Hamiltonian Simulation

Quantum Simulation

Chemistry & Materials Science

Quantum simulators leverage quantum mechanics itself to model molecular and material systems. Algorithms like Trotterization and Quantum Phase Estimation enable simulations that are intractable for even the most powerful classical supercomputers.

Why Quantum Simulation?

Simulating quantum systems classically is extraordinarily difficult. The state space of a quantum system grows exponentially with the number of particles, meaning a system with just a few dozen electrons might require more classical memory than exists on Earth.

Richard Feynman proposed in 1982 that quantum computers could naturally simulate quantum physics. Instead of tracking an exponential number of classical variables, a quantum computer uses qubits that directly map to the quantum system's degrees of freedom.

•Classical Bottleneck: Storing the wavefunction of N spin-1/2 particles requires 2^N complex numbers.
•Natural Mapping: N qubits can directly represent N spin-1/2 particles using the same Hilbert space.
•Applications: Drug discovery, catalyst design, battery materials, high-temperature superconductors.

Trotterization

Trotter-Suzuki decomposition is the workhorse of digital quantum simulation. It approximates the time evolution operator e^(-iHt) by breaking the Hamiltonian H into small, locally-interacting pieces that can be implemented as quantum gates.

While approximate, the error can be made arbitrarily small by increasing the number of Trotter steps. Higher-order Suzuki formulas reduce the number of steps needed, making simulations more efficient on near-term devices.

U(t) = e^{-iHt}

Time Evolution Operator

e^{-iHt} \approx \bigl(e^{-iH_1 t/n} e^{-iH_2 t/n} \cdots e^{-iH_k t/n}\bigr)^n

First-Order Trotter-Suzuki

\varepsilon = \mathcal{O}\!\left(\frac{t^2}{n}\right)

Trotter Error

•Hamiltonian Splitting: Decompose H into a sum of local terms H = H_1 + H_2 + ... + H_k.
•Trotter Step: Approximate e^(-iHt) ≈ (e^(-iH_1t/n) · e^(-iH_2t/n) · ...)^n for large n.
•Error Scaling: First-order Trotter error scales as O(t^2/n). Higher-order formulas improve this.

Quantum Phase Estimation (QPE)

QPE extracts eigenvalues of a Hamiltonian by encoding phase information into a register of ancilla qubits. Combined with a state preparation circuit, it can find ground-state energies to chemical precision — a task crucial for predicting reaction rates and molecular stability.

Variational Quantum Eigensolver (VQE) offers a NISQ-friendly alternative, using a classical optimizer to tune a quantum circuit that prepares low-energy states without the deep circuits required by QPE.

•Ground State Energy: QPE can estimate eigenvalues to arbitrary precision given sufficient ancilla qubits.
•VQE for NISQ: Shallow circuits + classical optimization make VQE suitable for current hardware.
•Dynamics: Trotterization simulates time evolution, revealing how systems change over time.

Code in Action

Runnable implementations you can copy and experiment with.

Trotterized Time Evolution

A simple 2-qubit Trotter step simulating the Heisenberg XXX model Hamiltonian.

trotter_xxx.py

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
import numpy as np

# Parameters
J = 1.0      # coupling strength
t = 0.5      # total time
n_steps = 4  # Trotter steps
dt = t / n_steps

qc = QuantumCircuit(2, 2)

# Initial state: |01>
qc.x(1)

# Trotter steps
for _ in range(n_steps):
    # e^(-i H_XX dt) via XX rotation
    qc.h([0, 1])
    qc.cx(0, 1)
    qc.rz(-2 * J * dt, 1)
    qc.cx(0, 1)
    qc.h([0, 1])
    
    # e^(-i H_YY dt) via YY rotation
    qc.rx(np.pi / 2, [0, 1])
    qc.cx(0, 1)
    qc.rz(-2 * J * dt, 1)
    qc.cx(0, 1)
    qc.rx(-np.pi / 2, [0, 1])
    
    # e^(-i H_ZZ dt) via ZZ rotation
    qc.cx(0, 1)
    qc.rz(-2 * J * dt, 1)
    qc.cx(0, 1)

qc.measure([0, 1], [0, 1])

simulator = AerSimulator()
job = simulator.run(qc, shots=4096)
print(job.result().get_counts())

Trotterized Evolution with PennyLane

A PennyLane implementation of Trotterized time evolution for the 2-qubit Heisenberg XXX model.

trotter_pennylane.py

import pennylane as qml
import numpy as np

n_qubits = 2
dev = qml.device('default.qubit', wires=n_qubits)

J = 1.0
t = 0.5
n_steps = 4
dt = t / n_steps

@qml.qnode(dev)
def trotter_circuit():
    qml.PauliX(wires=1)
    for _ in range(n_steps):
        # XX term
        qml.Hadamard(wires=0)
        qml.Hadamard(wires=1)
        qml.CNOT(wires=[0, 1])
        qml.RZ(-2 * J * dt, wires=1)
        qml.CNOT(wires=[0, 1])
        qml.Hadamard(wires=0)
        qml.Hadamard(wires=1)
        # YY term
        qml.RX(np.pi / 2, wires=0)
        qml.RX(np.pi / 2, wires=1)
        qml.CNOT(wires=[0, 1])
        qml.RZ(-2 * J * dt, wires=1)
        qml.CNOT(wires=[0, 1])
        qml.RX(-np.pi / 2, wires=0)
        qml.RX(-np.pi / 2, wires=1)
        # ZZ term
        qml.CNOT(wires=[0, 1])
        qml.RZ(-2 * J * dt, wires=1)
        qml.CNOT(wires=[0, 1])
    return qml.counts(wires=range(n_qubits))

counts = trotter_circuit()
print("Trotter counts:", counts)

Explore next

Variational Algorithms

Explore next

Quantum Machine Learning

from qiskit import QuantumCircuit from qiskit_aer import AerSimulator import numpy as np # Parameters J = 1.0 # coupling strength t = 0.5 # total time n_steps = 4 # Trotter steps dt = t / n_steps qc = QuantumCircuit(2, 2) # Initial state: |01> qc.x(1) # Trotter steps for _ in range(n_steps): # e^(-i H_XX dt) via XX rotation qc.h([0, 1]) qc.cx(0, 1) qc.rz(-2 * J * dt, 1) qc.cx(0, 1) qc.h([0, 1]) # e^(-i H_YY dt) via YY rotation qc.rx(np.pi / 2, [0, 1]) qc.cx(0, 1) qc.rz(-2 * J * dt, 1) qc.cx(0, 1) qc.rx(-np.pi / 2, [0, 1]) # e^(-i H_ZZ dt) via ZZ rotation qc.cx(0, 1) qc.rz(-2 * J * dt, 1) qc.cx(0, 1) qc.measure([0, 1], [0, 1]) simulator = AerSimulator() job = simulator.run(qc, shots=4096) print(job.result().get_counts())

import pennylane as qml import numpy as np n_qubits = 2 dev = qml.device('default.qubit', wires=n_qubits) J = 1.0 t = 0.5 n_steps = 4 dt = t / n_steps @qml.qnode(dev) def trotter_circuit(): qml.PauliX(wires=1) for _ in range(n_steps): # XX term qml.Hadamard(wires=0) qml.Hadamard(wires=1) qml.CNOT(wires=[0, 1]) qml.RZ(-2 * J * dt, wires=1) qml.CNOT(wires=[0, 1]) qml.Hadamard(wires=0) qml.Hadamard(wires=1) # YY term qml.RX(np.pi / 2, wires=0) qml.RX(np.pi / 2, wires=1) qml.CNOT(wires=[0, 1]) qml.RZ(-2 * J * dt, wires=1) qml.CNOT(wires=[0, 1]) qml.RX(-np.pi / 2, wires=0) qml.RX(-np.pi / 2, wires=1) # ZZ term qml.CNOT(wires=[0, 1]) qml.RZ(-2 * J * dt, wires=1) qml.CNOT(wires=[0, 1]) return qml.counts(wires=range(n_qubits)) counts = trotter_circuit() print("Trotter counts:", counts)