Chapter 33: Digital Quantum Simulation

33.1 Hamiltonian Simulation: The Core Problem

The central problem of digital quantum simulation is this: given a Hamiltonian $H$ acting on $n$ qubits, an initial state $|\psi(0)\rangle$, and a time $t$, implement the unitary operator $e^{-iHt}$ as a quantum circuit. The output state $|\psi(t)\rangle = e^{-iHt} |\psi(0)\rangle$ is the solution to the Schrodinger equation:

The Hamiltonian $H$ is a $2^n \times 2^n$ Hermitian matrix. Exponentiating it directly would require manipulating this exponentially large matrix - exactly the classical bottleneck we are trying to avoid. The key insight is that physically relevant Hamiltonians have structure that can be exploited.

Local Hamiltonians

Most physical Hamiltonians are local: they decompose as a sum of terms, each acting on a small number of qubits:

where each $H_j$ acts non-trivially on at most $k$ qubits (typically $k = 2$ for pairwise interactions). Examples include:

• Ising model: $H = J\sum_{\langle i,j \rangle} Z_i Z_j + h\sum_i X_i$. Each term acts on at most 2 qubits.
• Heisenberg model: $H = J\sum_{\langle i,j \rangle} (X_i X_j + Y_i Y_j + Z_i Z_j)$. Pairwise interactions between neighboring spins.
• Electronic structure Hamiltonians: After second quantization and a mapping to qubits (Jordan-Wigner or Bravyi-Kitaev), molecular Hamiltonians become sums of Pauli tensor products.

The locality condition is crucial: while $e^{-iHt}$ is an exponentially large matrix, each $e^{-iH_j t}$ acts on only a few qubits and can be implemented with $O(1)$ quantum gates. The challenge is that $e^{-i(H_1 + H_2)t} \neq e^{-iH_1 t} e^{-iH_2 t}$ in general, because $H_1$ and $H_2$ typically do not commute.

Complexity of Hamiltonian Simulation

For a Hamiltonian with $m$ local terms, the simulation cost depends on:

• $t$: the simulation time (longer times require more work).
• $m$: the number of terms in the Hamiltonian.
• $\|H\|$: the norm of the Hamiltonian (stronger interactions are harder).
• $\epsilon$: the desired accuracy (tighter precision costs more).

Different simulation methods achieve different trade-offs among these parameters. The theoretical optimum is $O(\alpha t + \log(1/\epsilon))$ queries, where $\alpha$ characterizes the access model (e.g., the 1-norm of coefficients). This optimal scaling is achieved by qubitization and QSVT (Section 33.3).

33.2 Trotterization

The Trotter-Suzuki decomposition is the simplest and most widely used method for Hamiltonian simulation. It is named after Hale Trotter (1959) and Masuo Suzuki (1976, 1990), who developed the mathematical theory of product formulas for operator exponentials. Suzuki's 1976 work generalized the Trotter product formula, while his 1990 paper introduced the recursive fractal decomposition that yields higher-order formulas.

First-Order Trotter Formula

The idea is straightforward. Although $e^{-i(A+B)t} \neq e^{-iAt} e^{-iBt}$ when $[A,B] \neq 0$, the Trotter formula says that the product becomes a good approximation when we divide time into small steps:

For finite $n$ (the number of Trotter steps), the error is:

The error scales as $O(t^2/n)$ for two terms. For a Hamiltonian with $m$ terms, $H = \sum_{j=1}^m H_j$, the error per step generalizes to $O(m^2 \|H\|^2 t^2/n)$. One Trotter step applies all $m$ individual evolutions, so to achieve accuracy $\epsilon$ the total gate count is $O(m^3 \|H\|^2 t^2/\epsilon)$.

Second-Order Suzuki-Trotter Formula

A better approximation symmetrizes the product. The second-order Suzuki-Trotter formula is:

This "palindrome" structure (forward then backward) cancels the leading-order error term, giving:

The error now scales as $O(t^3/n^2)$ - a significant improvement. To achieve accuracy $\epsilon$, we need $n = O(t^{3/2}/\sqrt{\epsilon})$ steps, which is fewer than the first-order formula requires.

Higher-Order Formulas

Suzuki showed how to recursively construct product formulas of any even order $2k$:

where $p_k = (4 - 4^{1/(2k-1)})^{-1}$. The $2k$-th order formula has error $O(t^{2k+1}/n^{2k})$. Higher orders converge faster with $n$ but require more gate operations per step (the formulas become longer). In practice, the optimal order depends on the specific Hamiltonian and desired accuracy.

Implementing Trotter Steps for the Ising Model

Let us implement Trotterization for the transverse-field Ising model on two qubits:

One Trotter step requires:

1. $e^{-iJ Z_0 Z_1 \Delta t}$: This is the $\text{Rzz}(2J\Delta t)$ gate, which applies a $ZZ$ rotation. In QASM: rzz(2*J*dt) q[0], q[1].
2. $e^{-ihX_0 \Delta t}$: This is the $\text{Rx}(2h\Delta t)$ gate on qubit 0. In QASM: rx(2*h*dt) q[0].
3. $e^{-ihX_1 \Delta t}$: Same rotation on qubit 1. In QASM: rx(2*h*dt) q[1].

Try the sandbox below. It simulates the Ising model with $J = 1$, $h = 0.5$, total time $t = \pi/2$, using a variable number of Trotter steps. The code uses 4 steps (each with $\Delta t = \pi/8$).

The circuit above starts with qubit 0 in state $|1\rangle$ and qubit 1 in state $|0\rangle$, then simulates the Ising Hamiltonian for time $t = \pi/2$. The $ZZ$ interaction causes energy exchange between the qubits, while the transverse field $X$ induces spin flips. You should see a non-trivial probability distribution over all four basis states, reflecting the quantum dynamics.

Try modifying the code: reduce to 2 Trotter steps (double the angles per step) or increase to 8 steps (halve the angles). With more steps, the simulation becomes more accurate, but the circuit becomes deeper. This trade-off is at the heart of Trotterization.

The Heisenberg Model: Three Interaction Terms

The Heisenberg model $H = J(X_i X_j + Y_i Y_j + Z_i Z_j)$ requires three interaction gates per bond per Trotter step:

• $e^{-iJ X_i X_j \Delta t}$: the $\text{Rxx}(2J\Delta t)$ gate
• $e^{-iJ Y_i Y_j \Delta t}$: the $\text{Ryy}(2J\Delta t)$ gate
• $e^{-iJ Z_i Z_j \Delta t}$: the $\text{Rzz}(2J\Delta t)$ gate

This sandbox simulates the Heisenberg spin chain: starting with qubit 0 in $|1\rangle$ and qubit 1 in $|0\rangle$, the $XX + YY + ZZ$ interaction causes the spin excitation to hop between sites. The Heisenberg model conserves total spin, so you should see significant probability only on $|10\rangle$ and $|01\rangle$ (the one-excitation sector), with very little probability on $|00\rangle$ or $|11\rangle$ (which would indicate Trotter error).

Trotter Error in Practice

The theoretical error bounds for Trotterization are often pessimistic. Childs, Su, Tran, Wiebe, and Zhu (2021) developed a refined "commutator scaling" theory that exploits the commutativity structure of the Hamiltonian terms. Their key insight: the Trotter error depends not on the norms of individual terms but on the norms of their commutators. If many terms commute (or nearly commute), the actual error is much smaller than the worst-case bound predicts.

For instance, in a lattice Hamiltonian where terms on distant sites commute exactly, the commutator-based bound gives errors that are tight to within a factor of approximately 5, compared to the naive bound which can overestimate by orders of magnitude.

33.3 Advanced Techniques: Qubitization and QSP

While Trotterization is intuitive and practical, it is not asymptotically optimal. The gate count for first-order Trotter scales as $O(m^2 t^2/\epsilon)$, which has suboptimal dependence on both $t$ and $\epsilon$. Two advanced techniques achieve the theoretically optimal scaling: qubitization (Low and Chuang, 2019) and quantum signal processing (Low and Chuang, 2017), unified under the QSVT framework (Chapter 31).

Linear Combinations of Unitaries (LCU)

Before describing qubitization, we need the LCU lemma (Childs and Wiebe, 2012). Suppose the Hamiltonian can be written as a linear combination of unitaries:

where each $U_j$ is a unitary operator. For Pauli Hamiltonians (which arise from quantum chemistry), this decomposition is natural: each $U_j$ is a tensor product of Pauli operators, and $\alpha_j$ is the corresponding coefficient.

The LCU technique uses two operations:

• PREPARE: A unitary that maps $|0\rangle$ to $\sum_j \sqrt{\alpha_j/\lambda} |j\rangle$, where $\lambda = \sum_j \alpha_j$ is the 1-norm.
• SELECT: A controlled unitary that maps $|j\rangle|\psi\rangle$ to $|j\rangle U_j |\psi\rangle$.

Using these, one can implement a block encoding of $H/\lambda$: the unitary $U = \text{PREPARE}^\dagger \cdot \text{SELECT} \cdot \text{PREPARE}$ satisfies $\langle 0 | U | 0 \rangle = H/\lambda$ (projecting the ancilla register). This block encoding is the gateway to qubitization and QSVT.

Qubitization

Qubitization, developed by Low and Chuang (2019), converts a block encoding of $H/\lambda$ into a quantum walk operator whose eigenvalues encode the eigenvalues of $H$. Specifically, the "qubitized walk operator" $W$ has eigenvalues $e^{\pm i \arccos(\lambda_k/\lambda)}$, where $\lambda_k$ are the eigenvalues of $H$.

The key property of qubitization is that it embeds the Hamiltonian's eigenvalues into an $SU(2)$ subspace - a two-dimensional space for each eigenvalue. This embedding is exactly what quantum signal processing needs: by applying a sequence of signal processing rotations to $W$, one can implement any polynomial transformation of the eigenvalues.

Achieving Optimal Simulation

To simulate $e^{-iHt}$ using qubitization and QSP:

1. Construct a block encoding of $H/\lambda$ using PREPARE and SELECT oracles.
2. Build the qubitized walk operator $W$.
3. Find QSP angles $\phi_0, \ldots, \phi_d$ such that the polynomial $P(x) \approx e^{-i\lambda t \cdot x}$ (a polynomial approximation to the complex exponential).
4. Apply the QSP sequence: alternating $W$ with single-qubit rotations determined by the angles $\phi_k$.

The degree of the polynomial $P$ needed to approximate $e^{-i\lambda t x}$ to accuracy $\epsilon$ is $d = O(\lambda t + \log(1/\epsilon)/\log\log(1/\epsilon))$. Each degree requires one use of $W$ (which uses one PREPARE and one SELECT), giving the total query complexity.

Comparison of Simulation Methods

The progression from Trotterization to qubitization represents a journey from $O(t^2/\epsilon)$ to $O(t + \log(1/\epsilon))$ - a dramatic improvement in asymptotic complexity. However, the constant factors and ancilla requirements of advanced methods can be substantial, and for many practical instances (moderate $t$, modest $\epsilon$), Trotterization with optimized ordering remains competitive.

33.4 Simulating Chemistry

Quantum chemistry is the flagship application of quantum simulation. The goal is to compute the electronic structure of molecules: the energies, wavefunctions, and properties of electrons in the Coulomb potential of atomic nuclei. This is the foundation of all of chemistry, materials science, and much of biology.

The Electronic Structure Problem

Under the Born-Oppenheimer approximation (nuclei are stationary), the electronic Hamiltonian for a molecule is:

The three terms represent: kinetic energy of electrons, electron-nuclear attraction, and electron-electron repulsion. It is the electron-electron repulsion (the last term) that makes the problem hard: it couples all electrons together, creating the "correlation" that classical methods struggle to capture.

Second Quantization and Qubit Mapping

To simulate this on a quantum computer, we first express $H_{\text{elec}}$ in second quantization using a finite basis of molecular orbitals:

where $a_p^\dagger$ and $a_p$ are fermionic creation and annihilation operators, and $h_{pq}$, $h_{pqrs}$ are one- and two-electron integrals computed classically. The operators satisfy anticommutation relations $\{a_p, a_q^\dagger\} = \delta_{pq}$.

We then map fermionic operators to qubit operators using one of several encodings:

• Jordan-Wigner transformation: Maps each orbital to one qubit. The occupation $n_p = a_p^\dagger a_p$ becomes $(I - Z_p)/2$. Creation operators become strings of $Z$ operators (to enforce antisymmetry). The result is a Hamiltonian expressed as a sum of Pauli strings: $H = \sum_k \alpha_k P_k$, where each $P_k$ is a tensor product of Pauli operators.
• Bravyi-Kitaev transformation: An alternative mapping that reduces the length of Pauli strings from $O(n)$ to $O(\log n)$, at the cost of more complex encoding.

For a molecule described by $M$ molecular orbitals, the qubit Hamiltonian has $M$ qubits and $O(M^4)$ Pauli terms (from the two-electron integrals). This is the Hamiltonian we simulate.

A Simple Example: Hydrogen Molecule

The simplest non-trivial molecule is H$_2$. In a minimal basis (STO-3G), H$_2$ requires 4 spin-orbitals and thus 4 qubits. After applying symmetries and simplifications, the qubit Hamiltonian in the Jordan-Wigner encoding reduces to:

where the coefficients $g_k$ depend on the internuclear distance. The $XX$ and $YY$ terms come from electron hopping, while the $ZZ$ term comes from electron-electron repulsion. We can simulate this with Trotter steps using the gates we already know:

This sandbox simulates a simplified H$_2$ Hamiltonian starting from the Hartree-Fock state $|01\rangle$. The $ZZ$ term represents electron-electron interaction, while $XX$ and $YY$ terms represent electron hopping (correlation). Run the circuit and observe the output: the presence of $|10\rangle$ probability indicates electron correlation beyond the mean-field (Hartree-Fock) approximation. This is precisely the quantum advantage - capturing correlation effects that classical mean-field methods miss.

Resource Estimates for Real Molecules

For industrially relevant molecules, the resource requirements are substantial:

The FeMoCo active site of nitrogenase - the enzyme responsible for biological nitrogen fixation - is a prime target for quantum simulation. Understanding how FeMoCo catalyzes the conversion of $\text{N}_2$ to $\text{NH}_3$ at ambient temperature and pressure (a reaction that industrially requires extreme conditions in the Haber-Bosch process) could revolutionize fertilizer production and reduce its enormous carbon footprint.

33.5 Simulating Materials and Condensed Matter

Beyond molecules, quantum simulation can address fundamental questions in condensed matter physics - the study of solids, liquids, and other many-body systems. Many of the most important open questions in physics involve quantum many-body systems that are beyond the reach of classical simulation.

The Hubbard Model

The Hubbard model is perhaps the most important unsolved model in condensed matter physics:

Despite its simplicity (only two parameters, $t$ and $U$), the 2D Hubbard model is believed to capture the essential physics of high-temperature superconductivity in cuprate materials. Yet its phase diagram remains unknown in important parameter regimes. Classical methods - quantum Monte Carlo, DMRG, dynamical mean-field theory - each have limitations (the sign problem, dimensionality restrictions, or mean-field approximations) that prevent a definitive solution.

A quantum computer could simulate the Hubbard model directly. Using the Jordan-Wigner transformation on a 2D lattice, the Hamiltonian maps to a qubit Hamiltonian with nearest-neighbor hopping terms (long Pauli strings due to the Jordan-Wigner encoding) and on-site interaction terms ($ZZ$ type). Trotterization then proceeds as before, with $\text{Rxx}$, $\text{Ryy}$ gates for hopping and $\text{Rzz}$ gates for interactions.

Topological Phases

Some of the most exotic quantum phases of matter are topological phases - states characterized not by local order parameters (like magnetization) but by global topological invariants. Examples include:

• Topological insulators: Materials that are insulating in their bulk but conduct on their surface through topologically protected edge states.
• Fractional quantum Hall states: Two-dimensional electron gases in strong magnetic fields that exhibit fractionalized excitations (anyons) with exotic exchange statistics.
• Quantum spin liquids: Magnetic materials where quantum fluctuations prevent ordering even at zero temperature, potentially hosting anyonic excitations useful for topological quantum computing.

Quantum simulation can probe these phases by preparing ground states of topological Hamiltonians and measuring topological invariants (e.g., entanglement entropy, string order parameters). This is an area where analog quantum simulators (particularly Rydberg atom arrays) have already made significant experimental progress.

Lattice Gauge Theories

Quantum field theories - the theoretical framework of particle physics - can be discretized on a lattice for numerical simulation. Lattice gauge theories describe the strong nuclear force (quantum chromodynamics) and electroweak interactions. Classical lattice Monte Carlo methods are powerful for equilibrium properties but fail for real-time dynamics (due to the sign problem) and for systems at finite density (relevant to neutron stars and heavy-ion collisions).

Quantum simulation of lattice gauge theories is a rapidly growing field. The Hamiltonian formulation of lattice gauge theories maps naturally to qubit systems (with gauge constraints implemented via energy penalties or symmetry-preserving encodings). Early quantum simulations have demonstrated proof-of-concept for 1+1 dimensional models (Schwinger model, lattice QED), and scaling to higher dimensions is an active area of research.

Near-Term Prospects

What can we realistically expect from quantum simulation in the near term? The answer depends on the hardware:

• NISQ era (now - ~2030): Noisy quantum computers with 50-1000 qubits and limited circuit depth. Variational algorithms (VQE, QAOA) and shallow Trotter circuits can address small molecules and lattice models. Results must be carefully benchmarked against classical methods. Some scientific insights (e.g., qualitative phase diagrams) may be achievable even without full error correction.
• Early fault-tolerant era (~2030+): Logical qubits with low error rates enabling deeper circuits. First applications where quantum simulation provides unambiguous advantages over all classical methods - likely in quantum chemistry (molecules with 50-100 orbitals) or small lattice models.
• Full-scale fault-tolerant era: Thousands of logical qubits with arbitrary circuit depth. Industrial-scale quantum chemistry, materials design, and simulation of quantum field theories. This is when the vision of Feynman's 1982 proposal is fully realized.

Interactive: Hamiltonian Expectation Value

A key output of digital quantum simulation is the expectation value of an observable on the time-evolved state. Adjust the evolution parameter below to see how the $\langle ZZ \rangle$ expectation value changes as the simulated Ising interaction evolves the system.