Sign-problem free quantum stochastic series expansion algorithm on a quantum computer

The study of strongly correlated many-body systems often requires numerical methods because standard analytical techniques fail. Quantum Monte Carlo (QMC) methods, including the Stochastic Series Expansion (SSE), are powerful for simulating a broad range of Hamiltonians at all temperatures and dimensions. With recent advances in quantum computing, there is promise of quantum algorithms achieving speedups by directly exploiting superposition. This work proposes a quantum version of the SSE algorithm (quantum SSE) and contrasts it with the classical SSE approach. The central question is whether quantum hardware can remove key limitations of classical SSE—namely the no-branching requirement that leads to sign problems—and offer complexity advantages. The authors argue that lifting the no-branching condition on a quantum computer enables measuring more general observables and guarantees nonnegative configuration weights sampled from quantum circuit measurement probabilities that scale linearly with system size. They position quantum SSE against methods like Quantum Metropolis Sampling, which avoids sign problems but relies on deep phase estimation circuits and Trotterization, potentially introducing systematic errors that SSE avoids.

The paper reviews several QMC approaches and related methods: classical SSE, which samples the Taylor expansion of e^{-βH} and yields numerically exact results up to statistical error; world line QMC methods and continuous-time variants; and DMRG for low-dimensional systems. It highlights limitations due to the sign problem in classical QMC and the no-branching condition necessary for efficient classical SSE. Quantum Metropolis Sampling avoids sign problems via phase estimation but introduces Trotterization and deep circuits. Continuous-time world line methods avoid Trotterization but yield more complex operator structures than SSE, making SSE more suitable for near-term quantum implementations. Prior complexity results establish that some instances (e.g., certain lattice configurations for Ising-type models) are NP-hard and that k-local ground state energy estimation is QMA-complete, framing expectations for what quantum SSE can and cannot achieve.

The authors construct a quantum SSE algorithm that estimates configuration weights for SSE sampling using quantum circuits. They first present a special case where Hamiltonian terms H_b mutually commute (e.g., decompositions into I and σ^x for a nearest-neighbor σ^xσ^x model). By shifting the Hamiltonian to ensure positive semidefiniteness, they define operator blocks H^b = h0 I + h0 σ^x, guaranteeing nonnegative contributions. A controlled-unitary U acting on system and ancilla qubits is designed so that the expectation value ⟨ψ_n|U|ψ_n⟩/2^n equals ⟨α|H_{b1}…H_{bn}|α⟩/2^n, where |ψ_n⟩ encodes the system state and ancilla superpositions. Measuring the probability of all-zero outcomes over repeated circuit preparations provides an unbiased estimator of the relative weight q(n,b,α), with sampling variance scaling as 1/t; amplitude estimation can further improve scaling to 1/t^2 with respect to the number of unitary calls. These weights feed a Metropolis-Hastings sampler over configurations C=(n,b,α), with acceptance probabilities determined by ratios of estimated q values; updates to n, b, and α are performed separately. For general Hamiltonians where H_{b1}…H_{bn} may be non-Hermitian and sign issues arise, they only require the real part of ⟨α|H_{b1}…H_{bn}|α⟩. They add a sufficiently large constant shift to each term to ensure nonnegativity of the real part. Specifically, given cutoff M ≥ n, they define H_b := |h_b|[sgn(h_b)+2M I_A], representing an unequal superposition of 2^M unitaries. With an extended state |ψ_n⟩ = |α_a⟩ ⊗ |β_b⟩ ⊗ ... ⊗ |β_b⟩ ⊗ |c⟩, where |β_b⟩ = sqrt(2M/(2M+1))|0⟩ + 1/sqrt(2M+1)|1⟩ and |+_c⟩, they define a controlled unitary V_{A,B,C} such that ⟨ψ_n|V_{A,B,C}|ψ_n⟩ yields the sum of the string and its Hermitian-conjugate contributions normalized by (2M+1)^n. The estimator q(n,b,a) = |Re⟨α|H_b…H_b|α⟩|^2/(2M+1)^{2n} is then sampled via circuit measurements, ensuring nonnegative weights. The Metropolis procedure proceeds with acceptance ratios formed from these q estimates. They implement a proof-of-principle simulation for a 1D antiferromagnetic spin-1/2 chain with periodic boundary conditions using Qiskit. After adding identity components to make individual bond operators positive semidefinite, the effective terms are H_i = 1 − σ_i^x σ_{i+1}^x (J=1). A controlled unitary U_B maps ancilla-system states so that ⟨α|U_B…U_B|α⟩ = 2^n ⟨α|H_0…H_0|α⟩. The circuit uses ancilla preparation in |−⟩, prepares system basis states via H and non-Clifford T gates, applies layers of controlled-Pauli operations to realize the operator string, and measures in the computational basis; the probability of all-zero outcomes yields the square of the desired expectation value. Metropolis sampling initializes with an operator string, runs 10^4 burn-in steps, and then accumulates ⟨n⟩ to compute the energy via E = ⟨n⟩/β + N at β=1 for N=3,4,5.

- Quantum SSE lifts the classical SSE no-branching constraint, enabling sampling of configuration weights as nonnegative probabilities derived directly from quantum circuit measurements, thereby avoiding the sign problem. - In sign-problem instances where classical SSE has exponential cost due to negative weights, a single Monte Carlo update in quantum SSE requires circuits whose size scales linearly with system size N (O(N)), with depth in the QNC class (O(log N)) for relevant controlled-Pauli layers and efficiently preparable basis states. - Even when classical SSE is efficient, quantum SSE permits measurement of more general observables by choosing bases that diagonalize the observable and can be prepared efficiently on quantum hardware, including non-Clifford state preparations that are classically hard to simulate. - Numerical simulations with Qiskit for a 1D antiferromagnetic spin-1/2 chain (N=3,4,5; β=1) show the mean energy from quantum SSE converges to exact diagonalization results, validating the algorithm. - The method avoids Trotterization and thus systematic Trotter errors; errors are purely statistical from sampling and can be improved via amplitude estimation.

The findings address the central limitation of classical SSE—the sign problem stemming from the no-branching condition—by leveraging quantum superposition to directly estimate nonnegative configuration weights without tracking an exponentially large superposition. This removes a key bottleneck and enables efficient Monte Carlo updates with O(N) gate counts and polylogarithmic depth under suitable parallelization of controlled-Pauli layers. The ability to select bases that diagonalize desired observables without being constrained by no-branching expands the set of measurable quantities (e.g., state projectors for overlap with thermal states). Compared to classical SSE, which faces exponential overhead in sign-problem regimes and is restricted in observable choice, quantum SSE provides both complexity and capability advantages. The demonstration on small spin chains confirms convergence to exact thermodynamics, supporting the practical validity of the approach. Nonetheless, while quantum SSE accelerates configuration weight estimation, it does not guarantee efficient solutions to QMA-complete ground-state problems or exponential convergence of Monte Carlo sampling overall.

The paper introduces a quantum implementation of the SSE Monte Carlo algorithm that guarantees nonnegative configuration weights and avoids the sign problem by lifting the no-branching constraint inherent in classical SSE. The quantum approach estimates configuration weights via shallow quantum circuits with O(N) gate counts (and QNC-class depth) and uses these estimates in a Metropolis sampler. Simulations for small 1D antiferromagnetic chains show convergence to exact energies. The method enables measurement of more general observables than classical SSE and eliminates Trotter errors. While this removes a major classical bottleneck and yields potential exponential advantages in sign-problem settings, it does not imply efficient solutions to all many-body problems (e.g., QMA-complete cases). Future directions include optimizing the minimal constant shifts required to ensure positivity for specific models, extending implementations to broader classes of Hamiltonians and observables, and deploying on NISQ devices where shallow, parallel circuits are advantageous.

- The general construction may require adding large constant shifts to Hamiltonian terms to ensure nonnegative real contributions, which could be suboptimal; model-specific optimization is needed. - Although configuration weight estimation is polynomial-time on quantum hardware, this does not guarantee rapid overall Monte Carlo convergence or polynomial-time solutions to QMA-complete problems. - Practical performance depends on the ability to efficiently prepare basis states |α⟩ and implement controlled-Pauli layers with low error; noise on NISQ devices and finite sampling may impact accuracy. - The demonstrated results are on small systems; scaling performance on larger, noisy hardware remains to be validated experimentally.