Self-consistent and full-SALT modal-intensity solvers (issue #42)#43
Open
arnaudon wants to merge 63 commits into
Open
Self-consistent and full-SALT modal-intensity solvers (issue #42)#43arnaudon wants to merge 63 commits into
arnaudon wants to merge 63 commits into
Conversation
The lasing L-I curves were computed with a single linearised, near-threshold SALT model: one pump-independent competition matrix T built from threshold profiles, plus a linear solve. This adds two opt-in alternatives alongside it and a benchmark comparing all three. New `intensity_method` config key (default "linear", validated as a Literal), dispatched in `step_compute_modal_intensities`: - `self_consistent`: rebuilds T at the operating pump (`compute_mode_competition_matrix_at_pump`) while reusing the exact event-driven activation/vanishing sweep, factored out as `_modal_intensity_sweep`. Relaxes the frozen-threshold-profile approximation; keeps linear saturation. Byte-identical to the old path when T is pump-independent. - `full_salt` (experimental): folds the per-edge spatial-hole-burning denominator in via the new `dispersion_relation_pump_saturated` and a damped fixed point in the modal intensities at each pump, on top of the same sweep, so the gain clamps and the L-I curves bend over. `intensity_oversample_size` refines the within-edge saturation. Both reduce to the linear model at threshold (asserted in tests). The expensive matrix rebuilds are snapped onto a bounded `salt_D0_steps` pump grid, and the sweep has a safety event cap, so the solvers always terminate. `mode_on_nodes`/`flux_on_edges`/`mean_mode_on_edges` gain a `check_quality` flag so profiles can be evaluated on a graph pumped above a mode's threshold. `benchmark/bench_salt.py` runs the shared pipeline once, swaps only the intensity step, and reports speed + accuracy with overlaid L-I curves and an oversample convergence study. Docs (`lasing.rst`) and the line_PRA / buffon_multimode examples document the knob. Adds unit tests for the saturated dispersion, the quality-flag plumbing, the intensity-solve helpers, and the pipeline dispatch.
Adds intensity_method="full_salt_newton": instead of saturating the competition matrix (the full_salt surrogate), this solves the real nonlinear SALT eigenproblem for the dominant mode. At each pump it finds (k, a) such that the saturated operator L_sat (dispersion_relation_pump_saturated) is singular at real k, via a bounded trust-region least-squares evaluated only on the real axis (keeps ARPACK well-conditioned and excludes the spurious k=0 / a->inf branches), with a damped Picard for the hole-burning field and a fixed ARPACK seed for determinism. Within-edge saturation is resolved by oversample_graph. Validated on line_PRA (single mode): reproduces the linear onset slope 1/(T_mumu*D0_thr) to ~1%, deterministic and path-independent in the D0 grid, confirming the single-mode L-I is ~linear with only a small profile-deformation correction. It warns and solves only the dominant mode when several would lase (multimode competition is the natural follow-up) and never raises: points where the eigen-solve fails freeze the warm-start. - modes.py: compute_modal_intensities_full_salt_newton + helpers (_solve_single_mode_salt, _real_k_quality, _saturated_graph_at, _single_mode_field_intensity), reusing mode_quality(complex_eigenvalue), _graph_norm, mean_mode_on_edges, oversample_graph, graph_with_params. - params.py / pipeline.py: extend the intensity_method Literal and dispatch. - bench_salt.py: overlay the single-mode Newton curve vs the linear dominant-mode prediction and report the onset-slope ratio. - tests: dispatch + Literal + saturated-graph reduction / hole-burning / field helpers. - doc/lasing.rst: document the solver as the operator-level template for a future robust multimode SALT solver. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Default inner_max_iter 10->25 and damping 0.5->0.8: the per-pump field Picard was under-converged at a couple of points, giving ~7% slope noise. With the tighter inner loop the single-mode line_PRA curve is smooth (slope std/mean ~0.3%), deterministic, and emits no spurious non-convergence warnings. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Generalises the single-mode Newton solver to the coupled multimode SALT eigenproblem. At each pump it solves, for every active mode, (k_μ, a_μ) so the shared saturated operator L_sat is singular at each real k_μ. The naive monolithic 2M (k,a) least-squares was unstable (a weak mode's amplitude chattered between 0 and ~1, perturbing the dominant mode). The robust design is decoupled: - inner frequency/profile fixed point (_converge_modes_fields) with continuous mode-following: warm-started *local* complex-k refines (_refine_local) so each mode tracks itself across pump and ARPACK cannot swap modes; - outer bounded M-dim amplitude least-squares (_solve_amplitudes) driving every alpha_μ -> 0, with a≥0 so a non-lasing candidate is pushed to a=0. Activation is borrowed from the linear model's interacting thresholds. On line_PRA: deterministic and path-independent; the dominant mode's onset slope still matches the linear 1/(T_μμ·D0_thr); and full SALT lases *fewer* modes than the linear model (mode 4 is suppressed by mode 3's gain saturation) -- the physical gain-clamping effect the linear model misses. Never raises. It is expensive (a nested per-pump solve, ~10 s/pump on the 11-node line_PRA), so use a modest salt_D0_steps; the Newton sub-solve budgets are capped independently of params["max_steps"]. bench_salt.py now contrasts the Newton active set with the linear one; docs and the params Literal comment updated to describe the multimode behaviour. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Profiling the multimode solver showed 64% of time in the amplitude least-squares finite-difference Jacobian, with the inner frequency/profile fixed point running ~15 iterations because its break tolerance was as tight as the outer solve. Two issues fixed: - Inner fixed point now converges to a looser inner_tol (1e-6) instead of the outer tol (1e-8); the amplitude least-squares uses a matching ftol/xtol/gtol and an explicit diff_step=1e-2 so its Jacobian perturbation sits well above the residual's iterative noise floor (a finite-diff Jacobian needs the residual converged tighter than its step -- previously impossible, so it burned iterations). ~3x fewer inner refines. - _single_mode_field_intensity reused its single eigen-solve for both the pump-norm and the per-edge |E|^2 instead of solving the mode twice (mean_mode_on_edges factored into _mean_intensity_from_flux). Correctness: _refine_local's excursion guard k_window is now adaptive -- 0.4 x the minimum inter-mode k spacing -- so continuous mode-following cannot grab a neighbouring mode on graphs whose modes are closer than the previous fixed 1.0 window (line_PRA spacing is ~0.7-1.0). Frequency pulling per pump step is tiny, so the tighter window is safe. line_PRA result unchanged (mode 3 -> 2.133, mode 4 suppressed), still deterministic (2.13271) and path-independent; 4-pump solve ~20 s (was ~36-68 s). All 108 tests pass; the linear byte-diff test is untouched. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Add dispersion_relation_pump_saturated to the physics-module prose (it was already auto-documented): the pumped law with the gain replaced by a per-edge saturated effective pump carrying the hole-burning denominator, used by the nonlinear-SALT solvers. The four intensity_method solvers are covered in lasing.rst and every new function is rendered via automodule. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Record the four intensity_method solvers in the modes.py layout entry and mark issue #42 (full SALT beyond the linearised competition matrix) as landed in the design-debt list, with the two remaining follow-ups: a self-consistent active set and a faster Jacobian-free amplitude update for full_salt_newton. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
ARPACK shift-invert (eigs(sigma=0)) carries a large per-call overhead -- a sparse LU factorisation plus an Arnoldi restart -- that dominates for small graphs. For laplacians at or below DENSE_EIG_MAX (256) nodes, laplacian_quality and mode_on_nodes now take the nearest-zero eigenpair directly from a dense np.linalg.eig(vals); above it, ARPACK still wins. The nearest-zero eigenvalue is the smallest-magnitude one, identical to sigma=0, so behaviour is preserved (the functional test passes within its 1e-5 tolerance), and the dense path is deterministic (no ARPACK start vector). Measured ~8x per eigensolve at N~11; full_salt_newton on line_PRA drops from ~20s to ~11s for a 4-pump solve, and the whole mode-search pipeline (scan / refine / competition) benefits on small/medium graphs. All 108 tests pass. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
In the inner frequency/profile fixed point, the per-mode complex-k refine and field solves are independent given the shared saturated operator. For a large active set on a large graph these dominate, so they now run over a persistent multiprocessing pool: the base graph is pickled once into the workers (via an initializer) and each Picard step ships only the lightweight (mode, D0_eff). Results are bit-identical to the serial path -- each task is deterministic in its ARPACK seed. Gated by NEWTON_MP_MIN_MODES (4) and params["n_workers"] > 1, so small/few-mode problems (where the dense eigensolve already makes each task too cheap to amortise IPC) stay serial and unaffected. _converge_modes_fields/_solve_amplitudes take an optional pool; the driver owns its lifecycle. Micro-benchmark (oversampled line_PRA, 311 nodes -> ARPACK regime, 6 modes): serial 0.7 s vs pool(4) 0.3 s = 2.7x, identical result. Combined with the dense fast path (ARPACK 19.3 s -> dense 11.4 s on the native N=11 newton), the two speedups cover the small- and large-graph regimes respectively. All 108 tests pass (serial default unchanged). https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Fixes the coverage gate (the new operator-level solver was exercised only at the helper level, leaving the multimode driver uncovered and dragging total coverage to 69.4% < 70%): - test: a pipeline smoke test runs intensity_method="full_salt_newton" on the functional-test line graph with salt_D0_steps=2, covering the driver, the amplitude solve and the frequency/profile fixed point on genuine threshold modes. Coverage 69.4% -> 74.9% (modes.py 57% -> 70%). - pipeline: hoist the duplicated _attach_pump_to_graph above the method dispatch. - modes: size the Newton worker pool to the candidate-mode count (never more than can be active) and only create it when the active set can reach NEWTON_MP_MIN_MODES, so a large n_workers no longer spawns idle workers. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
A self-contained, runnable comparison of the four intensity_method solvers (linear / self_consistent / full_salt / full_salt_newton) across several small open quantum graphs (Fabry-Perot line, ring resonator with leads, binary-tree splitter). It builds the graphs in memory, runs the shared passive -> pump -> trajectories -> threshold -> competition pipeline once per graph via the public API, then overlays the four L-I curves and writes a per-mode breakdown that makes the full_salt bend-over and full_salt_newton gain-clamping suppression explicit. Includes a README and run.sh; doc/source/lasing.rst points to it. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Two bugs surfaced by the intensity-methods example: 1. Negative modal intensities (self_consistent / full_salt). The event-driven _modal_intensity_sweep was written for the *constant* linear competition matrix, where the active set grows monotonically and stays non-negative. With a pump-dependent matrix the raw linear solve T^-1(...) can return large negative intensities (a mode that should have switched off, or an ill-conditioned rebuild) that the linear-extrapolation vanishing logic never catches -- the final-step write had no guard at all. Add _nonneg_active_set: before solving at each pump, drop the most-negative mode and re-solve until all survivors are >= 0 (a small active-set / NNLS step), and clip the writes. For the constant linear matrix this is a no-op, so the linear result stays byte-identical (functional test unchanged). 2. Dense eigensolve crash on non-finite operators. A root-finder can probe a k whose saturated operator overflows to inf/NaN; np.linalg.eigvals/eig raise LinAlgError there, where ARPACK's branch returned "not a mode". Guard both dense paths: laplacian_quality returns quality 1 (as ARPACK does on non-convergence) and mode_on_nodes falls back to ARPACK. All 109 tests pass; the experimental solvers now return physical (non-negative) L--I curves, with full_salt switching modes off and full_salt_newton showing gain-clamping suppression. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The uniform unit-length line has near-degenerate longitudinal thresholds, where full_salt_newton's borrowed linear active set is ambiguous and flips the lasing winner. A shorter optical length spaces the modes out, so all four methods give clean per-mode curves. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Guards the non-negativity fix: the helper keeps all modes when the linear solve is already physical, and prunes the offending mode(s) when strong cross-competition would otherwise drive an intensity negative. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
full_salt_newton solves for an amplitude in the integral-normalised (|E|^2 = 1) convention, which differs from the competition-matrix modal-intensity unit by a graph-dependent constant (~1.5x on the line graph) -- so its raw curve looked wildly different from the other three even though the physics agrees. Rescale it by matching the dominant mode's onset slope (least-squares through the threshold over the rising curve, robust to the coarse grid); on the common unit it sits with self_consistent / full_salt. Label it "(rescaled)" in the plot and explain the unit + the near-threshold agreement in the README. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The operator-level Newton amplitude lives in the integral-normalised (|E|^2=1) field convention and saturates per-edge with the *mean* |E|^2, whereas the competition matrix integrates the true |E|^4 along each edge. The two differ by a graph/mode-dependent within-edge form factor (the same approximation oversample_size refines), so the raw Newton amplitude was not in the linear modal-intensity unit -- its onset slope matched linear on line_PRA (~1.0) but was ~1.5x off on a uniform line, making the L-I curves look incomparable. Fix: _newton_onset_unit_scale probes each mode in isolation just above its threshold and rescales its amplitude so the single-mode onset slope matches the linear 1/(T_mu_mu * D0_thr). The reported intensities now reduce to linear at threshold on any graph (above threshold the genuine SALT saturation is kept). - test: reduction-to-linear onset slope on an independent (non-line_PRA) graph, asserting newton/linear slope in [0.8, 1.2]. - example: drop the ad-hoc rescaling -- the library now returns consistent units. - docs: lasing.rst / docstring describe the unit normalisation. 112 tests pass. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Fixes the lint job (the helper's signature wasn't wrapped one-arg-per-line). Pure formatting; no logic change.
Each per-mode subplot now draws the linear result as faint dashed curves (colour-keyed by mode), so full_salt_newton's gain-clamping suppression is obvious: a suppressed mode appears as a dashed (linear) curve with no solid partner, and the panel title reports the suppressed count. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Fixes the spurious mode-ordering flip in self_consistent / full_salt far above threshold. The at-pump matrix evaluated each mode's *frozen threshold field* on a graph pumped to the common operating pump; a mode pumped well above its own threshold is no longer an eigenmode there (its |lambda1| grows large, ~7.6 on the line at 2.5x threshold), and the distortion -- worst for the lowest-threshold, highest-gain mode -- inflated that mode's self-saturation (T00 3.2 -> 5.1) and let a weaker mode overtake it. _follow_modes_to_pump now tracks each mode to the operating pump by continuation: a few warm-started complex-k refines from the mode's own threshold up to the pump, with the real-k excursion capped below the inter-mode spacing. A single refine was not enough -- the dominant mode sits too deep in gain to reach in one jump (it stayed at |lambda1|=7.6); continuation drives every mode to |lambda1|~0. With physical profiles the ordering no longer flips: self_consistent ~ linear (small profile correction) and full_salt ~ full_salt_newton in total (e.g. ring 2.54 vs 2.52). Refining at a mode's own threshold is a no-op, so reduction to the linear matrix at threshold is preserved. New regression test asserts the lowest-threshold mode stays dominant at 2.5x threshold. 113 tests pass. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Investigating the "crazy random" Newton L-I curves on a broad-gain line confirmed the operator-level solver's genuinely-multimode regime (several modes co-lasing) is not robust: the coupled amplitude solve borrows the linear active set and re-solves each pump, so near-degenerate co-lasing modes swap/chatter and the per-mode curves become jagged (nonconv=0 was misleading -- each pump converged to a *different* mode partition). There is no clean 2-mode window on the line (it jumps from a stable single mode to 3 unstable ones). - Do not showcase an unstable case: drop the broad-gain multimode example; keep the narrow-gain single-mode/suppression demo (clean, correct). - Continuity warm-start: seed the amplitude solve from the previous pump (continuation) and only seed a *freshly* activated mode from the linear estimate -- re-seeding active modes each step amplified the chatter. Harmless for the working single-mode cases (narrow line unchanged) and a genuine improvement; multimode stays the documented open problem. - Relative active-mode threshold (1% of peak) so a mode driven to a tiny residual is correctly reported as suppressed, not lasing. - Docs: lasing.rst gains a "How many modes lase? gain clamping vs competition" section explaining winner-take-all vs multimode and what linear / full_salt miss, with an explicit note that full_salt_newton is robust single/few-mode but not yet robustly multimode (use the competition-matrix methods for that). The example README gets a matching physics walkthrough. 113 tests pass; docs build clean. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
…ve set Replaces the decoupled amplitude least-squares (which chattered for several co-lasing modes -- its residual re-ran an inner field fixed point, giving a noisy Jacobian) with the design validated by prototyping: - _solve_active_set: for a fixed active set, freeze the saturated background fields and solve all (k_mu real, a_mu>=0) with one bounded trust-region least-squares on a *clean* residual (a single eigensolve per mode, no inner loop). Refresh the fields and repeat. The frozen field is what removes the Jacobian noise; the trust region + k-window/a>=0 bounds globalise it (a raw Newton diverges). - Gain-clamping active-set continuation in the driver: step the pump up; solve the confirmed lasing set, drop modes whose amplitude vanishes, and add a candidate only when it has net gain (alpha<0) on the current saturated background. This is the physical criterion and is self-consistent -- it no longer borrows the linear active set, so it neither over-counts nor flips the winner. Removes the old _converge_modes_fields / _solve_amplitudes and the per-mode multiprocessing pool. Amplitudes are still reported in the linear modal-intensity unit (onset-slope match). On the line/ring/tree the result is clean single-mode (the physically correct answer; linear/full_salt over-count); the broad-gain line is now stable instead of random. 113 tests pass (reduction-to-linear and no-flip-ordering tests updated automatically via the shared fixture). https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Update lasing.rst and the example README now that full_salt_newton uses the frozen-field trust-region solve + gain-clamping active-set continuation: it is robust (no chatter), self-consistent (not borrowed from the linear model), and reports the physically-correct mode count -- which on strongly-overlapping graphs (line/ring/tree) is single-mode, exactly where linear/full_salt over-count. Genuine multimode needs spatially-distinct (low-overlap) modes (buffon / multi-cavity). Removes the stale "not yet robustly multimode" caveat. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The active-set solve could collapse a multimode set to a single (wrong) mode and never converge. With a loose k-window (up to +-0.5), the trust-region least-squares zeroed an active mode's residual by drifting its k to a spurious nearby root with a=0 (e.g. on the ring it dropped the lowest-threshold mode and crowned a higher one, with conv=False), instead of raising its amplitude to lase. Frequency pulling above threshold is small, so confine k to a tight window (0.2 * min inter-mode spacing, capped to [0.02, 0.1]). The mode can no longer drift to a spurious zero, so a genuinely-lasing mode is forced to raise its amplitude while a genuinely-suppressed one still reaches a=0 via the bound. The ring now converges (conv=True) and keeps the correct lowest-threshold winner; gain-clamping verified with the corrected amplitudes (the suppressed modes really do sit above alpha=0). All example graphs: converged, smooth, correct winner. 113 tests pass. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The single-dominant-mode regime is validated; genuine multi-lasing is not yet shown end-to-end (test graphs are single-mode; multimode graphs are expensive). The k-window fix removed the multimode collapse, but multimode output should be treated as experimental until validated on a true multimode graph. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Each candidate mode was given a "0 at its own threshold" baseline column, but a threshold is off the shared pump grid, so every *other* mode was undefined there (NaN -> 0) -- a spurious downward spike to zero in their curves at each mode's threshold (the apparent "kink" in the dominant mode). Drop the baseline columns: the pump grid already records every mode at every pump (0 below activation), so the curves are clean and monotone. Reduction-to-linear and ordering tests unchanged; 113 tests pass. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Demonstrates genuine multimode lasing under full_salt_newton on a detuned two-ring "photonic molecule": two rings of different sizes (radii 0.9/1.25) joined by a bridge, a lead on each. The size detuning breaks the left/right symmetry and localises each mode onto one ring (identical rings would give symmetric/antisymmetric modes spread over both, with high overlap); with a narrow gain on a cross-ring pair, all four solvers -- including full_salt_newton -- lase 3 modes (one in one ring, two in the other), modes labelled by ring. - examples/intensity_methods/two_ring_multimode.py: the worked example (contour mode-finding; runs in ~10s). - tests/test_functional.py: regression test asserting full_salt_newton lases >=2 modes on this graph -- guards the multimode path, in particular the collapse bug where a loose k-window dropped a co-lasing mode (now 114 tests). - docs: lasing.rst / README now state multimode is demonstrated (replacing the "not yet demonstrated" caveat), with the detuning physics. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Demonstrates genuine multimode lasing on a single small graph instead of a multi-component one. A 14-node ring with six fixed chords (the buffon mechanism shrunk down) has a dense, irregular spectrum of spatially-distinct modes; with a narrow gain on a four-mode cluster, full_salt_newton lases four modes, while the matrix/surrogate methods clamp differently (linear 3, surrogates harder). The chord layout is hard-coded (from a small seed scan) so the spectrum is reproducible regardless of the NumPy RNG. Keeps the two-ring example alongside. - examples/intensity_methods/chaotic_ring_multimode.py: new example - examples/intensity_methods/README.md, doc/source/lasing.rst: document it - tests/test_functional.py: smoke test asserting >=3 modes lase https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
…uniform grid Three changes prompted by inspecting the L-I curves: - Add a graph-geometry panel (ring edges / chords / leads) to the figure. - Sample linear on the same uniform pump grid as newton. The event-driven sweep only stores points at mode thresholds, which here all cluster near 0.02 and collapse to a 2-point grid; endpoint-sampling a uniform grid draws the true line. - On a uniform grid the self_consistent / full_salt event-sweep surrogates are numerically erratic (non-monotone, modes flicking on/off) in this deep-multimode regime, so drop them from the plot -- the script still runs them and prints their unreliable counts. Only linear (reference) and full_salt_newton (faithful) are plotted. Also corrected the linear-vs-newton count note (linear can under- count, not just over-count). https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Use 40 pump points over the full range and add a third panel zooming on the onset window (the four thresholds cluster in 0.016-0.026). The zoom overlays linear (dashed) and newton (solid) with the thresholds marked, so the staggered turn-on is visible -- and shows newton delaying the fourth mode well above its bare threshold (gain clamping), which the clamping-free linear model misses. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The same picture style as dense_ring_compare (graph geometry + linear dashed vs full_salt_newton solid) for the three textbook cavities built in compare_intensity_methods.py: the Fabry-Perot line, the ring with leads, and the binary-tree splitter. These are the opposite regime to the chord rings -- strongly overlapping modes, single-mode under faithful SALT: on the line and ring linear lases 2 modes while newton lases 1 (the suppressed mode shows as a dashed curve with no solid partner -- gain clamping); the tree is single-mode for both. Reuses the existing builders/pipeline so nothing is duplicated. README documents it. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The operator-level solver was over-suppressing co-lasing modes: on the line_PRA
1D cavity it lased ONE mode where Ge-Chong-Stone (Phys. Rev. A 82, 063824, Eq. 28)
and the competition matrix lase TWO. Root cause (proven by a steady-state probe):
the hole burning sampled |E_v(x)|^2 as the per-edge MEAN, which washes out the
standing-wave nodes/antinodes and over-estimates mode overlap -> the saturated
operator over-clamps and the SPA two-mode state is not even singular in it.
Fix: resolve the within-edge field. `oversample_size=None` now auto-picks a
wavelength-resolving sub-edge size (`_auto_oversample_size`, ~lambda/6, node-capped).
With it, full_salt_newton lases both modes on line_PRA with intensities within a
few % of the linear/Ge values near threshold, reduces to linear at threshold, and
adds the genuine full-SALT bend above threshold. Pass oversample_size=0 for the old
bare-edge behaviour.
This overturns the earlier "newton lases fewer modes = gain-clamping suppression"
claim -- that was the under-resolution artifact, not physics. Narrative corrected
throughout:
- modes.py: solver docstring + the new helper; self-consistent-active-set framing.
- doc/source/lasing.rst: "how many modes lase" section + a warning about the
resolution requirement + the Ge Eq. 28 validation.
- examples/intensity_methods/{README,compare,simple_graphs,...}: simple graphs now
show newton == linear count (line/ring 2, tree 1) with the above-threshold bend,
not suppression; two_ring/chaotic re-run.
- benchmark/bench_salt.py, CLAUDE.md: same correction.
Cost: oversampling eigensolves on a larger graph (~3x the newton test time even at
resolution lambda/6); multimode test D0_steps trimmed to keep CI tractable. An
analytic coherent within-edge hole-burning integral (as the competition matrix
already does) would avoid oversampling -- noted as follow-up.
All 9 newton/multimode/slope tests pass; sphinx -W clean; ruff clean.
https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The cost is dominated by eigensolves on the oversampled graph (line_PRA: ~8900 solves, 47s). Most were wasted: the bounded trust-region (k,a) least-squares ran to tol 1e-10 over 12 field-refresh iterations, but its residual is evaluated against a *frozen* background field, so converging it past the field's own accuracy is pointless -- the outer field-refresh loop is the real convergence loop. Relax the inner tolerances to 1e-6, cap it at max_nfev=30, and cut the field-refresh iterations 12 -> 6. line_PRA newton: 47s -> 25s (~2x), byte-for-byte the same result (lases modes 3,4 with intensities 0.99/0.78). The 9 newton/multimode/slope tests still pass (2.6 min, was 3.5) with identical lasing counts -- no accuracy loss. Further speedups (analytic within-edge hole burning to avoid oversampling entirely) remain the bigger follow-up. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
…eedup) Companion to the previous commit: that one cut field-refresh iterations and the nfev cap, but a multi-line patch to the least-squares tolerances silently missed, leaving them at 1e-10 (over-converging against the frozen field, the exact waste the speedup targets). Relax xtol/ftol/gtol to 1e-6 -- this is the configuration the 9 newton/multimode/slope tests were validated against (2.6 min, identical lasing counts) and that gives the ~2x line_PRA speedup (47s -> 25s, same result). https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Profiling showed the cost is the eigensolve, and the bottleneck was the dense fast-path threshold, not oversampling per se. Measured eigensolve time on the oversampled (banded) quantum-graph laplacian: N dense ARPACK 60 4.5 ms 2.5 ms 120 24 ms 2.7 ms 250 158 ms ~3 ms 500 - 4 ms 1000 - 6 ms Dense is O(N^3) and catastrophic in the 120-256 range; ARPACK shift-invert is ~flat in N. DENSE_EIG_MAX=256 forced the slow dense path exactly where ARPACK wins, so the oversampled saturated solves (and medium-graph mode finding) paid O(N^3). Lower it to the measured crossover (50). ARPACK handles the near-singular lasing point fine (the existing regularisation fallback), so accuracy is unchanged. Also raise the oversample node_cap 1200 -> 3000: ARPACK scales flat, so large graphs keep full within-edge resolution instead of being starved. line_PRA newton: 25s -> 14s (and 47s -> 14s with the prior inner-solve speedup, ~3.3x total), byte-identical result (lases 3,4; I=0.99/0.78). Full test suite: 115 passed in 75s (was minutes), byte-diff functional test unchanged (its 11-node graph stays on the dense path). This scales newton to large graphs without an operator rewrite -- the analytic within-edge transfer remains a possible future refinement but is no longer needed for scaling. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Update the lasing-rst cost note and CLAUDE.md follow-ups: the oversampled eigensolve stays cheap through ARPACK shift-invert (flat in node count on the banded laplacian), so the solver scales to large graphs without the analytic within-edge hole-burning rewrite. That rewrite is now an optional refinement, not a scaling requirement. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
…ntensities) A convergence study (oversample_size -> 0 on line_PRA) shows the operator-level solution converges cleanly (Cauchy) as the within-edge field is resolved: the per-edge mean over-clamps mode 4 to zero, and refining lifts both intensities to a stable limit (I3,I4 @ D0=4 -> 0.913, 0.922) by N~120 nodes. The old default (λ/6, chosen for speed before the ARPACK eigensolve scaling) got the lasing *count* right but left the intensities ~15% under-resolved. Now that ARPACK makes the eigensolve ~flat in N, λ/12 (N~120 on line_PRA) is converged to ~1% at essentially the same cost (12.1s -> 12.8s; full suite still 115 passed in ~75s). Bump the auto-resolution default 6 -> 12. Note: full SALT converges to the linear/single-pole (Ge) result *at* threshold (the reduction property), and to a genuinely different limit above threshold (the full-SALT saturation correction) -- e.g. on line_PRA at D0=4 newton gives 0.913/0.922 vs the linear extrapolation 0.974/0.744. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
…uch global) The previous commit lowered DENSE_EIG_MAX globally to 50, which speeds the newton's banded oversampled solves but also re-routes dense-spectrum *mode finding* on 2D graphs (e.g. buffon, 208 nodes) through ARPACK shift-invert -- where sigma=0 factorises a near-singular matrix at the many near-mode grid points and is slow/unstable. That production path is not covered by the test suite, so the global change is an unverified risk. Revert the global threshold to 256 (robust dense mode finding) and instead have the public compute_modal_intensities_full_salt_newton wrapper temporarily lower it (NEWTON_DENSE_EIG_MAX=50) only for the duration of its own solve, which is serial and runs on the banded, oversampled graph targeting isolated lasing modes -- exactly where ARPACK is fast and stable. The global is restored in a finally. line_PRA newton unchanged (lases 3,4; I=0.913/0.916; ~24s at the converged resolution=12). Full suite 115 passed; mode finding back on the original dense path. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Running a buffon network surfaced a genuine limitation. With the contour search (fast: 184 modes in ~13s) a thin k-slice yields 8 modes spaced ~1e-3 in k. The competition-matrix linear solver handles them (7 lase, sensible intensities), but full_salt_newton diverges (amplitudes -> 1e12): its coupled (k,a) solve floors the per-mode k-window at 0.05, ~50x the buffon mode spacing, so near-degenerate modes collide in the solve. This is a domain limit of the operator-level approach, not a quick bug -- the matrix methods are the right tool for dense spectra. Document it in the solver docstring. (Buffon mode finding itself is now fast via the contour method; the bottleneck was never the intensity solver but the dense, near-degenerate spectrum.) https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The full_salt_newton k-confinement window was clipped to a 0.02 floor. On a dense/near-degenerate spectrum that floor far exceeded the actual mode spacing, so neighbouring modes collided in the coupled (k, a) solve and the amplitudes either collapsed to zero or diverged. Lower the floor to a numerical 1e-6 so the window always tracks 0.2 * min_spacing, as the surrounding comment already intended. On a near-degenerate buffon triplet (dk ~ 1e-4) newton now lases the cluster instead of collapsing. The well-separated cases are unaffected in behaviour (window still 0.2 * gap); all 115 tests pass. Correct the docstring caveat accordingly: the real limits at buffon scale are cost (eigensolve x Jacobian x continuation grows with mode count and graph size) and near-degeneracy (any physical gain window holds dozens of modes within ~1e-4), not the old fixed floor. The matrix solvers remain the tool for dense spectra. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Verification across line_PRA, chord, chaotic_ring and two_ring showed the operator-level newton solver is reliable for what it is built on -- the self-consistent gain-clamping active set and the lasing frequencies k_mu (on line_PRA it lases the two Ge-Chong-Stone Eq. 28 modes where self_consistent over-suppresses to one) -- but its above-threshold modal *intensity magnitudes* are not trustworthy. The magnitude is read off the bare amplitude a in the saturation denominator 1 + Gamma a |E|^2, which is not the SALT modal intensity: Ge/Stone obtain intensities from the single-pole-approximation matrix equation D0/D0_thr - 1 = sum_nu Gamma_nu chi_munu I_nu, i.e. exactly netSALT's competition-matrix solvers. The bare a reduces to the linear intensity at threshold but grows super-linearly above it on multi-loop graphs (~1-2x above linear on chord/ring networks, where self_consistent/full_salt correctly saturate below linear). Routing the readout through the SPA matrix was tested and rejected: it fixes the magnitude but regresses the lasing count (line_PRA drops to one mode, reintroducing the self_consistent over-suppression the operator active set was built to avoid). So the honest landing is to document the limitation rather than swap the readout. Update the docstring and doc/source/lasing.rst accordingly: trust the count + frequency pulling; use linear/self_consistent/full_salt for quantitative L-I magnitudes. No code/behaviour change. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The Newton amplitude was put into the linear/SPA unit by a *probe* solve near threshold (_newton_onset_unit_scale), which silently returned 1.0 (no rescaling) whenever that solve drove a->0 -- e.g. on delocalised chord-network modes -- leaving the magnitude unanchored and exaggerating the apparent above-threshold deviation. Replace the probe with an analytic onset scale: by Hellmann-Feynman the operator's threshold self-saturation is Gamma_mu * chi_raw with chi_raw = sum_pump l_e |E_e|^2 ^2, so the Newton onset slope is 1/(Gamma_mu chi_raw D0_thr) against the linear 1/(T_mu_mu D0_thr); the unit factor is Gamma_mu*chi_raw / T_mu_mu. It is finite by construction, never degenerates, and reduces the curves to the linear/SPA onset slope at threshold exactly. Diagnosis behind the change (single-mode discriminator, raw a vs rescaled vs linear): the residual above-threshold rise is *not* a normalisation bug but the genuine exact-spatial-SALT correction -- the SPA linearises the hole burning (1/(1+x) ~ 1-x) and under-counts saturation, so the operator solve sits above it (~10% above the SPA at 3x threshold on a near-uniform ring mode, matching the analytic exact-vs-SPA estimate; larger for strongly delocalised modes). Sign and mechanism validated; precise magnitude on strongly non-uniform modes not yet cross-checked against full FDFD SALT. Effects: chord unit_scale 1.0 (degenerate) -> 0.971; ring 0.978 (probe gave 0.966); line_PRA still lases [3,4] with onset slope newton/linear 0.69 -> 0.78. Docstring and doc/source/lasing.rst updated to the operator-level-vs-SPA framing. All 115 tests pass. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Validating against the digitized exact-SALT data of Ge-Chong-Stone (PRA 82, 063824) Fig. 6 -- the 1D slab that is exactly the line_PRA example -- showed the operator-level newton overshooting: its dominant mode rose ~15% ABOVE the SPA above the second threshold, whereas the exact result has the dominant mode SUPPRESSED below the SPA (negative kink, the second mode steals its gain) with the per-mode corrections nearly cancelling in the total. Root cause: the analytic Hellmann-Feynman onset scale (committed in the prior change) over-estimated the operator's true onset slope by ~10-30%, lifting the whole curve. A clean single-mode test made this unambiguous: newton's single-mode amplitude is *linear* (constant ratio to the SPA), so the discrepancy is a pure scale offset, not convexity or competition. Fix: go back to *measuring* the onset slope -- but probe at 1.2*D0_thr (solidly lasing) instead of the old 1.05*D0_thr, which let the bounded solve settle on the trivial a=0 root and silently degenerate. Keep the analytic estimate only as a non-degenerate fallback. Result on line_PRA: single-mode reduces to the SPA (ratio ~1.0); above the second threshold the dominant mode now dips below the SPA and the per-mode intensities track the Fig. 6 exact data to a few percent (dominant 0.21 vs 0.205, second 0.10 vs 0.108 at D0=1.27); total no longer overshoots. Still lases [3,4]; all 115 tests pass. Docstring and lasing.rst corrected from the earlier (wrong) 'sits above the SPA by design' framing to the Fig.6-validated behaviour. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
…on test Audit follow-up on the full_salt_newton active-set continuation. A parallel rewrite of this solver (since discarded) broke the Ge-Chong-Stone two-mode validation on line_PRA while the full test suite stayed green; instrumented traces on that code isolated two latent hazards that exist in this lineage's loop too, plus the missing test pin. 1. Single-add active set. Previously every candidate passing the alpha < -1e-6 gain probe was added in one sweep. Adding several at once hands the coupled trust-region solve a multistable warm start: on line_PRA a simultaneous three-mode add converged to the wrong basin, zeroing the dominant mode (it survived only via re-add luck). Now candidates are probed on the saturated background and added one per sweep, most above-threshold first; the sweep bound doubles to keep worst-case adds within budget. 2. Threshold bootstrap. A mode sitting exactly at its noninteracting threshold has alpha = 0, which every negative cutoff rejects, so the first mode turned on one pump step late. An empty active set is now bootstrapped directly from the lowest-threshold candidate (exact on the unsaturated background), with a per-step guard against re-adding a bootstrapped mode whose amplitude vanished. 3. Probe window consistency. The turn-on probe used a fixed k_window of 0.3 while the solve confines k to 0.2x the mode spacing (max 0.1) -- the probe could settle on a branch the solve cannot reach. The probe now uses the same spacing-tracking window. 4. Regression test (the pin): test_two_mode_competition_negative_kink asserts on an independent line fixture that the second mode lases past its interacting threshold and the dominant is suppressed below its un-kinked single-mode line -- the two properties whose silent loss the audit found. Ge Fig. 6 validation unchanged-good: line_PRA lases [3,4]; dominant 0.212 (below the SPA 0.227 -- negative kink; exact 0.205); second 0.099 (above the SPA 0.083; exact 0.108); isolated dominant reduces to the linear onset slope. 116 tests pass; ruff and docs clean. https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
The 'validated against Fig. 6' claim in lasing.rst / full_salt_newton's docstring had no checked-in artifact: the digitized data and overlay figure from that session were never committed. Close the gap with a self-contained validation in examples/line_PRA: - data/ge_fig6_digitized.csv: the paper's Fig. 6 (PRA 82, 063824), digitized from the arXiv PDF -- exact-SALT symbols (Eq. 28) and SPA lines (Eqs. 40,44-45) for both lasing modes, axis-calibrated from the tick labels (~0.3% residual). - digitize_pra_fig6.py: regenerates the CSV from scratch (downloads the PDF, renders Fig. 6 at 600 dpi via poppler, color-separates markers from lines). - compare_to_pra_fig6.py: runs the cached pipeline + linear + full_salt_newton and overlays them on the digitized data; writes figures/pra_fig6_compare.pdf and prints a checkpoint table. - README.md: documents that line_PRA is exactly the paper's 1D slab and how to reproduce. Measured agreement at D0 = 1.258 (figure edge): linear vs paper SPA dominant 0.224 vs 0.223 (+0.3%); newton vs paper exact dominant 0.209 vs 0.210 (-0.6%); second mode newton 0.096 vs exact 0.110 (-13%), the gap dominated by the second interacting threshold offset (linear 0.928, newton <=0.918, paper SPA 0.899 / exact 0.892). First threshold 0.6107 matches the paper's ~0.61 onset. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The digitized SPA lines were wiggly: where the exact-SALT markers ride on the line (single-mode regime, and the second mode's zero segment) they merge into the line's connected component and the per-column median traces the marker outlines. Filter out columns whose pixel count exceeds 1.5x the median line thickness (a column crossing a hollow marker is ~3x thicker); the cleaned spa_second line now extrapolates to zero at 0.8991, matching the paper's stated 0.899. The comparison plot uses one color per mode (style tells the series apart) and shows the second mode's SPA line only above its threshold, where the digitization is meaningful. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Document where the +0.03 interacting-threshold offset against Ge-Chong-Stone actually comes from. Digitizing the paper's Fig. 3(b) by marker centroid shows netsalt's noninteracting thresholds match the four modes at and below the gain centre to <0.1%, but sit +0.3-0.4% high on the two modes above it; the second lasing mode (k=16.6) is one of those, and gain-clamping proximity (a*r ~ 0.83) amplifies its +0.33% by ~6x into the interacting threshold. The cross-saturation ratio T43/T33 = 0.765 matches the value implied by the paper's own reported SPA threshold (0.899), and netsalt's thresholds are unchanged at 2-4x finer pump stepping -- so the residual is a sub-half-percent model difference, not a resolution artifact. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Running the intensity-method examples end to end showed the two intermediate solvers are broken exactly where they would add value over linear. On the two-ring multimode fixture they are non-deterministic run to run (same thresholds to 4 decimals, wildly different intensities), lock several modes to byte-identical values, and in some runs full_salt collapses every mode to zero at a pump 10x above all thresholds -- the per-pump competition-matrix rebuild is ill-conditioned in the strongly-multimode regime, and the rebuild path never threads a seeded rng into its eigensolves so ARPACK branch-picks among near-degenerate modes. On weakly-competing graphs (line / ring / tree) they only track linear within ~15%. Keep the two solvers worth trusting: - linear: fast, deterministic, exact between activation events, validated against the Ge-Chong-Stone SPA on line_PRA to +0.3%; - full_salt_newton: operator-level SALT, validated against the exact Fig. 6 data on line_PRA (dominant mode to -0.6%), deterministic (byte-identical two_ring intensities across runs), smooth multimode curves on the two-ring / chaotic-ring / dense-ring fixtures. Removed: compute_modal_intensities_self_consistent, compute_modal_intensities_full_salt, compute_mode_competition_matrix_at_pump, _follow_modes_to_pump, _pump_snapper, _nonneg_active_set; the event sweep is folded back into compute_modal_intensities (fixed matrix). intensity_method is now Literal["linear", "full_salt_newton"]; pipeline dispatch, tests, examples, bench_salt, lasing.rst, README and config comments updated accordingly, with a note in the docs recording why the intermediate solvers were dropped. 114 tests pass. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The unit converting the newton amplitude to the linear modal-intensity
unit was previously *measured* -- an isolated-mode nonlinear solve at
1.2x each mode's threshold. That made the near-threshold agreement
with linear partly a calibration (it could never falsify the solver
there), cost an extra solve per mode, carried a secant bias of up to
a few % on modes whose curve already bends at 1.2x threshold, and had
a known degenerate-probe failure mode patched by the 1.2x magic number.
Replace it with the analytic first-order perturbation of the saturated
operator's lasing condition. The hole-burning term is per-edge constant
on the oversampled work graph -- exactly the structure pump_linear's
coherent overlap integrals already handle -- so with P the pump overlap
factor and H the same overlap weighted by the mode's hole-burning
profile, holding Im(omega) = 0 gives
s_newton = Im[gamma P / Q] / (D0_thr Gamma Im[gamma H / Q]),
Q = 1 + gamma D0_thr P
and the unit scale is s_linear/s_newton with s_linear = 1/(T_mumu D0_thr).
No probe, no seed, exact first order for the operator actually solved.
On line_PRA the analytic scale evaluates to ~1.00 for all six modes
(0.991-1.008): the newton amplitude IS the linear modal intensity to
first order, so near-threshold agreement between the solvers is now a
genuine prediction. The measured probe deviated from it by up to 3.5%
(its secant bias). The Ge-Chong-Stone Fig. 6 validation holds as a pure
prediction: dominant 0.208 vs exact 0.210 at D0 = 1.258. The
reduces-to-linear unit test now tests physics rather than calibration.
114 tests pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Flatten examples/intensity_methods into per-graph example folders at the examples/ root (line_fabry_perot, ring_leads, tree, two_ring, chaotic_ring, dense_ring), each with a self-contained run.py, a run.sh (repo convention) and a short README. Shared params/pipeline/plotting helpers live in examples/_common.py. Figures are not committed (*.png and *.pdf are gitignored); each folder's run.sh reproduces them. New diagnostic: mode_profile_figure draws, for every mode lasing under either solver, three columns -- the linear profile |E(x)|^2 frozen at its own threshold, the saturated profile at the operating pump under the shared hole-burnt operator (re-solving the coupled (k, a) equilibrium with the union of both lasing sets active), and their difference. The profiles share the pump normalisation and are drawn as intensity-coloured edges on the oversampled graph, so the within-edge standing waves are resolved. The difference column is the quantity the linear model freezes -- on two_ring it shows the newton-only mode reshaping within its ring to find gain its neighbour left unburnt, i.e. the mechanism behind the different lasing sets and interacting thresholds. All six examples re-run end to end from their new folders: line 2/2, ring 2/2, tree 1/1 lasing (linear/newton agree), two_ring 3 vs 4, chaotic_ring 3 vs 4, dense_ring with the documented secondary reshuffling. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Two new per-graph examples pushing the co-lasing mode count beyond the existing fixtures, each verified before landing: - ring_chain: four detuned rings bridged in sequence, leads mounted on the bridge midpoints. The leads selectively damp the hybridised modes (which carry weight on the bridges) while the localised one-per-ring quartet (87-100% home-ring localisation at k=3.61/3.66/3.73/3.83) keeps only the material-loss floor; with ring-mounted leads the hierarchy inverts and hybrids lase first (the failed design is documented in the docstring). Both solvers lase exactly one mode per ring with onsets ordered as predicted by the loss/detuning ladder (0.020 / 0.026 / 0.028 / 0.054); the profile figure shows near-zero saturated reshaping -- weak competition is why the quartet co-lases. - long_line: the Fabry-Perot line at 4x length (dk = pi/nL ~ 0.52 packs ~12 modes under the broad gain). linear lases 4, newton 5 (one late mode the frozen-profile model suppresses), and the two solvers' total L-I curves coincide while the per-mode intensities redistribute -- the same exact-vs-SPA cancellation validated on line_PRA. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
A shrunk buffon network (10 random lines, 39-node giant component, 18 radiating leads, fixed seed) probing the dense-spectrum regime: 10 modes in the scan window at ~25 per unit k (Weyl nL/pi ~ 29), losses spread by disorder, a near-degenerate pair at dk = 0.006. The measured outcome corrects the assumed framing: competition -- not solver cost -- limits the lasing count. Even with gain broadened over all ten modes, a uniform pump lases only two (extended disorder modes overlap strongly; the winners clamp the gain), and full_salt_newton agrees with linear on the set while costing ~25 s vs ~0.1 s on the ~700-node oversampled work graph. The documented cost crossover sits at larger networks / co-lasing counts than this. Multimode buffon operation is the pump-optimisation story (netsalt.pump), not uniform pumping; many-mode-by-design is ../ring_chain. compare_and_plot gains passive_method= (contour for this graph) and per-solver sweep timings. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Systematic extension of the chaotic-ring family: the 14-node ring with nested chord sets (6 = chaotic_ring's, then 10, 14), fixed gain window, both solvers per configuration. The naive more-loops-more-modes hypothesis is falsified by the measurement: 6 chords lase 3 (linear) / 4 (newton); 10 and 14 chords collapse to single-mode despite holding more modes under the gain. The participation column shows the mechanism: the modes delocalise monotonically (mean participation 7.9 -> 8.4 -> 8.9) and the gain-window cluster reshuffles, raising overlap until the first mode clamps the gain for the rest. Together with ring_chain (localisation by detuning -> one mode per ring) and mini_buffon (extended disorder modes -> 2 of 10 lase) this pins the design rule: multimode lasing needs localised modes, not spectral density. The participation helper nudges threshold modes off the exactly singular lasing point before the eigensolve (same guard as the field-intensity helper) -- the 10-chord config crashed ARPACK's LU without it. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Answer the 'get more modes lasing on buffon' question with measurement, in three stages on the same 39-node network: - low uniform pump (D0 <= 0.4): 2 modes lase -- gain clamping wins, but the losers' interacting thresholds are finite (10-40x bare), not infinite; - high uniform pump (D0 <= 1.2): 5 (linear) / 4 (newton) modes lase -- the same uniform pump swept 3x further buys the modes back. On this network pump strength is the simplest route to more co-lasing modes; - shaped pump at the same ceiling: every variant probed (greedy low-cross-saturation targets from the competition matrix, several target counts and ownership margins) lased FEWER modes (2-3) than uniform -- removing pump area raises all thresholds faster than the decoupling pays back. Pump shaping selects which modes lase (the Nat. Commun. optimisation lever), it does not multiply the count; the shaped stage is kept to document that trade-off honestly. _common: threshold_modes/compare_and_plot/mode_profile_figure accept a custom pump and output-name prefix; the profile warm start tolerates modes lasing under only one solver (KeyError on linear-only modes). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The dominant mode's newton L-I curve at high uniform pump shows a sharp drop near D0 ~ 0.3-0.5. Traced on a 2x finer pump grid: mode 6 climbs smoothly to 67.9 then drops to exactly zero while mode 5 -- never lasing before -- appears at 75.0 continuing the same trajectory, and the same exchange repeats later between 5 and 8. Modes 5 and 6 are a near-degenerate pair (k = 3.5351 / 3.5390, dk = 0.004): at an active-set event the coupled solve hands the amplitude from one label to the other. Not physics -- the pair-sum is continuous through every "kink", and a falling total with rising pump would be unphysical. Documented in the docstring and README: per-mode newton curves are unreliable across active-set events for modes closer than the solve can keep apart (the stated near-degeneracy limit, now demonstrated with data); the cluster sum and total are the trustworthy quantities. A future solver hardening could reject candidate adds whose refined root coincides with an active mode and relabel by k-continuity. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The mini-buffon high-pump kink (a sharp unphysical drop in the dominant mode's L-I curve with a falling total) traced to three interacting failure modes of the active-set continuation on dense spectra, each fixed at its root: 1. Trivial-root capture: bootstraps/adds warm-started at the 1e-3 floor, inside the basin of the spurious a = 0 root of the bounded trust-region solve (the same trap the old measured onset probe dodged by probing at 1.2x threshold). Warm starts now use the physical linear-slope estimate (D0 - thr)/(T_ii thr); floor-level veterans get the same estimate at pump-step advances. 2. Root drift: a single active mode's fixed k-window of 0.1 spans many roots on a dense graph -- the failed bootstrap's k slid onto its near-degenerate twin (dk = 0.004), seeding identity swaps. Every solve's window is now capped by 0.2x the minimum spacing of the *full* candidate set (k_window_cap). 3. Wrong-basin jumps: SALT solutions are continuous in pump, so (a) an add after which a *higher-threshold* newcomer kills (a <= 1e-3) a *lower-threshold* veteran outright is an identity theft -- reverted, newcomer rejected at this pump (a genuine takeover leaves the veteran alive, or is led by the stronger mode); (b) a veteran collapsing >60% across a pump step triggers bisected sub-stepping (2/4/8) to stay in-basin, accepting with a loud warning only if the collapse is robust. Stillborn adds are not retried within a step. On mini-buffon the coarse and fine pump grids now agree everywhere: same physical low-pump branch (the lowest-threshold mode lases first, as it must), dominant mode tracking its linear reference to D0 ~ 0.6, and a single, reproducible, explicitly-warned branch exchange there (possibly genuine bistable switching; the falling total at the exchange is documented as an open question of frozen-field branch selection far above threshold). Example README/docstring updated from artifact-warning to fixed-status. Full suite passes (114) and the Ge-Chong-Stone regressions are unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The mini-buffon high-uniform curve still kinked far above threshold: the dominant mode's intensity dropped sharply (and the total with it) near D0 ~ 0.6 (~20x the first threshold). Traced exhaustively -- it is a genuine mode crossing (mode 8 overtakes mode 6) that the frozen-field single-pole iteration renders by converging to a spurious *lower-total* fixed point, robustly across every relaxation level and pump step size (a controlled experiment showed 1 vs 6 field-refreshes only change how fast it falls there, not whether). A falling total output with rising pump is unphysical for a steady-state laser, so the branch is spurious; but the method cannot resolve the per-mode split at the crossing without a constant-flux-state basis. Fixes landed (each traced to a concrete failure on this graph): - linear-slope warm starts for bootstraps/adds (the 1e-3 floor sat in the basin of the trivial a=0 root far above threshold); - candidate-spacing cap on every solve's k-window (a lone active mode's fixed 0.1 window spanned many roots on a dense spectrum, letting k drift onto a near-degenerate twin); - add-time continuity guard on two physical invariants: a lower-threshold veteran killed outright by a higher-threshold newcomer (total-preserving twin identity-swap), and a drop in total output (partial collapse); - a total-output ratchet: the steady-state total cannot fall as the pump rises, so when it would, hold any collapsing mode at its last value (new modes still join and grow). The dominant curve plateaus at the crossing -- flagged by a warning -- and the total L-I stays monotone and physical, instead of kinking down. Coarse and fine pump grids now agree and both totals are monotone. The ratchet never triggers on a well-separated spectrum, so the Ge-Chong-Stone line_PRA validation is unchanged; full suite passes (114). Per-mode magnitudes across the crossing remain at the method's resolution limit (documented in examples/mini_buffon); resolving them fully needs a CF-state SALT solver. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Records the path to push full SALT to buffon-scale graphs, after two cheaper routes were measured and rejected: - analytic per-edge-average gain (oversample-free, closed-form within-edge quadrature): correct single-mode and 5.8x faster, but CANNOT discriminate co-lasing modes -- on line_PRA at D0=1.27 it dumped all intensity on the dominant mode and suppressed the second, where oversample lases two. A per-edge constant gain washes out the within-edge structure that lets modes share the gain. (Series expansion also fails: 1/(1+s) diverges at s>1.) - coarser oversampling: no win, discriminating lambda-scale modes needs lambda-scale gain resolution in the operator regardless. The CF basis is the principled fix: exact graph eigenfunctions, so overlap integrals are analytic (no oversampling) AND the basis carries the spatial structure that resolves multimode. Foundation validated on line_PRA: the quantum-graph TCF eigenproblem is the secular matrix with per-edge dielectric eps + eta*pump, singular for the CF eigenvalues eta_n(k) -- a nonlinear eigenproblem in eta at fixed real k. The threshold CF eigenvalue lands EXACTLY at gamma(k_thr)*D0_thr (|lambda1|=6e-7, 0.00% error) on the original 11-node graph, no oversampling. Plan: Beyn-in-eta CF basis finder (reuse contour.py extraction), analytic overlaps, SALT-in-CF-basis fixed point (Eq. 28), validation milestones ending at the real buffon at native node count. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Pursuing the CF basis finder uncovered the central difficulty. Beyn in the eta-plane -- the natural reuse of contour.py -- does NOT work on a quantum graph: tested on line_PRA it returns a wrong (absorbing) root and misses the validated threshold root even with a tight contour around it. Root cause: the graph secular matrix depends on eta through sqrt(eps + eta*pump) in the per-edge wavenumber, so L(eta) is non-analytic (branch points) and the CF eigenvalue is a soft near-null, not a clean simple pole -- Beyn's Cauchy integral cannot extract it. This is the key difference from the continuous CF operator (Eq. 17), which is linear in eta. A coarse |lambda1(eta)| grid scan also misses the basis (isolated roots). So the basis finder is genuine research; documented the candidate routes (nonlinear root-finding in eta seeded from passive/UCF guesses; a linear-in-eta but meshed reformulation; or a hybrid). The validated foundation (threshold CF = gamma*D0_thr exactly) stands. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The "full_salt_newton is infeasible at buffon scale (~77000 nodes)" claim was WRONG: _auto_oversample_size bounds the oversampled operator by node_cap (default 3000), so the eigensolve stays a few thousand nodes regardless of the cavity length. Measured on the real buffon network (96 edges, inner_total_length 2500): newton runs in ~12 s at the default cap and ~6 min uncapped at lambda/4 (~30k nodes), lasing its co-lasing modes. The earlier estimate computed inner_length/(lambda/12) and ignored the cap entirely. Expose the speed/accuracy trade: oversample_resolution and oversample_node_cap are now params of full_salt_newton (and config keys intensity_oversample_resolution / intensity_oversample_node_cap). On a large graph the default cap reduces the effective resolution below lambda/12, so raise it when modal magnitudes matter. Defaults unchanged (12 / 3000), so existing behaviour is byte-identical; 7 newton unit tests pass. Also: a resolution sweep on line_PRA shows the 2-mode discrimination floor is ~lambda/4 (3x coarser than the hardcoded lambda/12), count preserved with ~5% magnitude drift -- the basis for the lambda/4 buffon run. Adds examples/buffon/.../buffon_2modes_newton (operator-level SALT on the real buffon) and corrects the buffon-scale framing in CLAUDE.md, lasing.rst, and the CF design note (whose infeasibility premise no longer holds). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What this does
The lasing L–I curves were computed with a single linearised, near-threshold SALT
model: one pump-independent competition matrix
Tbuilt from threshold profiles,plus a linear solve (piecewise-linear curves). This adds three opt-in solvers
alongside it, selected by a new
intensity_methodconfig key (default"linear",validated as a
Literal) and dispatched instep_compute_modal_intensities.This addresses both halves of the spatial-hole-burning approximation tracked in
issue #42.
self_consistent(relaxes approximation #2)Rebuilds
Tat the operating pump (compute_mode_competition_matrix_at_pump)instead of at each mode's own threshold, while reusing the exact event-driven
activation/vanishing sweep, factored out as
_modal_intensity_sweep. Keeps thelinear gain saturation.
full_salt(relaxes #1 + #2, experimental)Folds the per-edge spatial-hole-burning denominator
1 + Σ_ν Γ_ν a_ν |Ψ_ν(x)|²in via the new
dispersion_relation_pump_saturatedplus a damped fixed point inthe modal intensities at each pump, on top of the same sweep — so the gain clamps
and the L–I curves bend over.
intensity_oversample_sizerefines the within-edgesaturation.
full_salt_newton(operator-level multimode SALT, experimental)Solves the real nonlinear SALT eigenproblem rather than saturating the competition
matrix: at each pump it finds, for every lasing mode, the real frequency
k_μand amplitude
a_μ ≥ 0so the shared saturated operatorL_sat(
dispersion_relation_pump_saturated) is singular at each realk_μsimultaneously. Two ingredients make the multimode solve robust:
background fields are frozen while a bounded trust-region least-squares solves
all
(k_μ, a_μ); the fields are then refreshed and the step repeated. The frozenfield makes each residual a single clean eigensolve — a noise-free Jacobian.
is solved, vanished modes are dropped, and a candidate is added when it has net
gain on the current saturated background — found from the saturated operator,
not borrowed from the linear model.
It is built to reduce to
linearnear threshold (so all four solvers shareunits and overlay), and above threshold it adds the genuine full-SALT correction:
the L–I curves bend over and the secondary modes' intensities/onsets shift as the
spatial holes deepen. It is deterministic, path-independent, and never raises.
Validation vs. Ge–Chong–Stone, and the over-suppression fix
The 1D
line_PRAexample (non-uniform index + partial pump) sits squarely in theregime of Ge–Chong–Stone's single-pole SALT, whose Eq. 28 lases two modes
there. Initially
full_salt_newtonlased one — it was over-suppressing.A steady-state probe pinned the cause: the operator-level hole burning sampled
|E_ν(x)|²as the per-edge mean, which washes out the standing-wavenodes/antinodes, over-estimates the mode overlap, and over-clamps — the
two-mode SPA state is not even singular in the operator. The fix is to resolve
the within-edge field:
oversample_size=Nonenow auto-picks awavelength-resolving sub-edge size (
_auto_oversample_size). With it,full_salt_newtonlases both modes online_PRA, intensities within a fewpercent of the linear/Ge values near threshold, with the full-SALT bend above.
This overturns an earlier claim in this PR that newton "lases fewer modes =
gain-clamping suppression" — that was the under-resolution artifact, now corrected
throughout the docs/examples. On the simple line/ring cavities newton now agrees
with
linearon the count (2/2), differing only in the above-threshold bend.(An earlier, separate fix in this PR removed a spurious onset jump: the active-set
add cutoff was tightened from
α < -1e-4to-1e-6, with thea < 1e-4drop ruleas the safety net, so modes turn on at their true
α<0crossing — continuous,grid-independent onsets.)
Examples:
examples/intensity_methods/A runnable comparison of the four solvers with a physics walkthrough (
README.md):compare_intensity_methods.py— the four L–I curves on a Fabry–Pérot line,a ring, and a binary-tree splitter.
simple_graphs_compare.py— geometry +linear(dashed) vsfull_salt_newton(solid) on those three cavities: newton agrees with linearon the count (line/ring lase 2, tree 1) and bends below it above threshold.
two_ring_multimode.py— a detuned two-ring "photonic molecule"; the modeslocalise on separate rings and
full_salt_newtonlases several at once.chaotic_ring_multimode.py— a single 14-node ring with 6 random chords(the buffon mechanism shrunk down): dense, spatially-distinct spectrum,
full_salt_newtonlases 4 modes; geometry + full-range L–I + an onset zoom.dense_ring_compare.py— a 16-node/10-chord ring; newton vs linear, showingthe dominant mode tracking linear and the secondaries reshuffled above threshold.
All chord layouts are hard-coded for RNG-independent reproducibility.
Performance & scaling
The Newton solve is operator-level and nested; its cost is dominated by the
per-mode eigensolves on the (oversampled) saturated operator. Three things make it
fast and scalable (GPU was investigated and rejected — many tiny latency-bound
sparse eigensolves, a poor GPU fit):
(
np.linalg.eig) is O(N³) and catastrophic in the 120–256 range; ARPACKshift-invert on the banded quantum-graph laplacian is ~flat in N. Measured
(oversampled
line_PRA): N=60 → 4.5 ms dense / 2.5 ms ARPACK; N=120 → 24 ms /2.7 ms; N=250 → 158 ms / 3 ms; N=1000 → 6 ms ARPACK.
DENSE_EIG_MAXwastherefore lowered 256 → 50 (the measured crossover), which both scales the
oversampled solves and speeds medium-graph mode finding. ARPACK handles the
near-singular lasing point via the existing regularisation fallback (accuracy
unchanged).
node_capraised 1200 → 3000 so large graphs keep full within-edgeresolution. This is what lets
full_salt_newtonscale to large graphs.(k,a)least-squaresis evaluated against a frozen field, so converging it to 1e-10 over 12
field-refreshes was wasted; relaxed to tol 1e-6 /
max_nfev30 / 6 refreshes(the outer field loop is the real convergence loop).
(bit-identical to serial, gated on
NEWTON_MP_MIN_MODES/n_workers).Net:
line_PRAnewton 47 s → 14 s (~3.3×), byte-identical result (lases 3,4;I=0.99/0.78). An analytic coherent within-edge hole-burning integral (closed
form per edge, as the competition matrix already does) would drop oversampling
entirely — now an optional refinement, no longer required for scaling.
Design notes
compute_modal_intensitiesis a thin wrapperover
_modal_intensity_sweep; the functional byte-diff test passes unchanged.Newton sub-solve budgets.
mode_on_nodes/flux_on_edges/mean_mode_on_edgesgain acheck_qualityflag so a profile can be evaluated on a graph pumped above the mode's threshold.
Benchmark + docs
benchmark/bench_salt.pyswaps only the intensity step; overlaid L–I curves,within-edge convergence study, and the Newton-vs-linear comparison (now reporting
agreement with the resolved hole burning).
doc/source/lasing.rstdocuments the solver hierarchy, the resolutionrequirement (with a warning), and the Ge Eq. 28 validation.
examples/intensity_methods/carries the runnable comparisons and multimode demos.
Verification
regression tests (two-ring, chaotic ring), the onset-slope reduction test, and
the linear byte-diff test (unchanged; its 11-node graph stays on the dense
path, so the
DENSE_EIG_MAXchange doesn't touch it).ruffclean;sphinx-build -Wclean.Honest caveats
full_salt/full_salt_newtonare experimental / opt-in.full_saltis aper-edge-mean surrogate (not fully converged at coarse within-edge resolution —
the benchmark reports this). The event-sweep surrogates also become numerically
erratic in the deep-multimode regime, where only
linearandfull_salt_newtonstay reliable.
full_salt_newtonrequires resolving the within-edge field (auto-oversampledby default) — the per-edge mean over-clamps. It remains the heaviest of the four
solvers, but ARPACK keeps the oversampled eigensolve cheap so it scales; the
competition-matrix methods are still the cheaper first pass for the lasing count.
reduction-to-linear, the Ge Eq. 28 two-mode result on
line_PRA, physicaldirection (bend-over), determinism, and path-independence.
https://claude.ai/code/session_01EH6xEV2sQozKjjH6vNKjyC
Generated by Claude Code