cavsqueeze

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? Is my echo decay exponential, as the model assumes? And do my measured data actually follow the laws the fits assume?
• Does my detector see the projection noise of the spins, and how much noise does it add? How large is the squeezing once that noise is removed, with the error of the calibration included?
• 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). 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
• Install, requirements and units
• Examples (each with the output it prints)
• Tutorials (three scripts that go from data to a result)
• What is in the package
• When it refuses, and why
• How the results are checked
• Corrections in earlier versions
• Limits
• Where it comes from
• 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.19.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?

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

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

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

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

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

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

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

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

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

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.

Since 1.18 the fit also takes line="voigt" (a Gaussian and a Lorentzian part together, with T2 from the echo), and fit_stretched_decay tests whether the echo decay is exponential. The tutorial examples/tutorial_decays_to_model.py uses both (see Tutorials).

11. Do the tomography data follow the rotation law?

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

12. Does the detector see projection noise?

import numpy as np
from cavsqueeze import calibrate_projection_noise

# Synthetic, seeded calibration: coherent states of N spins (binomial shots)
# read by a detector with noise variance 9, at eight atom numbers.  In the
# second record every shot's pi/2 pulse also jitters by 0.02 rad (rms), which
# moves all spins together and adds noise that grows like N^2.
rng = np.random.default_rng(5)
N = np.arange(500.0, 4001.0, 500.0)
for jitter in (0.0, 0.02):
    shots = []
    for n in N:
        tip = rng.normal(0, jitter, 2000)                # common tipping angle
        up = rng.random((2000, int(n))) < (1 + np.sin(tip))[:, None] / 2
        shots.append(up.sum(axis=1) - n / 2 + rng.normal(0, 3.0, 2000))
    line = calibrate_projection_noise(N, shots)
    quad = calibrate_projection_noise(N, shots, quadratic=True)
    print(f"jitter {jitter} rad: line p = {line.p_value:.2g}, "
          f"4 c1 = {line.projection_ratio:.3f}; parabola p = {quad.p_value:.2f}, "
          f"4 c1 = {quad.projection_ratio:.2f} +/- {quad.projection_ratio_sigma:.2f}, "
          f"c2 = {quad.technical_coefficient:.1e} +/- {quad.technical_coefficient_sigma:.1e}")

jitter 0.0 rad: line p = 0.46, 4 c1 = 1.010; parabola p = 0.39, 4 c1 = 0.96 +/- 0.07, c2 = 3.3e-06 +/- 4.6e-06
jitter 0.02 rad: line p = 1.7e-22, 4 c1 = 2.297; parabola p = 0.55, 4 c1 = 1.12 +/- 0.11, c2 = 9.0e-05 +/- 8.7e-06

Before a squeezing measurement, an experiment checks its readout on coherent spin states. For N independent spins, each up or down with probability 1/2, J_z is a binomial count minus N/2, so its variance is N/4 (projection noise; W. M. Itano et al., Phys. Rev. A 47, 3554 (1993)). A detector adds a constant variance c0. Noise that moves all spins together adds a term in N^2: if every spin is tipped by the same random angle, the variance is (N/4) E[cos^2] + (N^2/4) Var(sin) of that angle, exactly. calibrate_projection_noise (new in 1.19) fits c0 + c1 N (or, with quadratic=True, c0 + c1 N + c2 N^2) to the measured variances and reports projection_ratio = 4 c1 (1 for pure projection noise) with its error bar and the p-value p_projection of c1 = 1/4, the detector variance detection_variance = c0 with its error bar, and chi2, dof, p_value for the fit.

Here the truth is 4 c1 = 1.000 and c2 = 0 without jitter, and 4 c1 = 1.000 and c2 = 1.0e-4 with it (from the formula above). Without jitter the straight line fits. With jitter the straight line is rejected (p = 1.7e-22) and its slope is wrong by more than a factor 2; the parabola fits, and its two numbers are within about one error bar of the truth. The fit weights each point by the error of its sample variance, read from the fitted curve and iterated, as in estimate_squeezing; shot_noise="empirical" uses each point's measured kurtosis instead of assuming Gaussian shots.

The detector variance and its error bar go into the squeezing estimate: estimate_squeezing(..., detection_variance=c0, detection_variance_sigma=...) subtracts it and adds its error to that of var_min, xi2_S and xi2_R (new in 1.19). The tutorial examples/tutorial_detection_noise.py does both steps (see Tutorials).

Tutorials

Three scripts in the examples/ folder of the source repository (and, since 1.19, of the source distribution on PyPI) go from data to a final result with error bars. Each uses only the package and seeded synthetic "measurements" drawn from known values, so the result can be compared with the truth. With your own data, replace the synthetic block with your measurements. Each runs in a few seconds, and a test checks that it still prints the numbers below.

From decays to a squeezing prediction (python examples/tutorial_decays_to_model.py). An echo decay gives T2, and fit_stretched_decay checks that the echo is exponential (stretch exponent x = 1, as the model assumes). A Ramsey decay, fitted with the echo's T2 and its error bar, gives a Voigt line: its total FWHM and Lorentzian fraction (truth: 150 Hz and 0.3, T2 = 5 ms). The fitted values go into the model, which predicts the squeezing after a 1 ms echo twist of 10^12 spins. The error bar of the prediction comes from changing each fitted input by its error bar.

echo:     T2 = 5.028 +/- 0.038 ms, exponent x = 1.002 +/- 0.016, p(x = 1) = 0.89
Ramsey:   FWHM = 148.8 +/- 1.9 Hz, Lorentzian fraction = 0.315 +/- 0.026, p = 0.61
model:    xi2_R after 1.0 ms = -4.55 +/- 0.06 dB (true line and T2: -4.56 dB)

From shot records to a clock's stability (python examples/tutorial_tomography_to_clock.py). A fringe file and a tomography file of a squeezed state of 1000 spins (truth: contrast 0.9989, xi2_R = 0.2508, that is 6.01 dB) are loaded, fitted and checked with their chi-square tests. The squeezing then sets the projection noise of an illustrative clock (429 THz, 0.4 s Ramsey time in a 1 s cycle), shown next to the Dick floor of a laser with white frequency noise h0 = 1e-35 (a placeholder).

contrast: 0.9988 +/- 0.0007 (fringe fit chi2 = 8.8 for 13 dof, p = 0.79)
xi2_R:    0.2473 +/- 0.0096 = 6.07 +/- 0.17 dB (tomography chi2 = 3.5 for 9 dof, p = 0.94)
clock:    projection 1.459e-17 +/- 2.8e-19 (no squeezing 2.933e-17), Dick 2.739e-18, total 1.484e-17 at 1 s

From a detector calibration to corrected squeezing (new in 1.19; python examples/tutorial_detection_noise.py). Coherent states at six atom numbers (binomial shots) are read by a detector that adds noise of variance 25 (spin units squared). calibrate_projection_noise checks the projection noise (4 c1 = 1) and measures the detector variance with its error bar. The tomography of a squeezed state of 1000 spins (truth: xi2_R = 0.2508), read by the same detector, is then fitted twice: without the correction the detector noise makes the state look less squeezed (0.343); with it the estimate agrees with the truth, and its error bar includes the error of the calibration. The last line plans the next run: the errors of the contrast and of the calibration set a floor under the error bar of xi2_R that more tomography shots cannot lower, and shots_for_squeezing gives the shots for a target above it.

calibration: detection variance 25.3 +/- 3.2, 4 c1 = 1.014 +/- 0.027 (p = 0.60), fit p = 0.89
squeezing:   xi2_R = 0.3430 +/- 0.0088 raw, 0.2417 +/- 0.0154 = 6.17 +/- 0.28 dB corrected (true 0.2508)
planning:    error-bar floor 0.0127 from contrast and calibration; +/- 0.015 needs 2396 shots per angle

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 (since 1.18 the times can be given in any order).
• 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). New in 1.18: line="voigt" fits a Voigt line and returns its total FWHM and lorentz_fraction, the two numbers lineshape("voigt", ...) takes (also fwhm_gaussian and fwhm_lorentzian).
• fit_stretched_decay (returns StretchedDecayFit, new in 1.18): fits an echo decay with A exp(-(t/T_M)^x) and tests x = 1, the exponential decay of the model (p_exponential).
• calibrate_projection_noise (returns ProjectionNoiseFit, new in 1.19): the readout noise of coherent spin states against the atom number, c0 + c1 N (or + c2 N^2), with the test of projection noise c1 = 1/4 and the detector variance c0 with its error bar (example 12).
• detection_variance_sigma (new in 1.19) in estimate_squeezing, plan_tomography and shots_for_squeezing: the error bar of the subtracted detector variance, added to that of var_min in quadrature. Subtracting a constant moves the fitted var_min by that constant and leaves the fit weights unchanged, so this needs no approximation.

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 (since 1.18, at least the top 100) 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.

New in 1.18:

• the output times t_eval of evolve and evolve_meanfield lie outside [0, t] (before, evolve with t = 0 returned the starting state for every requested time), or the times given to ramsey_cumulant and ramsey_meanfield are negative or not finite;
• save_shots_csv gets an angle with fewer than 2 shots, a non-finite angle or non-finite shots (before, an angle without shots was left out of the file without a word, and the others gave files that load_shots_csv could not read);
• sample_variance_sigma gets a kurtosis below 1 or not finite (no distribution has a kurtosis below 1; before, the result was nan);
• dick_allan_deviation gets a spectrum with nan or negative values at the harmonics (before, the floor came back as nan), a tail_tol that is not finite and positive, or dick_fourier_coefficients a cutoff that is not a whole number (2.7 used to mean 2). The cutoff check now looks at no fewer than the top 100 harmonics (the top half when there are fewer than 200); before, a window of one harmonic could fall on a zero of the coefficients and pass a result 2 % short;
• plan_tomography and shots_for_squeezing get a theta_min that is not finite (before, NumPy's LinAlgError), and cavsqueeze.ensemble.tail_resolved_classes gets class counts below 1 or impossible edges (before, ZeroDivisionError for M_core = 0);
• fit_coherence_decay(line="voigt") gets no T2, or one part of the line (Lorentzian or Gaussian) is not at least min_significance error bars above zero (the message names the simpler line to fit);
• fit_stretched_decay gets fewer than 4 points or 3 distinct times, times that do not sample the decay, or a T_M or exponent that is not at least min_significance error bars above zero.

New in 1.19:

• calibrate_projection_noise gets inputs of different lengths, atom numbers that are not finite and positive, fewer than 2 shots or non-finite shots at an atom number, identical shots, fewer atom numbers than fitted coefficients plus one, fewer distinct atom numbers than coefficients, or data whose fitted variance curve reaches zero or below at a measured atom number (it cannot then set the weights);
• detection_variance_sigma is negative or not finite, or a target of shots_for_squeezing lies below the floor set by the contrast error and the detection-variance error together (the message names that floor);
• pulse and twist_imperfect get a pulse duration that is negative or not finite (before, a negative duration gave an instantaneous pulse without a word), or pulse gets an n_steps that is not a whole number of at least 1 (before, 0 ended in ZeroDivisionError).

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

251 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.
• For many spins (new in 1.18): at the optimal twisting angle the smallest variance agrees with the same formulas evaluated with 60 decimal digits (Python's decimal module) to 1e-9 for N = 10^7, 10^8 and 10^9, and with exact evolution in the (N + 1)-dimensional Dicke basis at N = 10^6 to 1e-8. For N up to 1000 the values differ from those of 1.17 by less than 1e-11 of N/4.

Readout protocols (new in 1.18)

• Without decay, the twisting is a shear of the transverse noise by Q = chi N t. For N = 10^6 and Q = 0.5 and 2, twist_untwist gives the gain 1 + Q^2 (and (1 + Q^2)/(1 + 4 sigma_det^2/N) with detection noise) and the amplification sqrt(1 + Q^2) to 2e-5, and plain_squeezed_readout gives 1/xi2 with xi2 = 1 + Q^2/2 - Q sqrt(1 + Q^2/4) to 2e-5. An ideal twist_imperfect equals twist exactly.

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) 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.
• Voigt line (new in 1.18): decays built from the product formula are recovered to 1e-9 (FWHM, lorentz_fraction and both parts, for fractions 0.1, 0.3 and 0.6), and decays computed by numerical Fourier integration of the profile of lineshape("voigt", ...) give back its FWHM and fraction to 1e-7. Over 300 seeded datasets, in which the given T2 is itself drawn with its error bar, the scatter of the FWHM and of the fraction is between 0.8 and 1.2 times the reported error bars (1.00 and 1.00 in the test run).
• Stretched echo (new in 1.18): noise-free decays with x = 0.7, 1, 1.5 and 2.4 are recovered to 1e-9. The exponential reference chi2_exponential equals a separate fit by SciPy's curve_fit to 1e-7. For exponential data (200 seeds) the scatter of T_M and x is 0.8 to 1.2 times the error bars, and 1 % to 10 % of the datasets give p_exponential < 0.05 (3.5 % in the test run; the p-values pass a Kolmogorov-Smirnov test for uniformity). For x = 1.8 (100 seeds) every dataset rejects x = 1.

Projection-noise calibration (new in 1.19)

• Variances that follow c0 + c1 N (+ c2 N^2) exactly give back the coefficients to 1e-9, for atom numbers up to 10^6.
• At the weights the fit converges to, the coefficients and their covariance equal those of NumPy's polyfit (a separate implementation) to 1e-9 of the error bars, and the weights equal the fitted curve times sqrt(2/(M-1)).
• Binomial shots (pure projection noise) plus detection noise of variance 16, 200 seeded calibrations: no bias beyond 3 standard errors, the scatter of c0 and c1 is 0.8 to 1.2 times the reported error bars, and between 1 and 22 of the 200 tests of c1 = 1/4 and of the fit give p < 0.05 (10 expected).
• Shots tipped by a common random angle (0.02 rad rms): the parabola gives back c1 = (1 + e)/8 and c2 = (1 - e)/8, e = exp(-2 s^2), derived in the test, within 4 error bars, and the straight line is rejected (p < 1e-6).
• With 5 % of detector shots at 4 times the width, the Gaussian error bar of c0 is more than 1.3 times too small over 200 seeds, and the empirical one is within 0.8 to 1.2 times the scatter.
• The whole chain (calibration, then a tomography fit with detection_variance and detection_variance_sigma), 200 seeds: the scatter of xi2_S is 0.8 to 1.2 times the reported error bar, and without detection_variance_sigma the error bar is more than 1.5 times too small. On data, subtracting the detector variance moves var_min by that constant (to rounding) and leaves its error bar and the fit covariance unchanged, for both weightings and both shot_noise options; the planned error bar equals the estimate's to 1e-9.

Speed-ups (1.18 and 1.19)

• The 1.17 code of each routine made faster in 1.18 is kept in tests/test_speed_118.py as the reference: the Dick Fourier coefficients for nine cycle layouts (20000 harmonics, including exact resonances) agree to 1e-12 of the largest coefficient, the solver's right-hand side for 1, 3 and 16 classes to 1e-12 of its largest entry, and a whole solver run to 1e-9; optimal_squeezing gives exactly the result of a separate squeezing_after call at its best time. (Until 1.19 these tests demanded bit-identical results. The vectorized Dick coefficients call NumPy's array functions where the 1.17 code called scalar ones, whose last digits can differ between processors, so the tests now use these tolerances.)
• The 1.18 right-hand side is kept in tests/test_speed_119.py as the reference for the 1.19 one: for random states with 1, 3 and 16 classes, with and without the correlation feedback and with the continuous measurement, they agree to 1e-12 of the largest entry (5e-16 measured), and whole runs of README examples 1 and 2 agree to 1e-9 of the largest entry of the state (8e-14 measured) and to 1e-6 in xi2 (1e-9 measured).

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. Each bug fixed in 1.18 has a test (tests/test_fixes_118.py) that fails on 1.17.0, and each one fixed in 1.19 a test (tests/test_fixes_119.py) that fails on 1.18.0. The tutorials are run by tests/test_tutorials.py, which checks every printed number to its last digit. 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.19.0 fixes two input problems.

• The 1.18 notes said ramsey_cumulant and ramsey_meanfield take the times in any order. A time that appeared twice, even in a sorted list such as [1e-3, 1e-3], ended in SciPy's error "Values in t_eval are not properly sorted". evolve and evolve_meanfield now accept output times in any order and with repeats, and return each in the order given (for T2 = 10 ms and no coupling, [2e-3, 1e-3, 2e-3, 0.0] gives the contrasts 0.818731, 0.904837, 0.818731, 1). Results for strictly increasing times are unchanged (the same solver call).
• pulse with a negative duration returned the instantaneous pulse without an error (the same state as duration=0); it is now refused, as are a non-finite duration and an n_steps that is not a whole number of at least 1 (0 used to end in ZeroDivisionError).

A docstring that disagreed with the code was corrected: CavityParams.chi and CavityParams.Gamma_SR include g^2 (they are per unit product of the coupling weights, not "per unit g_j g_k").

1.18.0 fixes wrong numbers that came without an error.

• oat_closed_form lost its accuracy for many spins. The formulas raise cosines close to 1 to the power N, which multiplies their rounding error by about N. At the optimal twisting angle the smallest variance was off by -4.0e-4 of its value for N = 10^7 and by +3.9 % for N = 10^8, and for N = 10^9 it came out negative (-2251.75; the true value is 275.16). It is now computed in a form without that loss (see How the results are checked).
• ramsey_cumulant integrated up to the last time in the list, not the largest. With the times [1e-3, 0.0] it returned the contrast 1.0 at 1 ms, where the true value is 0.9048. The root cause, now refused: evolve with t = 0 returned the starting state for any requested output time. The Ramsey functions now accept the times in any order.
• ramsey_meanfield and evolve_meanfield failed with a TypeError for t = 0.
• save_shots_csv left an angle without shots out of the file without a word (3 angles saved, 2 loaded back).
• sample_variance_sigma returned nan for a kurtosis below 1, and dick_allan_deviation returned nan for a spectrum with nan or negative values.
• The harmonic-cutoff check of dick_allan_deviation could look at a single harmonic: with m_max = 10 and a Ramsey time of 0.4 of the cycle that harmonic is exactly zero, and the result (8.476e-16 for h0 = 1e-30) was 2.1 % below the converged 8.660e-16 without an error. For m_max of 1000 or more the check is unchanged.

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.

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, Lorentzian or (since 1.18) Voigt line. Echo decays in solids are often not exponential (spectral diffusion). fit_stretched_decay (new in 1.18) measures the stretch exponent and tests the exponential form, but the solver itself has only exponential (Markovian) dephasing: when the test fails, T2 is not a parameter that describes the measured decay. The fits' 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.
• calibrate_projection_noise treats the atom numbers as exact. An error in N that is common to all points scales c1 (and so projection_ratio), and point-to-point errors in N are not included in the error bars. Its p-values, like those of the tomography fit, assume nearly Gaussian sample variances (the tests use 1500 or more shots per point).

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, 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). The stretched echo form exp(-(2 tau/T_M)^x) is the one used, for example, by M. Le Dantec et al., arXiv:2106.14974 (2021), who cite T. Bottger, C. W. Thiel, Y. Sun and R. L. Cone, Phys. Rev. B 73, 075101 (2006).

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, which always resolves to the latest release. Citation metadata is in 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; see 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).