Transverse-Field Ising Model

Simulate the real-time dynamics of a 1D transverse-field Ising model (TFIM) by Trotterizing its time-evolution operator — a landmark model of condensed-matter physics, a standard benchmark for quantum simulation and variational algorithms. The reference implementation runs today on the CPU statevector target (qpp-cpu); the GPU backends referenced below are on the roadmap and not yet released.

What it is

The Ising model describes a row of spins (qubits) where neighbouring spins couple to one another, while the whole chain is immersed in a transverse magnetic field. Its Hamiltonian is:

H = -J · Σ<i> Z_i Z_{i+1}  -  h · Σ<i> X_i

The two terms play opposing roles:

The crux is that X and Z do not commute: the field term constantly tries to flip spins that the coupling has aligned, so order and fluctuation compete. The ratio h/J decides which wins.

Intuition — the tug-of-war between order and fluctuation

Think of h/J as a dial. When h/J ≈ 0 the coupling dominates and the ground state is a neatly aligned ferromagnetic order. When h/J is large the field dominates and each spin is independently pulled along X, the bulk magnetization washed out by fluctuations into a disordered paramagnetic state. Between the two lies a critical point (h/J)c where the system undergoes a quantum phase transition.

Why it matters — quantum phase transitions

Unlike classical phase transitions driven by heating, a quantum phase transition happens at absolute zero, driven by quantum (not thermal) fluctuations, triggered purely by tuning the Hamiltonian parameter h/J. The 1D TFIM is special because it is exactly solvable (mapped to free fermions via the Jordan–Wigner transformation), with the critical point precisely at (h/J)c = 1. Near the critical point the correlation length diverges, the energy gap closes, and the system exhibits universal critical scaling.

Because the answer is known yet the problem is genuinely "hard", the TFIM is an ideal touchstone for testing quantum algorithms:

• Real-time dynamics. Evolve under exp(−i·H·t) to watch excitations (such as domain walls) propagate — exactly what this page does.
• VQE / ground state. Approximate the ground-state energy with a variational quantum eigensolver and check it point by point against the exact solution.
• Locating the transition. Sweep h/J and recover the critical point (h/J)c = 1 from observables like the magnetization.

Why it works — Trotterized time evolution

We want to apply the time-evolution operator exp(−i·H·t), but the X and ZZ terms in H do not commute, so we cannot simply "exponentiate them separately". The Trotter decomposition gives a controlled approximation: slice the evolution into many small steps dt, and within each step apply the field term and then the coupling term separately:

exp(−i·H·dt) ≈ exp(i·h·dt·ΣXi) · exp(i·J·dt·ΣZiZi+1)

Each individual term maps to a simple gate: the field term is a single-qubit rotation about X (rx); the ZZ coupling can be implemented with a pair of CNOTs sandwiching an rz. Repeating this small step n_steps times evolves the state forward to time n_steps·dt. The smaller dt, the more accurate the approximation — at the cost of more gates.

What it is good at — and its limits

• A trustworthy benchmark, because there is an exact solution to check against. Being able to compare point by point with the analytic result makes it the gold standard for validating simulators and hardware.
• The statevector grows exponentially with the number of spins. N spins need 2^N complex amplitudes — precisely the motivation for the planned GPU statevector backend (and multi-GPU sharding) on the roadmap; today the reference implementation runs on the qpp-cpu CPU target.
• Trotter error must be managed. A first-order decomposition has O(dt²) error per step; switching to a second-order (Suzuki–Trotter) step substantially lowers the discretization error.
• For low entanglement, there are shortcuts. For very long, weakly entangled chains, tensor-network methods are often more efficient than a dense statevector.

Building the Hamiltonian

Read it against the model above: spin.z(i) * spin.z(i+1) accumulates the ZZ couplings along the open chain term by term, spin.x(i) accumulates the transverse field, and the overall minus signs correspond to −J and −h, matching the formula one-to-one.

import qalgora
from qalgora import spin
import numpy as np

n_spins = 8
J, h = 1.0, 1.0

def tfim_hamiltonian(n, J, h):
    H = 0
    for i in range(n - 1):                 # ZZ couplings (open chain)
        H -= J * spin.z(i) * spin.z(i + 1)
    for i in range(n):                     # transverse field
        H -= h * spin.x(i)
    return H

hamiltonian = tfim_hamiltonian(n_spins, J, h)

Trotterized time evolution

Read it against the "Why it works" section above: rx(-2.0 * h * dt, q[i]) implements the single-step evolution of the field term; the x.ctrl then rz(-2.0 * J * dt, ...) then x.ctrl sandwich is the standard implementation of the ZZ coupling term. evolve starts from a single flipped spin (a domain wall) and repeatedly applies the step to evolve it forward.

@qalgora.kernel
def trotter_step(q: qalgora.qview, dt: float, J: float, h: float):
    n = q.size()
    # transverse-field term: exp(i h dt X_i)
    for i in range(n):
        rx(-2.0 * h * dt, q[i])
    # coupling term: exp(i J dt Z_i Z_{i+1})
    for i in range(n - 1):
        x.ctrl(q[i], q[i + 1])
        rz(-2.0 * J * dt, q[i + 1])
        x.ctrl(q[i], q[i + 1])

@qalgora.kernel
def evolve(n: int, steps: int, dt: float, J: float, h: float):
    q = qalgora.qvector(n)
    x(q[0])                       # start from a single flipped spin (domain wall)
    for _ in range(steps):
        trotter_step(q, dt, J, h)

Tracking the magnetization

Take the observable to be the average ⟨Z⟩ over the spins and use qalgora.observe to measure its expectation value at each time slice; this lets you plot the magnetization evolving in time — the spreading domain wall pulls the curve away from its initial value over time.

qalgora.set_target("qpp-cpu")      # runnable today — CPU statevector reference implementation
# qalgora.set_target("gpu")        # planned: GPU backend not yet released

dt, n_steps = 0.05, 100
magnetization = spin.z(0)
for i in range(1, n_spins):
    magnetization += spin.z(i)
magnetization /= n_spins

times, mz = [], []
for step in range(n_steps):
    res = qalgora.observe(evolve, magnetization, n_spins, step, dt, J, h)
    times.append(step * dt)
    mz.append(res.expectation())

print("final ⟨Z⟩ =", mz[-1])

Scaling up & next steps

• Increase n_spins and watch the CPU qpp-cpu wall-clock grow (the planned gpu backend targets the 30+ spin regime).
• Sweep h/J to locate the quantum phase transition in the steady-state magnetization.
• Use a second-order (Suzuki–Trotter) step to cut the time-discretization error.
• For very long chains with low entanglement, a tensornet target is planned.