Hamiltonian Simulation with Trotterization in Qiskit

Simulating quantum systems (how a molecule evolves in time, how electrons interact in a lattice) is one of the most compelling uses of quantum computers. Richard Feynman’s original motivation for quantum computing was exactly this: classical computers require exponential resources to track all the amplitude combinations of an evolving quantum state. A quantum computer can represent those states directly.

The standard tool for converting Hamiltonian simulation into a quantum circuit is Trotterization. This tutorial explains how it works and how to implement it in Qiskit.

Setup

pip install qiskit qiskit-aer

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.quantum_info import SparsePauliOp, Statevector
from qiskit.synthesis import LieTrotter, SuzukiTrotter
from qiskit_aer import AerSimulator

The problem: simulating time evolution

Given a Hamiltonian H (the energy operator of a quantum system), the quantum state at time t evolves as:

|ψ(t)⟩ = e^(-iHt) |ψ(0)⟩

The operator e^(-iHt) is unitary (time evolution preserves probability). The challenge: for most interesting Hamiltonians, this operator cannot be implemented directly as a quantum circuit. It would require exponentially many terms.

Trotterization breaks the problem into manageable pieces. If H = A + B where A and B do not commute, then:

e^(-i(A+B)t) ≈ (e^(-iAt/n) e^(-iBt/n))^n

This approximation becomes exact as n → ∞. For finite n, the error is O(t²/n). The first-order Trotter formula, and its second-order Suzuki improvement, are the standard tools.

Example 1: A two-site Ising model

The transverse-field Ising model is a standard benchmark:

H = J * (Z₀Z₁) + h * (X₀ + X₁)

where J is the coupling strength and h is the transverse field.

# Define the Hamiltonian
J = 1.0  # coupling
h = 0.5  # transverse field

H = SparsePauliOp.from_list([
    ('ZZ', J),     # Z0 Z1 coupling
    ('XI', h),     # X on qubit 0, I on qubit 1
    ('IX', h),     # I on qubit 0, X on qubit 1
])
print("Hamiltonian:")
print(H)

# Time parameters
t = 1.0       # total time
n_steps = 4   # Trotter steps

# Build the Trotter circuit
evo_gate = PauliEvolutionGate(H, time=t, synthesis=LieTrotter(reps=n_steps))
qc = QuantumCircuit(2)
qc.h([0, 1])        # Start in |+>|+> state
qc.append(evo_gate, [0, 1])
qc = qc.decompose()  # Expand the evolution gate into native gates

print(f"\nTrotter circuit depth: {qc.depth()}")
print(f"Gate count: {qc.count_ops()}")

Simulating the time evolution

def evolve_state(H: SparsePauliOp, t: float, n_steps: int, initial_state: list) -> Statevector:
    """
    Evolve an initial state under Hamiltonian H for time t using n_steps Trotter steps.
    Returns the final statevector.
    """
    qc = QuantumCircuit(H.num_qubits)
    qc.initialize(initial_state)
    evo = PauliEvolutionGate(H, time=t, synthesis=LieTrotter(reps=n_steps))
    qc.append(evo, range(H.num_qubits))
    qc = qc.decompose()
    
    sv = Statevector(qc)
    return sv

# Initial state: |00>
n_qubits = 2
initial = [1, 0, 0, 0]  # |00> in computational basis

# Evolve for different times
times = np.linspace(0, 2*np.pi, 20)
z0_expectations = []

Z0 = SparsePauliOp.from_list([('IZ', 1.0)])  # Z on qubit 0

for t in times:
    sv = evolve_state(H, t, n_steps=8, initial_state=initial)
    exp_z0 = sv.expectation_value(Z0).real
    z0_expectations.append(exp_z0)

print("Time evolution of <Z0>:")
for t, val in zip(times[::4], z0_expectations[::4]):
    print(f"  t={t:.2f}: <Z0> = {val:.4f}")

First-order vs second-order Trotter

The first-order Trotter formula (LieTrotter) has error O(t²/n). The second-order Suzuki formula (SuzukiTrotter order=2) has error O(t³/n²) but requires slightly more gates per step. For deep circuits, the second-order formula often gives much better accuracy per total gate count.

# Compare accuracy for different Trotter formulas
from qiskit.synthesis import SuzukiTrotter
from scipy.linalg import expm

# Use the transverse-field Ising model: [ZZ, XI] != 0, so Trotter error is non-trivial.
H_small = SparsePauliOp.from_list([('ZZ', 1.0), ('XI', 0.5), ('IX', 0.5)])
t_test = 1.0
initial_test = [1, 0, 0, 0]

H_matrix = H_small.to_matrix()
U_exact = expm(-1j * H_matrix * t_test)
psi_exact = U_exact @ np.array(initial_test, dtype=complex)

obs = SparsePauliOp.from_list([('ZI', 1.0)])  # Z on first qubit
exact_val = np.real(psi_exact.conj() @ obs.to_matrix() @ psi_exact)

print("Trotter approximation error (|<Z0> - <Z0>_exact|):")
print(f"{'Steps':<8} {'LieTrotter':<15} {'SuzukiTrotter (2nd)':<20}")

for n in [1, 2, 4, 8, 16]:
    # First order
    sv1 = evolve_state(H_small, t_test, n_steps=n, initial_state=initial_test)
    err1 = abs(sv1.expectation_value(obs).real - exact_val)

    # Second order
    qc2 = QuantumCircuit(2)
    qc2.initialize(initial_test)
    evo2 = PauliEvolutionGate(H_small, time=t_test, synthesis=SuzukiTrotter(order=2, reps=n))
    qc2.append(evo2, [0, 1])
    qc2 = qc2.decompose()
    sv2 = Statevector(qc2)
    err2 = abs(sv2.expectation_value(obs).real - exact_val)

    print(f"{n:<8} {err1:<15.6f} {err2:<20.6f}")

The second-order formula converges much faster: errors shrink as O(1/n²) vs O(1/n) for the first-order formula. Use SuzukiTrotter(order=2) by default for production simulations.

Example 2: Heisenberg chain (XXX model)

The Heisenberg model describes interacting spins and is harder to simulate classically for large systems:

def heisenberg_hamiltonian(n_sites: int, J: float = 1.0) -> SparsePauliOp:
    """Build a 1D Heisenberg XXX chain Hamiltonian."""
    terms = []
    for i in range(n_sites - 1):
        # XX coupling between sites i and i+1
        xx_str = 'I' * i + 'XX' + 'I' * (n_sites - i - 2)
        terms.append((xx_str[::-1], J))  # Qiskit: rightmost qubit = index 0
        
        # YY coupling
        yy_str = 'I' * i + 'YY' + 'I' * (n_sites - i - 2)
        terms.append((yy_str[::-1], J))
        
        # ZZ coupling
        zz_str = 'I' * i + 'ZZ' + 'I' * (n_sites - i - 2)
        terms.append((zz_str[::-1], J))
    
    return SparsePauliOp.from_list(terms)

# 4-site Heisenberg chain
H_heis = heisenberg_hamiltonian(4, J=1.0)
print(f"Heisenberg Hamiltonian ({H_heis.num_qubits} sites):")
print(f"  Terms: {len(H_heis)}")
print(f"  Paulis: {H_heis.paulis}")

# Evolve for time t=0.5
qc_heis = QuantumCircuit(4)
qc_heis.x([0, 2])  # Start in |1010> -- Neel state
evo_heis = PauliEvolutionGate(H_heis, time=0.5, synthesis=SuzukiTrotter(order=2, reps=4))
qc_heis.append(evo_heis, range(4))
qc_heis = qc_heis.decompose()

print(f"\nHeisenberg circuit (4 sites, t=0.5):")
print(f"  Depth: {qc_heis.depth()}")
print(f"  Two-qubit gates: {qc_heis.count_ops().get('cx', 0)}")

Running on a noisy simulator

Trotter circuits are deep, and noise accumulates. This is why Hamiltonian simulation is one of the most demanding applications for NISQ hardware.

from qiskit_aer.noise import NoiseModel, depolarizing_error
from qiskit import transpile

# Build a realistic noise model (0.5% two-qubit gate error)
noise_model = NoiseModel()
noise_model.add_all_qubit_quantum_error(depolarizing_error(0.005, 2), ['cx'])
noise_model.add_all_qubit_quantum_error(depolarizing_error(0.001, 1), ['h', 'rx', 'ry', 'rz'])

backend = AerSimulator(noise_model=noise_model)

# Simple 2-site Ising -- compare ideal vs noisy
H_ising = SparsePauliOp.from_list([('ZZ', 1.0), ('XI', 0.5), ('IX', 0.5)])
t_sim = 1.5

# Ideal statevector
qc_ideal = QuantumCircuit(2)
qc_ideal.initialize([1, 0, 0, 0])
evo = PauliEvolutionGate(H_ising, time=t_sim, synthesis=SuzukiTrotter(order=2, reps=6))
qc_ideal.append(evo, [0, 1])
qc_ideal = qc_ideal.decompose()
sv_ideal = Statevector(qc_ideal)
exp_ideal = sv_ideal.expectation_value(SparsePauliOp.from_list([('ZI', 1.0)])).real

# Noisy shot-based simulation
qc_noisy = qc_ideal.copy()
qc_noisy.save_density_matrix()
qc_noisy_t = transpile(qc_noisy, backend)
result = backend.run(qc_noisy_t, shots=4096).result()
dm = result.data()['density_matrix']
obs_matrix = SparsePauliOp.from_list([('ZI', 1.0)]).to_matrix()
exp_noisy = np.real(np.trace(np.array(dm) @ obs_matrix))

print(f"Ideal <Z0>:  {exp_ideal:.4f}")
print(f"Noisy <Z0>:  {exp_noisy:.4f}")
print(f"Noise error: {abs(exp_ideal - exp_noisy):.4f}")

Practical guidance

Choosing n_steps (Trotter steps):

• Start with n = 1 and increase until the observable converges
• For first-order: error ~ O(t²/n). Double n, halve the error
• For second-order: error ~ O(t³/n²). Double n, quarter the error
• Check convergence by comparing two consecutive values

Choosing the Trotter formula:

• LieTrotter(reps=n): first order, fewer gates but larger error
• SuzukiTrotter(order=2, reps=n): second order, recommended default
• SuzukiTrotter(order=4, reps=n): fourth order, best accuracy per step, most gates

Circuit depth vs accuracy: Every Trotter step adds depth. On NISQ hardware, more depth means more noise. There is a fundamental tension: more steps = better approximation but more noise. Use error mitigation (ZNE) to partially compensate, or use fewer steps and accept the Trotter error.

What to simulate: Hamiltonian simulation is most useful for:

• Quantum chemistry: molecular ground state dynamics, reaction rates
• Condensed matter: spin chain dynamics, topological phases
• High-energy physics: lattice gauge theories

For static ground state energy (not dynamics), use VQE instead. Trotter is for time evolution problems.