Metadata-Version: 2.4
Name: cavsqueeze
Version: 1.17.0
Summary: Beyond-mean-field cavity-mediated spin squeezing for solid-state clock-transition ensembles
Author-email: Tanvir Mahmud Mahim <tanvir.mahim@bracu.ac.bd>
License: Apache-2.0
Project-URL: Homepage, https://github.com/TaN-MM-Org/cavsqueeze
Project-URL: Issues, https://github.com/TaN-MM-Org/cavsqueeze/issues
Project-URL: Paper-Companion-Repository, https://github.com/Tanvir-Mahmud-Mahim/yb-cawo4-cavity-squeezing
Keywords: spin squeezing,cavity QED,cumulant expansion,one-axis twisting,rare-earth ions,quantum metrology
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Requires-Dist: matplotlib>=3.7
Provides-Extra: exact
Requires-Dist: qutip>=5.0; extra == "exact"
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Requires-Dist: qutip>=5.0; extra == "test"
Dynamic: license-file

# cavsqueeze

[![Tests](https://github.com/TaN-MM-Org/cavsqueeze/actions/workflows/tests.yml/badge.svg)](https://github.com/TaN-MM-Org/cavsqueeze/actions/workflows/tests.yml)
[![PyPI](https://img.shields.io/pypi/v/cavsqueeze?label=PyPI&color=blue)](https://pypi.org/project/cavsqueeze/)
[![License](https://img.shields.io/badge/License-Apache_2.0-blue.svg)](LICENSE)
[![DOI](https://img.shields.io/badge/DOI-10.5281%2Fzenodo.22278034-blue)](https://doi.org/10.5281/zenodo.22278034)

`cavsqueeze` is a Python package for simulating and measuring **spin
squeezing** in a large group of spins (atoms, ions or defects in a
crystal) that all talk to one microwave or optical **resonator**
(a cavity). Squeezing is a way of preparing the spins so that they
measure a phase or a frequency more precisely than the same number of
independent spins could. The package was developed for
171Yb3+:CaWO4 (ytterbium ions in a calcium tungstate crystal) and
applies to any spin ensemble coupled to a resonator.

It answers these questions:

- How much squeezing does a given cavity and spin ensemble produce,
  and after how long, once the real imperfections of a solid are
  included: a spread of spin frequencies, decay into the cavity,
  thermal photons, loss of coherence, and uneven coupling?
- What does that squeezing buy a clock or a magnetometer?
- Given shot-by-shot data from an experiment, how much squeezing was
  actually there, and with what error bar?
- How many shots does a measurement need to reach a target error bar,
  before it is taken? And how large is the error bar when the
  detector's noise is not Gaussian?
- What coherence time `T2` and line width do my own echo and Ramsey
  measurements give, so that I can put them into the model? And do my
  measured data actually follow the laws the fits assume?
- When does the noise of the clock laser (the Dick effect) make
  squeezing useless?

The solver groups spins with the same frequency and coupling into
**classes**. Its memory and running time grow with the number of
classes, not with the number of spins: the state it stores has the
same size for 10^3 spins as for 10^15 (example 2 uses 10^15). The main
results are checked by automated tests against exact solutions, closed
formulas and an independent second solver (see
[How the results are checked](#how-the-results-are-checked)). When an
input cannot give a meaningful answer, the package stops with an error
that says why, rather than returning a number that looks fine but is
not.

## Contents

- [A short guide to the words used here](#a-short-guide-to-the-words-used-here)
- [Install, requirements and units](#install-requirements-and-units)
- [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 in earlier versions](#corrections-in-earlier-versions)
- [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

- **Spin, collective spin.** Each spin is a two-level system (up or
  down). Adding up all N spins gives one large **collective spin** `J`
  with components `J_x`, `J_y`, `J_z`. `J_z = (N_up - N_down)/2`
  is what a population measurement reads out.
- **Coherent spin state.** All spins point the same way (here along
  `+x`). Its measurement noise is the **standard quantum limit**
  (SQL): a phase uncertainty of `1/sqrt(N)` radians per shot. This
  noise, which comes from each spin landing randomly up or down when
  it is measured, is called **projection noise**.
- **Squeezing.** Correlations between the spins shrink the noise in
  one direction perpendicular to the mean spin, at the cost of more
  noise in the other. Two numbers measure it, both equal to 1 for a
  coherent spin state and smaller than 1 when squeezed:
  - `xi2_S` (**Kitagawa-Ueda parameter**): the smallest variance
    across the mean spin, divided by its coherent-state value `N/4`.
  - `xi2_R` (**Wineland parameter**): `xi2_S` divided by the contrast
    squared. It is the one that sets phase precision. The gain over
    the SQL in decibels is `-10 log10(xi2_R)`.
- **Contrast.** The length of the mean spin as a fraction of its
  largest possible value, `2|J|/N`. Decay and dephasing shrink it.
- **One-axis twisting (OAT).** The interaction that creates the
  squeezing here: each spin's precession speed depends on `J_z`, which
  shears the noise circle into an ellipse.
- **Cavity-mediated interaction.** The spins do not touch each other
  directly; they exchange energy through the cavity. With the cavity
  tuned away from the spins by `Delta` (the **detuning**), this gives
  a twisting strength `chi` and a collective decay rate `Gamma_SR`
  (formulas in `help(cavsqueeze.resonator)`). `g` is the coupling
  of one spin to the cavity and `kappa` the cavity linewidth.
- **Inhomogeneous line.** In a solid the spins do not all have the
  same frequency. The spread is the **line**, described by its shape
  (Gaussian, Lorentzian or Voigt) and its full width at half maximum
  (FWHM).
- **Cumulant expansion.** The solver follows averages and pairwise
  correlations of the spins and drops higher-order correlations
  (a "second-order", Gaussian approximation). Tests compare it with
  exact solutions where those exist.
- **Echo.** A pi pulse (a pulse that turns every spin by half a
  turn) half way through the twisting. It undoes the dephasing
  (spreading out of the spins' directions) caused by static frequency
  differences.
- **Ramsey measurement, tomography.** A Ramsey sequence turns an
  accumulated phase into a population difference. **Fringe**: the
  mean `J_z` as the phase of the last pulse is scanned. **Tomography**:
  the shot-to-shot variance of `J_z` as the noise ellipse is rotated
  by an angle `theta`; its smallest value is the squeezed variance.
- **Allan deviation.** The standard measure of a clock's fractional
  frequency instability after averaging for a time `tau`.
- **Dick effect.** A pulsed clock looks at its laser (the local
  oscillator) only part of the time, so laser noise leaks into the
  result. It sets a floor that no squeezing can lower.
- **Exact references.** QuTiP is a Python library for quantum
  simulations. A **master-equation** solution computes the full
  quantum state, including decay; it is exact but only feasible for a
  few spins. **PIQS** (part of QuTiP) solves the same problem exactly
  for many *identical* spins by using their symmetry (the **Dicke**
  basis).

## Install, requirements and units

```
pip install cavsqueeze          # core package
pip install cavsqueeze[exact]   # adds QuTiP, needed for cavsqueeze.exact and to_qutip
```

To run the tests from a copy of the source:

```
git clone https://github.com/TaN-MM-Org/cavsqueeze
cd cavsqueeze
pip install -e .[test]
pytest tests
```

It needs Python 3.10 or newer, NumPy 1.24 or newer, SciPy 1.10 or
newer and Matplotlib 3.7 or newer (Matplotlib is a declared
dependency, although no module of the package imports it). QuTiP 5.0
or newer is optional. QuTiP 5.0.0 imports `setuptools`, so with that
QuTiP version `setuptools` must also be installed.

Units and conventions:

- Inside the package every rate and frequency is an **angular**
  frequency in rad/s. `from_hz` takes ordinary frequencies in Hz and
  multiplies them by 2 pi. Line widths passed to `lineshape` are also
  in rad/s.
- Temperature `T` in kelvin, times in seconds, `T2` in seconds
  (`from_hz` sets the dephasing rate to `1/T2`).
- `J_z` and all variances are in spin units: `J_z = (N_up - N_down)/2`,
  variances in spin units squared. Angles and phases are in radians.
- The spin state variable is `<sigma^+> = (x + i y)/2`. A spin whose
  frequency is above the reference (detuning `delta > 0`) turns from
  `+x` towards `+y`.
- The protocols start from the coherent spin state along `+x`.

## Examples

Each example runs as written, and the output shown is what it printed
with cavsqueeze 1.17.0. All device values are illustrative
(chosen to show the behavior), not the values of a specific
experiment. Examples that draw random numbers use a fixed seed.

### 1. How much squeezing does a cavity give, and when?

```python
import numpy as np
from cavsqueeze import from_hz, homogeneous, optimal_squeezing

# Illustrative values, not a specific device.
N = 1e10                                   # number of spins
p = from_hz(g_hz=1e6 / np.sqrt(N),         # coupling of one spin (Hz)
            kappa_hz=1e4,                  # cavity linewidth (Hz)
            Delta_hz=30e6,                 # cavity-spin detuning (Hz)
            T=0.02,                        # temperature (K)
            T2=0.15)                       # spin coherence time (s)
best = optimal_squeezing(p, homogeneous(N), 1e-6, 1e-2)
print(f"best squeezing {10 * np.log10(best['xi2']):.2f} dB "
      f"after {best['t'] * 1e6:.1f} us, contrast {best['contrast']:.4f}")
```

```
best squeezing -23.59 dB after 124.3 us, contrast 0.9992
```

`optimal_squeezing` scans the twisting time between the two limits
given (here 1 us to 10 ms), with an echo pulse in the middle, and
returns the time with the smallest `xi2_R`. It first scans a coarse
grid of 12 times spaced evenly on a log scale, then a finer linear
grid between the neighbors of the best coarse point. `best` also
holds the whole scan (`scan_t`, `scan_xi2`).

### 2. A disordered line, with 10^15 spins

```python
import numpy as np
from cavsqueeze import (from_hz, homogeneous, lineshape,
                        equal_probability_classes, squeezing_after)

N = 1e15
p = from_hz(g_hz=1e6 / np.sqrt(N), kappa_hz=1e4, Delta_hz=30e6,
            T=0.02, T2=0.15)
# A Gaussian line 2 kHz wide (FWHM), cut into 16 classes of equal weight.
line = lineshape("gaussian", 2 * np.pi * 2e3)
for name, ens in [("one class, no disorder", homogeneous(N)),
                  ("16 classes, 2 kHz line", equal_probability_classes(line, 16, N))]:
    r = squeezing_after(p, ens, 2e-5)
    print(f"{name:24s} xi2 = {10 * np.log10(r['xi2']):6.2f} dB, "
          f"contrast {r['contrast']:.4f}")
```

```
one class, no disorder   xi2 = -12.86 dB, contrast 0.9999
16 classes, 2 kHz line   xi2 = -12.84 dB, contrast 0.9972
```

Both runs take seconds: the size of the problem is the number of
classes (1 or 16), not the number of spins. `squeezing_after` uses an
echo by default, which removes most of the effect of the static
frequency spread; the contrast still drops a little.

### 3. What the squeezed state buys a clock

```python
import numpy as np
from cavsqueeze import (Rates, from_hz, homogeneous, css_x, twist,
                        squeezing_parameters, metrological_gain_db,
                        clock_allan_deviation)

N = 1e8
p = from_hz(g_hz=1e6 / np.sqrt(N), kappa_hz=1e4, Delta_hz=30e6,
            T=0.02, T2=0.15)
ens = homogeneous(N)
rt = Rates.from_params(p, ens)            # the model's rates, per class
state = twist(css_x(rt.M), rt, 5e-5)      # 50 us of twisting, with an echo

m = squeezing_parameters(state, ens.n)
print(f"xi2_S = {m['xi2_S']:.5f}, xi2_R = {m['xi2_R']:.5f}, "
      f"contrast = {m['contrast']:.5f}")
print(f"gain over the standard quantum limit: "
      f"{metrological_gain_db(m['xi2_R']):.2f} dB")
print(f"phase noise per shot: {m['dphi']:.3e} rad "
      f"(standard quantum limit {1 / np.sqrt(N):.3e} rad)")
# Illustrative clock: 429 THz transition, 0.1 s Ramsey time, 1 s averaging.
print(f"projected Allan deviation: "
      f"{clock_allan_deviation(m['dphi'], nu0=4.29e14, T_ramsey=0.1, tau=1.0):.3e}")
```

```
xi2_S = 0.01013, xi2_R = 0.01013, contrast = 0.99967
gain over the standard quantum limit: 19.94 dB
phase noise per shot: 1.007e-05 rad (standard quantum limit 1.000e-04 rad)
projected Allan deviation: 1.181e-20
```

`twist` evolves the state; `squeezing_parameters` reads both squeezing
parameters, the contrast and the per-shot phase noise `dphi` from it.
`clock_allan_deviation` uses the standard projection-noise formula
`sigma_y = dphi / (2 pi nu0 T_ramsey) * sqrt(T_cycle / tau)`, with the
cycle time equal to the Ramsey time unless you give `T_cycle`.

### 4. The exact one-axis-twisting benchmark

```python
from cavsqueeze import oat_closed_form

for mu in (0.0, 0.05, 0.2):
    r = oat_closed_form(100, mu)
    print(f"mu = {mu:4.2f}: <Jx> = {r['Jx']:7.3f}, Vmin = {r['Vmin']:6.3f}, "
          f"xi2_S = {r['xi2_S']:.4f}, xi2_R = {r['xi2_R']:.4f}")
```

```
mu = 0.00: <Jx> =  50.000, Vmin = 25.000, xi2_S = 1.0000, xi2_R = 1.0000
mu = 0.05: <Jx> =  48.477, Vmin =  3.128, xi2_S = 0.1251, xi2_R = 0.1331
mu = 0.20: <Jx> =  30.453, Vmin =  3.946, xi2_S = 0.1578, xi2_R = 0.4255
```

`oat_closed_form(N, mu)` gives the Kitagawa-Ueda closed-form results for
N spins twisted by the angle `mu` with no decay and no disorder
(`mu = 2 chi t`). The dissipative solver can be compared against it.
At `mu = 0.2` the smallest variance has grown again and the contrast
has fallen, so `xi2_R` is much worse than `xi2_S`: twisting too long
costs precision.

### 5. Squeezing from measured shot data

```python
import os, tempfile
import numpy as np
from cavsqueeze import (oat_closed_form, save_shots_csv, load_shots_csv,
                        estimate_contrast, estimate_squeezing,
                        metrological_gain_db)

# Make a synthetic, seeded experiment: N = 100 spins, one-axis-twisted
# by mu = 0.05 (the exact answer is known from oat_closed_form).
N, rng = 100, np.random.default_rng(0)
oat = oat_closed_form(N, 0.05)
C_true = 2 * oat["Jx"] / N                         # true fringe contrast

# Fringe scan: mean J_z = (N C / 2) cos(phase), plus projection noise.
phases = np.linspace(0, 2 * np.pi, 12, endpoint=False)
fringe = [0.5 * N * C_true * np.cos(ph) + rng.normal(0, 5, 400)
          for ph in phases]
# Tomography: J_z variance at angle theta follows the squeezed ellipse.
a = oat["alpha_min"]
cov = np.array([[np.cos(a), -np.sin(a)], [np.sin(a), np.cos(a)]]) @ \
    np.diag([oat["Vmin"], oat["Vmax"]]) @ \
    np.array([[np.cos(a), np.sin(a)], [-np.sin(a), np.cos(a)]])
angles = np.linspace(0, np.pi, 12, endpoint=False)
tomo = [rng.normal(0, np.sqrt(np.array([np.cos(t), np.sin(t)]) @ cov
                              @ np.array([np.cos(t), np.sin(t)])), 4000)
        for t in angles]

folder = tempfile.mkdtemp()
save_shots_csv(os.path.join(folder, "fringe.csv"), phases, fringe)
save_shots_csv(os.path.join(folder, "tomography.csv"), angles, tomo)

ph, fr = load_shots_csv(os.path.join(folder, "fringe.csv"))
cest = estimate_contrast(ph, fr, N=N)
th, sh = load_shots_csv(os.path.join(folder, "tomography.csv"))
est = estimate_squeezing(th, sh, N=N, contrast=min(cest.contrast, 1.0),
                         contrast_sigma=cest.contrast_sigma)
print(f"contrast: {cest.contrast:.4f} +/- {cest.contrast_sigma:.4f} "
      f"(true {C_true:.4f})")
print(f"xi2_R:    {est.xi2_R:.4f} +/- {est.xi2_R_sigma:.4f} "
      f"(true {oat['xi2_R']:.4f})")
print(f"gain:     {metrological_gain_db(est.xi2_R):.2f} dB")
```

```
contrast: 0.9713 +/- 0.0020 (true 0.9695)
xi2_R:    0.1420 +/- 0.0043 (true 0.1331)
gain:     8.48 dB
```

The data here are synthetic, drawn from the exact twisted state of
example 4 so the true answer is known. With real data you would only
run the last block. The steps:

1. `save_shots_csv` / `load_shots_csv` read and write shot records in a
   plain-text format: a header `angle_rad,jz`, then one row per shot.
   Rows with the same angle form that angle's shots.
2. `estimate_contrast` fits the fringe. The mean follows
   `A cos(phi) + B sin(phi) + d` for any state, so the fit is a
   weighted linear least-squares fit: a standard fit whose answer is
   given by a formula, with no iterative search.
3. `estimate_squeezing` fits the tomography. The variance at angle
   `theta` follows `c + a cos(2 theta) + b sin(2 theta)` for any state
   (it is how a 2x2 covariance matrix rotates), so this fit is also
   linear. Each angle is weighted by how precisely its variance is
   known, which is proportional to the variance itself; since 1.16 the
   variance used for that weight is read from the fitted curve, and
   the fit is repeated until it no longer changes (see
   [Corrections in earlier versions](#corrections-in-earlier-versions)
   for why). By default the uncertainty of each sample variance of `M`
   shots is taken as `sigma^2 sqrt(2/(M-1))`, which holds for Gaussian
   shot noise; `shot_noise="empirical"` measures it from the shots
   instead (example 9). A known detection-noise variance can be
   subtracted with `detection_variance=...`.

In this run the estimate is about two error bars from the true value,
which is within normal scatter for one experiment. (With 1.15.1 this
example printed `xi2_R: 0.1421` and `gain: 8.47 dB`; the change is the
new fit weights.)

### 6. Planning the number of shots

```python
import numpy as np
from cavsqueeze import plan_tomography, shots_for_squeezing

# Illustrative plan: 480 spins, expected variances 90 and 160 (spin units^2),
# contrast 0.92 +/- 0.01, nine equally spaced angles.
angles = np.linspace(0.0, np.pi, 9, endpoint=False)
plan = plan_tomography(angles, shots_per_angle=50, var_min=90.0,
                       var_max=160.0, N=480, contrast=0.92,
                       contrast_sigma=0.01)
print(f"50 shots per angle: xi2_R = {plan['xi2_R']:.4f} "
      f"+/- {plan['xi2_R_sigma']:.4f}")

m, plan = shots_for_squeezing(0.02, angles, 90.0, 160.0, N=480,
                              contrast=0.92, contrast_sigma=0.01)
print(f"for +/- 0.02: {m} shots per angle (gives {plan['xi2_R_sigma']:.4f})")

try:
    shots_for_squeezing(0.015, angles, 90.0, 160.0, N=480,
                        contrast=0.92, contrast_sigma=0.01)
except ValueError as err:
    print("refused:", err)
```

```
50 shots per angle: xi2_R = 0.8861 +/- 0.1159
for +/- 0.02: 22116 shots per angle (gives 0.0200)
refused: target 0.015 is below the contrast floor 0.0193 = xi2_R * 2 sigma_C / C: no number of tomography shots reaches it -- improve the fringe (contrast) measurement instead
```

Because the tomography fit is linear, its error bars depend only on
the design (angles, shots per angle, the variances you expect), so
they can be computed before measuring. `plan_tomography` returns them;
`shots_for_squeezing` finds the smallest number of shots per angle
that meets a target error bar on `xi2_R`. The contrast error sets a
floor, `xi2_R * 2 sigma_C / C`, that more tomography shots cannot
lower, so a target below it is refused. Equally spaced angles over
half a turn are the efficient choice: they make the three fitted
numbers independent.

### 7. When the laser, not the atoms, limits a clock

```python
from cavsqueeze import (power_law_psd, total_clock_allan_deviation,
                        synchronized_comparison)

# Illustrative clock: 429 THz, 0.4 s Ramsey time in a 1 s cycle.
# The laser noise level h0 is a placeholder -- use your laser's measured value.
laser = power_law_psd(h0=1e-30)
for label, dphi in [("no squeezing", 1e-2), ("10 dB squeezing", 3.16e-3)]:
    r = total_clock_allan_deviation(dphi, nu0=4.29e14, T_ramsey=0.4,
                                    tau=1.0, T_cycle=1.0, S_y=laser)
    print(f"{label:16s} projection {r['qpn']:.2e}  Dick {r['dick']:.2e}  "
          f"total {r['total']:.2e}  Dick-limited: {r['dick_limited']}")

# Two ensembles read by the same laser: the Dick term cancels in the difference.
out = synchronized_comparison(3.16e-3, 3.16e-3, nu0=4.29e14, T_ramsey=0.4,
                              tau=1.0, T_cycle=1.0)
print(f"synchronized comparison: difference {out['differential']:.2e}, "
      f"per clock {out['per_clock']:.2e}")
```

```
no squeezing     projection 9.27e-18  Dick 8.66e-16  total 8.66e-16  Dick-limited: True
10 dB squeezing  projection 2.93e-18  Dick 8.66e-16  total 8.66e-16  Dick-limited: True
synchronized comparison: difference 4.14e-18, per clock 2.93e-18
```

With this laser noise the Dick floor is about 100 times the projection
noise, so 10 dB of squeezing leaves the total unchanged
(`dick_limited` is `True`). The Dick floor is computed from the
laser's noise spectrum through the Fourier coefficients of the Ramsey
timing, in closed form; `power_law_psd` builds the usual
`h0 + h_-1/f + h_-2/f^2` spectrum from your measured coefficients (the
package ships none).

Record clock comparisons avoid the Dick effect by reading two
ensembles with the same laser in the same cycles, so the laser noise
is common to both and cancels in their difference (Robinson et al.,
Nat. Phys. 20, 208 (2024): 1.9(2) dB of stability enhancement at the
10^-17 level; Yang et al., PRL 135, 193202 (2025): 2.0(2) dB beyond
the standard quantum limit at 1.1x10^-18). `synchronized_comparison`
models exactly that: the difference has the two projection noises
added in quadrature and no Dick term, and `per_clock` is the
difference divided by sqrt(2). Noise that is not common to both
ensembles (for example path-length differences or gradients) is not
modeled and adds on top.

### 8. Handing the state to QuTiP

```python
import numpy as np
from cavsqueeze import (Rates, from_hz, homogeneous, css_x, twist, bosonic_mode,
                        to_qutip)

N = 1e8
p = from_hz(g_hz=1e6 / np.sqrt(N), kappa_hz=1e4, Delta_hz=30e6,
            T=0.02, T2=0.15)
rt = Rates.from_params(p, homogeneous(N))
state = twist(css_x(rt.M), rt, 2e-5)

mode = bosonic_mode(state, rt.n)
print(f"r = {mode['r']:.4f}, thermal occupation = {mode['n_th']:.5f}, "
      f"purity = {mode['purity']:.5f}, squeezing = {mode['squeezing_db']:.2f} dB")
rho, mode = to_qutip(state, rt.n)            # needs QuTiP
print("QuTiP density matrix:", rho.shape, f"trace {rho.tr().real:.6f}")
```

```
r = 1.4828, thermal occupation = 0.00239, purity = 0.99525, squeezing = -12.86 dB
QuTiP density matrix: (216, 216) trace 1.000000
```

For a strongly polarized spin, the small fluctuations across the mean
spin behave like one bosonic mode (the Holstein-Primakoff picture).
`bosonic_mode` gives that mode's 2x2 covariance, scaled so that a
coherent spin state is the vacuum (variance 1/2), and splits it into
a squeezing strength `r`, an angle, a thermal occupation and a purity.
`to_qutip` builds the matching QuTiP density matrix, then checks that
the exported state reproduces the two extreme variances to 1 %, and
raises if the Fock-space cutoff (the number of photon-number levels
QuTiP keeps) was too small. It keeps second moments (variances and
covariances) only.

### 9. Error bars when the detector has outliers

```python
import numpy as np
from cavsqueeze import estimate_squeezing

# Synthetic tomography of a squeezed state of N = 100 spins (variances 5 and
# 40 in spin units squared, xi2_S = 0.2), read out by a detector that gives
# 5 % of shots with three times the usual spread (outliers).
rng = np.random.default_rng(1)
angles = np.linspace(0, np.pi, 6, endpoint=False)
var = 22.5 - 17.5 * np.cos(2 * (angles - 0.3))
z = rng.normal(0, 1, (6, 400))
z = np.where(rng.random((6, 400)) < 0.05, 3 * z, z) / np.sqrt(0.95 + 0.05 * 9)
shots = z * np.sqrt(var)[:, None]

for noise in ("gaussian", "empirical"):
    est = estimate_squeezing(angles, shots, N=100, contrast=1.0,
                             shot_noise=noise)
    print(f"shot_noise={noise!r:12s} xi2_S = {est.xi2_S:.4f} "
          f"+/- {est.xi2_S_sigma:.4f}")
```

```
shot_noise='gaussian'   xi2_S = 0.2202 +/- 0.0182
shot_noise='empirical'  xi2_S = 0.2202 +/- 0.0382
```

The error bar of a squeezing estimate comes from how much each
angle's sample variance scatters. For Gaussian shots that scatter is
fixed by the variance alone; with outliers (a heavy-tailed
distribution) it is larger. `shot_noise="empirical"` measures it from
the shots through the **kurtosis** (the fourth moment of the shots
divided by the squared variance; 3 for a Gaussian), using the exact
formula `Var(s^2) = sigma^4 [kurtosis/M - (M-3)/(M(M-1))]` for `M`
shots, and carries it through the fit. The estimate itself is the same
either way; only the error bar changes. Here the true value 0.2 is
1.1 Gaussian error bars away but 0.5 empirical ones. In the tests, over
200 experiments like this one the Gaussian error bar is too small by
more than a factor 1.5, and the empirical one matches the actual
scatter to within 20 %. `sample_variance_sigma(shots)` gives the same
error bar for a single shot array.

### 10. T2 and the line width from your own decays

```python
import numpy as np
from cavsqueeze import fit_coherence_decay, from_hz, lineshape

# Synthetic, seeded contrasts standing in for a measurement (for example
# estimate_contrast at each time), each with error 0.01.  Truth: T2 = 5 ms,
# a Gaussian line of FWHM 150 Hz, contrast 0.9 at t = 0.
rng = np.random.default_rng(2)
err = np.full(12, 0.01)
s = 2 * np.pi * 150.0 / (2 * np.sqrt(2 * np.log(2)))    # line std, rad/s
t_echo = np.linspace(1e-3, 15e-3, 12)                    # s, total echo time
echo = 0.9 * np.exp(-t_echo / 5e-3) + rng.normal(0, 0.01, 12)
t_ramsey = np.linspace(0.2e-3, 4e-3, 12)                 # s
ramsey = (0.9 * np.exp(-t_ramsey / 5e-3 - 0.5 * (s * t_ramsey) ** 2)
          + rng.normal(0, 0.01, 12))

e = fit_coherence_decay(t_echo, echo, err)                  # echo: T2 only
print(f"echo:   T2 = {e.T2 * 1e3:.2f} +/- {e.T2_sigma * 1e3:.2f} ms, "
      f"p = {e.p_value:.2f}")
r = fit_coherence_decay(t_ramsey, ramsey, err, line="gaussian",
                        T2=e.T2, T2_sigma=e.T2_sigma)       # Ramsey: the line
print(f"Ramsey: FWHM = {r.fwhm / (2 * np.pi):.1f} +/- "
      f"{r.fwhm_sigma / (2 * np.pi):.1f} Hz, p = {r.p_value:.2f}")
bad = fit_coherence_decay(t_ramsey, ramsey, err)            # wrong model
print(f"Ramsey fitted without a line: p = {bad.p_value:.1e}")

# The fitted values are the model's inputs (rates in rad/s, times in s):
params = from_hz(g_hz=1e-4, kappa_hz=1e4, Delta_hz=30e6, T=0.02, T2=e.T2)
line = lineshape("gaussian", r.fwhm)
```

```
echo:   T2 = 5.05 +/- 0.09 ms, p = 0.24
Ramsey: FWHM = 149.4 +/- 2.3 Hz, p = 0.98
Ramsey fitted without a line: p = 4.8e-29
```

The model needs the dephasing rate `1/T2` and the line of the spin
frequencies. `fit_coherence_decay` (new in 1.17) fits both from decays
you measure. Measure them with the interaction switched off or too weak
to matter on these time scales (for example with the cavity tuned far
away); otherwise the fitted rates include the interaction.

- An echo decays as `A exp(-t/T2)`, with `t` the total time of both
  halves: the pi pulse removes the effect of the frequency spread, so
  `line="none"` gives `T2`.
- A Ramsey decay also contains the line. For a Gaussian line of FWHM
  `F` (standard deviation `s = F / (2 sqrt(2 ln 2))`) it is
  `A exp(-t/T2 - s^2 t^2 / 2)`; for a Lorentzian line of FWHM `F` it is
  `A exp(-(1/T2 + F/2) t)`. The Ramsey fit takes `T2` and its error bar
  from the echo fit and propagates that error into the width.

The fit is a weighted least-squares fit of the contrast itself. Its
error bars come from the curvature of the fit (the standard linear
approximation); in the tests they match the scatter over 300 seeded
datasets to 20 %. `chi2`, `dof` and `p_value` say whether the decay
follows the formula within the error bars: the last line shows a
Ramsey decay fitted without its Gaussian line, which the p-value
rejects. The fit refuses a Lorentzian width without a known `T2`
(only the sum `1/T2 + F/2` can be measured from a Ramsey decay) and
any rate or width that is not at least two error bars above zero.
The fitted values (in rad/s and s) go straight into `from_hz` and
`lineshape`, as the last two lines show.

### 11. Do the tomography data follow the rotation law?

```python
import numpy as np
from cavsqueeze import estimate_squeezing

# Synthetic tomography of a squeezed state (variances 5 and 40 in spin
# units squared, N = 100), 12 angles, 2000 shots each.  In the second
# record the rotation actually applied was 5 % larger than recorded.
rng = np.random.default_rng(3)
angles = np.linspace(0, np.pi, 12, endpoint=False)
for label, applied in [("angles as recorded", angles),
                       ("angles 5 % too large", 1.05 * angles)]:
    var = 22.5 - 17.5 * np.cos(2 * (applied - 0.3))
    shots = rng.normal(0, 1, (12, 2000)) * np.sqrt(var)[:, None]
    est = estimate_squeezing(angles, shots, N=100, contrast=1.0)
    print(f"{label:21s} xi2_S = {est.xi2_S:.4f} +/- {est.xi2_S_sigma:.4f}, "
          f"chi2 = {est.chi2:4.1f} ({est.dof} dof), p = {est.p_value:.1g}")
```

```
angles as recorded    xi2_S = 0.2048 +/- 0.0048, chi2 =  5.6 (9 dof), p = 0.8
angles 5 % too large  xi2_S = 0.1979 +/- 0.0046, chi2 = 101.1 (9 dof), p = 9e-18
```

The variance at angle `theta` follows `c + a cos(2 theta) + b sin(2 theta)`
for any state, so a fit that does not describe the data points at the
measurement, not the state: wrong or drifting rotation angles, a state
that changed between angles, or detection noise that depends on the
angle. Since 1.17 `estimate_squeezing` reports `chi2`, its degrees of
freedom (angles minus 3) and the p-value, and `estimate_contrast` does
the same for the fringe. In the second record above the squeezing
estimate looks normal, but `p` is tiny. The test assumes nearly
Gaussian sample variances, which needs many shots per angle. In the
tests, with 200 shots per angle, between 2 % and 9 % of 400 correct
experiments give `p < 0.05` (5 % is expected).

## What is in the package

Every name below is importable from `cavsqueeze` directly. Each
docstring (`help(cavsqueeze.estimate_squeezing)`, for example) gives the
inputs, units and conventions.

**Cavity and spins**

- `CavityParams`, `from_hz` -- coupling `g`, cavity linewidth `kappa`,
  detuning `Delta`, spin frequency, temperature and dephasing; the
  twisting and collective-decay rates follow from them.
- `thermal_occupation` -- the thermal photon number at a frequency and
  temperature.
- `loop_gap_dispersive(N, T, T2)` -- the parameters of the dispersive
  loop-gap resonator of Fukumori et al. (arXiv:2604.26909):
  g/2pi = 15 mHz, kappa/2pi = 660 kHz, Delta/2pi = 22 MHz.
- `Ensemble`, `homogeneous` -- a set of classes (detuning, coupling
  weight, number of spins); one class with no disorder.
- `lineshape` -- Gaussian, Lorentzian or Voigt line from its FWHM. For
  the Voigt line (a Gaussian and a Lorentzian combined),
  `lorentz_fraction` sets the Lorentzian part's FWHM as a fraction of
  the total, and since 1.16 the Gaussian part is solved for so that the
  total FWHM is the one asked for, checked to 1e-9 (`split="olivero"`
  gives the approximate split used before).
- `equal_probability_classes` -- cuts a line into M classes of equal
  weight. `product_classes` combines this with a spread of coupling
  strengths, and `log_uniform_weights` makes such a spread.
  (`cavsqueeze.ensemble.tail_resolved_classes` resolves the far tails
  of a Lorentzian line and treats spins far outside it as free
  "spectator" spins.)

**Solver** (second-order cumulant expansion)

- `Rates` -- the model rates per class, from `Rates.from_params(params,
  ensemble)`. Setting `Rates.meas` (and `meas_eta`, the detection
  efficiency) adds a continuous measurement of `J_z` through the
  resonator and conditions the state on it.
- `State`, `product_state` -- the solver state: per-class averages and
  pair correlations.
- `evolve`, `evolve_meanfield` -- time evolution with and without
  correlations.
- `rotate` -- an instantaneous rotation (pulse) of all spins.
- `collective_moments`, `transverse_variances`, `wineland_xi2`,
  `coherence` -- mean collective spin and covariance; the variances
  across the mean spin; the Wineland parameter; the Ramsey contrast.

**Pulse sequences**

- `css_x` -- the starting state (coherent spin state along `+x`).
- `twist`, `twist_imperfect`, `pulse` -- twisting with or without an
  echo; an echo with a pulse of finite length or uneven angle; a
  pulse during which the interaction continues.
- `squeezing_after`, `squeezing_trace`, `optimal_squeezing` --
  squeezing after one time, over a list of times, and at the best
  time.
- `ramsey_cumulant`, `ramsey_meanfield` -- Ramsey decay of the
  contrast.
- `twist_untwist` -- the twist-untwist readout of Davis, Bentsen,
  Schleier-Smith, PRL 116, 053601 (2016), with detection noise.
- `plain_squeezed_readout` -- the gain of an ordinary squeezed Ramsey
  readout with detection noise.

**Metrology**

- `squeezing_parameters`, `phase_sensitivity`, `metrological_gain_db`
  -- `xi2_S`, `xi2_R`, contrast and per-shot phase noise; gain in dB.
- `clock_allan_deviation`, `magnetometer_sensitivity` -- the
  projection-noise-limited clock instability and Ramsey magnetometer
  field noise.
- `oat_closed_form` -- exact one-axis twisting without decay.

**Dick effect**

- `ramsey_sensitivity`, `dick_fourier_coefficients` -- the Ramsey
  sensitivity function over one cycle and its Fourier coefficients in
  closed form.
- `power_law_psd` -- the laser noise spectrum from its coefficients.
- `dick_allan_deviation`, `total_clock_allan_deviation` -- the Dick
  floor; the floor and the projection noise together, with a
  `dick_limited` verdict.
- `synchronized_comparison` -- two ensembles read by one laser.

**Measured data**

- `save_shots_csv`, `load_shots_csv` -- the shot-record file format.
- `estimate_contrast` (returns `ContrastEstimate`) -- fringe contrast
  and its error bar.
- `estimate_squeezing` (returns `SqueezingEstimate`) and the fit behind
  it, `variance_tomography` -- squeezing parameters and error bars
  from tomography shots. Options (new in 1.16): `weighting="model"`
  (default) or `"sample"` (the 1.15 fit weights), and
  `shot_noise="gaussian"` (default) or `"empirical"` (example 9).
  With few shots per angle the model weights cannot always be
  computed (the fitted curve can dip to zero or below at a measured
  angle). The default then returns the `"sample"` fit instead, warns
  (`RuntimeWarning`) and records why in `weighting_fallback`, so it
  never refuses data the 1.15 weights can fit. In a seeded script (not
  part of the test suite; setup in CHANGELOG 1.16.0) with 6 angles this
  happened for 27-29 % of the datasets with 3 shots per
  angle, 11-13 % with 5 and 1-2 % with 10.
- `sample_variance_sigma` -- the error bar of one sample variance,
  from the measured kurtosis or a given one.
- Fit quality (new in 1.17): `SqueezingEstimate` and `ContrastEstimate`
  carry `chi2`, `dof` and `p_value`, which test whether the data follow
  the rotation law or the fringe law (example 11); `variance_tomography`
  returns the same three numbers.
- `plan_tomography`, `shots_for_squeezing` -- predicted error bars and
  the shot count for a target.
- `fit_coherence_decay` (returns `DecayFit`, new in 1.17): `T2` from
  an echo decay, and the FWHM of a Gaussian or Lorentzian line from a
  Ramsey decay, with error bars and a fit-quality test (example 10).

**Bridge to QuTiP**

- `bosonic_mode`, `to_qutip` -- see example 8.

**Reference solvers** (not imported at the package root)

- `cavsqueeze.exact` (needs QuTiP) -- exact master-equation solutions
  for a few distinguishable spins (`full_hilbert`) and the
  permutation-invariant Dicke solver PIQS for identical spins
  (`dicke_piqs`).
- `cavsqueeze.cumulant_raw` -- the same cumulant equations written in
  raw moments, kept as a test reference.
- `cavsqueeze.dtwa` -- an independent solver based on the discrete
  truncated Wigner approximation (J. Schachenmayer, A. Pikovski and
  A. M. Rey, Phys. Rev. X 5, 011022 (2015)): many classical
  trajectories started from random samples of the initial state. It
  covers twisting and detunings (no decay).

## When it refuses, and why

`cavsqueeze` raises `ValueError` instead of guessing when:

- `lineshape` is given a name other than gaussian, lorentzian or voigt,
  a width that is not finite and positive, or a `lorentz_fraction`
  outside [0, 1] (new in 1.16; before, a negative width gave `nan`
  detunings without an error);
- an `Ensemble` has arrays of different lengths, non-finite values,
  negative occupations or no spins; `equal_probability_classes` gets
  fewer than 1 class; `product_classes` gets negative or all-zero
  probabilities (new in 1.16);
- cavity parameters are not finite, `kappa`, the temperature or the
  dephasing rate is negative, `kappa` and `Delta` are both zero, or
  `T2` is not positive (new in 1.16; before, a negative `T2` gave a
  negative dephasing rate without an error);
- an evolution time is negative, or `optimal_squeezing` gets limits
  that do not satisfy `0 < t_lo < t_hi` (new in 1.16; since 1.17
  `evolve_meanfield` refuses a negative or non-finite time too, where it
  used to integrate backwards and return a contrast above 1);
- a clock or magnetometer projection gets a phase noise, clock
  frequency, Ramsey time, averaging time or gyromagnetic ratio that is
  not finite and positive (only its size is used for the gyromagnetic
  ratio), `metrological_gain_db` gets `xi2_R <= 0`, or the Dick floor
  gets an averaging time that is not positive (new in 1.16; before,
  a negative `tau` returned `nan`);
- `oat_closed_form` is asked for fewer than 2 spins, or (new in 1.17)
  for a number of spins that is not a whole number (2.9 used to be
  truncated to 2 without a word) or a twisting angle that is not
  finite;
- a clock or magnetometer cycle time `T_cycle` is shorter than the
  Ramsey time, or (new in 1.17) not finite (before, `nan` or `inf` came
  back as the result);
- the Ramsey timing is impossible (`T_ramsey <= 0`, a negative pulse
  length, or a cycle shorter than the Ramsey time plus two pulses), any
  of its times is not finite (new in 1.17; before, the Dick floor came
  back as `nan`), or fewer than 1 harmonic is asked for;
- laser noise coefficients are negative or (new in 1.17) not finite;
- `log_uniform_weights` gets fewer than 1 weight or an infinite range
  (new in 1.17);
- the harmonic cutoff of `dick_allan_deviation` is too low: if the top
  10 % of the harmonics still carry more than `tail_tol` (default
  1e-3) of the sum, it asks for a larger `m_max` rather than returning
  a truncated answer;
- `synchronized_comparison` gets a per-shot phase noise that is not
  finite and positive;
- the mean spin is zero (there is then no direction to squeeze across),
  the covariance is not positive definite, or it breaks the
  uncertainty relation (`bosonic_mode`);
- the Fock-space cutoff distorts the state exported by `to_qutip`;
- a fit has fewer than 3 angles, or angles that give fewer than 3
  different directions (angles that differ by a half turn are the same
  direction for the tomography; by a full turn for the fringe);
- a shot array has fewer than 2 shots, non-finite values or zero
  scatter, or the number of shot arrays does not match the angles;
- `N` given to the measured-data functions is not an integer (a Python
  or NumPy integer; `100.0` is refused) of at least 2;
- the contrast given to `estimate_squeezing` or `plan_tomography` is
  not in (0, 1], or its error is negative or (new in 1.17) not finite;
- the fitted fringe contrast exceeds 1 by more than two error bars
  (usually a wrong `N` or a `J_z` calibration problem; the message
  says so). A smaller overshoot is returned unclipped;
- the detection variance is negative or not finite, or subtracting it
  leaves a smallest variance of zero or less;
- a shot file is empty, has the wrong header, a row without exactly two
  fields, no data rows, or an angle with fewer than 2 shots; saving
  refuses repeated angles;
- the planned variances do not satisfy `0 < var_min <= var_max`, the
  planned shot counts are not whole numbers (new in 1.17; 9.9 used to
  be planned as 9), or a target error bar is not positive, uses a
  degenerate angle set, or is below the contrast floor (example 6);
- `fit_coherence_decay` (new in 1.17) gets a Lorentzian line without a
  known `T2`, `T2` together with `line="none"`, inputs of different
  lengths, negative or non-finite times, non-finite contrasts, error
  bars that are not finite and positive, too few points (one more than
  the fitted parameters) or distinct times, or times that leave the
  parameters undetermined; or when a fitted rate or width is not at
  least `min_significance` (default 2) error bars above zero.

Also: `evolve` raises `RuntimeError` if SciPy's integrator reports a
failure, and a planned design with fewer than 3 distinct directions is
reported by `plan_tomography` as `identifiable = False` (and refused
by `shots_for_squeezing`).

## How the results are checked

170 tests run on every push to `main` and every pull request, on
Python 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.24.0, SciPy 1.10.0,
Matplotlib 3.7.0, QuTiP 5.0.0). A further job runs on Python 3.14
without QuTiP; there the tests that need QuTiP are skipped. 3 more
tests always skip here: they check scripts of the paper's companion
repository, which are not in this repository. The whole suite takes a
few minutes; most of it is one test of the trajectory solver.

**The cumulant solver**

- For 2 spins the second-order equations are complete, so the solver
  must equal the exact QuTiP master-equation solution: every pair
  correlation agrees to 2e-6 and every average to 2e-6 plus 1e-5 of
  its size, for three random sets of detunings and couplings, with
  decay, thermal absorption and dephasing switched on.
- For 4 spins it agrees with the exact solution to 2e-3 at short times
  and to 0.15 at the end of the run (the error of the approximation).
- For 40 identical spins with decay, thermal absorption and dephasing,
  its Wineland parameter agrees with the exact PIQS solution to 0.2 dB
  and its mean spin to 5 % of `N/2`, up to a twisting strength
  `chi N t = 1.5`.
- The production solver (correlations stored with the averages
  subtracted) agrees with the raw-moment form: averages to 1e-8 plus
  1e-5 of their size, correlations to 1e-7, the Wineland parameter to
  1e-7, for classes of one and of several spins.
- Rotations: a coherent spin state keeps `xi2 = 1` to 1e-12, and
  rotating back and forth returns the state to 1e-12 plus 1e-5 of
  each value's size. A pulse of vanishing length equals an
  instantaneous rotation to 1e-5.
- With the measurement switched on and nothing else, the variance of
  `J_z` follows the closed-form result `(N/4)/(1 + Gamma_m N t/4)` to
  1e-5 relative, and the uncertainty product `Var(J_y) Var(J_z)` stays
  at its smallest allowed value, `(N/4)^2`, to 1e-6 relative.
- The mean-field equations (averages only, no correlations) match the
  pendulum equations of the paper to 0.2 %.
- A leading-order formula for the growth rate of the `J_z` variance
  under collective emission (Sec. S7 of the paper) matches an exact
  6-spin master-equation solution to 0.2 %. This test uses QuTiP only;
  it checks the formula, not the solver.

**The independent trajectory solver**

- For 200 identical spins, its best Wineland parameter matches the
  exact one-axis-twisting result to 0.05 dB (4000 trajectories, seed
  1). This tolerance is smaller than the random scatter between
  seeds: with seeds 1 to 4 the difference is +0.03, -0.13, +0.03 and
  -0.22 dB, so the test passes for seed 1 but would not for every
  seed.
- Making its time step four times smaller (200 to 800 steps) changes
  the answer by less than 0.01 dB.
- Without interaction, a Gaussian line dephases as `exp(-t^2/2)` to
  0.02.
- For 4 spins with different detunings, its mean spin follows exact
  free precession when there is no interaction, and agrees with the
  cumulant solver when there is, to 0.08 in spin units (4000
  trajectories). This pins the sign of `J_y` (new in 1.15.1).
- The moments are returned at exactly the requested times: for spins
  with one common detuning and no interaction, the mean spin and the
  covariance at each time equal those at `t = 0` rotated by
  `delta t`, to 1e-9 (new in 1.16; 1.15.1 is off by 1.3e-4 here).

**Metrology and one-axis twisting**

- A coherent spin state gives `xi2_S = xi2_R = 1`, contrast 1 and
  `dphi = 1/sqrt(N)`.
- On a state squeezed by the solver, `xi2_R = xi2_S / contrast^2` to
  1e-10 and equals the solver's own Wineland parameter to 1e-12.
- The clock and magnetometer formulas follow the expected scalings
  with `tau`, atom number and dead time.
- `oat_closed_form` against brute-force exact evolution for N = 2, 5
  and 8 at random twisting angles: the mean spin, and the variance at
  the predicted best angle, to 1e-12; the smallest and largest
  variance over a 2001-point angle grid to 1e-4. `mu = 0` gives exactly
  `xi2_S = 1`.

**Dick effect**

- The closed-form Fourier coefficients equal independent adaptive
  numerical quadrature to 1e-9 (three cycle layouts, with and without
  finite pulses).
- Continuous interrogation (no dead time, no pulses) gives every
  coefficient below 1e-12 and a Dick deviation below 1e-18.
- Without pulses, the size of each coefficient divided by the zeroth
  one equals `|sin(pi m eta)/(pi m eta)|` to 1e-12, where `eta` is the
  Ramsey time divided by the cycle time.
- The floor scales as `1/sqrt(tau)` to 1e-12; doubling the harmonic
  cutoff changes it by less than 1e-3; too low a cutoff is refused.
- For white frequency noise `S_y = h0` and instantaneous pulses, the
  Allan variance equals the closed form `h0 (1 - eta) / (2 eta tau)`
  (from the Fourier series `sum sin^2(m x)/m^2 = x (pi - x)/2`) to
  6e-5 for `eta` = 0.1, 0.4 and 0.9, and the small shortfall equals the
  harmonics left out above the default cutoff, about
  `1 / (pi^2 eta (1 - eta) m_max)`, to 1e-6 of itself (new in 1.17).
- With a noisy laser, about 10 dB of squeezing improves the total by
  less than 10 %; with a quiet laser, by more than 3 times.
- `synchronized_comparison`: for two equal ensembles the difference is
  sqrt(2) times one clock (to 1e-18 relative) and `per_clock` equals
  one clock (to 1e-15 relative); the single-ensemble values equal
  `clock_allan_deviation` exactly.

**Measured data and planning**

- A noise-free tomography record returns the covariance matrix to
  1e-9, and the smallest and largest variances equal NumPy's
  eigenvalues to 1e-9.
- Shots drawn from the exact one-axis-twisting state return its
  `xi2_R` within 4 reported error bars.
- Coherent-spin-state shots return `xi2_R` compatible with 1 (within 4
  error bars). Through the full file-to-estimate path, with a fitted
  contrast, the test allows 4 error bars plus 0.1 (and the contrast
  4 error bars plus 0.02 from 1).
- Over 150 (squeezing) and 60 (contrast) seeded repeats, the scatter of
  the estimates is between 0.6 and 1.6 times (squeezing) and 0.5 and
  2 times (contrast) the reported error bar.
- A clean fringe returns contrast and phase to 1e-7.
- Shot files round-trip exactly.
- The formula for the scatter of a sample variance,
  `sigma^4 [kurtosis/M - (M-3)/(M(M-1))]`, equals the exact value
  obtained by listing every possible outcome of a three-valued
  distribution (2, 4 and 5 shots), to 1e-12 relative.
- With weights that are not the true errors, the fit covariance
  reported with `error_sigmas` matches 3000 simulated fits to 10 %;
  with the true errors as weights it equals the usual matrix to 1e-12.
- Over 400 seeded experiments (9 angles, 20 shots each), the fitted
  mean variance `(V1 + V2)/2` is 14 % low with `weighting="sample"`
  and within 3 standard errors of the truth with `weighting="model"`.
  The model-weighted covariance equals `plan_tomography` at the
  fitted values to 1e-9.
- With 3 shots per angle (200 seeded datasets), whenever the default
  falls back it returns exactly the `weighting="sample"` result with
  one `RuntimeWarning`, and it refuses only datasets that
  `weighting="sample"` also refuses.
- With 5 % outlier shots, over 200 experiments the Gaussian error bar
  is too small by more than 1.5 times, and the empirical one is within
  0.85 and 1.2 times the scatter.
- For a coherent state, the low bias of the fitted smallest variance
  (see [Limits](#limits)) is between 0.85 and 1.3 reported error bars
  over 400 experiments; the expected value for many shots is
  `sqrt(pi/3) = 1.02`.
- Fit quality (new in 1.17): the reported `chi2` equals the weighted
  sum of squared residuals recomputed in the test from the shots (to
  1e-9), for the tomography with both weightings and for the fringe.
  The general form used with `shot_noise="empirical"` equals the usual
  sum when the weights are the errors (to 1e-10), and for Gaussian data
  with other weights it follows the chi-square law (2000 seeded fits:
  mean within 3 standard errors of the degrees of freedom, p-values
  uniform by a Kolmogorov-Smirnov test). Over 400 seeded experiments
  with 200 shots per angle (both `shot_noise` options) and 400 fringe
  experiments, between 2 % and 9 % give `p < 0.05`. Rotation angles or
  fringe phases 5 % larger than recorded give `p < 0.01` in each of 20
  experiments.

**Decay fits** (new in 1.17)

- Noise-free echo, Gaussian-line and Lorentzian-line decays built from
  the formulas are recovered to 1e-9.
- The echo decay computed by the cumulant solver (a 2 kHz Gaussian
  line in 16 classes, no interaction) gives back the `T2` given to
  `from_hz` to 1e-7. The Ramsey decay computed by the mean-field solver
  on a Gaussian line cut into 20000 classes gives back the FWHM given
  to `lineshape` and the `T2` to 5e-5 (the limit set by the 20000
  classes). This checks that the fitted numbers mean what the model's
  inputs mean.
- Over 300 seeded noisy datasets, the scatter of the fitted FWHM and
  `1/T2` is between 0.8 and 1.2 times the reported error bar, with no
  bias beyond 3 standard errors; the same holds for the Lorentzian
  width with a known `T2`. The mean `chi2` is within 3 standard errors
  of its degrees of freedom and 2 % to 9 % of the p-values are below
  0.05.
- With `T2` known only to +/- 0.2 ms, passing `T2_sigma` gives a width
  error bar within 0.8 to 1.2 times the scatter; leaving it out makes
  the error bar more than 1.3 times too small. The propagated term
  equals a finite-difference refit to 2 %.
- A Gaussian-line Ramsey decay fitted without the line gives
  `p < 1e-6`.

**Line shapes and inputs**

- The Voigt line's total FWHM, checked by computing the profile as a
  direct numerical convolution of its Gaussian and Lorentzian parts
  (independently of SciPy's Voigt function), equals the requested
  width: the profile at half the width is half its peak to 1e-9, for
  Lorentzian fractions 0.1, 0.3, 0.6 and 0.9. At the same fractions the
  old split misses by between 1e-6 and 2.5e-4. Its tabulated quantiles
  agree with numerical integration of the profile to 1e-7, and its
  pure-Lorentzian limit with the exact Lorentzian quantiles to 1e-6.
- Each refusal added in 1.16 is triggered in the tests, and the
  accepted edge cases (no dephasing, zero temperature, a lossless
  cavity, a negative gyromagnetic ratio) still work.
- The planned covariance equals the fitted one to 1e-12, equally
  spaced angles give the diagonal design to 1e-12, and over 300 seeded
  experiments the scatter matches the planned error bar to 20 %.
- `bosonic_mode` maps a coherent spin state to the vacuum, and its
  squeezing agrees with `squeezing_parameters` to 1e-10; the QuTiP
  export reproduces the covariance at 25 angles to 1e-4; a deliberately
  small Fock cutoff is caught.

Each input that 1.16 accepted with a wrong or `nan` answer and 1.17
refuses has a test (`tests/test_fixes_117.py`) that fails on 1.16.0.
The remaining tests check refusals, further identities of the same
functions, properties of the paper's pendulum model (the locking
condition and a conserved quantity), and that the version number is
the same in the package, its metadata and CITATION.cff.

## Corrections in earlier versions

**1.17.0 refuses inputs that gave wrong numbers without an error.**
No result of a valid input changes. The inputs that are now refused:

- `oat_closed_form` truncated a number of spins that was not a whole
  number: `oat_closed_form(2.9, 0.3)` returned the result for 2 spins
  (`xi2_S = 0.8506`), and `oat_closed_form(100.7, 0.05)` the one for
  100.
- `plan_tomography` truncated a fractional shot count: 9.9 shots per
  angle were planned as 9 (`var_min_sigma` 28.72 in the test design,
  where 10 shots give 27.08).
- `evolve_meanfield` integrated backwards for a negative time: with
  `T2 = 1 ms` and `t = -1 ms` it returned a contrast of 2.718 (a spin
  longer than 1). `evolve` already refused this.
- A `nan` or infinite `contrast_sigma` gave a `nan` or infinite
  `xi2_R_sigma`; a `nan` or infinite `T_cycle`, or a `nan` Ramsey or
  pulse time, gave a `nan` or infinite clock instability or Dick floor;
  `nan` or infinite laser-noise coefficients were accepted;
  `log_uniform_weights` returned empty arrays for 0 weights and `nan`
  weights for an infinite range.

Docstrings that disagreed with the code were corrected: the fine scan
of `optimal_squeezing` runs between the neighbors of the best coarse
point (the docstring said `[t*/2, 2 t*]`); the gain formula of
`plain_squeezed_readout` includes the contrast,
`1/(xi2 + 4 sigma_det^2 / (S2 C^2))`; and the harmonic-cutoff check of
`dick_allan_deviation` looks at the top 10 % of the harmonics (it said
"last decade") and bounds the omitted part only roughly (see its
docstring).

**1.16.0 changed the tomography fit weights.** Each angle's sample
variance is weighted by how precisely it is known, and that precision
is proportional to the variance itself. Up to 1.15.1 the variance in
the weight was the measured one, so an angle whose variance came out
low by chance got a larger weight, and the fit leaned low. With few
shots per angle this is large: over 400 simulated experiments with 9
angles and 20 shots each, the fitted mean variance was 14 % low. Since
1.16.0 the weight uses the fitted curve, and the fit is repeated until
it no longer changes; in the same test the mean is off by
-0.2 % +/- 0.6 %, consistent with no bias. With many shots the
difference is small (example 5 moved from `xi2_R = 0.1421` to
`0.1420`). `weighting="sample"` gives the old result, and it is also
what the default falls back to (with a warning) when there are too few
shots per angle for the new weights.

**1.16.0 returns trajectory-solver moments at the requested times.**
`cavsqueeze.dtwa.evolve` stepped on one uniform grid and returned the
moments at the grid point nearest each requested time, up to half a
step away (for example `t = 0.37` was read at `0.370075`). Each
requested time is now reached exactly. In the one-axis-twisting test
the best squeezing moves from -14.1731 dB to -14.1736 dB (exact:
-14.2020 dB).

**1.16.0 made the Voigt width match the request.** The Gaussian and
Lorentzian parts of a Voigt line were split with the
Olivero-Longbothum approximation, which gives a total FWHM off by up
to about 2.3e-4 of the requested width. The split is now solved for
numerically, and the width is checked to 1e-9. For
`lorentz_fraction = 0.3` the class detunings of a 16-class
discretization move by at most 2.0e-4 of their value.

**1.15.1 fixed a sign in the trajectory solver.** `cavsqueeze.dtwa`
integrates the same equations as the cumulant solver, whose variable
is `<sigma^+> = (x + i y)/2`, but read the results out as if the
variable were `<sigma^->`. The `J_y` component of the mean spin, and
the `xy` and `yz` covariances, therefore came out with the wrong sign
(in the 4-spin case of the new test at `t = 1`, `<J_y> = -0.34`
where the cumulant solver gives `+0.33`).
Squeezing parameters did not change, because they do not depend on
that sign. The read-out is fixed, seeded runs give the same squeezing
parameters as before, and a new test checks the sign.

**Corrected claims (documentation only, 1.15.1).**

- Earlier documents said the Dick coefficients gave "exactly zero"
  noise for continuous interrogation and the "exact"
  rectangular-window ratio; the tests check these to 1e-12. The
  module docstring of `cavsqueeze.dick` also said the coefficients
  were checked against an FFT; the test uses adaptive quadrature.
- The 1.10.0 notes said every quantity returned by `oat_closed_form`
  is checked against exact evolution to 1e-12. The mean spin and the
  smallest variance are; the largest variance is checked to 1e-4 on
  an angle grid, and `xi2_S` and `xi2_R` are computed from these,
  not compared separately.
- "Machine precision" and "exact identity" for the planner mean
  agreement to 1e-12.
- The `dtwa` module said the cumulant solver was checked against 8
  distinguishable spins; the tests use 2 and 4.

The full history is in [CHANGELOG.md](CHANGELOG.md).

## Limits

- The cavity is eliminated: its field follows the spins instantly,
  which holds when the cavity is far detuned or much faster than the
  spin dynamics.
- The solver keeps averages and pair correlations only. Tests bound
  the resulting error for small systems (above); for strongly
  non-Gaussian states it grows.
- The trajectory solver covers the twisting and the detunings, not
  decay or thermal photons.
- The continuous-measurement option does not condition the spectator
  spins of `tail_resolved_classes`; use a discretization without
  spectators with it.
- `fit_coherence_decay` knows only the decays of this model: an
  exponential echo decay and a Gaussian or Lorentzian line. Echo decays
  in solids are often not exponential (spectral diffusion); a small
  p-value says so, but the package has no model for it. A Voigt line
  is not fitted. The fit's error bars are the standard linear
  approximation, checked in the tests for decays sampled well past
  their half point.
- The `chi2` tests of the tomography and fringe fits assume nearly
  Gaussian sample variances and means, which needs many shots per
  angle (the tests use 200 to 2000). With `shot_noise="empirical"` and
  few shots, small p-values are more common than stated.
- By default the error bars of `estimate_squeezing` assume Gaussian
  shot noise, and heavy-tailed detection noise makes them too small;
  use `shot_noise="empirical"`. That option estimates the kurtosis of
  each angle's shots from the shots themselves, which needs enough
  shots to see the tails (the tests use 400 per angle). The planning
  functions (`plan_tomography`, `shots_for_squeezing`) still assume
  Gaussian shots.
- The fitted smallest variance is the smaller of two noisy numbers and
  is biased low. For a coherent state (two equal variances) measured at
  equally spaced angles the bias is about one reported error bar, so
  such a state can look slightly squeezed. This bias is not
  corrected; compare with the error bar before claiming squeezing of
  about that size.
- The clock and magnetometer formulas use linear error propagation
  around the operating point (the same approximation that defines
  `xi2_R`).
- The Voigt line's quantiles come from a tabulated cumulative
  distribution (checked to 1e-7 against numerical integration).
- The bosonic export keeps second moments only.

## Where it comes from

The package was developed for 171Yb3+:CaWO4 and applies to any spin
ensemble coupled to a resonator. The physics, the conventions and the
validation are described in: T. M. Mahim, M. M. Rahman, A. S. M.
Mohsin, *Synchronization sets the coherence and the squeezing limit of
a spin ensemble in a cavity*. The paper's companion repository,
[yb-cawo4-cavity-squeezing](https://github.com/Tanvir-Mahmud-Mahim/yb-cawo4-cavity-squeezing),
contains the scripts, datasets and figures that reproduce the paper and
is archived on Zenodo; this repository is the software's home for
development, releases and support.

Other sources named in the code: Kitagawa and Ueda, Phys. Rev. A 47,
5138 (1993); Wineland et al., Phys. Rev. A 46, R6797 (1992); Sorensen
et al., Nature 409, 63 (2001); Itano et al., Phys. Rev. A 47, 3554
(1993); Ludlow et al., Rev. Mod. Phys. 87, 637 (2015); G. J. Dick,
Proc. PTTI 1987; G. Santarelli et al., IEEE Trans. Ultrason.
Ferroelectr. Freq. Control 45, 887 (1998); A. Quessada et al., J. Opt.
B 5, S150 (2003); Schulte et al., Nat. Commun. 11, 5955 (2020).

## Citing, support and license

Please cite the paper if you use this code. To cite the software, use
the concept DOI
[10.5281/zenodo.22278034](https://doi.org/10.5281/zenodo.22278034),
which always resolves to the latest release. Citation metadata is in
[CITATION.cff](CITATION.cff).

The package is maintained by Tanvir Mahmud Mahim (BRAC University),
who reviews issues and pull requests; the scientific content is
developed with the co-authors of the associated paper. Bug reports,
questions and pull requests are welcome through
[GitHub issues](https://github.com/TaN-MM-Org/cavsqueeze/issues); see
[CONTRIBUTING.md](CONTRIBUTING.md) for the development setup and the
testing requirements. Tagged releases are published to PyPI by CI.

Licensed under Apache-2.0 (see [LICENSE](LICENSE)).
