Metadata-Version: 2.4
Name: sqzcomb
Version: 0.15.0
Summary: Squeezed light in Kerr microcombs: Lugiato-Lefever steady states, linearized quantum noise, and detected output spectra
Author-email: "Tanvir M. Mahim" <tanvir.mahim@bracu.ac.bd>
License: Apache-2.0
Project-URL: Homepage, https://tanvir-mahmud-mahim.github.io/software/
Project-URL: Repository, https://github.com/TaN-MM-Org/sqzcomb
Project-URL: Issues, https://github.com/TaN-MM-Org/sqzcomb/issues
Project-URL: Changelog, https://github.com/TaN-MM-Org/sqzcomb/releases
Keywords: squeezed light,microcomb,Kerr,Lugiato-Lefever,quantum optics,homodyne
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.22
Requires-Dist: scipy>=1.10
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Provides-Extra: interop
Requires-Dist: qutip>=4.7; extra == "interop"
Dynamic: license-file

# sqzcomb

[![PyPI](https://img.shields.io/pypi/v/sqzcomb)](https://pypi.org/project/sqzcomb/) [![DOI](https://img.shields.io/badge/DOI-10.5281%2Fzenodo.22015375-blue)](https://doi.org/10.5281/zenodo.22015375) [![tests](https://github.com/TaN-MM-Org/sqzcomb/actions/workflows/ci.yml/badge.svg)](https://github.com/TaN-MM-Org/sqzcomb/actions)

`sqzcomb` is a Python package for predicting **squeezed light from
Kerr microcombs**: tiny optical ring resonators that, when pumped with
a laser, can produce light whose noise in one quadrature is below the
normal quantum noise floor. The package follows the light end to end:
from the classical field circulating in the ring, through the quantum
noise around it, to **the noise a detector on the output fiber would
actually report**. That last step is the point. Squeezing inside the
ring is not what an experiment sees; only what leaves through the
output port is, and the ring coupling that is best for one is not best
for the other.

It answers questions such as:

- How much squeezing will my ring deliver to a detector, given its
  linewidth, output coupling and pump power?
- How much do optical loss, detector electronics and local-oscillator
  phase jitter take away, and what loss or jitter budget does a target
  need, with all three together?
- At which local-oscillator angle is the light most squeezed, at each
  analysis frequency? (Found exactly, without scanning angles.) And
  which one fixed angle is best for a whole band of frequencies?
- Can a second, coupled ring get the light out better than one ring
  alone can?
- Are two comb lines entangled at the output, and would a detector be
  able to show it?
- From a measured pair of squeezed and antisqueezed levels, or from
  measured noise spectra, what did the source itself produce? And how
  do raw spectrum-analyzer traces become dB relative to shot noise,
  with error bars? How much does an uncertain shot-noise reference
  add to those error bars?
- Before measuring: will the planned frequencies pin down the
  parameters at all?
- What happens above threshold, or when single photons matter and the
  light is no longer Gaussian?
- How much does classical noise (a wandering resonance, a noisy pump
  laser) spoil the squeezing?
- How squeezed is a pulse of light made with a pulsed pump?

When a question falls outside what the model can answer, the package
stops with an error that says why, instead of returning a number that
looks fine but is not (see [When it refuses, and
why](#when-it-refuses-and-why)).

## Contents

- [A short guide to the words used here](#a-short-guide-to-the-words-used-here)
- [Install and requirements](#install-and-requirements)
- [Units and conventions](#units-and-conventions)
- [Examples](#examples) (each with the output it prints)
- [What is in the package](#what-is-in-the-package)
- [When it refuses, and why](#when-it-refuses-and-why)
- [How the results are checked](#how-the-results-are-checked)
- [Corrections](#corrections)
- [Limits](#limits)
- [Where it comes from](#where-it-comes-from)
- [Citing, support and license](#citing-support-and-license)

## A short guide to the words used here

- **Microring, Kerr microcomb** -- a small ring-shaped optical
  resonator. Its material has the **Kerr effect** (the refractive index
  changes slightly with light intensity), which lets a strong pump
  laser create light at new, evenly spaced frequencies (a **comb** of
  **lines**, or **modes**).
- **Linewidth (kappa)** -- how fast light leaks out of the ring; the
  width of its resonance. **Escape efficiency** `eta = kappa_ex /
  kappa` is the fraction of that leakage that goes into the useful
  output port; the rest is lost inside the ring. `eta = 1/2` is
  **critical coupling**.
- **Pump, detuning, threshold** -- the pump is the driving laser. Its
  **detuning** is how far it sits from the ring's resonance. Above the
  **threshold** pump power the ring starts to oscillate on new comb
  lines. Most of this package describes the noise around a stable
  steady state (below threshold, or a soliton); the above-threshold
  tools are described in examples 11 and 12. In the model,
  `F` is the pump strength and `alpha` the detuning, both in
  normalized units; `mu` is the pump strength seen by the quantum
  noise (the "parametric gain"), with `mu = 1` the threshold of a
  single squeezed mode.
- **Dispersion** -- how the ring's resonance spacing changes from line
  to line (`D2` in Hz, `d2` in normalized units). "Anomalous"
  dispersion is what bright solitons need.
- **Lugiato-Lefever equation (LLE)** -- the standard equation for the
  field inside such a ring. `sqzcomb` solves it in its usual
  dimensionless ("normalized") form.
- **Soliton** -- a single stable pulse circulating in the ring; a
  **soliton crystal** is several equally spaced pulses.
- **Quadrature** -- light's noise can be split into two components,
  like the cosine and sine parts of a wave, selected by a phase `phi`.
  **Vacuum** (no light) has equal noise in every quadrature; this is
  the **shot-noise** level.
- **Squeezing / antisqueezing** -- noise below / above vacuum in a
  quadrature, quoted in **dB relative to vacuum** (negative is
  squeezed). The squeezed and antisqueezed quadratures are 90 degrees
  apart.
- **Homodyne detection** -- the standard way to measure one
  quadrature: the light is mixed with a strong reference beam, the
  **local oscillator (LO)**, whose phase selects the quadrature. The
  noise is recorded against **analysis frequency** (offset from the
  carrier, `omega` in normalized units, `f` in Hz).
- **Detection efficiency** -- the fraction of the light that reaches
  and is counted by the detector. Any loss mixes in vacuum noise and
  washes squeezing out.
- **Covariance matrix** -- a table of all the quadrature noises and
  their correlations; it fully describes the (Gaussian) states used
  here. **xxpp** is the ordering (all `x` quadratures, then all `p`).
- **Entanglement** -- quantum correlations between two lines stronger
  than any classical ones. **Logarithmic negativity** `E_N` measures
  it (0 means not entangled); it comes from the **PPT test** (positive
  partial transpose), a standard mathematical check for entanglement
  of such states. The **Duan test** is a simpler check based on the
  noise of sums and differences of the two lines' quadratures; it can
  miss entanglement that the PPT test finds.
- **Symplectic eigenvalues** -- numbers computed from a covariance
  matrix that say how noisy (mixed) a state is; all equal to the
  vacuum value means a pure state.
- **Photonic molecule** -- two coupled rings. Here the second ring is
  passive and can act as the output channel for the first. `J` is the
  rate at which light hops between the rings and `gamma` the decay
  rate of the second ring (both in units of `kappa/2`).
- **Drift matrix** -- the matrix `M` of the linear equations that the
  small quantum fluctuations obey. Almost every function here takes or
  returns one. It must be **stable** (all eigenvalues with negative
  real part) for a steady noise level to exist.
- **Gaussian state, linearization** -- the usual approximation for
  bright light: small noise around a large classical field, whose
  statistics are then fully set by the covariance matrix. It fails when
  individual photons matter.
- **Master equation, Fock basis, Wigner function** -- the full quantum
  description of a few modes: the state is a matrix over photon
  numbers 0, 1, 2, ... (the **Fock basis**, cut off at some maximum).
  The **Wigner function** is a picture of that state over the two
  quadratures; where it goes **negative** the state is not Gaussian
  and has no classical counterpart.
- **Technical noise** -- classical noise added by the hardware: the
  resonance frequency wandering with temperature
  (**thermorefractive noise**), pump-laser intensity and phase noise. It
  is given as a **power spectral density (PSD)**: noise power per unit
  bandwidth.
- **Temporal mode** -- with a pulsed local oscillator the detector
  measures one pulse-shaped mode of the output light, set by the LO's
  pulse shape.
- **Purity** -- 1 for a pure quantum state; lower when the state also
  carries classical randomness (a squeezed state with excess noise).

## Install and requirements

```
pip install sqzcomb              # numpy + scipy
pip install sqzcomb[interop]     # adds QuTiP, for drift_from_qutip only
```

For development: clone the repository and `pip install -e .[test]`.

It needs Python 3.9 or newer, NumPy 1.22 or newer and SciPy 1.10 or
newer. QuTiP (4.7 or newer) is needed only for `drift_from_qutip`.

## Units and conventions

- **Model units.** The physics runs in the normalized units of the
  LLE: rates in units of `kappa/2` (so a single ring's field decays at
  rate 1 and a single parametric mode reaches threshold at `mu = 1`),
  time in units of `2/kappa`, analysis frequency `omega` in units of
  `kappa/2`. The dispersion sign follows the LLE module: `d2 < 0` is
  anomalous.
- **Lab units.** `RingSpec` and its helpers take ordinary Hz (the
  linewidth `kappa_hz` is the full width at half maximum, `g0_hz` is
  `g0/2pi`), watts and metres, and convert once. `hbar` and `c` come
  from `scipy.constants`.
- **Noise level.** A quadrature variance of **0.5 is vacuum**
  (`squeezing_db(v) = 10 log10(v / 0.5)`). With a real pump parameter
  `mu > 0`, the quadrature at `phi = pi/2` is squeezed and `phi = 0`
  antisqueezed.
- **Covariance matrices** are in xxpp order with vacuum `(hbar/2) I`.
  `covariance_xxpp` and the entanglement functions use `hbar = 2` by
  default (vacuum = identity); `output_covariance_xxpp` returns
  `hbar = 1` (vacuum = 0.5 I, the same scale as the variances above).
- **Measured traces** (`fit_noise_spectra`, `infer_source`, the CSV
  files) are in dB relative to shot noise, the way a spectrum analyzer
  trace is calibrated. `shot_noise_normalize` makes them from the raw
  traces in dBm.
- **Dark noise and the reference level.** By default
  `detected_variance` adds the dark (electronic) noise to the signal
  and quotes it against the pure vacuum level. A laboratory ratio to a
  shot-noise trace taken through the same receiver carries the dark
  noise in the reference too; `dark_in_reference=True` gives that
  number (example 19).

## Examples

Each example below runs as written, and the output shown is what it
printed with sqzcomb 0.15.0 (examples 1 to 19 print the same as with
0.14.0). The device and pump numbers are
illustrative values chosen for the example, not measured data, except
the measured pair in example 6, which is quoted from the paper named
there.

### 1. From lab numbers to model units

```python
import sqzcomb as sc

# Illustrative ring values (not a real device).
ring = sc.RingSpec(kappa_hz=100e6,      # loaded linewidth (FWHM), Hz
                   eta_esc=0.8,         # escape efficiency kappa_ex / kappa
                   g0_hz=10.0,          # single-photon Kerr shift g0 / 2 pi, Hz
                   lambda_pump_m=1.55e-6,
                   reference="illustrative values for the README")

print(f"comb threshold at alpha = 1: {sc.threshold_power(ring) * 1e3:.3f} mW")
print(f"pump F at 2 mW:              {sc.normalized_pump(ring, 2e-3):.4f}")
print(f"detuning alpha for 50 MHz:   {sc.normalized_detuning(ring, 50e6):.4f}")
print(f"d2 for D2/2pi = 100 kHz:     {sc.normalized_dispersion(ring, 1e5):.4f}")
print(f"omega for a 10 MHz offset:   {sc.normalized_frequency(ring, 10e6):.4f}")
print(f"photons in the ring at rho=1: {sc.intracavity_photons(ring, 1.0):.3e}")
```

```
comb threshold at alpha = 1: 0.126 mW
pump F at 2 mW:              3.9870
detuning alpha for 50 MHz:   1.0000
d2 for D2/2pi = 100 kHz:     -0.0020
omega for a 10 MHz offset:   0.2000
photons in the ring at rho=1: 5.000e+06
```

The 2 mW pump is well above threshold for this ring (`F` = 3.99,
against `F = 1` at threshold for `alpha = 1`). `threshold_power` is the power at which the
flat-state cubic of the LLE has the root `|psi|^2 = 1`, so
`F_th^2 = 1 + (alpha - 1)^2`. Here `psi` is the normalized field in
the ring and `rho = |psi|^2` its normalized intensity. Every
`RingSpec` needs a `reference`
string saying where its numbers came from; an empty one is refused.

### 2. One ring against two: the 3 dB limit

```python
import numpy as np
from sqzcomb import (single_mode_parametric, output_quadrature_variance,
                     squeezing_db, photonic_molecule, molecule_threshold,
                     output_variance_ports)

# One ring, pumped at 99 % of threshold (mu = 0.99), noise read at the
# centre of the spectrum (omega = 0) in the squeezed quadrature (phi = pi/2).
M = single_mode_parametric(0.99)
for eta in (0.5, 0.8, 1.0):
    v = output_quadrature_variance(M, eta, 0.0, mode_index=0, n_modes=1,
                                   phi=np.pi / 2)
    print(f"single ring, eta = {eta}: {squeezing_db(v):7.2f} dB")

# Two rings: the Kerr ring has NO output port; light leaves through a
# second, passive ring (J^2 / gamma = 3), read at its own port.
J, gamma = np.sqrt(27.0), 9.0
mu_th = molecule_threshold(J, gamma)
M2, gammas = photonic_molecule(0.99 * mu_th, J, gamma=gamma)
v2 = output_variance_ports(M2, gammas, eta=1.0, port_mode=1, omega=0.0,
                           phi=0.0)
print(f"threshold of the pair: mu = {mu_th:.4f}")
print(f"two rings, read through ring 2: {squeezing_db(v2):7.2f} dB")
```

```
single ring, eta = 0.5:   -3.01 dB
single ring, eta = 0.8:   -6.99 dB
single ring, eta = 1.0:  -45.98 dB
threshold of the pair: mu = 4.0000
two rings, read through ring 2:   -6.02 dB
```

A single ring at critical coupling (`eta = 0.5`) cannot deliver more
than about 3 dB of squeezing to the detector, however hard it is
pumped: half of the squeezed light is lost inside the ring and replaced
by vacuum. Stronger output coupling does better. In the two-ring
**photonic molecule**, the pumped Kerr ring has no output port of its
own; its light leaves through the second ring, which acts as an output
channel with effective efficiency `(J^2/gamma) / (1 + J^2/gamma) =
3/4` here. The result passes the 3 dB value. Note the reading phase:
hopping to the second ring turns the squeezed quadrature by a quarter
turn, so it is read at `phi = 0` there.

### 3. A pumped comb: steady state, then detected squeezing

```python
import numpy as np
from sqzcomb import (lle_evolve, homogeneous_steady_states,
                     fluctuation_matrix, output_quadrature_variance,
                     squeezing_db)

F, alpha = 0.9, 0.3                          # illustrative, below threshold
psi = lle_evolve(np.full(16, 0.05 + 0j), F=F, alpha=alpha, t_end=100.0)
print(f"|psi|^2 from the solver: {abs(psi[0]) ** 2:.6f}")
print(f"root of the cubic:       {homogeneous_steady_states(F, alpha)[0]:.6f}")

M, modes = fluctuation_matrix(psi, alpha)    # quantum noise around it
i0 = int(np.where(modes == 0)[0][0])         # the pumped line
phis = np.linspace(0.0, np.pi, 181)
for eta in (0.2, 0.5, 0.8, 1.0):
    best = min(squeezing_db(output_quadrature_variance(
        M, eta, 0.0, i0, modes.size, phi=p)) for p in phis)
    print(f"eta = {eta}: best squeezing {best:6.3f} dB")
```

```
|psi|^2 from the solver: 0.698836
root of the cubic:       0.698836
eta = 0.2: best squeezing -0.730 dB
eta = 0.5: best squeezing -2.126 dB
eta = 0.8: best squeezing -4.194 dB
eta = 1.0: best squeezing -6.460 dB
```

`lle_evolve` finds the classical field by time-stepping the LLE; for a
flat state it lands on the root of the cubic `rho (1 + (alpha -
rho)^2) = F^2`. `fluctuation_matrix` then builds the quantum-noise
equations around it, and `output_quadrature_variance` gives the noise
at the output port. The script `examples/squeezing_vs_coupling.py`
does the same for six couplings from 0.1 to 1.0.

### 4. What the detector takes away

```python
from sqzcomb import (detected_squeezing_db, dark_from_clearance_db,
                     required_efficiency_db, phase_noise_squeezing_db,
                     max_phase_noise, thermal_occupation)

v_source = 0.05                        # 10 dB below vacuum (vacuum = 0.5)
print(f"70 % efficiency:            {detected_squeezing_db(v_source, 0.7):.3f} dB")
dark = dark_from_clearance_db(15.0)    # receiver noise 15 dB below pure shot noise
print(f"70 % and 15 dB clearance:   {detected_squeezing_db(v_source, 0.7, dark):.3f} dB")
print(f"efficiency needed for -6 dB: {required_efficiency_db(-6.0, -10.0):.4f}")

# Local-oscillator phase jitter mixes in the antisqueezed quadrature.
for jitter in (0.001, 0.03, 0.1):       # RMS, radians
    db = phase_noise_squeezing_db(-10.0, 11.0, jitter)
    print(f"-10 dB / +11 dB source, {jitter * 1e3:5.0f} mrad jitter: {db:.3f} dB")
print(f"largest jitter that still gives variance 0.1: {max_phase_noise(0.1, 0.05, 6.3):.4f} rad")

print(f"thermal photons at 193 THz and 300 K: {thermal_occupation(193e12, 300.0):.2e}")
```

```
70 % efficiency:            -4.318 dB
70 % and 15 dB clearance:   -3.962 dB
efficiency needed for -6 dB: 0.8320
-10 dB / +11 dB source,     1 mrad jitter: -9.999 dB
-10 dB / +11 dB source,    30 mrad jitter: -9.538 dB
-10 dB / +11 dB source,   100 mrad jitter: -6.504 dB
largest jitter that still gives variance 0.1: 0.0898 rad
thermal photons at 193 THz and 300 K: 3.90e-14
```

Loss mixes in vacuum: `V_detected = eta V + (1 - eta)/2 + V_dark`.
Local-oscillator phase jitter mixes in the antisqueezed quadrature;
for Gaussian jitter of RMS size `sigma` the average is the closed form
`<cos 2 theta> = exp(-2 sigma^2)`. The better the source, the larger its
antisqueezing and the more jitter hurts. At optical frequencies and
room temperature the thermal photon number is tiny, which is why the
baths default to vacuum. `dark_from_clearance_db(15.0)` reads the 15 dB
from the pure shot-noise level; for the gap between an analyzer's
shot-noise trace (which carries the dark noise too) and its dark trace,
pass `trace_gap=True` (example 19). Here the dark noise is quoted
against the pure vacuum level; against a shot-noise trace that carries the same
dark noise the number differs slightly (example 19).

### 5. Entanglement and supermodes

```python
import numpy as np
from sqzcomb import (photonic_molecule, intracavity_covariance,
                     covariance_xxpp, entanglement_report,
                     principal_quadratures, symplectic_eigenvalues,
                     output_entanglement_spectrum)

# The two coupled rings of example 2, driven at mu = 0.8 (illustrative).
M, gammas = photonic_molecule(mu=0.8, J=1.0)
sigma = covariance_xxpp(intracavity_covariance(M, gammas))   # vacuum = identity
rep = entanglement_report(sigma, 0, 1)
print(f"log. negativity: {rep['log_negativity']:.4f} "
      f"(entangled: {rep['entangled']})")
print(f"Duan sum {rep['duan_sum']:.4f}, bound {rep['duan_bound']:.1f} "
      f"(Duan test fires: {rep['duan_violated']})")
w, v = principal_quadratures(sigma)
print("principal variances:", np.round(w, 4))
print("symplectic eigenvalues:", np.round(symplectic_eigenvalues(sigma), 4))

# Output entanglement of a twin-beam pair, as a detector would see it.
twin = np.array([[-1, 0, 0, 0.5], [0, -1, 0.5, 0],
                 [0, 0.5, -1, 0], [0.5, 0, 0, -1]], dtype=complex)
omegas = [0.0, 1.0, 5.0]
print("E_N at the output:", np.round(output_entanglement_spectrum(
    twin, 0.8, omegas, 0, 1, 2), 4))
print("... with 70 % detection:", np.round(output_entanglement_spectrum(
    twin, 0.8, omegas, 0, 1, 2, detection_efficiency=0.7), 4))
```

```
log. negativity: 0.0911 (entangled: True)
Duan sum 2.6440, bound 2.0 (Duan test fires: False)
principal variances: [0.5795 0.9307 1.241  2.5368]
symplectic eigenvalues: [1.0484 1.0484]
E_N at the output: [1.2417 0.6779 0.0605]
... with 70 % detection: [0.6887 0.4225 0.042 ]
```

Inside the driven two-ring molecule the rings are entangled (`E_N >
0`), but the simple symmetric Duan test does not detect it; the
sharper PPT test used for `E_N` does. `principal_quadratures` finds the
most-squeezed combination of all quadratures (0.5795 here, against 1 for
vacuum in these units). `symplectic_eigenvalues` above 1 show the state
is mixed. The second part uses a hand-written drift matrix of a
two-line twin-beam source (`da1/dt = -a1 + 0.5 a2*` and the same with 1
and 2 swapped) and gives the entanglement a pair of detectors would
see at each analysis frequency, with and without 70 % detection
efficiency.

### 6. Inferring the source from a measured pair

```python
from sqzcomb import infer_source

# Measured pair of Ulanov et al., Nat. Commun. 16, 10791 (2025):
# 1.71 dB squeezing and 5.54 dB antisqueezing, after all losses.
out = infer_source(sq_db=-1.71, anti_db=5.54, sigma_db=0.1)
print(f"implied total efficiency:  {out['eta']:.4f} +/- {out['eta_sigma']:.4f}")
print(f"inferred source squeezing: {out['sq_db_source']:.2f} "
      f"+/- {out['sq_db_source_sigma']:.2f} dB")
print(f"inferred source antisqueezing: {out['anti_db_source']:.2f} dB")

try:
    infer_source(sq_db=-3.0, anti_db=1.0)
except ValueError as err:
    print("refused:", err)
```

```
implied total efficiency:  0.3724 +/- 0.0204
inferred source squeezing: -8.99 +/- 0.25 dB
inferred source antisqueezing: 8.99 dB
refused: the measured uncertainty product S*A = 0.157739 is below the vacuum product 0.25: impossible for any source followed by loss (loss only grows the product). Check the shot-noise calibration
```

Under the standard model -- a pure squeezed state followed by loss --
the measured squeezed and antisqueezed levels fix the total efficiency
and the source in closed form. This is a model statement, not a
measurement: if the source is less pure than assumed (excess phase
noise, thermal noise), the inferred numbers flatter it. The error bars
come from numerical derivatives of the closed form. The paper's own
inference is not reproduced here (the test suite checks only that this
pair gives a physical efficiency and a stronger source).

### 7. Fitting measured noise spectra, and planning the measurement

```python
import numpy as np
from sqzcomb import noise_spectrum_db, fit_noise_spectra, plan_noise_measurement

# Synthetic "measured" traces: illustrative true values plus 0.05 dB noise.
mu, eta, kappa_hz = 0.55, 0.62, 12e6
f = np.linspace(0.5e6, 30e6, 25)                  # analysis frequencies, Hz
rng = np.random.default_rng(1)
sq = noise_spectrum_db(mu, eta, kappa_hz, f, "squeezed") + rng.normal(0, 0.05, f.size)
an = noise_spectrum_db(mu, eta, kappa_hz, f, "anti") + rng.normal(0, 0.05, f.size)

fit = fit_noise_spectra(f, sq_db=sq, anti_db=an, sigma_db=0.05)
print(f"mu    = {fit.mu:.4f} +/- {fit.sigma['mu']:.4f}")
print(f"eta   = {fit.eta:.4f} +/- {fit.sigma['eta']:.4f}")
print(f"kappa = {fit.kappa_hz / 1e6:.3f} +/- {fit.sigma['kappa_hz'] / 1e6:.3f} MHz")
print(f"chi^2 / dof = {fit.chi2:.1f} / {fit.chi2_dof}")

# Before measuring: would one quadrature alone be enough?
one = plan_noise_measurement(mu, eta, kappa_hz, f, sigma_db=0.05,
                             quadratures=("squeezed",))
both = plan_noise_measurement(mu, eta, kappa_hz, f, sigma_db=0.05)
print("squeezed trace only identifiable:", one["identifiable"])
print("both traces identifiable:", both["identifiable"],
      f"(predicted sigma of mu {both['sigma']['mu']:.4f})")
try:
    fit_noise_spectra(f, sq_db=sq)
except ValueError as err:
    print("refused:", err)
```

```
mu    = 0.5474 +/- 0.0021
eta   = 0.6181 +/- 0.0035
kappa = 12.109 +/- 0.073 MHz
chi^2 / dof = 35.1 / 47
squeezed trace only identifiable: False
both traces identifiable: True (predicted sigma of mu 0.0021)
refused: one quadrature is a single Lorentzian -- two shape numbers for three unknowns (mu, eta, kappa). Supply both quadratures, or pass the independently measured kappa_hz
```

One quadrature alone is a single bell-shaped (Lorentzian) curve: two shape numbers
(depth and width) for three unknowns (`mu`, `eta`, `kappa`). The fit
therefore needs both traces, or a separately measured `kappa_hz`.
`plan_noise_measurement` gives the same error bars the fit would report
(the same `(J^T W J)^-1` matrix, the standard least-squares error
formula) before any data exist, and
`design_noise_frequencies` picks the most informative analysis
frequencies from a candidate list.

### 8. Calibrating g0 from the threshold power

```python
import os, tempfile
import numpy as np
from sqzcomb import (ring_from_threshold, threshold_power,
                     save_spectra_csv, load_spectra_csv)

# Illustrative measurement: threshold 0.15 mW +/- 0.005 mW,
# linewidth 100 MHz +/- 2 MHz.
ring, sigma_g0 = ring_from_threshold(
    kappa_hz=100e6, eta_esc=0.8, threshold_w=0.15e-3, lambda_pump_m=1.55e-6,
    reference="illustrative threshold sweep",
    sigma_kappa_hz=2e6, sigma_threshold_w=0.005e-3)
print(f"g0/2pi = {ring.g0_hz:.3f} +/- {sigma_g0:.3f} Hz")
print(f"threshold_power of the new ring: {threshold_power(ring) * 1e3:.6f} mW")
print(ring.reference)

# Measured spectra travel in a plain CSV file.
path = os.path.join(tempfile.mkdtemp(), "spectra.csv")
f = np.array([1e6, 2e6, 5e6])
save_spectra_csv(path, f, [-3.1, -2.8, -1.9], [6.2, 5.9, 4.0], sigma_db=0.1)
f2, sq2, an2, sig2 = load_spectra_csv(path)
print(open(path).read().splitlines()[0], "| round trip exact:",
      np.array_equal(f, f2) and np.array_equal(an2, [6.2, 5.9, 4.0]))
```

```
g0/2pi = 8.388 +/- 0.437 Hz
threshold_power of the new ring: 0.150000 mW
g0 calibrated by sqzcomb.lab.ring_from_threshold from the measured threshold power 0.00015 W at alpha=1; illustrative threshold sweep
f_hz,sq_db,anti_db,sigma_db | round trip exact: True
```

The single-photon Kerr shift `g0` is hard to measure directly; the comb
threshold power is routine. `ring_from_threshold` inverts
`threshold_power` for `g0` and propagates the measurement errors
(`g0` scales as `kappa^2 / P`).

### 9. A soliton, and the marginal mode

```python
import numpy as np
from sqzcomb import (soliton_seed, newton_state, fluctuation_matrix,
                     output_quadrature_variance, squeezing_db)

F, alpha, d2 = 1.9, 3.0, -0.25                 # illustrative, anomalous d2
sol, info = newton_state(soliton_seed(128, F, alpha, (d2,)), F, alpha, (d2,))
print("converged:", info["residual"] < 1e-12,
      f"| peak |psi|^2 = {np.abs(sol).max() ** 2:.3f}")

M, modes = fluctuation_matrix(sol, alpha, (d2,))
print("largest growth rate below 1e-8 in size:",
      abs(np.linalg.eigvals(M).real.max()) < 1e-8)
i0 = int(np.where(modes == 0)[0][0])
try:
    output_quadrature_variance(M, 0.9, 0.5, i0, modes.size)
except ValueError as err:
    print("refused:", str(err)[:60], "...")
v = output_quadrature_variance(M, 0.9, 0.5, i0, modes.size,
                               allow_marginal=True)
print(f"with allow_marginal=True, pump line at phi = 0: {squeezing_db(v):+.3f} dB")
```

```
converged: True | peak |psi|^2 = 7.290
largest growth rate below 1e-8 in size: True
refused: drift matrix is marginally stable (an eigenvalue's real part ...
with allow_marginal=True, pump line at phi = 0: +2.342 dB
```

`newton_state` solves the stationary LLE and never returns a state
whose residual is above the tolerance. Around a soliton the noise
equations have an eigenvalue at zero: sliding the pulse around the ring
costs nothing. The spectra functions refuse such a marginal matrix
unless you pass `allow_marginal=True`; a truly unstable one is refused
either way.

### 10. Bringing a model from QuTiP

```python
import numpy as np
import qutip
from sqzcomb import drift_from_qutip, single_mode_parametric

a = qutip.destroy(12)
H = 0.5j * 0.4 * (a.dag() ** 2 - a ** 2)       # quadratic: accepted
M, _ = drift_from_qutip(H, [1.0])
print("same drift matrix as single_mode_parametric(0.4):",
      np.allclose(M, single_mode_parametric(0.4)))
try:
    drift_from_qutip(a.dag() ** 2 * a ** 2, [1.0])   # a Kerr term
except ValueError as err:
    print("refused:", err)
```

```
same drift matrix as single_mode_parametric(0.4): True
refused: H is not quadratic in the mode operators; refusing to linearize it silently
```

`drift_from_qutip` reads the coefficients of a quadratic Hamiltonian
(the energy operator that defines a quantum model; "quadratic" means
it has at most two mode operators per term)
and rebuilds it; anything it cannot rebuild (such as a Kerr term) is
refused rather than silently linearized.

### 11. Above threshold: bright states, and the exact answer

```python
import numpy as np
from sqzcomb import (kerr_parametric_states, kerr_parametric_master,
                     steady_state, master_moments, master_output_spectrum,
                     output_quadrature_variance)

mu, chi = 1.5, 0.1               # 50 % above threshold; Kerr shift per photon
for s in kerr_parametric_states(mu, delta=0.0, chi=chi):
    print(f"alpha = {s['alpha']:.3f}  photons = {s['n']:6.3f}  stable = {s['stable']}")

bright = kerr_parametric_states(mu, 0.0, chi)[1]
H, c_ops, a, dims = kerr_parametric_master(mu, 0.0, chi, cutoff=70)
rho = steady_state(H, c_ops, dims)            # exact, no linearization
print(f"exact mean photon number: {master_moments(rho, a)['n']:.3f}")
phi = np.angle(bright["alpha"]) + np.pi / 2   # at right angles to alpha
for om in (0.0, 2.0):
    lin = output_quadrature_variance(bright["drift"], 1.0, om, 0, 1, phi=phi)
    exact = master_output_spectrum(H, c_ops, a, 1.0, [om], phi)[0]
    print(f"omega = {om}: linearized {lin:.4f}, exact {exact:.4f} (vacuum 0.5)")
```

```
alpha = 0.000+0.000j  photons =  0.000  stable = False
alpha = 3.052+1.365j  photons = 11.180  stable = True
alpha = -3.052-1.365j  photons = 11.180  stable = True
exact mean photon number: 10.642
omega = 0.0: linearized 0.9000, exact 12.4750 (vacuum 0.5)
omega = 2.0: linearized 0.6176, exact 0.5078 (vacuum 0.5)
```

Above threshold a Kerr ring does not grow forever: each photon shifts
the resonance by `chi` (in units of `kappa/2`, `chi = 2 g0 / kappa`),
and the growth stops at two bright states `+alpha` and `-alpha` (the
two phases of a parametric oscillator). `kerr_parametric_states` finds
them by algebra and gives the drift matrix around each one, so all the
spectra functions work there. `kerr_parametric_master` and
`steady_state` solve the same model exactly in a photon-number basis,
with no linearization. The two agree better as `chi` gets smaller (the
tests check that the gap shrinks). At `omega = 0` they differ a lot:
now and then the ring jumps between `+alpha` and `-alpha`, which the
linearized theory leaves out. That adds noise in a narrow band around
`omega = 0` whose width equals the slowest decay rate of the master
equation (about `8e-4` here, in units of `kappa/2`). It shows up even
in this quadrature, at right angles to `alpha`, because the exact
lobes sit very slightly rotated from the classical `alpha`; in the
quadrature along `alpha` the extra noise is far larger. The Fock basis must be large enough:
`steady_state` refuses when the top levels hold more than `1e-6` of the
population.

### 12. A non-Gaussian state

```python
import numpy as np
from sqzcomb import (kerr_parametric_master, coherent_state, master_evolve,
                     wigner, wigner_negativity)

N, chi = 40, 1.0
H, _, a, dims = kerr_parametric_master(0.0, 0.0, chi, cutoff=N)  # Kerr only
rho0 = coherent_state(2.0, N)                 # a laser-like (Gaussian) state
rho_t = master_evolve(H, [], rho0, [np.pi / chi])[0]
x = np.linspace(-6, 6, 241)
print(f"negative volume, start:     {wigner_negativity(rho0, x, x):.4f}")
print(f"negative volume, t = pi/chi: {wigner_negativity(rho_t, x, x):.4f}")
print(f"Wigner function at the origin: {wigner(rho_t, [0.0], [0.0])[0, 0]:+.4f}")
```

```
negative volume, start:     0.0000
negative volume, t = pi/chi: 0.2932
Wigner function at the origin: +0.0001
```

`master_evolve` follows the master equation in time. Starting from a
coherent (laser-like, Gaussian) state, the Kerr effect alone turns it
into a superposition of two coherent states (a "cat state") at
`t = pi / chi`; its Wigner function has negative regions, which no
Gaussian state has. For this lossless case the test suite compares
every matrix element with the exact analytic answer. With the ring's
loss switched on (pass `c_ops` instead of `[]`), the cat only survives
if the Kerr shift per photon is much larger than the loss rate: the
negative volume at `t = pi / chi` is 0.29 without loss, and with loss
about 0 for `chi = 1`, 0.019 for `chi = 10` and 0.165 for `chi = 50`.
The steady states of the parametric oscillator were non-negative in
the six cases we tried.

### 13. Technical noise: a wandering resonance

```python
import numpy as np
from sqzcomb import (RingSpec, lle_evolve, fluctuation_matrix,
                     output_quadrature_variance, squeezing_db,
                     lle_mode_amplitudes, resonance_noise_drive,
                     classical_noise_variance, normalized_psd)

ring = RingSpec(kappa_hz=100e6, eta_esc=0.8, g0_hz=10.0,
                lambda_pump_m=1.55e-6, reference="illustrative values")
chi = 2 * ring.g0 / ring.kappa                  # Kerr shift per photon, units kappa/2
F, alpha = 0.9, 0.3
psi = lle_evolve(np.full(16, 0.05 + 0j), F=F, alpha=alpha, t_end=100.0)
M, modes = fluctuation_matrix(psi, alpha)
i0 = int(np.where(modes == 0)[0][0])
c = lle_mode_amplitudes(psi, modes, chi)        # mean field, photon units
b = resonance_noise_drive(c)                    # all lines shift together

# An illustrative resonance-frequency noise: flat 1 Hz^2/Hz at 5 MHz.
f = 5e6
om, S = normalized_psd(ring, f, 1.0, kind="frequency")
phis = np.linspace(0, np.pi, 181)
q = [output_quadrature_variance(M, 0.8, om, i0, modes.size, phi=p) for p in phis]
k = int(np.argmin(q))
extra = classical_noise_variance(M, 0.8, om, [b], [S], i0, phi=phis[k])
print(f"photons in the pumped line: {abs(c[i0]) ** 2:.3e}")
print(f"quantum only:        {squeezing_db(q[k]):.3f} dB")
print(f"with the added noise: {squeezing_db(q[k] + extra):.3f} dB")
```

```
photons in the pumped line: 3.494e+06
quantum only:        -4.191 dB
with the added noise: -3.159 dB
```

`lle_mode_amplitudes` converts the steady state (here a flat,
single-line state; a comb or soliton works the same way) to photon
units (the
classical noise competes with vacuum noise, one photon's worth, so the
photon number matters). `resonance_noise_drive` describes the whole
comb of resonances shifting together, `normalized_psd` converts a
laboratory PSD (here a made-up flat 1 Hz^2/Hz) to model units, and
`classical_noise_variance` gives the variance it adds at the output,
to be added to the quantum noise. `pump_noise_drive` and
`gain_noise_drive` do the same for pump-laser amplitude and phase
noise. The package carries the noise through the resonator exactly;
the PSD itself (thermorefractive, Raman, laser) is yours to supply,
with its source.

### 14. A pulsed pump

```python
import numpy as np
from sqzcomb import parametric_pulse_drift, temporal_mode_variance

# Pump pulse (in units of the cavity time 2/kappa) whose peak gain
# briefly exceeds threshold (|mu| = 1.6 > 1); the detector's LO pulse
# is a Gaussian centred slightly later.
mu = lambda t: 1.6 * np.exp(-((t - 5.0) / 1.5) ** 2)
lo = lambda t: np.exp(-((t - 7.0) / 1.5) ** 2)       # normalized for you
drift = parametric_pulse_drift(mu)
for phi in (np.pi / 2, 0.0):
    r = temporal_mode_variance(drift, (0.0, 14.0), lo, eta=0.85, phi=phi,
                               max_step=0.05)
    print(f"phi = {phi:.3f}: {r['squeezing_db']:+.2f} dB "
          f"(peak photons in the ring {r['max_photons']:.2f})")
```

```
phi = 1.571: -5.91 dB (peak photons in the ring 4.69)
phi = 0.000: +16.71 dB (peak photons in the ring 4.69)
```

With a pulsed pump the drift matrix changes in time, so there is no
steady state. `temporal_mode_variance` integrates the quantum noise
forward in time and returns the noise of the pulse-shaped mode that the
LO selects (vacuum = 0.5, as everywhere). For a short time the gain
here is above threshold (`|mu| = 1.6`), which is fine for this linear
model: the light is simply amplified. For a real Kerr ring the
linearization needs the photon number to stay small against `1/chi`,
which is why the peak photon number is reported. For a pulse train
that repeats once per round trip, see "Pulsed pumping" under [What is
in the package](#what-is-in-the-package).

### 15. Two rings with different line spacings

```python
import numpy as np
from sqzcomb import (lle_evolve, ring_line_frequencies,
                     vernier_molecule_fluctuation_matrix,
                     output_variance_ports, squeezing_db)

F, alpha, d2 = 0.9, 0.5, -0.15
psi = lle_evolve(np.full(64, 0.1 + 0j), F=F, alpha=alpha, dispersion=(d2,),
                 t_end=300.0)
fsr = 600.0                                  # main-ring FSR in units of kappa/2
aux = ring_line_frequencies(np.arange(-8, 9), offset=0.3,
                            fsr=fsr * 1.07)  # auxiliary ring: 7 % wider spacing
M, g, modes, pairs = vernier_molecule_fluctuation_matrix(
    psi, alpha, J=1.5, gamma_b=1.0, fsr=fsr, aux_frequencies=aux,
    dispersion=(d2,), modes=range(-3, 4))
for k, d in zip(pairs["partner"], pairs["aux_delta"]):
    print(f"main line {k:+d}: nearest auxiliary line detuned by {d:8.2f}")
m = modes.size
j = m + int(np.where(pairs["partner"] == 0)[0][0])   # aux partner of line 0
phis = np.linspace(0, np.pi, 181)
best = min(output_variance_ports(M, g, 1.0, j, 0.0, phi=p) for p in phis)
print(f"best squeezing leaving through line 0's partner: {squeezing_db(best):.3f} dB")
```

```
main line -3: nearest auxiliary line detuned by   125.70
main line -2: nearest auxiliary line detuned by    83.70
main line -1: nearest auxiliary line detuned by    41.70
main line +0: nearest auxiliary line detuned by    -0.30
main line +1: nearest auxiliary line detuned by   -42.30
main line +2: nearest auxiliary line detuned by   -84.30
main line +3: nearest auxiliary line detuned by  -126.30
best squeezing leaving through line 0's partner: -1.990 dB
```

When the second ring's line spacing differs from the main ring's
(here by 7 %), its lines slide past the comb like a Vernier scale: only
one main line (here line 0) has a partner close to resonance, and the
rest see theirs far off. `vernier_molecule_fluctuation_matrix` pairs
each auxiliary line with its nearest main line and refuses when a line
sits too close to midway for that pairing to hold. `ring_line_frequencies`
builds the line frequencies from an offset, a spacing and a dispersion
(convert Hz with `normalized_frequency`). Mind the sign: its
`dispersion` uses the physical sign (`D2 > 0` for anomalous
dispersion, the lines spreading apart), while the main ring's
`dispersion` argument uses the LLE sign (`d2 < 0` anomalous, the value
`normalized_dispersion` returns).

### 16. Inference and fitting beyond the simplest model

```python
import numpy as np
from sqzcomb import infer_source, fit_spectra_model, parametric_spectra_model

# The measured pair of example 6, with and without a known purity.
pure = infer_source(-1.71, 5.54)
impure = infer_source(-1.71, 5.54, purity=0.8)
print(f"pure source assumed: eta = {pure['eta']:.3f}, source {pure['sq_db_source']:.2f} dB")
print(f"purity 0.8:          eta = {impure['eta']:.3f}, source {impure['sq_db_source']:.2f} dB")

# A detuned squeezer seen at two fixed LO angles (synthetic, noise-free).
names, model = parametric_spectra_model(detuned=True)
true = dict(mu=0.6, eta=0.7, kappa_hz=80e6, delta=0.3, phi_sq=1.4)
f = np.linspace(2e6, 300e6, 60)
data = {k: 10 * np.log10(v / 0.5) for k, v in model(true, f).items()}
bounds = dict(mu=(0, 0.99), eta=(0, 1), kappa_hz=(1e6, 1e9),
              delta=(-2, 2), phi_sq=(0, np.pi))
start = dict(mu=0.5, eta=0.5, kappa_hz=60e6, phi_sq=1.5)
try:
    fit_spectra_model(f, data, model, names, p0=dict(start, delta=0.1),
                      bounds=bounds)
except ValueError as err:
    print("refused:", str(err)[:95], "...")
fit = fit_spectra_model(f, data, model, names, p0=start, bounds=bounds,
                        fixed=dict(delta=0.3))
print(f"with delta known: mu = {fit.params['mu']:.4f}, eta = {fit.params['eta']:.4f}, "
      f"kappa = {fit.params['kappa_hz'] / 1e6:.2f} MHz")
```

```
pure source assumed: eta = 0.372, source -8.99 dB
purity 0.8:          eta = 0.415, source -6.64 dB
refused: the data are fitted equally well by clearly different parameter sets: {mu=0.6, eta=0.7, kappa_h ...
with delta known: mu = 0.6000, eta = 0.7000, kappa = 80.00 MHz
```

`infer_source` can now take a known purity and a known phase jitter.
If the source is in fact impure, assuming it pure makes it look better
than it is (here -8.99 dB instead of -6.64 dB). One measured pair
cannot tell you the purity; it has to come from elsewhere. When a
known purity allows two answers, the function refuses unless you say
which one (`branch`). `fit_spectra_model` fits any spectrum model; with
a detuning, two fixed-angle traces are fitted exactly by more than one
parameter set (one is the mirror image with the opposite detuning, and
there is another with a different gain and efficiency), and the fit
says so rather than returning one of them. Knowing the detuning
removes the ambiguity.

### 17. g0 from the nonlinear refractive index

```python
from sqzcomb import kerr_shift_from_n2

# Illustrative inputs (not material data): n2 = 2.4e-19 m^2/W, n0 = 2.0,
# V_eff = 6e-16 m^3, 1550 nm. Use numbers from your own source.
print(f"g0/2pi = {kerr_shift_from_n2(2.4e-19, 2.0, 6e-16, 1.55e-6):.3f} Hz")
```

```
g0/2pi = 0.743 Hz
```

`kerr_shift_from_n2` evaluates `g0 = hbar omega0^2 c n2 / (n0^2
V_eff)`. The package still ships no material values: `n2` depends on
the material, how it was made and the wavelength, and `V_eff` on the
geometry and on how the mode volume is defined. A `g0` measured from
the comb threshold (example 8) is usually the better number.

### 18. The best quadrature angle, exactly

```python
import numpy as np
from sqzcomb import (single_mode_parametric, optimal_quadrature,
                     output_quadrature_variance, squeezing_db)

# A detuned squeezer (illustrative: gain 0.6, detuning 0.3, eta = 0.8).
M = single_mode_parametric(0.6, delta=0.3)
r = optimal_quadrature(M, 0.8, [0.0, 0.5, 2.0], 0, 1)
for w, db, ph in zip([0.0, 0.5, 2.0], r["squeezing_db"], r["phi_min"]):
    print(f"omega = {w}: best {db:.4f} dB at phi = {ph:.4f} rad")

# The same number from a scan of 181 angles is slightly worse.
scan = min(output_quadrature_variance(M, 0.8, 0.5, 0, 1, phi=p)
           for p in np.linspace(0.0, np.pi, 181))
print(f"181-angle scan at omega = 0.5: {squeezing_db(scan):.4f} dB")
```

```
omega = 0.0: best -5.8030 dB at phi = 1.7915 rad
omega = 0.5: best -4.9141 dB at phi = 1.7588 rad
omega = 2.0: best -1.5193 dB at phi = 1.6275 rad
181-angle scan at omega = 0.5: -4.9130 dB
```

Every quadrature variance depends on the angle as `c0 + c1 cos 2phi +
c2 sin 2phi`, so three evaluations give the smallest and largest
value and their angles exactly; a scan over angles can only come close
(here 0.001 dB short). With a detuning the best angle changes with the
analysis frequency, so a fixed local-oscillator angle cannot be
optimal everywhere. `quadrature_extremes` does the same for any
function of the angle, for example
`lambda p: output_variance_ports(M, g, eta, 1, w, phi=p)` for a
photonic molecule; it checks the form with a fourth evaluation and
refuses a function that does not have it (such as one returning dB).

### 19. The detection budget, and raw analyzer traces

```python
import numpy as np
from sqzcomb import (detected_squeezing_db, dark_from_clearance_db,
                     required_efficiency_db, max_phase_noise,
                     shot_noise_normalize)

# Dark noise 15 dB below shot noise, 70 % efficiency, a -10 dB source.
dark = dark_from_clearance_db(15.0)
print(f"against pure vacuum:        {detected_squeezing_db(0.05, 0.7, dark):.3f} dB")
print(f"against a shot-noise trace: "
      f"{detected_squeezing_db(0.05, 0.7, dark, dark_in_reference=True):.3f} dB")

# The whole budget at once: -10 dB / +11 dB source, 20 mrad jitter,
# the dark noise above. Which efficiency still gives -6 dB?
print(f"efficiency needed for -6 dB: "
      f"{required_efficiency_db(-6.0, -10.0, dark, anti_db=11.0, theta_rms=0.02):.4f}")
v_sq, v_an = 0.5 * 10 ** -1.0, 0.5 * 10 ** 1.1
print(f"jitter allowed for -6 dB at 90 % efficiency: "
      f"{max_phase_noise(0.5 * 10 ** -0.6, v_sq, v_an, 0.9, dark) * 1e3:.1f} mrad")

# Raw analyzer traces (made up for the example, in dBm): three sweeps
# of the squeezed trace, the shot-noise trace and the dark trace at
# two analysis frequencies.
sq = [[-84.02, -83.10], [-83.95, -83.21], [-84.10, -83.15]]
shot = [[-80.00, -80.02], [-79.97, -80.05], [-80.03, -79.98]]
dark_tr = [[-95.1, -95.0], [-94.9, -95.2], [-95.0, -95.1]]
out = shot_noise_normalize(sq, shot, dark_tr)
print("dB re shot noise:", np.round(out["db"], 3),
      "+/-", np.round(out["sigma_db"], 3))
print("dark noise, vacuum units:", np.round(out["dark_variance"], 5))
```

```
against pure vacuum:        -3.962 dB
against a shot-noise trace: -4.097 dB
efficiency needed for -6 dB: 0.8720
jitter allowed for -6 dB at 90 % efficiency: 51.4 mrad
dB re shot noise: [-4.245 -3.286] +/- [0.05 0.04]
dark noise, vacuum units: [0.01633 0.01601]
```

A spectrum analyzer compares the squeezed trace with a shot-noise
trace taken through the same receiver, so both carry the dark noise.
Their ratio is the same as an extra loss of efficiency `(1/2) / (1/2 +
V_dark)` (the equivalence of electronic noise and loss shown by J.
Appel et al., Phys. Rev. A 75, 035802 (2007)); `dark_in_reference=True`
gives that number, and the default keeps the convention of earlier
versions. `required_efficiency` and `max_phase_noise` now take loss,
jitter and dark noise together and invert the budget exactly.
`shot_noise_normalize` turns raw traces (one row per repeated sweep)
into dB relative to shot noise: sweeps are averaged in linear power,
the dark trace is subtracted if given, and the error bar is the
standard error from the sweep-to-sweep scatter. The traces here are
made-up numbers for the example. It does not correct for an analyzer
that averaged in dB (log) mode. The dark variance it reports (0.01633
at the first frequency, where the shot and dark traces are 14.9993 dB
apart) is what `dark_from_clearance_db(14.9993, trace_gap=True)` gives;
without `trace_gap=True` the same 14.9993 dB would be read from the
pure shot-noise level and give 0.01581.

### 20. One local-oscillator angle for a whole band

```python
import numpy as np
from sqzcomb import single_mode_parametric, optimal_band_quadrature

# The detuned squeezer of example 18 (gain 0.6, detuning 0.3, eta = 0.8),
# recorded from omega = 0 to 2 with one fixed local-oscillator angle.
M = single_mode_parametric(0.6, delta=0.3)
band = np.linspace(0.0, 2.0, 21)
r = optimal_band_quadrature(M, 0.8, band, 0, 1)
print(f"best fixed angle: {r['phi']:.4f} rad "
      f"(band average {r['squeezing_db_band']:.3f} dB)")
for k in (0, 10, 20):
    print(f"omega = {band[k]:.1f}: {r['squeezing_db'][k]:.3f} dB at that angle, "
          f"{r['penalty_db'][k]:.3f} dB short of the best angle there "
          f"({r['phi_opt'][k]:.4f} rad)")
```

```
best fixed angle: 1.7506 rad (band average -3.253 dB)
omega = 0.0: -5.532 dB at that angle, 0.271 dB short of the best angle there (1.7915 rad)
omega = 1.0: -3.330 dB at that angle, 0.053 dB short of the best angle there (1.7000 rad)
omega = 2.0: -1.449 dB at that angle, 0.070 dB short of the best angle there (1.6275 rad)
```

In a measurement the local oscillator usually stays at one angle while
the analyzer sweeps many frequencies. With a detuning the best angle
moves with frequency (example 18), so no single angle is best
everywhere. `optimal_band_quadrature` finds the one angle that gives
the lowest average noise over the band, exactly and without a scan:
every variance has the form `c0 + c1 cos 2phi + c2 sin 2phi`, and so
does any weighted average of them. It also reports, at each frequency,
how much the fixed angle loses against the best angle there. Pass
`weights` (for example the width of each frequency bin) to weigh the
band the way your detector does. The average is taken over the
variance, which is the noise power a detector integrating over the
band would see.

### 21. An uncertain shot-noise reference

```python
import numpy as np
from sqzcomb import (noise_spectrum_db, fit_noise_spectra,
                     plan_noise_measurement, infer_source)

# The synthetic traces of example 7, divided by a shot-noise reference
# that is 0.15 dB too low (illustrative): both traces move up together.
mu, eta, kappa_hz = 0.55, 0.62, 12e6
f = np.linspace(0.5e6, 30e6, 25)
rng = np.random.default_rng(1)
sq = noise_spectrum_db(mu, eta, kappa_hz, f, "squeezed") + rng.normal(0, 0.05, f.size)
an = noise_spectrum_db(mu, eta, kappa_hz, f, "anti") + rng.normal(0, 0.05, f.size)
sq, an = sq + 0.15, an + 0.15

plain = fit_noise_spectra(f, sq_db=sq, anti_db=an, sigma_db=0.05)
cal = fit_noise_spectra(f, sq_db=sq, anti_db=an, sigma_db=0.05,
                        calib_sigma_db=0.2)          # reference known to 0.2 dB
print(f"offset ignored:   eta = {plain.eta:.4f} +/- {plain.sigma['eta']:.4f}, "
      f"chi^2 / dof = {plain.chi2:.1f} / {plain.chi2_dof}")
print(f"offset allowed:   eta = {cal.eta:.4f} +/- {cal.sigma['eta']:.4f}, "
      f"chi^2 / dof = {cal.chi2:.1f} / {cal.chi2_dof}")
print(f"fitted reference offset: {cal.calib_offset_db:+.3f} "
      f"+/- {cal.sigma['calib_offset_db']:.3f} dB")

# Before measuring: the error bar the 0.2 dB calibration leaves on eta.
for c in (None, 0.2):
    p = plan_noise_measurement(mu, eta, kappa_hz, f, sigma_db=0.05,
                               calib_sigma_db=c)
    print(f"planned sigma of eta, calibration error {c}: {p['sigma']['eta']:.4f}")

# One measured pair (example 6): the same calibration error moves both
# levels together, which matters more than two independent errors.
for c in (None, 0.1):
    out = infer_source(-1.71, 5.54, sigma_db=0.1, calib_sigma_db=c)
    print(f"source {out['sq_db_source']:.2f} +/- {out['sq_db_source_sigma']:.2f} dB, "
          f"eta {out['eta']:.4f} +/- {out['eta_sigma']:.4f} (calibration error {c})")
```

```
offset ignored:   eta = 0.5821 +/- 0.0036, chi^2 / dof = 339.5 / 47
offset allowed:   eta = 0.6173 +/- 0.0040, chi^2 / dof = 35.5 / 47
fitted reference offset: +0.147 +/- 0.008 dB
planned sigma of eta, calibration error None: 0.0035
planned sigma of eta, calibration error 0.2: 0.0040
source -8.99 +/- 0.25 dB, eta 0.3724 +/- 0.0204 (calibration error None)
source -8.99 +/- 0.43 dB, eta 0.3724 +/- 0.0300 (calibration error 0.1)
```

Every point of a squeezing trace is divided by the same shot-noise
reference. If that reference is off, every point of both traces moves
by the same number of dB. More frequencies do not average this error
away, and a fit that treats all errors as independent can be wrong by
many of its own error bars without saying so (here `eta` comes out
0.5821 +/- 0.0036 against a true 0.62; only the large chi-square hints
at it). With `calib_sigma_db`, `fit_noise_spectra` fits a common
offset of the reference along with the parameters, constrained by the
stated calibration error. This is the same as a fit with the full
correlated error matrix: `sigma_db^2` on the diagonal, plus
`calib_sigma_db^2` in every entry. The error bars are the ones that
matrix gives.
Here both traces approach shot noise at high frequency, so the data
pin the offset down well. `plan_noise_measurement` gives the same error
bars before measuring, and `infer_source` propagates the calibration
error into the inferred source, taking into account that it moves the
squeezed and antisqueezed levels together. The offset of 0.15 dB is
made up for the example.

## What is in the package

**Classical field**

- `lle_evolve(psi0, F, alpha, dispersion, t_end, dt, d1, t0)` -- time-steps
  the LLE (split-step method; the Kerr half-steps and the linear step
  are each solved exactly). The pump `F` can be a number, a profile
  around the ring, or a function of time (see "Pulsed pumping" below).
- `homogeneous_steady_states(F, alpha)` -- the flat-state intensities,
  roots of `rho (1 + (alpha - rho)^2) = F^2`.
- `newton_state`, `soliton_seed`, `continuation` -- solitons and
  soliton crystals by Newton's method, a pulse-shaped (sech) starting guess, and a
  sweep over detuning; unconverged states are never returned.

**Quantum noise and output spectra**

- `fluctuation_matrix(psi_s, alpha, dispersion, modes, d1)` -- the
  drift matrix around a steady state.
- `single_mode_parametric(mu, delta)` -- the drift matrix of one
  below-threshold parametric mode (the textbook squeezer).
- `output_quadrature_variance` -- the noise of one output quadrature
  (or the joint quadrature of two lines) at frequency `omega`, with
  optional thermal baths; `squeezing_db` converts to dB.
- `optimal_quadrature`, `quadrature_extremes` -- the most squeezed and
  most antisqueezed quadrature and their angles, exactly, at one
  frequency or over a spectrum (example 18).
- `optimal_band_quadrature` (new in 0.15): the one fixed angle that
  gives the lowest (optionally weighted) average noise over a band of
  frequencies, and what it costs at each frequency (example 20).
- `output_covariance_xxpp` -- the full output covariance matrix;
  `output_entanglement`, `output_entanglement_spectrum` -- `E_N`
  between two output lines, optionally after detection loss.

**Photonic molecules**

- `photonic_molecule(mu, J, delta_a, delta_b, gamma)` -- two coupled
  rings; `molecule_threshold(J, gamma)` -- its exact threshold on
  resonance; `output_variance_ports` -- the detected noise when rings
  have different decay rates and the output port sits on one (or
  several) of them, with optional thermal baths and `allow_marginal`
  as for a single ring.
- `molecule_fluctuation_matrix` -- the multimode version: every
  retained comb line coupled to a matching line of an auxiliary ring.
- `vernier_molecule_fluctuation_matrix`, `ring_line_frequencies` --
  the same when the two rings' line spacings differ: each auxiliary
  line is paired with its nearest main line.

**Above threshold and non-Gaussian states**

- `kerr_parametric_states`, `kerr_parametric_drift` -- the bright
  steady states of a Kerr parametric oscillator and the drift matrix
  around each.
- `kerr_parametric_master`, `fock_operators`, `liouvillian`,
  `steady_state`, `master_evolve`, `master_moments`,
  `master_output_spectrum` -- the exact master equation of one or two
  modes in a photon-number basis: steady state, time evolution,
  moments, and the output noise spectrum.
- `fock_state`, `coherent_state`, `wigner`, `wigner_negativity` --
  states to start from, the Wigner function, and its negative volume.

**Technical noise**

- `classical_noise_variance`, `classical_noise_variance_ports` -- the
  variance a classical noise adds at the output (one port, or the
  molecule's ports).
- `resonance_noise_drive`, `pump_noise_drive`, `gain_noise_drive` --
  how a wandering resonance, pump amplitude or phase noise (on one
  line, or on every line a pulsed pump feeds), or noise on the
  parametric gain enters; `lle_mode_amplitudes` -- the comb's
  mean field in photon units; `normalized_psd` -- laboratory PSD to
  model units.

**Pulsed pumping**

- `covariance_evolution`, `temporal_mode_variance`,
  `parametric_pulse_drift` -- the quantum noise followed in time for a
  pump that changes in time, and the noise of a pulse-shaped output
  mode.
- A pulse train repeating once per round trip enters the LLE as a pump
  profile `F(theta)`; `d1` is the drift when the pulse repetition rate
  `f_rep` differs from the ring's free spectral range `f_FSR`,
  `d1 = 4 pi (f_rep - f_FSR) / kappa`. `newton_state` and
  `fluctuation_matrix` accept both, so solitons driven by pulses and
  their noise spectra work as for a steady pump.

**Gaussian states and entanglement**

- `intracavity_covariance`, `covariance_xxpp`,
  `symplectic_eigenvalues`, `principal_quadratures` -- the state
  inside the ring, its covariance matrix, how mixed it is, and its
  most-squeezed collective quadratures.
- `drift_from_qutip` -- a drift matrix from a quadratic QuTiP
  Hamiltonian.
- `two_mode_reduction`, `ppt_symplectic_eigenvalue`,
  `logarithmic_negativity`, `duan_epr_sum`, `entanglement_report` --
  two-mode entanglement measures.
- `thermal_occupation(frequency_hz, temperature_k)` -- the
  thermal (Bose-Einstein) photon number, from the exact SI values of `h` and
  `k_B`.

**Imperfect detection**

- `detected_variance`, `detected_squeezing_db`,
  `dark_from_clearance_db`, `lossy_channel_xxpp` -- loss and
  electronic noise, on single variances or on covariance matrices
  (`trace_gap=True` for a clearance read as the gap between the
  shot-noise and dark traces);
  `dark_in_reference=True` and `dark_equivalent_efficiency` for a
  result quoted against a shot-noise trace that carries the dark noise.
- `required_efficiency`, `required_efficiency_db` -- the efficiency a
  target needs, optionally with dark noise and phase jitter too.
- `phase_noise_variance`, `phase_noise_squeezing_db`,
  `max_phase_noise` -- LO phase jitter, and the jitter a target
  allows (optionally after loss and dark noise).

**Lab units, fitting and planning**

- `RingSpec`, `normalized_pump`, `threshold_power`,
  `normalized_detuning`, `normalized_dispersion`,
  `normalized_frequency`, `physical_frequency`, `intracavity_photons`
  -- conversions between lab and model units.
- `fit_noise_spectra`, `NoiseFit`, `noise_spectrum_db` -- fit
  `(mu, eta, kappa)` and optionally a dark-noise floor to measured
  spectra, and the model curve itself. With `calib_sigma_db` (new in
  0.15) the fit also carries the uncertainty of the shot-noise
  reference, an error shared by every point (example 21).
- `fit_spectra_model`, `ModelFit`, `parametric_spectra_model`,
  `molecule_spectra_model`, `jitter_average` -- fit any spectrum model
  (ready-made: the single-mode model with optional detuning, LO jitter
  and dark floor; the two-ring molecule, optionally detuned), with
  checks that the data determine the parameters.
- `infer_source` -- the source behind a measured squeezed /
  antisqueezed pair, optionally with a known purity and phase jitter,
  and (new in 0.15, `calib_sigma_db`) a shot-noise calibration error
  that moves both levels together.
- `kerr_shift_from_n2` -- `g0` from `n2`, the refractive index and the
  mode volume that you supply.
- `plan_noise_measurement`, `design_noise_frequencies`,
  `ring_from_threshold`, `save_spectra_csv`, `load_spectra_csv` --
  measurement planning (with `calib_sigma_db` since 0.15), `g0`
  calibration, and a checked CSV format for spectra.
- `shot_noise_normalize` -- raw analyzer traces in dBm (repeated
  sweeps allowed) to dB relative to shot noise, with optional dark
  subtraction and standard errors.

Each function's docstring (`help(sqzcomb.output_quadrature_variance)`,
for example) gives its inputs, units and conventions.

## When it refuses, and why

`sqzcomb` raises an error instead of guessing when:

- a drift matrix is unstable (above threshold), where a linearized
  noise level does not exist -- in the spectra, the port spectra, the
  output covariance and the intracavity covariance;
- a drift matrix is only marginally stable (an eigenvalue with zero
  real part, as around a soliton), unless `allow_marginal=True` says
  that direction is understood;
- a thermal occupation or dark noise is negative, a decay rate
  (`gamma`, `gamma_b`) is not positive, a detection efficiency is
  outside (0, 1], a port's escape fraction `eta` is outside [0, 1]
  (in every spectra function; before 0.14 the single-ring spectra
  returned NaN, and `output_entanglement` 0), a variance or a dark
  clearance is negative, or `thermal_occupation` gets a non-positive
  frequency or a negative temperature;
- a mode index is not an integer in `[0, n_modes)` (a negative index
  used to read the wrong quadrature), `n_modes` does not match the
  drift matrix, or a joint quadrature names the same mode twice;
- a detector is placed on a line that never enters its bus
  (`output_variance_ports`);
- a matrix is not a valid covariance matrix (odd size, not symmetric,
  not positive semidefinite, not real, invalid for the entanglement
  formulas);
- a QuTiP Hamiltonian is not quadratic, or a mode has fewer than 3
  Fock (photon-number) levels;
- Newton's method does not reach its tolerance (also inside
  `continuation`, naming the failing detuning), or a soliton seed is
  asked for with normal dispersion (`d2 >= 0`);
- a `RingSpec` has no `reference`, a non-positive linewidth, `g0` or
  wavelength, or an escape efficiency outside (0, 1]; a pump power or
  intensity is negative;
- a loss or jitter budget is asked for a target the source cannot
  reach, or one that loss or jitter could never produce;
- `fit_noise_spectra` gets only one trace and no known `kappa_hz`,
  an antisqueezed trace whose mean is more than 0.25 dB below shot
  noise, fewer than 5 frequencies, or non-positive `sigma_db`; it also
  stops if the fit fails or the data do not determine the parameters;
- `infer_source` gets a non-finite value, a non-positive `sigma_db`,
  an antisqueezed level at or below shot noise, a
  squeezed level at or above it, a squeezed-times-antisqueezed product
  below the vacuum value (impossible after loss), or a pair implying an
  efficiency outside (0, 1];
- `plan_noise_measurement` and `design_noise_frequencies` get `mu`
  outside (0, 1), `eta` outside (0, 1], a non-positive `kappa_hz`,
  frequency or `sigma_db`, unknown quadratures, a bad `n_pick`, or
  (design only) candidate frequencies that cannot determine the
  parameters at all;
- a spectra CSV file is empty, has the wrong header, no data rows,
  rows of the wrong length, a non-numeric value or a non-positive
  frequency, or the columns to be saved differ in length;
- `steady_state` or `master_evolve` (given `dims`) find more than
  `trunc_tol` (default `1e-6`) of the population in the top levels of
  the photon-number basis; `master_output_spectrum` is given a thermal
  bath; `kerr_parametric_states` is given `chi = 0` (nothing would stop
  the growth);
- `infer_source` is given a purity outside (0, 1], a negative jitter, a
  jitter too large for the measured pair, a purity that no efficiency
  can match, or a purity that allows two answers without `branch`;
- `fit_spectra_model` finds that the data do not determine the
  parameters (it names the combination), or that clearly different
  parameter sets fit equally well;
- `vernier_molecule_fluctuation_matrix` finds an auxiliary line too
  close to midway between two main lines for the pairing to hold;
- `newton_state` does not converge, for example because a pulse drift
  `d1` is too large for the soliton to lock, so that no steady state
  exists (a refusal can also just mean the starting guess was too far
  off: let `lle_evolve` settle the state first);
- a pump profile has the wrong length or is not finite; a PSD is
  negative; a mode function vanishes on the time window;
- `quadrature_extremes` is given a function of the angle that is not
  of the form `c0 + c1 cos 2phi + c2 sin 2phi` (for example one that
  returns dB);
- `required_efficiency` is given jitter without the antisqueezed
  level, or a target that even efficiency 1 cannot reach with the
  stated dark noise and jitter; `max_phase_noise` a target that the
  loss and dark noise alone already miss;
- `shot_noise_normalize` gets traces of different lengths, a
  non-finite value, or a dark trace that is not below the others;
  `load_spectra_csv` a non-finite value or a non-positive `sigma_db`;
- (new in 0.15) the modes read by `output_variance_ports` or
  `classical_noise_variance_ports` name the same mode twice (a repeated
  port alone is fine); `noise_spectrum_db` gets a quadrature name other
  than `"squeezed"` or `"anti"` (`"antisqueezed"` also works);
  `lossy_channel_xxpp` or `output_entanglement` gets a NaN efficiency,
  and the entanglement functions a covariance matrix with NaN or
  infinite entries; `soliton_seed` gets a detuning that is not
  positive or fewer than one pulse;
- (new in 0.15) `calib_sigma_db` is not finite and positive, or is
  given to `fit_noise_spectra` without `sigma_db`;
  `optimal_band_quadrature` gets an empty or non-finite frequency
  list, or weights that are negative, not finite, all zero or of the
  wrong length.

## How the results are checked

195 automated tests run on every push and pull request, on Python 3.9,
3.10, 3.11, 3.12, 3.13 and 3.14, and once more on Python 3.10 with the
oldest versions the package allows (NumPy 1.22.0, SciPy 1.10.0, QuTiP
4.7.0). The five QuTiP tests are skipped when QuTiP is not installed,
which is the case in the main matrix; they run in the oldest-versions
job. The numerical checks compare the package with a closed-form
result, an exact identity, or a second, independent calculation; none
compares against a number stored from an earlier run. The main ones,
with the tolerances the tests use:

**Output spectra**

- A passive cavity returns vacuum (0.5) at 10 random couplings,
  frequencies and phases, to 1e-12; with thermal baths at `n_bar`,
  exactly `(2 n_bar + 1)/2` on a grid of couplings, frequencies and
  phases, to 1e-12.
- The single parametric mode matches the closed form
  `0.5 (1 - 4 eta mu / ((1 + mu)^2 + omega^2))` to a relative 1e-10.
- At critical coupling and `mu = 0.9999` the detected squeezing lies
  between -3.011 and -2.99 dB; with full extraction at `mu = 0.99` it
  is below -20 dB.
- With thermal inputs the parametric spectrum is `(2 n_bar + 1)` times
  the vacuum one, to 1e-12; a hot-loss / cold-port case matches a
  hand-derived formula to 1e-12; the Bose function equals 1 exactly
  where `h f = k_B T ln 2` (1e-12).
- The LLE solver's free decay matches `exp(-t)` (relative 1e-3), and
  its flat state matches the cubic root (relative 1e-4).

**Best quadrature angle (new in 0.14)**

- `optimal_quadrature` equals the eigenvalues of the 2x2 output
  covariance block from `output_covariance_xxpp` (a separate code
  path) to 1e-12, and its angle equals the eigenvector's to 1e-9, for
  a detuned squeezer with a complex gain and for the pumped line of a
  comb; through the molecule's port with `J = 0` it equals the
  single-mode closed forms (relative 1e-12). An angle scan never goes
  below it, and misses it by no more than the grid allows
  (`2 r (dphi/2)^2`). For the twin beam, `E_N = -ln(2 V_min)` now
  holds to 1e-10 (the scan-based test needed 1e-6).

**New in 0.15**

- `optimal_band_quadrature` equals the smallest eigenvalue, and its
  angle the eigenvector, of the band average of the 2x2 output
  covariance blocks from `output_covariance_xxpp` (a separate code
  path), to 1e-12 and 1e-9, with equal and with random weights; no
  angle of a 181-angle scan does better on the band. With a real gain
  and no detuning it gives `pi/2` and the closed form at every
  frequency (1e-12), and with one frequency, or all weight on one, it
  gives that frequency's `optimal_quadrature` (1e-14).
- `optimal_quadrature` now builds the output matrix once per frequency;
  its results equal the per-angle calls of 0.14 exactly (tested with
  `==`).
- With `calib_sigma_db`, the error bars of `fit_noise_spectra` and
  `plan_noise_measurement` equal `(J^T C^-1 J)^-1` built in the test
  from the closed-form spectra and the explicit matrix
  `C = sigma_db^2 I + calib_sigma_db^2 (all ones)` (planner 1e-4, fit
  2 %); the fitted offset and chi-square equal the closed-form optimum
  offset and `r^T C^-1 r` (relative 1e-6). In 30 seeded trials with a
  common 0.2 dB offset, the plain fit scatters more than 3 times its
  error bar, and the calibrated fit between 0.6 and 1.6 times its own.
- `infer_source` with `calib_sigma_db` matches the closed-form
  gradients of `eta = -2 s a / (s + a)` and of the source level
  `10 log10(-s / a)` (relative 1e-6), and 1500 seeded trials with a
  common offset (within 10 %); with `sigma_db` alone it is unchanged.
- The four fixes under [Corrections](#corrections) each have a test
  that fails on 0.14.0.

**Photonic molecules**

- With `J = 0` the molecule reproduces the single-mode closed form to
  1e-12; undriven, it returns vacuum at 25 random settings of every
  parameter, to 1e-12.
- The undriven supermodes split by `2 J` (1e-12); the threshold
  formula agrees with the eigenvalues just below and just above it
  (0.1 % either side) for four `(J, gamma)` pairs.
- At `omega = 0` the output through the second ring equals a single
  mode with decay `1 + J^2/gamma` and efficiency
  `eta (J^2/gamma)/(1 + J^2/gamma)`, to 1e-12, for four parameter sets;
  the `J^2/gamma = 3` molecule gives between -6.03 and -6.01 dB.
- With thermal baths the port spectra of a passive molecule are
  exactly `(2 n_bar + 1)/2` at 12 random settings (1e-12); with
  `J = 0` a hot loss bath and a colder port give the same result as
  the single-ring path (1e-12), and a resonant molecule with every bath
  at `n_bar` gives `(2 n_bar + 1)` times its vacuum spectrum (relative
  1e-12). At threshold (`allow_marginal=True`) the port path equals the
  single-ring path (1e-13) and the closed form (relative 1e-12).
- The multimode molecule equals the two-ring model at one line
  (1e-14), reduces to the plain comb machinery at `J = 0` (1e-12), and
  satisfies the same `omega = 0` equivalence for the twin-beam
  quadrature of a line pair (1e-12).

**Gaussian states and entanglement**

- Intracavity covariances of vacuum and of the parametric mode match
  closed forms (1e-12 and 1e-10); in the tested molecule and
  lossy-channel cases every symplectic eigenvalue stays at or above
  `hbar/2`.
- For the two-mode squeezed vacuum (the textbook entangled pair, with
  squeezing strength `r`), `E_N = 2r` and the Duan sum
  `= hbar e^{-2r}` to a relative 1e-12; the principal variances are
  `(hbar/2) e^{-/+2r}` to 1e-12.
- The fast PPT formula equals an explicit partial transposition on 10
  random states (relative 1e-9); `E_N` is unchanged by local
  operations (relative 1e-10).
- Every quadrature read off `output_covariance_xxpp` equals the
  separate `output_quadrature_variance` result to 1e-13; for a twin
  beam, `E_N(omega) = -ln(2 V_min)` to 1e-6 (limited by the phase grid
  used to find `V_min`); lowering the detection efficiency step by
  step (0.9, 0.5, 0.2) never increases `E_N`.
- The QuTiP adapter reproduces the parametric mode and the two-ring
  molecule, including the detuning sign, to 1e-9.

**Detection and phase noise**

- Vacuum is unchanged by loss (1e-15); two losses equal one of the
  product efficiency (1e-15 scalar, 1e-13 matrix); the scalar and
  matrix forms agree (1e-13).
- With `dark_in_reference=True` the result equals plain loss with
  efficiency `eta (1/2)/(1/2 + V_dark)` and the literal ratio of the
  two traces (1e-15). `required_efficiency` inverts loss, jitter and
  dark noise together on a grid of 2 sources, 3 jitters, 3 dark levels,
  both conventions and 3 efficiencies (1e-12; grid points whose
  target is not below shot noise are skipped), and `max_phase_noise`
  inverts jitter after loss and dark noise (relative 1e-9); their
  defaults are unchanged bit for bit.
- The jitter formula equals direct numerical integration over the
  Gaussian distribution to 1e-10; zero, static and very large jitter
  give the known limiting values (to 1e-12 or better); `max_phase_noise` inverts it to 1e-12; loss
  and jitter give the same result in either order (1e-12).

**Solitons**

- From a flat seed, Newton lands on the cubic root (1e-10); a soliton
  converges below 1e-12 and moves by less than 5e-3 when handed to the
  independent time-stepper for 20 time units.
- The translation mode (sliding the pulse around the ring) is a zero
  mode of the drift matrix, `M z = 0`, to a relative 1e-9, and the largest eigenvalue real part is below 1e-8 in
  size.
- The two-pulse crystal equals the single soliton at four times the
  dispersion, resampled, to 1e-8 in amplitude.

**Lab units, fitting, planning, inference**

- The frequency conversion round-trips (relative 1e-15), and `F`
  scales as the square root of power (relative 1e-12); at `threshold_power`,
  `F^2 = 1 + (alpha - 1)^2` (relative 1e-12) and the cubic has the
  root 1 (1e-10).
- The fit model equals the closed form to 1e-10 dB; noise-free spectra
  give back `mu`, `eta` and `kappa` to 1e-6; a dark floor is recovered
  to 1e-4. With 0.1 dB noise, the scatter of 40 seeded fits lies
  between 0.4 and 2 times the reported error bars.
- The planner's error bars equal the fit's within 2 %, also at
  `eta = 1`; one quadrature alone gives a scaled singular value below
  1e-8 of the largest. The greedy frequency choice was never beaten
  by any of 30 random subsets of the same size.
- `shot_noise_normalize` on noise-free traces built from
  `detected_variance` gives back `detected_squeezing_db` with
  `dark_in_reference=True` (no subtraction) or with no dark noise
  (subtraction) to 1e-12 dB, and the dark variance to a relative
  1e-12, which also equals `dark_from_clearance_db(gap,
  trace_gap=True)` for the shot-to-dark trace gap (relative 1e-12);
  over 400 seeded repetitions of 8 noisy sweeps, the scatter of
  the result is between 0.85 and 1.15 times the reported error bar.
- `ring_from_threshold` round-trips through `threshold_power` to a
  relative 1e-12; its error bar matches finite differences within
  0.1 %. The CSV round trip is bit-exact.
- `infer_source` inverts `detected_variance` on a grid of 4 sources
  times 4 efficiencies (efficiency to 1e-9, squeezed variance to 1e-12);
  its error bars match the scatter of 600 seeded Monte Carlo trials
  (repeated simulated measurements with random noise) within 20 %.
  With a known purity and jitter it inverts loss plus jitter on a grid
  of 3 squeezing levels, 4 purities, 4 efficiencies and 3 jitters
  (efficiency to 1e-9, source to a relative 1e-9), and the grid
  contains two-answer cases, which it flags; the pure default is
  unchanged digit for digit; assuming purity where there is none gives
  a lower efficiency and a source more than 1 dB better than the truth.
- `fit_spectra_model` with the single-mode model equals
  `fit_noise_spectra` (relative 1e-6); noise-free data give back the
  jitter model's and the molecule's parameters to a relative 1e-6; the
  scatter of 30 seeded fits matches the reported jitter error bar
  within 35 %; one trace is refused as unidentifiable; the detuned
  pair is refused as ambiguous (its mirror image fits to 1e-10 dB) and
  fitted to 1e-6 once the detuning is fixed. The jitter law equals
  direct numerical averaging to a relative 1e-10.
- `kerr_shift_from_n2` gives the package's threshold power the
  textbook form `kappa^2 n0^2 V_eff / (8 eta omega0 c n2)` (relative
  1e-12).

**Above threshold and non-Gaussian states**

- With no Kerr term, the exact master-equation spectrum equals the
  linearized one at four quadrature angles, three frequencies, and
  also for a complex pump with detuning (1e-8); the intracavity
  covariance matches too (1e-8), and so does the closed form of the
  parametric mode (1e-9).
- The steady state and the Wigner function agree with QuTiP's
  independent solvers (1e-8 and 1e-10). The Wigner function of one
  photon is `-1/pi` at the origin (1e-14), integrates to 1 (1e-9), has
  the exact position marginal (1e-9) and negative volume
  `2 exp(-1/2) - 1` (relative 1e-4).
- Lossless Kerr evolution equals the analytic phase pattern element by
  element (1e-10) and makes a cat state with negative volume above 0.2.
- Bright states solve the classical equation (1e-9), their drift
  matrix equals a finite-difference Jacobian (1e-6), and above
  threshold the exact photon number is within 12 % of the bright
  state's; the exact spectrum approaches the linearized one as `chi`
  shrinks, and the extra low-frequency noise has the width of the
  slowest master-equation rate (2 %). A too-small photon-number basis
  is refused, also by `master_output_spectrum`.

**Technical noise**

- A passive mode with mean field `c` gets exactly
  `4 eta |c|^2 S / (1 + omega^2)` in one quadrature and nothing in the
  other (relative 1e-12). A time-domain calculation (stationary
  covariance of the resonator plus a coloured classical noise) equals
  the integral of the added spectrum (relative 1e-7), and a seeded
  simulation of 12 million time steps matches the spectrum within 4 %.
  A slow shift of the resonance or the gain moves the exact bright
  state of the parametric oscillator as the zero-frequency response
  predicts (relative 1e-4). The unit conversions keep the variance
  (Parseval, relative 1e-9).

**Pulsed pumping**

- A passive cavity returns exactly the vacuum value 0.5 for a pulsed
  LO mode (1e-12); with a steady pump the time-domain result converges
  to the steady covariance (1e-9) and a long pulse gives the steady
  spectrum weighted by the pulse's spectrum (relative 1e-8, also for
  a mode on two sidebands). A seeded simulation of 12 000 noisy
  trajectories with a pump pulse briefly above threshold matches the
  pulse's variance within four standard errors.
- A uniform pump profile equals a scalar pump (1e-13); with a uniform
  pump, `d1` only moves the pattern (1e-12); a moving pump profile seen
  from its own frame equals the `d1` equation (2e-5). That checks the
  change of frame behind `d1`; its conversion to `f_rep`, `f_FSR` and
  `kappa` is derived in the `lle` module's docstring. With a pulse-shaped pump, the soliton on the pulse
  peak is unstable along its translation direction, and the stable one
  sits beside the peak with every growth rate below -0.01 (so the
  spectra need no `allow_marginal`); its mirror image is also a steady
  state (1e-11).

**Two rings with different line spacings**

- With equal spacings the new builder equals the matched molecule
  (1e-12). For passive rings its eigenvalues equal those of the exact
  model with every line coupled to every line to within the neglected
  couplings (below `1e-2` here), and every port returns vacuum
  (1e-12).

## Corrections

**0.15.0 (this release)** fixes four inputs that gave a wrong answer
without an error. Each fix has a test that fails on 0.14.0:

- `output_variance_ports` and `classical_noise_variance_ports` counted
  a mode twice when it was named twice among the modes read
  (`mode_index=[1, 1]`, or `port_mode=[1, 1]` with the default
  `mode_index`), and so returned twice the variance: 3.8444 instead
  of 1.9222 for the twin beam of example 5 at `eta = 0.8`, `omega = 0`,
  `phi = 0.3`. This is now refused. A repeated port with distinct
  modes read still works and gives the same number as before.
- `noise_spectrum_db` returned the antisqueezed curve for every
  quadrature name other than `"squeezed"`. A typo such as
  `"Squeezed"` gave +8.141 dB at 1 MHz where the squeezed trace is
  -5.210 dB (`mu = 0.5`, `eta = 0.8`, `kappa_hz = 10e6`). Unknown
  names are now refused; `"anti"` and `"antisqueezed"` work as before.
- A NaN detection efficiency passed `lossy_channel_xxpp`, and the
  entanglement functions turned the resulting NaN matrix into
  `E_N = 0`, "not entangled" (the twin beam above has `E_N = 1.2417`,
  or 0.6887 at 70 % detection). Both are now refused.
- `soliton_seed` returned an all-NaN seed for a negative detuning
  (with only a NumPy warning) and raised `ZeroDivisionError` at zero
  detuning. It now refuses both with a message, and also refuses
  fewer than one pulse.

The `soliton` module docstring gave the seed's width as
`sqrt(|d2| / alpha)`; the code uses, correctly, `sqrt(|d2| / (2
alpha))`, the exact width of the sech solution of the conservative
part of the equation in this package's convention (derived in the
docstring). Only the docstring changed.

**0.14.0** fixed inputs that gave wrong answers without an error:

- `output_quadrature_variance`, `output_covariance_xxpp` and the
  output-entanglement functions accepted an escape fraction `eta`
  outside [0, 1]. The variance and the covariance came out NaN (with
  only a NumPy warning), and `output_entanglement_spectrum` then
  reported `E_N = 0`, "not entangled" (for example for the twin beam
  of example 5 at `eta = 1.5`). They now refuse.
- A negative mode index (`mode_index=-1`) was accepted by
  `output_quadrature_variance`, `output_variance_ports` and the
  technical-noise functions and read the last mode at the angle `-phi`
  instead of `phi` (for `photonic_molecule(0.5, 1.0, delta_a=0.4)`,
  `eta = 0.5`, `omega = 0`, `phi = 0.4`: 0.5098, which belongs to
  `-phi`, instead of 0.4347 for mode 1). Naming the same mode twice
  for a joint quadrature returned twice the single-mode variance. Both are now refused, as is an `n_modes`
  that does not match the drift matrix.
- `output_variance_ports` refused a marginally stable drift matrix
  (for example a soliton's) as "above threshold" and had no
  `allow_marginal`; it now uses the same stability test as the
  single-ring spectra, so a largest growth rate within 1e-6 of zero
  needs `allow_marginal=True` (before, any negative value, even
  `-4.9e-32` from rounding exactly at threshold, was accepted).
- `load_spectra_csv` accepted `nan` or `inf` in the dB columns and a
  zero or negative `sigma_db`; it now refuses them.

**0.13.0** fixed one bug: `output_quadrature_variance`
and `output_variance_ports` read the quadrature at angle `-phi`
instead of the documented `phi` (`X_phi = cos phi x + sin phi p`, the
convention of `output_covariance_xxpp`). The two readings agree at
`phi = 0` and `pi/2`, for the best angle found by scanning all angles,
and at every angle when the pump parameter is real and there is no
detuning, so examples 1 to 10 and `fit_noise_spectra` are unchanged.
Otherwise (a complex pump parameter or a detuning, read at some other
angle) a result at `phi` from 0.12.1 or earlier belongs to `-phi`.
Every Lugiato-Lefever comb state has a detuning, so comb results read
at a fixed angle other than 0 or `pi/2` did change. The exact
master equation of 0.13.0 confirmed the fix independently.

**0.12.1** fixed one bug:
`plan_noise_measurement` computed its sensitivities with a step that
went past `eta = 1`. At `eta = 1`, which the function accepts, the
model then returned NaN, the planner reported a measurable design as
"not identifiable", and `design_noise_frequencies` refused it. It now
steps backward at that edge. The fit itself was not affected.

That release also corrected the documentation. The old README's short
prediction example stopped with an "unstable" error (its 256-line
state is above threshold); it said the xxpp export uses `hbar = 2`,
which is true of `covariance_xxpp` but not of `output_covariance_xxpp`
(vacuum 0.5 I); it listed Python 3.9-3.14 in CI, but 3.10 was not
tested; and several "every" or "exact" statements were stronger than
the tests (see [How the results are checked](#how-the-results-are-checked)
for the real tolerances). The full history is in
[CHANGELOG.md](CHANGELOG.md).

## Limits

What 0.15.0 adds to the picture:

- **One angle for a band.** The 0.14 limit "the package does not yet
  choose the best single angle for a band" is closed for the single
  ring (example 20), and for the molecule's ports through
  `quadrature_extremes` (see the docstring of
  `optimal_band_quadrature`). The angle minimizes the average
  variance over the band. Other goals, such as the best worst case in
  dB, are not offered.
- **Calibration error.** Only an error of the shot-noise reference is
  modelled, as one offset in dB shared by every point of both traces
  (both traces divided by the same reference). An error that differs
  between the traces, or that changes with frequency, is not. Its
  size, `calib_sigma_db`, has to come from your own calibration
  record; the package cannot know it.

What 0.14.0 added:

- **Quadrature angle.** The best and worst angle are now found exactly
  (example 18), so results no longer carry the error of an angle scan.
  In an experiment the local oscillator usually sits at one angle for
  all frequencies; with a detuning that angle cannot be the best one
  everywhere, and 0.14 did not choose the best single angle for a
  band (0.15 does; see above).
- **Dark noise.** Two conventions are available (example 19); the
  default is still the one of earlier versions. Which one matches a
  published number depends on how that number was normalized, which
  the package cannot know.
- **Raw traces.** `shot_noise_normalize` assumes the traces were
  averaged in linear power, taken with the same analyzer settings, and
  independent from sweep to sweep; its error bar is first order in the
  scatter and needs at least two sweeps of every trace.

What 0.13.0 changed about the limits of earlier versions, and what is
left:

- **Above threshold and non-Gaussian states.** The Kerr parametric
  oscillator's bright states and the exact master equation (examples
  11 and 12) cover these for one or two modes. The exact model's cost
  grows as `cutoff^(2n)` for `n` modes, so a whole comb above
  threshold is still only available linearized around a stable steady
  state (such as a soliton); the linearized spectra still refuse an
  unstable drift matrix.
- **Technical noise.** Any classical noise with a known spectrum can
  now be carried through the resonator (example 13). The spectrum
  itself (thermorefractive, Raman, laser noise) is not computed from
  material physics: you supply it, with its source.
- **Pulsed pumping.** A pump that changes over many round trips
  (example 14) and a pulse train repeating once per round trip (pump
  profile and `d1`) are covered. A pump that changes during a single
  round trip in any other way is not.
- **Two rings with different line spacings** are covered (example 15)
  within the pairing approximation, which the builder checks.
- **Inference and fitting.** `infer_source` takes a known purity and
  jitter, and `fit_spectra_model` fits other models. One measured pair
  still cannot determine the purity, and a detuned squeezer's two
  fixed-angle traces cannot determine the detuning; the functions say
  so rather than guess.
- **Material constants.** Still none ship, on purpose: `n2`, loss and
  dispersion depend on the material, its fabrication and the
  wavelength, and a published value for one device is not a value for
  yours. `kerr_shift_from_n2` computes `g0` from numbers you supply;
  every `RingSpec` still requires a `reference`.

## Where it comes from

The methods were developed for the associated paper:

> T. M. Mahim, M. M. Rahman and A. S. M. Mohsin, "Overcoming the 3 dB
> squeezing extraction limit in silicon carbide microcombs with a
> photonic molecule," Optics Express 34(18), 34822-34834 (2026),
> https://doi.org/10.1364/OE.612248 (open access); code for the paper:
> https://github.com/Tanvir-Mahmud-Mahim/sic-molecule-squeezer

This package is the general-purpose engine; the paper repository
reproduces the specific published study. Other literature the methods
rely on is cited in the docstrings and in [CHANGELOG.md](CHANGELOG.md).

## Citing, support and license

If `sqzcomb` helps your work, please cite it with the concept DOI
[10.5281/zenodo.22015375](https://doi.org/10.5281/zenodo.22015375),
which always resolves to the latest release; every release is archived
on Zenodo. [CITATION.cff](CITATION.cff) has the details.

Written and maintained by Tanvir Mahmud Mahim (Department of
Electrical and Electronic Engineering, BRAC University), who reviews
every change and takes the final decision on scope and releases.
Design questions are discussed in the open in issues and pull
requests, and the standing rule of [CONTRIBUTING.md](CONTRIBUTING.md)
binds the maintainer exactly as it binds contributors: a change that
touches physics arrives with a test, and a constant arrives with its
source.

Support runs through the
[issue tracker](https://github.com/TaN-MM-Org/sqzcomb/issues). Usage
questions are welcome alongside bug reports; a docstring that left a
unit or a sign convention unclear is treated as a documentation bug,
not user error. While the version is below 1.0 the API may still move
between minor versions; such changes are called out in the release
notes. The normalized-unit conventions above are stable: any change
to them would be a breaking change named in the release notes, never
a quiet renormalization.

Licensed under Apache-2.0.
