Simulating Quantum Systems

Write down the state of $n$ interacting spins. By the tensor product of Chapter 2 it lives in ${\mathbb{C}}^{2^n}$: a list of $2^n$ complex amplitudes. Fifty spins is already $2^{50}\approx 10^{15}$ numbers — sixteen petabytes of amplitudes for a system you could hold on a fingertip. To evolve that state forward in time a classical computer must store and update every one of them. This is not a failure of cleverness; it is the dimension of the space. Richard Feynman's observation in 1982 was that the difficulty is physical: nature evolves those $2^n$ amplitudes effortlessly, so perhaps a computer built out of quantum parts could too.

The evolution in question is fixed by a single operator. A closed quantum system with Hamiltonian $H$ (a Hermitian matrix, the energy observable) evolves by the Schrödinger equation, whose solution is the one-parameter family of unitaries

That $U(t)$ is unitary because $H$ is Hermitian — the matrix exponential of $-iHt$ preserves length exactly, as a legal gate must. So “simulate the system” means apply the gate $e^{-iHt}$. The catch is that $H$ is a $2^n\times 2^n$ matrix; we cannot form its exponential directly on a real machine any more than we could form the state. We need to build $e^{-iHt}$ out of gates that each touch only a few qubits — and that is exactly what the structure of physical Hamiltonians lets us do.

Real Hamiltonians are local: they are sums of terms, each acting on only a few qubits,

If the pieces all commuted, we would be done — $e^{-iHt}$ would factor into ${\prod }_j{e}^{-i{H}_{j}t}$, and each $e^{-i{H}_{j}t}$ is a small gate we know how to build. But they do not commute. The exponential of a sum equals the product of exponentials only when the exponents commute; in general the Baker–Campbell–Hausdorff formula leaves a residue,

so $e^{-i(A+B)}\ne e^{-iA}{e}^{-iB}$ whenever the commutator $[A,B]$ is nonzero. The Trotter–Suzuki idea is to make that residue small by taking small steps. Slice the evolution into $r$ pieces of duration $t/r$ and approximate each slice by the product formula. For two terms $H=A+B$, the first-order (Lie–Trotter) formula is

with a total error of order $t^2/r$: each of the $r$ steps carries an $O((t/r{)}^2)$ error from the discarded $\frac{1}{2}[A,B](t/r{)}^2$ term, and the errors add. Halving the step size halves the error. We can do better for free by symmetrizing the step. The second-order (Strang) formula splits one term in half around the other,

and the symmetry cancels the leading error term exactly (any symmetric ordering gives second order; the lab, matching (9.5) with $A$ the interaction layer and $B$ the field layer of §9.3, symmetrizes the field layer), leaving a total error of order $t^3/r^2$: now doubling the steps cuts the error by four. These two error orders — error $\propto 1/r$ for first order, $\propto 1/r^2$ for second — are the central quantitative claim of the chapter, and §9.5 verifies both as clean slopes on a log-log plot. The generalization to many non-commuting terms, and to arbitrarily high order, is due to Masuo Suzuki; the price of order $2k$ is more sub-steps per slice, a trade the practitioner tunes.

Make it concrete on the workhorse model of quantum magnetism, the transverse-field Ising model on a chain of $n$ qubits:

The first sum is an interaction: neighbouring spins prefer to align along $Z$ (coupling $J$). The second is a transverse field of strength $h$ that tips each spin toward the $X$ axis and so fights the alignment. The two pieces do not commute — $[Z_i{Z}_{i+1},X_i]\ne 0$ — which is what makes the dynamics rich and the naive product formula wrong. Group the terms into two layers: $A=-J{\sum }_i{Z}_i{Z}_{i+1}$ and $B=-h{\sum }_i{X}_i$. Within a layer the terms commute (all the $Z_i{Z}_{i+1}$ are simultaneously diagonal; all the $X_i$ act on different qubits), so each layer's exponential factors exactly into local gates.

And those local gates are elementary. A single $ZZ$ term is diagonal with eigenvalues $\pm 1$, so its exponential is just a phase stamped on each basis state: $e^{i\theta Z_i{Z}_j}$ multiplies amplitude $k$ by $e^{+i\theta }$ when bits $i,j$ agree and $e^{-i\theta }$ when they differ. A single $X$ term exponentiates by the identity $e^{i\varphi X}=\cos \varphi I+i\sin \varphi X$, a one-qubit rotation. So one Trotter step is a $ZZ$-phase layer followed by an $X$-rotation layer — a shallow, entirely local circuit.

One bookkeeping note, since “bit $i$ of amplitude $k$” only means something once the labelling is fixed. The lab numbers the chain sites $i=0,\dots ,n-1$ from one end and gives site $i$ bit $i$ of the amplitude index — site $0$ is the least-significant bit, $k={\sum }_i{s}_i{2}^i$. That is the reverse of the ket labelling used in Chapters 6 to 8, where qubit $0$ is written leftmost and is the most-significant bit. Nothing in this chapter turns on the choice: the open chain with uniform $J$ and $h$ is unchanged by the site reversal $i\leftrightarrow n-1-i$, so the two labellings build the identical matrix and $⟨Z_0⟩=⟨Z_{n-1}⟩$ to round-off — a number the lab prints rather than assumes. It is stated here so that no one has to infer it from a shift operator. With the labelling fixed, here is the whole step as a circuit:

The payoff is a real physical prediction: the magnetization $⟨Z_0(t)⟩$ of the first spin as the chain evolves from all-up. The exact curve comes from diagonalizing $H$ (the lab's referee oracle); the Trotter curve is the circuit above, run $r$ times. Slide $r$ up and watch the approximation converge — and switch between first and second order to feel the difference between $1/r$ and $1/r^2$:

At a handful of steps the Trotter curve lags and overshoots the true magnetization — the discarded commutators show up as a visible phase error. A dozen second-order steps already track the exact curve to the width of the line. That is the whole method: trade circuit depth for accuracy, on a schedule you can predict in advance.

Time evolution is not only interesting for watching a system move; it is the gateway to energies. Recall from Chapter 7 that phase estimation reads the eigenphase of a unitary. Feed it the unitary $U=e^{-iHt}$ — built by Trotter — on an eigenstate $|E_k⟩$ of $H$. Since $e^{-iHt}|E_k⟩=e^{-i{E}_{k}t}|E_k⟩$, $|E_k⟩$ is an eigenvector of $U$ and the algorithm returns its phase. Be careful with the convention: (7.4) writes the eigenvalue as $U|u⟩=e^{2\pi i\phi }|u⟩$ with $\phi \in [0,1)$, so the phase register does not hold $E_{k}t$ but

Both the $2\pi$ and the minus sign are real: the sign is the one in $e^{-iHt}$, and the $\mathrm{mod}1$ is a genuine constraint — the evolution time $t$ must be short enough that every $|E_k|t/2\pi$ of interest stays inside one period, or the energies alias, exactly the wraparound Chapter 7 managed with its register-sizing rule. With that fixed, reading ${\phi }_k$ off the phase register yields the energy $E_k$. This composition — Trotter to build the evolution, phase estimation to read its phase — is the textbook quantum algorithm for eigenvalue estimation, and the ground-state energy of a molecule or material is exactly the quantity chemistry and materials science most want. It is the direct ancestor of the variational quantum eigensolver of Chapter 13, which trades the deep, coherent phase-estimation circuit for many shallow ones.

Now the honest accounting. Hamiltonian simulation is the clearest case of an exponential quantum advantage we have: the classical cost is the dimension $2^n$, the quantum cost is polynomial in $n$ and in the evolution time, and — unlike Grover's quadratic edge or the fragile promises of some learning tasks — there is no known classical algorithm that dequantizes it for generic local Hamiltonians. Chemistry and materials are its first real customers. But the advantage is asymptotic, and the useful instances are large: an accurate energy for an industrially interesting molecule needs many well-controlled qubits and deep circuits, which means it needs fault tolerance — the error correction of Part III. Today's noisy machines can simulate small chains like ours beautifully; the molecules that would change an industry are still waiting on the hardware. That gap between the clean theory and the demanding practice is precisely why the next Part exists.

The lab builds the transverse-field Ising chain from scratch and evolves it two independent ways. The exact path diagonalizes the dense $2^n\times 2^n$ Hamiltonian (with nalgebra, the only place a library touches the physics) and forms $e^{-iHt}$ from its eigenvalues. The Trotter path never builds a big matrix: it moves amplitudes by hand, one local exponential at a time. A $ZZ$ term is a diagonal phase; an $X$ term is a one-qubit rotation:

1/// Apply e^{i*theta * Z_i Z_j} to the state: a diagonal phase, +theta on basis
2/// states whose bits i and j agree, -theta where they differ (the +/-1
3/// eigenvalues of Z_iZ_j). This is one factor of the ZZ interaction layer.
4fn apply_zz_phase(state: &mut [C], i: usize, j: usize, theta: f64) {
5    let bi = 1usize << i;
6    let bj = 1usize << j;
7    let ep = C::from_polar(1.0, theta);
8    let em = C::from_polar(1.0, -theta);
9    for (k, amp) in state.iter_mut().enumerate() {
10        let equal = (k & bi != 0) == (k & bj != 0);
11        *amp *= if equal { ep } else { em };
12    }
13}
14
15/// Apply e^{i*phi * X_q} = cos(phi) I + i sin(phi) X to qubit q: a 2x2 rotation
16/// mixing the amplitude pairs that differ in bit q. One factor of the X layer.
17fn apply_x_rot(state: &mut [C], q: usize, phi: f64) {
18    let bit = 1usize << q;
19    let co = C::new(phi.cos(), 0.0);
20    let si = C::new(0.0, phi.sin()); // i sin(phi)
21    for i in 0..state.len() {
22        if i & bit == 0 {
23            let j = i | bit;
24            let a = state[i];
25            let b = state[j];
26            state[i] = co * a + si * b;
27            state[j] = si * a + co * b;
28        }
29    }
30}

One Trotter step assembles those factors into a layer of $ZZ$ phases and a layer of $X$ rotations. The first-order step lays the two layers back to back; the second-order step symmetrizes, wrapping the $ZZ$ layer in two half-strength $X$ layers so the leading error cancels:

1/// One first-order (Lie–Trotter) step of size dt: the ZZ interaction layer
2/// followed by the X field layer. For H = -J ZZ - h X the exponents pick up a
3/// sign: e^{-i(-J ZZ)dt} = e^{+iJ dt ZZ}, e^{-i(-h X)dt} = e^{+ih dt X}.
4fn trotter1_step(state: &mut [C], bonds: &[(usize, usize)], n: usize, dt: f64) {
5    for &(i, j) in bonds {
6        apply_zz_phase(state, i, j, J * dt);
7    }
8    for q in 0..n {
9        apply_x_rot(state, q, H_FIELD * dt);
10    }
11}
12
13/// One second-order (symmetric Strang) step: half an X layer, a full ZZ layer,
14/// then the other half X layer. The symmetry cancels the leading error term.
15fn trotter2_step(state: &mut [C], bonds: &[(usize, usize)], n: usize, dt: f64) {
16    for q in 0..n {
17        apply_x_rot(state, q, 0.5 * H_FIELD * dt);
18    }
19    for &(i, j) in bonds {
20        apply_zz_phase(state, i, j, J * dt);
21    }
22    for q in 0..n {
23        apply_x_rot(state, q, 0.5 * H_FIELD * dt);
24    }
25}

The headline referee measures the error $\Vert {\psi }_{\text{Trotter}}-{\psi }_{\text{exact}}\Vert$ as a function of the step count $r$ and fits its slope on log-log axes. Both fitted exponents are refereed against the predicted $1/r$ and $1/r^2$ laws with a tolerance of $10^{-2}$ on the deviation — so “within a hundredth of $1$ and $2$” is not a description of one lucky run but the condition under which the lab refuses to finish. The achieved values are live in the panel:

1// --- The error-vs-steps study (first & second order) ------------------
2let r_list: Vec<usize> = vec![2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128, 256, 512, 1024];
3let exact_conv = exact_evolve(&evecs, &evals, &psi0_real, T_CONV);
4let mut error_curve: Vec<ErrPoint> = Vec::new();
5for &r in &r_list {
6    let p1 = trotter_evolve(&psi0, &bonds, n, T_CONV, r, 1);
7    let p2 = trotter_evolve(&psi0, &bonds, n, T_CONV, r, 2);
8    let err1 = l2_diff(&p1, &exact_conv);
9    let err2 = l2_diff(&p2, &exact_conv);
10    assert!(err1.is_finite() && err2.is_finite());
11    error_curve.push(ErrPoint { r, err1, err2 });
12}
13
14// --- Referees 3 & 4: the scaling exponents ----------------------------
15// Fit slope over a clean asymptotic window r in [8, 64] (small enough that
16// the second-order error stays well above the round-off floor).
17let fit_r_min = 8usize;
18let fit_r_max = 64usize;
19// …
20let slope1 = loglog_slope(&rs, &e1);
21let slope2 = loglog_slope(&rs, &e2);
22let dev1 = (slope1.abs() - 1.0).abs();
23let dev2 = (slope2.abs() - 2.0).abs();
24assert!(dev1.is_finite() && dev2.is_finite());
25referees.push(Referee {
26    name: "First-order error scaling  |slope| ≈ 1".into(),
27    value: dev1,
28    // Achieved 2.12e-4. Tightened from 0.1.
29    tol: 1e-3,
30    pass: dev1 < 1e-3,
31});

Seven referees run in all and write their verdicts — value and tolerance — to JSON, which the panel below reads live; no figure in this section is transcribed by hand. Two of the seven confirm that both evolutions preserve the norm (every local exponential is unitary by construction, so these are checks on the arithmetic, not on the physics). The five that can catch a wrong Hamiltonian all compare the hand-rolled Trotter path against the eigendecomposition: the two scaling exponents, the second-order error at the largest step count in the sweep, the magnetization curve, and the energy drift of the symmetric splitting. Every tolerance sits within a factor of five of the value achieved — read the two columns against each other in the panel.

That margin is deliberate and it is the honest part. A tolerance set orders of magnitude above what a lab actually achieves is not a referee, it is decoration: it will wave through a real regression and never notice. Tightening the second-order slope check to $10^{-2}$ has teeth — corrupting the $ZZ$ coupling by a few parts in a hundred thousand, far too little to see in any plot on this page, is enough to fail it, where the tolerance it replaced would have shrugged.

Run it yourself with cargo run --release in Rust-QML/ch09-trotter. This is the last lab of Part II, and it closes a loop opened in Chapter 1: the exponential state space that made a qubit such a poor safe is, for simulation, the whole point — the machine holds the $2^n$ amplitudes for free and pushes them forward with a shallow circuit.