Source-linked AI summary

Exponentially Accelerated Sampling of Pauli Strings for Nonstabilizerness

Zhenyu Xiao, Shinsei Ryu

arXiv:2601.00761v2quant-phcond-mat.mes-hallcond-mat.stat-mech

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.

Abstract

from arXiv · show

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.

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.

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.

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.

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.Source: §III · pp. 12–13
  • 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.Source: §III · p. 13
  • 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.Source: §III · p. 13
  • 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.Source: §III · p. 13
  • 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.Source: §III · p. 13

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.Source: A. Single-shot protocol · pp. 13–14
  • 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.Source: A. Single-shot protocol · p. 14

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.Source: A. Single-shot protocol · p. 14

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.Source: A. Single-shot protocol · p. 14
  • 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.Source: A. Single-shot protocol · pp. 14–15

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.Source: B. Multi-round protocol · p. 15

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.Source: B. Multi-round protocol · pp. 15–16
  • 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.Source: B. Multi-round protocol · p. 15
  • 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.Source: B. Multi-round protocol · p. 16

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.Source: B. Multi-round protocol · p. 16
  • 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.Source: B. Multi-round protocol · p. 16

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.Source: Figure S4 · p. 17

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 η.Source: §V · p. 17

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.Source: §V · p. 17

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).Source: §VI · p. 17

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.Source: §VI · p. 18

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.Source: §VI · p. 18

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.Source: §VI · p. 18
  • 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.Source: §VI · p. 18

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.Source: §VII · p. 19

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.Source: §VII · p. 19

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.Source: §VII · p. 19

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.Source: §VII · p. 19

Comparison

  • Brute-force enumeration of ⟨ψ|P|ψ⟩ for every Pauli string supported on A costs O(2^N+2N_A).Source: §VII · p. 19
Loading 2601.00761v2…