Hamiltonian Simulation with Qiskit: Pauli Evolution and Trotterization

Hamiltonian simulation is one of the most compelling near-term applications of quantum computers. Richard Feynman’s original motivation for quantum computing was precisely this: a quantum system can simulate another quantum system exponentially more efficiently than any classical computer can. Simulating the time evolution of a quantum many-body system is classically intractable beyond a few dozen spins, yet it is a natural task for a quantum processor.

The fundamental challenge is that the time evolution operator U(t) = exp(-iHt) cannot typically be implemented as a single quantum gate. For a many-body Hamiltonian, H is a sum of non-commuting terms. The Trotter-Suzuki decomposition converts this sum into a sequence of individually implementable exponentials, at the cost of a controllable approximation error.

The Transverse-Field Ising Model

The transverse-field Ising model (TFIM) is the workhorse of quantum simulation research. It is simple enough to write down but rich enough to exhibit quantum phase transitions and non-trivial dynamics. For N spins on a 1D chain with periodic boundary conditions:

H = -J * sum_{i} Z_i Z_{i+1}  -  h * sum_{i} X_i

The first term couples neighboring spins via ZZ interactions. The second term is the transverse magnetic field, with strength h, that drives quantum fluctuations. When h/J < 1 the system is in a ferromagnetic phase; when h/J > 1 it is in a paramagnetic phase. The phase transition at h/J = 1 is a canonical example of a quantum phase transition.

Time evolution under this Hamiltonian mixes the ZZ coupling and X field terms. Because Z_i Z_{i+1} and X_i do not commute (they share qubit i), we cannot simply exponentiate each term separately and multiply. We need Trotterization.

First-Order Trotterization

The first-order Lie-Trotter formula splits the evolution into small time steps of size dt:

exp(-iHt) ≈ [exp(-i H_ZZ dt) * exp(-i H_X dt)]^(t/dt)

where H_ZZ = -J * sum Z_i Z_{i+1} and H_X = -h * sum X_i. Each factor is easy to implement:

• exp(-i H_ZZ dt) is a product of ZZ rotation gates (Rzz) on neighboring pairs.
• exp(-i H_X dt) is a product of Rx rotations on each qubit.

The approximation error per step scales as O(dt^2), so total error for time t scales as O(dt * t) (linear in dt, first order).

Second-Order Suzuki-Trotter

The Suzuki-Trotter product formula of second order is significantly more accurate at the same cost:

exp(-iHt) ≈ [exp(-i H_ZZ dt/2) * exp(-i H_X dt) * exp(-i H_ZZ dt/2)]^(t/dt)

This symmetric decomposition cancels the leading error term, leaving O(dt^3) error per step and O(dt^2 * t) total error. For the same total simulation time, second-order Trotter needs far fewer steps to reach the same accuracy.

Implementation with Qiskit

Qiskit’s quantum_info and synthesis modules provide SparsePauliOp for defining Hamiltonians and PauliEvolutionGate with SuzukiTrotter for building the evolution circuit:

import numpy as np
from qiskit.quantum_info import SparsePauliOp, Statevector
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import SuzukiTrotter, LieTrotter
from qiskit import QuantumCircuit
import matplotlib.pyplot as plt

# Parameters
N = 4          # number of spins
J = 1.0        # ZZ coupling strength
h = 0.5        # transverse field strength
total_time = 3.0

def build_tfim_hamiltonian(N, J, h):
    """Build the transverse-field Ising Hamiltonian as a SparsePauliOp."""
    pauli_list = []

    # ZZ coupling terms: -J * Z_i Z_{i+1}
    for i in range(N - 1):
        zz_string = ['I'] * N
        zz_string[i] = 'Z'
        zz_string[i + 1] = 'Z'
        pauli_list.append((''.join(reversed(zz_string)), -J))

    # Transverse field terms: -h * X_i
    for i in range(N):
        x_string = ['I'] * N
        x_string[i] = 'X'
        pauli_list.append((''.join(reversed(x_string)), -h))

    return SparsePauliOp.from_list(pauli_list)

H = build_tfim_hamiltonian(N, J, h)
print("Hamiltonian terms:")
for term, coeff in zip(H.paulis, H.coeffs):
    print(f"  {coeff.real:+.2f} * {term}")

Now build and compare the Trotter circuits:

def build_trotter_circuit(H, total_time, n_steps, order=1):
    """
    Build a Trotter evolution circuit.
    order=1: LieTrotter, order=2: SuzukiTrotter second order.
    """
    N = H.num_qubits
    dt = total_time / n_steps

    if order == 1:
        synthesizer = LieTrotter(reps=1)
    else:
        synthesizer = SuzukiTrotter(order=2, reps=1)

    qc = QuantumCircuit(N)

    # Initialize in all-up state |0...0>
    # (no gates needed; |0...0> is the default)

    # Apply n_steps Trotter steps
    evo_gate = PauliEvolutionGate(H, time=dt, synthesis=synthesizer)
    for _ in range(n_steps):
        qc.append(evo_gate, range(N))

    qc = qc.decompose(reps=3)
    return qc

# Build circuits with different step counts
qc_1st_4  = build_trotter_circuit(H, total_time, n_steps=4,  order=1)
qc_1st_16 = build_trotter_circuit(H, total_time, n_steps=16, order=1)
qc_2nd_4  = build_trotter_circuit(H, total_time, n_steps=4,  order=2)
qc_2nd_8  = build_trotter_circuit(H, total_time, n_steps=8,  order=2)

print(f"1st order, 4 steps  - depth: {qc_1st_4.depth()}")
print(f"1st order, 16 steps - depth: {qc_1st_16.depth()}")
print(f"2nd order, 4 steps  - depth: {qc_2nd_4.depth()}")
print(f"2nd order, 8 steps  - depth: {qc_2nd_8.depth()}")

Computing Energy Expectation Values

We track the energy over time to see how Trotter error accumulates:

def energy_vs_time(H, total_time, n_steps, order=1):
    """
    Simulate Trotter evolution and compute <H> at each step.
    Uses Statevector simulation for exactness.
    """
    N = H.num_qubits
    dt = total_time / n_steps
    H_matrix = H.to_matrix()

    if order == 1:
        synthesizer = LieTrotter(reps=1)
    else:
        synthesizer = SuzukiTrotter(order=2, reps=1)

    energies = []
    times = []

    # Start in |0...0>
    sv = Statevector.from_label('0' * N)

    for step in range(n_steps + 1):
        # Energy expectation value
        energy = sv.expectation_value(H_matrix).real
        energies.append(energy)
        times.append(step * dt)

        if step < n_steps:
            # Evolve by one Trotter step
            evo_gate = PauliEvolutionGate(H, time=dt, synthesis=synthesizer)
            step_qc = QuantumCircuit(N)
            step_qc.append(evo_gate, range(N))
            step_qc = step_qc.decompose(reps=3)
            sv = sv.evolve(step_qc)

    return np.array(times), np.array(energies)

# Exact energy (many steps = high accuracy reference)
t_exact, E_exact = energy_vs_time(H, total_time, n_steps=200, order=2)

# Coarse Trotter approximations
t_1st, E_1st = energy_vs_time(H, total_time, n_steps=8, order=1)
t_2nd, E_2nd = energy_vs_time(H, total_time, n_steps=8, order=2)

# Plot
plt.figure(figsize=(10, 5))
plt.plot(t_exact, E_exact, 'k-', linewidth=2, label='Exact (200 steps, 2nd order)')
plt.plot(t_1st, E_1st, 'r--o', markersize=5, label='1st order Trotter (8 steps)')
plt.plot(t_2nd, E_2nd, 'b--s', markersize=5, label='2nd order Trotter (8 steps)')
plt.xlabel('Time')
plt.ylabel('Energy <H>')
plt.title('TFIM Energy vs Time: Trotter Error Comparison')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('trotter_energy.png', dpi=150)
plt.show()
print("Plot saved to trotter_energy.png")

Trotter Error vs. Step Size

We can directly measure how the accumulated Trotter error depends on the number of steps:

def trotter_error(H, total_time, n_steps, order):
    """Compute norm of difference between Trotter and exact evolution."""
    N = H.num_qubits

    # Exact statevector (via matrix exponentiation)
    H_matrix = H.to_matrix()
    U_exact = np.linalg.matrix_power(
        np.eye(2**N), 1  # placeholder
    )
    from scipy.linalg import expm
    U_exact = expm(-1j * H_matrix * total_time)
    psi0 = np.zeros(2**N)
    psi0[0] = 1.0
    psi_exact = U_exact @ psi0

    # Trotter statevector
    if order == 1:
        synthesizer = LieTrotter(reps=1)
    else:
        synthesizer = SuzukiTrotter(order=2, reps=1)

    sv = Statevector.from_label('0' * N)
    dt = total_time / n_steps
    for _ in range(n_steps):
        evo_gate = PauliEvolutionGate(H, time=dt, synthesis=synthesizer)
        step_qc = QuantumCircuit(N)
        step_qc.append(evo_gate, range(N))
        step_qc = step_qc.decompose(reps=3)
        sv = sv.evolve(step_qc)

    psi_trotter = sv.data
    return np.linalg.norm(psi_exact - psi_trotter)

step_counts = [2, 4, 8, 16, 32, 64]
errors_1st = [trotter_error(H, total_time=2.0, n_steps=n, order=1) for n in step_counts]
errors_2nd = [trotter_error(H, total_time=2.0, n_steps=n, order=2) for n in step_counts]
dt_values  = [2.0 / n for n in step_counts]

plt.figure(figsize=(8, 5))
plt.loglog(dt_values, errors_1st, 'r-o', label='1st order (LieTrotter)')
plt.loglog(dt_values, errors_2nd, 'b-s', label='2nd order (SuzukiTrotter)')

# Reference lines: slope 1 and slope 2
dt_ref = np.array(dt_values)
plt.loglog(dt_ref, dt_ref * errors_1st[0] / dt_values[0], 'r--', alpha=0.4, label='O(dt^1)')
plt.loglog(dt_ref, dt_ref**2 * errors_2nd[0] / dt_values[0]**2, 'b--', alpha=0.4, label='O(dt^2)')

plt.xlabel('Step size dt')
plt.ylabel('State vector error (L2 norm)')
plt.title('Trotter Error Scaling: 1st vs 2nd Order')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('trotter_error_scaling.png', dpi=150)
plt.show()

The log-log plot will show two straight lines with slopes 1 and 2, confirming the theoretical scaling. For the same accuracy target, second-order Trotter needs roughly sqrt(N_steps_1st) steps, a significant circuit depth saving on hardware.

Why Quantum Simulation Matters for Condensed Matter Physics

The transverse-field Ising model is a toy example, but the same techniques apply to:

• Hubbard model: models strongly correlated electron systems, relevant to high-temperature superconductivity. Classical simulation requires exponential resources in system size.
• Quantum chemistry: molecular Hamiltonians expressed as sums of Pauli operators. VQE and Hamiltonian simulation together tackle the electronic structure problem.
• Lattice gauge theories: QED and QCD on a lattice, relevant to particle physics. The qubit encoding maps gauge fields to quantum registers.

The utility threshold for Hamiltonian simulation, where a quantum computer provides a meaningful advantage over the best classical algorithm, is estimated to require 100-1000 logical qubits with moderate circuit depth. For the TFIM and related models on chains of 50-100 sites, classical methods based on tensor networks (DMRG) remain competitive. Quantum advantage kicks in for 2D systems with strong frustration, long evolution times, or dynamical quantities that tensor networks cannot efficiently compute.

Qiskit’s PauliEvolutionGate and SuzukiTrotter provide a clean interface for prototyping these simulations and testing them at small scale before targeting real hardware.