9 Lecture 9: Hamiltonian Simulation and Trotterization
Class date: October 12, 2026
Hamiltonian simulation asks us to implement the time evolution generated by a Hamiltonian \(H\). When \(H\) is a sum of simple terms, each term may be easy to simulate even though their sum is not. Product formulas convert this decomposition into a quantum circuit. Their accuracy is controlled by the noncommutativity of the terms, while their cost is determined by the number and order of the time steps.
9.1 The Hamiltonian-simulation problem
For a time-independent Hamiltonian \(H\), the Schrödinger equation
\[ i\frac{d}{dt}\lvert\psi(t)\rangle=H\lvert\psi(t)\rangle \tag{9.1}\]
has solution
\[ \lvert\psi(t)\rangle=U(t)\lvert\psi(0)\rangle, \qquad U(t)=e^{-iHt}. \tag{9.2}\]
The same unitary evolves a density operator according to
\[ \rho(t)=U(t)\rho(0)U^{\dagger}(t). \tag{9.3}\]
Hamiltonian dynamics appears directly in quantum chemistry and many-body physics, where it describes spectra, transport, reactions, phase transitions, and thermalization. It also appears as a subroutine in other quantum algorithms. For example, phase estimation applied to \(e^{-iHt}\) can reveal eigenvalues of \(H\).
The algorithmic task is to construct a circuit for a unitary \(\widetilde U(t)\) satisfying
\[ \left\lVert \widetilde U(t)-e^{-iHt}\right\rVert\leq\epsilon, \tag{9.4}\]
where \(\lVert\cdot\rVert\) denotes the operator norm. This guarantee controls the error for every normalized input state because
\[ \left\lVert\left(\widetilde U(t)-e^{-iHt}\right)\lvert\psi\rangle\right\rVert_2 \leq \left\lVert\widetilde U(t)-e^{-iHt}\right\rVert. \tag{9.5}\]
9.2 Decomposing the Hamiltonian
Suppose
\[ H=\sum_{j=1}^{m}H_j, \tag{9.6}\]
where each \(H_j\) is Hermitian and \(e^{-iH_j\tau}\) can be implemented efficiently. Local Hamiltonians and Hamiltonians expanded as sums of Pauli strings naturally have this form.
If all terms commute pairwise, then
\[ e^{-iHt} =e^{-i(\sum_{j=1}^{m}H_j)t} =\prod_{j=1}^{m}e^{-iH_jt}. \tag{9.7}\]
In this case the decomposition is exact. For noncommuting terms, however,
\[ [H_j,H_k]\neq0 \quad\Longrightarrow\quad e^{-i(H_j+H_k)t}\neq e^{-iH_jt}e^{-iH_kt} \tag{9.8}\]
in general. Product formulas approximate the evolution by alternating short evolutions under the individual terms.
9.3 First-order product formula
Divide the total time into \(r\) steps of duration
\[ \Delta=\frac{t}{r}. \tag{9.9}\]
One first-order step is
\[ S_1(\Delta) =e^{-iH_1\Delta}e^{-iH_2\Delta}\cdots e^{-iH_m\Delta}. \tag{9.10}\]
Repeating the step gives the Lie–Trotter approximation
\[ e^{-iHt}\approx S_1(t/r)^r. \tag{9.11}\]
Each Trotter step evolves under every component Hamiltonian for a short time. Increasing \(r\) reduces the digital approximation error but increases the circuit depth.
9.3.1 The error comes from noncommutativity
For \(H=A+B\), the Baker–Campbell–Hausdorff expansion gives
\[ e^{-iA\Delta}e^{-iB\Delta} =\exp\left(-i(A+B)\Delta-\frac12[A,B]\Delta^2+O(\Delta^3)\right) \tag{9.12}\]
as \(\Delta\to0\), for bounded \(A\) and \(B\). The leading correction vanishes when \(A\) and \(B\) commute. In operator norm, the one-step error therefore scales as
\[ \left\lVert e^{-i(A+B)\Delta}-e^{-iA\Delta}e^{-iB\Delta}\right\rVert =O\left(\lVert[A,B]\rVert\Delta^2\right). \tag{9.13}\]
The following bound explains how this local error accumulates.
Lemma 9.1 (Error accumulation over repeated steps) For unitaries \(U\) and \(V\),
\[ \lVert U^r-V^r\rVert\leq r\lVert U-V\rVert. \tag{9.14}\]
Proof. Using the telescoping identity
\[ U^r-V^r=\sum_{q=0}^{r-1}U^{r-1-q}(U-V)V^q, \]
unitary invariance of the operator norm and the triangle inequality give Equation 9.14.
For a fixed Hamiltonian decomposition, the first-order global error consequently has the form
\[ \left\lVert e^{-iHt}-S_1(t/r)^r\right\rVert =O\left(\frac{t^2}{r}\sum_{1\leq j<k\leq m}\lVert[H_j,H_k]\rVert\right), \tag{9.15}\]
up to constants depending on the Hamiltonian norms and the chosen ordering. Suppressing those Hamiltonian-dependent constants, achieving error \(\epsilon\) requires
\[ r=O\left(\frac{t^2}{\epsilon}\right). \tag{9.16}\]
9.4 Second-order product formula
For \(H=A+B\), the symmetric or Strang splitting is
\[ S_2(\Delta) =e^{-iA\Delta/2}e^{-iB\Delta}e^{-iA\Delta/2}. \tag{9.17}\]
The full simulation uses
\[ e^{-i(A+B)t}\approx S_2(t/r)^r. \tag{9.18}\]
The symmetry removes the second-order term in the local error. The leading contribution is governed by nested commutators, giving
\[ \left\lVert e^{-i(A+B)\Delta}-S_2(\Delta)\right\rVert =O\left(\left(\lVert[A,[A,B]]\rVert+\lVert[B,[B,A]]\rVert\right)\Delta^3\right). \tag{9.19}\]
After \(r\) steps, the global error scales as
\[ O\left(r\Delta^3\right)=O\left(\frac{t^3}{r^2}\right), \tag{9.20}\]
with the nested-commutator factors understood. Thus the step count needed for fixed error scales as
\[ r=O\left(\frac{t^{3/2}}{\sqrt{\epsilon}}\right) \tag{9.21}\]
when the Hamiltonian-dependent constants are fixed. Second-order splitting often provides a useful balance: it reduces the number of time steps substantially while only modestly increasing the number of exponentials per step.
For \(m\) terms, a symmetric second-order step can be written as
\[ S_2(\Delta) =e^{-iH_1\Delta/2}\cdots e^{-iH_{m-1}\Delta/2}e^{-iH_m\Delta} e^{-iH_{m-1}\Delta/2}\cdots e^{-iH_1\Delta/2}. \tag{9.22}\]
9.5 Higher-order Suzuki formulas
Higher even-order formulas are obtained recursively from symmetric lower-order formulas. For \(k\geq2\), define
\[ S_{2k}(\Delta) =S_{2k-2}(p_k\Delta)^2 S_{2k-2}((1-4p_k)\Delta) S_{2k-2}(p_k\Delta)^2, \tag{9.23}\]
where
\[ p_k=\frac{1}{4-4^{1/(2k-1)}}. \tag{9.24}\]
The coefficient \(1-4p_k\) is negative, so a higher-order step contains both forward and backward evolutions under the component Hamiltonians. This composition cancels all error terms through order \(2k\), leaving local error
\[ O(\Delta^{2k+1}). \tag{9.25}\]
After \(r\) steps, the corresponding global error has the scaling
\[ O\left(\frac{t^{2k+1}}{r^{2k}}\right). \tag{9.26}\]
More generally, a product formula of order \(p\) requires
\[ r=O\left(\frac{t^{1+1/p}}{\epsilon^{1/p}}\right) \tag{9.27}\]
when Hamiltonian-dependent constants are suppressed. As \(p\) increases, the exponent \(1+1/p\) approaches \(1\), and the dependence on simulation time approaches linear scaling. This does not make arbitrarily high order automatically optimal: the number of exponentials within one Suzuki step grows recursively, and the commutator-dependent constants also affect the best choice of order.
9.5.1 Why linear dependence on time is a natural limit
For general sparse Hamiltonians accessed through an oracle, the no-fast-forwarding theorem gives a worst-case lower bound proportional to the simulated time:
\[ \text{query complexity}=\Omega(t) \tag{9.28}\]
for fixed constant accuracy and a normalized Hamiltonian family. Thus the nearly linear time dependence approached by higher-order formulas is close to the best possible generic scaling. The statement is a worst-case result; special Hamiltonians with additional structure can sometimes be simulated more rapidly.
9.6 Randomized product formulas
Standard Trotterization uses the same ordering of the terms in every step. For
\[ H=\sum_{\ell=1}^{L}H_\ell, \tag{9.29}\]
choose a permutation \(\sigma\) of \(\{1,\ldots,L\}\) and define
\[ W_1(\Delta;\sigma) =\prod_{\ell=1}^{L}e^{-i\Delta H_{\sigma(\ell)}}. \tag{9.30}\]
A symmetric step is
\[ W_2(\Delta;\sigma) =\left(\prod_{\ell=1}^{L}e^{-i\Delta H_{\sigma(\ell)}/2}\right) \left(\prod_{\ell=L}^{1}e^{-i\Delta H_{\sigma(\ell)}/2}\right). \tag{9.31}\]
A randomized product formula samples a new ordering at each time step. Different orderings produce different leading error operators. With a fixed ordering, these coherent errors may repeatedly align; randomization prevents the same alignment from occurring systematically. Averaging over the sampled orderings can therefore improve the error bound, particularly its dependence on the number of Hamiltonian terms. The improvement is an expected or averaged statement and does not imply exact cancellation in every individual realization.
9.7 Implementing Pauli-string evolutions
Suppose a Hamiltonian term is a Pauli string
\[ P=\sigma_1\otimes\sigma_2\otimes\cdots\otimes\sigma_n, \qquad \sigma_j\in\{I,X,Y,Z\}. \tag{9.32}\]
Because \(P^2=I\), its exponential has the closed form
\[ e^{-i\theta P}=\cos(\theta)I-i\sin(\theta)P. \tag{9.33}\]
The corresponding circuit is constructed in four steps:
- Change basis on every nonidentity factor so that each \(X\) or \(Y\) becomes \(Z\). Use \(\mathrm{H}\) for \(X\), and apply \(S^{\dagger}\) followed by \(\mathrm{H}\) for \(Y\).
- Use a controlled-NOT ladder to compute the parity of all supported qubits onto one target qubit.
- Apply \(\mathrm{R}_z(2\theta)\) to that target, using the convention \(\mathrm{R}_z(\phi)=e^{-i\phi Z/2}\).
- Reverse the controlled-NOT ladder and the basis changes.
The two-qubit \(ZZ\) interaction is the basic example.
9.8 Example: transverse-field Ising model
Consider the one-dimensional transverse-field Ising Hamiltonian
\[ H=-J\sum_{j=1}^{n-1}Z_jZ_{j+1}-h\sum_{j=1}^{n}X_j. \tag{9.34}\]
Split it into
\[ H_Z=-J\sum_{j=1}^{n-1}Z_jZ_{j+1}, \qquad H_X=-h\sum_{j=1}^{n}X_j. \tag{9.35}\]
All terms within \(H_Z\) commute, and all terms within \(H_X\) commute. Their separate evolutions therefore factorize exactly:
\[ e^{-iH_Z\tau}=\prod_{j=1}^{n-1}e^{iJ\tau Z_jZ_{j+1}}, \qquad e^{-iH_X\tau}=\prod_{j=1}^{n}e^{ih\tau X_j}. \tag{9.36}\]
The two blocks do not commute with each other. A second-order step is
\[ e^{-iH\Delta} \approx e^{-iH_Z\Delta/2}e^{-iH_X\Delta}e^{-iH_Z\Delta/2}. \tag{9.37}\]
Each \(ZZ\) factor is implemented using Figure 9.1, with the sign incorporated into the rotation angle, and each \(X\) factor is a single-qubit rotation. Under the convention \(\mathrm{R}_x(\phi)=e^{-i\phi X/2}\),
\[ e^{ih\tau X_j}=\mathrm{R}_{x,j}(-2h\tau). \tag{9.38}\]
For the two-qubit Hamiltonian
\[ H=-JZ_1Z_2-h(X_1+X_2), \tag{9.39}\]
a first-order step may be written as
\[ e^{-iH\Delta} \approx e^{iJ\Delta Z_1Z_2}e^{ih\Delta X_1}e^{ih\Delta X_2}. \tag{9.40}\]
Thus one time step becomes a short sequence of controlled-NOT gates and single-qubit rotations, and the same sequence is repeated \(r\) times.
9.9 Accuracy and implementation cost
For fixed total time, reducing \(\Delta=t/r\) reduces the product-formula error. It also increases the number of gates. The leading scalings are
\[ \begin{aligned} \text{first order:}&\qquad O(t^2/r),\\ \text{second order:}&\qquad O(t^3/r^2), \end{aligned} \tag{9.41}\]
with commutator-dependent constants suppressed. On an ideal fault-tolerant computer, one selects the formula order and step count to minimize logical resources at the target accuracy. On a noisy device, increasing \(r\) also creates more opportunities for gate error and decoherence. The best observed accuracy can therefore occur at a finite step count, where digital simulation error and hardware error are balanced.
Trotterization is direct and broadly applicable, but it has three main limitations:
- strongly noncommuting terms can produce large error constants;
- long-time or high-precision simulation can require deep circuits; and
- a decomposition that is mathematically convenient may not be optimal for the available hardware.
These limitations motivate simulation methods based on truncated Taylor series and linear combinations of unitaries, as well as qubitization and problem-tailored analog or digital constructions.
9.10 Summary
Hamiltonian simulation reduces the continuous evolution \(e^{-iHt}\) to gates implementing short evolutions under simpler terms. Commuting terms factorize exactly, while noncommuting terms generate product-formula errors controlled by commutators. First-order, symmetric second-order, and higher-order Suzuki formulas trade longer steps for lower approximation error, and randomized orderings can reduce coherent accumulation. For Pauli Hamiltonians, each short evolution compiles into basis changes, a parity computation, and a single-qubit rotation. The useful step count is ultimately set by the balance between digital approximation error and implementation cost.