Noise

qliff.noise simulates noisy circuits inside the stabilizer formalism. Each channel is a quasiprobability mixture of stabilizer (Clifford) channels:

$$\mathcal{E}(\cdot )=\sum _{\mu }{q}_{\mu }{\mathcal{S}}_{\mu }(\cdot ).$$

A sampler runs many pure-Clifford trajectories and reweights each by $\mathrm{sign}(q_{\mu })\gamma$, where $\gamma =\sum _{\mu }|q_{\mu }|$ is the per-location overhead. Pauli channels have $q_{\mu }\ge 0$ summing to one, so $\gamma =1$. Coherent and non-unitary channels carry negative weights and $\gamma >1$. This is the method of arXiv:2512.07304.

You rarely touch a channel directly. Add noise to a Circuit, then call estimate or sample -- the classes below are what those methods build.

Channel

Channel is an abstract base; subclasses expose a stabilizer-channel decomposition. A branch is a pair (weight, ops), where ops is a list of (gate, targets). The identity (no-fault) branch comes first.

Property/Method	Description
is_pauli	True if all weights are probabilities (records are sampleable)
is_unitary	True for a coherent unitary channel, whose one Kraus operator carries interfering amplitudes
arity	number of qubits the channel acts on
gamma	the sampling overhead $\gamma =\sum _{\mu }|q_{\mu }|$; 1.0 for a Pauli channel
branches(targets)	the (weight, ops) decomposition for the given qubits
sample(targets, rng)	draw one branch as (sign(w) * gamma, ops)
kraus_branches(targets)	Kraus operators, each a sum of (amplitude, ops) Pauli terms
pauli_twirl(targets)	(px, py, pz), the lossy Pauli approximation; single-qubit only

branches is the sampling decomposition. The last two rows are the other views a decoder may need: Kraus operators (exact for coherent rotation, whose single operator kraus_branches(targets)[0] holds the interfering amplitude terms, and exact for damping via the doubled process-matrix (chi) decoder), and the twirl (lossy, and only where a Pauli-only decoder is asked for).

Channel catalog

Every noise instruction resolves to a Channel through make_channel(name, arg), which reads the CHANNEL_META registry; build the object directly when you need one outside a circuit. The first seven rows are PauliChannels (weights are probabilities); the last three are general (is_pauli = False) and carry signed quasiprobabilities. Rotation decomposes over the Cliffords diagonal in its own axis, $\{I,Z,S,S^{†}\}$; the damping channels decompose over $\{I,Z,\text{reset}\}$, plus a reset to $|1⟩$ when $\tau <1-\mu$.

Channel	Circuit method	Arguments	Description
make_channel("DEPOLARIZE1", p)	DEPOLARIZE1(q, p)	$p$, dict or vector	$X,Y,Z$ each with probability $p/3$
make_channel("DEPOLARIZE2", p)	DEPOLARIZE2(pair, p)	$p$, dict or vector	the 15 two-qubit Paulis, each $p/15$
PauliChannel({"X": p})	X_ERROR(q, p)	$p$	$X$ with probability $p$
PauliChannel({"Y": p})	Y_ERROR(q, p)	$p$	$Y$ with probability $p$
PauliChannel({"Z": p})	Z_ERROR(q, p)	$p$	$Z$ with probability $p$
PauliChannel({"X": .., "Z": ..})	PAULI_CHANNEL_1(q, w)	$p$, dict or vector	the same channel as DEPOLARIZE1
PauliChannel({"ZZ": .., "XI": ..})	PAULI_CHANNEL_2(pair, w)	$p$, dict or vector	the same channel as DEPOLARIZE2
Rotation(axis, theta)	RZ(q, theta) / RX(q, theta) / RY(q, theta)	axis, $\theta$	coherent rotation $e^{-i\theta P/2}$
AmplitudeDamping(p)	AMPLITUDE_DAMP(q, p)	$p$	energy decay $|1⟩\to |0⟩$; returns a GeneralizedDamping
GeneralizedDamping(lam, mu, tau)	c.noise(ch, q)	$\lambda ,\mu ,\tau$	combined amplitude + phase damping; build from $T_1,T_2$ with from_times

Anisotropy

Isotropic is the default. A bare scalar spreads evenly: $p/3$ over $X,Y,Z$ for the one-qubit channels, $p/15$ over the fifteen pairs for the two-qubit ones. That is what almost every circuit wants, and it is unchanged.

When you want anisotropy, pass a {label: rate} dict instead of the scalar. Anything unlisted is zero:

c.DEPOLARIZE1(q, {"Z": 0.02, "X": 0.001})   # dephasing-dominated
c.DEPOLARIZE2([a, b], {"ZZ": 0.02})         # pure ZZ crosstalk
PauliChannel({"ZZ": 0.02})                  # the channel object directly

DEPOLARIZE1 and PAULI_CHANNEL_1 are the same channel under two names, as are DEPOLARIZE2 and PAULI_CHANNEL_2; the depolarizing name reads better where the rate is isotropic, the Pauli name where it is not. Both take a scalar, a dict, or the dense vector in canonical order

$$IXIYIZXIXXXYXZYIYXYYYZZIZXZYZZ$$

(the dense form exists for interop; prefer the dict, since nobody can read fifteen positional floats and a mis-ordered one is silent). A mistyped, wrong-case or wrong-arity label raises naming the problem rather than landing as a zero rate.

PauliChannel(rates) reads the arity off the key width: {"X": ..} is one qubit, {"ZZ": ..} is two. Mixed widths are refused rather than guessed at.

X_ERROR, Y_ERROR and Z_ERROR carry a single Pauli each, so they take a scalar only and refuse a vector. A 3-tuple handed to X_ERROR used to be re-read as the whole channel, silently turning a bit flip into a general Pauli channel.

For the common single-axis case the code builders take bias, the ratio $\eta =p_{\text{axis}}/(p_b+p_c)$ of the biased axis to the other two, and bias_axis, instead of a dict. The split itself is qliff.noise.channel.bias_split(p, eta, axis="Z"), returning (px, py, pz) for axis in "X", "Y", "Z":

$$p_{\text{axis}}=p\frac{\eta }{1+\eta },p_{\text{other}}=\frac{p}{2(1+\eta )}\text{each},$$

so $\eta =\frac{1}{2}$ is depolarizing on any axis and $\eta \to \mathrm{\infty }$ is pure single-axis noise. A negative $\eta$ raises. The axis is a knob rather than a constant because $X$ and $Z$ are duals everywhere else: swapping bias_axis and memory together gives a bit-identical detector error model, which a $Z$-locked bias made impossible. bias / bias_axis are accepted by every code builder -- see the noise bias.

bias shapes the 1-qubit Pauli channels only. There is no canonical two-qubit lift of a one-qubit $\eta$ (the tensor square of a biased 1Q channel does not reduce to DEPOLARIZE2 at $\eta =\frac{1}{2}$), so asking for it on a two-qubit channel raises instead of quietly meaning something else; pass a dict.

Amplitude damping

AmplitudeDamping(p) decomposes exactly over $\{I,Z,\text{reset}\}$:

$$q_I=\frac{(1-p)+\sqrt{1-p}}{2},q_Z=\frac{(1-p)-\sqrt{1-p}}{2}<0,q_R=p.$$

The overhead $\gamma =p+\sqrt{1-p}\approx 1+\frac{p}{2}$ sits barely above the Pauli value of $1$, so this non-unitary channel costs nearly as little to sample as Pauli noise. It is not a Pauli channel -- estimate it with expect (estimate reweights automatically).

GeneralizedDamping(lam, mu, tau) is the $T_1$/$T_2$ generalisation, built from times with GeneralizedDamping.from_times(t1, t2, t). It supplies its own Kraus decomposition, so it works with the doubled process-matrix (chi) decoder and with pauli_twirl(); the twirl is $p_x=p_y=(1-\mu )/4$, $p_z=(1+\mu -2\lambda )/4$. AmplitudeDamping(p) is this channel at $\lambda =\sqrt{1-p}$, $\mu =1-p$, $\tau =p$.

Rotation

Rotation(axis, theta) folds the nearest Clifford power $S^k$ in exactly and mixes only the residual angle $r$, $|r|\le \pi /4$, so $\gamma =\cos r+|\sin r|\le \sqrt{2}$ at any $\theta$. At a multiple of $\pi /2$ it is the Clifford gate itself (see Circuit).

Rates are validated

Channel rates are checked on construction: negative rates, rates above 1, and a rate set summing above 1 raise with the offending value named, as do non-CP $(\lambda ,\mu ,\tau )$. A dense vector of the wrong length is refused naming the count it wanted, an unknown dict label naming the labels it allows, and a vector handed to a single-Pauli channel as not a number. An out-of-range rate used to produce a non-physical map that stayed trace-preserving, so nothing downstream noticed and the sampler returned a confident wrong number.

This applies to rates, not to the signed decompositions: Rotation and GeneralizedDamping in its signed regime legitimately carry negative quasiprobabilities and are unaffected.

Leakage

Leakage is the one hardware mechanism here that is not a channel. A leaked qubit has no qubit density matrix, so there is nothing to decompose: no CHANNEL_META entry, no stabilizer-branch decomposition, no LEAKAGE instruction and no detector error model. It is configured with a LeakageModel object handed to a LeakySimulator, which a WeightedDetectorSampler runs once per shot, never with a circuit.append(...) noise instruction.

What makes it tractable anyway is that coherence between the computational levels and the leaked level decays on the anharmonicity timescale, far faster than a syndrome round. Which qubits are leaked is therefore a classical fact you can sample, and conditioned on it everything else is an ordinary qubit process again. One bit per qubit suffices even though the transition involves a pair, because the post-transition state is a product state.

The mechanism

The leak rides the same $|11⟩\to |02⟩$ avoided crossing the two-qubit gate itself uses, so it fires once per two-qubit gate with probability leak times the pair's $|11⟩$ occupation. That occupation is read exactly off the tableau, without disturbing it, as

$$P_{11}=\frac{⟨(I-A)(I-B)⟩}{4}=\frac{1-⟨A⟩-⟨B⟩+⟨AB⟩}{4},$$

which on a stabilizer state is one of $0$, $\frac{1}{4}$, $\frac{1}{2}$, $1$. The hardware gate is a CZ, so under CX the target sits in the H-conjugated frame and its excitation operator is $X$ while the control's is $Z$.

On the event, one qubit leaves the subspace and the other drops to its ground state. The qubit that does not leak therefore takes a real error despite never leaving the subspace. site selects which member leaks: "target" for $|02⟩$, "control" for $|20⟩$. It is a fixed property of the gate, set by whose 1-2 transition the pulse uses, not a random choice.

While a qubit is leaked

Operation	What happens
single-qubit gate	skipped, the leaked qubit is decoupled from the tableau
two-qubit gate	a rotation is applied to the healthy partner instead
measurement	a biased random bit, governed by readout
reset, or measure-and-reset	the flag clears

readout is P(measure 1 given leaked). 0.5 models a discriminator that cannot separate the leaked level from the code states, 1.0 one that reports every leaked qubit as excited. The record still gets exactly one bit per measured qubit either way.

The reset row is the one that matters in practice. Ancillas are reset every round, so ancilla leakage self-clears. Data qubits are not, so data leakage persists.

Seeping back happens with probability seep per two-qubit gate the leaked qubit takes part in, and returns the qubit in $|1⟩$, which is where the leaked level decays to. The onward $|1⟩\to |0⟩$ relaxation is ordinary $T_1$ and belongs to AMPLITUDE_DAMP, not here.

The partner kick is coherent, not stochastic

phi is the per-gate angle the healthy partner is rotated by. Under kick="COHERENT" (the default) it is applied as the Rotation it is, decomposed through the same quasiprobability path as every other non-Pauli channel, so no twirl enters anywhere and each trajectory carries an importance weight. That weight is exactly 1.0 on any shot where nothing leaked, so negativity is paid only where leakage actually happened. kick="TWIRL" replaces the rotation with its Pauli approximation and gives unit weights throughout. It is kept for comparison with the literature, not because it is right.

At distance 3 the erasure of the leaked qubit accounts for the effect, and the rotation's contribution was not resolvable at that shot budget. No difference in logical error rate is claimed between the two arms here.

Argument	Default	Meaning
leak	0.0	P(leak) per two-qubit gate, times the pair's $|11⟩$ occupation
seep	0.0	P(return to $|1⟩$) per two-qubit gate a leaked qubit joins
phi	0.0	per-gate rotation angle imprinted on the healthy partner
readout	0.5	P(measure 1 given leaked)
kick	"COHERENT"	"COHERENT" or "TWIRL"
site	"target"	which member of the pair leaks: "target" or "control"

Rates outside $[0,1]$ and unknown kick / site names are refused at construction. LeakySimulator is the per-shot object underneath, a Simulator carrying one leaked flag per qubit, and sim.leaked lists the qubits currently outside the computational subspace.

Sampling

WeightedDetectorSampler(circuit, simulator) takes any simulator(num_qubits, seed) factory and runs one per shot. Its sample(shots, seed) returns three values where DetectorSampler.sample returns two: detection events, observable flips, and the per-trajectory importance weights.

from functools import partial

from qliff.noise import LeakageModel, LeakySimulator
from qliff.qec import WeightedDetectorSampler, rotated_surface_code
from qliff.qec.decoder import make

circ = rotated_surface_code(3, 3, 0.001)
model = LeakageModel(leak=0.02, seep=0.05, phi=0.3)
leaky = WeightedDetectorSampler(circ, partial(LeakySimulator, model=model))

dets, obs, weights = leaky.sample(300, seed=0)

# the circuit's own leakage-free error model: what hardware forces on a decoder
decoder = make("bposd", circ)
wrong = (decoder.decode_batch(dets) != obs).any(axis=1)

ler = float((weights * wrong).sum() / weights.sum())

The weighted mean is unbiased under the true leakage process. Under kick="TWIRL" every weight is 1.0 and it degenerates to the plain failing fraction. At leak=0 it reproduces DetectorSampler(circuit, Simulator), the plain per-shot Python sampler, bit for bit.

Why no detector error model can represent it

A DetectorErrorModel can only produce syndromes in the column space of its check matrix $H$ over ${\mathbb{F}}_2$. Leakage produces syndromes outside it. Measured on a distance-3 rotated surface code memory over 3 rounds, where $\mathrm{rank}(H)=12$ of 16 detector dimensions, at 4000 shots per row:

noise	fraction of syndromes outside $\mathrm{col}(H)$
DEPOLARIZE1 $p=0.001$	0
DEPOLARIZE1 $p=0.05$	0
DEPOLARIZE1 $p=0.15$	0
leakage, leak=0.005	3.3%
leakage, leak=0.02	12.6%
leakage, leak=0.05	28.0%
leakage, leak=0.10	45.8%

The leakage rows used seep=0.05, phi=0.3 and the coherent kick.

This is not a statement about miscalibrated priors. The syndrome is not in the image of the model at any rates, so no reweighting and no fitted detector error model reaches it. Minimum-weight matching does not misdecode these syndromes, it fails to return a matching at all, and pymatching raises. Use bposd for leakage work, which returns a best-effort answer instead of raising. The twirled and coherent arms give the same escape fraction to within sampling error, which locates the effect in the classical flag rather than in the rotation.

Limitations

• No coherence between the computational and the leaked levels. This is the load-bearing assumption.
• The branch in which the gate does not leak carries a back-action $\sqrt{I-p|11⟩⟨11|}$, which is not a stabilizer operation. It is neglected, an $O(p)$ effect per gate.
• Leakage is injected only at two-qubit gates. Real devices also leak on single-qubit gates, while idling, and at readout.
• Seepage is evaluated per two-qubit gate rather than per unit time.
• Only the first leaked level is modelled.
• A SWAP involving a leaked qubit does not transport the flag.
• This is the per-shot Python path with no batched core behind it, so it is the slowest sampler in the package.

Sampler

Sampler wraps a circuit and runs trajectories. It has two methods.

Method	Handles	Notes
expect(obs, shots, seed=None, stratify=False)	any channel	importance estimate, unbiased for any noise
sample(shots, seed=None)	Pauli only	uint8 array (shots, measurements); raises on a general channel

expect draws one branch per noise location and reweights each trajectory by $\mathrm{sign}(q_{\mu })\gamma$. On Pauli noise every weight is $1$, so it is plain Monte-Carlo. With stratify=True it rewrites the estimate as $F=\sum _{k}P(k)F_k$: $P(k)$ is a Poisson-binomial probability of $k$ faulty locations, and $F_k$ the conditional estimate. Since $F_k$ varies slowly with $k$, this cuts variance sharply at the same shot budget.

from qliff import Circuit
from qliff.noise import Sampler

c = Circuit(1)
c.H(0).RZ(0, 0.3)

Sampler(c).expect("X", 20000)
Sampler(c).expect("X", 20000, stratify=True)

sample is Pauli-only: a measured bitstring cannot be reweighted by a negative quasiprobability, so a non-Pauli channel raises.

TIP

Prefer c.estimate(observable, shots). It uses flat importance sampling for Pauli circuits and stratifies otherwise, so you rarely build a Sampler by hand.

Rare-event splitting

Direct Monte-Carlo needs on the order of $1/\text{LER}$ shots to see one failure, so a logical error rate of ${10}^{-15}$ is unreachable at any shot budget. SplittingEstimator spends a fixed decode budget across a ladder of physical error rates instead: it factors the tail probability into a product of conditional level-crossing ratios, each estimated by a Metropolis walk over the decoder's failing set. The method is Bravyi-Vargo multilevel splitting (arXiv:1308.6270); the estimate is unbiased.

from qliff.noise import SplittingEstimator
from qliff.qec import rotated_surface_code

est = SplittingEstimator(lambda p: rotated_surface_code(5, 5, p), 0.002)
res = est.estimate(p_top=0.04, p_target=0.002, levels=9)   # (ler, rel_err, decodes)

res.ler        # logical error rate at the deepest ladder rate
res.rel_err    # relative standard error, in quadrature over the levels
res.decodes    # decoder calls spent

Argument	Meaning
circuit_fn(p)	builder returning the memory Circuit at physical rate p
decoder_ref	rate at which the fixed decoder is built
decoder_cls	any batch decoder; MwpmDecoder by default
ladder	explicit strictly-descending rate list, in place of p_top/p_target/levels
top_shots	shots for the direct-MC measurement at the top of the ladder
burn, samples	length of each per-level Metropolis chain

The decoder is built once, at decoder_ref, and held fixed down the ladder. The failing set it induces then depends on the code and the decoder alone and not on $p$, which is what makes splitting in $p$ valid. Build it at a rate whose error model contains every mechanism the ladder will use; the top rate is always safe. A signature absent from the reference model raises rather than misaligning the priors. est.priors_at(p) gives the per-mechanism prior vector at any rate, in that model's mechanism order. est.direct_mc(p, shots) runs plain Monte-Carlo with the same fixed decoder and returns (failures, last failing fault vector or None), for cross-checking wherever MC is affordable.

Measured on rotated-surface patches under the amplitude damping channel's Pauli twirl, the log-log slope of LER against the damping rate reproduces $\lceil d_{\text{fail}}/2\rceil$ exactly. $d_{\text{fail}}$ is the patch dimension the protected logical does not run along (see surface variants):

$d_{\text{fail}}$	2	3	4	6	8	12
$\lceil d_{\text{fail}}/2\rceil$	1	2	2	3	4	6
measured slope	1.00	1.99	1.91	2.92	3.97	6.01

The last column is a $2\times 12$ patch, which reaches $\sim {10}^{-23}$ in seconds. A 2-million-shot signed importance estimate of the same quantity under amplitude damping returns noise instead, and can come back negative.

The generator is a local sweep script, not tracked in the repository, which builds each circuit as PAULI_CHANNEL_1 carrying AmplitudeDamping(gamma).pauli_twirl((0,)). It runs the twirl rather than the channel because the estimator is Pauli-only, for the reason below. The twirl is validated rather than assumed: against well-resolved direct Monte Carlo on the real amplitude damping channel it agrees to within 3-8%.

WARNING

SplittingEstimator is Pauli-only, and the twirl it falls back on is a modelling approximation. It works from a DetectorErrorModel, which has nowhere to put a signed or complex branch weight, so circuit_fn must return a Pauli circuit. Under a non-Pauli channel what is being estimated is the Pauli twirl of the circuit, not the circuit. The twirl gets the decoding decisions right but underestimates the logical error rate: true / twirled is ~1.35x under amplitude damping at $d=3$ and ~2.78x ($25/9$) under coherent rotation, growing with distance (~4.23x at $d=5$). Measured against well-resolved direct MC on the real channel at 15 points, the twirled splitting estimate agreed within 0.91x to 1.31x. Do not read a threshold or a $\mathrm{\Lambda }$ straight off it.

Custom channels

Subclass Channel, set is_pauli (and arity = 2 for a two-qubit channel), and return the branches: an identity branch first, then (weight, ops) faults (weights may be negative quasiprobabilities). No Rust or recompilation needed. Drop it into a circuit with c.noise(channel, q).

from qliff import Circuit
from qliff.noise import Channel


class Dephase(Channel):
    is_pauli = True

    def __init__(self, p):
        self.p = p

    def branches(self, targets):
        q = targets[0]

        return [(1.0 - self.p, []), (self.p, [("Z", (q,))])]


c = Circuit(1)
c.H(0).noise(Dephase(0.1), 0)

A custom channel is classified by what it actually emits, not by its name. If every live branch commutes with a memory basis, qliff.noise.channel.blind_axis(channel, arg) reports that basis and the QEC builders warn -- see the memory basis.