Zero-Noise Extrapolation

Zero-noise extrapolation is the most disarmingly simple idea in the near-term quantum toolbox: you cannot switch the noise off, but you can reliably turn it up, and a curve you can follow uphill you can also run back downhill to where it began. Every gate on a real superconducting or trapped-ion device is a slightly leaky CPTP map, so the expectation value ⟨O⟩\langle O\rangle your circuit reports is systematically biased away from the ideal answer. ZNE treats that biased number not as a fixed defect but as one sample of a smooth function E(λ)E(\lambda) of a tunable noise strength λ\lambda, where λ=1\lambda=1 is the machine's native noise and λ=0\lambda=0 is the perfect circuit you can never actually run. Measure EE at several amplified settings λ=1,3,5,…\lambda=1,3,5,\dots, fit a model, and extrapolate the fit to λ=0\lambda=0. This article is the deep-tier synthesis of the error-mitigation thread. It takes the operator machinery from aq-qci-cptp-maps-quantum-operations and aq-qci-kraus-representation-q9, the continuous-time dissipator from aq-qci-lindblad-gksl-master-equation, and the single-parameter noise archetype from aq-qci-depolarizing-channel. Three consequences organize everything below. First, the whole method rests on analyticity: if the dominant noise is weak and Markovian, a Dyson expansion gives E(λ)=E0+∑kakλkE(\lambda)=E_0+\sum_k a_k\lambda^k, and extrapolation is just reading off E0E_0. Second, the amplification map is free to choose — analog pulse-stretching or the digital trick of unitary folding, U↦U(U†U)nU\mapsto U(U^\dagger U)^n, which changes nothing logically yet piles up decoherence. Third, ZNE is mitigation, not correction: it lowers the systematic bias of an estimator while inflating its variance, buying accuracy with extra shots and no extra qubits — a NISQ-era bargain, not a road to fault tolerance.

Key takeaways
  • Zero-noise extrapolation (ZNE) estimates a noise-free expectation value E0=⟨O⟩λ=0E_0=\langle O\rangle_{\lambda=0} by deliberately running a circuit at several amplified noise levels λ≥1\lambda\ge1 and extrapolating the measured curve E(λ)E(\lambda) back to λ=0\lambda=0.
  • It works because weak Markovian noise makes E(λ)E(\lambda) analytic: E(λ)=E0+∑k≥1akλkE(\lambda)=E_0+\sum_{k\ge1}a_k\lambda^k, so a polynomial, Richardson, or exponential fit inverts to recover E0E_0.
  • Noise is amplified either by analog pulse-stretching or by digital unitary folding U↦U(U†U)nU\mapsto U(U^\dagger U)^n, which leaves the computed unitary unchanged but multiplies the physical error by λ=2n+1\lambda=2n+1.
  • Order-nn Richardson extrapolation cancels the first nn error terms, leaving bias O(λn+1)O(\lambda^{n+1}), but the Lagrange weights γi\gamma_i grow so the variance ∑iγi2 Var[E(λi)]\sum_i\gamma_i^2\,\mathrm{Var}[E(\lambda_i)] balloons — the universal ZNE trade-off.
  • For a single global depolarizing rate E(λ)=(1−p)λE0E(\lambda)=(1-p)^\lambda E_0 exactly, so exponential extrapolation is bias-free; coherent or non-Markovian noise breaks analyticity and ZNE fails.
  • ZNE needs no extra qubits and no noise model, but it neither protects quantum information nor scales to deep circuits — it postpones, rather than replaces, quantum error correction.
  • Key idea: treat the measured expectation value as a smooth function of an artificially tunable noise strength λ\lambda, sample it at several λ≥1\lambda\ge1, and extrapolate the fitted curve to λ=0\lambda=0 — buying reduced bias at the price of increased variance and no extra qubits.
The core insight

Treat the measured expectation value as a smooth function of an artificially tunable noise strength λ\lambda, sample it at several λ≥1\lambda\ge1, and extrapolate the fitted curve to λ=0\lambda=0 to estimate the noiseless result — buying reduced bias at the price of increased variance and no extra qubits.

Where this lives today

Zero-noise extrapolation graduated from theory to production faster than almost any other quantum-information idea. The first hardware demonstration came in 2019, when Kandala and colleagues at IBM combined pulse-stretched noise amplification with Richardson extrapolation to sharpen the variational-eigensolver energy of small molecules on a superconducting processor, measurably extending its computational reach. Digital folding then made the method hardware-agnostic: the open-source Mitiq library implements global and local folding plus a menu of extrapolators, so any circuit can be mitigated in a few lines of code with no pulse access. The most-cited recent milestone is IBM's 2023 "utility" experiment (Kim and coworkers), which ran a 127-qubit kicked-Ising simulation beyond the reach of brute-force classical simulation and relied on a ZNE variant — probabilistic error amplification followed by extrapolation — to pull trustworthy expectation values out of deeply noisy circuits, reframing the debate about near-term quantum advantage. Today ZNE is a standard layer in the mitigation stack, routinely composed with readout-error mitigation and, where a noise model is available, probabilistic error cancellation. It remains a near-term mitigation strategy rather than a replacement for fault tolerance — but it can extend the useful regime of NISQ-era experiments.

The physical picture

You are handed a machine with a noise dial you cannot turn below 11, and a readout whose expectation value is biased by the device noise. The move ZNE makes is to turn the dial the wrong way on purpose — up to λ=1,2,3\lambda=1,2,3 — and watch how the answer drifts. If the three readings fall on a smooth, gently decaying curve, that curve carries a memory of where it started: extend it backwards, through the region you were never allowed to visit, to the point λ=0\lambda=0 where the dial would read "no noise at all." The noiseless answer is never measured; it is inferred from the shape of the decay. The first figure makes this concrete with a circuit whose noise is a single depolarizing rate, so the curve is a pure exponential and the extrapolation lands exactly on target.

Figure
Zero-noise extrapolation of a depolarizing-dominated observable E(λ) = (0.9)^λ with ideal value E₀ = 1: three measured points at λ = 1, 2, 3 read 0.900, 0.810, 0.729, and running the fitted curve back to λ = 0 recovers E₀ = 1.000, undoing the native 10 percent bias
The punchline: measure at λ = 1, 2, 3 and run the fit back to λ = 0 to recover E₀ = 1.000 1.00 0.90 0.81 0.729 0 1 2 3 noise scale λ (λ = 1 is native, λ = 0 is ideal) observable E(λ) 0.900 0.810 0.729 E₀ = 1.000 (extrapolated) native bias 0.100

The elegance of the digital version is that you never touch a pulse or trust a noise model. Inserting U†U=1U^\dagger U=\mathbb{1} is a logical no-op — the computed unitary is untouched — but every physical gate you add is another noisy CPTP map, so accumulated decoherence grows in proportion to the gate count. Fold once (n=1n=1) and you traverse the gate set three times, so λ=3\lambda=3; fold twice and λ=5\lambda=5. The second figure shows this ladder of identity-padded circuits and the linearly growing error that feeds the extrapolation.

Figure
Digital unitary folding amplifies noise without changing the logic: replacing UU by U(U†U)nU(U^\dagger U)^n leaves the computed unitary identical while multiplying the physical gate count, and hence the accumulated error, by λ = 2n + 1 — giving the discrete noise scales λ = 1, 3, 5 fed into the extrapolation
Unitary folding: U(U†U)ⁿ = U logically, but physical error ∝ λ = 2n + 1 n = 0 U λ = 1 error ∝ 1 n = 1 U U† U λ = 3 error ∝ 3 n = 2 U U† U U† U λ = 5 each inserted U†U pair adds two passes: λ grows 1 → 3 → 5, logic unchanged

The core equations

The object ZNE manipulates is a single number: the expectation value of a Hermitian observable OO on the state prepared by the noisy circuit. Making the noise scale λ\lambda tunable turns that expectation value into a function E(λ)E(\lambda), which ZNE extrapolates to the zero-noise limit λ=0\lambda=0.

E(λ)=Tr[O Eλ(ρ0)](1)E(\lambda)=\mathrm{Tr}\big[O\,\mathcal{E}_\lambda(\rho_0)\big] \qquad (1)

where: E(λ)E(\lambda) — measured expectation of OO at noise scale λ\lambda; OO — target Hermitian observable; Eλ\mathcal{E}_\lambda — the effective noisy channel of the entire circuit at scale λ\lambda; ρ0\rho_0 — ideal input state; λ\lambda — dimensionless noise scale, λ=1\lambda=1 native and λ=0\lambda=0 ideal. Reading it physically: E(1)E(1) is what the machine produces at its native noise level, and E(0)E(0) is the zero-noise value that cannot generally be sampled directly. The justification for extrapolation is that, in the weak-noise regime and under suitable noise-scaling assumptions, E(λ)E(\lambda) can be expanded smoothly around λ=0\lambda=0. A Dyson expansion of a scaled Lindblad evolution (equation 7) in powers of the small error rate produces such a series, whose constant term is the ideal value.

E(λ)=E0+∑k=1∞ak λk(2)E(\lambda)=E_0+\sum_{k=1}^{\infty}a_k\,\lambda^k \qquad (2)

where: E0E_0 — the ideal noiseless expectation, the extrapolation target; aka_k — expansion coefficients fixed by the circuit, the observable, and the noise generators. The whole game is to estimate E0E_0 — the λ→0\lambda\to0 intercept — from measurements taken only at λ≥1\lambda\ge1, without ever knowing the aka_k. Richardson extrapolation does exactly this: use a degree-nn polynomial through n+1n+1 measured points and evaluate it at λ=0\lambda=0. Because polynomial interpolation is linear in the data, the estimator is a fixed linear combination of the measurements with Lagrange weights.

E^0=∑i=0nγi E(λi),γi=∏j≠iλjλj−λi(3)\hat E_0=\sum_{i=0}^{n}\gamma_i\,E(\lambda_i),\qquad \gamma_i=\prod_{j\neq i}\frac{\lambda_j}{\lambda_j-\lambda_i} \qquad (3)

where: E^0\hat E_0 — the extrapolated zero-noise estimate; λi\lambda_i — the chosen noise scales, i=0,…,ni=0,\dots,n; γi\gamma_i — Lagrange interpolation weights evaluated at λ=0\lambda=0; E(λi)E(\lambda_i) — measured expectation at scale λi\lambda_i. The weights obey ∑iγi=1\sum_i\gamma_i=1, which is the requirement that a noiseless constant is returned unbiased. The amplification map that generates the scales λ\lambda can be purely digital. Global unitary folding replaces the circuit unitary UU by U(U†U)nU(U^\dagger U)^n: logically the identity, physically 2n+12n+1 traversals of the gate set.

fold(U)=U (U†U)n=U  (logically),λ=2n+1(4)\mathrm{fold}(U)=U\,(U^\dagger U)^{n}=U\ \ (\text{logically}),\qquad \lambda=2n+1 \qquad (4)

where: UU — ideal circuit unitary; nn — number of inverse-forward fold pairs; λ=2n+1\lambda=2n+1 — effective noise scale, the number of times the gate set is executed. Per-gate or partial folding folds only a random subset of gates, reaching non-integer λ\lambda so the curve can be sampled more finely. When the noise-scaling protocol is modeled by a depolarizing semigroup, the relevant traceless observable component can have an exponential dependence on the noise scale, E(λ)∝e−cλE(\lambda)\propto e^{-c\lambda}. In the idealized model used here, fitting that exponential form recovers the zero-noise value exactly.

E(λ)=A+B e−Cλ→ depolarizing E(λ)=(1−p)λ E0(5)E(\lambda)=A+B\,e^{-C\lambda}\quad\xrightarrow{\ \text{depolarizing}\ }\quad E(\lambda)=(1-p)^{\lambda}\,E_0 \qquad (5)

where: BB — amplitude of the decaying contribution; CC — effective decay rate; pp — global depolarizing strength per pass; E0E_0 — ideal expectation. For the generic exponential form E(λ)=A+Be−CλE(\lambda)=A+B e^{-C\lambda}, the zero-noise value is E0=A+BE_0=A+B. For the traceless-observable depolarizing model used here, A=0A=0 and E(λ)=(1−p)λE0=e−λ[−ln⁡(1−p)]E0E(\lambda)=(1-p)^\lambda E_0=e^{-\lambda[-\ln(1-p)]}E_0, so C=−ln⁡(1−p)C=-\ln(1-p). The bias reduction is not free. Because the estimator is a signed sum of noisy measurements, its variance is a weighted sum of the individual variances with the squares of the weights — and those squares grow fast with the extrapolation order.

Var[E^0]=∑i=0nγi2 Var[E(λi)](6)\mathrm{Var}[\hat E_0]=\sum_{i=0}^{n}\gamma_i^{2}\,\mathrm{Var}\big[E(\lambda_i)\big] \qquad (6)

where: Var[E^0]\mathrm{Var}[\hat E_0] — variance of the extrapolated estimator; γi2\gamma_i^2 — squared Lagrange weights, which can be ≫1\gg1. This is the cost sheet: every order of bias you cancel is paid for in extra shots to keep the error bars from exploding. Underneath all of it is the continuous-time picture. Analog amplification stretches the gate duration, which is equivalent to multiplying the dissipator of the Lindblad generator by λ\lambda while the Hamiltonian part is held fixed by construction.

dρdt=−iℏ[H,ρ]+λ∑k(LkρLk†−12{Lk†Lk,ρ})(7)\frac{\mathrm{d}\rho}{\mathrm{d}t}=-\frac{i}{\hbar}[H,\rho]+\lambda\sum_k\Big(L_k\rho L_k^\dagger-\tfrac12\{L_k^\dagger L_k,\rho\}\Big) \qquad (7)

where: HH — the ideal gate Hamiltonian; LkL_k — Lindblad jump operators encoding the noise channels; λ\lambda — scaling factor on the dissipator, set physically by pulse-stretching. Expanding the solution in the small dissipative rate is precisely the Dyson series that produces the analytic expansion (2).

Derivation AExponential extrapolation when one depolarizing rate dominates

The goal of this derivation is to prove that a single global depolarizing channel makes the noise-scaled expectation exactly exponential in λ\lambda, so that a two-point exponential fit recovers the ideal value E0E_0 with no residual bias at all — the one regime where ZNE is not an approximation but an identity. We adopt the traceless-observable convention Tr[O 1/2N]=0\mathrm{Tr}[O\,\mathbb{1}/2^N]=0, model the amplified circuit as λ\lambda successive applications of the depolarizing channel with fixed strength pp, and track how the observable rides the coherent part of the state as it is dragged toward the maximally mixed point.

Step 1. Write one pass of the channel, Dp(ρ)=(1−p)ρ+p 1/2N\mathcal{D}_p(\rho)=(1-p)\rho+p\,\mathbb{1}/2^N: the state is linearly mixed toward the maximally mixed state with weight pp. Why: depolarizing is the isotropic, structureless noise model, so it acts by a single scalar contraction of the density matrix toward the center of state space — the cleanest possible one-parameter noise.

Step 2. Take the expectation after one pass, Tr[O Dp(ρ)]=(1−p) Tr[Oρ]+p Tr[O 1/2N]\mathrm{Tr}[O\,\mathcal{D}_p(\rho)]=(1-p)\,\mathrm{Tr}[O\rho]+p\,\mathrm{Tr}[O\,\mathbb{1}/2^N], and drop the second term by the traceless convention, leaving (1−p) Tr[Oρ](1-p)\,\mathrm{Tr}[O\rho]. Why: the maximally mixed component contributes nothing to a traceless observable, so only the coherent part (1−p)ρ(1-p)\rho survives the trace — the entire effect is a multiplicative shrink of EE.

Step 3. Compose λ\lambda passes. Because the mixed component is a fixed point, the channels compose multiplicatively on the coherent part, Dp∘λ(ρ)=(1−p)λρ+(1−(1−p)λ)1/2N\mathcal{D}_p^{\circ\lambda}(\rho)=(1-p)^{\lambda}\rho+\big(1-(1-p)^{\lambda}\big)\mathbb{1}/2^N, so E(λ)=(1−p)λE0E(\lambda)=(1-p)^{\lambda}E_0. Why: iterating a scalar contraction just raises the scalar to a power, which is what turns "one noisy pass" into a clean exponential decay of the observable in λ\lambda.

Step 4. Invert to recover the target, E^0=E(λ) (1−p)−λ\hat E_0=E(\lambda)\,(1-p)^{-\lambda}, so any two scales pin down both pp and E0E_0 exactly: the ratio gives (1−p)λ2−λ1=E(λ2)/E(λ1)(1-p)^{\lambda_2-\lambda_1}=E(\lambda_2)/E(\lambda_1) and then E0E_0 follows. Why: two unknowns need two equations, and because the model is the true channel here, the solution is exact rather than a fit — this is why exponential ZNE is favoured whenever global depolarizing dominates. *As a check: at λ=0\lambda=0, E(0)=E0E(0)=E_0 as required; as λ→∞\lambda\to\infty, E→0E\to0, the maximally mixed limit. For the numeric case p=0.1p=0.1, the ratio E(3)/E(1)=0.81=(0.9)2E(3)/E(1)=0.81=(0.9)^2, confirming the exponential law, and E0=0.9/0.9=1.000E_0=0.9/0.9=1.000.*

Go deeper · Derivation B · 1 minLinear extrapolation kills the leading-order bias
Derivation BLinear extrapolation kills the leading-order bias

The aim here is to show that the simplest nontrivial ZNE — two measurements and a straight line through them, evaluated at zero — already removes the entire O(λ)O(\lambda) error term, leaving a residual bias of only O(λ2)O(\lambda^2), and to expose the fixed factor of variance inflation that comes with it. We start from the analytic expansion E(λ)=E0+a1λ+a2λ2+⋯E(\lambda)=E_0+a_1\lambda+a_2\lambda^2+\cdots, measure at the two smallest convenient scales λ0=1\lambda_0=1 and λ1=2\lambda_1=2, and form the specific linear combination that annihilates a1a_1.

Step 1. Evaluate the expansion at the two scales, E(1)=E0+a1+a2+⋯E(1)=E_0+a_1+a_2+\cdots and E(2)=E0+2a1+4a2+⋯E(2)=E_0+2a_1+4a_2+\cdots. Why: writing both measurements in terms of the same unknown coefficients is what lets us design a combination that cancels the one we do not want, a1a_1, while keeping the one we do, E0E_0.

Step 2. Take the combination 2E(1)−E(2)2E(1)-E(2) with weights γ0=2, γ1=−1\gamma_0=2,\ \gamma_1=-1; note γ0+γ1=1\gamma_0+\gamma_1=1, the consistency condition, and compute 2E(1)−E(2)=2E0+2a1+2a2−E0−2a1−4a22E(1)-E(2)=2E_0+2a_1+2a_2-E_0-2a_1-4a_2. Why: the weights are forced — summing to one preserves E0E_0, and the ratio is fixed by demanding the a1a_1 terms 2a1−2a12a_1-2a_1 cancel; there is no freedom left once the two scales are chosen.

Step 3. Collect terms to get E^0=2E(1)−E(2)=E0−2a2−6a3−⋯\hat E_0=2E(1)-E(2)=E_0-2a_2-6a_3-\cdots: the linear error is gone and the leading residual is second order. Why: a straight line cannot exactly follow a curved E(λ)E(\lambda), so the leading leftover −2a2-2a_2 is the curvature the line misses — reducing the noise-scale spacing suppresses this term quadratically, while adding a third point can cancel the quadratic term.

Step 4. Read off the variance cost from equation (6): Var⁡[E^0]=4Var⁡[E(1)]+Var⁡[E(2)]\operatorname{Var}[\hat E_0]=4\operatorname{Var}[E(1)]+\operatorname{Var}[E(2)]. If the two measurements have equal variance, this is a fivefold variance relative to one measurement. Why: the weight 22 enters squared, so cancelling a1a_1 comes at the cost of increased statistical variance — the simplest instance of the ZNE bias–variance trade-off. *As a check: the weights follow the general Richardson formula with λ0=1,λ1=2\lambda_0=1,\lambda_1=2, giving γ0=λ1/(λ1−λ0)=2\gamma_0=\lambda_1/(\lambda_1-\lambda_0)=2 and γ1=λ0/(λ0−λ1)=−1\gamma_1=\lambda_0/(\lambda_0-\lambda_1)=-1; on the numeric case 2(0.9)−0.81=0.992(0.9)-0.81=0.99, a bias of 0.010.01 versus the native 0.100.10.*

Go deeper · Derivation C · 2 minOrder-n Richardson weights are the zero-noise Lagrange interpolant
Derivation COrder-n Richardson weights are the zero-noise Lagrange interpolant

This derivation generalizes the linear case to arbitrary order, proving that with n+1n+1 scales the Richardson estimator annihilates the first nn error coefficients a1,…,ana_1,\dots,a_n and leaves bias O(λn+1)O(\lambda^{n+1}), while making explicit the price — a sum of squared weights that grows without bound as the order climbs. The key realization is that the whole estimator is nothing but the value at λ=0\lambda=0 of the unique degree-nn polynomial passing through all the measured points, so its coefficients are the Lagrange weights evaluated at the origin.

Step 1. Construct the interpolating polynomial P(λ)=∑i=0nE(λi) ℓi(λ)P(\lambda)=\sum_{i=0}^{n}E(\lambda_i)\,\ell_i(\lambda) with Lagrange basis ℓi(λ)=∏j≠i(λ−λj)/(λi−λj)\ell_i(\lambda)=\prod_{j\ne i}(\lambda-\lambda_j)/(\lambda_i-\lambda_j), which satisfies ℓi(λk)=δik\ell_i(\lambda_k)=\delta_{ik} and therefore passes exactly through the data points. Why: if the true E(λ)E(\lambda) is analytic, a degree-nn interpolant reproduces all polynomial terms through order nn, so the first uncancelled contribution to the zero-noise estimate is of order O(λn+1)O(\lambda^{n+1}).

Step 2. Evaluate at the forbidden point, E^0=P(0)=∑iE(λi) ℓi(0)\hat E_0=P(0)=\sum_i E(\lambda_i)\,\ell_i(0), and identify the weights γi=ℓi(0)=∏j≠i(0−λj)/(λi−λj)=∏j≠iλj/(λj−λi)\gamma_i=\ell_i(0)=\prod_{j\neq i}(0-\lambda_j)/(\lambda_i-\lambda_j)=\prod_{j\neq i}\lambda_j/(\lambda_j-\lambda_i), where the two sign flips per factor cancel. Why: extrapolation to zero is just interpolation evaluated outside the data range, so the same Lagrange machinery that fits the points also delivers the intercept — no new idea, only a substitution λ=0\lambda=0.

Step 3. Substitute the Taylor series (2). Because PP reproduces the first n+1n+1 coefficients, E^0=E0+O(λn+1)\hat E_0=E_0+O(\lambda^{n+1}): the terms a1,…,ana_1,\dots,a_n are all absorbed by the polynomial and cancelled, and only the first uncaptured order survives. Why: a degree-nn polynomial has n+1n+1 free coefficients, exactly enough to match E0E_0 and the first nn error terms, so the residual is the first coefficient it cannot represent, an+1a_{n+1}.

Step 4. Note the stability caveat: ∑iγi=1\sum_i\gamma_i=1 always, but ∑iγi2\sum_i\gamma_i^2 grows steeply with nn — for λ={1,2,3}\lambda=\{1,2,3\} the weights are {3,−3,1}\{3,-3,1\} and ∑iγi2=19\sum_i\gamma_i^2=19. Why: higher order buys geometrically smaller bias but pays geometrically larger variance, so in practice n≤3n\le3 is the working regime and exponential fits are preferred when the noise justifies them. *As a check: the general formula reduces to {2,−1}\{2,-1\} for n=1n=1 and gives exact recovery when E(λ)E(\lambda) is itself a degree-nn polynomial; the λ={1,2,3}\lambda=\{1,2,3\} weights {3,−3,1}\{3,-3,1\} sum to 11 and yield E^0=0.999\hat E_0=0.999 on the numeric case, a bias of 0.0010.001.*

Limiting cases

CaseResultMeaning
Weak global depolarizingE(λ)=E0 (1−p)λE(\lambda)=E_0\,(1-p)^{\lambda}one decay rate; exponential extrapolation is exact and recovers E0E_0 with a single fit parameter.
Perturbative MarkovianE(λ)=E0+∑k=1dakλk+O(λd+1)E(\lambda)=E_0+\sum_{k=1}^{d}a_k\lambda^k+O(\lambda^{d+1})order-dd Richardson cancels the first dd error terms; residual bias is O(λd+1)O(\lambda^{d+1}).
Coherent / non-MarkovianE(λ)E(\lambda) not analytic in λ\lambdascaling does not multiply one rate; the smooth ansatz fails and ZNE leaves uncontrolled bias.
Deep circuit, large λ\lambdaE(λ)→⟨O⟩1/2NE(\lambda)\to\langle O\rangle_{\mathbb{1}/2^N}amplified circuits saturate to the maximally mixed value; there is no slope left to extrapolate.
High extrapolation orderVar[E^0]∼(∑iγi2) Var[E(λi)]\mathrm{Var}[\hat E_0]\sim(\sum_i\gamma_i^2)\,\mathrm{Var}[E(\lambda_i)]weights grow, so variance inflation eventually overwhelms the bias reduction.
No amplificationλ≡1\lambda\equiv1, E^0=E(1)\hat E_0=E(1)with a single point there is nothing to extrapolate; you keep the native bias.

Worked examples

1 · Linear versus Richardson on a depolarizing circuit (fully worked)

A circuit is dominated by global depolarizing noise with per-pass strength p=0.1p=0.1, and the traceless observable has ideal value E0=1E_0=1, so E(λ)=(0.9)λE(\lambda)=(0.9)^{\lambda}. You fold to obtain λ=1,2,3\lambda=1,2,3 and measure E(1)=0.900E(1)=0.900, E(2)=0.810E(2)=0.810, E(3)=0.729E(3)=0.729. Compare the native value, linear ZNE, and order-2 Richardson, and quote the variance cost of each.

Step 1. Establish the baseline. The native reading is E(1)=0.900E(1)=\mathbf{0.900}, a systematic bias of 1−0.900=0.1001-0.900=\mathbf{0.100} — a full 10 percent low. This is the number a naive run reports.

Step 2. Apply linear ZNE with scales λ={1,2}\lambda=\{1,2\} and weights {2,−1}\{2,-1\}: E^0=2E(1)−E(2)=2(0.900)−0.810=1.800−0.810=0.990\hat E_0=2E(1)-E(2)=2(0.900)-0.810=1.800-0.810=\mathbf{0.990}. The residual bias is 1−0.990=0.0101-0.990=\mathbf{0.010}, a tenfold reduction. The variance cost from equation (6) is ∑γi2=22+12=5×\sum\gamma_i^2=2^2+1^2=\mathbf{5\times} a single measurement.

Step 3. Apply order-2 Richardson with scales λ={1,2,3}\lambda=\{1,2,3\}. The Lagrange weights are γ0=λ1λ2(λ1−λ0)(λ2−λ0)=2⋅31⋅2=3\gamma_0=\tfrac{\lambda_1\lambda_2}{(\lambda_1-\lambda_0)(\lambda_2-\lambda_0)}=\tfrac{2\cdot3}{1\cdot2}=3, γ1=λ0λ2(λ0−λ1)(λ2−λ1)=1⋅3(−1)(1)=−3\gamma_1=\tfrac{\lambda_0\lambda_2}{(\lambda_0-\lambda_1)(\lambda_2-\lambda_1)}=\tfrac{1\cdot3}{(-1)(1)}=-3, γ2=λ0λ1(λ0−λ2)(λ1−λ2)=1⋅2(−2)(−1)=1\gamma_2=\tfrac{\lambda_0\lambda_1}{(\lambda_0-\lambda_2)(\lambda_1-\lambda_2)}=\tfrac{1\cdot2}{(-2)(-1)}=1, and indeed ∑iγi=3−3+1=1\sum_i\gamma_i=3-3+1=1. Then E^0=3(0.900)−3(0.810)+1(0.729)=2.700−2.430+0.729=0.999\hat E_0=3(0.900)-3(0.810)+1(0.729)=2.700-2.430+0.729=\mathbf{0.999}, a residual bias of just 0.001\mathbf{0.001}.

Step 4. Weigh the trade-off. Going from linear to order-2 dropped the bias from 0.0100.010 to 0.0010.001 (another tenfold), but the variance factor climbed from ∑γi2=5\sum\gamma_i^2=5 to 32+32+12=193^2+3^2+1^2=\mathbf{19}. Order-2 Richardson recovers E0E_0 to 0.9990.999 at roughly 19×19\times the single-shot variance, versus 0.9900.990 at 5×5\times for linear. Because this noise is purely depolarizing, an exponential fit does better still: E(3)/E(1)=0.729/0.900=0.81=(1−p)2⇒1−p=0.9E(3)/E(1)=0.729/0.900=0.81=(1-p)^2\Rightarrow 1-p=0.9, hence E^0=E(1)/(1−p)=1.000\hat E_0=E(1)/(1-p)=\mathbf{1.000}, exact. The lesson: match the model to the physics before spending shots on high order.

2 · Unitary folding: how many gates, and what noise scale?

A hardware-efficient ansatz has 200200 two-qubit gates, each with an average error of 0.5%0.5\%. You want ZNE at noise scales λ=1,3,5\lambda=1,3,5 by global folding. How many physical gates does each setting execute, what is the approximate accumulated two-qubit error at each, and where does the method break down?

Step 1. Convert the fold count to a scale with λ=2n+1\lambda=2n+1. Native is n=0⇒λ=1n=0\Rightarrow\lambda=1; one fold n=1⇒λ=3n=1\Rightarrow\lambda=3; two folds n=2⇒λ=5n=2\Rightarrow\lambda=5. Each fold appends a full U†UU^\dagger U pair, so the gate count is λ×200\lambda\times200: 200200, 600600, and 10001000 two-qubit gates respectively. The computed unitary is identical in all three; only the physical length grows.

Step 2. Estimate accumulated error as λ×(per-pass error)\lambda\times(\text{per-pass error}). One native pass carries roughly 200×0.005=1.0200\times0.005=1.0 expected two-qubit faults, so the amplified circuits carry about 1.01.0, 3.03.0, and 5.05.0 expected faults. The effective error rate scales linearly with λ\lambda, which is exactly the single-parameter behaviour the extrapolation assumes.

Step 3. Locate the breakdown. Once the accumulated error becomes large, the state can be driven toward a highly mixed state, so E(3)E(3) and E(5)E(5) may have little signal left and the analytic expansion (2) may cease to be a useful local model — the curve flattens and the extrapolation loses its lever arm. With ∼1\sim1 fault already at λ=1\lambda=1, folding to λ=5\lambda=5 leaves the observable near its noise floor, so this circuit is at the edge of where ZNE helps: mitigation works when the native circuit is mildly noisy, not when it is already scrambled. In practice one would fold more gently (partial folding for non-integer λ\lambda near 1) or shorten the ansatz first.

Physics of the portrait

Zero-noise extrapolation looks like a quantum-hardware invention of the late 2010s, but its mathematical heart is a century older — a numerical-analysis trick for wringing an exact answer out of a sequence of approximate ones, first written down to solve partial differential equations by hand.

Origins · 2 minLewis Fry Richardson and how this was discovered
LF
Lewis Fry Richardson
1881–1953 · numerical extrapolation to the limit, a century later repurposed as quantum error mitigation

Lewis Fry Richardson was an English mathematician, physicist, and meteorologist whose ideas kept arriving decades ahead of the machinery needed to use them. In a 1911 paper on the finite-difference solution of physical problems, and more fully in the 1927 study "The deferred approach to the limit" with J. Arthur Gaunt, he introduced what is now called Richardson extrapolation: compute a quantity at several finite step sizes hh, recognize that the answer depends analytically on hh with a known leading power, and combine the estimates so the leading error term cancels — extrapolating to the unreachable limit h→0h\to0. He used it to squeeze accuracy out of the crude arithmetic of the day, and famously imagined a "forecast factory" of thousands of human computers integrating the weather equations, a vision realized only when electronic computers arrived. Richardson's method became a staple of numerical analysis (it is the engine inside Romberg integration and the Bulirsch-Stoer ODE solver), and his name also attaches to the Richardson number in turbulence and to his pioneering, Quaker-pacifist-motivated statistical study of the "deadly quarrels" of war. The quantum connection came in 2017, when Kristan Temme, Sergey Bravyi, and Jay Gambetta — and independently Ying Li and Simon Benjamin — realized that a noisy quantum circuit's expectation value depends analytically on the noise strength exactly as Richardson's finite-difference answer depends on the step size hh. Replace hh by the noise scale λ\lambda, run the circuit at several amplified noise levels, and Richardson's 1911 combination extrapolates to the noiseless limit. A hand-computation trick for weather forecasting became the first practical error-mitigation method for quantum computers.

History
1911
— Richardson introduces finite-difference extrapolation to improve approximate solutions of physical problems
1927
— Richardson and Gaunt formalize the deferred approach to the limit
2017
— Temme, Bravyi, and Gambetta introduce zero-noise extrapolation and probabilistic error cancellation as error-mitigation methods for short-depth quantum circuits
2017
— Li and Benjamin independently propose error extrapolation via stretched evolution
2019
— Kandala and coworkers at IBM demonstrate zero-noise extrapolation on a superconducting quantum processor, extending the computational reach of variational quantum eigensolvers
2020
— Giurgica-Tiron and colleagues introduce digital unitary folding in the Mitiq toolkit

Going further

Texts. Nielsen and Chuang, Quantum Computation and Quantum Information, §8.3 (quantum operations) and Ch. 10 (error-correction context) · Preskill, Lecture Notes on Quantum Computation, Ch. 3 (noise and master equations) · Watrous, The Theory of Quantum Information, Ch. 4 (channels) · Wilde, Quantum Information Theory, Ch. 4 · Kitaev, Shen, and Vyalyi, Classical and Quantum Computation (channels and errors) · original sources: Temme, Bravyi, and Gambetta, Phys. Rev. Lett. 119, 180509 (2017) and Li and Benjamin, Phys. Rev. X 7, 021050 (2017). The key numerical anchors are: the Lagrange/Richardson weights γi=∏j≠iλj/(λj−λi)\gamma_i=\prod_{j\neq i}\lambda_j/(\lambda_j-\lambda_i) satisfy ∑iγi=1\sum_i\gamma_i=1; the linear weights are {2,−1}\{2,-1\} for λ={1,2}\lambda=\{1,2\} and the order-2 weights are {3,−3,1}\{3,-3,1\} for λ={1,2,3}\lambda=\{1,2,3\}, with ∑γi2=5\sum\gamma_i^2=5 and 1919 respectively; the depolarizing composition Dp∘λ\mathcal{D}_p^{\circ\lambda} gives E(λ)=(1−p)λE0E(\lambda)=(1-p)^\lambda E_0; and the worked-example arithmetic for p=0.1p=0.1 yields native 0.9000.900, linear 0.9900.990, Richardson 0.9990.999, exponential 1.0001.000. The folding gate counts (200200, 600600, 10001000) and fault estimates are illustrative of the stated 200200-gate, 0.5%0.5\%-error setup.

Prerequisites. Quantum Operations and CPTP Maps · Kraus Representation of Quantum Channels · Lindblad / GKSL Master Equation · The Depolarizing Channel. Leads to. Markovian vs Non-Markovian Dynamics · Master-Equation Derivations: Born-Markov and Secular Approximations · Decoherence.

Ask the tutor: Ask me to work a ZNE fit on your own measured points, to derive the Richardson weights for any set of noise scales, or to judge whether your circuit's noise is analytic enough for extrapolation to be trustworthy.

Isn't ZNE just quantum error correction without the overhead?

No, and the distinction is the whole point. Error correction removes errors: it encodes logical information across many physical qubits, detects errors with syndrome measurements, and corrects them, so a single logical result comes out clean — at large qubit and gate cost. ZNE corrects nothing on any individual run and uses no extra qubits or syndromes. Instead it treats the biased expectation value statistically, models how that bias grows with noise, and extrapolates it away. So ZNE can only sharpen an average over many shots, never protect a fragile state mid-computation, and it does not scale to arbitrarily deep circuits. It is a near-term expedient that buys accuracy today; error correction is the long-term solution ZNE explicitly does not replace.

How do you actually turn the noise up on real hardware?

Two families of methods. Analog amplification stretches the control pulses: run each gate more slowly so the qubit spends proportionally longer exposed to decoherence, which multiplies the Lindblad dissipator by λ\lambda while leaving the intended rotation unchanged. This needs low-level pulse access and careful calibration. Digital amplification uses unitary folding: insert identity-equivalent sequences U†UU^\dagger U into the circuit, either globally (U↦U(U†U)nU\mapsto U(U^\dagger U)^n, giving λ=2n+1\lambda=2n+1) or gate-by-gate. Folding requires only the native gate set and no hardware model, which is why it dominates in practice and is built into libraries like Mitiq. Partial or "local" folding of a random subset of gates reaches non-integer λ\lambda close to 11, letting you sample the curve finely near where it matters most.

Which extrapolation model should I choose — linear, Richardson, or exponential?

Match the model to the noise, not to how tight the fit looks. If the dominant error is close to global depolarizing — often a decent approximation for deep, twirled circuits — the true curve is E(λ)=(1−p)λE0E(\lambda)=(1-p)^\lambda E_0, so an exponential fit is unbiased and should win. If the noise is weak and Markovian but not obviously depolarizing, a low-order Richardson extrapolation (order 1 or 2) cancels the leading error terms with modest variance. Reserve higher-order Richardson for when you have the shot budget to absorb the growing weights. The trap is over-fitting: a flexible model hugs the data points, but a wrong functional form extrapolates badly to λ=0\lambda=0. Start simple, check stability against adding a point, and stop before the error bars overtake the bias you removed.

Why do the error bars explode at high extrapolation order?

Because the estimator is a signed sum of noisy measurements, and its variance uses the squares of the weights: Var[E^0]=∑iγi2 Var[E(λi)]\mathrm{Var}[\hat E_0]=\sum_i\gamma_i^2\,\mathrm{Var}[E(\lambda_i)]. The Richardson weights must sum to one to keep the estimate unbiased, but to cancel more error terms they take large positive and negative values that nearly cancel — for λ={1,2,3}\lambda=\{1,2,3\} they are {3,−3,1}\{3,-3,1\}, so ∑γi2=19\sum\gamma_i^2=19, already a nineteenfold variance inflation over a single measurement. Add another scale and the weights grow again. So each extra order shrinks the bias geometrically but grows the variance geometrically too. Past a certain order the added statistical noise is larger than the systematic bias you cancelled, and the estimate gets worse. This is why practitioners rarely go beyond order 3 and lean on exponential fits when the physics allows.

Does ZNE fix readout (measurement) error?

Not on its own, and assuming it does is a common way to get a wrong answer. Unitary folding amplifies gate noise by lengthening the circuit, but a measurement happens once at the end regardless of how many gates precede it, so its confusion-matrix error does not scale with λ\lambda the way the extrapolation assumes. Feeding unscaled readout bias into a ZNE fit corrupts the intercept. The standard fix is to layer methods: apply readout-error mitigation (invert the measurement confusion matrix) first, then apply ZNE to the residual gate noise, and add probabilistic error cancellation where a noise model is available. ZNE is one tool in a stack, not a universal solvent.

Can the noise scale λ be non-integer, or less than 1?

Non-integer above 1: yes, and it is useful. Global folding gives only odd integers λ=2n+1\lambda=2n+1, but partial folding — inserting U†UU^\dagger U around a random subset of gates — produces any real λ≥1\lambda\ge1, letting you sample the curve densely near λ=1\lambda=1 where the extrapolation lever is strongest. Below 1: essentially no. Reaching λ<1\lambda<1 would mean running with less noise than the hardware's native level, which folding cannot do (you can add identity padding but not subtract physical error) and pulse-stretching cannot do either. Effective noise reduction is the province of probabilistic error cancellation, a different technique; ordinary ZNE lives entirely in λ≥1\lambda\ge1 and extrapolates outward to 00.

When does ZNE fail outright, even with unlimited shots?

Whenever its two structural assumptions break. First, if the noise is not faithfully described by a single scalable scalar λ\lambda — coherent over-rotations, non-Markovian memory, cross-talk, or leakage — then amplifying "the noise" does not simply multiply one rate, E(λ)E(\lambda) is not analytic in λ\lambda, and the smooth ansatz extrapolates to the wrong intercept no matter how many shots you take. Second, if the native circuit is already so deep that even λ=1\lambda=1 has driven the observable near its maximally mixed floor, there is no slope left: the amplified points carry almost no signal. Unlimited statistics cannot repair a biased model or resurrect a scrambled signal — the failures are systematic, not statistical. ZNE is trustworthy in the mild-noise, weakly-Markovian regime and should be sanity-checked against a smaller classically simulable instance whenever you leave it.

Before this

Where this leads