Every item of the original roadmap is now implemented, each pinned to exact anchors (66 tests total, up from 47; CI on Python 3.9 through 3.13).
Solitons and soliton crystals by verified Newton continuation (newton_state, soliton_seed, continuation): localized steady states of the LLE, solved with the exact discrete form of the package's own fluctuation matrix as the Jacobian -- the operator that damps a perturbation and the operator Newton inverts are the same linearization. The solver refuses to return anything whose stationary residual is not verified below tolerance, and the tests close the loop through independent code paths: a flat seed lands on the exact cubic root; the converged soliton barely moves under the independent split-step evolver, sits on the exact lower-branch background, and carries the exact translation (Goldstone) zero mode with eigenvector d psi / d theta; a two-pulse crystal equals, grid point for grid point, the single soliton of the equation with dispersion scaled by four (the exact rescaling of the periodic domain); and a truncated stable mode set around the soliton shows sideband squeezing below vacuum.
Marginal vs unstable, distinguished: the Goldstone mode makes a soliton's drift matrix marginally rather than asymptotically stable, so the spectra now refuse a marginal matrix with an explanation unless allow_marginal=True says the translation mode is understood -- while a genuinely positive growth rate stays refused regardless of the flag.
Supermode decomposition (principal_quadratures): the eigendecomposition of any multimode covariance matrix, whose smallest eigenvalue is -- exactly, by linear algebra -- the deepest squeezing any generalized quadrature of the state attains, with its supermode as eigenvector. Pinned to the two-mode squeezed vacuum closed form ((hbar/2) e^(-/+2r), EPR supermodes), flat on vacuum, and verified to lower-bound a scan of quadratures on the photonic molecule. Deliberately distinct from the Williamson spectrum, which measures mixedness.
Thermal input noise: Bose occupations for the extraction-port and loss baths in output_quadrature_variance, per-mode occupations in intracavity_covariance, and thermal_occupation computing the physical number from the exact SI constants (h and k_B are exact by definition since 2019). Anchors are closed forms: a passive cavity with baths at n_bar emits exactly (2 n_bar + 1)/2 at every frequency, coupling and phase; the parametric oscillator's whole spectrum scales by exactly (2 n_bar + 1); a passive mode holds exactly n_bar photons; the hot-loss/cold-port mixture matches its hand-derived form; and n_bar crosses 1 exactly at hbar omega = k_B T ln 2.
Future-proofing: a CHANGELOG covering every release, and the CI matrix widened to Python 3.9, 3.11, 3.12 and 3.13.
The scope that remains -- non-Gaussian quantum dynamics, thermorefractive/Raman technical noise, pulsed pumping, Bloch-Messiah for pure states -- is named in the README as deliberate boundaries with their reasons, so a designed-out feature can be told from a gap.
Install: pip install sqzcomb -- archived on Zenodo (DOI badge in the README).