absnoise

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)
• Tutorials (three end-to-end scripts in examples/)
• What changed in 0.14.0
• What changed in 0.13.0
• 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.
• 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): 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

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.

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

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

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?

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

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

Corrected (details and before/after numbers under 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).

Corrected (details and before/after numbers under 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):

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

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