Introduction

In 1982, Richard Feynman posed a deceptively simple question at a physics conference: if you want to simulate a quantum system — say, a molecule or a chain of magnetic spins — why not use a machine that is itself quantum?

The reason to ask is that classical computers are terrible at it. The state of a quantum system with n particles lives in a space of 2n2^{n} complex numbers. Double the particles and the memory requirement squares. A 50-qubit system already needs more classical memory than every computer on Earth combined. Nature runs these calculations instantly; we cannot afford to watch.

Hamiltonian simulation is the task of computing how a quantum system evolves in time given its energy function — its Hamiltonian, H. The evolution of the quantum state |ψ⟩ over time t is given by Schrödinger's equation, whose formal solution is the unitary operator e−iHte^{-iHt}. The challenge is to implement that operator efficiently on a quantum computer.

This was the application Feynman originally envisioned, and it remains the clearest example of quantum advantage: a quantum computer can track the state naturally, qubit by qubit, while a classical machine drowns in exponential overhead.

Watch Two Spins Evolve

The demo below simulates a 2-spin system under the Heisenberg XX interaction: two qubits coupled so that they exchange energy. The Hamiltonian is H = J(X⊗X + Y⊗Y), where X and Y are the Pauli matrices and J sets the coupling strength.

The Trotter product formula decomposes e−iHte^{-iHt} into short, alternating steps — one for each term in H — so the evolution can be built from simple quantum gates. Smaller steps mean higher accuracy but more gates; larger steps are faster but introduce error.

<p class="hint">Simulating a 2-spin system: H = J(X⊗X + Y⊗Y). The Trotter formula decomposes e<sup>−iHt</sup> into short alternating steps. Watch probability flow between the four basis states as the two spins exchange energy.</p>
<div class="controls">
  <label>Coupling J: <span id="jVal">1.0</span>
    <input id="jSlider" type="range" min="0.2" max="3.0" step="0.1" value="1.0">
  </label>
  <label>Trotter step dt: <span id="dtVal">0.10</span>
    <input id="dtSlider" type="range" min="0.02" max="0.5" step="0.02" value="0.10">
  </label>
</div>
<canvas id="chart" width="520" height="200"></canvas>
<div class="legend">
  <span class="dot" style="background:#1d6fa5"></span>|00⟩
  <span class="dot" style="background:#e06c1b"></span>|01⟩
  <span class="dot" style="background:#1a8a3f"></span>|10⟩
  <span class="dot" style="background:#9b3cc0"></span>|11⟩
</div>
<div class="btns">
  <button id="playBtn" type="button">Pause</button>
  <button id="resetBtn" type="button" class="ghost">Reset</button>
</div>
<div class="info" id="info">t = 0.00 | Trotter steps: 0</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .87rem; color: #444; margin: 0 0 .6rem; line-height: 1.45; }
.controls { display: flex; flex-direction: column; gap: .4rem; margin-bottom: .6rem; }
label { font-size: .9rem; display: flex; align-items: center; gap: .5rem; flex-wrap: wrap; }
input[type=range] { width: 160px; }
canvas { display: block; width: 100%; max-width: 520px; border: 1px solid #d0d7df;
         border-radius: 8px; background: #f8fafc; }
.legend { font-size: .82rem; display: flex; gap: .8rem; margin: .4rem 0; flex-wrap: wrap; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; margin-right: 3px; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-top: .4rem; }
button { font: 600 14px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
.info { font-size: .85rem; color: #555; margin-top: .35rem; font-family: ui-monospace, monospace; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// ── Hamiltonian Simulation Demo ────────────────────────────────────────────
// System: 2 qubits (4-dimensional Hilbert space).
// H = J*(XX + YY), where X,Y are Pauli matrices.
// State: complex vector of 4 amplitudes [c00, c01, c10, c11].
// Trotter: e^{-iHdt} ≈ e^{-iJ*XX*dt} * e^{-iJ*YY*dt}
//
// e^{-i*J*XX*t} on the full 4-dim space:
//   XX = [[0,0,0,1],[0,0,1,0],[0,1,0,0],[1,0,0,0]]
//   eigenvalues ±1, so e^{-i*theta*XX} = cos(theta)*I - i*sin(theta)*XX
//
// e^{-i*J*YY*t} on the full 4-dim space:
//   YY = [[0,0,0,-1],[0,0,1,0],[0,1,0,0],[-1,0,0,0]]
//   same eigenvalue structure (±1), so e^{-i*theta*YY} = cos(theta)*I - i*sin(theta)*YY

const COLORS = ['#1d6fa5','#e06c1b','#1a8a3f','#9b3cc0'];
const LABELS = ['|00⟩','|01⟩','|10⟩','|11⟩'];
const HIST   = 200; // history points to plot

// Complex number as [re, im]
const C = (r,i=0) => [r,i];
const cadd = ([ar,ai],[br,bi]) => [ar+br, ai+bi];
const cmul = ([ar,ai],[br,bi]) => [ar*br-ai*bi, ar*bi+ai*br];
const cscale = ([r,i],s) => [r*s, i*s];
const cneg  = ([r,i]) => [-r,-i];

// Apply e^{-i*theta*P} where P is XX or YY.
// Both XX and YY have the 4x4 structure:
//   Diagonal = I (multiplied by cos(theta))
//   Off-diag: XX swaps |00⟩<->|11⟩ and |01⟩<->|10⟩ with -i*sin factor
//             YY swaps |00⟩<->|11⟩ with +sin and |01⟩<->|10⟩ with -sin (times -i)
// Let's derive directly:
// XX matrix (indices 00=0,01=1,10=2,11=3):
//   XX|00⟩=|11⟩, XX|01⟩=|10⟩, XX|10⟩=|01⟩, XX|11⟩=|00⟩
// YY matrix: YY = (Y⊗Y)
//   Y = [[0,-i],[i,0]]
//   Y⊗Y |00⟩ = (-i)(−i)|11⟩ = -|11⟩
//   Y⊗Y |01⟩ = (-i)(i)|10⟩  = +|10⟩
//   Y⊗Y |10⟩ = (i)(-i)|01⟩  = +|01⟩
//   Y⊗Y |11⟩ = (i)(i)|00⟩   = -|00⟩
// So: e^{-i*theta*XX}: new[k] = cos(theta)*old[k] + -i*sin(theta)*(XX*old)[k]
//     e^{-i*theta*YY}: new[k] = cos(theta)*old[k] + -i*sin(theta)*(YY*old)[k]

function applyExpXX(state, theta) {
  const c = Math.cos(theta);
  const s = Math.sin(theta);
  const [c00,c01,c10,c11] = state;
  const isin = C(0, -s);
  // (-i*sin)*XX*state: XX swaps 0<->3 and 1<->2
  return [
    cadd(cscale(c00,c), cmul(isin,c11)),
    cadd(cscale(c01,c), cmul(isin,c10)),
    cadd(cscale(c10,c), cmul(isin,c01)),
    cadd(cscale(c11,c), cmul(isin,c00)),
  ];
}

function applyExpYY(state, theta) {
  const c = Math.cos(theta);
  const s = Math.sin(theta);
  const [c00,c01,c10,c11] = state;
  const isin = C(0, -s);
  // (-i*sin)*YY*state: YY|00⟩=-|11⟩, YY|01⟩=+|10⟩, YY|10⟩=+|01⟩, YY|11⟩=-|00⟩
  return [
    cadd(cscale(c00,c), cmul(isin, cneg(c11))),
    cadd(cscale(c01,c), cmul(isin, c10)),
    cadd(cscale(c10,c), cmul(isin, c01)),
    cadd(cscale(c11,c), cmul(isin, cneg(c00))),
  ];
}

function trotterStep(state, J, dt) {
  // e^{-i J XX dt} then e^{-i J YY dt}
  state = applyExpXX(state, J * dt);
  state = applyExpYY(state, J * dt);
  return state;
}

function prob(amp) { return amp[0]*amp[0] + amp[1]*amp[1]; }

// ── State ───────────────────────────────────────────────────────────────────────────
let J  = 1.0;
let dt = 0.10;
// Initial state |01⟩: c01 = 1, rest 0
let state = [C(0),C(1),C(0),C(0)];
let t = 0, steps = 0;
let running = true;
let history = [state.map(prob)]; // array of [p00,p01,p10,p11]
let animId = null;

const canvas = document.getElementById('chart');
const ctx    = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

function draw() {
  ctx.clearRect(0,0,W,H);
  // background
  ctx.fillStyle = '#f8fafc'; ctx.fillRect(0,0,W,H);
  // gridlines
  ctx.strokeStyle = '#e0e7ef'; ctx.lineWidth = 1;
  [0.25,0.5,0.75,1.0].forEach(v => {
    const y = H - v*H*0.88 - H*0.06;
    ctx.beginPath(); ctx.moveTo(0,y); ctx.lineTo(W,y); ctx.stroke();
    ctx.fillStyle='#aab'; ctx.font='10px system-ui';
    ctx.fillText((v*100).toFixed(0)+'%', 3, y-2);
  });

  const n = Math.min(history.length, HIST);
  if(n < 2) return;
  const slice = history.slice(-n);

  COLORS.forEach((col,ci) => {
    ctx.beginPath();
    ctx.strokeStyle = col; ctx.lineWidth = 2;
    slice.forEach((probs, xi) => {
      const x = (xi/(n-1))*W;
      const y = H - probs[ci]*H*0.88 - H*0.06;
      xi===0 ? ctx.moveTo(x,y) : ctx.lineTo(x,y);
    });
    ctx.stroke();
  });
}

function step() {
  state = trotterStep(state, J, dt);
  t += dt; steps++;
  history.push(state.map(prob));
  document.getElementById('info').textContent =
    fmt('t = %t% | Trotter steps: %steps%', { t: t.toFixed(2), steps: steps });
  draw();
}

function loop() {
  if(running) {
    step();
    animId = setTimeout(loop, 40);
  }
}

function resetSim() {
  clearTimeout(animId);
  state   = [C(0),C(1),C(0),C(0)];
  t = 0; steps = 0;
  history = [state.map(prob)];
  document.getElementById('info').textContent = 't =0.00| Trotter steps:0';
  draw();
  if(running) { animId = setTimeout(loop, 40); }
}

document.getElementById('jSlider').addEventListener('input', function() {
  J = parseFloat(this.value);
  document.getElementById('jVal').textContent = J.toFixed(1);
  resetSim();
});
document.getElementById('dtSlider').addEventListener('input', function() {
  dt = parseFloat(this.value);
  document.getElementById('dtVal').textContent = dt.toFixed(2);
  resetSim();
});
document.getElementById('playBtn').addEventListener('click', function() {
  running = !running;
  this.textContent = running ? 'Pause' : 'Play';
  if(running) { animId = setTimeout(loop, 40); }
});
document.getElementById('resetBtn').addEventListener('click', resetSim);

// Verify: at t = pi/(4J), energy should transfer fully from |01⟩ to |10⟩
// (complete swap for XX+YY at half period)
console.log(fmt('[hamiltonian-simulation] demo loaded, J=%J%, dt=%dt%', { J: J, dt: dt }));
draw();
animId = setTimeout(loop, 100);
Result
Console

Adjust the coupling J and Trotter step size dt and watch the four basis-state probabilities (|00⟩, |01⟩, |10⟩, |11⟩) oscillate over time. Start in the |01⟩ state and you will see energy flowing back and forth between the two spins — exactly the quantum phenomenon Feynman wanted to simulate.

The Real Complexity

How hard is Hamiltonian simulation?

  • Classically: exponentially hard in general. The state vector has 2n2^{n} amplitudes. Exact classical simulation of n qubits takes O(2n)O(2^{n}) memory and time. For even n = 50, this exceeds planetary-scale computing budgets.

  • The Trotter–Suzuki formula (1959/1991). If H = H1H_{1} + H2H_{2} + ··· + HkH_{k}, each term easy to exponentiate separately, then:

    e−iHt≈(e−iH1t/r e−iH2t/r ⋯ e−iHkt/r)re^{-iHt} \approx (e^{-iH_{1}t/r}\, e^{-iH_{2}t/r}\, \cdots\, e^{-iH_{k}t/r})^{r}

    for large r. The error shrinks as 1/r. This decomposes time evolution into a sequence of small, implementable quantum gates.

  • Seth Lloyd's 1996 result. The first proof that quantum computers can simulate local Hamiltonians efficiently — in polynomial time and gate count — using the Trotter formula. This established Hamiltonian simulation as the flagship quantum-over-classical speedup.

  • BQP-complete. The problem sits in BQP (quantum polynomial time) and is believed to be outside classical polynomial time (P and BPP), placing it at the heart of the quantum advantage debate.

  • Trotter error. Each Trotter step introduces a commutator error of order O((Δt)2)O((\Delta t)^2). Modern methods — product formulas of higher order, qubitization, and quantum signal processing — reduce this to near-optimal gate counts.

The key insight is that quantum hardware natively represents the state with n qubits, so there is no exponential overhead in storage. The simulation runs in the same Hilbert space as the real system.

Where It Matters

Hamiltonian simulation is not an abstract exercise — it is the reason quantum computers are being built:

  • Quantum chemistry. Simulating the electronic structure of molecules (their Hamiltonians are Fermionic) predicts reaction rates, bond energies, and catalytic properties. Classical methods fail for large or strongly correlated molecules; quantum simulation scales naturally.
  • Drug discovery. Proteins fold and bind under quantum mechanical forces. Simulating these interactions at atomic resolution could cut years off the drug-design pipeline.
  • Materials science. High-temperature superconductors, battery materials, and topological phases of matter are all governed by many-body Hamiltonians that are intractable classically.
  • Nuclear physics. Simulating the strong interaction between quarks inside protons requires lattice QCD — a huge Hamiltonian problem perfectly suited to quantum hardware.
  • Quantum error correction. Understanding noise in quantum devices themselves requires simulating open-system Hamiltonians, a stepping stone on the path to fault tolerance.

See also Shor's Algorithm and Grover Search for other BQP problems, and quantum error correction for the hardware context.

Conclusion

Hamiltonian simulation is perhaps the purest distillation of why quantum computing exists. Feynman's 1982 insight — that only a quantum machine can efficiently track a quantum system — translates directly into the Trotter product formula: break the time evolution into small steps, each implementable as a quantum gate, and let the hardware do what nature does naturally.

The problem is BQP-complete: any quantum computation can be recast as a Hamiltonian simulation, so solving it efficiently is equivalent to having a general quantum computer. For chemistry, materials, and biology, this is the speedup that matters most — not breaking encryption, but understanding atoms.

As fault-tolerant quantum processors scale up, Hamiltonian simulation will likely be the first task where they outpace every classical supercomputer in a scientifically meaningful way. Feynman asked the right question.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/hamiltonian-simulation/Content licensed under CC BY-NC 4.0.