Metadata-Version: 2.4
Name: vacspin
Version: 0.4.0
Summary: Group-IV colour-centre spin-photon interfaces: strained spin Hamiltonians, cyclicity, transition-resolved Purcell budgets and exact single-shot readout statistics
Author: Tanvir Mahmud Mahim
License: Apache-2.0
Project-URL: Homepage, https://github.com/TaN-MM-Org/vacspin
Project-URL: Issues, https://github.com/TaN-MM-Org/vacspin/issues
Project-URL: Changelog, https://github.com/TaN-MM-Org/vacspin/blob/main/CHANGELOG.md
Keywords: quantum,color-center,SnV,SiV,diamond,cavity-QED,Purcell,spin-photon-interface,single-shot-readout
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.22
Provides-Extra: test
Requires-Dist: pytest>=7; extra == "test"
Dynamic: license-file

# vacspin

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

`vacspin` is a Python package for working out how well a single
electron spin trapped in diamond can be read out with light and linked
to another spin. The spin sits in a **group-IV vacancy colour centre**:
a tin, silicon, germanium or lead atom next to a missing carbon atom
(written SnV-, SiV-, GeV-, PbV-). Such a centre can store one quantum
bit and send it out as a photon, which is why it is studied as a node
of a quantum network.

Starting from the numbers that describe one centre, and the cavity and
optics around it, the package answers questions such as:

- At a given magnetic field, what is the qubit frequency, and how many
  photons can the spin emit before it flips by accident?
- How much does an optical cavity speed up the emission, and what
  fraction of the photons reach the detector?
- With that detector, how reliably can one measurement tell the two
  spin states apart, and in how short a time window?
- What detection efficiency, or what window, does a target reliability
  need?
- How often would two such nodes become entangled?
- From your own measured frequencies, what are the parameters of *your*
  centre, with error bars -- and will a planned set of measurements be
  able to determine them at all?
- With those error bars, how uncertain is any number predicted from
  the fitted centre (a qubit frequency at a new field, a cyclicity)?

The whole chain runs on NumPy alone. Results are checked by automated
tests against exact formulas, published measurements or a second
calculation done a different way (see
[How the results are checked](#how-the-results-are-checked), which
also gives each tolerance). When the physics behind a formula no longer
applies, or a question cannot be answered from the data given, 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, units and conventions](#install-units-and-conventions)
- [Examples](#examples) (each with the output it prints)
- [Cited parameter sets](#cited-parameter-sets)
- [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

- **Ground and excited manifold** -- the centre has a group of four
  low-energy states (ground) and four states one optical photon higher
  (excited). In each group the four states form two pairs.
- **Kramers doublet** -- one such pair. At zero magnetic field the two
  states of a pair have exactly the same energy; a field splits them.
  The lowest ground pair is the **qubit**, and its splitting is the
  **qubit frequency**.
- **Zero-field splitting** -- the energy gap between the two pairs of
  a manifold at zero field. It comes from **spin-orbit coupling**
  (`lam`, a property of the atom) and **strain** (`ups`, how the
  crystal around the centre is squeezed; it varies from sample to
  sample). The two combine as `sqrt(lam^2 + 4 ups^2)`.
- **Spin Hamiltonian** -- the 4 x 4 matrix whose eigenvalues are the
  energies of one manifold. `vacspin` builds it from spin-orbit
  coupling, strain and the magnetic field (the Zeeman terms).
- **Symmetry axis** -- each centre points along one of diamond's
  `<111>` crystal directions. Its angles in your lab frame are
  `theta_rad` and `phi_rad`.
- **Transition strength** -- how strongly light drives one
  ground-to-excited transition. There are 4 x 4 = 16 transitions.
- **Cyclicity** (`Lambda`) -- the strength of the transition that keeps
  the spin unchanged divided by the strength of the one that flips it:
  roughly, how many photons the spin emits before it flips. High
  cyclicity makes readout easy. It depends strongly on the angle
  between the field and the symmetry axis.
- **Orbital quenching** -- the crystal weakens how strongly a magnetic
  field acts on the electron's orbital motion. The published reduction
  factors enter the Zeeman factors `f_g`, `f_e`, `delta_g`, `delta_e`.
- **Emission budget** -- where the light from a decay goes. The
  **lifetime** is how long the excited state lasts; the **quantum
  efficiency** is the share of decays that give light at all; the
  **zero-phonon line** (ZPL) is the sharp colour emitted without
  shaking the crystal, and the **Debye-Waller factor** is the share of
  the light in it; the **branching ratio** is the share of that light on
  the one transition the cavity is tuned to. Their product is the
  `radiative_fraction`.
- **Purcell factor** (`F`) -- how much a small optical cavity speeds up
  the emission into its mode. It grows with the cavity's **quality
  factor** `Q` (how long light stays in it) and shrinks with its
  **mode volume** `V`.
- **Cavity linewidth** (`kappa`) and **coupling rate** (`g`) -- `kappa`
  is how fast light leaks out of the cavity, given as a frequency width
  `f_cav / Q`; `g` is how fast the centre and the cavity trade energy.
  In the **bad-cavity, weak-coupling** regime this package covers, both
  `g` and the centre's own decay are much slower than `kappa`, so the
  cavity simply speeds up the emission. The **cooperativity**
  `4 g^2 / (kappa gamma)` measures the same speed-up; here it equals the
  Purcell factor times the radiative fraction.
- **Detection efficiency** (`eta`) -- the fraction of emission events
  that end as a click on the detector, after every loss on the way.
- **Single-shot readout fidelity** -- the probability that one
  measurement names the spin state correctly. You count photons in a
  time **window** `tau`; the spin is called "bright" if the count
  reaches a **threshold**. 0.5 is a coin toss, 1 is perfect.
- **Dark counts / leak / switch-on** -- ways the dark spin state can
  still produce clicks: detector background (`noise_rate`, which adds
  clicks to both spin states), weak off-resonant scattering (`leak`),
  and an off-resonant excitation that flips the dark spin into the
  bright state during the window (switch-on).
- **Heralded entanglement** (Barrett-Kok scheme) -- two nodes each
  send a photon; a joint detection announces that the two spins are
  now entangled.

## Install, units and conventions

```
pip install vacspin
```

It needs Python 3.9 or newer and NumPy 1.22 or newer, and nothing else.
For development: clone the repository and `pip install -e .[test]`,
then run `pytest`.

Units and conventions (each function's docstring states its own):

- **Energies and spin frequencies in GHz** (ordinary frequency, that is
  angular frequency divided by 2 pi): spin-orbit and strain constants,
  energies, qubit frequency, splittings, Rabi rates and fitted
  observables.
- **Magnetic fields in tesla**, as `(Bx, By, Bz)` vectors in your
  **lab frame**. The package rotates them into the centre's own frame
  using `theta_rad` and `phi_rad`. `field_on_circle(angle_deg, B)`
  gives a field of size `B` at an angle in degrees in the lab x-z
  plane (or the x-y or y-z plane).
- **Cavity quantities in Hz**: `omega_q_hz` (the qubit frequency in Hz,
  despite the name not an angular frequency), `kappa_hz` (the cavity
  linewidth, `f_cav / Q_L`), `g_hz` (the coupling rate, likewise given
  as an ordinary frequency: `g / (2 pi)` in Hz).
- **Decay and count rates in 1/s**: the optical decay rate `gamma`
  (`1/lifetime`), `gamma_cav`, `noise_rate`. Times in seconds.
- **Mode volume** `v_rel` in units of `(wavelength / refractive
  index)^3`.
- The **gyromagnetic ratio** is `GAMMA_GHZ_PER_T = 28.0` GHz/T.

## Examples

Each example below runs as written, and the output shown is what it
printed with vacspin 0.4.0. The cavity, detector, attempt-rate and
noise values are illustrative, not a recommended design; the centre's
parameters come from the cited sets described in
[Cited parameter sets](#cited-parameter-sets).

### 1. The spin model: splittings, qubit frequency, cyclicity

```python
import vacspin as vs

spin = vs.snv_rosenthal2023()          # cited SnV- parameter set
print(spin.reference)

# Zero-field splittings of the ground and excited manifolds (GHz)
print(f"ground splitting:  {vs.zero_field_splitting(spin.lam_g, spin.ups_g):.2f} GHz")
print(f"excited splitting: {vs.zero_field_splitting(spin.lam_e, spin.ups_e):.2f} GHz")

# Qubit frequency at 125 mT, 147 degrees in the lab x-z plane
b = vs.field_on_circle(147.0, 0.125)   # lab-frame field vector in tesla
print(f"qubit frequency:   {vs.qubit_frequency(spin, b):.3f} GHz")

# Cyclicity (photons per spin flip, roughly) at two field angles, 180 mT
for angle in (53.0, 147.0):
    lam = vs.cyclicity(spin, vs.field_on_circle(angle, 0.180))
    print(f"cyclicity at {angle:.0f} deg: {lam:.1f}")
```

```
Rosenthal et al., PRX 13, 031022 (2023), Table I (device values); quenching: Thiering & Gali, PRX 8, 021063 (2018)
ground splitting:  903.00 GHz
excited splitting: 2999.99 GHz
qubit frequency:   3.640 GHz
cyclicity at 53 deg: 7.9
cyclicity at 147 deg: 158.5
```

The measured values for this device are a 902.98 GHz ground splitting,
a 3.677 GHz qubit frequency at this field, and a cyclicity of 8.6 at
53 degrees. The model gets these with no fitted numbers; the tests
hold it to 0.2 %, 2 % and 10 % of them respectively. The measured
cyclicity near alignment is 2244. At the nominal 147-degree angle the
model gives only 158.5; the test suite reproduces 2244 (within 10 %)
for a field tilted 9 degrees from the symmetry axis, so the
difference is consistent with the field in the experiment being about
10 degrees off the axis rather than exactly at the nominal angle. That is why example 2 uses the measured
2244 as its input.

`transition_table(spin, b)` gives the energies and strengths of all 16
transitions, and `rabi_rate(spin, b, b_drive)` the microwave Rabi rate
in GHz.

### 2. The cavity budget and a fast readout

```python
import vacspin as vs

spin = vs.snv_rosenthal2023()
f_q = vs.qubit_frequency(spin, vs.field_on_circle(147.0, 0.125))   # GHz

ci = vs.CavityInterface(
    q_loaded=500, v_rel=0.68,          # cavity: loaded Q, mode volume in (lam/n)^3
    eta_wg=0.9996,                     # share of cavity light leaving by the useful port
    budget=vs.snv_emission(),          # cited SnV- emission budget
    lambda0=2244.0,                    # bare cyclicity (the measured value)
    omega_q_hz=f_q * 1e9,              # qubit frequency in Hz
    xi_pol=vs.xi_pol_overlap(), xi_pos=0.5,
    t_link=0.99, eta_chip=0.8, eta_det=0.9).require_valid()

print(f"Purcell factor F_C:      {ci.f_c:.2f}")
print(f"lifetime shortened by:   {ci.zeta:.2f}x  ({ci.tau_cav_s * 1e9:.3f} ns)")
print(f"photons into the mode:   {ci.beta:.3f}")
print(f"detected per decay, eta: {ci.eta:.3f}")
print(f"cavity linewidth kappa:  {ci.kappa_hz / 1e9:.1f} GHz")
print(f"coupling rate g:         {ci.g_hz / 1e9:.2f} GHz")
print("regime checks:", ci.validity)

out = vs.fidelity(ci.eta, ci.lambda_cav, ci.gamma_cav, tau=50e-9, s=2.0)
print(f"fidelity in 50 ns: {out['fidelity']:.5f} "
      f"(threshold {out['threshold']} photon, mean bright count {out['n_bright']:.2f})")
```

```
Purcell factor F_C:      19.29
lifetime shortened by:   7.26x  (0.620 ns)
photons into the mode:   0.909
detected per decay, eta: 0.648
cavity linewidth kappa:  968.6 GHz
coupling rate g:         7.52 GHz
regime checks: {'spin_selectivity': True, 'bad_cavity': True, 'weak_coupling': True}
fidelity in 50 ns: 0.99966 (threshold 1 photon, mean bright count 17.31)
```

`xi_pol` and `xi_pos` say how well the centre's dipole lines up with
the cavity field in direction and position (each from 0 to 1);
`xi_pol_overlap()` gives 2/3 for a centre in a `<100>`-oriented
diamond film. `eta` multiplies every loss from the emitter to the
detector: the share of emission into the cavity mode (`beta`), the
out-coupling port, the link, the chip and the detector.

This fidelity counts only the bright spin's own statistics: `leak` and
`noise_rate` are left at zero here, so the dark state never produces a
click. The test suite runs nearly the same design point (link 0.991,
out-coupling `1 - 500/1.39e6`) with a computed
off-resonant leak and a 1000 counts/s background and asserts a
fidelity above 98 % at the best of four windows (10 to 50 ns).

`require_valid()` is not decoration: the Purcell formulas assume that
the cavity loses light much faster than the emitter decays (bad
cavity), that the coupling is weak, and that the cavity picks out one
spin state (spin selectivity). If any of these fails the interface
refuses and names the failed condition:

```python
import vacspin as vs

ci = vs.CavityInterface(q_loaded=2e7, v_rel=0.68, eta_wg=0.999,
                        budget=vs.snv_emission(), lambda0=2244.0,
                        omega_q_hz=3.677e9, xi_pol=vs.xi_pol_overlap(),
                        xi_pos=0.5)
print(ci.validity)
try:
    ci.require_valid()
except ValueError as err:
    print("refused:", err)
```

```
{'spin_selectivity': False, 'bad_cavity': False, 'weak_coupling': False}
refused: interface outside its regime of validity: spin_selectivity, bad_cavity, weak_coupling failed. The Purcell/rate formulas assume the bad-cavity, weak-coupling, spin-selective regime; change Q_L, V or the operating point.
```

### 3. Readout without a cavity, and the inverse questions

```python
import vacspin as vs

gamma = vs.snv_emission().gamma0       # optical decay rate, 1/s (1 / 4.5 ns)

# Slow confocal readout: eta = 0.2 %, cyclicity 2244, 50 us window,
# background 0.2 counts per window.
n_b, n_d = vs.readout_counts(0.002, 2244.0, gamma, 50e-6, s=10.0,
                             noise_rate=0.2 / 50e-6)
print(f"mean counts: bright {n_b:.2f}, dark {n_d:.2f}")
print(f"threshold-1 fidelity: {vs.fidelity_threshold1(n_b, n_d):.4f}")

# How long a window reaches 85 %?
tau = vs.required_window(0.85, eta=0.002, lam=2244.0, gamma=gamma, s=10.0,
                         noise_rate=0.2 / 50e-6, tau_max=1.0)
print(f"window for 85 %: {tau * 1e6:.1f} us")

# What detection efficiency reaches 99 % in a 10 us window?
eta = vs.required_efficiency(0.99, lam=2244.0, gamma=gamma, tau=10e-6, s=10.0)
print(f"efficiency for 99 % in 10 us: {eta:.4f}")

# A target the model cannot reach is refused, with the best value named.
try:
    vs.required_efficiency(0.999, lam=3.0, gamma=gamma, tau=1e-6)
except ValueError as err:
    print("refused:", err)

# The same readout, but the dark spin also scatters light off-resonantly
# (leak = 1 % of the resonant rate), which also switches it on.
out = vs.fidelity(0.002, 2244.0, gamma, 50e-6, s=10.0, leak=0.01,
                  noise_rate=0.2 / 50e-6)
print(f"with leak: fidelity {out['fidelity']:.3f}, mean counts "
      f"bright {out['n_bright_total']:.2f}, dark {out['n_dark_total']:.2f}")
```

```
mean counts: bright 4.22, dark 0.20
threshold-1 fidelity: 0.9020
window for 85 %: 8.8 us
efficiency for 99 % in 10 us: 0.0218
refused: target fidelity 0.999 is unreachable at this operating point: even eta = 1 gives 0.9000. Increase the cyclicity, the window, or reduce the dark counts.
with leak: fidelity 0.500, mean counts bright 4.22, dark 4.30
```

These inputs are the confocal case the test suite uses for the
measured readout of Rosenthal et al., arXiv:2403.13110, which reports
about 4 bright counts and at most 0.2 dark counts in a 50 us window at
cyclicity 2244; the test asserts a bright-minus-dark mean between 3
and 5. The efficiency 0.2 % is chosen so that this package's
mean-count formula gives those 4 counts; the paper's own fitted
efficiency is about 0.1 %, so treat 0.2 % as a matched input, not a
quoted one. `s` is the **saturation parameter** of the optical drive
(how hard it is driven); `fidelity_threshold1` is the closed form for
a threshold of one photon. `required_window` and `required_efficiency`
search the full model (`fidelity`) and refuse when no value reaches the
target. Compared with example 2, the cavity turns a window of
microseconds into one of tens of nanoseconds; the test suite checks
that the required window shrinks by more than 100 times.

In the last case the leak switches the dark spin on within about a
microsecond, so for most of the window both spin states look bright
and one measurement cannot tell them apart (fidelity 0.5). `fidelity`
evaluates this switch-on channel exactly (since 0.4.0; 0.3.1 gave 0.72
here, see [Corrections](#corrections-in-earlier-versions)).

### 4. Remote entanglement

```python
import vacspin as vs

p = vs.barrett_kok_success(0.65)            # both nodes: eta = 0.65
print(f"success per attempt: {p:.4f}")
print(f"ceiling with perfect nodes: {vs.barrett_kok_success(1.0)}")
rate = vs.entanglement_rate(1e6, 0.65)      # 1 MHz attempt rate (your number)
print(f"entangled pairs per second: {rate:.0f}")
```

```
success per attempt: 0.2113
ceiling with perfect nodes: 0.5
entangled pairs per second: 211250
```

The success probability per attempt is `(1/2) eta_A eta_B`. The factor
1/2 is a ceiling of the two-photon scheme itself, not a loss any
detector upgrade removes. The attempt rate depends on your reset and
communication times, so you supply it; the package does not guess it.

### 5. Your sample, your numbers: plan, measure, fit

Every sample is strained differently, so the numbers that describe
your centre should come from your own measurements. The `lab` tools
work in the order a lab does: plan, measure, then calibrate.

```python
import dataclasses
import numpy as np
import vacspin as vs

p0 = vs.snv_rosenthal2023()                      # starting point

# Plan: 6 field points, alternating the two observable kinds
kinds = ["qubit", "orbital_g"] * 3
fields = [vs.field_on_circle(z, 0.15) for z in (20, 60, 100, 140, 160, 80)]

# 1. Before measuring: can these points pin down lam_g and ups_g,
#    and how large would the error bars be at 0.05 GHz per point?
info = vs.parameter_information(kinds, fields, p0,
                                vary=("lam_g", "ups_g"), sigmas_ghz=0.05)
print("identifiable:", info["identifiable"])
print({k: round(v, 2) for k, v in info["sigma"].items()})

# Zero-field splittings alone cannot separate the two:
zf = vs.parameter_information(["orbital_g"] * 4, np.zeros((4, 3)), p0,
                              vary=("lam_g", "ups_g"))
print("zero-field only identifiable:", zf["identifiable"])

# 2. Pick the 4 most informative of the 6 planned points
pick = vs.design_fields(kinds, fields, 4, p0, vary=("lam_g", "ups_g"))
print("chosen points:", pick["indices"])

# 3. "Measure": illustrative synthetic data from a sample whose
#    ups_g is 10 % higher than the cited value, plus 0.05 GHz noise.
truth = dataclasses.replace(p0, ups_g=1.1 * p0.ups_g)
rng = np.random.default_rng(1)
measured = [vs.qubit_frequency(truth, b) if k == "qubit"
            else vs.solve(truth, b)[0][2] - vs.solve(truth, b)[0][0]
            for k, b in zip(kinds, fields)]
measured = np.array(measured) + 0.05 * rng.standard_normal(len(kinds))

fit = vs.fit_spin_parameters(kinds, fields, measured, p0,
                             vary=("lam_g", "ups_g"), sigmas_ghz=0.05)
for name in ("lam_g", "ups_g"):
    print(f"{name}: {fit.values[name]:.2f} +/- {fit.sigma[name]:.2f} GHz "
          f"(true {getattr(truth, name):.2f})")
print(f"chi2 = {fit.chi2:.2f} for {fit.chi2_dof} degrees of freedom")
my_sample = fit.params                           # usable everywhere
print(my_sample.reference[:60] + "...")
```

```
identifiable: True
{'lam_g': 8.24, 'ups_g': 9.63}
zero-field only identifiable: False
chosen points: [1, 0, 5, 3]
lam_g: 824.35 +/- 8.88 GHz (true 830.15)
ups_g: 201.49 +/- 9.09 GHz (true 195.44)
chi2 = 3.19 for 4 degrees of freedom
fitted by vacspin.lab.fit_spin_parameters to 6 measured poin...
```

The observable kinds are `"qubit"` (the qubit frequency), `"orbital_g"`
(the gap between the two ground pairs) and `"orbital_e"` (the same gap
in the excited manifold), all in GHz. Zero-field splittings alone can
never separate spin-orbit coupling from strain, because they depend on
both only through `sqrt(lam^2 + 4 ups^2)`; the planner says so, and
`fit_spin_parameters` and `design_fields` refuse such a design with an
explanation. The error bars come from the standard weighted
least-squares covariance. A `chi2` close to the number of degrees of
freedom (measurements minus fitted parameters) means the fit matches
the data within the stated measurement errors. `design_fields` adds points one at a time,
each time taking the one that most increases the determinant of the
information matrix (greedy D-optimal design); this is a good
heuristic, not a proof of the best possible subset. Observation
records can be saved and reloaded with `save_observations_csv` and
`load_observations_csv`; the round trip returns identical values.

### 6. Error bars on anything you predict from the fit

Continuing example 5 (run it first; `fit` and `truth` come from
there):

```python
# The qubit frequency at a field you did not measure, with its error bar
b_new = vs.field_on_circle(35.0, 0.30)
q = vs.propagate_uncertainty(fit, lambda p: vs.qubit_frequency(p, b_new))
print(f"qubit frequency at 300 mT, 35 deg: {q['value']:.3f} +/- {q['sigma']:.3f} GHz"
      f" (true {vs.qubit_frequency(truth, b_new):.3f})")

# Any function of the parameters works, including several outputs at once
z = vs.propagate_uncertainty(
    fit, lambda p: [vs.zero_field_splitting(p.lam_g, p.ups_g),
                    vs.cyclicity(p, b_new)])
print("ground splitting: {:.2f} +/- {:.2f} GHz".format(z["value"][0], z["sigma"][0]))
print("cyclicity:        {:.1f} +/- {:.1f}".format(z["value"][1], z["sigma"][1]))
```

```
qubit frequency at 300 mT, 35 deg: 4.126 +/- 0.141 GHz (true 4.032)
ground splitting: 917.57 +/- 0.03 GHz
cyclicity:        5.8 +/- 0.4
```

`propagate_uncertainty` takes the slopes of your function with
respect to the fitted parameters and combines them with the fit's
covariance matrix (first-order error propagation). It is exact for a
quantity that depends linearly on the fitted parameters, and a good
approximation when the error bars are small compared with how quickly
the quantity curves; the tests compare it with the spread obtained by
drawing many parameter sets from the fit's uncertainty (within 5 %).
The ground splitting is known far better than `lam_g` and `ups_g`
separately, because the fit pins down the combination
`sqrt(lam^2 + 4 ups^2)` directly; `cov` in the result carries such
correlations.

## Cited parameter sets

No physical number is made up, and none is accepted without a source:
the `reference` field of both parameter classes (`SpinParameters`,
`EmissionBudget`) is mandatory, and the package refuses an empty or
too-short one.

- `snv_rosenthal2023()` + `snv_emission()`: the SnV- device values of
  Rosenthal et al., PRX 13, 031022 (2023) (Table I) with the
  Thiering-Gali quenching factors (Thiering & Gali, PRX 8, 021063
  (2018)), and the emission budget of arXiv:2403.13110 (lifetime
  4.5 ns, quantum efficiency 0.8) / Goerlitz et al., New J. Phys. 22,
  013048 (2020) (Debye-Waller factor 0.57) / Lee et al.,
  arXiv:2511.05740 (branching ratio 0.75, zero-phonon line 619 nm). Checked against measurement in the tests
  with no fitted numbers: the 902.98 GHz splitting, the 3.677 GHz qubit
  frequency, the cyclicity landscape (8.6 at 53 degrees, 2244 near
  alignment), and a Rabi rate on the MHz scale (tolerances in
  [How the results are checked](#how-the-results-are-checked)).
- `siv_hepp2014()`: the unstrained SiV- splittings (50 and 260 GHz,
  zero-phonon line 737 nm) of Hepp et al., PRL 112, 036405 (2014).
  Strain and orbital quenching are sample-specific and deliberately
  not shipped (set to zero); the reference string says so.
- `gev_bhaskar2017()`: the measured GeV- splittings (152 and
  981 GHz, zero-phonon line 602 nm) of Bhaskar et al., PRL 118,
  223603 (2017), with the same deliberate omissions as SiV- (Siyushev
  et al., PRB 96, 081201(R) (2017) measured 170 GHz on a strained
  emitter): strain and quenching are your sample's numbers, and
  `vacspin.lab` exists to fit them.

For SiV- and GeV- the default orientation is the `<111>` axis of a
`(001)`-oriented sample; pass `theta_rad` and `phi_rad` to change it.
No emission budget ships for SiV- or GeV-; build an `EmissionBudget`
from your own cited numbers.

The lead-vacancy (PbV-) centre is documented but not shipped, on
purpose: its ground splitting is measured (about 3900 GHz; Wang et
al., ACS Photonics 8, 2947 (2021), confirmed by Chen et al.,
arXiv:2605.27841 (2026)), but no verified measured excited-state
splitting was available, and this package does not ship half a
parameter set. For PbV-, or your own sample of any centre, fill in
`SpinParameters` and `EmissionBudget` from your measurements or the
literature; the provenance travels with every prediction.

## What is in the package

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

**Parameter sets** (`vacspin.params`)

- `SpinParameters` -- spin-orbit constants `lam_g`, `lam_e`, strain
  `ups_g`, `ups_e`, orbital and spin Zeeman factors `f_g`, `f_e`,
  `delta_g`, `delta_e`, the axis angles `theta_rad`, `phi_rad`, and a
  required `reference`.
- `EmissionBudget` -- lifetime `tau0_s`, quantum efficiency `eta_q`,
  zero-phonon fraction `eta_dw`, branching ratio `eta_br`, zero-phonon
  wavelength `zpl_nm`, and a required `reference`; `gamma0` and
  `radiative_fraction` are derived from them.
- `snv_rosenthal2023`, `snv_emission`, `siv_hepp2014`,
  `gev_bhaskar2017` -- the cited sets above.

**Spin Hamiltonian** (`vacspin.hamiltonian`; Rosenthal et al., PRX 13,
031022 (2023), Appendix B)

- `h_manifold` -- the 4 x 4 Hamiltonian of one manifold (spin-orbit +
  strain + spin and orbital Zeeman), with the field in the centre's
  frame.
- `lab_to_spin` -- rotates a lab-frame vector into the centre's frame.
- `solve` -- energies and eigenvectors of both manifolds at a
  lab-frame field; `strain_scale` multiplies the strain (a
  strain-tuning knob).
- `qubit_frequency` -- the splitting of the lowest ground pair.
- `zero_field_splitting` -- the closed form `sqrt(lam^2 + 4 (Ux^2 +
  Uy^2))`.
- `qubit_frequency_perpendicular` -- the approximate formula
  `2 gamma B_perp ups / Delta` for a purely sideways field (Eq. B9),
  valid for small fields; use `qubit_frequency` for anything
  quantitative.
- `GAMMA_GHZ_PER_T` -- the gyromagnetic ratio, 28.0 GHz/T.

**Transitions** (`vacspin.transitions`)

- `transition_strength`, `transition_table` -- one or all 16
  transition strengths, and their energies.
- `cyclicity` -- the cyclicity of the lowest ground state on the lower
  (`manifold="lower"`) or upper excited pair; infinite when the
  spin-flip line is exactly dark.
- `rabi_rate` -- the microwave Rabi rate (GHz) for a drive field.
- `field_on_circle` -- a field of given size at an angle in one lab
  plane.

**Cavity** (`vacspin.cavity`; conventions of Lee et al.,
arXiv:2511.05740)

- `purcell_max` -- the ideal Purcell factor `(3 / 4 pi^2) Q / V`.
- `lorentzian_suppression` -- how much a transition detuned from the
  cavity is enhanced less: `1 / (1 + (2 Q delta / f_cav)^2)`.
- `xi_pol_overlap` -- the polarisation overlap of a `<111>` dipole
  with an in-plane cavity field.
- `CavityInterface` -- the full budget: `f_max`, `f_c` (Purcell factor
  of the coupled line), `f_c2` (of the detuned spin-flipping partner),
  `zeta` (lifetime shortening), `gamma_cav`, `tau_cav_s`, `beta`,
  `eta`, `lambda_cav` (cavity-boosted cyclicity), `kappa_hz`, `g_hz`,
  `cooperativity`, the `validity` checks, `require_valid()` and
  `summary()`. An optional `delta_partner_hz` adds the Purcell factor
  of a competing same-spin line and a check that the cavity is
  narrower than their spacing.

**Readout** (`vacspin.readout`; Rosenthal et al., arXiv:2403.13110,
Appendix A)

- `polarization_rate` -- how fast the bright record ends through spin
  flips.
- `readout_counts` -- mean bright and dark counts in a window.
- `fidelity_threshold1` -- the closed-form fidelity at a threshold of
  one photon (Eq. A20).
- `fidelity` -- the fidelity at the best whole-number threshold, from
  the full count distributions of both spin states: the bright one
  (spin flips end its record) plus background, and the dark one
  (background, `leak`, and the switch-on channel, evaluated exactly).
  Besides `fidelity` and `threshold` it returns `n_bright` (mean signal
  count), `n_dark` (mean background plus leak count),
  `n_bright_total` and `n_dark_total` (the means of the two full
  distributions) and `truncation_bound` (a guaranteed bound, at most
  1e-9, on the error from tracking a finite number of counts).
- `required_efficiency`, `required_window` -- the smallest efficiency
  or shortest window that reaches a target fidelity.
- `poisson_pmf`, `geometric_pmf` -- the two limiting count
  distributions (no spin flips; flips always end the record).

**Remote entanglement** (`vacspin.remote`; Barrett and Kok, PRA 71,
060310(R) (2005))

- `barrett_kok_success` -- success probability per attempt,
  `(1/2) eta_A eta_B`.
- `entanglement_rate` -- attempt rate times that probability.

**Your own sample** (`vacspin.lab`)

- `parameter_information` -- before measuring: is the plan
  identifiable, and what error bars would it give?
- `design_fields` -- the most informative subset of candidate points.
- `fit_spin_parameters` -- fit any of the ten spin parameters to your
  measured frequencies; returns a `SpinFit` with `params` (a ready
  `SpinParameters` whose `reference` records the fit and the base set),
  `values`, `sigma`, `cov`, `chi2`, `chi2_dof` and more.
- `propagate_uncertainty` -- after fitting: the value and error bar
  (and covariance) of any quantity computed from the fitted
  parameters, by first-order error propagation.
- `save_observations_csv`, `load_observations_csv` -- a plain CSV
  format for observation records (columns `kind, bx_t, by_t, bz_t,
  value_ghz` and optionally `sigma_ghz`).

`__version__` gives the installed version.

## When it refuses, and why

`vacspin` raises an error instead of guessing when:

- a parameter set has no real `reference` (fewer than 8 characters), a
  number that is not finite (NaN or infinity), a spin-orbit constant
  that is not positive, a negative strain, a lifetime or wavelength
  that is not positive and finite, or an efficiency outside (0, 1];
- the magnetic field is not finite, `strain_scale` is negative, or the
  sideways field given to `qubit_frequency_perpendicular` is negative
  or not finite;
- a manifold or field-plane name is unknown;
- a cavity's `Q` or `V` is not positive and finite, an efficiency or
  overlap is outside [0, 1], the bare cyclicity is not finite and
  positive, or the qubit frequency or `delta_partner_hz` is not
  positive and finite;
- `require_valid()` finds the interface outside the bad-cavity,
  weak-coupling or spin-selective regime (and, when
  `delta_partner_hz` is given, the partner-selectivity condition); it
  names every failed condition;
- a readout input is out of range: `eta` outside [0, 1], a cyclicity
  that is not finite and positive, a decay rate, window or saturation
  parameter that is not positive and finite, a `leak` or `noise_rate`
  that is negative or not finite, `f0` outside (0, 1], or mean counts
  that are negative or not finite;
- `fidelity` would need count distributions wider than 2^18 = 262144
  counts to stay within its 1e-9 accuracy bound (a very long window
  with a very high background);
- a target fidelity is outside (0.5, 1), or cannot be reached (by any
  efficiency up to 1, or any window up to `tau_max`, which must be
  finite and longer than 1 ns); the message gives the best value found;
- a detection efficiency is outside [0, 1] or the attempt rate is not
  positive and finite (remote entanglement);
- in `vacspin.lab`: an unknown observable kind or parameter name, a
  repeated or empty `vary`, fields of the wrong shape or not finite,
  measurement errors that are not positive, too few observations for
  the parameters (and, without `sigmas_ghz`, for an error scale), a
  design that cannot tell the parameters apart, an `n_pick` outside
  its range, a fit that does not converge (`RuntimeError`), a CSV
  file with a wrong header, a short row, a non-numeric value or no
  data, or (in `propagate_uncertainty`) a function that returns a
  non-finite value or slope.

## How the results are checked

84 automated tests run on every push and pull request, on Python 3.9,
3.10, 3.11, 3.12, 3.13 and 3.14, and once more on Python 3.10 with the
oldest versions the package allows (NumPy 1.22.0, pytest 7.0.0). The
numerical checks compare the package with an exact formula, a
published measurement, or a second calculation done a different way.
The rest check that the refusals fire, plus one check that the version
numbers agree and one that every exported name exists. The main
checks, with the tolerances the tests actually use:

**Spin Hamiltonian and transitions**

- The Hamiltonian is Hermitian (equal to its own conjugate transpose,
  as any energy matrix must be) for 20 random inputs (to 1e-12).
- At zero field each manifold is two degenerate pairs (to 1e-9 GHz),
  and their gap equals the closed-form zero-field splitting (to 1 part
  in 10^9), for 10 random parameter sets with strain.
- The frame rotation keeps vector lengths (to 1e-12) and is
  orthonormal (a pure rotation, to 1e-12).
- The approximate transverse-field qubit frequency agrees with full
  diagonalisation within 2 % at 50 mT.
- Without strain, a field along the axis splits the qubit by exactly
  `gamma B` (to 1e-9 GHz), and a sideways field splits it by less than
  1 % of that.
- Scaling the strain moves the ground splitting along the closed form
  (to 1 part in 10^9).
- The 16 transition strengths from each ground state add up to 3 (to
  1e-10) for SnV- and SiV- at three fields, and for GeV- at two.
- For a field exactly along the symmetry axis the cyclicity is above
  10^12 for SnV-, SiV- and GeV- (for SnV- the spin-flip strength is
  also checked to be below 1e-20).
- The Rabi rate is zero without a drive, exactly doubles when the
  drive doubles (to 1 part in 10^12), and lies between 1 and 20 MHz for
  a 0.6 mT drive (measured: about 6.25 MHz).
- The Hamiltonian built from precomputed matrices (since 0.4.0, for
  speed) equals the term-by-term construction for 50 random inputs
  (to 1e-12 GHz).

**Against published measurements (SnV-, no fitted numbers)**

- Ground splitting within 0.2 % of 902.98 GHz; excited within 0.2 % of
  3000 GHz.
- Qubit frequency within 2 % of the measured 3.677 GHz at 125 mT,
  147 degrees.
- Cyclicity within 10 % of 8.6 at 53 degrees, 180 mT, and within 10 %
  of 2244 for a 180 mT field tilted 9 degrees from the axis.
- The confocal readout of arXiv:2403.13110: a bright-minus-dark mean
  count between 3 and 5 (measured about 4), the dark count equal to
  the stated 0.2 (to 1e-9), and a threshold-1 fidelity between 0.82
  and 0.92 with a preparation factor `f0 = 0.935` (measured 0.874).
- The GeV- and SiV- values equal their sources exactly, and so do
  their closed-form zero-field splittings; for GeV- the full solver
  also gives 152 and 981 GHz (to 1e-9 GHz).

**Cavity**

- `purcell_max` gives exactly 1 for `Q/V = 4 pi^2 / 3` (to 1e-12),
  and the coupling
  rate satisfies `F_C = 4 g^2 / (kappa gamma_rad)` with every rate in Hz
  (to 1 part in 10^9); `g_hz` also equals the angular-unit value
  divided by 2 pi and the cooperativity equals `F_C` times the
  radiative fraction (to 1 part in 10^12).
- The Lorentzian (the cavity's response versus detuning) is 1 on
  resonance, 1/2 at `delta = kappa/2` (to
  1e-12), and below 1e-9 far away; the polarisation overlap is 2/3 (to
  1e-6) at the default angles and 0 for a dipole along the normal.
- `beta` and `eta` stay in (0, 1] for Q from 10 to 10^7, `beta` rises
  with Q and passes 0.99, and `zeta` follows its formula (to 1e-12).
- The cavity-boosted cyclicity is at least the bare one, follows its
  formula (to 1 part in 10^9), and grows by less than 5 % for a very
  broad cavity (Q = 10).
- `require_valid()` passes at Q = 500 and refuses at Q = 2 x 10^7.

**Readout statistics**

- The Poisson and geometric distributions add up to 1 (to 1e-12 and
  1e-9) and have the right means (to 1e-9 and 1e-6 of `1 + mean`).
- With almost no spin flips, the bright-count distribution equals the
  Poisson distribution (to 1e-10 per value).
- When flips always end the record, it equals the geometric
  distribution with mean `eta (Lambda + 1)` (to 1e-10); the
  per-emission version with mean `eta Lambda` differs by more than
  1e-4 but less than `10 / Lambda^2`.
- Where the best threshold is 1, `fidelity` equals `1 - P(no bright
  photon)/2` (to 1e-9) and agrees with the closed form of Eq. A20
  within 0.02.
- The fidelity stays in [0.5, 1] and rises with `eta`.
- `polarization_rate` reaches `gamma/2` at strong drive (within
  `1e-3 gamma`), falls below 1e-4 of its value at a large detuning, and
  scales as `1/(1 + Lambda)` (to `1e-12 gamma`).
- `required_efficiency` returns an efficiency that reaches the target
  (and half of it does not); `required_window` returns a window that
  reaches it while a window 0.1 % shorter does not; both refuse
  unreachable targets.
- `required_window` below its 1 ns search grid: for an emitter without
  spin flips or background the shortest window has the closed form
  `ln(1 / (2 (1 - F))) / (eta R)`; the result agrees to 1 part in 10^5.
- Background in the bright state: without spin flips and leak both
  states are Poisson, and `fidelity` (value and threshold) equals the
  best-threshold Poisson closed form (to 1e-9), in two cases including
  150 signal counts on 200 background counts.
- The mean bright count for very rare spin flips (cyclicity 10^9 to
  10^15) matches its series limit `eta R tau (1 - Gp tau / 2)` (to 1
  part in 10^12).
- The dark switch-on distribution matches a direct double numerical
  integration over the switch and flip times (to 1e-6 per value, two
  cases covering both evaluation branches), sums to 1 (to 1e-12), and
  its mean equals the numerically integrated survival function times
  the count rate (to 1 part in 10^9), in four cases.
- Monte Carlo: 400 000 simulated shots per spin state of the counting
  process; the fidelity at the returned threshold and the mean bright
  and dark counts agree with `fidelity` within 5 standard errors, in
  three cases
  (slow switch-on, switch-on within the first microsecond, and the
  cavity design point).
- `readout_counts` and `fidelity` differ, as documented, only by the
  leak counts `eta leak R tau` in the bright mean (to 1e-12 relative;
  zero without leak).

**The source study's design point** (Q = 500 cavity, from the study's
stated inputs, with the qubit frequency computed by the package)

- Purcell factor within 5 % of 19 and detection efficiency within 5 %
  of 0.65, inside the valid regime.
- With a computed off-resonant leak and a 1000 counts/s background:
  fidelity above 98 % at the best of four windows (10, 20, 30,
  50 ns); above 94 % at 100 ns with `f0 = 0.99`; and lower when `eta`
  is halved.
- The window needed for 85 % fidelity is more than 100 times shorter
  with the cavity than confocally.

**Your own sample** (`vacspin.lab`)

- The observables equal `qubit_frequency` and the closed-form
  splittings (to 1e-9 GHz).
- A fit to noise-free synthetic data, started 8 % and 15 % off,
  recovers `lam_g` and `ups_g` to 1 part in 10^5.
- In a linear case the error bar equals `sigma / sqrt(n)` (to 1e-6
  GHz), from both the fit and the planner, and the two agree to 0.1 %
  in a nonlinear case.
- Over 400 seeded synthetic experiments, the scatter of the fitted
  values matches the reported error bars within 15 %.
- For zero-field data only, the information matrix has rank one (it
  carries information about only one combination of the two
  parameters: its second singular value is below 1e-10 of the first);
  the planner reports the design as not identifiable, and the fit and
  the design tool refuse it.
- The greedy design's chosen set is at least as informative as each
  of 30 random subsets of the same size. (The rank-one determinant
  identity behind the greedy rule is checked on random matrices.)
- A CSV round trip returns identical values.
- The faster eigenvalue-only path used by the fits equals the full
  solver for all three observables, 20 random fields and strain
  scales (to 1e-9 GHz).
- `propagate_uncertainty` is exact where it should be: a fitted
  parameter returns its own error bar (to 1e-9 relative), the
  zero-field splitting at zero strain returns `sigma / sqrt(n)` (to
  1e-9 GHz), and a fixed linear combination `A x` returns
  `A C A^T` (to 1e-8 relative). For a nonlinear prediction (the qubit
  frequency at a new field) it agrees within 5 % with the spread of
  4000 parameter sets drawn from the fit's uncertainty. At a
  parameter boundary (zero strain) the slope is taken one-sided and
  is still right for a linear quantity (to 1e-6 relative), and an
  error raised by your own function reaches you instead of being
  hidden.

**Inputs**

- NaN and infinite values are refused across the parameter classes,
  the cavity, readout and remote-entanglement functions (with array
  input to `polarization_rate` still accepted).

## Corrections in earlier versions

**0.4.0 fixed four problems in the readout model and the inputs.**

- `fidelity` left the detector background (`noise_rate`) out of the
  bright spin's counts, although detector dark counts and stray light
  do not depend on the spin, and this package's own `readout_counts`
  includes them in the bright mean.
  Background counts can only help the bright state pass the
  threshold, so the old fidelity was too low, sometimes badly: with
  150 signal counts on 200 background counts it returned 0.500
  instead of 0.999998.
- The dark switch-on channel was a 24-point sum over switch times and
  was added to the background error and capped at 1. When the switch
  happens early in the window the sum missed most of it: for the
  confocal readout with `leak = 0.01` (example 3) 0.3.1 gave 0.718,
  while the model's exact value and a Monte Carlo simulation give
  0.500. It is now evaluated exactly and combined with the background
  exactly.
- The mean bright count `eta (Lambda + 1)(1 - exp(-Gp tau))` lost
  digits for very large cyclicity (relative error 6e-5 at 10^15); it
  now uses `expm1`.
- NaN or infinite inputs could pass silently and give NaN or infinite
  results (for example a NaN spin-orbit constant, `purcell_max(nan,
  1)`, `entanglement_rate(inf, 0.5)`), and `CavityInterface` with
  `Q = 0` stopped with a `ZeroDivisionError`. These are now refused
  with a clear message. `required_window` returned 1 ns whenever the
  target was already reached at 1 ns; it now finds the shorter window.

What changed in numbers: see the "Behaviour changes" note in
[CHANGELOG.md](CHANGELOG.md). With no background and no leak every
result is unchanged.

**0.3.1 fixed three problems found in a review of 0.3.0.**

- `CavityInterface.g_hz` mixed units: it multiplied the cavity
  linewidth in Hz by the decay rate in 1/s, which made the reported
  coupling rate (`g_hz`, and `g_mhz` in `summary()`) too large by
  `sqrt(2 pi)`, about 2.5 times, and made the weak-coupling check
  refuse some designs it should have accepted. It now uses the decay
  linewidth in Hz, `gamma / (2 pi)`, as the bad-cavity check already
  did. The Purcell factor, efficiencies, lifetime, cyclicity and
  cooperativity are unchanged. If you used `g_hz` or the
  `weak_coupling` verdict from 0.3.0 or earlier, please re-check them.
- `readout_counts` accepted a negative `leak` or `noise_rate` and
  returned a negative dark count; it now refuses them, as `fidelity`
  already did.
- The readout functions accepted an infinite window, decay rate or
  saturation parameter (for example `fidelity` returned NaN for
  `tau = inf`); they now refuse non-finite values with a clear message.

Also in 0.3.1: the CI matrix now includes Python 3.10, and a CI job
tests the oldest allowed dependencies.

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

## Limits

- **No photonic device design.** The cavity's Q, mode volume and
  out-coupling are inputs you simulate or measure; the companion
  study's pipeline shows how.
- **Weak-coupling, bad-cavity regime only.** The readout model is the
  two-level bad-cavity treatment. `require_valid()` refuses the
  strong-coupling regime rather than describing it wrongly. Its
  thresholds are fixed factors: decay linewidth below a tenth of the
  cavity linewidth, `g` below a tenth of the cavity linewidth, and
  decay linewidth below a fifth of the qubit frequency.
- **The readout model is a two-state counting model.** Both count
  distributions are evaluated exactly within it, but it leaves out
  some second-order effects: once a bright spin has flipped to dark
  it produces no leak counts and cannot switch back on, and a dark
  spin's leak counts are counted for the whole window even after it
  has switched on (the leak is small whenever the readout is good).
  The `nt_switch` argument of `fidelity` is ignored since 0.4.0.
- **Leak counts belong to the dark spin only.** In `fidelity` the
  bright spin gets the background (`noise_rate`) but no `leak`
  counts, because the leak is off-resonant scattering of the dark
  spin. `readout_counts` keeps the simpler bookkeeping of its source
  and adds the whole dark mean, leak included, to the bright mean, so
  with `leak > 0` its bright mean exceeds `fidelity`'s
  `n_bright_total` by exactly `eta leak R tau` (example 3 with
  `leak = 0.01`: 4.3176 against 4.2166). With `leak = 0` they agree.
- **The measured cyclicity near alignment depends on alignment.** At
  the nominal angles the model gives a much lower cyclicity than the
  measured 2244; see example 1.
- **Spin coherence times** (T1, T2) are sample-dependent measurements
  and are not part of the model.
- **SiV- and GeV- ship without strain, quenching or emission budget**;
  supply or fit your own.
- **The `lab` fits** use a local Levenberg-Marquardt search (a standard
  step-by-step least-squares method) from your starting set; the error
  bars assume independent Gaussian measurement errors, and the greedy
  design is a heuristic. `propagate_uncertainty` is first-order: when
  an error bar is large compared with the scale on which the predicted
  quantity curves (or reaches a boundary such as zero strain), sample
  the fit's uncertainty instead.

## Where it comes from

`vacspin` is the general-purpose engine distilled from the
SnV-on-TFLN (tin-vacancy centre on thin-film lithium niobate) interface
study below, and generalised to any group-IV
colour centre and any cavity. The study's repository holds the device
design (FEM/EME/GME electromagnetic simulations) and reproduces the specific study.

> T. M. Mahim, M. M. Rahman and A. S. M. Mohsin, "Fast single-shot
> readout of tin-vacancy spins with an overcoupled diamond nanocavity
> on thin-film lithium niobate" (submitted to Optics Express, 2026);
> pipeline: https://github.com/Tanvir-Mahmud-Mahim/snv-tfln-cavity-interface

Formalism: Rosenthal et al., PRX 13, 031022 (2023) and
arXiv:2403.13110; transition-resolved cavity conventions: Lee et al.,
arXiv:2511.05740; heralding: Barrett and Kok, PRA 71, 060310(R)
(2005); optimal design: F. Pukelsheim, Optimal Design of Experiments,
SIAM (2006).

## Citing, support and license

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

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

Support runs through the
[issue tracker](https://github.com/TaN-MM-Org/vacspin/issues). Usage
questions are welcome alongside bug reports; a docstring that left a
unit or a sign convention unclear is treated as a documentation bug,
not user error. While the version is below 1.0 the API may still move
between minor versions; such changes are called out in the release
notes.

Licensed under Apache-2.0.
