Source-linked AI summary
Quantum algorithm for simulating real time evolution of lattice Hamiltonians
Jeongwan Haah, Matthew B. Hastings, Robin Kothari, Guang Hao Low
TL;DR
The paper asks how efficiently quantum computers can simulate time evolution under geometrically local lattice Hamiltonians. It decomposes the evolution using Lieb-Robinson bounds and achieves quasilinear gate complexity, while proving matching worst-case lower bounds even under weaker output requirements.
Problem
Simulating time evolution of geometrically local lattice Hamiltonians is important for modeling physical systems, but efficient gate-complexity guarantees and lower bounds are needed.
Method
The algorithm decomposes full time evolution into products of smaller local time-evolution unitaries using Lieb-Robinson bounds.
Results
O(nT polylog(nT/ε)) geometrically local two-qubit gates with depth O(T polylog(nT/ε)) suffice, and Ω̃(nT) gates are necessary in worst-case one-dimensional instances.
Takeaways & Limitations
The algorithm is quasilinear in spacetime volume up to polylogarithmic factors, and the gate count is optimal up to those factors for worst-case instances.
Takeaways & Limitations
The stated lower bound applies only for T ≤ 4n, while the main simulation guarantee uses spectral-norm, hence strong, unitary approximation.
Abstract
from arXiv · showhide
We study the problem of simulating the time evolution of a lattice Hamiltonian, where the qubits are laid out on a lattice and the Hamiltonian only includes geometrically local interactions (i.e., a qubit may only interact with qubits in its vicinity). This class of Hamiltonians is very general and is believed to capture fundamental interactions of physics. Our algorithm simulates the time evolution of such a Hamiltonian on $n$ qubits for time $T$ up to error $ε$ using $\mathcal O( nT \mathrm{polylog} (nT/ε))$ gates with depth $\mathcal O(T \mathrm{polylog} (nT/ε))$. Our algorithm is the first simulation algorithm that achieves gate cost quasilinear in $nT$ and polylogarithmic in $1/ε$. Our algorithm also readily generalizes to time-dependent Hamiltonians and yields an algorithm with similar gate count for any piecewise slowly varying time-dependent bounded local Hamiltonian. We also prove a matching lower bound on the gate count of such a simulation, showing that any quantum algorithm that can simulate a piecewise constant bounded local Hamiltonian in one dimension to constant error requires $\tilde Ω(nT)$ gates in the worst case. The lower bound holds even if we only require the output state to be correct on local measurements. To our best knowledge, this is the first nontrivial lower bound on the gate complexity of the simulation problem. Our algorithm is based on a decomposition of the time-evolution unitary into a product of small unitaries using Lieb-Robinson bounds. In the appendix, we prove a Lieb-Robinson bound tailored to Hamiltonians with small commutators between local terms, giving zero Lieb-Robinson velocity in the limit of commuting Hamiltonians. This improves the performance of our algorithm when the Hamiltonian is close to commuting.
1 Introduction
The paper studies efficient quantum simulation of geometrically local lattice Hamiltonians, whose dynamics model fundamental physical interactions. It gives near-optimal algorithms and gate-complexity lower bounds for this problem.
- Problem setting: The simulation problem takes an initial n-qubit state, a Hamiltonian, time T, and error ε, and outputs an approximation to the evolved state or unitary.Approximation is measured using spectral norm distance between the ideal and implemented unitaries.
- Problem setting: Geometrically local Hamiltonians act on adjacent qubits arranged in a constant-dimensional Euclidean lattice and model fundamental physical forces and many condensed matter systems.The paper distinguishes geometrically local interactions from non-geometrically local k-local interactions.
- Algorithmic results: O(nT polylog(nT/ε)) local two-qubit gates and depth O(T polylog(nT/ε)) suffice for lattice Hamiltonian simulation.The algorithm uses geometrically local gates and only polylog(nT/ε) ancilla qubits.
- Lower bounds: A piecewise constant bounded one-dimensional Hamiltonian requires Ω̃(nT) two-qubit gates for constant-error simulation in the worst case.This lower bound remains valid with unlimited ancillas, nonlocal gates, and possibly infinite gate sets.
- Lower bounds: The Ω̃(nT) lower bound persists when correctness is required only for local observables, including a computational-basis measurement on the first output qubit.For T ≤ n, the local-measurement lower bound is Ω̃(T^2).
- Algorithmic results: The algorithm decomposes full time evolution into many smaller unitaries acting on logarithmic-size Hamiltonian instances, using Lieb-Robinson bounds.The smaller simulations require exponentially small error, enabling polylogarithmic dependence on 1/ε.
2 Algorithm and analysis
The algorithm recursively decomposes local Hamiltonian evolution into smaller overlapping block evolutions using Lieb-Robinson bounds, then simulates those blocks with existing methods. This achieves near-linear spacetime gate complexity and polylogarithmic accuracy dependence, including for suitable time-dependent Hamiltonians.
- Complexity: O(Tn polylog(Tn/ε)) two-qubit local gates and O(T polylog(Tn/ε)) depth suffice for accuracy ε.The overall gate count follows when block algorithms are polynomial in block size and polylogarithmic in accuracy.
- Decomposition: The decomposition reduces full time evolution to products of smaller unitaries supported on overlapping regions, with exponentially small error in the overlap width.The construction repeatedly applies the decomposition lemma to large regions, producing small block evolutions that can be implemented separately.
- Lieb-Robinson analysis: Lieb-Robinson bounds provide exponentially decaying influence across separated regions, while slower-decaying interactions require larger overlaps.For strictly local interactions, the relevant error is O(e^-µℓ) at fixed time.
- Decomposition: For a one-dimensional chain, the method applies the decomposition across blocks and repeats it over O(T) time slices for total evolution time T.Each time slice has constant-scale duration, and the full propagator is written as a product of slice propagators.
- Complexity: Choosing overlap width ℓ = O(log(TL/ε)) yields m = O(TL^D/ℓ^D) blocks, each acting on O(ℓ^D) qubits.The block count and overlap width control the accumulated approximation error and the cost of implementing the local evolutions.
- Time dependence: The same gate complexity extends to piecewise slowly varying bounded time-dependent Hamiltonians when suitable time-dependent simulation algorithms are used.The required slowly varying condition supports the underlying block-simulation routines; the decomposition itself is applied across time slices.
3 Optimality
The optimality proof constructs hard one-dimensional Hamiltonians by embedding local quantum circuits into piecewise constant Hamiltonian evolution. Counting functions and distinguishable unitaries shows that simulation requires near-linear gate count in the spacetime volume, even under relaxed circuit and output assumptions.
- Hamiltonian embedding: Every depth-T circuit of local two-qubit gates can be represented exactly as evolution under a piecewise constant bounded one-dimensional Hamiltonian for time T.The construction assigns a two-qubit Hamiltonian to each circuit layer and concatenates the resulting time intervals.
- Counting argument: Local depth-T circuits on n qubits can compute at least 2^eΩ(Tn) distinct Boolean functions when 2 ≤ n ≤ T ≤ 2n.The proof groups qubits into blocks, computes many functions per block, and combines outputs using parity.
- Counting argument: Circuits using G arbitrary non-local two-qubit gates compute at most 2^eO(G log n) Boolean functions, even with unlimited ancillas.This upper bound applies to any gate set, including possibly infinite sets.
- Unitary counting: A parallel unitary-counting argument constructs 2^eΩ(Tn) pairwise distinguishable unitaries, while G-gate circuits can realize at most 2^eO(G log n) such unitaries.This establishes the corresponding lower bound for approximating time-evolution operators, not only measured output functions.
4 Discussion
The discussion extends the decomposition to fermionic systems and compares its geometric assumptions with alternative simulation choices. It also notes that smaller Lieb-Robinson velocity improves accuracy, while ancilla-free implementations retain non-polylogarithmic accuracy dependence.
- Discussion: The decomposition is more accurate when local Hamiltonian terms have smaller commutators and therefore smaller Lieb-Robinson velocity.Lieb-Robinson bounds with exponential separation dependence extend beyond Euclidean space, but suitable decompositions may be graph-specific.
- Discussion: Without ancillas, high-order Suzuki formulas give arbitrarily small power-law accuracy dependence, but cannot achieve polylogarithmic dependence on Tn/ϵ.The obstruction is the exponential prefactor in the resulting gate complexity.
- Fermionic systems: The method applies to fermionic Hamiltonians because even fermion parity preserves Lieb-Robinson bounds without modification.Majorana operators provide a convenient representation for reducing fermionic systems to qubit Hamiltonians.
- Fermionic systems: In one dimension, Jordan-Wigner strings cancel when fermion ordering follows the chain, yielding a local qubit Hamiltonian.Periodic boundaries can be handled by viewing the chain as a doubled finite line, at the cost of doubling qubit density.
- Fermionic systems: In higher dimensions, auxiliary fermions cancel Jordan-Wigner strings and preserve the same gate and depth complexity with O(1) auxiliaries per target fermion.Auxiliary fermions can be reused for interaction terms involving fermions that are close in the chosen ordering.
A Heisenberg model benchmark
A numerical Heisenberg-model benchmark tests the Lieb-Robinson block decomposition and its error scaling. The implementation combines blockwise quantum signal processing and qubitization, with explicit Clifford+T gate counts.
- Model: The antiferromagnetic one-dimensional Heisenberg model with inhomogeneous fields provides a numerical benchmark for quantum simulation algorithms.With non-uniform field coefficients, the model generally requires numerical treatment.
- Error behavior: The decomposition error is almost independent of overlap position and decreases exponentially with overlap size ℓ.For the tested small systems, errors were evaluated over a range of evolution times up to ℓ = 9.
- Complexity: The full evolution decomposes into m = O(Tn/ℓ) blocks, with total error bounded by m(ϵLR + ϵ□), requiring ℓ = O(log(Tn/ϵ)).The two contributions arise from decomposing the global evolution and approximately simulating each block.
- Error behavior: Choosing ℓ = O(t + log(1/ϵLR)) is consistent with the Lieb-Robinson error bound.The numerical fit is used to estimate overlap sizes needed for a target decomposition error.
- Complexity: An alternative recursive block schedule reduces the number of blocks to m = 3Tn/(2tℓ), potentially improving efficiency depending on the block simulator.The schedule uses 2m/3 blocks of size ℓ for time t and m/3 blocks of size 2ℓ for time 2t.
- Benchmark: Figure 3 reports Clifford+T gate counts for T = n, ϵ = 10^-3, and uniformly random hj ∈ [−1, 1], using overlap sizes ℓ = 7, 8, 9.The comparison includes blockwise QSP/qubitization, undecomposed QSP with complexity ˜O(n^3), and optimized Lie-Trotter-Suzuki simulations.
B Further algorithmic improvements
The decomposition can accommodate a single interaction term much stronger than the remaining local terms without subdividing the entire system’s time step. A cut isolates the strong term so its error dependence is controlled.
- Inhomogeneous interactions: A strong term with norm J ≫ 1 need not force global time-step subdivision into ⌈J⌉ pieces.The analysis instead localizes the treatment of the strong interaction around a cut.
- Inhomogeneous interactions: Cutting at the strong term places it in Hbd, making the truncation error linear in J while retaining exponential decay in overlap size ℓ.This isolates the large interaction strength from the rest of the chain in the error analysis.
B.2 Reducing number of layers in higher dimensions
Higher-dimensional decompositions use colored tessellations whose same-color cells are well separated, reducing the number of evolution layers. The general construction requires 2α − 1 layers for an α-colorable tessellation.
- Higher dimensions: The hyperplane construction uses 3^D layers per unit time, which grows rapidly with spatial dimension D.A more uniform tessellation, such as a 4-colorable BCC construction in three dimensions, can reduce this layer count.
- Two dimensions: A 3-colorable hexagonal tiling decomposes two-dimensional unit evolution through five forward and backward block-evolution stages.The stages alternate evolution on color unions, overlaps, rings, and enlarged cells.
- Two dimensions: The five-stage construction uses separated rings formed by the ℓ-neighborhood of one color’s cells to complete evolution over the full lattice.Backward evolution acts on overlap rings, followed by forward evolution on enlarged hexagons.
- General construction: For an α-colorable tessellation with well-separated same-color cells, unit evolution requires 2α − 1 layers.The construction proceeds inductively by evolving α − 1 colors, then correcting overlaps with backward and forward evolutions.
- General construction: A regular D-dimensional lattice admits a (D + 1)-colorable tessellation obtained from a triangulation and subsequent fattening.This provides a general coloring scheme for regular lattices.
- Higher dimensions: In the BCC construction, lattice points are colored by c = x + y + z mod 4, and same-color cells are separated by the tessellation geometry.The associated Voronoi cells form a valid four-colored tessellation for the decomposition.
C.1 Introduction and Assumptions
The appendix develops Lieb-Robinson bounds for Hamiltonians with small inter-term commutators and exponentially decaying interactions. Under these assumptions, the resulting velocity scales as √η and the bounds extend to piecewise continuous time dependence.
- Motivation: Small commutators can reduce the Lieb-Robinson velocity, improving the method’s accuracy; the velocity vanishes in the commuting limit.The appendix contrasts this behavior with methods that do not generally improve under small commutators.
- Assumptions: The exponential-decay assumptions bound interaction strength through Eq. (15), without imposing ∥hX∥≤1.The assumptions use the cardinality |X| and are stronger than some earlier exponential-decay conditions.
- Main bound: For operators supported on separated sets, the appendix proves a Lieb-Robinson bound under assumptions (14,15).The result is stated as Lemma 12 for operators supported on sets X and Y.
- Time dependence: The bound remains valid for piecewise continuous time-dependent Hamiltonians when the assumptions hold at every time.The proof notation is time-independent, but the lemma uses a time-ordered exponential in the time-dependent case.
C.3 Lieb-Robinson Velocity
This section derives the exponential-decay Lieb-Robinson estimate by recursively bounding operator-growth quantities. The resulting lemma follows after summing the bounds over recursion order.
- Recursive bounds: The proof uses assumption (15) to bound sums over linked site sets with exponentially decaying distance dependence.The argument repeatedly applies propositions controlling sums involving interaction norms, diameters, and distances.
- Recursive bounds: The dependence on |X|2 in assumption (15) is needed when recursively bounding the kth line of the expansion.The proof alternates bounds for odd and even recursion levels before applying the linked-set estimate repeatedly.
- Conclusion: Summing the recursion over k yields Lemma 12.The preceding estimates are combined to establish the main exponential-decay Lieb-Robinson result.
C.4 Proof of Lemma 5
The proof extends the earlier argument to obtain a variant of Lemma 5 under exponential-decay locality. Strict locality supplies a stronger bound because sufficiently short interaction paths vanish.
- Exponential decay: A modification of the preceding exponential-decay argument proves a variant of Lemma 5.The proof bounds contributions outside a region by changing the distance measure to the complement of that region.
- Strict locality: For strictly local interactions, hX=0 when diam(X)>1, which implies the exponential-decay assumption and suffices for Theorem 1.The strictly local condition permits arbitrary µ>0 after estimating ζ as a function of µ.
- Strict locality: Strict locality eliminates the first ℓ−1 terms of the expansion when no interaction path of that length can connect X and Y.The proof uses the absence of nonzero links between separated supports.
- Conclusion: The proof concludes the stated strictly local variant of Lemma 5.The final step uses η≤1 after applying the modified bounds.
C.5 Higher Order Commutators
The appendix sketches stronger Lieb-Robinson bounds under assumptions on higher-order commutators. It generalizes the recursive proof by introducing additional operator-growth quantities.
- Higher-order bounds: Bounding higher-order commutators can produce a better Lieb-Robinson bound for sufficiently small η′.The paper proposes extending the two-commutator analysis to third and higher orders.
- Proof strategy: The generalized proof introduces EB(X,Y,t) to control the evolution of DB(X,t) through higher-order commutator terms.The quantities CB, DB, and EB are bounded recursively in succession.
- Assumptions: For commutator order β, the corresponding locality assumption replaces Eq. (15) with a higher-order condition.The paper identifies β=2 with Eq. (14) and β=3 with Eq. (50).
C.6 Strictly Local Hamiltonians
The section develops bounds for strictly local Hamiltonians under local commutator assumptions, deriving Lieb-Robinson behavior and its implications for analytic approximation and simulation circuits.
- Assumptions: Bounded-range interactions can be represented using a rescaled metric, while local commutators are assumed to be bounded by a constant K.The analysis does not impose a uniform numerical bound on individual term norms, but each term remains bounded for finite systems.
- Lieb-Robinson bound: The commutator expansion bounds distant influence by sums over interaction-graph sequences whose lengths must reach the separation between regions.The derivation repeatedly substitutes recursive bounds and relaxes sequence restrictions to obtain a tractable positive sum.
- Lieb-Robinson bound: For bounded-degree interaction graphs, the Lieb-Robinson velocity is bounded by a constant times dK.A simpler finite-series argument gives vLR = O(dK), though it can be slightly worse than the infinite-series bound on expanding graphs.
- Polynomial approximation: Functions analytic on a Bernstein ellipse admit Chebyshev approximations with degree O(log(1/ϵ)) for accuracy ϵ.The result follows from exponentially bounded Chebyshev coefficients and supports polylogarithmic-precision approximation.
- Hamiltonian simulation: Hamiltonian simulation is outlined in three steps: encode H, construct a qubiterate, and use quantum signal processing to transform its eigenphases.The qubiterate eigenphases are related to Hamiltonian eigenvalues, while the signal-processing sequence uses controlled-W operations and single-qubit rotations.
- Hamiltonian simulation: Jacobi-Anger truncation approximates e^(-iαt sin θ), with error controlled by the tail sum of Bessel-function coefficients.The required phase angles can be precomputed classically given α ≥ ∥H∥ and t.
E.1 Encoding coefficients in reflections
The section constructs block-encoding oracles for Hamiltonians expressed as Pauli sums and modifies coefficient encoding so reflection-based qubitization remains applicable.
- Oracle construction: A general Hamiltonian is represented as a linear combination of M Pauli operators, whose coefficients are encoded through oracle construction.The resulting oracles O and G satisfy the required block-encoding relation and have asymptotically similar gate complexity.
- Coefficient encoding: Arbitrary state preparation can be excessively costly when many Hamiltonian coefficients are identical.Uniform or structured coefficient patterns can instead use Hadamards, arithmetic with amplitude amplification, added identity terms, or unary control logic.
- Coefficient encoding: Encoding coefficients in unitary operators avoids putting all coefficient information in the preparation state, but a direct construction may violate O^2 = 1_as.The direct exponential encoding therefore does not meet the prerequisite of the simple qubitization circuit.
- Coefficient encoding: A two-qubit circuit modification encodes coefficient information in unitary operators while preserving the reflection condition O^2 = 1_as.This supplies the structural property needed by the simple qubitization construction.
- Circuit optimization: When coefficient-dependent rotation angles β_j vanish, the corresponding Q_j operation can be replaced by the identity.This yields a situational constant-factor gate reduction for sparse deviations from unit coefficients.