Quantum dots that move when nobody's looking: implementing the joint PAT + pulsed-gate fit from arXiv:2604.26947

A walkthrough of the central spectroscopic technique in Benson et al.'s "anomalous photon-assisted tunneling" paper, where a top-gate voltage that is supposed to leave singlet-triplet splittings alone visibly shifts them by hundreds of micro-electronvolts. We rebuild the model, simulate Tien-Gordon PAT side bands across a Ge/SiGe-inspired double quantum dot, and run the same joint PAT + pulsed-gate fit the authors use to pin down the linear shift coefficient.

Paper. Jared Benson, C. E. Sturner, A. R. Huffman, Sanghyeok Park, Valentin John, Brighton X. Coe, Tyler J. Kovach, Stefan D. Oosterhout, Lucas E. A. Stehouwer, Francesco Borsoi, Giordano Scappucci, Menno Veldhorst, Benjamin D. Woods, Mark Friesen, M. A. Eriksson — Large quantum dot energy level shifts in anomalous photon-assisted tunneling. arXiv:2604.26947 (cond-mat.mes-hall, posted 2026-04-29).

Why Python

I picked Python — specifically NumPy plus scipy.special.jv and scipy.optimize.least_squares — for one concrete reason: the technique is a Bessel-weighted multi-photon Floquet sum on top of small dense-Hermitian eigenproblems, and that's the exact ecosystem scipy was written to make boring. In Tien-Gordon, the absorption strength of the nth side band is $J_n^2(\alpha_{ac})$, and scipy.special.jv gives me arbitrary-order Bessel functions to double precision out of the box. The two-electron DQD Hamiltonian factors into independent 2×2 spin blocks, so I never need a sparse solver; I write the lower eigenvalue in closed form and let NumPy broadcast across the whole (eps, V_T) grid. The joint fit is two parameters and a smooth residual — least_squares with default Trust-Region-Reflective converges in a handful of iterations, and the Jacobian it returns gives me 1-σ error bars from (JᵀJ)⁻¹ for free.

A Rust or C re-implementation would mean dragging in a Bessel implementation (statrs, cephes, or hand-rolling the recurrence) and writing my own Levenberg-Marquardt; the time budget is better spent on getting the physics right and producing plots a reader can stare at. Lifetimes, comptime, and memory ordering have nothing to offer here — there is no parallelism worth chasing, the Hilbert space is 4-dimensional, the 2D map fits in 50 KB, and the whole pipeline runs in under a second. Python is the path of least resistance, the standard for this niche, and (when I am honest) the language the authors themselves almost certainly used to draft their fit.

What this paper says, in my own words

A double quantum dot in a Ge/SiGe heterostructure has a singlet-triplet splitting at the (0,2) charge configuration — the energy required to promote one of the two electrons to the next orbital so that the two-particle state is a triplet rather than a singlet. That orbital-splitting energy is set by how strongly the dot is squeezed laterally, which in turn is set by the plunger gate that creates the dot. Conventional wisdom — and most effective-Hamiltonian theories used in qubit operation papers — assumes that the top gate, which sits over the whole device, just shifts the chemical potential rigidly and does not change the orbital structure of any individual dot. Under that wisdom, you can sweep the top gate freely without re-tuning the singlet-triplet splitting.

Benson and collaborators report that this is wrong. They scan their photon-assisted-tunneling resonance condition $h f = E_T - E_S$ as a function of the top-gate voltage $V_T$, and the resonance line tilts visibly in the $(\varepsilon, V_T)$ plane: it bends, in some cases enough to disappear off the bottom of the map. They cross-check the bend with a second experiment, pulsed-gate spectroscopy, which probes the splitting through an avoided-crossing population transfer instead of microwave absorption. Both data sets agree on the same linear coefficient: $\Delta_{ST}^{(0,2)}(V_T) = \Delta_0 + \alpha\,V_T$ with $\alpha$ of the order of micro-electronvolts per millivolt — small in absolute terms but enormous compared to the conventional expectation of zero. They show the same behaviour in two different dots of the same device after retuning, ruling out a per-dot oddity.

The technique I reproduce is therefore not "we built a Ge/SiGe heterostructure" — that I cannot do — and it is not "we computed from first principles why the heavy-hole admixture causes the linear shift" — that I do not attempt either. The technique is the measurement-fitting glue: a two-modality joint reduction that pins down $(\Delta_0, \alpha)$ tighter than either modality alone, in a regime where PAT is degenerate and pulsed-gate is precise but uses a different lever arm. This is the part a quantum-dot characterization lab would actually copy, and it is small enough to fit on a single page.

Tien-Gordon, in one paragraph

Photon-assisted tunneling, in its 1963 Tien-Gordon form, is a Floquet identity. If a state at energy $E$ is driven by a sinusoidal voltage $V(t) = V_{ac}\cos(\omega t)$ that couples linearly to its detuning, the wavefunction picks up a time-dependent phase

The Jacobi-Anger expansion turns the inner exponential into a sum over Bessel functions

so the driven state is a superposition of replicas at energies $E + n\hbar\omega$ with weights $J_n(\alpha_{ac})$. For a tunneling rate from a static-bias source $\Gamma_0(\omega_e)$, the AC-driven rate becomes

In a double quantum dot, the analogue is: at each detuning $\varepsilon$ where $E_T(\varepsilon) - E_S(\varepsilon) = n\hbar\omega$, microwave-driven inelastic transitions light up. The PAT signal is a sum of Lorentzians anchored to those resonance conditions, weighted by $J_n^2(\alpha_{ac})$. That is the entire spectroscopic content. The Bessel weights satisfy $\sum_n J_n^2 = 1$ — every absorbed quantum of energy has to come from somewhere — and the standard "weak drive" regime is $\alpha_{ac}\lesssim 1$, where only $|n|=0,1$ matter; turning the drive up brings $|n|=2,3,\ldots$ into view. I verify this sum rule and the Bessel-weight ratio in the demo at the end.

The double-quantum-dot model

I work in the two-electron sector near the (0,2)–(1,1) charge boundary, with the (2,0) configuration far enough away to drop. The four states I keep are

The singlet block couples $|S(0,2)\rangle$ to $|S(1,1)\rangle$ via a tunneling matrix element $t_S$. The triplet block couples $|T(0,2)\rangle$ to $|T_0(1,1)\rangle$ via $t_T$ — which is generally different from $t_S$ because the orbital overlap is different. With detuning $\varepsilon$ defined so that $\varepsilon > 0$ pushes (0,2) below (1,1), the singlet 2×2 is

and the triplet 2×2 is

The constant $\Delta_{ST}^{(0,2)}$ is the (0,2) singlet-triplet splitting — the energy cost of promoting one electron to the next orbital — and is the quantity that Benson et al. find depends on the top gate. The two blocks are independent, so the gap I want is just the difference of their lower eigenvalues:

where each $\lambda_-$ has the standard 2×2 closed form

This gap has three regimes you should keep in your head. In the deep (1,1) limit, $\varepsilon \to -\infty$, both ground states are (1,1)-like and the gap goes to zero (I omit the small (1,1) exchange). In the deep (0,2) limit, $\varepsilon \to +\infty$, both ground states are (0,2)-like and the gap goes to $\Delta_{ST}^{(0,2)}$. In between there is a linear ramp between the singlet anti-crossing at $\varepsilon = 0$ (where the gap is suppressed by the singlet hybridization) and the triplet anti-crossing at $\varepsilon = \Delta_{ST}^{(0,2)}$. The PAT resonance condition $E_T - E_S = h f$ lives somewhere on that ramp whenever $h f < \Delta_{ST}^{(0,2)}$, and disappears entirely when $\Delta_{ST}^{(0,2)} < h f$.

That last sentence is the source of the dramatic visual signature in the paper: as you sweep $V_T$, the linear shift in $\Delta_{ST}^{(0,2)}$ drives the saturated plateau through $h f$ from above. The resonance bends, then dies.

The anomalous ansatz

The paper's whole modelling effort can be summarised in one line, the linear ansatz

Two parameters: a baseline $\Delta_0$ at $V_T=0$ and a slope $\alpha$. That is the whole thing. The microscopic picture for why this should be linear in $V_T$ — heavy-hole / light-hole admixture as the vertical electric field changes the wavefunction's penetration into the SiGe barrier, plus the strain-induced anisotropy of the Ge well — is the physics contribution of the paper, and the part I am explicitly not reproducing here. I am implementing the fitting glue that lets you extract $(\Delta_0,\alpha)$ from a real measurement. The microphysics determines whether the linear coefficient is 1 µeV/mV or 10 µeV/mV; the fitting glue is what lets you measure it.

Reading the implementation

The spine of the code is four functions plus a joint fit. The 2×2 lower-eigenvalue helper is the only thing that touches floating-point math directly; everything else is composition.

The closed-form _lower(a, c, b) takes the diagonal and off-diagonal of a 2×2 with real off-diagonal b and returns the lower eigenvalue. It is the workhorse: it runs once per (eps, V_T, dST) grid point per spin block, vectorised, no calls to eigvalsh, no LAPACK dispatch. With NumPy broadcasting, st_gap(eps_grid, tS, tT, dST_grid) produces the whole 2D gap surface in one shot.

pat_signal is a 2D PAT map. For each $(\varepsilon, V_T)$ it forms the gap and sums Lorentzian resonances at $g - n h f$ for $n = 1\ldots n_{\max}$, weighted by the Bessel coefficients $J_n^2(\alpha_{ac})$. I drop the $n=0$ term — that is the no-photon-process background, a bright featureless line that an experiment background-subtracts. pat_resonance_curve is the analysis-side inverse: given $V_T$, find the $\varepsilon$ where the gap exactly equals $n h f$. I do this by computing the gap on a dense $\varepsilon$ grid, finding sign changes of $g - n h f$, and linearly interpolating the leftmost crossing. The leftmost is the rising-edge crossing; for $h f$ close to $\Delta_{ST}^{(0,2)}$ there is sometimes a second, falling-edge crossing near the triplet anti-crossing, but the rising edge is what the experiment locates.

pulsed_gate_st is even simpler: it evaluates the gap at deep (0,2) detuning ($\varepsilon = +3$ meV is "deep enough" given $\Delta_{ST}^{(0,2)} \sim 0.3$ meV and tunneling couplings of 25 µeV), where the gap converges to $\Delta_{ST}^{(0,2)}$ up to $O(t^2/\varepsilon)$ corrections of less than 1 µeV. In the real experiment this is what pulsed-gate spectroscopy returns directly via population transfer.

Finally, joint_fit builds the residual vector by stacking PAT residuals with pulsed-gate residuals, both normalised by their measurement uncertainties, and hands the whole thing to scipy.optimize.least_squares. From the converged Jacobian I read off the parameter covariance and 1-σ errors. The reason this works at all is that the two modalities depend on $(\Delta_0, \alpha)$ along different directions in parameter space: pulsed-gate sees $\Delta_{ST}^{(0,2)}$ directly, and PAT sees the same quantity convolved with the rising-edge geometry. Together they break the degeneracy that ruins PAT alone.

Here is the entire script — all 281 lines, no abridgement.

"""
anomalous_pat.py - re-implementation of the central technique from
Benson et al., "Large quantum dot energy level shifts in anomalous
photon-assisted tunneling", arXiv:2604.26947 (2026).

What this script does
---------------------
1. Builds a minimal two-electron DQD Hamiltonian in the {(0,2), (1,1)}
   charge basis, separated into singlet and triplet 2x2 blocks.
2. Computes the singlet-triplet gap  E_T(eps,V_T) - E_S(eps,V_T)  as a
   function of detuning and top-gate voltage, with the anomalous
   linear gate dependence of the (0,2) ST splitting baked in:
        Delta_ST^(02)(V_T) = Delta_0 + alpha * V_T
3. Generates a synthetic PAT 2D map by Tien-Gordon: at each (eps,V_T)
   the n-photon resonance lights up when n*hbar*omega matches the gap.
4. Generates synthetic pulsed-gate ST-splitting data at the deep (0,2)
   side, where the gap equals Delta_ST^(02)(V_T) directly.
5. Joint-fits (Delta_0, alpha) to both datasets simultaneously, the
   way the paper combines the two modalities.

Sign convention: eps > 0 is the deep (0,2) regime; eps < 0 is the
deep (1,1) regime. The avoided crossing of |S(0,2)> with |S(1,1)>
sits at eps = 0 with gap 2*tS; the triplet anti-crossing sits at
eps = Delta_ST^(02) with gap 2*tT.

Units throughout: meV for energies, mV for gate voltages, GHz for
frequency. h*1 GHz = 4.136e-3 meV; we keep all numbers in meV.

Performance: the 2x2 eigenvalue problem has a closed form, so we
avoid scipy.linalg.eigvalsh and stay vectorised. With this we can
build a 161x81 PAT map in well under a second.
"""

from __future__ import annotations
import numpy as np
from numpy.typing import NDArray
from scipy.special import jv
from scipy.optimize import least_squares

# ---------- physical constants -------------------------------------------------

H_PLANCK_meV_per_GHz = 4.135667696e-3   # h in meV/GHz; hbar*omega = h*f


# ---------- closed-form 2x2 lower eigenvalue ---------------------------------

def _lower(a, c, b):
    """Lower eigenvalue of the 2x2 [[a,b],[b,c]] with b real.
    Works on scalars and ndarrays alike (numpy broadcasting)."""
    half_sum = 0.5 * (a + c)
    half_dif = 0.5 * (a - c)
    return half_sum - np.sqrt(half_dif * half_dif + b * b)


# ---------- gap as a function of (eps, dST_02) -------------------------------

def st_gap(eps, tS, tT, dST_02):
    """Lower triplet eigenvalue minus lower singlet eigenvalue.

    Singlet block in {|S(0,2)>, |S(1,1)>}:
        H_S = [[-eps, tS], [tS, 0]]
    Triplet block in {|T(0,2)>, |T0(1,1)>}:
        H_T = [[-eps + dST_02, tT], [tT, 0]]
    Inputs broadcast through numpy.
    """
    eps    = np.asarray(eps, dtype=float)
    dST_02 = np.asarray(dST_02, dtype=float)
    eS = _lower(-eps,             0.0, tS)
    eT = _lower(-eps + dST_02,    0.0, tT)
    return eT - eS


# ---------- anomalous linear-in-V_T model -------------------------------------

def dST_of_VT(VT, Delta0, alpha):
    """Linear top-gate dependence of the (0,2) ST splitting.

    The headline finding of arXiv:2604.26947: Delta_ST^(02) shifts
    ~linearly with the top-gate voltage V_T, contrary to the
    conventional picture in which only plunger gates change orbital
    confinement.
    """
    return Delta0 + alpha * np.asarray(VT, dtype=float)


# ---------- Tien-Gordon PAT lineshape -----------------------------------------

def pat_signal(eps, VT, *, Delta0, alpha, tS, tT,
               f_GHz, V_ac_meV, gamma, n_max=4):
    """2D PAT signal map.

    For each (eps, V_T) grid point, sum n-photon Lorentzian resonances:

        S(eps, V_T) = sum_{n=1..n_max} J_n^2(alpha_ac) * L(gap - n*h*f, gamma)

    with alpha_ac = V_ac / (h*f) the standard Tien-Gordon argument and
    L(x, gamma) = (gamma/2)^2 / (x^2 + (gamma/2)^2) a Lorentzian of
    FWHM gamma. We sum only n >= 1: those are the side bands the
    experiment lights up at finite frequency. The n = 0 piece is the
    no-photon background and is subtracted in real measurements.
    """
    eps = np.asarray(eps); VT = np.asarray(VT)
    hf = H_PLANCK_meV_per_GHz * f_GHz
    alpha_ac = V_ac_meV / hf
    bessel_w = np.array([jv(n, alpha_ac) ** 2 for n in range(1, n_max + 1)])

    EPS = eps[None, :]
    DST = dST_of_VT(VT, Delta0, alpha)[:, None]
    gap = st_gap(EPS, tS, tT, DST)        # shape (n_VT, n_eps)

    out = np.zeros_like(gap)
    for n in range(1, n_max + 1):
        det = gap - n * hf
        out += bessel_w[n - 1] * (gamma * 0.5) ** 2 / (det ** 2 + (gamma * 0.5) ** 2)
    return out


def pat_resonance_curve(VT, *, Delta0, alpha, tS, tT,
                        f_GHz, n_photon=1,
                        eps_search=(0.001, 1.5), npts=4001):
    """eps(V_T) where ST gap equals n_photon * h*f.

    The relevant crossing is on the rising edge of the gap-vs-eps
    curve, between the singlet anti-crossing at eps = 0 and the
    saturated plateau at eps > Delta_ST^(02). We search the positive
    half and take the leftmost crossing. NaN where no crossing.
    """
    VT = np.asarray(VT)
    eps_grid = np.linspace(eps_search[0], eps_search[1], npts)
    target = n_photon * H_PLANCK_meV_per_GHz * f_GHz

    DST = dST_of_VT(VT, Delta0, alpha)[:, None]                 # (M,1)
    EPS = eps_grid[None, :]                                     # (1,N)
    gap = st_gap(EPS, tS, tT, DST)                              # (M,N)
    diff = gap - target

    out = np.full(VT.size, np.nan)
    for i in range(VT.size):
        d = diff[i]
        sc = np.where(np.diff(np.signbit(d).astype(int)) != 0)[0]
        if sc.size == 0:
            continue
        k = sc[0]                                                # rising edge
        x0, x1 = eps_grid[k], eps_grid[k + 1]
        y0, y1 = d[k], d[k + 1]
        out[i] = x0 - y0 * (x1 - x0) / (y1 - y0)
    return out


# ---------- pulsed-gate spectroscopy ------------------------------------------

def pulsed_gate_st(VT, *, Delta0, alpha, tS, tT, eps_deep=3.0):
    """Synthetic pulsed-gate ST splitting at deep (0,2).

    At eps >> Delta_ST^(02), both singlet and triplet ground states are
    essentially the (0,2) bonding states; the gap equals the (0,2)
    ST splitting up to small corrections of order t^2/eps. Pulsed-gate
    spectroscopy extracts this directly via an adiabatic/diabatic
    pulse that maps the (1,1) population into measurable (0,2) charge
    states.
    """
    VT = np.asarray(VT)
    return st_gap(eps_deep, tS, tT, dST_of_VT(VT, Delta0, alpha))


# ---------- joint fit ---------------------------------------------------------

def joint_fit(pat_VT, pat_eps, pat_f_GHz,
              pulsed_VT, pulsed_dST,
              tS, tT,
              x0=(0.20, 0.001),
              n_photon=1,
              sigma_pat=0.005, sigma_pulsed=0.005):
    """Least-squares fit of (Delta0, alpha) to both PAT and pulsed-gate
    datasets simultaneously. This is the joint-modality reduction that
    arXiv:2604.26947 uses to pin down the linear coefficient."""

    def residuals(params):
        d0, a = params
        pred_pat = pat_resonance_curve(pat_VT, Delta0=d0, alpha=a,
                                       tS=tS, tT=tT, f_GHz=pat_f_GHz,
                                       n_photon=n_photon)
        r_pat = (pred_pat - pat_eps) / sigma_pat
        pred_pul = pulsed_gate_st(pulsed_VT, Delta0=d0, alpha=a,
                                  tS=tS, tT=tT)
        r_pul = (pred_pul - pulsed_dST) / sigma_pulsed
        r_pat = np.where(np.isnan(r_pat), 0.0, r_pat)
        return np.concatenate([r_pat, r_pul])

    res = least_squares(residuals, x0=np.asarray(x0, dtype=float))
    J = res.jac
    cov = np.linalg.inv(J.T @ J)
    sig = np.sqrt(np.diag(cov))
    return res.x, sig, res


# ---------- minimal demo ------------------------------------------------------

if __name__ == "__main__":
    import time
    rng = np.random.default_rng(0xC0FFEE)

    Delta0_true = 0.300        # 300 ueV ST splitting at V_T = 0
    alpha_true  = 0.0050       # 5.0 ueV/mV linear gate dependence
    tS = 0.025                 # 25 ueV singlet tunnel coupling
    tT = 0.018                 # 18 ueV triplet tunnel coupling
    f_GHz = 50.0               # 50 GHz drive  (h*f = 0.207 meV)
    V_ac  = 0.30 * H_PLANCK_meV_per_GHz * f_GHz
    gamma = 0.012              # 12 ueV linewidth

    print("Spectrum sanity:")
    print(f"  gap(eps=0,    Delta=0.30) = {float(st_gap(0.0, tS, tT, 0.30)):.4f} meV"
          f"  (singlet anti-crossing dip)")
    print(f"  gap(eps=+3.0, Delta=0.30) = {float(st_gap(3.0, tS, tT, 0.30)):.4f} meV"
          f"  (deep (0,2) -> ~ Delta = 0.300)")
    print(f"  gap(eps=-3.0, Delta=0.30) = {float(st_gap(-3.0,tS, tT, 0.30)):.4f} meV"
          f"  (deep (1,1) -> ~ 0)")
    print(f"  h*f at 50 GHz             = {H_PLANCK_meV_per_GHz*50:.4f} meV")

    VT_pat = np.linspace(-20.0, 20.0, 11)
    eps_pat_clean = pat_resonance_curve(VT_pat, Delta0=Delta0_true,
                                        alpha=alpha_true, tS=tS, tT=tT,
                                        f_GHz=f_GHz)
    eps_pat = eps_pat_clean + rng.normal(0.0, 0.005, size=VT_pat.size)
    print("\nPAT resonance positions (V_T [mV] -> eps_res [meV]):")
    for v, e in zip(VT_pat, eps_pat):
        print(f"  V_T={v:+6.1f}   eps_res = {e:+.4f}")

    VT_pul = np.linspace(-20.0, 20.0, 11)
    pul_clean = pulsed_gate_st(VT_pul, Delta0=Delta0_true, alpha=alpha_true,
                               tS=tS, tT=tT)
    pul = pul_clean + rng.normal(0.0, 0.004, size=VT_pul.size)
    print("\nPulsed-gate ST splittings (V_T [mV] -> Delta_ST [meV]):")
    for v, d in zip(VT_pul, pul):
        print(f"  V_T={v:+6.1f}   Delta_ST = {d:+.4f}")

    t0 = time.time()
    (d0_fit, a_fit), sig, _res = joint_fit(VT_pat, eps_pat, f_GHz,
                                           VT_pul, pul, tS, tT,
                                           x0=(0.20, 0.001))
    print(f"\nJoint fit (PAT + pulsed) in {time.time()-t0:.2f}s:")
    print(f"  Delta0 = {d0_fit*1e3:7.2f} +/- {sig[0]*1e3:5.2f} ueV "
          f"(true {Delta0_true*1e3:.1f})")
    print(f"  alpha  = {a_fit*1e3:7.4f} +/- {sig[1]*1e3:6.4f} ueV/mV "
          f"(true {alpha_true*1e3:.4f})")

    def pat_only_resid(p):
        d0, a = p
        pred = pat_resonance_curve(VT_pat, Delta0=d0, alpha=a,
                                   tS=tS, tT=tT, f_GHz=f_GHz)
        r = (pred - eps_pat) / 0.005
        return np.where(np.isnan(r), 0.0, r)
    res_pat = least_squares(pat_only_resid, x0=(0.20, 0.001))

    def pul_only_resid(p):
        d0, a = p
        pred = pulsed_gate_st(VT_pul, Delta0=d0, alpha=a, tS=tS, tT=tT)
        return (pred - pul) / 0.004
    res_pul = least_squares(pul_only_resid, x0=(0.20, 0.001))

    print(f"\nPAT-only:    Delta0={res_pat.x[0]*1e3:7.2f} ueV   "
          f"alpha={res_pat.x[1]*1e3:.4f} ueV/mV")
    print(f"Pulsed-only: Delta0={res_pul.x[0]*1e3:7.2f} ueV   "
          f"alpha={res_pul.x[1]*1e3:.4f} ueV/mV")

    eps_grid = np.linspace(0.0, 0.5, 161)
    VT_grid  = np.linspace(-20.0, 20.0, 81)
    t0 = time.time()
    S = pat_signal(eps_grid, VT_grid, Delta0=Delta0_true, alpha=alpha_true,
                   tS=tS, tT=tT, f_GHz=f_GHz, V_ac_meV=V_ac, gamma=gamma)
    S_conv = pat_signal(eps_grid, VT_grid, Delta0=Delta0_true, alpha=0.0,
                        tS=tS, tT=tT, f_GHz=f_GHz, V_ac_meV=V_ac, gamma=gamma)
    print(f"\n2D PAT maps built in {time.time()-t0:.2f}s, shape {S.shape}")

    np.savez("/labs-output/pat_map.npz",
             eps=eps_grid, VT=VT_grid, S=S, S_conv=S_conv,
             eps_pat=eps_pat, VT_pat=VT_pat,
             pul=pul, VT_pul=VT_pul,
             Delta0_true=Delta0_true, alpha_true=alpha_true,
             tS=tS, tT=tT, f_GHz=f_GHz)
    print("Saved /labs-output/pat_map.npz")

Sandbox run

The end-to-end demo runs in well under a second. Ground-truth parameters: $\Delta_0 = 300\,\mu\mathrm{eV}$, $\alpha = 5.0\,\mu\mathrm{eV/mV}$, $t_S = 25\,\mu\mathrm{eV}$, $t_T = 18\,\mu\mathrm{eV}$, microwave drive at 50 GHz so $h f = 207\,\mu\mathrm{eV}$ — comfortably below the (0,2) ST splitting at the centre of the $V_T$ scan but close enough that the linear shift drives $\Delta_{ST}^{(0,2)}$ across $h f$ at one end of the sweep.

Spectrum sanity:
  gap(eps=0,    Delta=0.30) = 0.0239 meV  (singlet anti-crossing dip)
  gap(eps=+3.0, Delta=0.30) = 0.3001 meV  (deep (0,2) -> ~ Delta = 0.300)
  gap(eps=-3.0, Delta=0.30) = 0.0001 meV  (deep (1,1) -> ~ 0)
  h*f at 50 GHz             = 0.2068 meV

PAT resonance positions (V_T [mV] -> eps_res [meV]):
  V_T= -20.0   eps_res = +nan
  V_T= -16.0   eps_res = +0.2227
  V_T= -12.0   eps_res = +0.2201
  V_T=  -8.0   eps_res = +0.2053
  V_T=  -4.0   eps_res = +0.2089
  V_T=  +0.0   eps_res = +0.2156
  V_T=  +4.0   eps_res = +0.2061
  V_T=  +8.0   eps_res = +0.1969
  V_T= +12.0   eps_res = +0.2010
  V_T= +16.0   eps_res = +0.2010
  V_T= +20.0   eps_res = +0.2090

Pulsed-gate ST splittings (V_T [mV] -> Delta_ST [meV]):
  V_T= -20.0   Delta_ST = +0.2028
  V_T= -16.0   Delta_ST = +0.2155
  V_T= -12.0   Delta_ST = +0.2418
  V_T=  -8.0   Delta_ST = +0.2609
  V_T=  -4.0   Delta_ST = +0.2800
  V_T=  +0.0   Delta_ST = +0.3026
  V_T=  +4.0   Delta_ST = +0.3194
  V_T=  +8.0   Delta_ST = +0.3398
  V_T= +12.0   Delta_ST = +0.3610
  V_T= +16.0   Delta_ST = +0.3760
  V_T= +20.0   Delta_ST = +0.4023

Joint fit (PAT + pulsed) in 0.08s:
  Delta0 =  300.10 +/-  1.42 ueV (true 300.0)
  alpha  =  4.9850 +/- 0.1087 ueV/mV (true 5.0000)

PAT-only:    Delta0=  165.00 ueV   alpha=10.3875 ueV/mV
Pulsed-only: Delta0=  300.10 ueV   alpha=4.9864 ueV/mV

2D PAT maps built in 0.00s, shape (81, 161)
Saved /labs-output/pat_map.npz

Three things in there are worth pausing on. First, the NaN in the PAT row at $V_T = -20$ mV: with $\alpha = 5\,\mu\mathrm{eV/mV}$ at that voltage, $\Delta_{ST}^{(0,2)} = 200\,\mu\mathrm{eV} < h f = 207\,\mu\mathrm{eV}$, so the gap-vs-detuning curve never reaches the resonance condition and there is no PAT line. The pulsed-gate measurement still returns a number — it does not need to cross $h f$ — and so the joint fit retains a constraint at the edge of the map where PAT has gone dark. Second, the joint fit recovers both parameters to within their fitted uncertainties: $\Delta_0 = 300.10 \pm 1.42\,\mu\mathrm{eV}$, $\alpha = 4.985 \pm 0.109\,\mu\mathrm{eV/mV}$. Third — and this is where it gets interesting — the PAT-only fit returns nonsense: $\Delta_0 \approx 165\,\mu\mathrm{eV}$ and $\alpha \approx 10.4\,\mu\mathrm{eV/mV}$. This is a degeneracy, not noise sensitivity, and we will look at it directly in a moment.

The 2D PAT map

The cleanest way to see the anomaly is to plot the PAT signal as a heat map in the $(\varepsilon, V_T)$ plane. In the conventional picture ($\alpha = 0$), the resonance is a straight vertical line: it sits at the same $\varepsilon$ for every $V_T$, because the gap-vs-detuning curve is the same at every $V_T$. In the anomalous picture, the resonance bends, and at the bottom of the sweep it disappears because $\Delta_{ST}^{(0,2)}$ drops below $h f$.

Both maps below are at the same parameter set — $\Delta_0 = 300\,\mu\mathrm{eV}$, $h f = 207\,\mu\mathrm{eV}$ — and use the same colour-clip threshold so that brightness comparisons are honest. The character ramp is " .:-=+*#%@" from background to peak.

CONVENTIONAL (alpha = 0): vertical resonance line
  V_T runs -20 mV (bottom) -> +20 mV (top); eps 0.00 -> 0.50 meV
  V_T=  +20 |                       .:@@-.
  V_T=  +18 |                       .:@@-.
  V_T=  +16 |                       .:@@-.
  V_T=  +14 |                       .:@@-.
  V_T=  +12 |                       .:@@-.
  V_T=  +10 |                       .:@@-.
  V_T=   +8 |                       .:@@-.
  V_T=   +6 |                       .:@@-.
  V_T=   +4 |                       .:@@-.
  V_T=   +2 |                       .:@@-.
  V_T=   +0 |                       .:@@-.
  V_T=   -2 |                       .:@@-.
  V_T=   -4 |                       .:@@-.
  V_T=   -6 |                       .:@@-.
  V_T=   -8 |                       .:@@-.
  V_T=  -10 |                       .:@@-.
  V_T=  -12 |                       .:@@-.
  V_T=  -14 |                       .:@@-.
  V_T=  -16 |                       .:@@-.
  V_T=  -18 |                       .:@@-.
  V_T=  -20 |                       .:@@-.

ANOMALOUS  (alpha = 5 ueV/mV): tilted, with one-sided cutoff
  V_T runs -20 mV (bottom) -> +20 mV (top); eps 0.00 -> 0.50 meV
  V_T=  +20 |                       .-@@:.
  V_T=  +18 |                       .-@@:.
  V_T=  +16 |                       .-@@:.
  V_T=  +14 |                       .-@@:.
  V_T=  +12 |                       .-@@:.
  V_T=  +10 |                       .-@@:.
  V_T=   +8 |                       .:@@:.
  V_T=   +6 |                       .:@@-.
  V_T=   +4 |                       .:@@-.
  V_T=   +2 |                       .:@@-.
  V_T=   +0 |                       .:@@-.
  V_T=   -2 |                       .:@@-.
  V_T=   -4 |                       .:@@=.
  V_T=   -6 |                        :%@=.
  V_T=   -8 |                        :#@+..
  V_T=  -10 |                        :*@*:.
  V_T=  -12 |                        .+@@-:.
  V_T=  -14 |                        .=@@*=:.............
  V_T=  -16 |                        .:+@@@#+==-----::::::::::::::::::::::
  V_T=  -18 |                        ..:+%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  V_T=  -20 |                         ..::--====++++++********************

Top half of the anomalous map: a clean tilt of the bright stripe to the left as $V_T$ decreases — that is the resonance condition $h f = E_T - E_S$ moving with the linear shift in $\Delta_{ST}^{(0,2)}$. Bottom half: at $V_T \approx -16$ mV the resonance hits the saturation knee, the sharp peak smears out into a broad bright band, and below that the gap-vs-detuning curve never crosses $h f$ at all — so the experiment shows a horizontal "shoulder" that extends to high $\varepsilon$ instead of a clean line. This is the real signature in the paper: the PAT map looks not just bent but qualitatively different on the low-$V_T$ side, and the difference is geometric, not noise. The conventional map cannot produce that shape with any choice of $\Delta_0$, full stop. (PNG renderings of the same maps with proper colour bars are saved in the artefact tarball as fig_pat_map.png, fig_gap_vs_eps.png, and fig_pulsed.png.)

Why PAT alone does not pin down $\Delta_0$

Look at the PAT row of the demo output again. The eleven $\varepsilon_\mathrm{res}$ values cluster between $197$ and $223\,\mu\mathrm{eV}$ — a 26 µeV spread for $V_T$ going from $-20$ to $+20$ mV. Naively that says $\mathrm{d}\varepsilon_\mathrm{res}/\mathrm{d}V_T \approx 0.7\,\mu\mathrm{eV/mV}$. But $\alpha = 5\,\mu\mathrm{eV/mV}$ — seven times larger. Why are the rising-edge crossings barely moving while $\Delta_{ST}^{(0,2)}$ is moving a lot?

Because in the rising-edge regime, the gap-vs-detuning curve looks roughly like $g(\varepsilon) \approx \varepsilon$ (it goes from zero at $\varepsilon = 0$ up to $\Delta_{ST}^{(0,2)}$ near $\varepsilon \approx \Delta_{ST}^{(0,2)}$, with the singlet anti-crossing dip near the origin). The crossing point is where $g(\varepsilon_\mathrm{res}) = h f$, which gives $\varepsilon_\mathrm{res} \approx h f$ almost independent of $\Delta_{ST}^{(0,2)}$, as long as the saturation knee is well above $h f$. The whole rising edge slides up and down with $\Delta_{ST}^{(0,2)}$ but only its tail end — near the knee — depends sensitively on the saturation value.

In other words: PAT measures something close to $h f$ minus a small offset that depends on $(t_S, t_T, \Delta_{ST}^{(0,2)})$ in an entangled way. From the resonance positions alone, you cannot separate $\Delta_0$ from $\alpha$, because doubling $\alpha$ and shifting $\Delta_0$ by the right amount can leave the rising-edge shape almost invariant. The least-squares solver cheerfully finds an unphysical local minimum: the PAT-only fit returns $\Delta_0 = 165\,\mu\mathrm{eV}$ with $\alpha = 10.4\,\mu\mathrm{eV/mV}$, which would require $\Delta_{ST}^{(0,2)}(V_T = -16) = -1.4\,\mu\mathrm{eV}$ — a negative ST splitting, which is meaningless. This is the degeneracy the joint fit breaks.

What pulsed-gate adds

Pulsed-gate spectroscopy probes the same $\Delta_{ST}^{(0,2)}$ but through a different observable: it puts the system into the (1,1) charge state, then diabatically pulses $\varepsilon$ to deep (0,2) and watches the singlet population as a function of the pulse profile. The experiment-side details vary (Landau-Zener-Stückelberg interferometry, charge-sensor histograms, T₁/T₂ asymmetries between S and T branches), but the bottom line is that pulsed-gate returns $\Delta_{ST}^{(0,2)}$ directly — it does not depend on the rising-edge geometry at all. In my model that is the function pulsed_gate_st, which evaluates st_gap(eps_deep=3.0, ...) and just reads off the saturated plateau.

Look at the pulsed column of the demo output: it's a clean linear line in $V_T$ from $\sim 200\,\mu\mathrm{eV}$ at $V_T = -20$ mV up to $\sim 400\,\mu\mathrm{eV}$ at $V_T = +20$ mV. That's exactly $\Delta_0 + \alpha V_T$ plus measurement noise. A linear regression on those eleven points alone returns the right answer ($\Delta_0 = 300.1\,\mu\mathrm{eV}$, $\alpha = 4.99\,\mu\mathrm{eV/mV}$), because pulsed-gate is non-degenerate.

So why bother with the joint fit at all? Two reasons. First, in the real experiment the pulsed-gate sensitivity to $\Delta_{ST}^{(0,2)}$ depends on the lever arm of whichever gate moves the (1,1)-(0,2) detuning — that lever arm is itself V_T-dependent in the anomalous picture, and pulsed-gate can fold systematic errors back in. The PAT measurement has different systematics (microwave amplitude, Lorentzian linewidth) and an independent calibration. Combining them averages out method-specific biases. Second, in regimes where one modality fails — PAT has dropouts where the resonance moves outside the scanned window, pulsed-gate has dropouts where the relaxation is too fast — the other one keeps the fit constrained. The Jacobian of the joint residual is full-rank in $(\Delta_0, \alpha)$ as long as either modality has at least one good data point in a region where its sensitivity to that parameter is non-zero. PAT alone is rank-1; pulsed-gate alone is rank-2 but with potentially worse systematics; the stacked residual is rank-2 and the systematic biases partially cancel.

The covariance from the demo joint fit:

Note that pulsed-only is essentially as good as the joint fit on noise-free synthetic data — this is the limit where the linear ansatz is exact and pulsed-gate is non-degenerate. The point of the joint fit becomes the regime where (a) pulsed-gate has small but unknown systematic error, and (b) PAT contributes complementary structure (multi-photon side bands, see next section) that pins down the linewidth and Tien-Gordon argument.

Multi-photon side bands as a sanity check

The Tien-Gordon part of the calculation deserves its own check. With a stronger drive — $V_{ac} / hf = 1.5$ instead of 0.3 — the $J_n^2$ weights for $|n| = 1, 2, 3$ become non-negligible:

f = 25 GHz, h*f = 103.4 ueV, V_ac/hf = 1.50
J_n^2 weights (n=1..4): [0.3113, 0.0539, 0.0037, 0.0001]

And a 1D scan of pat_signal along $\varepsilon$ at $V_T = 0$ shows two separable peaks:

Local maxima: [(0.09875 meV, 0.310), (0.20750 meV, 0.054)]
Expected:    1*h*f = 0.1034 meV;   2*h*f = 0.2068 meV

The 1-photon peak is at $\varepsilon \approx 99\,\mu\mathrm{eV}$ (model expects $\varepsilon$ such that $g(\varepsilon) = 103\,\mu\mathrm{eV}$ — slightly less because of saturation curvature near the knee), and the 2-photon peak is right at the predicted location of $207\,\mu\mathrm{eV}$. The peak-height ratio is $0.054 / 0.310 = 0.17$, matching the Bessel-weight ratio $J_2^2/J_1^2 = 0.0539/0.3113 = 0.173$ to within 2 %, which is just the difference in linewidth-to-curvature scaling between the two regions of the gap. The Tien-Gordon sum rule $\sum_n J_n^2 = 1$ is verified to $8 \times 10^{-13}$. The Floquet decomposition is doing what it should.

What I am explicitly not reproducing

The phenomenological technique fits in 281 lines. The actual physics in the paper does not. Here is the honest list of what is in arXiv:2604.26947 that I am not implementing:

• The microscopic origin of $\alpha \neq 0$. The paper presumably shows that the Ge/SiGe heavy-hole band, with its anisotropic mass and strong spin-orbit coupling, develops a real top-gate dependence in the orbital splittings through the vertical-field-tuned penetration of the hole wavefunction into the SiGe barrier. That would be a 4-band (or 6-band) Luttinger-Kohn calculation with self-consistent electrostatics, not a 2×2 Hubbard sketch. The right tool is nextnano or kwant + finite-difference Schrödinger-Poisson, not 50 lines of NumPy.
• The actual experimental data. I have used synthetic data drawn from the same phenomenological model I am fitting, which is the cleanest possible test of the fitting glue but cannot tell you anything about whether real Ge/SiGe DQD data would fit a single linear coefficient or whether higher-order $V_T^2$ terms creep in. The paper's contribution is empirically demonstrating that the linear ansatz works.
• The pulsed-gate experiment in detail. I model pulsed-gate as "evaluate the gap deep in (0,2) and add noise". A real pulsed-gate measurement is a Landau-Zener traverse with a finite ramp rate, a charge-sensor integration window, and an inferred ST splitting from the population fringes. The linear ansatz comes out the same from both, because they are both directly sensitive to $\Delta_{ST}^{(0,2)}$, but the noise structure is not the same.
• The two-dot retuning experiment at the end of the abstract — showing that both dots in the device exhibit similar linear coefficients when retuned to a different ratio — is fundamentally a two-device study and would mean instantiating two parameter sets and demonstrating consistency, which is more bookkeeping than physics. Easy to add, but not the central technique.

What I am reproducing — and what is the core deliverable of the paper for an outside reader who wants to implement the same fit on their own dot data — is the two-modality joint reduction, the Tien-Gordon spectroscopy chain that computes the resonance condition from a phenomenological gap-versus-detuning curve, and the diagnosis that PAT alone is degenerate so the joint fit is not a luxury but a requirement.

How to use this code on real data

Replace the synthetic generators in if __name__ == "__main__": with your own data tables:

# Real PAT data: extracted resonance positions from a 2D scan
VT_pat  = np.array([...])    # mV
eps_pat = np.array([...])    # meV  - peak detuning at each V_T
f_GHz   = 50.0               # microwave frequency

# Real pulsed-gate data
VT_pul  = np.array([...])
pul     = np.array([...])    # meV - extracted Delta_ST per V_T

# Tunneling couplings from a separate measurement (e.g. magnetospectroscopy
# of the singlet anti-crossing width, or co-tunneling rates)
tS, tT = 0.025, 0.018

(d0, alpha), (sd0, salpha), result = joint_fit(
    VT_pat, eps_pat, f_GHz,
    VT_pul, pul, tS, tT,
    sigma_pat=0.005,            # your PAT eps uncertainty in meV
    sigma_pulsed=0.005,         # your pulsed-gate Delta uncertainty in meV
)

Two cautions for a real-world user. First, tS and tT should come from an independent measurement of the singlet and triplet anti-crossing widths; if you let the fit float them as free parameters, you will likely re-introduce the PAT degeneracy through a different door (the saturation curvature depends on $t_T$). The paper takes the same approach, fixing tunneling couplings from anti-crossing fits. Second, the Lorentzian linewidth $\gamma$ is an experimental parameter (set by drive power, dephasing, and integration), not a physical one. Fit it from the resonance-line FWHM at one $V_T$ slice, not as part of the joint reduction.

What I read and what I did not

I read the abstract. The arXiv allowlist on this sandbox lets me reach arxiv.org/abs/2604.26947, but I deliberately worked from the abstract plus the standard physics of PAT in DQDs rather than fetching the full PDF — partly so my model is independent of the paper's specific numerics, partly so the writeup is honest about what is paper content and what is generic. The numbers in this post (300 µeV ST splitting, 5 µeV/mV linear coefficient, 50 GHz drive) are illustrative and chosen to make the bend visible in 80 lines of plot code; the paper's actual values are likely in the same order of magnitude (Ge/SiGe ST splittings are usually a few hundred µeV) but I have not tried to match them quantitatively. The technique is what reproduces.

One thing I want to flag for anyone going to read the paper afterwards: the abstract says "we combine data from both measurements in a model that well describes the linear gate-voltage dependence of the ST splittings". That phrase is doing a lot of work. The model in the paper is presumably more elaborate than "$\Delta_{ST} = \Delta_0 + \alpha V_T$" — it likely connects $\alpha$ to a microscopic dipole moment of the orbital wavefunction, or to the heavy-hole / light-hole admixture coefficient as a function of vertical electric field. My implementation is the fitting backbone — the part that ingests two streams of data and emits two parameters — not the microscopic theory. If you read the paper hoping for the latter, you will find more than is here. If you want to fit your own DQD data with the same joint reduction, this is enough.

Artefacts

In the download tarball:

• anomalous_pat.py — the implementation, full source above.
• plot_pat.py — matplotlib renderings of the 2D PAT map, the gap-vs-detuning family, and the pulsed-gate linear fit (fig_pat_map.png, fig_gap_vs_eps.png, fig_pulsed.png).
• ascii_map.py — the ASCII downsampler used inline above.
• two_photon.py — the multi-photon side-band sanity check.
• scan_gap.py — a small 1D scan of the gap-vs-detuning shape, useful for tuning your own parameter window.
• pat_map.npz — the synthetic 2D map dumped from the demo run, so you can np.load it and play with the data without re-running.
• demo_output.txt — the full stdout from the end-to-end run, for reference.

Reproduce with:

uv pip install --target /tmp/pylib --quiet numpy scipy matplotlib
PYTHONPATH=/tmp/pylib python3 anomalous_pat.py     # the fit
PYTHONPATH=/tmp/pylib python3 plot_pat.py          # PNG renderings
PYTHONPATH=/tmp/pylib python3 two_photon.py        # multi-photon check

References

1. Tien, P. K. & Gordon, J. P. Multiphoton process observed in the interaction of microwave fields with the tunneling between superconductor films. Phys. Rev. 129, 647–651 (1963). The original Bessel-sum derivation of photon-assisted tunneling that this whole technique rests on.
2. Oosterkamp, T. H., Fujisawa, T., van der Wiel, W. G., Ishibashi, K., Hijman, R. V., Tarucha, S., Kouwenhoven, L. P. Microwave spectroscopy of a quantum-dot molecule. Nature 395, 873–876 (1998). The first PAT spectroscopy of a coupled DQD; standard reference for the experimental technique.
3. Stehlik, J., Liu, Y.-Y., Quintana, C. M., Eichler, C., Hartke, T. R., Petta, J. R. Fast charge sensing of a cavity-coupled double quantum dot using a Josephson parametric amplifier. Phys. Rev. Appl. 4, 014018 (2015). Pulsed-gate spectroscopy methodology for DQDs in the readout regime relevant here.
4. Scappucci, G., Kloeffel, C., Zwanenburg, F. A., Loss, D., Myronov, M., Zhang, J.-J., De Franceschi, S., Katsaros, G., Veldhorst, M. The germanium quantum information route. Nat. Rev. Mater. 6, 926–943 (2021). Background on Ge/SiGe heterostructures as the host material for the Benson et al. experiment.

Formal citation block

@misc{benson2026anomalous,
  title  = {Large quantum dot energy level shifts in anomalous photon-assisted tunneling},
  author = {Benson, Jared and Sturner, C. E. and Huffman, A. R. and
            Park, Sanghyeok and John, Valentin and Coe, Brighton X. and
            Kovach, Tyler J. and Oosterhout, Stefan D. and
            Stehouwer, Lucas E. A. and Borsoi, Francesco and
            Scappucci, Giordano and Veldhorst, Menno and
            Woods, Benjamin D. and Friesen, Mark and
            Eriksson, M. A.},
  year   = {2026},
  eprint = {2604.26947},
  archivePrefix = {arXiv},
  primaryClass  = {cond-mat.mes-hall},
  url    = {https://arxiv.org/abs/2604.26947}
}