Source-linked AI summary
Exponentially Accelerated Sampling of Pauli Strings for Nonstabilizerness
Zhenyu Xiao, Shinsei Ryu
TL;DR
The paper presents a FWHT-based procedure for computing stabilizer Rényi entropy and nullity. For two-qubit Haar-random-gate doping, magic increases with Clifford scrambling, reaches a higher saturation value than in the T-gate case, and approaches saturation faster.
Problem
The paper addresses computation of stabilizer Rényi entropy and nullity for an N-qubit wavefunction.
Method
The method uses fast Walsh–Hadamard transforms to sample Pauli strings and compute stabilizer Rényi entropy and nullity.
Results
At large Clifford depth, two-qubit Haar-random gates reach E[M2]/N ≈0.875, exceeding the T-gate value E[M2]/N ≈0.415, while reaching ΔM2 = 10^-5 at N = 20 and NC = 4 with around 40 gates versus around 120 T gates.
Takeaways & Limitations
Two-qubit Haar-random gates generate magic more efficiently than T gates, though they have higher experimental cost.
Abstract
from arXiv · showhide
Quantum magic, quantified by nonstabilizerness, measures departures from stabilizer structure and underlies potential quantum speedups. We introduce an efficient classical framework for computing stabilizer Rényi entropies and stabilizer nullity of generic $N$-qubit wavefunctions. The method combines the fast Walsh-Hadamard transform with an exact partition of Pauli operators, reducing the average cost per sampled Pauli string from $\mathcal{O}(2^N)$ to $\mathcal{O}(N)$. We further develop a Monte Carlo estimator with Clifford preconditioning and find that the required number of samples shows no visible growth with $N$ in our benchmarks. Applying the method to $T$-doped random Clifford circuits, we identify the scrambling ratio $η$ (Clifford gates per $T$ gate) as the key parameter governing magic growth. Each $T$ gate approaches its dilute-limit nonstabilizerness power with only modest Clifford scrambling. Our approach enables quantitative studies of magic in highly entangled states and long-time nonequilibrium dynamics.
Pseudocode for FWHT-based Pauli sampling
Algorithm 1 computes stabilizer Rényi entropy and nullity by partitioning Pauli strings into x-families and processing each family with a Walsh–Hadamard transform.
- The algorithm takes an N-qubit wavefunction and tolerance ϵ, then outputs stabilizer Rényi entropy Mα and nullity ν.
- For every x ∈ F_2^N, it defines f_x(b) from pairs of wavefunction amplitudes separated by bitwise XOR.
- The pseudocode loops over all x-families, applying the same family-wise calculation needed to assemble the global estimators.
- The algorithm accumulates each family’s partial moment mα;x and nullity count ν_x before returning Mα and ν.
Doping with two-qubit Haar-random gates
Two-qubit Haar-random gates generate magic more efficiently than single-qubit T gates, with Clifford scrambling controlling saturation and growth.
- The protocol interleaves NC layers of random brick-wall Clifford gates with N/2 independent Haar-random two-qubit gates on neighboring pairs.
- E[M2/N] saturates at NC ≳15, reaching approximately 0.875 at large NC.
- Without Clifford scrambling, E[M2] = N log2(7/4)/2 ≈0.403N, comparable to the T-gate protocol’s approximately 0.415N saturation.
- The gap ΔM2 decays exponentially with the total number of injected two-qubit Haar gates, faster than for T gates.
- At N = 20 and NC = 4, reaching ΔM2 = 10^-5 requires around 40 two-qubit Haar gates versus around 120 T gates.
A. Analytic fluctuation ratio for the product magic state
For an unscrambled product magic state, the fluctuation ratio grows exponentially because different x-families have sharply different Pauli-support profiles.
- The analysis considers |ψm⟩ = |T⟩^⊗N/2 ⊗ |0⟩^⊗N/2, assuming even N without changing the large-N scaling.
- The Pauli spectrum factorizes over the |T⟩ and |0⟩ sectors after writing x = (x_A,x_B).
- Families with x_B ≠ 0 have m2;x = 0, while families with x_B = 0 depend on the Hamming weight of x_A.
- The first two moments over x ≠ 0 are evaluated analytically, and the large-N limit is governed by the leading term under the square root.
- The resulting fluctuation ratio produces the exponential sample scaling quoted before Clifford preconditioning.
B. Brick-wall circuits and Ising ground states
Supplementary benchmarks test Clifford-preconditioned Monte Carlo sampling on brick-wall circuit states and mixed-field Ising ground states, finding bounded fluctuations across the studied sizes.
- The benchmarks include Haar-random brick-wall circuits of depths dH = 1, 2, 4, 6 and open-chain mixed-field Ising-model ground states.
- Figure S1 compares normalized fluctuation r2 before and after preconditioning at depth NC = 2N as system size varies.
- After preconditioning, brick-wall states have r2 below 3, so N ∼10^4 samples suffice for accuracy ϵ ∼0.03.
- For N ≥16, preconditioned brick-wall states show no clear increasing trend in r2, and deeper circuits with dH ≥2 reduce r2 at larger sizes.
- For mixed-field Ising ground states, r2 remains below 1 for N = 10, 12, …, 22 after preconditioning.
- Across all three state classes, r2 is close to its saturated value by NC ≈2N, supporting NC = 2N as a practical default.
III. DIRECT PAULI-STRING SAMPLING WITH METROPOLIS–HASTINGS UPDATES
The Metropolis–Hastings sampler targets Pauli strings according to squared expectation values, but each proposal costs O(2^N), and low acceptance can produce strongly correlated samples.
- The target distribution is π(P) := |⟨ψ|P|ψ⟩|2/d over Pauli strings.Samples from its stationary Markov chain estimate the Pauli expansion sum and then the stabilizer Rényi entropy.
- Each update changes one local Pauli operator to one of the other three choices using a symmetric proposal.The acceptance probability therefore reduces to the Metropolis–Hastings ratio for the target distribution.
- O(2^N) operations are required for each proposed move because it evaluates a full-state Pauli expectation value.The MC+FWHT scheme instead has amortized cost O(N) per Pauli string.
- Low acceptance makes consecutive samples highly correlated, reducing the effective number of independent samples below the recorded sample count.The chain uses N_burn = 10^4 burn-in steps before recording samples.
- For C|ψ_m⟩ at N = 16, the MH acceptance ratio is ≈0.03 and estimates fluctuate substantially more despite 10× more CPU time.This state contains T magic states scrambled by a random Clifford circuit of depth N_C = 2N.
A. Single-shot protocol
The single-shot protocol prepares an entangled stabilizer state, applies T gates, and analyzes the resulting nonstabilizerness using a replica-based stabilizer linear entropy and Clifford averaging.
- Protocol: The protocol starts from single-qubit Clifford states, applies N_C layers of random two-qubit Clifford gates, and then applies T to every qubit.The Clifford circuit creates the entangled stabilizer state |ψ(N_C)⟩ before T injection.
- Protocol: At N_C = 0, the T-doped state remains a product state, so M_2 is additive across qubits.This provides the no-scrambling baseline for the single-shot analysis.
- No-scrambling baseline: For each site, T|ϕ_i⟩ has M_2 = 0 with probability 1/3 and M_2 = log_2(4/3) otherwise.The distinction is whether the initial Clifford state is a Z eigenstate.
- Analytical framework: At sufficiently large N_C, the random brick-wall circuit is approximated by a uniform average over the N-qubit Clifford group.This replacement is the stated analytical assumption behind the large-scrambling calculation.
- Analytical framework: The four-replica Clifford-invariant subspace is two-dimensional, allowing the averaged quantity to be represented using Π_sym and QΠ_sym.This reduction supports a closed-form single-shot calculation.
B. Multi-round protocol
The multi-round protocol alternates independent Clifford scrambling with T-gate layers and reduces the averaged dynamics to a scalar recursion whose asymptotic behavior determines convergence toward the Haar value.
- Protocol: The m-round unitary alternates T layers τ with independent random Clifford circuits C_1,…,C_m, starting from |Φ⟩.The total number of injected T gates is N_T = mN.
- Recursion: After each Clifford twirl and T layer, the four-copy state evolves as ρ_m+1 = τ^⊗4 T(ρ_m) τ†⊗4.The twirling channel preserves the relevant permutation-symmetric Clifford-invariant structure.
- Recursion: The averaged dynamics close into a one-round scalar recursion because the relevant invariant subspace is spanned by Π_sym and QΠ_sym.The coefficients are fixed by trace preservation and the Q moment.
- Recursion: The recursion uses the initial condition x_0 = 1 because the initial product state |Φ⟩ is stabilizer with M_lin = 0.This condition is inserted into the recurrence to obtain the closed-form expression for x_m.
- Asymptotics: For small fluctuations, the expected second-order stabilizer Rényi entropy is approximated from the averaged stabilizer linear entropy.The resulting expression depends on the recursion parameter λ_N and the number of rounds m.
- Asymptotics: Because λ_N scales as (3/4)N and N_T = mN, the gap to the Haar value decays exponentially with total T-gate count.The asymptotic decay rate per T gate follows from this large-N scaling.
- Numerical validation: Figure S4 tests the predicted decay using two schedules at η = 2 and three system sizes at η = 5.Each data point averages 80 random realizations, with exact evaluation for N ≤ 16 and Monte Carlo estimates for larger systems.
V. CONVERGENCE OF M2 TO THE HAAR VALUE AT ADDITIONAL SCRAMBLING RATIOS AND SYSTEM SIZES
Additional simulations show that the late-time convergence rate depends on the scrambling ratio η rather than the temporal distribution of T gates, and is already near saturation at η = 5 across tested sizes.
- Schedule dependence: At η = 2 and N = 20, uniform “5T/1L” and bursty “20T/4L” schedules give s = 0.178 ± 0.007 and s = 0.170 ± 0.014.The agreement supports the conclusion that the growth rate is insensitive to temporal T-gate distribution at fixed η.
- System-size dependence: At η = 5, the fitted slopes are s_N=16 = 0.295±0.020, s_N=20 = 0.284±0.010, and s_N=24 = 0.290±0.010.These values are close to the analytical asymptote s_∞ = ln(4/3) ≈ 0.288.
VI. NONSTABILIZERNESS POWER OF EXTENSIVE HAAR-RANDOM TWO-QUBIT GATES
The section analyzes a layer of independent Haar-random two-qubit gates applied to a Clifford-scrambled stabilizer state and derives its averaged nonstabilizerness. The resulting extensive layer generates extensive nonstabilizerness and approaches the reported saturation values for studied system sizes.
- Setup: The final T layer is replaced by N/2 independent Haar-random two-qubit gates acting on neighboring qubit pairs.Each gate acts on (2j − 1, 2j) and is sampled independently from the Haar measure on U(4).
- Averaging method: The calculation averages over the product of independent two-qubit Haar unitaries using a two-qubit Haar-twirl decomposition.The derivation groups permutations by cycle type after factorizing the final layer across qubit pairs.
- Approximation: A fluctuation-smallness assumption permits approximating the relevant nonstabilizerness quantity from its averaged linear measure.The section then evaluates the resulting expression for finite system sizes.
- Result: One extensive Haar-random two-qubit layer generates extensive nonstabilizerness on top of a Clifford-scrambled stabilizer state.The exact expression is used to evaluate −log2(1 − E_U E_C[Mlin]) for the system sizes studied in the main text.
- Result: The reported values are approximately 0.799 for N = 12, 0.846 for N = 16, and 0.875 for N = 20.These values agree well with the saturation values shown in the main text.
VII. STABILIZER R´ENYI ENTROPY OF MIXED STATES
This section extends stabilizer Rényi entropy calculations to a mixed-state subsystem by restricting attention to Pauli correlators supported on that subsystem. The method uses FWHTs to organize these correlators, while brute-force enumeration has cost O(2^N+2N_A).
- Definition: The stabilizer Rényi entropy is defined for the reduced density matrix ρ_A of a subsystem A containing N_A qubits.The subsystem dimension is d_A = 2^N_A.
- Subsystem restriction: Evaluating the stabilizer Rényi quantities for ρ_A requires only Pauli correlators supported on A.The subsystem is taken to be the first N_A qubits, with the remaining qubits forming its complement.
- FWHT evaluation: For each fixed x_A, an N_A-qubit FWHT obtains the full set of correlators indexed by z_A.The transform acts on a length-2^N_A vector constructed from amplitudes summed over the complement configuration.
- Complexity: Computing the complement sums for all subsystem configurations costs O(2^N), followed by an FWHT cost of O(N_A 2^N_A).The stated complexity is obtained before sweeping over all x_A values.
- Comparison: Brute-force enumeration of ⟨ψ|P|ψ⟩ for every Pauli string supported on A costs O(2^N+2N_A).