absnoise

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). 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
• Install
• Units and conventions
• Examples (each with the output it prints)
• What changed in 0.12.0
• What changed in 0.11.0
• What is in the package
• When it refuses, and why
• How the results are checked
• Corrections in earlier versions
• Limits
• Where it comes from
• Citing, support and license

A short guide to the words used here

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.

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): 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.12.0 (Python 3.13, NumPy 2.5.3, SciPy 1.18.1). Values not taken from the shipped recipes -- correlation times, geometries, noise levels -- are illustrative, not measurements.

1. A detector budget for a shipped junction

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

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

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

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

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

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

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

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

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?

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.

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

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

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

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).
• matched_filter_sigma(P, Q, W, tau_th, tauA) -- the closed-form energy resolution behind energy_resolution, from the three noise levels of nep_spectrum (phonon, occupation, readout) at zero frequency.
• 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, 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.

How the results are checked

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

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

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

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.
• Files: Ic files round-trip exactly; trace files round-trip the samples exactly and dt within 1e-18 s.

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.12.0 (this release) 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 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.

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).
• 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; a real readout with a finite bandwidth does worse.
• 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 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). 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, which always resolves to the latest version; every release is archived on Zenodo. 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 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. 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.