Metadata-Version: 2.4
Name: absnoise
Version: 0.15.0
Summary: Andreev-bound-state occupation noise and detector budgets for proximity Josephson junctions
Author-email: "Tanvir M. Mahim" <tanvir.mahim@bracu.ac.bd>
License-Expression: Apache-2.0
Project-URL: Homepage, https://tanvir-mahmud-mahim.github.io/software/
Project-URL: Repository, https://github.com/TaN-MM-Org/absnoise
Project-URL: Issues, https://github.com/TaN-MM-Org/absnoise/issues
Project-URL: Changelog, https://github.com/TaN-MM-Org/absnoise/releases
Keywords: Andreev bound states,Josephson junction,graphene,thermal detector,bolometer,calorimeter,occupation noise,superconductivity
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
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.24
Requires-Dist: scipy>=1.10
Provides-Extra: plot
Requires-Dist: matplotlib; extra == "plot"
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Dynamic: license-file

# absnoise

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

`absnoise` is a Python package about one question: **how well can a
superconducting junction measure temperature, or detect a single
photon?** The junctions in question are proximity Josephson junctions:
two superconductors joined by a short piece of normal (not
superconducting) conductor, in the shipped examples graphene. Such a
junction carries its supercurrent (a current that flows with zero
resistance) through a few discrete energy levels, the Andreev levels. Each level is
randomly filled and emptied, even in perfect thermal equilibrium. That
flicker -- not the readout electronics -- sets the final limit on how
well the junction can sense temperature.

The package follows the whole chain, from the levels to a detector:

- Where are the Andreev levels, and how do they move with phase and
  temperature?
- How large and how fast is the occupation flicker, and what
  temperature resolution does it allow?
- For a real device read out by a microwave resonator: what are the
  frequency-noise spectrum, the noise-equivalent power and the energy
  resolution for a single photon?

It also works backward from measured data: a critical-current curve, a
readout time trace, or a noise spectrum. And it helps plan those
measurements before the fridge time is spent.

The main results are checked by automated tests against exact formulas
or a second, independent calculation (see
[How the results are checked](#how-the-results-are-checked)). When a
question cannot be answered from the data or parameters given, the
package stops with an error that says why, instead of 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](#install)
- [Units and conventions](#units-and-conventions)
- [Examples](#examples) (each with the output it prints)
- [Tutorials](#tutorials) (three end-to-end scripts in `examples/`)
- [What changed in 0.14.0](#what-changed-in-0140)
- [What changed in 0.13.0](#what-changed-in-0130)
- [What changed in 0.12.0](#what-changed-in-0120)
- [What changed in 0.11.0](#what-changed-in-0110)
- [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

**The junction**

- **Superconducting gap** (`Delta`) -- the energy a superconductor
  needs to break a pair of electrons. It shrinks as the temperature
  rises and closes at the **critical temperature** `Tc`. Here the gap
  is the one the superconductor induces in the junction (the
  "proximity gap"), taken as `Delta = 1.764 kB Tc` at zero temperature.
- **Phase** (`phi`) -- the phase difference across the junction, in
  radians. The supercurrent depends on it (the **current-phase
  relation**); the largest supercurrent is the **critical current**
  `Ic`.
- **Andreev levels** -- discrete energy levels inside the gap that
  carry the supercurrent. Their energies depend on the phase.
- **Transparency** (`tau`, between 0 and 1) -- how easily electrons
  pass the contacts. `tau = 1` is a perfectly clean contact.
- **Channel** -- one of the separate paths (quantum modes) that
  electrons can take across the junction; a wide junction has many.
- **Short junction** -- a junction much shorter than the
  superconducting **coherence length** `xi` (roughly, how far the
  electron pairing of the superconductors reaches into the normal
  conductor). Then each channel has one level, `E = Delta sqrt(1 - tau sin^2(phi/2))`. For a longer junction
  the package solves the exact levels.
- **Recipe** -- one junction design: materials, `Tc`, size, gate
  voltage, transparency, measured critical current and resistance.
- **Matsubara sum**: a standard way to write a thermal quantity as a
  sum over the frequencies `(2n+1) pi kB T / hbar` (n = 0, 1, 2, ...).
  For a junction it gives the free energy of the levels and of the
  continuum together, without finding the levels.

**The noise**

- **Occupation** -- whether a level holds an extra electron-like
  excitation (a quasiparticle). It flips between filled and empty at
  random: a **telegraph** signal.
- **Correlation time** (`tauA`) -- how long an occupation typically
  lasts before it changes. The package does not predict `tauA`; you
  supply it, or measure it with the fitting tools.
- **Power spectral density (PSD)** -- how the noise power is spread over
  frequency. A telegraph signal has a **Lorentzian** spectrum: flat up
  to a **knee** at `1/(2 pi tau)`, then falling. **White noise** has
  the same density at all frequencies. A **periodogram** is the PSD
  estimated from one recorded trace; averaging many periodograms
  reduces their scatter.
- **Monte Carlo** -- a simulation driven by random numbers, used here
  to generate example noise traces and check the formulas against
  them.
- **Heat capacity** -- energy needed to warm something by 1 K. `C_e`
  is that of the junction's electrons; `C_A` is the part that comes
  from the Andreev occupations.
- **Thermal conductance** `G` and **thermal time** `tau_th = C_e / G`
  -- how fast the electrons lose heat to the crystal lattice (the
  **phonons**). Random heat exchange with the phonons causes
  **phonon thermal-fluctuation noise** (TFN), the classic noise limit
  of a thermal detector.
- **Cauchy-Schwarz bound** -- the package's lower limit on the
  temperature error from occupation noise,
  `var(T) >= 2 kB T^2 tauA / (C_A t)` for averaging time `t`. A
  uniform short junction reaches it exactly.

**The detector**

- **Readout** -- the junction ends a microwave resonator. A change in
  temperature changes the junction's **Josephson inductance** `L_J`,
  which moves the resonance frequency. `y` is the fractional frequency
  shift. Two readouts are named by `which`: `"L"` (the inductance at
  phase near 0; the default of the device budget) and `"I"` (the
  supercurrent at the phase where it is largest).
- **Noise-equivalent power (NEP)** -- the input power that would give a
  signal as large as the noise in a 1 Hz bandwidth (half a second of
  averaging, with the convention below), in W/sqrt(Hz). Smaller is
  better.
- **Energy resolution** -- the uncertainty (one standard deviation)
  with which a deposited energy, for example one photon, can be
  measured, using the best ("matched") filter. Smaller is better.
- **Allan variance** -- a measure of how much a reading changes from
  one averaging window to the next. It shows how long averaging still
  helps.
- **Hidden Markov model (HMM)** -- a statistical model of a hidden
  two-state signal seen through noise. Used here to recover the
  occupation from a noisy readout trace.

## Install

```
pip install absnoise
```

It needs Python 3.9 or newer, NumPy 1.24 or newer and SciPy 1.10 or
newer. The optional extra `absnoise[plot]` installs matplotlib; no part
of the package imports it. For development: clone the repository and
run `pip install -e .[test]`, then `pytest`.

## Units and conventions

- SI units throughout: kelvin, joules, seconds, amperes, hertz, metres.
  Phases are in radians.
- Spectra are **single-sided** (positive frequencies only). With this
  convention, the variance of a `t`-second average of white noise with
  spectral density `S` is `S / (2 t)`, and a Lorentzian with variance
  `var` and correlation time `tau` has plateau `S(0) = 4 var tau`.
- The shipped constants (`KB`, `HBAR`, ...) are CODATA 2018 values.
- Graphene thermal parameters follow Lee et al. (see
  [Where it comes from](#where-it-comes-from)): by default the
  electron-phonon cooling power is `P = Sigma A (Te^3 - Tp^3)` with
  `Sigma = 2.0 W m^-2 K^-3`.

## Examples

Each example below runs as written, and the output shown is what it
printed with absnoise 0.14.0 (Python 3.13, NumPy 2.5.3, SciPy 1.18.1).
Examples 1 to 12 print the same as with 0.13.0.
Values not taken from the shipped recipes -- correlation times,
geometries, noise levels -- are illustrative, not measurements.

### 1. A detector budget for a shipped junction

```python
import numpy as np
from absnoise import RECIPES, SensorBudget, H_PLANCK

budget = SensorBudget(RECIPES[1])        # the Ti/Al/Au junction of Jung et al.
T = 0.3 * budget.recipe.Tc               # operating temperature (K)
tauA = 1e-6                              # occupation correlation time (s), illustrative

dT_A, sums = budget.dT_andreev(T, tauA, t=1.0)    # 1 s of averaging
dT_ph = budget.dT_phonon(T, t=1.0)
print(f"T = {T:.3f} K, thermal time tau_th = {budget.tau_th(T) * 1e9:.2f} ns")
print(f"temperature resolution in 1 s: Andreev {dT_A * 1e6:.2f} uK, "
      f"phonon {dT_ph * 1e6:.2f} uK")
print(f"tauA/C_A = {tauA / sums['C_A']:.3e} s K/J, "
      f"tau_th/C_e = {budget.tau_th(T) / budget.Ce(T):.3e} s K/J")

sigE = budget.energy_resolution(T, tauA)
print(f"energy resolution: {sigE:.3e} J "
      f"(the energy of a {sigE / H_PLANCK / 1e9:.0f} GHz photon)")

f = np.array([0.0, 1e7, 1e8])
nep = budget.nep_spectrum(T, tauA, f)
for fi, a, p in zip(f, nep["andreev"], nep["phonon"]):
    print(f"NEP at {fi:.0e} Hz: Andreev {a:.2e}, phonon {p:.2e} W/sqrt(Hz)")
```

```
T = 0.225 K, thermal time tau_th = 0.78 ns
temperature resolution in 1 s: Andreev 40.48 uK, phonon 3.58 uK
tauA/C_A = 1.172e+15 s K/J, tau_th/C_e = 9.145e+12 s K/J
energy resolution: 1.752e-22 J (the energy of a 264 GHz photon)
NEP at 0e+00 Hz: Andreev 6.26e-18, phonon 5.53e-19 W/sqrt(Hz)
NEP at 1e+07 Hz: Andreev 6.27e-18, phonon 5.53e-19 W/sqrt(Hz)
NEP at 1e+08 Hz: Andreev 6.97e-18, phonon 5.53e-19 W/sqrt(Hz)
```

At this operating point the occupation noise, not the phonon noise,
sets the temperature resolution. That matches the criterion stated in
the code: occupation noise dominates when `tauA/C_A` is larger than
`tau_th/C_e`. In the NEP, the phonon part is flat at every frequency
(`sqrt(4 kB T^2 G)`), and the occupation part rises as
`sqrt(1 + (2 pi f tau_th)^2)`. The energy resolution is computed in
closed form over all frequencies (since 0.11.0; 0.10.x printed
1.810e-22 J here, 3.3 % too high, because its numerical integral
stopped at a finite frequency).

### 2. The temperature bound, finite length, and pair processes

```python
import numpy as np
from absnoise import RECIPES, ShortJunction, FiniteLJunction, sigma_spectrum

r = RECIPES[3]                           # Ti/Al(thick): L = 0.3 um, xi = 3.6 um
T = 0.3 * r.Tc

sj = ShortJunction(r)                    # every channel has the same level
sj.calibrate()
achieved, bound, _ = sj.temperature_bound(T, tauA=1e-6, t_int=1.0)
print(f"short-junction model: achieved / bound = {achieved / bound:.6f}")

fl = FiniteLJunction(r)                  # exact levels at the real length
print(f"finite-length model at phi = 1: achieved / bound = "
      f"{fl.saturation_ratio(1.0, T):.4f}")

# pair processes can only shorten the correlation time
for Gp in (0.0, 1.0, 10.0):
    _, S0, var, tau_eff = sigma_spectrum(0.25, 1.0, Gp, np.array([0.0]))
    print(f"Gp = {Gp:4.1f} Gs: tau_eff = {tau_eff:.4f} / Gs")
```

```
short-junction model: achieved / bound = 1.000000
finite-length model at phi = 1: achieved / bound = 1.0392
Gp =  0.0 Gs: tau_eff = 1.0000 / Gs
Gp =  1.0 Gs: tau_eff = 0.7500 / Gs
Gp = 10.0 Gs: tau_eff = 0.3750 / Gs
```

A short junction whose channels all have the same transparency sits
exactly on the bound. At the real length of this junction the levels
move differently with phase, and the resolution is a little worse than
the bound (never better). `sigma_spectrum` shows that pair processes
(two quasiparticles entering or leaving together, at rate `Gp`) can
only shorten the effective correlation time, so the single-quasiparticle
case is the worst one. Here `Gs` is the single-quasiparticle rate.

### 3. Characterize your junction from a critical-current curve

```python
import numpy as np
from absnoise import ic_model, fit_ic_curve, validate_ic_data, SensorBudget

# Synthetic "measurement": 10 temperatures, 1 % noise, fixed seed.
# The true values (illustrative) are Ic0 = 2 uA, Tc = 0.75 K, tau = 0.78.
T = np.linspace(0.05, 0.65, 10)
rng = np.random.default_rng(1)
sigma = 0.02e-6
Ic = ic_model(T, 2.0e-6, 0.75, 0.78) + rng.normal(0.0, sigma, T.size)

print(validate_ic_data(T, Ic))                   # nothing suspicious
print(validate_ic_data(T * 1e3, Ic)["flags"][:2])  # millikelvin typed as kelvin?
fit = fit_ic_curve(T, Ic, sigma_Ic=sigma)
print(f"Ic0 = {fit.Ic0 * 1e6:.3f} +- {fit.sigma['Ic0'] * 1e6:.3f} uA")
print(f"Tc  = {fit.Tc:.4f} +- {fit.sigma['Tc']:.4f} K")
print(f"tau = {fit.tau:.3f} +- {fit.sigma['tau']:.3f}")
print(f"chi2 = {fit.chi2:.2f} for {fit.chi2_dof} degrees of freedom")

# Package the fit (plus your own measured geometry) for the budgets.
rec = fit.to_recipe(name="my junction", label="mine", xi=5e-6,
                    L=0.2e-6, W=2e-6, Vbg=30.0, Rn=50.0)
b = SensorBudget(rec)
print(f"Josephson inductance at 0.2 K: {b.LJ(0.2) * 1e9:.3f} nH")
```

```
{'n_points': 10, 'flags': []}
[(0, 'T_large', 50.0), (1, 'T_large', 116.66666666666667)]
Ic0 = 2.007 +- 0.012 uA
Tc  = 0.7530 +- 0.0057 K
tau = 0.794 +- 0.032
chi2 = 3.79 for 7 degrees of freedom
Josephson inductance at 0.2 K: 0.226 nH
```

`fit_ic_curve` fits the package's own short-junction model for three
numbers: `Ic0` (the critical current at zero temperature), `Tc` and
`tau`. The error bars are the standard linearized estimate at the
solution (from the slopes of the model at the best fit; accurate when
the errors are small). `chi2` is given when you pass measurement
errors; a value near the number of degrees of freedom (the number of
points minus the three fitted numbers) means model and errors are
consistent. `validate_ic_data` never rescales anything: it returns flags
for values that look like a unit mix-up (above 30 K, or above 1 mA, by
default). `to_recipe` needs your measured geometry, because an Ic(T)
curve does not determine it.

### 4. Plan the measurements before the cooldown

```python
import numpy as np
from absnoise import (plan_ic_measurement, design_ic_temperatures,
                      psd_band_for_tau, plan_psd_measurement,
                      averages_for_tau)

# Expected junction (illustrative): Ic0 = 2 uA, Tc = 0.75 K, tau = 0.78.
Ic0, Tc, tau = 2.0e-6, 0.75, 0.78
sig = 0.02e-6                                  # expected error per point (A)

for top in (0.2, 0.65):                        # highest temperature reached (K)
    T = np.linspace(0.05, top, 8)
    p = plan_ic_measurement(T, sig, Ic0, Tc, tau)
    print(f"run up to {top} K: sigma(Tc) = {p['sigma']['Tc'] * 1e3:.1f} mK, "
          f"covers_tc_knee = {p['covers_tc_knee']}")

cand = np.linspace(0.03, 0.70, 15)             # temperatures the fridge can hold
pick = design_ic_temperatures(cand, 5, sig, Ic0, Tc, tau)
print("measure at:", np.round(np.sort(cand[pick["indices"]]), 3), "K")

f_lo, f_hi, f_knee = psd_band_for_tau(3e-4)    # expect tau near 0.3 ms
print(f"PSD band {f_lo:.1f} to {f_hi:.0f} Hz (knee {f_knee:.0f} Hz)")
f = np.geomspace(f_lo, f_hi, 48)
plan = plan_psd_measurement(f, 64, S0=1e-3, tau_s=3e-4, floor=1e-5)
print(f"64 averages: sigma(tau) = {plan['sigma']['tau_s'] * 1e6:.2f} us")
n, _ = averages_for_tau(3e-6, f, S0=1e-3, tau_s=3e-4, floor=1e-5)
print(f"averages needed for sigma(tau) <= 3 us: {n}")
```

```
run up to 0.2 K: sigma(Tc) = 1031.1 mK, covers_tc_knee = False
run up to 0.65 K: sigma(Tc) = 6.1 mK, covers_tc_knee = True
measure at: [0.03  0.078 0.365 0.413 0.7  ] K
PSD band 132.6 to 2122 Hz (knee 531 Hz)
64 averages: sigma(tau) = 15.56 us
averages needed for sigma(tau) <= 3 us: 1709
```

The planned error bars are the same matrix the fit reports, computed
from the package's own model before any data exist. A run that stays
far below `Tc` barely feels `Tc`; the plan flags this as
`covers_tc_knee = False` (the fit warns about the same case).
`design_ic_temperatures` picks the most informative temperatures one at
a time (a greedy method, not a proof of the best set).
`psd_band_for_tau` returns a band from `f_knee/4` to `4 f_knee` that
satisfies the rules `fit_telegraph_psd` enforces. Error bars fall as
`sqrt(psi'(n_avg))`, which is close to `1/sqrt(n_avg)` for many
averages (see [Example 8](#8-a-noise-spectrum-with-few-averages)), and
`averages_for_tau` inverts that. (0.10.x used `1/sqrt(n_avg)` and
printed 15.50 us and 1708.)

### 5. Occupation noise from a spectrum, and the Allan variance

```python
import numpy as np
from absnoise import (telegraph_traces, psd_single_sided, fit_telegraph_psd,
                      allan_variance, avar_exponential)

# Simulated occupation noise: 64 channels, f = 0.3, tau = 0.2 ms, 5 us steps.
rng = np.random.default_rng(0)
f_occ, tau, dt = 0.3, 2e-4, 5e-6
tr = telegraph_traces(f_occ, tau, dt, 40000, 64, rng)

# Channel-averaged periodogram, then averaged into 40 log-spaced bins.
S = np.mean([psd_single_sided(tr[:, c].astype(float), dt)[1]
             for c in range(64)], axis=0)[1:]
fr = np.fft.rfftfreq(40000, dt)[1:]
edges = np.geomspace(fr[0], fr[-1], 41)
fb, Sb = [], []
for lo, hi in zip(edges[:-1], edges[1:]):
    m = (fr >= lo) & (fr < hi)
    if m.sum() >= 3:
        fb.append(fr[m].mean())
        Sb.append(S[m].mean())

fit = fit_telegraph_psd(np.array(fb), np.array(Sb))
print(f"fitted tau = {fit.tau_s * 1e6:.0f} +- {fit.sigma['tau_s'] * 1e6:.0f} us "
      f"(true {tau * 1e6:.0f} us)")
print(f"variance under the fitted Lorentzian {fit.variance:.4f}, "
      f"trace variance {tr.astype(float).var():.4f}")

# Allan variance of the 64-channel sum against the closed form.
x = tr.astype(float).sum(axis=1)
taus, avar = allan_variance(x, dt, m_list=[20, 80, 320])
closed = avar_exponential(64 * f_occ * (1 - f_occ), tau, taus)
for T, a, c in zip(taus, avar, closed):
    print(f"T = {T * 1e3:.1f} ms: measured {a:.2e}, closed form {c:.2e}")
```

```
fitted tau = 201 +- 2 us (true 200 us)
variance under the fitted Lorentzian 0.2116, trace variance 0.2107
T = 0.1 ms: measured 3.12e+00, closed form 3.13e+00
T = 0.4 ms: measured 5.17e+00, closed form 5.12e+00
T = 1.6 ms: measured 2.33e+00, closed form 2.73e+00
```

`fit_telegraph_psd` fits `S0 / (1 + (2 pi f tau)^2) + floor` in
logarithmic units. The variance under the fitted Lorentzian,
`S0/(4 tau)`, is compared here with the variance of the simulated
traces. The Allan variance of the simulated signal follows the closed
form `avar_exponential`; the longest window has the fewest independent
samples and scatters most (the test suite checks agreement within 12 %
on a longer record).

### 6. Decode a noisy readout trace, with error bars

```python
import os, tempfile
import numpy as np
from absnoise import telegraph_traces, save_trace_csv, load_trace_csv, fit_hmm

# A simulated readout: occupation f = 0.3, tau = 1 ms, sampled every 50 us,
# with Gaussian readout noise 0.3 (levels 0 and 1). Illustrative values.
rng = np.random.default_rng(0)
dt = 5e-5
n = telegraph_traces(0.3, 1e-3, dt, 20000, 1, rng)[:, 0]
y = n + rng.normal(0.0, 0.3, n.size)

path = os.path.join(tempfile.mkdtemp(), "readout.csv")
save_trace_csv(path, dt, y)              # header t_s,y
dt2, y2 = load_trace_csv(path)           # refuses a non-uniform time grid

model, logliks = fit_hmm(y2)
f_hat, tau_hat = model.rates(dt2)
err = model.uncertainties(y2, dt2)["sigma"]
print(f"f = {f_hat:.3f} +- {err['f']:.3f}, "
      f"tau = {tau_hat * 1e3:.3f} +- {err['tauA'] * 1e3:.3f} ms")
print(f"levels {model.means[0]:.3f} and {model.means[1]:.3f}, "
      f"noise {model.sigma:.3f}, EM iterations: {len(logliks)}")
path_hat = model.viterbi(y2)
print(f"samples decoded correctly: {(path_hat == n).mean() * 100:.2f} %")

# The same process read out only every 0.5 ms (half of tau):
m = telegraph_traces(0.3, 1e-3, 5e-4, 20000, 1, rng)[:, 0]
coarse, _ = fit_hmm(m + rng.normal(0.0, 0.3, m.size))
print(f"sampled every 0.5 ms: tau = {coarse.rates(5e-4)[1] * 1e3:.3f} ms "
      f"(first-order formula of 0.10: "
      f"{5e-4 / (coarse.p_up + coarse.p_dn) * 1e3:.3f} ms)")
```

```
f = 0.290 +- 0.020, tau = 0.993 +- 0.056 ms
levels -0.001 and 1.003, noise 0.302, EM iterations: 6
samples decoded correctly: 99.70 %
sampled every 0.5 ms: tau = 1.014 ms (first-order formula of 0.10: 1.284 ms)
```

`fit_hmm` finds the occupation fraction, the correlation time, the two
readout levels and the noise from the trace alone, by
expectation-maximization (a standard method that improves the
estimates step by step until they stop changing). `viterbi` gives the most probable filled/empty
sequence and `posterior` the probability of "filled" at every sample.
The trace goes through the documented file format (`t_s,y`) on the way.

`uncertainties` (new in 0.11.0) gives one-standard-deviation error
bars. They come from how sharply the likelihood (the probability of
the recorded trace under the model) falls off around the best fit --
the standard large-sample estimate, good when the trace holds many
flips. Here the fitted `f` is 0.47 error bars and `tau` 0.12 error
bars away from the true values (0.3 and 1 ms).

The last lines read the same process only every 0.5 ms. Since 0.11.0
the decoder treats the samples as snapshots of a process that runs
in continuous time, which is exact for any sampling interval. The
0.10.x conversion was correct only for samples much closer together
than `tau`, and here it would have reported 1.284 ms instead of
1.014 ms.

### 7. A single photon: matched design, click, and a refusal

```python
from absnoise import (matched_Tc, matched_recipe, SensorBudget, click_template,
                      click_monte_carlo, two_temperature_click, H_PLANCK)

T0 = 0.05                                 # bath temperature (K)
print(f"matched junction Tc for T0 = 50 mK: {matched_Tc(T0):.4f} K")
b = SensorBudget(matched_recipe(T0, W=5.3e-6, L=0.5e-6))   # illustrative geometry
E = H_PLANCK * 26e9                       # one 26 GHz photon

tpl = click_template(b, E, tauA=1e-7, T0=T0)
print(f"peak electron temperature {tpl['Te_peak'] * 1e3:.1f} mK, "
      f"occupation dip {(1 - tpl['ms'].min() / tpl['m0']) * 100:.1f} %")

mc = click_monte_carlo(b, E, tauA=1e-7, T0=T0, n_trials=200, seed=3)
print(f"Monte Carlo: SNR {mc['snr']:.1f}, efficiency {mc['efficiency']:.2f}, "
      f"dark fraction {mc['dark_fraction']:.2f}")

try:
    two_temperature_click(b, E, 1e-7, T0)        # no phonon parameters given
except ValueError as err:
    print("refused:", err)
```

```
matched junction Tc for T0 = 50 mK: 0.0776 K
peak electron temperature 121.8 mK, occupation dip 1.5 %
Monte Carlo: SNR 1.9, efficiency 0.82, dark fraction 0.21
refused: c_ph, kappa_pb and delta_pb carry no vetted defaults: the phonon heat capacity and boundary escape law of a real stack are measurements of that stack. Supply all three from your device's thermal characterization (or a cited equivalent stack).
```

`matched_Tc` chooses the junction `Tc` so that the gap at the operating
temperature is `2.3994 kB T0` (the `MATCHED_RATIO` in the code).
`click_template` follows the electron temperature and the occupations
after the photon is absorbed, without noise; `click_monte_carlo` adds
occupation, phonon and readout noise and detects clicks with a matched
filter. Its numbers are statistics-limited; the tests only check
ordering relations. `two_temperature_click` also lets the local phonons
heat up, but it has no default phonon parameters and refuses to run
without them.

### 8. A noise spectrum with few averages

```python
import numpy as np
from absnoise import (telegraph_traces, psd_single_sided, fit_telegraph_psd)

# One long simulated record: 16 channels summed, f = 0.3, tau = 0.2 ms,
# 5 us steps (illustrative values).
rng = np.random.default_rng(4)
f_occ, tau, dt = 0.3, 2e-4, 5e-6
x = telegraph_traces(f_occ, tau, dt, 4 * 2**15, 16, rng).sum(axis=1)

# Cut it into 4 segments and average their periodograms: n_avg = 4.
segs = x.reshape(4, -1).astype(float)
fr = np.fft.rfftfreq(segs.shape[1], dt)
S = np.mean([psd_single_sided(s, dt)[1] for s in segs], axis=0)
keep = np.arange(1, 3000)                    # every bin up to 18 kHz
f_b, S_b = fr[keep], S[keep]

S0_true = 4 * 16 * f_occ * (1 - f_occ) * tau
plain = fit_telegraph_psd(f_b, S_b)
told = fit_telegraph_psd(f_b, S_b, n_avg=4)
print(f"true S0 = {S0_true:.3e}")
print(f"n_avg not given: S0 = {plain.S0:.3e}, tau = {plain.tau_s * 1e6:.0f} us")
print(f"n_avg = 4:       S0 = {told.S0:.3e} +- {told.sigma['S0']:.1e}, "
      f"tau = {told.tau_s * 1e6:.0f} +- {told.sigma['tau_s'] * 1e6:.0f} us")
print(f"chi2 = {told.chi2:.1f} for {told.chi2_dof} degrees of freedom")
```

```
true S0 = 2.688e-03
n_avg not given: S0 = 2.473e-03, tau = 206 us
n_avg = 4:       S0 = 2.816e-03 +- 1.6e-04, tau = 206 +- 6 us
chi2 = 2961.5 for 2996 degrees of freedom
```

A periodogram scatters a lot, and its logarithm (which the fit uses)
is on average too low: for an average of `n` periodograms by
`psi(n) - ln n`, where `psi` is the digamma function (a standard
special function; `psi'`, the trigamma function, is its derivative
and gives the scatter). For `n = 4` that makes S0 and the floor 12 %
too low. Given `n_avg`, the fit (since 0.11.0) removes this offset
and uses the exact scatter, so `chi2` close to the degrees of freedom
confirms the noise model. `tau` does not depend on the offset. Here
the uncorrected S0 lies 1.3 error bars below the true value and the
corrected one 0.8 error bars above it; with fewer averages the offset
is larger (44 % for `n_avg = 1`). `n_avg` must be the number of
independent periodograms behind every point: do not also average
neighbouring frequencies, or pass that larger number.

### 9. Error bars on the budget from your Ic(T) fit

```python
import numpy as np
from absnoise import ic_model, fit_ic_curve, SensorBudget

# The synthetic Ic(T) measurement of Example 3, fitted the same way.
T = np.linspace(0.05, 0.65, 10)
rng = np.random.default_rng(1)
sigma = 0.02e-6
Ic = ic_model(T, 2.0e-6, 0.75, 0.78) + rng.normal(0.0, sigma, T.size)
fit = fit_ic_curve(T, Ic, sigma_Ic=sigma)

geometry = dict(xi=5e-6, L=0.2e-6, W=2e-6, Vbg=30.0, Rn=50.0)  # yours
T_op, tauA = 0.2, 1e-6                     # operating point (illustrative)

lj = fit.propagate(lambda r: SensorBudget(r).LJ(T_op), **geometry)
print(f"L_J = {lj['value'] * 1e9:.4f} +- {lj['sigma'] * 1e9:.4f} nH")
dT = fit.propagate(lambda r: SensorBudget(r).dT_andreev(T_op, tauA, 1.0)[0],
                   **geometry)
print(f"Andreev resolution in 1 s: {dT['value'] * 1e6:.2f} "
      f"+- {dT['sigma'] * 1e6:.2f} uK")
sE = fit.propagate(lambda r: SensorBudget(r).energy_resolution(T_op, tauA),
                   **geometry)
print(f"energy resolution: {sE['value']:.3e} +- {sE['sigma']:.1e} J "
      f"(nonlinearity {sE['nonlinearity']:.3f})")

try:                                       # 10 mK below the fitted Tc
    fit.propagate(lambda r: SensorBudget(r).LJ(0.743), **geometry)
except ValueError as err:
    print("refused:", err)
```

```
L_J = 0.2264 +- 0.0046 nH
Andreev resolution in 1 s: 44.23 +- 0.79 uK
energy resolution: 1.782e-22 +- 3.2e-24 J (nonlinearity 0.007)
refused: func is too nonlinear over the fit's error bars: the second-order change over one standard deviation is 0.387 times the first-order error bar (limit 0.1). A first-order error bar would mislead; sample (Ic0, Tc, tau) from `cov` and look at the spread of func instead
```

`propagate` (new in 0.12.0) carries the error bars of the Ic(T) fit
into any number the package computes from the fitted junction. You
pass a function of a `Recipe` and your measured geometry. The geometry
stays fixed; only `Ic0`, `Tc` and `tau` vary. The method moves the
three parameters by one standard deviation, in both directions, along
each principal direction of the fit's covariance, and combines the
changes of the result. For a result that depends linearly on the
parameters, this is exactly the usual first-order error bar. The tests
check it against that exact result, against a quadratic closed form,
and against 400 parameter sets drawn at random from the covariance.

`nonlinearity` is the second-order change over one standard deviation,
in units of the error bar. Above 0.1 (the default limit) the call is
refused, because one error bar would then hide a skewed or shifted
spread. Here that happens 10 mK below the fitted `Tc`, where `L_J`
rises steeply. The phonon parameters (`sigma`, `delta`) and the
readout settings of `SensorBudget` carry no error bars here; only the
fit's do.

### 10. Is the readout trace a single-rate telegraph process?

```python
import numpy as np
from absnoise import telegraph_traces, fit_hmm, markov_dwell_test

rng = np.random.default_rng(2)
dt = 5e-5                                  # 50 us samples, noise 0.3

# A Markov occupation: f = 0.3, tau = 1 ms (illustrative values).
n = telegraph_traces(0.3, 1e-3, dt, 20000, 1, rng)[:, 0]
y = n + rng.normal(0.0, 0.3, n.size)
model, _ = fit_hmm(y)
res = markov_dwell_test(y, model, n_boot=500, seed=0)
print(f"Markov trace:    p = {res['p_value']:.2f}, dispersion "
      f"{res['dispersion'][0]:.2f} (empty), {res['dispersion'][1]:.2f} (filled)")

# Filled stays with two time scales: 0.25 ms or 3 ms, equal odds.
parts, s = [], 0
while sum(len(p) for p in parts) < 20000:
    mean = 47 if s == 0 else (5 if rng.random() < 0.5 else 60)
    parts.append(np.full(rng.geometric(1 / mean), s))
    s ^= 1
n2 = np.concatenate(parts)[:20000]
y2 = n2 + rng.normal(0.0, 0.3, n2.size)
model2, _ = fit_hmm(y2)
res2 = markov_dwell_test(y2, model2, n_boot=500, seed=0)
print(f"two time scales: p = {res2['p_value']:.3f}, dispersion "
      f"{res2['dispersion'][0]:.2f} (empty), {res2['dispersion'][1]:.2f} (filled)")
print("complete stays:", res2["n_dwells"])
```

```
Markov trace:    p = 0.50, dispersion 0.96 (empty), 1.11 (filled)
two time scales: p = 0.008, dispersion 0.87 (empty), 2.09 (filled)
complete stays: {0: 207, 1: 207}
```

`fit_hmm` assumes one Markov process: every stay in a state ends with
the same probability per sample, so the stays follow an exponential
law (geometric, counted in samples) with one correlation time. Several
trap states, slow drifts of the rates or a dead time break this, and
then the fitted `f` and `tau` describe no single process.
`markov_dwell_test` (new in 0.12.0) decodes the trace with `viterbi`,
measures every complete stay, and computes the dispersion index
`var / (mean (mean - 1))` of the stay lengths in samples, which is 1
for a geometric law. It compares this number with 500 traces simulated
from the fitted model and decoded in the same way, so decoding errors
(short stays missed in the noise) affect data and reference alike.

In the first trace nothing is found (p = 0.50). In the second, the
filled stays are far too spread out (dispersion 2.09), and the test
rejects the single-rate model. 0.008 is the smallest p-value that 500
simulated traces can give (4/501); ask for more if you need to resolve
smaller values. A large p-value means that no departure was found, not
that the model is proven. `dwell_times` returns the stays themselves,
for a histogram.

### 11. Josephson inductance and Kerr term of a finite-length junction

```python
import numpy as np
from absnoise import RECIPES, JunctionModel, junction_properties, gap_bcs
from absnoise.constants import HBAR, E_CHARGE

r = RECIPES[1]                               # Ti/Al/Au: L = 0.2 um, xi = 5.1 um
out = junction_properties(r, np.array([0.02, 0.225]))
for T, Ic, a1, a3 in zip(out["T"], out["Ic"], out["dIdphi0"], out["a3"]):
    LJ = HBAR / (2 * E_CHARGE) / a1          # Josephson inductance at phi = 0
    print(f"T = {T:.3f} K: Ic = {Ic * 1e6:.3f} uA, L_J = {LJ * 1e9:.4f} nH, "
          f"a3/a1 = {a3 / a1:.4f}")

# the current at phi = pi/2, calculated two independent ways
m = JunctionModel(r)                         # 121 phases; m.phis[60] = pi/2
m.Delta = gap_bcs(0.225, r.Tc, r.Delta)
I_levels = m.current_phase(0.225)[60]        # bound levels + continuum
I_sum = m.cpr_matsubara(0.225, m.phis[60])   # Matsubara sum
print(f"I(pi/2): levels + continuum {I_levels * 1e6:.5f} uA, "
      f"Matsubara sum {I_sum * 1e6:.5f} uA (scale 1)")
```

```
T = 0.020 K: Ic = 2.110 uA, L_J = 0.2196 nH, a3/a1 = -0.0594
T = 0.225 K: Ic = 2.007 uA, L_J = 0.2220 nH, a3/a1 = -0.0640
I(pi/2): levels + continuum 3.67599 uA, Matsubara sum 3.67775 uA (scale 1)
```

The full model of a recipe (`JunctionModel`) counts every transverse
channel at its real length, with the states above the gap (the
continuum). Since 0.13.0 the slope `a1 = I'(0)`, which sets the
Josephson inductance `L_J = (hbar/2e) / a1`, and the cubic term `a3`
of the current-phase relation `I = a1 phi + a3 phi^3 + ...` are exact
for this model (`cpr_series`). They come from the free energy written
as a sum over Matsubara frequencies (P. W. Brouwer and C. W. J.
Beenakker, Chaos, Solitons & Fractals 8, 1249 (1997), Eqs. (6) and (8)). That
sum needs no level search and does not split bound levels from the
continuum, so it stays accurate near zero phase, where bound levels
leave the gap and the two parts nearly cancel. `cpr_matsubara` gives
the whole current-phase relation the same way. The last line compares
it with the level-based `current_phase` at `phi = pi/2`; the two
calculations differ by 5e-4, the size of the energy and phase grid
errors of `current_phase`. 0.12.1 fitted `a1` and `a3` to the
level-based current and printed `a3/a1 = -0.0679` at 0.225 K here (6 %
off) and `L_J = 0.2198 nH` at 20 mK; its `JunctionModel.dIdphi0` was
57 % low for this recipe (see the corrections below).

### 12. The electron-phonon cooling law from a heating curve

```python
import numpy as np
from absnoise import RECIPES, SensorBudget, ep_temperature, fit_ep_coupling

r = RECIPES[1]                               # Ti/Al/Au, area 0.2 um x 1.8 um
P = np.concatenate([[0.0, 0.0], np.geomspace(1e-17, 3e-14, 12)])  # heater (W)
rng = np.random.default_rng(5)               # illustrative "measurement":
Te = ep_temperature(P, 0.03, r.area, 2.5, 3.0) + rng.normal(0.0, 1e-3, P.size)

fit = fit_ep_coupling(P, Te, r.area, Tp_K=0.03, sigma_Te=1e-3, fit_Tp=True)
print(f"Sigma = {fit.Sigma:.2f} +- {fit.sigma['Sigma']:.2f} W m^-2 K^-delta, "
      f"delta = {fit.delta:.3f} +- {fit.sigma['delta']:.3f}, "
      f"Tp = {fit.Tp * 1e3:.1f} +- {fit.sigma['Tp'] * 1e3:.1f} mK")
print(f"chi2 = {fit.chi2:.1f} for {fit.chi2_dof} degrees of freedom")

T = 0.225                                    # operating temperature (K)
G = fit.propagate(lambda s, d: SensorBudget(r, sigma=s, delta=d).Gep(T))
print(f"G_ep at {T} K: {G['value'] * 1e15:.1f} +- {G['sigma'] * 1e15:.1f} fW/K")

low = fit_ep_coupling(P[:5], Te[:5], r.area, Tp_K=0.03, sigma_Te=1e-3)
print(f"five lowest powers only: delta = {low.delta:.1f} +- {low.sigma['delta']:.1f}")
```

```
Sigma = 2.50 +- 0.06 W m^-2 K^-delta, delta = 3.002 +- 0.015, Tp = 29.5 +- 0.6 mK
chi2 = 12.0 for 11 degrees of freedom
G_ep at 0.225 K: 136.2 +- 1.0 fW/K
five lowest powers only: delta = 2.1 +- 1.5
```

Every detector number of `SensorBudget` depends on the cooling law
`P = Sigma A (Te^delta - Tp^delta)`, whose default (`Sigma = 2.0`,
`delta = 3`) follows Lee et al. `fit_ep_coupling` (new in 0.13.0)
measures it on your own device, the way Lee et al. did: from the
electron temperature at several heater powers. It fits `Te` (the
measured quantity) with the closed form
`Te = (Tp^delta + P / (Sigma A))^(1/delta)`, can fit the bath
temperature too, and gives error bars and a chi-square. `propagate`
carries them into any number computed from `(Sigma, delta)`, as
`IcFit.propagate` does for the Ic(T) fit. The last line shows why the
heating range matters: close to the bath the law is nearly linear, and
only the conductance `G` is determined, not `Sigma` and `delta`
apart. The fit warns when every point stays within 10 % of the bath
(in `Te^delta`), and refuses when the exponent runs to the end of its
search range (1 to 8) or the parameters cannot be told apart.

### 13. How much readout bandwidth does the energy resolution need?

```python
from absnoise import RECIPES, SensorBudget

b = SensorBudget(RECIPES[1])                 # Ti/Al/Au
T, tauA = 0.225, 1e-6                        # operating point (illustrative)
S_ro = 5e-15                                 # white readout floor of y (1/Hz)

ideal = b.energy_resolution(T, tauA, S_ro_y=S_ro)
print(f"unlimited band: sigma_E = {ideal:.3e} J")
for f_max in (1e5, 1e6, 1e7, 1e8):           # readout bandwidth (Hz)
    s = b.energy_resolution(T, tauA, S_ro_y=S_ro, f_max=f_max)
    print(f"f_max = {f_max:.0e} Hz: sigma_E = {s:.3e} J ({s / ideal:.3f} x)")
for floor in (S_ro, 0.0):
    f_bw = b.bandwidth_for_resolution(T, tauA, factor=1.1, S_ro_y=floor)
    print(f"bandwidth for 1.1 x the unlimited value, floor {floor:.0e}: "
          f"{f_bw / 1e6:.2f} MHz")
```

```
unlimited band: sigma_E = 2.629e-21 J
f_max = 1e+05 Hz: sigma_E = 1.010e-20 J (3.842 x)
f_max = 1e+06 Hz: sigma_E = 3.639e-21 J (1.384 x)
f_max = 1e+07 Hz: sigma_E = 2.705e-21 J (1.029 x)
f_max = 1e+08 Hz: sigma_E = 2.633e-21 J (1.001 x)
bandwidth for 1.1 x the unlimited value, floor 5e-15: 3.29 MHz
bandwidth for 1.1 x the unlimited value, floor 0e+00: 732.60 MHz
```

`energy_resolution` assumes a matched filter that uses every
frequency. A real readout passes only a band: the resonator
linewidth, the amplifier or the digitizer. With `f_max` (new in
0.14.0) the filter uses only the frequencies from 0 to `f_max`, as an
ideal low-pass readout allows, and the integral
`sigma_E^-2 = 4 int_0^f_max df / NEP^2` still has a closed form.
`bandwidth_for_resolution` returns the bandwidth at which the
resolution comes within a chosen factor of the unlimited value. A
readout floor makes the highest frequencies worthless, so it lowers
the bandwidth you need: 3.3 MHz here with the floor, against 733 MHz
without it, where the thermal time of 0.78 ns sets the scale. The
tests compare the closed form with numerical integration of the NEP
of `nep_spectrum`.

### 14. The exact critical current of the full model

```python
from absnoise import RECIPES, JunctionModel

m = JunctionModel(RECIPES[1])                # Ti/Al/Au, full model, scale 1
Ic_grid = m.Ic(0.02)                         # maximum on the phase grid
Ic_exact, phi_c = m.Ic_matsubara(0.02, return_phase=True)
print(f"Ic(20 mK), scale 1: phase grid {Ic_grid * 1e6:.5f} uA, "
      f"exact {Ic_exact * 1e6:.5f} uA at phi = {phi_c:.4f} rad "
      f"({(Ic_grid / Ic_exact - 1) * 1e4:.1f}e-4)")
scale = m.calibrate(method="matsubara")      # exact Ic(20 mK) = measured Ic20
print(f"scale {scale:.6f} (levels: {RECIPES[1].Ic20 / Ic_grid:.6f}), "
      f"Ic(20 mK) = {m.Ic_matsubara(0.02) * 1e6:.4f} uA")
```

```
Ic(20 mK), scale 1: phase grid 4.07482 uA, exact 4.07668 uA at phi = 1.9698 rad (-4.5e-4)
scale 0.517578 (levels: 0.517814), Ic(20 mK) = 2.1100 uA
```

`JunctionModel.Ic` takes the largest value of `current_phase` on its
phase grid, and `calibrate` scales the model so that this value at
20 mK equals the measured `Ic20`. `Ic_matsubara` (new in 0.14.0)
finds the maximum of the Matsubara-sum current instead (no energy or
phase grid; see Example 11), and `calibrate(method="matsubara")` uses
it. For Ti/Al/Au the two differ by 4.5e-4. For MoRe the level-based
current scatters by up to 3 % between neighboring grid phases near its
maximum, and its `Ic` lies 1.7 % above the exact maximum. Since 0.15.0
the exact maximum is the default of `calibrate` and of
`junction_properties`, so a calibrated model reproduces the measured
`Ic20` exactly. `method="levels"` gives the calibrations of earlier
releases.

## Tutorials

Three end-to-end scripts in the `examples/` folder go from synthetic
"measured" data (fixed seeds, illustrative true values) to a final
result with error bars. Each prints a short summary; the test suite
runs all three and checks every printed number against the output
below.

**`examples/tutorial_detector_budget.py`** (a few seconds): an Ic(T)
curve and a heating curve, fitted with `fit_ic_curve` and
`fit_ep_coupling`, give the junction and its cooling law; the detector
budget built from them gives the thermal time, the phonon NEP and the
energy resolution. The error bars of both fits are carried through
with `propagate` and added in quadrature (the two curves are
independent measurements).

```
Ic(T) fit:   Ic0 = 2.006 +- 0.006 uA, Tc = 0.7514 +- 0.0022 K, tau = 0.793 +- 0.013
heating fit: Sigma = 2.51 +- 0.06 W m^-2 K^-delta, delta = 2.991 +- 0.020, chi2 = 6.6 for 13 dof
budget at T = 0.200 K, tauA = 1 us (illustrative):
  thermal time tau_th = 0.69 +- 0.01 ns
  phonon NEP = 5.19e-19 +- 0.02e-19 W/sqrt(Hz)
  energy resolution = 1.998e-22 J +- 0.014e-22 (Ic fit) +- 0.009e-22 (heating fit)
                    = 1.998e-22 +- 0.016e-22 J, the energy of a 301 +- 2 GHz photon
```

**`examples/tutorial_readout_trace.py`** (about 10 s): a noisy readout
trace of one Andreev level. The decoder gives the occupation and the
correlation time with error bars, `markov_dwell_test` checks the
single-rate model, the noise spectrum of the same record gives a
second estimate of the correlation time, and the occupation-noise
temperature resolution of Ti/Al/Au at 0.3 `Tc` follows for that
correlation time. It scales exactly as `sqrt(tauA)`: Example 1 gives
40.48 uK for 1 us, and 40.48 uK times `sqrt(998)` is 1279 uK. Both
estimates of the correlation time come from the same record, so the
"difference" line is a loose consistency check, not a test.

```
decoder:  f = 0.301 +- 0.011, tauA = 0.998 +- 0.031 ms (EM iterations: 6)
model check: p = 0.92 (dispersion 0.90 empty, 0.93 filled; 637 and 638 stays)
spectrum: tauA = 0.974 +- 0.033 ms (chi2 = 1993 for 2045 dof)
difference decoder - spectrum: 0.5 error bars
Ti/Al/Au at 0.225 K: dT_A in 1 s = 1279.01 +- 20.03 uK for tauA = 0.998 ms
```

**`examples/tutorial_readout_bandwidth.py`** (new in 0.14.0, a few
seconds): how much readout bandwidth a junction calorimeter needs. A
synthetic noise record of the resonator frequency (occupation noise of
Ti/Al/Au at 0.3 `Tc` with `tauA` = 1 us, plus a white readout floor)
is averaged into a spectrum of 32 segments and fitted with
`fit_telegraph_psd(..., n_avg=32)`. The fit gives `tauA`, the floor
and the plateau; the plateau is checked against the budget's own
prediction. The energy resolution then follows for an unlimited band
and for readout bandwidths of 1, 3 and 10 MHz, together with the
bandwidth that brings it within 10 % of the unlimited value. The error
bars come from the fit's covariance of `tauA` and the floor
(`TelegraphPSDFit.cov`, new in 0.14.0), carried through along its
principal axes.

```
spectrum fit: tauA = 1.023 +- 0.021 us, floor = 5.58e-15 +- 0.02e-15 1/Hz, chi2 = 2820 for 2946 dof
plateau S0 = 1.747e-13 +- 0.056e-13 1/Hz, budget prediction 1.682e-13 (1.2 error bars apart)
energy resolution, unlimited band: 2.751e-21 +- 0.043e-21 J
readout bandwidth  1 MHz: 3.735e-21 +- 0.046e-21 J (1.36 times the unlimited value)
readout bandwidth  3 MHz: 3.035e-21 +- 0.044e-21 J (1.10 times the unlimited value)
readout bandwidth 10 MHz: 2.826e-21 +- 0.043e-21 J (1.03 times the unlimited value)
bandwidth for 1.1 times the unlimited value: 3.08 +- 0.03 MHz
```

## What changed in 0.14.0

New:

- `energy_resolution(..., f_max=...)` and `matched_filter_sigma(...,
  f_max=...)`: the matched-filter energy resolution for a readout of
  finite bandwidth, in closed form; and
  `SensorBudget.bandwidth_for_resolution`: the bandwidth that brings
  it within a chosen factor of the unlimited value (Example 13).
- `JunctionModel.Ic_matsubara`: the exact critical current of the full
  finite-length model, from the Matsubara-sum current;
  `calibrate(method="matsubara")` and
  `junction_properties(..., method="matsubara")` use it (Example 14).
  In 0.14.0 the default calibration stayed as it was; since 0.15.0 it
  is the exact one.
- `TelegraphPSDFit.cov`: the covariance of the fitted spectrum
  parameters, with their correlations.
- A third tutorial, on the readout bandwidth (see [Tutorials](#tutorials)).

Corrected (details and before/after numbers under
[Corrections in earlier versions](#corrections-in-earlier-versions)):

- `FiniteLJunction` within `1e-3` rad of zero phase (its finite
  difference step), where the two levels of an Andreev doublet cross:
  the inductance-readout sums were wrong by a factor of up to about
  1000 (saturation ratio 853.79 instead of 1.00747 for Ti/Al(thick) at
  `phi = 1e-4`).

Also new: `FiniteLJunction` refuses an unknown `which` (it was taken as
`"L"`) and a temperature at or above `Tc`; `SensorBudget` refuses a
non-positive `sigma`, `delta`, `L_r` or `nu_r`. The source
distribution now contains `examples/`, so the tutorial tests run from
it too.

Faster, with the same results: `markov_dwell_test` decodes its
simulated traces in batches (identical paths, checked exactly; 2.8 to
3.0 times faster for the 500 traces of Example 10). No function was
renamed or removed, and Examples 1 to 12 print what they printed with
0.13.0.

## What changed in 0.13.0

New:

- `JunctionModel.cpr_series(T)` and `JunctionModel.cpr_matsubara(T)`:
  the exact slope and cubic term of the current-phase relation, and the
  whole relation, of the finite-length model with its continuum, from
  the Matsubara-sum free energy (Example 11). `junction_properties` and
  `JunctionModel.dIdphi0` now use them.
- `fit_ep_coupling(P, Te, area, Tp)` and `EPCouplingFit`: the
  electron-phonon cooling law (`Sigma`, `delta`, optionally the bath
  temperature) from a measured heating curve, with error bars and
  `propagate` (Example 12). `ep_temperature` is its model curve.
- Two tutorials in `examples/` (see [Tutorials](#tutorials)).

Corrected (details and before/after numbers under
[Corrections in earlier versions](#corrections-in-earlier-versions)):

- `JunctionModel.dIdphi0` (44 % low for a zero-length junction, 57 %
  low for Ti/Al/Au) and the `a3` of `junction_properties` (6 % off for
  Ti/Al/Au).
- `abs_energies` at zero phase: it returned one level of each
  degenerate pair, so `FiniteLJunction` at `phi = 0` counted half of
  the occupation heat capacity.
- `ShortJunction.dlnIc_dT` and `dlnI1_dT` far below and just below
  `Tc` (inf, NaN or -inf).
- `sigma_spectrum` for tiny occupations (it stopped with an
  `AssertionError`) and for small rates (it returned zero).
- `continuum_share` below about 0.03 `Tc` (`ZeroDivisionError`, for
  example at 20 mK for Ti/Al/Au) and below 0.002 `Tc` (evaluated at a
  negative temperature).
- `SensorBudget.nep_spectrum` accepted a negative readout floor.

Faster, with the same results: `telegraph_traces` (identical traces;
4 times faster for the 64 channels of Example 5, about 350 times for
one channel), and the level search and free energy of `JunctionModel`
(levels about 2.8 times faster, free energy 2.4 times, `calibrate` 2.4
times; free energy equal to 1e-13). No function was renamed or removed. `JunctionModel.dIdphi0`
ignores its `k` argument (with a `DeprecationWarning`).

## What changed in 0.12.0

New:

- `IcFit.propagate(func, ...)`: error bars on any number computed from
  a fitted junction, from the covariance of the Ic(T) fit
  (Example 9). It refuses when the result is too nonlinear over the
  error bars for one error bar to describe it.
- `markov_dwell_test(y, model)` and `dwell_times(path, dt)`: a check
  that a decoded readout trace follows the single-rate telegraph model
  the decoder assumes, calibrated by simulating the fitted model
  (Example 10).

Corrected (details and before/after numbers under
[Corrections in earlier versions](#corrections-in-earlier-versions)):

- `gap_bcs` between 0.99 `Tc` and `Tc` (it returned 0 from 0.9995 `Tc`
  on, and was up to 0.7 % off just below).
- `avar_exponential` for averaging times much shorter than the
  correlation time (it lost all accuracy, and turned negative at
  `T = 1e-6 tau`).
- `ic_model` and `ShortJunction.Ic` at and above `Tc` (NaN, now 0), and
  so `plan_ic_measurement` with a planned temperature above `Tc`.
- The sign of `absnoise.levels.continuum_free_energy`.
- `plan_psd_measurement` with a non-integer `n_avg`.

Also new: `click_template` refuses a negative deposited energy, a
non-positive `tauA` or `T0`, and accepts a zero deposit (it raised
`ZeroDivisionError`); `avar_exponential` refuses non-finite inputs.
No function was renamed or removed, and results that were correct
before are unchanged (all eight earlier examples print exactly what
they printed with 0.11.0).

## What changed in 0.11.0

New:

- `SensorBudget.energy_resolution` is computed in closed form over all
  frequencies (`matched_filter_sigma`); 0.10.x integrated numerically
  over a band that ended too early and overestimated the resolution
  (see "Corrections" below). `method="grid"` reproduces the old
  numbers.
- `TelegraphHMM.uncertainties(y, dt)`: error bars for everything
  `fit_hmm` returns (Example 6); also `score` (the exact gradient of
  the likelihood) and `loglik`.
- `flip_probabilities(f, tauA, dt)`: the exact per-sample flip
  probabilities, now used by both `telegraph_traces` and the decoder.
- `fit_telegraph_psd(..., n_avg=n)` corrects the known offset of a
  log-averaged periodogram and uses its exact scatter; the PSD planner
  uses the same scatter (Example 8).
- `fit_hmm` starts from the best two-group split of the samples, which
  works for traces the old median split could not handle, and the
  decoder runs about 10 times faster.

The corrections, with before/after numbers, are listed under
[Corrections in earlier versions](#corrections-in-earlier-versions).

## What is in the package

**Junctions and levels**

- `Recipe`, `RECIPES` -- a junction design; the six measured recipes of
  Jung et al. (index 0 to 5: Ta/Ti/Au, Ti/Al/Au, Ti/Al(thin),
  Ti/Al(thick), Ti/Nb/NbN, MoRe).
- `abs_energies(phi, tau, c, Delta)` -- the exact Andreev levels of one
  channel of finite length, found by solving a closed-form equation.
- `gap_bcs(T, Tc, Delta0)` -- the gap `Delta(T)`, from the gap equation
  of BCS theory (the standard Bardeen-Cooper-Schrieffer theory of
  superconductivity) solved numerically, not the common interpolation
  formula. Since 0.12.0 it follows the square-root closing of the gap
  right up to `Tc`.
- `ShortJunction` -- the closed-form short-junction model: levels,
  current-phase relation, `Ic`, `dIdphi0`, occupation sums and
  `temperature_bound`.
- `FiniteLJunction` -- the same sums from the exact finite-length
  levels, and `saturation_ratio` (1 means on the bound).
- `JunctionModel`, `junction_properties` -- the full mode-summed model
  of a recipe, including states above the gap (the continuum),
  calibrated to the measured `Ic` at 20 mK. Since 0.13.0,
  `JunctionModel.cpr_series` (exact `I'(0)` and phase-cubic term) and
  `JunctionModel.cpr_matsubara` (the current-phase relation) evaluate
  the same model as a sum over Matsubara frequencies; since 0.14.0,
  `JunctionModel.Ic_matsubara` gives its exact critical current, and
  `calibrate` and `junction_properties` take `method="matsubara"`.
- `continuum_share`, `occupation_heat_capacities` -- how much of the
  current, and of the occupation heat capacity, comes from the
  continuum.

**Graphene thermal model** (`materials`)

- `carrier_density`, `fermi_energy`, `dos_ef`, `n_modes` -- electron
  density from the gate voltage, and quantities derived from it.
- `heat_capacity`, `ep_power`, `gth` -- electron heat capacity,
  electron-phonon cooling power and thermal conductance.
- `steady_temperature` -- the electron temperature at which cooling
  balances a steady absorbed power (readout self-heating).
- `fit_ep_coupling`, `EPCouplingFit`, `ep_temperature` (0.13.0): fit
  the cooling law to a measured heating curve; the fit result (with
  `G` and `propagate`); the model curve.

**Occupation noise**

- `channel_generator`, `sigma_spectrum` -- the four-state master
  equation of one channel (the equations for how the probabilities of
  its four occupation states change in time), with single and pair
  processes, and its exact spectrum.
- `tau_activated`, `noneq_penalty` -- an energy-dependent (thermally
  activated) correlation time, and the resolution penalty of an excess,
  non-equilibrium occupation.
- `telegraph_traces`, `psd_single_sided`, `flip_probabilities` --
  Monte Carlo telegraph signals (exact at any time step), their
  periodogram, and the exact per-sample flip probabilities.
- `allan_variance`, `avar_exponential`, `avar_white` -- the Allan
  variance of data, and its closed forms for exponentially correlated
  and white noise (`avar_exponential` is accurate for any averaging
  time since 0.12.0).

**Device budget and single photons**

- `SensorBudget` -- one recipe as a detector: `dT_andreev`,
  `dT_phonon`, `freq_noise_spectrum`, `energy_resolution`,
  `energy_resolution_analytic_A`, `nep_spectrum`, plus `Ce`, `Gep`,
  `tau_th`, `LJ` and `participation` (the share of the resonator's
  inductance that the junction supplies). Since 0.14.0
  `energy_resolution` takes a readout bandwidth `f_max`, and
  `bandwidth_for_resolution` returns the bandwidth a target
  resolution needs.
- `matched_filter_sigma(P, Q, W, tau_th, tauA, f_max=None)` -- the
  closed-form energy resolution behind `energy_resolution`, from the
  three noise levels of `nep_spectrum` (phonon, occupation, readout)
  at zero frequency, for an unlimited band or up to `f_max`.
- `matched_Tc`, `matched_recipe` -- the matched-level design.
- `click_template`, `click_monte_carlo` -- the response to one photon,
  noiseless and with noise.
- `two_temperature_click`, `total_energy` -- the same with a heated
  local phonon bath, and its conserved energy.

**From measured data**

- `fit_ic_curve`, `IcFit`, `ic_model` -- the Ic(T) fit, its result
  (with `to_recipe`, and `propagate` for error bars on anything
  computed from the fitted junction), and the model curve.
- `fit_telegraph_psd`, `TelegraphPSDFit`, `telegraph_psd_model` -- the
  noise-spectrum fit, its result (with the covariance `cov` since
  0.14.0), and the model spectrum.
- `fit_hmm`, `TelegraphHMM` -- the readout-trace decoder
  (`posterior`, `viterbi`, `forward_backward`, `rates`, `from_rates`,
  and `uncertainties`, `score`, `loglik` for error bars).
- `markov_dwell_test`, `dwell_times` -- does a decoded trace follow the
  single-rate telegraph model? And the stays in each state.
- `load_ic_csv`, `save_ic_csv` (header `T_K,Ic_A` or
  `T_K,Ic_A,sigma_Ic_A`), `load_trace_csv`, `save_trace_csv` (header
  `t_s,y`), `validate_ic_data` -- the file formats and unit checks.

**Planning**

- `plan_ic_measurement`, `design_ic_temperatures` -- planned Ic(T)
  error bars, and which temperatures to measure.
- `psd_band_for_tau`, `plan_psd_measurement`, `averages_for_tau` -- the
  frequency band, the planned spectrum-fit error bars, and the number
  of averages for a target error on `tau`.

**Constants** -- `HBAR`, `KB`, `E_CHARGE`, `H_PLANCK`, `PHI0`, `EPS0`
(CODATA 2018) and `V_F` (graphene Fermi velocity, 1.0e6 m/s).

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

## When it refuses, and why

`absnoise` raises an error instead of guessing when:

- `two_temperature_click` is called without all three phonon parameters
  (`c_ph`, `kappa_pb`, `delta_pb`): they are measurements of a real
  device, and none are shipped. It also refuses a non-positive `c_ph`
  or `delta_pb`, a negative `kappa_pb` or a negative deposited energy.
- `fit_telegraph_psd` gets a spectrum that changes by less than a
  factor of 4 across the band, or whose fitted knee lies outside
  `[2 f_min, f_max/2]`: such a band cannot determine `tau`. It also
  needs at least 8 points, positive frequencies (drop the zero-frequency
  point), a positive PSD and, if given, `n_avg >= 1`.
- `plan_psd_measurement` and `averages_for_tau` get a planned band the
  fit would refuse (the same rules, applied in advance), `n_avg < 1`
  or not finite, or a non-positive target; `psd_band_for_tau` gets a margin below 4.
- `fit_ic_curve` gets fewer than 4 points, fewer than 5 points without
  measurement errors (the noise level must then be estimated from the
  residuals), non-finite values, non-positive temperatures or currents,
  or non-positive errors. It raises `RuntimeError` if the fit does not
  converge or the data cannot separate the three parameters. If the
  data stay below 0.3 `Tc`, it warns (it does not refuse) that `Tc`
  and `tau` are weakly determined.
- `plan_ic_measurement` and `design_ic_temperatures` get `tau` outside
  (0, 1], a non-positive `Ic0`, `Tc`, temperature or error, a number of
  picks outside 3 to the number of candidates, or candidates that
  cannot determine the parameters even all together.
- a data file has the wrong header or number of fields, fewer than 2
  rows, temperatures that do not strictly increase, or (for a trace) a
  time grid that is not uniform; the file is refused, not repaired.
  `validate_ic_data` refuses non-finite values, non-positive
  temperatures and negative currents.
- `allan_variance` gets fewer than 4 samples, non-finite samples, a
  non-positive `dt` or a window outside `1 <= m <= len(x)//2`; the
  closed forms get a negative variance or PSD, or non-positive times
  (`avar_exponential` also refuses non-finite inputs).
- `TelegraphHMM` gets flip probabilities outside (0, 1), a non-positive
  noise or not exactly two levels.
- `steady_temperature` gets a negative power; `nep_spectrum` gets a
  negative or non-finite frequency; `click_template` and
  `two_temperature_click` get a scenario other than `"C"` or `"T"`.
- a `SensorBudget` method gets a temperature that is not positive, or
  one at or above `Tc` where the gap has closed (only `Ce`, `Gep`,
  `tau_th` and `dT_phonon` work above `Tc`), a non-positive `tauA` or
  averaging time, a `which` other than `"L"` or `"I"`, or a negative
  readout floor. (Up to 0.10.x these gave a `ZeroDivisionError` or a
  silent NaN.) `matched_Tc` gets a non-positive `T0`.
- `telegraph_traces` and `flip_probabilities` get `f` outside [0, 1]
  or a non-positive `tauA` or `dt`. (Up to 0.10.x `telegraph_traces`
  instead refused coarse time steps, with a Python `assert`.)
- `IcFit.propagate` finds that the result is too nonlinear over the
  error bars (second-order change over one standard deviation above
  `max_nonlinearity` times the error bar, default 0.1), that a
  parameter would leave its range within one standard deviation
  (`Ic0 > 0`, `Tc > 0`, `0 < tau <= 1`), or that `func` returns a
  non-finite value.
- `markov_dwell_test` gets fewer than `min_dwells` (default 30)
  complete stays in a state, a mean stay below 2 samples (sampling too
  coarse), fewer than 19 simulated traces, or a model that is not a
  `TelegraphHMM`.
- `click_template` gets a negative or non-finite deposited energy, or
  a non-positive `tauA` or `T0`.
- `fit_hmm` gets fewer than 4 samples, non-finite samples or a
  constant trace; `forward_backward` and `viterbi` get an empty trace; `rates` gets flip probabilities adding up to 1 or
  more (successive samples anti-correlated: no telegraph process does
  that); `uncertainties` finds that the likelihood is not at a
  maximum or does not determine all five parameters.
- `fit_ep_coupling` (0.13.0) gets arrays of different lengths,
  non-finite values, a negative power, a non-positive temperature,
  area or error, too few points (the number of fitted parameters plus
  1, or plus 2 without errors), no heated point, or no point above the
  bath. It refuses a fit whose exponent ends at the edge of its search
  range (1 to 8) and data that cannot separate the parameters, and it
  warns when every point stays within 10 % of the bath in `Te^delta`.
  `EPCouplingFit.propagate` refuses in the same cases as
  `IcFit.propagate`. `ep_temperature` refuses a negative power and
  non-positive parameters.
- `sigma_spectrum` (0.13.0) gets `f` outside [0, 1], a non-positive or
  non-finite `Gs`, or a negative `Gp`.
- `ShortJunction.dlnIc_dT` and `dlnI1_dT` (0.13.0) get a temperature
  that is not positive or is at or above `Tc`; `continuum_share` gets a
  temperature outside (0, `Tc`). `SensorBudget.nep_spectrum` refuses a
  negative or non-finite readout floor (0.13.0), as `energy_resolution`
  already did.
- `energy_resolution` and `matched_filter_sigma` (0.14.0) get an
  `f_max` that is not finite and positive, or an `f_max` with
  `method="grid"`; `bandwidth_for_resolution` gets a `factor` that is
  not finite or not above 1.
- `JunctionModel.Ic_matsubara` (0.14.0) gets a temperature that is not
  positive, a model without a gap, fewer than 5 grid phases or a phase
  tolerance outside (0, 0.1); `calibrate` and `junction_properties` get
  a `method` other than `"levels"` or `"matsubara"`.
- `FiniteLJunction.andreev_sums` and `saturation_ratio` (0.14.0) get a
  `which` other than `"I"` or `"L"` (it was taken as `"L"`) or a
  temperature that is not positive or is at or above `Tc` (a NaN ratio
  up to 0.13.0); `SensorBudget` gets a `sigma`, `delta`, `L_r` or
  `nu_r` that is not finite and positive.

## How the results are checked

267 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.9 with the
oldest NumPy (1.24.0) and SciPy (1.10.0) the package allows. The
numerical checks compare the package with an exact formula, a second
calculation done a different way, or a simulation. A few checks
compare with stored numbers on purpose, for example the recipe table
(locked to its source) and one tabulated value of the BCS gap. Others
check that refusals fire. The main checks:

**Levels and gap**

- Short-junction limit: the exact solver gives
  `Delta sqrt(1 - tau sin^2(phi/2))` to 1e-12 (relative to `Delta`),
  for 6 transparencies and 4 phases.
- Fully transparent, finite length: the levels satisfy the known
  exact level condition for that case (the Kulik condition) to
  1e-10 rad.
- One clean channel at 1 mK: `Ic = 2 e Delta / hbar` within 0.5 %.
- BCS gap: equals `Delta0` at 0.01 `Tc` to 1e-10, is below 0.05
  `Delta0` at 0.99999 `Tc`, equals the tabulated 0.956887 at `Tc/2`
  within 2e-4, and falls monotonically.
- BCS gap near `Tc` (0.12.0): between 0.99 and 0.9998 `Tc` it equals
  direct solutions of the gap equation in integral form to 1e-9, and
  closer to `Tc` direct solutions of its Matsubara-sum form; the two
  forms agree to 1e-9 from 0.1 to 0.999 `Tc`. At 0.99999 and 0.999999
  `Tc` it matches the Ginzburg-Landau limit
  `u = (pi t / A) sqrt(8 ln(t1/t) / (7 zeta(3)))`, derived from the
  sum, to 2e-5. It is continuous at 0.99 `Tc` (1e-6), and below it
  returns the values of 0.11.0 (two checked to 1e-12).
- The free energy of the continuum (0.12.0): `F_bound + F_cont` stays
  continuous (within 2 % of the bound-level jump) where a bound level
  leaves for the continuum, and `continuum_free_energy` agrees with the
  `JunctionModel` code path to 5e-4.
- Short junction: the continuum carries no phase dependence (below
  1e-14).
- Levels at zero phase (0.13.0): the degenerate pair is returned twice,
  equal to the levels at `phi = 1e-6` within 1e-6; near zero phase
  (1e-9 to 1e-5) the levels satisfy the level condition with the angle
  `2 asin(sqrt(tau) sin(phi/2))` computed in the test to 1e-12 of that
  angle, and their splitting equals the first-order closed form within
  1e-5; `FiniteLJunction` gives the same `C_A` at `phi = 0` and
  `1e-7` (1e-9); away from the degeneracy the levels are identical to
  those of 0.12.1. For `tau = 1` the level near zero energy close to
  `phi = pi` is found (first-order closed form, 1e-6).
- Matsubara-sum current (0.13.0): for zero length it equals the
  closed-form short-junction current at every phase to 1e-9 of its
  maximum, at 0.027, 0.3 and 0.9 `Tc`, and its slope and cubic term
  equal their closed forms (derived in the test) to 1e-9 and 1e-8. At
  finite length (Ti/Al(thick)) it agrees with the level-plus-continuum
  `current_phase`, an independent calculation, within 2e-3 of the
  maximum away from zero phase, and its slope with an odd polynomial
  fitted to that current within 3e-3. The Legendre chi function used
  for the tail equals its series (1e-13) and `pi^2/8` at 1.
  `junction_properties` returns the closed-form slope and cubic term
  (1e-9, 1e-8; 0.12.1 was more than 1 % off in the cubic term), and
  `dIdphi0` the closed-form slope (0.12.1: 44 % low).
- Free energy speed-up (0.13.0): `free_energy_multi` and
  `free_energy_components` equal the 0.12.1 phase-by-phase loop (kept
  in the tests) to 1e-13 relative.
- Finite-length sums near zero phase (0.14.0): for `|phi|` below the
  difference step of 1e-3 rad, the inductance-readout responsivity and
  saturation ratio and the current-readout noise sum equal a
  calculation from the analytic level slopes and curvatures (implicit
  differentiation of the level condition along each branch, coded in
  the test) to 1e-6, and `C_A` to 1e-12. The 0.13.0 ratio was more
  than 10 times too large there. From `|phi| = 1e-3` on, the level data
  are identical to those of 0.13.0.
- Exact critical current (0.14.0): for zero length `Ic_matsubara`
  equals the closed-form maximum of the short-junction current at
  0.027 `Tc` to 1e-9 (its phase to 1e-6 rad), and the maximum over
  200 001 phases at 0.5 `Tc` to 1e-9. At finite length it agrees with
  the level-based `Ic` within 2e-3 and is never below the Matsubara
  current on a 40-phase grid. After `calibrate(method="matsubara")` it
  returns the measured `Ic20` (1e-12). In the `junction_properties`
  table with that method, the slope and the cubic term differ from the
  default only by the calibration scale (1e-12), and `Ic` is the
  `Ic_matsubara` value (1e-12).

**Occupation noise and the bound**

- A uniform short junction reaches the Cauchy-Schwarz bound to 1e-9;
  a finite-length junction (Ti/Al(thick), two phases) never goes below
  it (to 1e-9).
- The occupation responsivity (how much the signal changes per
  kelvin) matches numerical differentiation to 1e-5, and `dIdphi0` matches its closed form to 1e-12.
- Master equation: without pair processes the variance is `2f(1-f)`
  and the correlation time `1/Gs`, both to 1e-12; the variance does not
  depend on the pair rate (1e-12); the correlation time never rises
  with the pair rate; probability is conserved (1e-14); the spectrum
  integrates to the variance within 5e-4.
- Telegraph Monte Carlo: the plateau within 10 %, half the plateau at
  the knee within 0.12, and the variance of 50-second averages equals
  `S(0)/(2t)` within 25 %.
- Exact sampling: the flip probabilities equal the matrix exponential
  of the two-state rate matrix (computed by SciPy) to 1e-9 relative,
  for occupations 0.02 to 0.9 and time steps from 1e-6 to 7 `tau`, and
  reduce to the first-order values for tiny steps (1e-5); traces
  sampled at `dt = tau` have autocorrelation `e^-k` at lags 1 to 3
  within 0.01 (the first-order rule would give 0).
- Periodogram: `sum(S) df` equals the variance of the data to 1e-12
  (Parseval), for even and odd lengths.
- `sigma_spectrum` (0.13.0): for `f` down to 1e-200 and at 0 and 1
  the variance, plateau and correlation time equal the closed form of
  two independent spins (1e-12; 0.12.1 stopped with an
  `AssertionError`); scaling all rates by 1e-20 to 1e15 rescales the
  result exactly (1e-12; 0.12.1 returned 0 below 1e-14); for ordinary
  inputs it equals the 0.12.1 code (kept in the tests) to 1e-12.
- Telegraph speed-up (0.13.0): the traces are identical to those of the
  0.12.1 step-by-step loop (kept in the tests) for eight sets of
  occupations, time steps and shapes, including block boundaries.
- `occupation_heat_capacities`: the bound part equals the level sum
  `C_A` of `FiniteLJunction` within 1e-5 (two code paths); the
  continuum part is exactly 0 at zero length; for Ti/Al/Au at 0.3 `Tc`
  it is positive and below 10 % of the total.

**Allan variance**

- Linear drift: the estimator gives `(c T)^2 / 2` to 1e-12.
- The closed form matches the defining integral, rewritten without
  cancellation and integrated adaptively, to 1e-9 from `T = 1e-8 tau`
  to `3 tau`, and its Taylor series `2r/3 - r^2/2 + 7r^3/30 - r^4/12`
  (`r = T/tau`) to 1e-12 at `r = 1e-3` (0.12.0).
- The closed form matches numerical integration of the defining
  integrals to 1e-6, and both limits to 1e-3.
- Telegraph Monte Carlo matches the closed form within 12 %, white
  noise matches `S0/(2T)` within 5 %.

**Device budget**

- The frequency-noise spectrum is exactly half the plateau at the knee
  and 1/101 of it at ten times the knee (1e-12).
- The phonon resolution equals its closed form (1e-18 K); the
  occupation resolution respects the bound; the full energy resolution
  is never better than 0.99 times the occupation-only closed form.
- NEP: the phonon part equals `sqrt(4 kB T^2 G)` at every frequency to
  1e-14; the occupation part has exactly the shape
  `sqrt(1 + (2 pi f tau_th)^2)` (1e-12); the zero-frequency ratio of
  the two equals the documented criterion (1e-12) and doubles when
  `tauA` doubles (1e-12); the readout part rises faster.
- Energy resolution: the closed form equals adaptive numerical
  integration of `4 / NEP^2` over all frequencies to 1e-8 (two code
  paths), with and without a readout floor, for both readouts `"L"`
  and `"I"` and two correlation times; the bare closed form matches
  the same integration to 1e-9 for five sets of noise levels, and two
  exact limits (occupation noise only; readout floor only) to 1e-14.
  The 0.10.x band-limited integral (`method="grid"`) still equals the
  NEP integral on its own grid to 1e-10, and lies 3.3 % and 10.7 %
  above the closed form for `tauA` = 1 us and 1 ns.
- Out-of-range temperatures, correlation times and averaging times
  are refused, and (0.13.0) a negative readout floor in `nep_spectrum`.
- Readout bandwidth (0.14.0): the closed-form band integral equals
  adaptive quadrature to 1e-10 for real, complex, equal and nearly
  equal roots of the quartic, and without the quartic term. The
  band-limited resolution equals quadrature of `4 / NEP^2` from
  `nep_spectrum` to 1e-8, for both readouts with and without a floor,
  and falls as the bandwidth grows; at 1 Hz it equals
  `NEP(0) / sqrt(4 f_max)` (1e-12). The 0.10.x grid integral equals the
  closed form over its own band (1e-6). Without a floor the bandwidth
  for a factor `F` equals the closed form
  `sqrt(a / b) tan(pi / (2 F^2)) / (2 pi)` (1e-9); with a floor the
  resolution at the returned bandwidth is `F` times the unlimited value
  (1e-10).
- Log slopes (0.13.0): `ShortJunction.dlnIc_dT` at 0.99, 0.999 and
  0.9999 `Tc` equals a five-point derivative of `ic_model` (a separate
  implementation of Ic(T)) within 5e-3 (0.12.1: -inf above 0.998
  `Tc`); at 10 mK for MoRe both slopes are finite and below 1e-12
  (0.12.1: inf and NaN); at 0.5 `Tc` the value is the 0.12.1 one,
  bit for bit.
- `continuum_share` at 0.027 `Tc` returns a finite current share and
  NaN for the derivative share (0.12.1: `ZeroDivisionError`).

**Recipes and thermal model**

- The label, `Tc`, `tau`, `Ic20` and `Rn` of all six recipes equal
  Table I of Jung et al. exactly.
- Gate voltage 20 V gives the quoted density 1.7e16 m^-2 within 5 %;
  the heat capacity at 1 um^2, 1.7e16 m^-2 and 0.1 K lies between 4 and
  8 `kB`.
- `steady_temperature` round trip: `ep_power(steady_temperature(P))`
  returns `P` within `max(1e-12 P, 1e-30 W)`.
- Cooling-law fit (0.13.0): noise-free data return `Sigma`, `delta`
  (1e-8) and the bath temperature (1e-7); a Nelder-Mead search on a
  chi-square coded in the test finds the same minimum (1e-6); with
  `delta` and `Tp` fixed the error bar of `Sigma` equals the closed form
  `[sum (dTe/dSigma / sigma)^2]^(-1/2)` (1e-6); over 150 noisy fits the
  scatter of `Sigma` and `delta` is between 0.8 and 1.2 times the
  reported error bars and their means are within 4 standard errors of
  the true values; `G` equals `SensorBudget.Gep` with the fitted
  parameters (1e-14) and a numerical derivative of `ep_power` (1e-8);
  `propagate` is exact for a linear function (1e-10) and for `Sigma^2`
  (1e-9), and uses the marginal block when `Tp` is fitted; refusals
  and the near-linear warning fire.

**Single photons**

- Matched design: `Delta(T0) / (kB T0) = 2.3994` to 1e-9, both from
  the gap equation and inside the `SensorBudget` of the matched recipe
  (0.10.x gave 2.39954 there).
- The starting electron temperature satisfies
  `gamma (T_pk^2 - T0^2) / 2 = E_gamma` to 1e-12 (this checks the
  starting point, not the time integration).
- The noiseless template dips and returns within 0.5 % of its start;
  the Monte Carlo has signal-to-noise above 1 and efficiency above the
  dark fraction.
- Two-temperature model: with nothing deposited, both temperatures
  stay exactly at the start value (not even rounding drift); with no
  escape, energy is conserved within 1e-5 and both temperatures end at
  the independently computed common temperature within 1e-5; with a
  very large phonon heat capacity it reproduces `click_template`
  within 1e-3; with escape on, both temperatures shed more than 99 %
  of their excursion.

**Fits and decoder**

- Ic(T) model: the closed-form zero-temperature maximum equals a dense
  grid maximum within 1e-8 (transparencies 0.05, 0.3, 0.78, 0.99 and
  1); `ic_model`
  equals `ShortJunction.Ic` within 1e-12 for all six recipes; both are
  exactly 0 at and above `Tc`, and a planned point above `Tc` leaves
  the planned error bars unchanged (1e-12).
- Error propagation (`IcFit.propagate`): exact for a linear function
  of the parameters (1e-10, also the fit's own error bars and
  covariance block); for `Tc^2` the closed form `2 Tc sigma(Tc)`
  (1e-9) and the expected range of the nonlinearity measure; for the
  Josephson inductance, the spread of 400 parameter sets drawn from
  the covariance within 12 % (the standard error of that spread is
  3.5 %); refusals for a nonlinear function, `tau` within one error
  bar of 1, and non-finite results.
- Dwell-time check: stays counted by hand on a short path; the
  dispersion index of 200 000 geometric samples within 0.02 of 1; the
  simulated reference paths have the closed-form occupation and flip
  rate (1 % and 3 %); on 12 Markov traces (`telegraph_traces` plus
  noise, each fitted with `fit_hmm`) no p-value is below 0.01 and their
  mean lies between 0.3 and 0.75; a trace whose filled stays have two
  time scales gives dispersion above 1.4 and p-value at most 0.01.
- Ic(T) fit: noise-free data return `Ic0` and `Tc` within 1e-6 and
  `tau` within 1e-5; over 60 noisy fits the scatter is between 0.5
  and 1.5 times the reported error bars; a fitted recipe fed back into
  `ShortJunction` reproduces the curve within 1e-4.
- Spectrum fit: noise-free data return `S0` and `tau` within 1e-6 and
  the floor within 1e-4; on telegraph Monte Carlo, `tau` within 15 %,
  the variance within 15 % of the trace variance, `S0` within 20 %.
- Spectrum fit with `n_avg`: the offset correction changes `tau` by
  less than 1e-7 and scales `S0` and the floor by exactly
  `exp(ln n - psi(n))` (1e-7); on 300 spectra with the exact scatter
  of an average of two periodograms, the mean of `ln(S0_fit/S0)` is
  within 3 standard errors of 0 (0.10.x: -0.27), the floor within
  10 % (it keeps a bias of about -5 % from the fit's nonlinearity),
  the mean `chi2/dof` within 0.06 of 1 (0.10.x: 1.29), and the scatter
  of `tau` within 20 % of the planned error bar.
- Decoder: above 99.9 % of samples correct at noise 0.15; at noise 0.5
  the stated confidence matches the real accuracy within 0.01; the fit
  finds `f` and `tau` within 4 statistical errors and the levels and
  noise within 0.02; the log-likelihood never decreases (within 1e-6
  relative). `rates` inverts `from_rates` to 1e-9, also at
  `dt = 2 tau`.
- Decoder speed-up: the new recursions give the same posteriors,
  transition probabilities and log-likelihood as the 0.10.1 matrix
  code (kept in the tests) to 1e-12, and the same Viterbi path.
- Decoder error bars: the exact gradient equals finite differences of
  the log-likelihood to 1e-5; at high signal-to-noise the error bars
  of the flip probabilities equal the binomial closed form
  `sqrt(p (1 - p) / N)` within 0.5 %; over 30 simulated traces the
  errors of `tau` and `f`, divided by each trace's own error bar,
  have a standard deviation between 0.65 and 1.4 and a mean within
  0.6 of zero; at `dt = tau/2` the fit finds `tau` within 4 error bars
  (the 0.10.x conversion is more than 20 % off on the same fit).
- Decoder robustness: with one 200-sigma outlier (0.10.1 returned NaN)
  the log-likelihood equals a log-space forward recursion to 1e-10 and
  the fit runs with a non-decreasing likelihood; a single-sample trace
  matches its closed form to 1e-14.
- Decoder start: a noise-free trace filled 90 % of the time (which
  0.10.x could not start on) returns both levels to 1e-6 and the exact
  state sequence; the two-group split equals a brute-force search over
  all thresholds to 1e-12.
- Spectrum-fit covariance (0.14.0): over 400 spectra with the exact
  scatter of an average of four periodograms, the observed
  correlations of `S0`, `tau` and the floor agree with `cov` within
  0.12, and the observed scatter with its diagonal within 15 %;
  `sqrt(diag(cov))` equals `sigma` (1e-12).
- Dwell test speed-up (0.14.0): the batched decoder returns the same
  paths as `viterbi` (compared exactly, also with a far outlier and
  with exact ties), and `markov_dwell_test` the same p-values and
  dispersions as the 0.13.0 loop (kept in the tests).
- Files: Ic files round-trip exactly; trace files round-trip the
  samples exactly and `dt` within 1e-18 s.
- Tutorials (0.13.0, the third since 0.14.0): every script in
  `examples/` runs in the tests,
  and every printed line matches the output given in the README and
  in the script, each number within two units of its last printed
  digit.

**Planning**

- Planned Ic(T) error bars equal the fit's within 5 %; 250 simulated
  cooldowns scatter in `Tc` within 20 % of the plan; a run below 0.3
  `Tc` is flagged and its `Tc` error is more than 10 times larger; the
  greedy design is never worse than 20 random choices.
- Planned spectrum-fit error bars equal the fit's within 0.1 %, also
  for non-integer `n_avg` (1.5, 2.5, 7.3; 0.12.0); the
  `sqrt(psi'(n_avg))` scaling holds to 1e-9; `averages_for_tau` returns the
  smallest sufficient number; 300 simulated spectra scatter in `tau`
  within 20 % of the plan; a band from `psd_band_for_tau` is accepted
  by the fit and returns `tau` within 1e-4.

## Corrections in earlier versions

**0.15.0 (this release) corrects the default calibration.** Up to
0.14.0, `JunctionModel.calibrate` and `junction_properties` set the
scale so that the level-based `Ic`, the largest value of
`current_phase` on its phase grid, equals the measured `Ic20`. For MoRe
that value lies 1.7 % above the exact maximum of the same model, so
the exact critical current of a calibrated MoRe model at 20 mK was
1.7 % below the measured value it was calibrated to (0.98344 `Ic20`).
Both now use the exact maximum (`Ic_matsubara`) by default, and the
exact critical current equals `Ic20` (1e-9 in the test). The scales
of the shipped recipes, before and after:

| Recipe | 0.14.0 (`method="levels"`) | 0.15.0 (default) | Change |
|---|---|---|---|
| Ta(10)/Ti(60)/Au(5) | 0.383553 | 0.382872 | -0.18 % |
| Ti(6)/Al(60)/Au(5) | 0.517814 | 0.517578 | -0.05 % |
| Ti(6)/Al(70) | 0.554064 | 0.553802 | -0.05 % |
| Ti(6)/Al(200) | 0.658375 | 0.658100 | -0.04 % |
| Ti(6)/Nb(5)/NbN(50) | 0.133312 | 0.133331 | +0.01 % |
| MoRe(50) | 0.492841 | 0.501142 | +1.68 % |

Every quantity of a calibrated `JunctionModel` scales with it. Of the
README examples only Example 11 changes, in the last digit of `L_J`
(0.2195 to 0.2196 nH and 0.2219 to 0.2220 nH). Pass
`method="levels"` to reproduce 0.14.0. The calibration of MoRe now
takes about 27 s instead of 4 s (the cost of the exact maximum grows
as `1/T`); for the other recipes it took 1 to 10 s either way.
`ShortJunction.calibrate` was always exact and is unchanged. Tests:
`tests/test_v015.py`.

**0.14.0 corrects one result of 0.13.0.** It was found
by a test that fails on 0.13.0 and passes now (`tests/test_v014.py`).
It does not change the output of Examples 1 to 12.

- **Finite-length sums near zero phase.** `FiniteLJunction` takes the
  slopes and curvatures of the levels by central differences with a
  step of 1e-3 rad, and matched the levels at the three phases by
  their order in energy. At `phi = 0` the two levels of an Andreev
  doublet cross. For `|phi|` below the step, the stencil therefore
  switched from one level to the other in the middle: the slope of a
  doublet level came out as `-alpha phi / dphi` instead of `-alpha`,
  and its curvature as about `-2 alpha / dphi` (`alpha` is the slope
  of the doublet splitting). For Ti/Al(thick) at 0.3 `Tc` and
  `phi = 1e-4`, `saturation_ratio(..., "L")` was 853.79 instead of
  1.00747 (the value just outside the step). Each level is now followed
  along its own branch there. `C_A` was not affected, and from
  `|phi| = 1e-3` on nothing changed.

Also: the module docstring of `absnoise.levels` said that the junction
current counts bound states only; the code has included the continuum
since 0.1.0, and the docstring now says so.

**0.13.0 corrects six results of 0.12.1.** Each was
found by a test that fails on 0.12.1 and passes now
(`tests/test_v013.py`). None of them changes the output of Examples 1
to 10.

- **The slope of the full model at zero phase.**
  `JunctionModel.dIdphi0` fitted a cubic to the first 8 points of
  `current_phase`. The first point is a one-sided difference, about 14
  times too large, and the next points scatter by up to 17 % because
  bound levels leave the gap there and the bound and continuum parts
  nearly cancel. For a junction of zero length, where the exact slope
  is known in closed form, the result was 44 % low; for Ti/Al/Au at
  0.3 `Tc` it was 1.229e-6 instead of 2.865e-6 A/rad (before the
  calibration scale). `junction_properties` used a two-term fit over
  0 to 0.76 rad (`cpr_coeffs`) instead: its slope was within 5e-4,
  but its cubic term `a3` was 2.3 % off for zero length and 6.1 % for
  Ti/Al/Au at 0.3 `Tc` (-1.0077e-7 instead of -9.495e-8 A/rad^3; slope
  1.484052e-6 instead of 1.483368e-6 A/rad). Both now come from the
  exact Matsubara sum (`cpr_series`). `cpr_coeffs` itself is unchanged.
- **Levels at zero phase.** At `phi = 0` the two levels of each
  Andreev pair coincide, and `abs_energies` returned them once. The
  same happened for every phase below about `2e-8 / sqrt(tau)`,
  because the angle `arccos(1 - tau + tau cos(phi))` was taken of a
  number already rounded to 1. `FiniteLJunction` at `phi = 0` then
  counted half of the occupation heat capacity (Ti/Al(thick) at 0.3
  `Tc`: 1.7417e-21 instead of 3.4835e-21 J/K). The angle now comes from
  `2 asin(sqrt(tau) sin(phi/2))` near zero phase (and from a matching
  form near `pi`, which also finds the level near zero energy of a
  transparent channel close to `phi = pi`), and equal levels are kept.
  Away from those points the levels are identical to 0.12.1.
- **Log slopes of the short junction.** `ShortJunction.dlnIc_dT` and
  `dlnI1_dT` used a fixed step of 2e-3 `Tc`. Below that temperature
  the lower point was negative (MoRe at 10 mK: inf and NaN instead of
  0); within 2e-3 `Tc` of `Tc` the upper point had no supercurrent
  (Ti/Al/Au at 0.999 `Tc`: -inf instead of -1320 1/K). The step is
  now at most `T/2` and `(Tc - T)/10`; between 0.004 and 0.98 `Tc` it
  is unchanged. Between 0.98 and 0.998 `Tc` the values change because
  the old step was too coarse there (at 0.997 `Tc`: -532.6 before,
  -443.6 now, against -442.1 from a fine derivative of `ic_model`; the
  new step leaves a difference error below 0.4 % up to `Tc`).
- **Master equation with tiny occupations or slow rates.**
  `sigma_spectrum` divided by the square root of the stationary
  probabilities, which underflows for `f` below about 1e-162, and
  stopped with an `AssertionError` there and at `f = 0` or 1. It also
  dropped every mode with a rate below 1e-14 in the chosen unit, so
  with `Gs = Gp = 1e-15` it returned `S0 = 0` and `tau_eff = 0`
  (instead of `0.75 / Gs`). It now builds the symmetrized rate matrix
  from its closed-form entries and removes only the stationary mode.
- **`continuum_share` at low temperature.** Below about 0.03 `Tc`
  `dI/dT` is below rounding, and the share of `dI/dT` was 0/0: 0.12.1
  stopped with `ZeroDivisionError` (Ti/Al/Au at 20 mK, for example).
  Below 0.002 `Tc` its temperature step reached negative temperatures.
  It now returns NaN for that share, and the step is limited as for
  the log slopes.
- **Negative readout floor in `nep_spectrum`.** It was accepted and
  gave a NaN readout channel and a total below the other channels;
  it is refused now, as in `energy_resolution`.

**0.12.0 corrects five results of 0.11.0.** Each was
found by a test that fails on 0.11.0 and passes now
(`tests/test_v012.py`). None of them changes the output of an example
in this README.

- **The gap just below `Tc`.** `gap_bcs` read the gap from a cubic
  spline in `T/Tc` that ended at 0.9995 `Tc` and returned 0 above it.
  The gap closes as a square root, which such a spline cannot follow.
  Between 0.9995 `Tc` and `Tc` the result was 0 instead of 0.039 to
  0.006 `Delta0` (at 0.9997 `Tc`: 0 instead of 0.0307 `Delta0`), and
  `SensorBudget` refused those temperatures as "at or above Tc".
  Below, the spline was off from the solution of its own gap equation
  by +0.13 % at 0.998 `Tc` (0.077956 instead of 0.077854 `Delta0`)
  and by -0.67 % at 0.999 `Tc` (0.054884 instead of 0.055252; see
  [Limits](#limits) for how this equation differs from exact BCS
  theory this close to `Tc`). From 0.99 `Tc` on, the gap now comes from a table of
  `u / sqrt(1 - t/t1)`, built from the Matsubara-sum form of the same
  gap equation. It agrees with direct solutions to better than 1e-9.
  Below 0.99 `Tc` nothing changed.
- **Allan variance for short averaging times.** The four terms of the
  closed form of `avar_exponential` cancel to order `(T/tau)^3`, so
  for `T` much shorter than `tau` the result was mostly rounding
  error. With `var = tau = 1`: at `T = 1e-4` 6.66578e-5 instead of
  6.66617e-5, at `T = 1e-5` 7.77e-6 instead of 6.67e-6 (17 % high),
  and at `T = 1e-6` -1.1e-4 (negative) instead of 6.67e-7. Below
  `T = tau/2` the bracket is now summed as its Taylor series; the
  tests compare it with the defining integral, rewritten so that it
  does not cancel, to 1e-9.
- **Critical current above `Tc`.** `ic_model` and `ShortJunction.Ic`
  returned NaN at and above `Tc` (0/0 with a closed gap); the
  supercurrent there is 0. Through NaN sensitivities,
  `plan_ic_measurement` then called any plan with one temperature
  above `Tc` "not identifiable" (identifiable False, no error bars).
  Now such a point adds no information, and the plan equals the plan
  without it.
- **Sign of the continuum free energy.** The function
  `absnoise.levels.continuum_free_energy` returned the opposite sign
  of the continuum free energy that `JunctionModel` (and so the
  current, `continuum_share` and `occupation_heat_capacities`)
  uses. Its docstring said the sign was fixed by continuity of the
  total free energy where a bound level leaves for the continuum. With
  the old sign the total jumped there by twice the free energy of the
  departing level (2.013 `Delta` instead of 0.007 `Delta` across
  0.06 rad, in the test case). `JunctionModel` was right and is
  unchanged.
- **Planned spectrum error bars for a non-integer `n_avg`.**
  `plan_psd_measurement` truncated `n_avg` to an integer, although
  the fit accepts an effective number of averages. For the band of
  Example 4 with `n_avg = 1.5`, it planned `sigma(tau)` = 159.0 us
  (the value for `n_avg = 1`) while the fit reports 119.9 us. It now
  uses `n_avg` as given.

**0.11.0 corrects five results of 0.10.1.** Each
change below alters numbers that earlier versions printed.

- **Energy resolution.** `SensorBudget.energy_resolution` integrated
  numerically between two cut-off frequencies. Above the upper cut-off
  the integrand still carries weight (it falls only as `1/f^2` while
  occupation noise dominates, and stays flat where phonon noise
  dominates), so the result was too large: by 3.3 % in Example 1
  (1.810e-22 J before, 1.752e-22 J now), by 10.7 % for Ti/Al/Au at
  0.3 `Tc` with `tauA` = 1 ns, and by 18 % with `which="I"` there.
  The integral now has a closed form. `method="grid"` gives the old
  numbers.
- **Correlation time from sampled traces.** `TelegraphHMM.rates`
  turned the per-sample flip probability `s = p_up + p_dn` into
  `tau = dt / s`, which is right only for samples much closer than
  `tau`. For a process running in continuous time and read every
  `dt`, the exact relation is `tau = -dt / ln(1 - s)`. The old value
  was too long by 2.5 % at `dt = 0.05 tau` and by 27 % at
  `dt = 0.5 tau`. `f` was not affected. The simulator
  `telegraph_traces` used the matching first-order flip
  probabilities, so the package agreed with itself but not with a
  real, continuously running process. It now samples the continuous
  process exactly, and any time step is allowed. Examples 5 and 6
  therefore print slightly different numbers (Example 6 printed
  `f = 0.297, tau = 0.977 ms` before).
- **Spectrum fit with `n_avg`.** The fit treated the logarithm of an
  averaged periodogram as unbiased with variance `1/n_avg`. Its mean
  is in fact low by `psi(n) - ln n` and its variance is `psi'(n)`, so
  `S0` and the floor were low by a factor 0.56 (n = 1), 0.76 (n = 2),
  0.88 (n = 4) or 0.992 (n = 64), and `chi2` high by 64 % (n = 1)
  or 0.8 % (n = 64). `tau` is unchanged. The planner now uses the
  same scatter: Example 4 prints 15.56 us and 1709 averages (15.50 us
  and 1708 before).
- **Matched design.** `matched_Tc` solved the gap equation with
  `Delta0 = 1.7639 kB Tc`, while `Recipe` (and so the budget of
  `matched_recipe`) uses 1.764. The budget then saw
  `Delta(T0)/(kB T0) = 2.39954` instead of the stated 2.3994;
  `matched_Tc(0.05)` moves from 0.0776088 K to 0.0776059 K.
- **Periodogram.** `psd_single_sided` doubled the Nyquist bin (the
  highest frequency, for an even number of samples), so the spectrum
  held slightly more than the variance of the data.

Also: `SensorBudget` methods refuse temperatures at or above `Tc`
(they raised `ZeroDivisionError`), zero or negative temperatures, and
non-positive `tauA` (a negative value gave NaN); and the decode module
docstring promised likelihood-based error bars that did not exist
until `uncertainties`.

**0.10.1 fixes four problems in 0.10.0.**

- **NumPy 1.x.** The package allows NumPy 1.24 and newer, but four
  integrals called `numpy.trapezoid`, which exists only from NumPy 2.0.
  On NumPy 1.x, `energy_resolution`, the current and free-energy
  methods of `JunctionModel` (`Ic`, `calibrate`, ...),
  `junction_properties`, `continuum_share` and
  `occupation_heat_capacities` stopped with an `AttributeError`. They
  now work on both, and NumPy 1.24.0 with SciPy 1.10.0 is tested.
- **Current readout in `energy_resolution`.** With `which="I"` it used
  the phase 1e-4 rad instead of the phase of maximum current that
  `freq_noise_spectrum`, `nep_spectrum` and
  `energy_resolution_analytic_A` use. Its result disagreed with the
  NEP-based value (it was 66 % larger without a readout floor in the
  example setup of the tests). The default `which="L"` was not affected.
- **Full transparency.** `ic_model` returned NaN at `tau = 1`, which
  the fit's bounds and the planner allow; `plan_ic_measurement` then
  reported such a junction as not identifiable. It now uses the
  correct limit, `Delta0/2`.
- **`n_avg` below 1.** `fit_telegraph_psd` accepted `n_avg = 0` or a
  negative value and returned infinite or NaN error bars. It now
  refuses, as the planner already did.

Also: `TelegraphHMM` and `fit_hmm` were importable from `absnoise` but
missing from `__all__`; they are now listed. Docstrings that described
tests that do not exist (a Monte Carlo check of the pair-process
model, a convergence-order check of the two-temperature integrator, a
continuity check of the continuum sign, a comparison of the telegraph
simulation with the occupation sums) now say what is and is not
tested, and the spectrum-fit notes now give the floor's real test
tolerance (1e-4, not 1e-6).

**What the 0.10.0 README overstated.** It said the checks run on
Python 3.9-3.14 (Python 3.10 was not in the CI matrix), gave "every
spectrum" a Monte-Carlo check (the master-equation spectrum has none),
called the Ic(T) covariance "exact" (it is the linearized estimate),
said the tests lock "every recipe number" (they lock `Tc`, `tau`,
`Ic20`, `Rn` and the label, not the geometry and gate voltage), said
no check uses a stored number (the recipe table and one gap value do),
described a machine-precision "energy conservation of the click
integrator" (the check covers the starting temperature only), and
said the file loaders flag unit mix-ups (that is `validate_ic_data`).

No earlier release records a correction. The full history is in
[CHANGELOG.md](CHANGELOG.md).

## Limits

- The correlation time `tauA` is an input, not a prediction.
- The device budget (`SensorBudget`) uses the short-junction model and
  converts the signal to resonator frequency through the Josephson
  inductance at zero phase. `which="I"` is accepted and evaluated at
  the phase of maximum current, but the frequency conversion is still
  that inductance.
- The finite-length occupation sums (`FiniteLJunction`) count bound
  levels only; `occupation_heat_capacities` gives the size of the
  continuum part at one phase (the tests check its sign and size only
  for Ti/Al/Au).
- `JunctionModel.current_phase` (and so `Ic`, `calibrate` and the
  `Ic` of `junction_properties`) still uses bound levels plus the
  continuum on an energy grid, differentiated on the phase grid. It
  differs from the exact Matsubara sum (`cpr_matsubara`) by about 5e-4
  of `Ic` away from zero phase, and its two end points (`phi = 1e-3`
  and `pi - 1e-3`, one-sided differences) are off by 1 to 2 % of `Ic`
  (zero-length test case). Its maximum lies inside the grid for the
  shipped recipes. Use `cpr_matsubara` where these errors matter. At
  20 mK `Ic` lies 4.2e-4 to 1.8e-3 below the exact maximum
  (`Ic_matsubara`) for four recipes and 1.5e-4 above it for Ti/Nb/NbN.
  For MoRe (the largest `L/xi`, 0.43) the level-based current scatters
  by up to 3 % between neighboring grid phases near its maximum, and
  `Ic` lies 1.7 % above the exact maximum. Since 0.15.0 `calibrate` and
  `junction_properties` use the exact maximum by default, so this
  error no longer enters the calibration; with `method="levels"` it
  does.
- The cost of `Ic_matsubara` grows as `1/T`. At 20 mK it took about
  1.5 s for Ti/Al/Au (the level-based `Ic` about 4 s) and about 20 s
  for MoRe (2.1e5 Matsubara terms for each of 159 modes; the
  level-based `Ic` 3.4 s) on the development machine.
- `FiniteLJunction` differentiates the levels with a phase step of
  1e-3 rad. Near `phi = pi` for a fully transparent channel
  (`tau = 1`, where a level touches zero energy) the levels have a
  kink, and the derivatives there are not meaningful; none of the
  shipped recipes has `tau = 1`.
- `JunctionModel.cpr_coeffs` is a two-term fit over 0 to 0.76 rad and
  keeps the bias of the neglected higher terms (2 to 6 % in `a3`);
  `cpr_series` gives the exact coefficients.
- The Matsubara sum runs to 2000 times the gap and adds the rest of
  its leading term in closed form; the remainder is about 1e-10
  relative. Its cost grows as `1/T`.
- `fit_ep_coupling` fits one power law with one bath temperature. It
  takes the absorbed power as exact (put the uncertainty of your
  heater calibration into the result yourself), and the error bars are
  first order: inspect `cov` when `delta` is poorly determined.
- The continuum occupation *noise* is not modelled: it would need a
  correlation-time model with material parameters this package does
  not invent (its thermodynamic weight is computed).
- The phonons are one lumped temperature; a non-thermal phonon model has
  no cited parameters at these scales. The two-temperature model needs
  your measured phonon parameters.
- `click_monte_carlo` results are statistics-limited; the tests check
  only ordering relations.
- `fit_hmm` models Gaussian readout noise with the same width in both
  states, and assumes the upper readout level means "filled". Its
  error bars (`uncertainties`) are the large-sample estimate: the
  tests check their calibration on traces with a few hundred flips;
  with only a handful of flips they are rough. A far outlier no longer
  breaks the computation, but the Gaussian model widens its fitted
  noise to absorb it (one 200-sigma sample in 4000 raised the fitted
  noise from 0.20 to 0.70 in a test trace): remove glitches first.
- `fit_telegraph_psd` can remove the periodogram offset only when you
  give `n_avg`, one number for all points (without it, `S0` and the
  floor are low by `exp(psi(n) - ln n)` for `n` averages). Even with
  it, the fitted floor kept a bias of about -5 % at `n_avg = 2` in the
  tests (-1 % at 8). The fit uses the continuous-time Lorentzian; the
  spectrum of a sampled trace also contains folded-back (aliased)
  noise near the highest frequency, which the fit does not model.
- `SensorBudget.energy_resolution` assumes an ideal matched filter over
  all frequencies, or (`f_max`, 0.14.0) over the frequencies up to
  `f_max` with nothing above: an ideal low-pass readout. A real
  readout attenuates gradually and may add noise of its own shape;
  model that with `S_ro_y` only where its noise is white.
- `gap_bcs` solves the gap equation with `Delta0 / (kB Tc) = 1.7639`,
  the BCS value `pi e^-gamma = 1.76388` rounded. Its solution would
  close at 1.000013 `Tc` instead of at `Tc`. `gap_bcs` returns 0 from
  `Tc` on, so the gap steps from 0.0063 `Delta0` to 0 there. Close to
  `Tc` this small shift matters: the gap lies above the curve of the
  exact constant by 0.65 % at 0.999 `Tc`, 1.3 % at 0.9995 `Tc` and
  6 % at 0.9999 `Tc` (and by less than 3e-6 below 0.5 `Tc`). The
  constant is kept so that every result below 0.99 `Tc` stays as in
  earlier releases.
- `IcFit.propagate` is first order. It carries only the uncertainty of
  `Ic0`, `Tc` and `tau`; your geometry, the phonon parameters and the
  readout settings enter as exact numbers.
- `markov_dwell_test` looks at one property of the stays, their spread
  (the dispersion index). It detects several time scales and dead
  times, but not, for example, a correlation between successive
  stays. The simulated reference uses the fitted model as it is (the
  model is not refitted to each simulated trace), and the first and
  last stay of a trace are left out. With `n_boot` simulated traces
  the smallest p-value it can report is `4 / (n_boot + 1)`.
- The greedy temperature design is a good heuristic, not a proof of
  the best subset.

## Where it comes from

Methodological basis:

> T. M. Mahim, A. S. M. Mohsin and M. M. Rahman, "Andreev occupation
> noise sets the sensitivity limit of proximity Josephson thermal
> detectors"; code for the paper:
> https://github.com/Tanvir-Mahmud-Mahim/andreev-occupation-noise

This package is the general-purpose engine; the paper repository
reproduces the specific study.

Six junction recipes ship with full provenance (Table I of W. Jung et
al., Phys. Rev. Applied 26, 014078 (2026), arXiv:2503.06850);
graphene electron-phonon cooling follows G.-H. Lee et al., Nature
586, 42 (2020); physical constants are CODATA 2018. The phonon NEP is
compared with J. C. Mather, Appl. Opt. 21, 1125 (1982). The
matched-filter energy resolution is the frequency-domain optimal-filter
result `(int_0^inf 4 df / NEP^2)^(-1/2)` of S. H. Moseley, J. C. Mather
and D. McCammon, J. Appl. Phys. 56, 1257 (1984); its band-limited form
(0.14.0) stops the integral at the readout bandwidth. The Allan
variance follows D. W. Allan, Proc. IEEE 54, 221 (1966) and W. J.
Riley, Handbook of Frequency Stability Analysis, NIST Special
Publication 1065 (2008). The temperature design follows F. Pukelsheim,
Optimal Design of Experiments, SIAM (2006). The Matsubara-sum free
energy of `cpr_series` and `cpr_matsubara` follows from Eqs. (6) and (8) of P. W. Brouwer
and C. W. J. Beenakker, Chaos, Solitons & Fractals 8, 1249 (1997); the
tests check it against the closed-form short-junction current. The
heating-curve fit of `fit_ep_coupling` is the analysis of Lee et al.
(their Fig. 4c and Table I: `delta = 3`, `Sigma` = 2.04 to 3.30
W m^-2 K^-3). Near `Tc` the gap equation
is solved in its Matsubara-sum form, which follows from the integral
form through the partial-fraction expansion
`tanh(x/2) / x = sum_n 4 / (x^2 + (2n+1)^2 pi^2)`; the tests check the
two forms against each other. A replacement recipe number
must arrive with a new source. No phonon-stack parameters are shipped
at all -- they are measurements of a real device, and calls without
them raise.

## Citing, support and license

If `absnoise` helps your work, please cite it with the concept DOI
[10.5281/zenodo.22048608](https://doi.org/10.5281/zenodo.22048608),
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/absnoise/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.
