Metadata-Version: 2.4
Name: cavsqueeze
Version: 1.16.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?
- 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.16.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.

## 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.
- `plan_tomography`, `shots_for_squeezing` -- predicted error bars and
  the shot count for a target.

**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);
- 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;
- a clock or magnetometer cycle time `T_cycle` is shorter than the
  Ramsey time;
- the Ramsey timing is impossible (`T_ramsey <= 0`, a negative pulse
  length, or a cycle shorter than the Ramsey time plus two pulses), or
  fewer than 1 harmonic is asked for;
- laser noise coefficients are negative;
- 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` is not a whole number of at least 2;
- the contrast given to `estimate_squeezing` or `plan_tomography` is
  not in (0, 1], or its error is negative;
- 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`, or a
  target error bar is not positive, uses a degenerate angle set, or is
  below the contrast floor (example 6).

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

123 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.
- 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`.

**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.

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.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.
- 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)).
