A quantum evolution can preserve probability perfectly and still be the wrong evolution. Hamiltonian simulation is the task of reproducing the right dynamics, with an error that has been defined, bounded and checked.
Suppose a system is governed by two interactions. You know how to implement each interaction separately, but not how to implement both simultaneously. Can a sequence of short, separate actions reproduce their combined effect? The answer is often yes—but the order of the actions matters, and the approximation must be controlled.
This is the central problem behind Trotter–Suzuki product formulas. A complicated exponential is replaced by a product of simpler exponentials. The mathematical responsibility is then to understand what the replacement preserves, what it changes, how the error accumulates and how much computation is required to make that error acceptable. Modern analyses express this error using commutators rather than treating every Hamiltonian as an undifferentiated large matrix. [1,2]
This is Guide 13 in the BTT Quantum Mathematics series. It is university-level enrichment, not a statement of school examination requirements. The starting tools are matrix multiplication, complex exponentials, operator norms and elementary series. The earlier guide to unitary evolution and the Schrödinger equation supplies the physical meaning; this article supplies the approximation machinery.
Reading route: the simulation task → why splitting introduces error → a genuine error bound → symmetric and higher-order formulas → a numerical comparison → worked practice.
1. Define the simulation task before choosing an algorithm
For a time-independent Hermitian Hamiltonian H, the target closed-system evolution is U(t) = exp(−iHt/ℏ). A simulator implements another transformation, V(t), intended to approximate U(t). To make the question precise, specify H, the initial state or allowed input states, the evolution time, the observable to be predicted and the allowed error.
These are different tasks. Approximating U on every normalised input is stronger than approximating one particular evolved state. Approximating that state is stronger than predicting one expectation value accurately. A simulation that gives the correct magnetisation for one initial state has not thereby reproduced every amplitude of the full evolution.
Throughout the finite-dimensional examples below, we set ℏ = 1 and express H in compatible inverse-time units. This is a units convention, not a claim that the physical constant has disappeared. With energy-valued H, restore ℏ in every exponential and in the corresponding error bounds. Time t and H must combine into a dimensionless exponent.
2. Split the Hamiltonian into implementable pieces
Write H = A + B, where A and B are Hermitian. Perhaps A contains interactions between neighbouring spins while B contains local fields. The aim is to exploit easy implementations of exp(−iAδ) and exp(−iBδ), with δ a small time step.
The decomposition is part of the algorithm. Two decompositions of the same H can lead to different circuit depths, different commutator bounds and different opportunities for parallel execution. Merely counting the number of summands is not enough. Their algebraic relationships and the hardware operations needed to implement them also matter.
A useful first question is therefore not “How many Trotter steps shall we use?” It is “Which parts can be performed exactly, which parts commute, and which noncommuting boundaries actually need approximation?” Avoid approximating a relationship that can already be implemented exactly.
3. Why exp(A + B) is not generally exp(A)exp(B)
For ordinary numbers, exponentials convert addition into multiplication. For matrices, the familiar identity requires commutation. The commutator is [A,B] = AB − BA. If it is zero, the two actions are compatible with this factorisation. If it is nonzero, changing the order changes the transformation.
Expand both sides for a small δ. The exact evolution through second order is I − i(A+B)δ − (A²+AB+BA+B²)δ²/2, plus terms of order δ³. The split product exp(−iAδ)exp(−iBδ) instead has second-order contribution −(A²+B²)δ²/2 − ABδ².
Subtracting gives a particularly revealing result:
exp(−iAδ)exp(−iBδ) − exp(−i(A+B)δ) = −δ²[A,B]/2 + O(δ³).
The discrepancy starts with the commutator. It is not caused by a loss of normalisation, an inaccurate matrix multiplication or random measurement noise. It is caused by replacing simultaneous noncommuting dynamics with a sequence of separate dynamics.
4. Worked example: a commuting split is exact
Let A = Z and B = 2Z, where Z has diagonal entries 1 and −1. Because both matrices are multiples of the same operator, AB = BA. Their sum is 3Z, so exp(−it(A+B)) = exp(−i3tZ).
The split evolution is exp(−itZ)exp(−i2tZ). Its upper diagonal entry is exp(−it)exp(−i2t) = exp(−i3t); its lower entry is exp(it)exp(i2t) = exp(i3t). It is exactly the target evolution for every t, not merely for small t.
This example explains why a bound based only on the sizes of A and B can be wasteful. Both operators may be large while the splitting error is identically zero. A commutator-sensitive description recognises that special structure.
5. Worked example: X and Z do not commute
Now let A = X and B = Z. The Pauli matrices satisfy X² = Z² = I and XZ = −ZX. In the convention used here, [X,Z] = −2iY. Since the operator norm of Y is one, the norm of the commutator is two.
The exact Hamiltonian is H = X + Z. Squaring it gives H² = X² + XZ + ZX + Z² = 2I. Separate the even and odd terms in the exponential series to obtain
U(t) = cos(√2t)I − i sin(√2t)(X+Z)/√2.
This closed form will serve as an independent reference for our numerical checks. A trustworthy teaching example should have a result that can be checked in more than one way: here by a matrix exponential, by the Pauli identity and by direct evaluation of the two-by-two matrix.
6. The first-order Lie–Trotter formula
Divide t into r equal steps, δ = t/r. Define S₁(δ) = exp(−iAδ)exp(−iBδ). The approximation is V₁(t;r) = S₁(t/r)r. For bounded finite-dimensional operators, this converges to exp(−i(A+B)t) as r increases. [1]
Read matrix products from right to left when describing chronological action on a column vector: in S₁ as written, the B evolution acts first, followed by A. A circuit drawn left to right in time can therefore look reversed relative to its algebraic product. Declare the convention rather than silently switching it.
Each individual factor is unitary. Their product is unitary. Consequently, every finite-r approximation preserves the norm exactly in ideal arithmetic. Increasing r is about bringing the dynamics closer to the intended dynamics, not repairing probability conservation.
7. A useful error bound, not just an asymptotic slogan
Let ‖M‖ denote the operator norm induced by Euclidean vector length: the largest value of ‖Mv‖ over normalised vectors v. For bounded Hermitian A and B, the one-step first-order error obeys
‖S₁(δ) − exp(−i(A+B)δ)‖ ≤ δ²‖[A,B]‖/2.
One way to see the bound is to differentiate the split product. Its instantaneous generator contains exp(−iAs)Bexp(iAs) instead of B. The difference between these operators is an integral of a unitarily conjugated commutator. Its norm is at most s‖[A,B]‖. Integrating that mismatch from zero to δ gives δ²‖[A,B]‖/2.
The unitary factors surrounding the mismatch have norm one, which is why no artificial exponential growth factor is needed in this finite-dimensional real-time bound. This derivation also identifies its assumptions: bounded Hermitian generators and unitary real-time propagation. Do not transfer the same estimate unchanged to an unbounded differential operator or a dissipative evolution. [1,4]
8. From local error to global error
A local error concerns one step. A global error concerns the full time interval. For unitary S and T, expand Sr − Tr as a telescoping sum of r terms, each containing one factor S−T surrounded by unitary factors. Taking norms gives ‖Sr−Tr‖ ≤ r‖S−T‖.
Substituting δ = t/r into the one-step estimate yields
‖V₁(t;r) − U(t)‖ ≤ t²‖[A,B]‖/(2r).
Notice the power of r. The local error is proportional to δ², but there are r steps. The global error therefore decreases like 1/r at fixed t. Calling a method first order refers to this global step-size scaling, not to the first nonzero term in every intermediate expansion.
9. Worked example: choose a conservative step count
For H = X+Z at t = 1, the commutator norm is two. The bound simplifies to ‖V₁−U‖ ≤ 1/r. To certify an operator error no larger than 0.01 using this bound alone, choose r ≥ 100.
This is a sufficient step count, not a claim that 99 steps must fail. The actual error may be smaller because a worst-case bound does not exploit every cancellation. A good report distinguishes a certified upper bound from an observed error and from an empirical estimate extrapolated from smaller instances.
If t doubles while the tolerance stays fixed, this particular bound asks for four times as many steps. If A and B are both multiplied by a factor c while time stays fixed, their commutator is multiplied by c². The strength of the interactions belongs in the calculation; the number of steps is not determined by qubit count alone.
10. Symmetric splitting cancels the leading error
The second-order symmetric formula is S₂(δ) = exp(−iAδ/2)exp(−iBδ)exp(−iAδ/2). Its global approximation is V₂(t;r) = S₂(t/r)r. This construction is often called Strang splitting or second-order Suzuki–Trotter splitting. [2,3]
The symmetry is doing real mathematical work. Reversing the sign of δ produces the inverse sequence, so S₂(−δ) = S₂(δ)−1. In a local logarithmic expansion around δ = 0, this symmetry removes even powers. The leading correction to the target generator is therefore cubic rather than quadratic.
The local error is O(δ³), controlled by nested commutators such as [A,[A,B]] and [B,[A,B]]. Over a fixed interval, the global error is O(t³/r²), with constants depending on those operators. Doubling r should eventually reduce the error by about a factor of four, once the calculation has entered its asymptotic convergence regime.
11. Second order is not simply three times the work
A written S₂ step contains three exponentials, compared with two for S₁. But adjacent half-A steps can be merged. In S₂(δ)S₂(δ), the two middle factors exp(−iAδ/2)exp(−iAδ/2) combine exactly into exp(−iAδ).
Across r steps, this leaves r B exponentials and r+1 A exponentials: 2r+1 elementary exponentials rather than 3r, assuming the implementation permits the merger. The full physical gate count still depends on how each exponential is synthesised.
At a fixed target accuracy, second order may need far fewer steps than first order. Comparing one step against one step therefore answers the wrong resource question. Compare complete implementations that achieve the same tolerance for the same Hamiltonian, time and output task.
12. Fourth order through a controlled composition
A standard fourth-order Suzuki formula builds on S₂. Define p = 1/(4−41/3), approximately 0.41449077. Then
S₄(δ) = S₂(pδ)² S₂((1−4p)δ) S₂(pδ)².
The coefficients add to one overall time step: 4p+(1−4p)=1. Their cubic combination satisfies 4p³+(1−4p)³=0, which cancels the leading cubic contribution from the symmetric second-order approximation. The resulting local error is fifth order and the fixed-time global error is fourth order. [3]
The middle coefficient is negative, approximately −0.65796309. For an implementable unitary rotation, a negative evolution interval means using the inverse rotation. It does not mean reversing the passage of laboratory time. Nor does it license treating an irreversible noise channel as though it had a physically valid inverse.
13. Higher order has a cost
Recursive Suzuki constructions increase accuracy order while increasing the number of factors. For order 2k, the familiar symmetric recursion uses pk = 1/(4−41/(2k−1)) and five scaled copies of the preceding even-order formula. This is a systematic way to cancel successive error terms, not a promise that the highest available order is always best. [3]
Longer formulas may require more imperfect gates. Negative coefficients may complicate implementation. The constants multiplying the asymptotic error can also become important. A method with superior asymptotic order can lose at the finite precision, finite time and finite gate budget of a particular task.
Order is one design variable. The others include decomposition, term ordering, number of steps, permitted gate set and the error criterion. Choosing a formula is a constrained comparison, not a ranking by a single integer.
14. Numerical check: the same noncommuting problem three ways
The following values are direct numerical calculations for the explicitly stated two-by-two example H = X+Z, t = 1. The error is the operator norm of V−U. The first-order factor order is exp(−iXδ)exp(−iZδ); second and fourth order use the definitions above. These are classical matrix calculations, not quantum-hardware measurements or a performance benchmark.
| Steps r | First order | Second order | Fourth order |
|---|---|---|---|
| 1 | 0.799214 | 0.313666 | 0.0127194 |
| 2 | 0.362410 | 0.0711304 | 0.000688402 |
| 4 | 0.176261 | 0.0172728 | 0.0000420565 |
| 8 | 0.0875127 | 0.00428607 | 0.00000261582 |
| 16 | 0.0436792 | 0.00106951 | 0.000000163300 |
At sufficiently fine steps, doubling r roughly halves the first-order error, quarters the second-order error and divides the fourth-order error by sixteen. The coarse first row does not have to obey those ratios exactly: an asymptotic scaling statement concerns the small-step regime.
All three approximations remain unitary to numerical precision. That makes this table a direct counterexample to the claim that norm preservation proves simulation accuracy. A unitary answer can still be far from the intended unitary answer.
15. Pauli strings make elementary exponentials manageable
A Hamiltonian on qubits is often represented as H = Σj cjPj, where each Pj is a tensor product of Pauli matrices and identities. A Hermitian Pauli string satisfies P² = I. Its exponential therefore has the simple exact form exp(−iaP) = cos(a)I − i sin(a)P.
This identity follows from separating even and odd powers in the exponential series, just as in the X+Z example. No approximation is needed inside that identity. Product-formula error enters when exponentials of different noncommuting strings are multiplied to approximate the exponential of their sum.
Keep another factor-of-two convention visible: the usual rotation RZ(θ) equals exp(−iθZ/2). To implement exp(−icδZ), the rotation angle is θ = 2cδ, not cδ. A perfect product formula with every gate angle halved will faithfully simulate the wrong Hamiltonian.
16. Group commuting terms before splitting them
Consider a spin chain with HZ = JΣ ZjZj+1 and HX = hΣ Xj. Every ZZ term commutes with every other ZZ term, even when adjacent terms share one qubit, because their overlapping factors are both Z. The X terms also commute with one another.
Within each group, exponentiation can be factored exactly. The approximation concerns the separation between HZ and HX, which generally do not commute. This is the structure used in familiar transverse-field Ising simulations. [2]
Commutation and parallel execution are not identical. Two commuting operations that share hardware may still be impossible to execute simultaneously. Algebra determines whether reordering is valid; architecture determines whether that reordering saves depth. Both descriptions are needed for a credible resource estimate.
17. Translate state error into observable error
Suppose ‖U−V‖ ≤ ε and the input state ψ is normalised. Set a = Uψ and b = Vψ. Then ‖a−b‖ ≤ ε. For a bounded observable O, insert and subtract one mixed term to obtain
|⟨a|O|a⟩ − ⟨b|O|b⟩| ≤ 2‖O‖ε.
The estimate follows because each of the two difference terms contains a−b once and a normalised state once. For an observable with norm one, an operator error of 0.005 gives a conservative expectation-value error bound of 0.01.
This is a useful connection between an abstract operator bound and the quantity a student may actually want to predict. The converse does not follow: one expectation value can agree even when the underlying states disagree substantially. A single observable is a partial test.
18. Keep the error budget separated
A simulation can contain model error, basis-truncation error, product-formula error, gate-synthesis error, physical noise, numerical roundoff and finite-sample measurement error. They arise at different points in the calculation and are not interchangeable.
Increasing r reduces ideal splitting error but usually increases the number of physical operations. Taking more measurement shots reduces sampling uncertainty but does not correct a biased Hamiltonian model. Increasing floating-point precision reduces roundoff but does not improve a poor initial-state preparation.
Do not add unrelated quantities simply because all are called “errors”. Operator norms, infidelities, probability deviations and standard errors need a common interpretation before a combined bound is meaningful. A useful report names each contribution and the route by which it affects the final observable.
19. Time-dependent and unbounded problems need new assumptions
If H changes with time, the exact evolution is a time-ordered propagator. Freezing H at one time introduces a separate approximation, even before splitting it into terms. Noncommutation between Hamiltonians at different times also matters.
For differential operators, such as the kinetic-energy operator on an infinite-dimensional function space, the operator may be unbounded. Domain questions and regularity of the evolving state become relevant. The finite-dimensional norm argument above cannot simply be pasted onto that setting. Work on time-dependent unbounded Hamiltonians treats precisely these additional issues. [4]
A finite matrix obtained by truncation is bounded, but that only certifies the matrix problem you constructed. You must still assess whether the truncation represents the original physical problem accurately. Solving a discretised problem very precisely does not remove its discretisation error.
20. The mathematical checks that deserve to travel with the answer
Begin with the zero-time check: V(0) must equal I. Next check Hermiticity of the supplied Hamiltonian, unitarity of each real-time factor and consistency of units and rotation conventions. A sign error in exp(−iHt) versus exp(+iHt) reverses the intended dynamics while preserving norm.
Use a commuting example to verify exactness, a small noncommuting example with a known solution to verify correctness, and a refinement sequence r, 2r, 4r to inspect convergence. Record the input state and observable rather than reporting an unexplained list of decimal outputs.
Then ask a transfer question: does the method still behave correctly when the initial state changes, the interaction strength changes or the order of the split factors is reversed? These changes expose errors hidden by one convenient example. Verification should probe the structure, not merely reproduce the first successful run.
21. Worked practice: twelve questions with explanations
1. Why is δ = t/r rather than r/t?
There are r equal intervals whose total duration is t, so rδ = t. Therefore δ = t/r. The units are time. The expression r/t has inverse-time units and cannot be substituted for a time interval inside a dimensionless product Hδ when H already has inverse-time units.
2. Compute the leading error for A = X and B = Z
Since [X,Z] = −2iY, the leading split-minus-exact error is −δ²(−2iY)/2 = iδ²Y. Its leading norm is δ². This is consistent with the local bound δ²‖[X,Z]‖/2 = δ². The sign depends on the stated ordering of the split factors.
3. What does reversing A and B change?
The leading commutator changes from [A,B] to [B,A] = −[A,B]. Thus the leading signed error changes direction, while its norm bound stays the same. Reversing the sequence does not generally make the first-order formula exact. It gives a different approximation with the same nominal order.
4. Choose r when t = 2, ‖[A,B]‖ = 3 and ε = 0.02
The first-order global bound is t²‖[A,B]‖/(2r) = 4×3/(2r) = 6/r. Requiring 6/r ≤ 0.02 gives r ≥ 300. This is a sufficient count under the stated bounded-Hermitian assumptions, not an exact prediction of the smallest successful r.
5. A second-order method has error 0.008 at r = 20. What should r = 40 suggest?
In the asymptotic regime, error scales as 1/r², so the expected value is about 0.008/4 = 0.002. Treat this as a convergence prediction, not a certificate. Another error source, a coarse step or numerical saturation can prevent the observed ratio from being four.
6. How many elementary exponentials remain after merging r = 10 symmetric steps?
There are ten B factors and eleven A factors, giving twenty-one exponentials. The boundary A factors are half steps; interior A factors are full steps. This count assumes adjacent A evolutions can actually be merged and says nothing by itself about the physical gate cost of each factor.
7. Implement exp(−i0.3Z) using a standard Z rotation
Because RZ(θ) = exp(−iθZ/2), choose θ = 0.6 radians. Using 0.3 as the rotation angle would implement only half the intended evolution. This is a convention error, not a Trotter error, and refining the number of steps will not repair it.
8. Is a norm-preserving simulation necessarily correct?
No. For example, I preserves every norm, but it is not generally equal to exp(−iHt) at nonzero time. Likewise, every product formula in our comparison is unitary, yet their errors differ. Norm preservation is an invariant worth checking, but it is not a substitute for comparison with the intended dynamics.
9. What observable error follows from ε = 0.003 and ‖O‖ = 2?
The derived bound gives 2‖O‖ε = 2×2×0.003 = 0.012. The units are those of O. The estimate is conservative and uniform over normalised input states because we began from an operator-norm bound rather than a single-state test.
10. Why is the fourth-order middle coefficient negative?
The composition coefficients are chosen to preserve the total time while cancelling the cubic error term. With p = 1/(4−41/3), the central coefficient 1−4p is negative. Its implementation requires inverse elementary unitary evolutions, not an inverse of whatever uncontrolled noise happens during the experiment.
11. Why can more steps worsen an experimental answer?
More steps can reduce ideal splitting error while adding physical gates, decoherence exposure and calibration-sensitive operations. The total experimental discrepancy reflects all these contributions. An optimal practical step count therefore depends on the device and error budget, not just the mathematical convergence order.
12. What should be checked before comparing two simulation methods?
Use the same Hamiltonian, units, time interval, input task and accuracy criterion. Include state preparation, gate synthesis, measurement and classical post-processing where relevant. Otherwise an apparent advantage may come from comparing different problems or from leaving an expensive part of one method outside the accounting.
22. The learning bridge back to familiar Mathematics
This topic recombines ideas a learner may already recognise. Algebra preserves the distinction between AB and BA. Taylor series locate the first term where two constructions disagree. Trigonometry solves exponentials of Pauli operators. Inequalities convert local errors into global guarantees. Numerical tables test whether an asymptotic prediction is visible in finite calculations.
A useful teaching sequence is to solve the commuting case first, compute a nonzero commutator second, build the first-order formula third and only then introduce higher order. Starting with a sophisticated circuit diagram can hide the reason the algorithm exists. Starting with the failed scalar exponential identity makes the need for controlled approximation visible.
The durable lesson is broader than quantum computation: preserve the right invariants, distinguish an approximation from its target, and attach a meaningful error statement to the result. A method is ready to use when the learner can explain not only how it works but also when its guarantee stops applying.
Sources and further study
[1] Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe and Shuchen Zhu, A Theory of Trotter Error, published as Theory of Trotter Error with Commutator Scaling. The source develops commutator-sensitive product-formula analysis beyond the elementary two-term derivation given here.
[2] IBM Quantum Learning, Quantum simulation. A practical treatment of spin-Hamiltonian decomposition, first- and second-order product formulas and comparison with exact dynamics.
[3] Chi-Fang Chen and Fernando G. S. L. Brandão, Average-Case Speedup for Product Formulas. The preliminary product-formula definitions include the standard even-order Suzuki recursion; average-case results should not be confused with the worst-case norm criterion used in this guide.
[4] Dong An, Di Fang and Lin Lin, Time-dependent unbounded Hamiltonian simulation with vector norm scaling. This source addresses why time dependence, unbounded operators and state-dependent bounds require additional care.
Continue through Quantum Mathematics
Guide 14: Variational Quantum Eigensolvers, Rayleigh Quotients and Optimisation asks how to approximate a low-energy state rather than a time evolution. Guide 15: Quantum Metrology, Fisher Information and Precision Limits asks how dynamics encode a measurable parameter. Guide 16: Tensor Networks, Matrix Product States and Entanglement Compression asks when the resulting states admit an efficient classical representation.