Hamiltonian dynamics on digital quantum computers without discretization error

The paper addresses the problem of simulating real-time Hamiltonian dynamics on digital quantum computers without incurring discretization (Trotter) errors, while maintaining finite—indeed bounded on average—circuit depth independent of target precision. Conventional product formulas (Trotterization) approximate continuous dynamics with stepwise unitaries, requiring vanishing step sizes (and hence unbounded depth) to eliminate errors. Randomized methods like qDRIFT improve certain scalings for non-sparse Hamiltonians but still exhibit discretization error unless depth grows. Advanced techniques (LCU, quantum signal processing, quantum walks) reduce asymptotic error but entail large resource overheads. The authors seek an algorithm that provides exact, unbiased estimation of time-evolved expectation values with finite average circuit depth, applicable to general and time-dependent Hamiltonians, and practical for near-term noisy devices.

Prior work in Hamiltonian simulation spans: (a) product formulas (Trotter–Suzuki) and their refinements, including error analyses and randomized variants (qDRIFT, randomized product formulas) that can improve performance for non-sparse Hamiltonians but still require depth scaling with precision to suppress Trotter error; (b) sophisticated but resource-intensive approaches such as LCU, truncated Taylor series methods, and quantum signal processing/qubitization, which achieve better asymptotic precision but demand substantial overheads; (c) noise-aware or randomized compilation techniques to mitigate errors; (d) analog simulators which naturally implement continuous-time dynamics without Trotter error but lack digital flexibility. The paper also compares against the URCC algorithm (post-publication related work) that uses a different expansion (Dyson) and has bounded gate count per circuit O(t^2 μ), yet the present work achieves lower average gate count at matched variance by optimizing gate angles.

Core idea: Replace small-angle rotations by randomly applying larger-angle rotations with a known attenuation factor so that averaging over random circuits yields unbiased estimates of exact expectation values.

1. Small-angle from large-angle identity: For O with O^2=I, and parameters 0 < p ≤ 1, 0 ≤ τ' ≤ τ ≤ π/2, one has (1−p)I + p e^{i τ O} = λ e^{i τ' O}, with λ an attenuation factor and p determined by τ and τ'. Thus, sampling between identity and a fixed large-angle rotation reproduces a smaller rotation in expectation, up to known scaling.
2. Circuit-level extension: For a circuit with G such rotations, deleting a subset of rotations with probability 1−p and using angle τ (instead of τ') leads to an exact relation between averaged expectation values and the original circuit, with an overall attenuation λ^{2G}. The attenuation increases the shot cost but does not introduce bias.
3. Continuous-time limit and Poissonization: For a Hamiltonian H = Σ_n c_n O_n with O_n^2=I and c_n>0, the exact propagator e^{-iHt} is expressed as the limit of a product of small rotations (Trotter formula). Applying the random replacement to each term yields a sequence of random unitaries where, in the τ'→0 limit, events (non-identity gates) occur according to independent Poisson processes with rates c_n/ sin τ_n (τ_n are freely chosen fixed gate angles per term). Ordering all sampled gates by their random occurrence times produces a random circuit whose expectation value matches the exact time-evolved observable up to a known attenuation factor A = exp(−t Σ_n c_n tan(τ_n/2)).
4. Algorithm (time-independent case): Given H = Σ_n c_n O_n, unitary observable M, time t, initial state |ψ(0)>, and angles 0<τ_n≤π/2:

• Prepare |ψ(0)>.
• For each n, draw m_n ~ Poisson(t c_n / sin τ_n); draw m_n times uniformly in [0,t]; merge all times, sorting to form a sequence of gate indices k_1,...,k_M.
• Apply rotations e^{i τ_{k_m} O_{k_m}} in increasing time order.
• Apply M, then perform a backward evolution using the same sampled structure (with the same τ_n) to construct the overlap required for ⟨ψ(t)| M |ψ(t)⟩.
• Measure via a Hadamard test, dividing by A^2 = exp(−2 t Σ_n c_n tan(τ_n/2)). Repeat and average over circuit samples (and shots).
• For the Loschmidt amplitude ⟨ψ(0)| e^{iHt} |ψ(0)⟩, skip insertion of M and one of the evolutions, dividing by A.

5. Sampling and variance: With single-shot per circuit optimality (S=1), the estimator is unbiased with variance bounded by A^{-4}. Precision ε requires O(A^4 ε^{-2}) total shots; runtime equals (shots) × (average gates per shot).
6. Optimal gate angle (noiseless): Gate angles τ_n can be chosen to minimize runtime. At large t, a uniform angle choice τ_n=τ yields an optimal scaling with finite attenuation (A finite) and an average gate count per circuit E[N] ~ 4 t^2 μ^2 (μ = Σ_n c_n), independent of precision. Thus, average gate count is O(t^2 μ^2).
7. Noise-aware optimization: Under a simple depolarizing-noise-per-gate model with damping per two-qubit gate e^{−r_g}, the total signal attenuation becomes A^2 q, with q = exp(−2 r_g E[N]). Optimizing τ for noisy hardware modifies the optimal angle and preserves polynomial dependence on ε. The runtime remains polynomial in 1/ε, unlike Trotter where gate count must scale inversely with ε, leading to exponential shot cost in 1/ε under similar noise.
8. Distribution of gate counts: N_rot = Σ_n m_n is Poisson with mean E[N_rot] (scaling as O(t^2 μ^2) at optimal τ). Although unbounded, it is sharply concentrated (Chernoff bounds). If a hard maximum gate count per circuit is imposed, the discard probability can be made negligible; the maximum needed gate count to guarantee precision η scales as O(τ^2 t^2) + O(τ t log(1/η)). Average gate count remains independent of ε.
9. Time-dependent Hamiltonians: For H(t)=Σ_n c_n(t) O_n with c_n(t)≥0, replace constant-rate Poisson processes by time-dependent rates c_n(t)/sin τ_n. A time reparameterization z_τ(u)=∫_0^u c_n(s) ds maps to a standard Poisson process, yielding exact unbiased estimation with attenuation Λ = exp(−∫_0^t ds Σ_n c_n(s) tan(τ_n/2)). Steps are modified to sample times uniformly in [0, z_τ(t)] then invert z_τ to get physical times.
10. Background Hamiltonian trick: For a commuting subset O of terms, take τ_j→0 (for j∈O) first to turn those frequent commuting rotations into continuous evolution by H_background = Σ_{j∈O} τ_j O_j between non-commuting events. This reduces the attenuation by excluding commuting terms from A: A = exp(−t Σ_{j∉O} c_j tan(τ_j/2)).
11. Classical simulation of Clifford+T: The same probabilistic angle-reduction identity decomposes T into a randomized mixture of Clifford operations (I and S), enabling Monte Carlo averaging of Clifford circuits. The number of samples needed for precision ε scales as ~ 1/(ε^2 cos(π/8)^{2t}), matching the best known exponential-in-T classical simulation cost.

• Unbiased, zero-discretization-error estimator: Exact continuous-time expectation values ⟨ψ(t)|M|ψ(t)⟩ are obtained by averaging over random finite-depth circuits with a known attenuation factor. No Trotter error remains, even at finite average depth.
• Bounded average gate count independent of precision: At optimal gate angle (noiseless case), the average gate count per circuit scales as E[N] ≈ 4 t^2 μ^2, i.e., O(t^2 μ^2), independent of target precision ε. The attenuation A remains finite in this limit.
• Runtime with shot noise: With single-shot-per-circuit optimality and variance scaling, total runtime (shots × gates per shot) to reach precision ε scales as O(t^2 μ^2 ε^-2).
• Advantage for non-sparse Hamiltonians: Complexity depends on μ = Σ_n c_n (1-norm of coefficients), not on the number of terms, benefiting non-sparse cases.
• Noise-aware optimality and polynomial ε-dependence: Under a simple depolarizing gate-noise model, optimal τ modifies the total attenuation to A^2 q but maintains polynomial scaling in 1/ε, in contrast to Trotter where gate count must shrink with ε, leading to exponential shot cost under noise.
• Concentration of random gate counts: Although N_rot is unbounded (Poisson), its tail is exponentially suppressed (Chernoff bound). A practical cutoff can be set with negligible discard probability, and the maximal needed gate count to guarantee precision η scales as O(τ^2 t^2) + O(τ t log(1/η)).
• Time-dependent Hamiltonians: The method generalizes exactly via time-reparameterized Poisson processes, preserving unbiasedness and eliminating time-discretization error. Attenuation becomes Λ = exp(−∫_0^t ds Σ_n c_n(s) tan(τ_n/2)).
• Background (commuting) terms: Treating commuting subsets as a continuous background evolution reduces the attenuation factor by excluding those terms from A.
• Comparison with URCC: Matching variances implies our average gate count is half that of URCC (N_ours = N_URCC/2 on average) at the same shot overhead, due to gate-angle optimization.
• Numerical demonstrations: (i) Stretched H2O (STO-6G, 14 orbitals) Loschmidt amplitude with depolarizing noise per 2-qubit gate (p≈0.002) shows very good agreement for a ratio observable R(t) robust to depolarizing noise. (ii) 2D Ising model adiabatic state preparation matches exact energies, with error bars broadening with T as expected from attenuation. (iii) Against Trotter and qDRIFT on 2D Ising (h=2, 3×4), equal runtime (10 circuits × 10 shots): average absolute error 0.009 (ours) vs 0.064 (Trotter) and 0.135 (qDRIFT).
• Classical simulation of Clifford+T: Monte Carlo decomposition yields a sampling cost ∼ 1/(ε^2 cos(π/8)^{2t}), matching leading exponents of best-known classical simulators.

The work addresses the central challenge of obtaining exact real-time expectation values on digital quantum computers without discretization errors despite finite-depth circuits. By replacing small-angle evolutions with randomized large-angle gates and compensating via a known attenuation, the algorithm eliminates Trotter errors at the estimator level. Mapping the limit to Poisson processes provides a rigorous foundation and enables practical sampling procedures with adjustable gate angles that trade entanglement and attenuation. This directly benefits applications requiring high precision (e.g., quantum chemistry) and adiabatic state preparation, where Trotter errors can lead to heating. The dependence on μ rather than sparsity makes it attractive for dense Hamiltonians common in chemistry. The method’s runtime is polynomial in 1/ε and exhibits bounded average depth, making it particularly suitable for NISQ devices with moderate coherence. The extension to time-dependent Hamiltonians preserves exactness and avoids time-discretization errors. Empirical tests confirm practical gains over Trotter and qDRIFT at equal runtime, and comparison with URCC highlights lower average gate counts via gate-angle optimization. The framework also connects quantum evolution to a path-integral-like limit (τ→π/2) and suggests classical simulation avenues (Clifford+T) with optimal exponents, further underscoring conceptual breadth.

The paper presents an algorithm that computes exact time-evolved expectation values on digital quantum hardware with finite average circuit depth independent of the desired precision, eliminating discretization (Trotter) errors. It achieves average per-circuit gate count O(t^2 μ^2) and total runtime O(t^2 μ^2 ε^-2), is well-suited for non-sparse Hamiltonians, and naturally extends to time-dependent cases. Numerical experiments showcase strong performance and robustness under noise, outperforming standard Trotter and qDRIFT at equal runtime. The approach bridges regimes from fully quantum (small τ) to path-integral-like classical trajectories (τ≈π/2), and supports efficient classical simulation strategies for Clifford+T circuits with optimal scaling. Future directions include deeper investigation of chemistry and adiabatic applications, refined noise models and mitigation, optimized angle selection strategies (including background commuting subsets), and devising related algorithms for sampling tasks rather than expectation values.

• The number of gates per circuit is random and unbounded (Poisson); although tightly concentrated, practical implementations may impose cutoffs that can introduce bias if not managed with small discard probabilities.
• The attenuation factor, while known, increases the required number of shots; optimal angle choices balance gate count against shot overhead and device noise, but high precision still demands many total measurements.
• The noise model analyzed (global depolarizing per gate) is a simplified approximation; real devices exhibit more complex noise, potentially affecting performance and optimal parameters.
• The algorithm targets expectation values of unitaries (e.g., via Hadamard tests); adapting to general (non-unitary) observables or to sampling output distributions is nontrivial. The authors note it is not currently adapted for sampling-based tasks (e.g., classical optimization sampling).
• Resource estimates depend on Hamiltonian coefficient 1-norm μ; while favorable for dense Hamiltonians, practical compilation overheads for multi-qubit Pauli rotations and measurement schemes may impact performance on specific platforms.