Understanding Fault-tolerant Threshold Theorem in Practice

Fault-tolerant Threshold Theorem

The threshold theorem is the mathematical bedrock of scalable quantum computing. Intuitively, it states that a fault-tolerant quantum computation of any size can be accurately executed on imperfect hardware, provided that the base error rate of the physical operations, \(p\), remains strictly below a specific, non-zero constant known as the threshold, \(p_{th}\), the rightmost crossing point in the schematic above.

More rigorously, assuming a local stochastic error model where \(p < p_{th}\), we can take any ideal circuit \(\mathcal{C}\) of size \(N\) and construct a corresponding fault-tolerant circuit \(\mathcal{C}^{\prime}\). Even when subjected to continuous noise, this corrected circuit is guaranteed to yield an output that deviates from the ideal outcome by no more than an arbitrarily small tolerance, \(\epsilon > 0\). Crucially, the theorem ensures that this correction is practically achievable and the required hardware overhead, i.e., the total number of physical qubits and time steps needed for \(\mathcal{C}^{\prime}\), grows at most by a polylogarithmic factor, \(\mathcal{O}(\log^{c}(N/\epsilon))\) for some positive constant \(c\) [1].

In other words, operating in the green region (where \(p < p_{th}\)) of the schematic is not merely better, it is the condition under which arbitrarily large and reliable quantum circuits become physically constructable. Although the original theoretical framework relied on specific assumptions like independent stochastic noise, the threshold theorem is robust enough to apply to highly realistic, correlated noise environments as well. Moreover, it assures that the required number of physical qubits grows non-exponentially with the size of the computation, which means there is no fundamental physical barrier standing in the way of large-scale quantum computers. At least, theoretically!

To test the threshold theorem in practice, we use one of the initial candidates for quantum error correction, the surface code, which is a topological code where qubits are arranged on a 2D grid, with stabilizer measurements locally checking for parity errors among nearest neighbors. If you had like to deepen your understanding of surface codes before proceeding, our Introducing Lattice Surgery and A Game of Surface Codes demos are a great starting point. Here, we specifically look at its rotated variant, which requires only \(d^2\) data qubits to achieve the same distance \(d\). This gives a 50% reduction in qubit overhead when compared to the standard surface codes.

As shown above for \(d=3\) and \(d=5\) rotated surface codes, the bulk plaquettes alternate between \(Z\)- and \(X\)-type stabilizers in a checkerboard pattern, while the boundaries host weight-2 stabilizers: \(X\)-type on the top and bottom edges, \(Z\)-type on the left and right edges [3]. The logical \(X\) operator runs top-to-bottom along the left column, and the logical \(Z\) operator runs left-to-right along the top row, each forming a string of length \(d\) that cannot be detected by any stabilizer.

The function below encodes this geometry directly, returning the qubit indices involved in each stabilizer and logical operator for any distance \(d\) rotated surface code:

import numpy as np

def rotated_surface_code(d: int):
    """Generate stabilizers and logical operators for a rotated surface code."""
    # Create a grid of row and column indices for the plaquettes
    grid_idxs = np.arange(d - 1)
    row, col = np.meshgrid(grid_idxs, grid_idxs, indexing='ij')
    plaquette_qubits = np.stack([
        row * d + col, row * d + col + 1, # Top-left/right
        (row + 1) * d + col, (row + 1) * d + col + 1# Bottom-left/right
     ], axis=-1)

    # Assign X and Z stabilizers in a checkerboard pattern
    is_x_plaquette = (row + col) % 2 == 1
    x_stabilizers = plaquette_qubits[is_x_plaquette].tolist()
    z_stabilizers = plaquette_qubits[~is_x_plaquette].tolist()

    # Top boundary X-stabilizers sit on the even indices
    top_edges = np.stack([grid_idxs, grid_idxs + 1], axis=-1)
    x_stabilizers += top_edges[(grid_idxs % 2) == 0].tolist()

    # Bottom boundary X-stabilizers sit on the odd indices
    bottom_edges = np.stack([(d-1)*d + grid_idxs, (d-1)*d + grid_idxs + 1], axis=-1)
    x_stabilizers += bottom_edges[(grid_idxs % 2) != 0].tolist()

    # Left boundary Z-stabilizers sit on the odd indices
    left_edges = np.stack([grid_idxs * d, (grid_idxs + 1) * d], axis=-1)
    z_stabilizers += left_edges[(grid_idxs % 2 != 0)].tolist()

    # Right boundary Z-stabilizers sit on the even indices
    right_edges = np.stack([grid_idxs*d + (d-1), grid_idxs*d + 2*d - 1], axis=-1)
    z_stabilizers += right_edges[(grid_idxs % 2 == 0)].tolist()

    # Logical X runs top-to-bottom along the left-most column
    # Logical Z runs left-to-right along the top-most row
    logical_x, logical_z = list(range(0, d * d, d)), list(range(d))

    return x_stabilizers, z_stabilizers, logical_x, logical_z

dist = 3
x_stabs, z_stabs, log_x, log_z = rotated_surface_code(dist)
print(f"X stabilizers: {x_stabs}")
print(f"Z stabilizers: {z_stabs}")

X stabilizers: [[1, 2, 4, 5], [3, 4, 6, 7], [0, 1], [7, 8]]
Z stabilizers: [[0, 1, 3, 4], [4, 5, 7, 8], [3, 6], [2, 5]]

Let’s verify our construction for a distance \(d=3\) code, for which we expect 4 \(X\)-type stabilizers, 4 \(Z\)-type stabilizers, and 1 logical qubit encoded in 9 data qubits (\(9 - 8 = 1\)).

nx, nz = len(x_stabs), len(z_stabs)
print(f"Number of X / Z stabilizers: {nx} | {nz}")
print(f"Logical X / Z operator: {log_x} | {log_z}")

n_qubits = max(map(max, (x_stabs + z_stabs))) + 1
print(f"\nNumber of data qubits: {n_qubits}")
print(f"Number of logical qubits: {n_qubits - nx - nz}")

Number of X / Z stabilizers: 4 | 4
Logical X / Z operator: [0, 3, 6] | [0, 1, 2]

Number of data qubits: 9
Number of logical qubits: 1

The Pseudo-Threshold

Before sweeping over many code distances to locate the true threshold, it is instructive to first ask a simpler question: Is this code worth using at all for my hardware? The answer is given by the pseudo-threshold, the leftmost crossing point in the schematic we saw earlier. For a single code of distance \(d\), it is referred to as \(p_{\text{pseudo}}^{(d)}\), the physical error rate at which the encoded logical error rate (LER) equals the unencoded physical error rate.

Below the pseudo-threshold (green region), the code actively suppresses errors, i.e., the LER sits beneath the unencoded line, and encoding is immediately beneficial. Above it but below $p_text{th}$ (the amber region), the $d=3$ code becomes a net liability, as it introduces more overhead than it corrects. This makes $p_text{pseudo}^{(d)}$ the break-even point for a specific code distance \(d\) and gives the lower bound on the true threshold. Therefore, encoding is only worthwhile if the hardware’s physical error rate is below its pseudo-threshold. Moreover, if \(p_{\text{pseudo}}^{(d)}\) decreases with increasing \(d\), we can assess that the code is scalable, even before computing the more expensive asymptotic threshold.

For the rotated surface code, we compute the pseudo-threshold by evaluating the logical error rate of the minimum-distance code (\(d=3\)) and comparing it against the raw physical noise level. To do so, we need syndrome extraction circuits that measure all stabilizers and return a syndrome for the decoder. For simplicity, these circuits assume instantaneous, noiseless syndrome extraction. In real hardware, each stabilizer measurement requires typically 4 CNOTs followed by a noisy readout, each introducing additional error sources. This places our simulation in the code-capacity regime, which yields a higher threshold than what is achievable in practice, a key difference we will quantify at the end of the Simulating The Threshold section.

The syndrome_extraction function below uses the stabilizers and logical operators from the previous section to build these circuits, which are then executed on the default.clifford device with depolarizing noise at a specified number of shots.

import pennylane as qp

def syndrome_extraction(stabilizers, logical_ops, num_wires, noise_param, n_shots):
    """Extract the syndrome from the stabilizers and logical operators."""
    x_stabs, z_stabs = stabilizers
    x_lg_op, z_lg_op = logical_ops

    # Build the measurement operators for the X and Z stabilizers
    z_stab_ops = [qp.prod(*[qp.Z(q) for q in s]) for s in z_stabs]
    x_stab_ops = [qp.prod(*[qp.X(q) for q in s]) for s in x_stabs]
    z_logic_op = qp.prod(*[qp.Z(q) for q in z_lg_op])
    x_logic_op = qp.prod(*[qp.X(q) for q in x_lg_op])
    z_meas_ops = z_stab_ops + [z_logic_op]
    x_meas_ops = x_stab_ops + [x_logic_op]

    # Build the syndrome extraction circuit
    @qp.set_shots(n_shots)
    @qp.qnode(qp.device("default.clifford", wires=num_wires))
    def syndrome_circuit(meas_ops, x_basis=False):
        for q in range(num_wires):
            if x_basis:
                qp.H(wires=q)
            qp.DepolarizingChannel(noise_param, wires=q)
        return [qp.sample(op) for op in meas_ops]

    # Circuit 1: |+...+⟩ -> noise -> measure X-stabs + X-logical  (Z-error syndrome)
    # Circuit 2: |0...0⟩ -> noise -> measure Z-stabs + Z-logical  (X-error syndrome)
    z_stab_res = np.column_stack(syndrome_circuit(x_meas_ops, x_basis=True))
    x_stab_res = np.column_stack(syndrome_circuit(z_meas_ops, x_basis=False))
    return z_stab_res, x_stab_res

The results from the above syndrome extraction circuits are then decoded using the minimum weight perfect matching (MWPM) decoding algorithm from the PyMatching [4] library, as shown below, to give the corrected syndromes \(\vec{c}\).

from pymatching import Matching

def syndrome_decoding(stabilizers, syndrome_results, num_wires, noise_param):
    """Decode the syndrome using PyMatching and compute corrections."""
    x_stabs, z_stabs = stabilizers
    nx, nz = len(x_stabs), len(z_stabs)

    # Z-error and X-error syndromes from the results
    z_stab_res, x_stab_res = syndrome_results
    z_syndrome = ((1 - z_stab_res[:, :nx]) // 2).astype(np.uint8)
    x_syndrome = ((1 - x_stab_res[:, :nz]) // 2).astype(np.uint8)

    # Build the parity check matrices for the X and Z stabilizers
    H_x, H_z = np.zeros((nx, num_wires)), np.zeros((nz, num_wires))
    for ix, sx in enumerate(x_stabs):
        H_x[ix, sx] = 1
    for iz, sz in enumerate(z_stabs):
        H_z[iz, sz] = 1

    # Decode the syndrome using PyMatching and compute corrections
    q_eff = 2 * noise_param / 3
    wt = np.log((1 - q_eff) / q_eff) if 0 < q_eff < 1 else 1.0
    z_corr = Matching(H_x, weights=wt).decode_batch(z_syndrome)
    x_corr = Matching(H_z, weights=wt).decode_batch(x_syndrome)
    return z_corr, x_corr

So overall, we model the surface code such that all the qubits independently suffer a depolarizing error with probability (noise_param). We then compute the \(Z\)/ \(X\)-stabilizer syndromes from the \(X\)/\(Z\)-stabilizer measurements, which are then decoded using the MWPM decoder defined above [2]. We additionally measure the logical operators to obtain the values of logical observables that are used tocompute the logical errors as a final step. More specifically, the ler_eval function below does this by computing the residual \(\vec{p} = \vec{e} \oplus \vec{c} \cdot \vec{l}\). It reconstructs the \(\vec{e} \cdot \vec{l}\) from the circuit output and uses the corrected syndromes from the function above to compute the LER as \(p_{L} = 1 - (1 - p_{x}) * (1 - p_{z})\).

def ler_eval(stabilizers, logical_ops, noise_param, num_shots=10_000):
    """Evaluate LER for a given set of stabilizers and logical operators."""
    num_wires = 2 * max(len(stabilizers[0]), len(stabilizers[1])) + 1

    # Extract the syndrome measurements
    z_stab_res, x_stab_res = syndrome_extraction(
        stabilizers, logical_ops, num_wires, noise_param, num_shots
    )

    # Decode the syndrome and compute the logical corrections
    syndrome_results = (z_stab_res[:, :-1], x_stab_res[:, :-1])
    z_corr, x_corr = syndrome_decoding(
        stabilizers, syndrome_results, num_wires, noise_param
    )

    # Build the logical operators for the X and Z stabilizers
    log_x_vec, log_z_vec = np.zeros((2, num_wires), dtype=np.uint8)
    log_x_vec[logical_ops[0]], log_z_vec[logical_ops[1]] = 1, 1

    # Pauli frame tracking: compute error from the circuit output
    x_log_meas, z_log_meas = (x_stab_res[:, -1], z_stab_res[:, -1])
    lx_raw = ((1 - x_log_meas) // 2).astype(np.uint8)
    lz_raw = ((1 - z_log_meas) // 2).astype(np.uint8)
    p_lx = (lx_raw ^ (x_corr @ log_z_vec % 2)).mean()
    p_lz = (lz_raw ^ (z_corr @ log_x_vec % 2)).mean()
    return (1 - (1 - p_lx) * (1 - p_lz))

We can now evaluate the pseudo-threshold for a given set of stabilizers and logical operators by sweeping over a range of noise parameters and evaluating the logical error rate.

from matplotlib import pyplot as plt

lerror_rates = []
p_noise = np.geomspace(0.025, 0.25, 21)
for p in p_noise:
    ler = ler_eval((x_stabs, z_stabs), (log_x, log_z), p)
    lerror_rates.append(ler)

# Approximating the pseudo-threshold by linear interpolation
diff = np.array(lerror_rates) - p_noise
p_idx = np.where(diff <= 0)[0][-1]
p0, p1 = p_noise[p_idx], p_noise[p_idx + 1]
d0, d1 = diff[p_idx], diff[p_idx + 1]
pseudo_threshold = p0 - d0 * (p1 - p0) / (d1 - d0)

plt.figure(figsize=(6, 3))
plt.axvline(x=pseudo_threshold, color="black", linestyle="--", linewidth=1)
plt.text(pseudo_threshold, 0.04, r" p$_{pseudo}$="+f"{pseudo_threshold:.3f}")
plt.plot(p_noise, lerror_rates, marker="o", label=f"Surface Code (d={dist})", color="b")
plt.plot(p_noise, p_noise, linestyle="--", label="Unencoded", color="r")
plt.xlabel("Physical Error Rate (p)", fontsize=10)
plt.ylabel(r"Logical Error Rate (p$_{L}$)", fontsize=10)
plt.yscale("log")
plt.xscale("log")
plt.grid(True, which="both", ls="--", alpha=0.6)
plt.legend(fontsize=10)
plt.tight_layout()
plt.show()

The red dashed line is the baseline, i.e., the error rate you would see with no error correction at all. On the right side of the graph (high noise), the blue curve sits above the red line, meaning QEC is making things worse because the extra circuit operations introduce more noise than they correct. Moving leftward to lower physical error rates, the blue curve eventually dips below the red line. That crossing point is the pseudo-threshold for our distance \(d=3\) code.

Simulating The Threshold

While the pseudo-threshold we computed previously tells us when a specific code distance starts helping, the true asymptotic fault-tolerant threshold tells us something deeper for the entire code family, i.e., the physical error rate below which we can keep improving by increasing the code distance. Here, we sweep over distances \(d = 3, 5, 7\) and a range of depolarizing noise strengths.

def eval_threshold(distances, p_noise, num_shots):
    """Evaluates the threshold for a given set of distances and noise levels."""
    results = {d: [] for d in distances}
    for dist in distances:
        x_stabs, z_stabs, log_x, log_z = rotated_surface_code(dist)
        for p in p_noise:
            ler = ler_eval((x_stabs, z_stabs), (log_x, log_z), p, num_shots)
            results[dist].append(ler)
    return results

distances = [3, 5, 7]
p_noise = np.geomspace(0.036, 0.36, 15)
results = eval_threshold(distances, p_noise, num_shots=20_000)

We visualize the results on a log-log plot. Below the threshold, errors are suppressed exponentially with increasing distance, so the curves fan out dramatically on a logarithmic scale.

plt.figure(figsize=(6, 4))

ler_vals = []
for d in distances:
    result = np.array(results[d])
    plt.plot(p_noise, result, label=f"Distance {d}", marker="o")
    ler_vals.append(result)

# Approximation of the threshold by linear interpolation
diff = np.diff(np.array(ler_vals), axis=0)
idxs = np.argmax(np.diff(np.sign(diff), axis=1) != 0, axis=1)[0]
p0, p1 = p_noise[idxs], p_noise[idxs + 1]
d0, d1 = diff[np.arange(len(diff)), idxs], diff[np.arange(len(diff)), idxs + 1]
p_th  = np.mean(p0 - d0 * (p1 - p0) / (d1 - d0))
plt.axvline(x=p_th, color="black", linestyle="--", linewidth=1)
plt.text(p_th, 0.04, r" p$_{th}$="+f"{p_th:.3f}")

plt.xlabel("Physical Error Rate (p)", fontsize=10)
plt.ylabel("Logical Error Rate (p$_{L}$)", fontsize=10)
plt.yscale("log")
plt.xscale("log")
plt.grid(True, which="both", ls="--", alpha=0.6)
plt.legend(fontsize=10)
plt.tight_layout()
plt.show()

The curves for different distances cross at a single point, the threshold. To the right of the crossing (high noise), larger codes perform worse, i.e., they have more qubits for errors to strike, but cannot correct them all. To the left (low noise), larger codes perform better, and the improvement is exponential with distance.

Note that this code-capacity threshold (~15%) is considerably higher than the circuit-level threshold (~0.8%) reported in the pseudo-threshold section. This gap arises directly from our code-capacity assumption of instantaneous, perfect syndrome extraction; the real hardware noise would, in fact, drive the threshold lower. Nevertheless, the qualitative picture remains the same: logical error rate curves for increasing code distances cross at a single threshold point.

Conclusion

In this tutorial, we demonstrated the threshold theorem in practice by simulating both the pseudo-threshold and the asymptotic threshold for the rotated surface codes. We saw that below the threshold point $p_{th}$, increasing the code distance $d$ leads to an exponential suppression of logical errors. Visually, this qualitative behavior gets captured by the curves for different code distances crossing at a single distinct point $p_{th}$, which is the hallmark of the threshold theorem. For topological codes like the surface code, this suppression follows the power-law relationship $p_L propto C cdot (p / p_{th})^{(d + 1)/2}$, where $C > 0$ is a constant. Since the base is $p / p_{th} < 1$, when operating below the threshold point $p < p_{th}$, the logical error rate $p_L$ drops exponentially as $d$ increases. Therefore, by simply engineering a larger code, we can suppress the logical error rate to arbitrarily low levels.

Transitioning from the above mathematical scaling to physical hardware depends heavily on how we define and reach that threshold in a laboratory setting. While in our simulations, we targeted the code-capacity threshold, which represents a theoretical upper bound by assuming perfect, instantaneous syndrome extraction. In reality, syndromes are extracted using noisy multi-qubit gates and measurements, which would push the circuit-level threshold lower. However, the qualitative picture is still preserved, guaranteeing that we are fighting a winnable battle despite the significant engineering challenges that remain, notably in scaling up the number of physical qubits and executing fast, efficient logical operations. This means that by engineering hardware that keeps physical error rates below the threshold, we can unlock the path to arbitrarily complex, reliable quantum computations.

References