Randomized Benchmarking via Clifford Twirling

Applications of Scrambling, numerical solution.

Randomized Benchmarking via Clifford Twirling

---

Background

Randomized benchmarking reads the average error per gate off the way a survival probability falls as circuits get longer, sidestepping full process tomography, which scales badly and is itself corrupted by errors in preparing and measuring the qubit. A run applies a sequence of random Clifford gates and then a single recovery gate chosen to invert the whole sequence, so with perfect gates the qubit comes back exactly to where it started. Any drop in the return probability is therefore due to the noise, and averaging over the random Cliffords reshapes that noise into a standard form.

The reshaping is the Clifford twirl. Averaging an error channel over the Clifford group, which is a unitary 2-design, yields a map that commutes with every Clifford, and Schur’s lemma then forces such a map to act as a single scalar \(\alpha\) on the traceless part of any state while leaving the identity fixed, the signature of a depolarizing channel. The intuition is that the random gates scramble any structure the noise had, leaving only its overall strength.

Stacking \(m\) of these depolarizing channels multiplies the traceless part by \(\alpha^m\), so the survival probability follows \(\bar p(m) = A\alpha^m + B\), with the preparation and measurement errors swept entirely into the constants \(A\) and \(B\). One exponential fit returns \(\alpha\), and from it the average gate fidelity \(F_{\mathrm{avg}} = \alpha + (1-\alpha)/D\) and the error rate \(r = (D-1)(1-\alpha)/D\).

A clean single exponential assumes noise that is the same from gate to gate and from step to step, so correlated or coherent errors reveal themselves as deviations from that form. The calculation runs the protocol under depolarizing noise and fits \(\alpha\), then confirms that twirling an amplitude-damping channel, which is neither depolarizing nor even unital to begin with, lands on a depolarizing channel.

Motivation

Randomized benchmarking (RB) is the standard method for characterizing gate fidelity in quantum hardware. It exploits a deep structural result: averaging any CPTP error channel over the Clifford group (a unitary 2-design) twirls it into a depolarizing channel. This twirling theorem is a direct consequence of Schur’s lemma and the Weingarten calculus, the same machinery that produces barren plateaus also produces the clean exponential decay of RB.

The protocol measures only a single scalar quantity (survival probability) as a function of circuit depth, yet extracts the average gate fidelity with no need for full process tomography.

Preliminaries / Toolbox

1. Standard Randomized Benchmarking (RB): A protocol to estimate the average gate fidelity. Sequences of length \(m\) random Clifford gates \(C_i\) are followed by an exact inverse \(C_m^{-1} \cdots C_1^{-1}\) to yield an identity operation. Survival probability decays as \(A p^m + B\).

2. Twirling over the Clifford Group: By Schur’s lemma, averaging a channel \(\mathcal{E}\) over the Clifford group produces a depolarizing channel: \(\frac{1}{|\mathcal{C}|} \sum_{C \in \mathcal{C}} C^\dagger \mathcal{E}(C \rho C^\dagger) C = \mathcal{E}_{\mathrm{depol}}(\rho) = p \rho + (1-p) \frac{I}{d}\).

3. Depolarizing Parameter \(p\) vs Fidelity: The average gate fidelity is related to \(p\) (also written as \(\alpha\)) by \(\bar{F} = p + \frac{1-p}{d}\).

4. Sequence fidelity: The expected return probability is \(\mathrm{Tr}(E \cdot \mathcal{E}_{\mathrm{depol}}^m(\rho_0)) = A p^m + B\) where \(A, B\) absorb state preparation and measurement (SPAM) errors.

Exercise

(a) Prove the twirling theorem: for any CPTP map \(\Lambda\), the Clifford twirl

\[ \widetilde{\Lambda}(\rho) = \frac{1}{|\mathcal{C}|}\sum_{C \in \mathcal{C}} C^\dagger \Lambda(C \rho C^\dagger) C \]

is a depolarizing channel: \(\widetilde{\Lambda}(\rho) = \alpha\, \rho + (1-\alpha)\frac{I}{D}\).

(b) Show the depolarizing parameter is \(\alpha = \frac{D\, F_{\mathrm{avg}} - 1}{D - 1}\), where \(F_{\mathrm{avg}}\) is the average gate fidelity.

(c) Derive the RB decay: \(\bar{p}(m) = A\,\alpha^m + B\).

(d) Implement and run the full RB protocol numerically; fit \(\alpha\) and extract \(F_{\mathrm{avg}}\).

(e) Verify numerically that the Clifford twirl of the amplitude-damping channel (a non-depolarizing, non-unital channel) produces a depolarizing channel.

Solution

Part (a): The twirling theorem

Claim: If \(\mathcal{C}\) forms a unitary 2-design, then \(\widetilde{\Lambda}(\rho) = \alpha\,\rho + (1-\alpha)\frac{I}{D}\) for some \(\alpha\).

Proof: The twirled channel \(\widetilde{\Lambda}\) commutes with all Clifford unitaries: \(C \widetilde{\Lambda}(\rho) C^\dagger = \widetilde{\Lambda}(C\rho C^\dagger)\) for all \(C \in \mathcal{C}\).

By Schur’s lemma (applied to the irreducible decomposition of the adjoint representation), any CPTP map that commutes with all elements of a 2-design must act as a scalar on each irreducible sector. For qubits, the space of \(D \times D\) Hermitian matrices decomposes into:

• The trivial sector: \(\{c \cdot I \mid c \in \mathbb{R}\}\) (dimension 1).
• The traceless sector: \(\{A \mid \mathrm{Tr}(A) = 0\}\) (dimension \(D^2 - 1\)).

The twirled channel must act as:

\[ \widetilde{\Lambda}(I) = I \quad \text{(trace preservation)}, \qquad \widetilde{\Lambda}(A) = \alpha\, A \quad \text{(traceless sector)}. \]

For general \(\rho = \frac{I}{D} + (\rho - \frac{I}{D})\), this gives:

\[ \widetilde{\Lambda}(\rho) = \frac{I}{D} + \alpha\left(\rho - \frac{I}{D}\right) = \alpha\,\rho + (1-\alpha)\frac{I}{D}. \quad \checkmark \]

Part (b): Connecting \(\alpha\) to average gate fidelity

The average gate fidelity of a channel \(\Lambda\) (relative to the identity) is:

\[ F_{\mathrm{avg}} = \int d\psi\, \langle\psi|\Lambda(|\psi\rangle\langle\psi|)|\psi\rangle. \]

For the depolarizing channel \(\widetilde{\Lambda}(\rho) = \alpha\rho + (1-\alpha)I/D\):

\[ F_{\mathrm{avg}} = \int d\psi\, \left[\alpha\, \langle\psi|\psi\rangle^2 + (1-\alpha)\frac{\langle\psi|\psi\rangle}{D}\right] = \alpha + \frac{1-\alpha}{D}. \]

Solving for \(\alpha\):

\[ \boxed{\alpha = \frac{D\, F_{\mathrm{avg}} - 1}{D - 1}.} \]

At \(F_{\mathrm{avg}} = 1\): \(\alpha = 1\) (identity). At \(F_{\mathrm{avg}} = 1/D\): \(\alpha = 0\) (fully depolarizing).

Part (c): The exponential decay

In the RB protocol, a sequence of \(m\) random Clifford gates is applied, each followed by noise \(\Lambda\). The effective channel for the full sequence is \(m\) concatenated twirled channels:

\[ \widetilde{\Lambda}^m(\rho) = \alpha^m \rho + (1 - \alpha^m) \frac{I}{D}. \]

(Composing depolarizing channels: each application multiplies the traceless component by \(\alpha\), giving \(\alpha^m\) after \(m\) steps.)

The survival probability after the recovery gate is:

\[ \bar{p}(m) = \langle 0| \widetilde{\Lambda}^{m+1}(|0\rangle\langle 0|) |0\rangle = \alpha^{m+1}\underbrace{\langle 0|0\rangle\langle 0|0\rangle}_{=1} + (1-\alpha^{m+1})\underbrace{\frac{1}{D}}_{} \]

\[ = \alpha^{m+1} + \frac{1 - \alpha^{m+1}}{D} = \left(1 - \frac{1}{D}\right)\alpha^{m+1} + \frac{1}{D}. \]

This is the characteristic RB form \(\bar{p}(m) = A\,\alpha^m + B\) with \(A = (1 - 1/D)\alpha\) and \(B = 1/D\). A single exponential fit extracts \(\alpha\).

---

Symbolic Verification (SymPy)

Code

import sympy as sp

d, p, m = sp.symbols('d p m', positive=True)

# Sequence fidelity: F(m) = A*p^m + B
# Average gate fidelity: F_avg = p + (1-p)/d
F_avg = p + (1 - p) / d
print(f'Average gate fidelity: F_avg = {F_avg}')

# Solve for p in terms of F_avg
p_solved = sp.solve(F_avg - sp.Symbol('F'), p)
print(f'p = {p_solved}')

# Average gate error rate
r = (1 - p) * (d - 1) / d
print(f'Error rate r = {sp.simplify(r)}')

# For d=2 (single qubit)
print(f'\nSingle qubit (d=2):')
print(f'  F_avg = {F_avg.subs(d, 2)}')
print(f'  r = {sp.simplify(r.subs(d, 2))}')

Average gate fidelity: F_avg = p + (1 - p)/d
p = [(F*d - 1)/(d - 1)]
Error rate r = -(d - 1)*(p - 1)/d

Single qubit (d=2):
  F_avg = p/2 + 1/2
  r = 1/2 - p/2

---

Parts (d) and (e): Numerical Implementation

Code

import numpy as np
from scipy.optimize import curve_fit

np.random.seed(42)

# Single-qubit Clifford group (24 elements)
H = np.array([[1,1],[1,-1]], dtype=complex) / np.sqrt(2)
S = np.array([[1,0],[0,1j]], dtype=complex)
I2 = np.eye(2, dtype=complex)

def gen_cliffords():
    result = []
    seen = set()
    queue = [I2.copy()]
    while queue and len(result) < 24:
        U = queue.pop(0)
        for e in U.flatten():
            if abs(e) > 1e-10:
                U = U * np.exp(-1j*np.angle(e))
                break
        key = tuple(np.round(U.flatten(), 6))
        if key in seen: continue
        seen.add(key); result.append(U)
        for g in [H, S]:
            queue.append(g @ U)
    return result

cliffords = gen_cliffords()
print(f"Generated {len(cliffords)} single-qubit Cliffords")

Generated 24 single-qubit Cliffords

Code

# Noise model: depolarizing channel
sx = np.array([[0,1],[1,0]], dtype=complex)
sy = np.array([[0,-1j],[1j,0]])
sz = np.array([[1,0],[0,-1]], dtype=complex)

def depolarize(rho, p):
    return (1-p)*rho + (p/3)*(sx@rho@sx + sy@rho@sy + sz@rho@sz)

def noisy_gate(U, rho, p):
    return depolarize(U @ rho @ U.conj().T, p)

# RB protocol
p_error = 0.02
seq_lengths = [1, 2, 4, 8, 16, 32, 64, 128]
psi0 = np.array([1,0], dtype=complex)
rho0 = np.outer(psi0, psi0.conj())

survivals = []
for m in seq_lengths:
    p_list = []
    for _ in range(300):
        gates = [cliffords[i] for i in np.random.randint(0, len(cliffords), m)]
        U_total = np.eye(2, dtype=complex)
        for C in gates: U_total = C @ U_total
        gates.append(U_total.conj().T)  # recovery

        rho = rho0.copy()
        for C in gates:
            rho = noisy_gate(C, rho, p_error)
        p_list.append((psi0.conj() @ rho @ psi0).real)
    survivals.append(np.mean(p_list))
    print(f"m={m:4d}: p̄ = {survivals[-1]:.4f}")

m=   1: p̄ = 0.9737
m=   2: p̄ = 0.9611
m=   4: p̄ = 0.9368
m=   8: p̄ = 0.8920
m=  16: p̄ = 0.8158
m=  32: p̄ = 0.7049
m=  64: p̄ = 0.5863
m= 128: p̄ = 0.5153

Code

# Fit exponential decay
def rb_model(m, A, alpha, B):
    return A * alpha**m + B

ms = np.array(seq_lengths, dtype=float)
popt, _ = curve_fit(rb_model, ms, survivals, p0=[0.5, 0.98, 0.5])
A_fit, alpha_fit, B_fit = popt

D = 2
F_fit = alpha_fit + (1-alpha_fit)/D
alpha_th = 1 - 4*p_error/3
F_th = 1 - 2*p_error/3

print(f"\nFit: A={A_fit:.4f}, α={alpha_fit:.4f}, B={B_fit:.4f}")
print(f"α_fit = {alpha_fit:.6f}  (theory {alpha_th:.6f})")
print(f"F_avg = {F_fit:.6f}  (theory {F_th:.6f})")

Fit: A=0.4867, α=0.9733, B=0.5000
α_fit = 0.973333  (theory 0.973333)
F_avg = 0.986667  (theory 0.986667)

Code

# Part (e): Twirling the amplitude-damping channel
gamma = 0.1
K0 = np.array([[1,0],[0,np.sqrt(1-gamma)]], dtype=complex)
K1 = np.array([[0,np.sqrt(gamma)],[0,0]], dtype=complex)

def amp_damp(rho):
    return K0 @ rho @ K0.conj().T + K1 @ rho @ K1.conj().T

def twirl(rho, channel, cliff_list):
    out = np.zeros_like(rho)
    for C in cliff_list:
        out += C.conj().T @ channel(C @ rho @ C.conj().T) @ C
    return out / len(cliff_list)

print("\nTwirling amplitude-damping channel:")
for label, rho in [("  |0⟩⟨0|", np.array([[1,0],[0,0]],dtype=complex)),
                    ("  |+⟩⟨+|", np.array([[.5,.5],[.5,.5]],dtype=complex)),
                    ("  mixed ", np.array([[.7,.2j],[-.2j,.3]],dtype=complex))]:
    out = twirl(rho, amp_damp, cliffords)
    diff_out = out - I2/2
    diff_in = rho - I2/2
    mask = np.abs(diff_in) > 1e-10
    if mask.any():
        p_vals = diff_out[mask] / diff_in[mask]
        p_est = np.mean(p_vals.real)
        print(f"{label}: α = {p_est:.4f} (depolarizing confirmed)")

print("\nClifford twirl: any channel → depolarizing.")

Twirling amplitude-damping channel:
  |0⟩⟨0|: α = 0.9325 (depolarizing confirmed)
  |+⟩⟨+|: α = 0.9325 (depolarizing confirmed)
  mixed : α = 0.9325 (depolarizing confirmed)

Clifford twirl: any channel → depolarizing.

---

Visual Pedagogical Proof

Code

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

# Decay curve from the actual simulated RB run (cells above): data and fit.
m_grid = np.linspace(min(seq_lengths), max(seq_lengths), 200)

plt.figure(figsize=(7, 4))
plt.scatter(seq_lengths, survivals, color='royalblue', zorder=3,
            label='RB data (depolarizing noise)')
plt.plot(m_grid, rb_model(m_grid, *popt), color='crimson', lw=2,
         label=fr'fit $A\,\alpha^m+B$, $\alpha={alpha_fit:.3f}$')
plt.axhline(B_fit, color='gray', ls=':', label=r'asymptote $B=1/D$')
plt.xlabel('sequence length $m$')
plt.ylabel(r'survival probability $\bar p(m)$')
plt.title(r'Randomized benchmarking decay $F(m)=A\,p^m+B$')
plt.legend()
plt.grid(True, ls='--', alpha=0.4)
plt.tight_layout()
plt.show()

Takeaway

The Clifford twirl reduces an arbitrary CPTP error to a single depolarizing parameter \(\alpha\), so the survival probability follows \(F(m)=A\,\alpha^m+B\) and one exponential fit returns \(\alpha\) and the average fidelity. State-preparation and measurement errors fall entirely into \(A\) and \(B\), leaving the fidelity read from the decay rate unbiased.

Code

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Simulated RB decay data
m_lengths = np.linspace(1, 50, 20, dtype=int)
p_true = 0.92
exact_decay = 0.5 * p_true**m_lengths + 0.5
np.random.seed(42)
noisy_decay = exact_decay + np.random.normal(0, 0.02, size=len(m_lengths))

def decay_func(m, A, p, B):
    return A * p**m + B

popt, _ = curve_fit(decay_func, m_lengths, noisy_decay, p0=[0.5, 0.9, 0.5])

plt.rcParams['figure.dpi'] = 150
plt.figure(figsize=(7, 4))
plt.scatter(m_lengths, noisy_decay, color='royalblue', label='Simulated RB Data')
plt.plot(m_lengths, decay_func(m_lengths, *popt), 'crimson', linewidth=2,
         label=f'Exponential Fit (p={popt[1]:.3f})')
plt.plot(m_lengths, exact_decay, 'k--', alpha=0.5, label='Ideal Decay')

plt.xlabel("Sequence Length $m$")
plt.ylabel("Sequence Fidelity")
plt.title("Randomized Benchmarking Decay Curve")
plt.legend()
plt.grid(True, ls="--", alpha=0.5)
plt.tight_layout()
plt.show()

Connections

The Clifford twirl that reduces complicated gate noise to a single depolarizing parameter is the second Haar moment again, the same average that builds the classical shadows channel and gives the Page formula. The exponential decay you fit here is the experimental shadow of that one calculation, which is why benchmarking and shadow tomography rest on the same theory.

References

• Emerson et al., Scalable noise estimation with random unitary operators, Journal of Optics B: Quantum and Semiclassical Optics (2005).
• Magesan et al., Scalable and robust randomized benchmarking of quantum processes, Physical Review Letters (2011).
• Wallman and Flammia, Randomized benchmarking with confidence, IOP Publishing (2014). doi
• Magesan et al., Characterizing quantum gates via randomized benchmarking, Physical Review A (2012).
• Płodzień, Information Scrambling, Quantum Chaos, and Haar-Random States, lecture notes (2025). arXiv:2511.14397