Simulating the Transverse-Field Ising Model with Qiskit

The transverse-field Ising model (TFIM) is one of the simplest yet richest models in condensed matter physics. It exhibits a quantum phase transition driven entirely by quantum fluctuations, making it a canonical testbed for quantum simulation. On near-term quantum hardware, it is also tractable enough to implement with shallow Trotter circuits, which makes it an ideal first foray into Hamiltonian simulation with Qiskit.

The Model

The 1D TFIM on $n$ sites is governed by the Hamiltonian:

$H = -J \sum_{i=0}^{n-2} Z_i Z_{i+1} - h \sum_{i=0}^{n-1} X_i$

The first term, with coupling $J$, energetically prefers neighboring spins aligned along the $z$-axis. The second term, with transverse field $h$, introduces quantum fluctuations that tend to align spins along the $x$-axis. At $h/J \ll 1$, the ground state is ferromagnetic (ordered). At $h/J \gg 1$, it is paramagnetic (disordered). The transition occurs exactly at $h/J = 1$ in the thermodynamic limit.

Trotterization

To simulate time evolution $e^{-iHt}$ on a quantum circuit, we use the first-order Trotter decomposition. We split $H = H_{ZZ} + H_X$ and approximate:

$e^{-iHt} \approx \left(e^{-iH_{ZZ}\Delta t} \cdot e^{-iH_X \Delta t}\right)^r$

where $\Delta t = t / r$ and $r$ is the number of Trotter steps. The error scales as $O(r \Delta t^2)$, so more steps give higher fidelity at the cost of deeper circuits.

Each $e^{-iH_{ZZ}\Delta t}$ layer applies $\text{ZZPhase}(\theta) = e^{-i\theta Z_i Z_{i+1}}$ gates with $\theta = J \Delta t$, implemented as an $R_{ZZ}$ gate. Each $e^{-iH_X\Delta t}$ layer applies $R_X(2h\Delta t)$ rotations on every qubit.

Setup

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.primitives import StatevectorEstimator
from qiskit.quantum_info import SparsePauliOp, Statevector
import matplotlib.pyplot as plt

n_qubits = 4
J = 1.0

Building the Trotter Circuit

n_qubits = 4

def build_trotter_circuit(n, J, h, dt, n_steps):
    """Build a first-order Trotter circuit for the 1D TFIM."""
    qc = QuantumCircuit(n)

    # Start in the all-|+> state (eigenstate of X)
    qc.h(range(n))

    for _ in range(n_steps):
        # ZZ interaction layer
        for i in range(n - 1):
            qc.rzz(2 * J * dt, i, i + 1)

        # Transverse field layer
        for i in range(n):
            qc.rx(2 * h * dt, i)

    return qc

# Inspect the circuit for a few steps
qc_example = build_trotter_circuit(n_qubits, J=1.0, h=0.5, dt=0.1, n_steps=3)
print(qc_example.draw(output='text', fold=-1))
print(f"Depth: {qc_example.depth()}, Gate count: {qc_example.size()}")

The circuit grows linearly with the number of Trotter steps. For $n = 4$ qubits and 3 steps, you will see alternating layers of RZZ and RX gates.

Measuring Magnetization vs Time

We use the StatevectorEstimator to compute the expectation value of the total $Z$ magnetization, $M_z = \frac{1}{n}\sum_i \langle Z_i \rangle$, as a function of time.

def measure_magnetization(n, J, h, total_time, n_steps):
    """Compute <Mz> vs time via statevector simulation."""
    dt = total_time / n_steps
    estimator = StatevectorEstimator()

    # Build observable: average Z magnetization
    paulis = []
    for i in range(n):
        pauli_str = 'I' * (n - 1 - i) + 'Z' + 'I' * i
        paulis.append((pauli_str, 1.0 / n))
    Mz_op = SparsePauliOp.from_list(paulis)

    times = []
    magnetizations = []

    for step in range(1, n_steps + 1):
        qc = build_trotter_circuit(n, J, h, dt, step)
        job = estimator.run([(qc, Mz_op)])
        result = job.result()
        mz = result[0].data.evs
        times.append(step * dt)
        magnetizations.append(float(mz))

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

# Run for ferromagnetic and paramagnetic cases
times, mz_ferro = measure_magnetization(n_qubits, J, h=0.3, total_time=5.0, n_steps=50)
_, mz_para = measure_magnetization(n_qubits, J, h=2.0, total_time=5.0, n_steps=50)

plt.figure(figsize=(8, 4))
plt.plot(times, mz_ferro, label='h/J = 0.3 (ferromagnetic)')
plt.plot(times, mz_para, label='h/J = 2.0 (paramagnetic)')
plt.xlabel('Time (J=1 units)')
plt.ylabel('<Mz>')
plt.title('Magnetization Dynamics under TFIM')
plt.legend()
plt.tight_layout()
plt.savefig('tfim_dynamics.png', dpi=150)

Starting from the all-$|+\rangle$ state, the ferromagnetic case shows oscillations as the system explores its ordered phase, while the paramagnetic case quickly equilibrates near zero net magnetization.

Scanning the Phase Transition

The quantum phase transition at $h/J = 1$ is visible in the ground-state magnetization. We prepare the ground state approximately by imaginary time evolution (or simply use exact diagonalization) and scan $h/J$.

def exact_ground_state_mz(n, J, h):
    """Compute <Mz> in the ground state via exact diagonalization."""
    from scipy.sparse import csr_matrix, kron, eye
    from scipy.sparse.linalg import eigsh

    sx = csr_matrix(np.array([[0, 1], [1, 0]], dtype=complex) / 2)
    sz = csr_matrix(np.array([[1, 0], [0, -1]], dtype=complex) / 2)
    I2 = csr_matrix(eye(2))

    def kron_op(op, i, n):
        ops = [I2] * n
        ops[i] = op
        result = ops[0]
        for o in ops[1:]:
            result = kron(result, o, format='csr')
        return result

    H = csr_matrix((2**n, 2**n), dtype=complex)
    for i in range(n - 1):
        H -= J * 2 * kron_op(sz, i, n) @ kron_op(sz, i + 1, n) * 4
    for i in range(n):
        H -= h * 2 * kron_op(sx, i, n) * 2

    energies, vecs = eigsh(H, k=1, which='SA')
    psi = vecs[:, 0]

    Mz_total = sum(kron_op(sz, i, n) for i in range(n)) * (2 / n)
    mz = np.real(psi.conj() @ (Mz_total @ psi))
    return abs(mz)

h_values = np.linspace(0.1, 3.0, 30)
mz_exact = [exact_ground_state_mz(n_qubits, J, h) for h in h_values]

plt.figure(figsize=(7, 4))
plt.plot(h_values / J, mz_exact, 'o-', markersize=4)
plt.axvline(1.0, color='red', linestyle='--', label='Critical point h/J=1')
plt.xlabel('h / J')
plt.ylabel('|<Mz>| (ground state)')
plt.title('Phase Transition in the 1D TFIM (n=4)')
plt.legend()
plt.tight_layout()
plt.savefig('tfim_phase_transition.png', dpi=150)

The magnetization drops sharply near $h/J = 1$, with finite-size rounding smoothing the transition for small $n$. Larger system sizes sharpen the feature and push it closer to the thermodynamic critical point.

Comparing Trotter to Exact Evolution

To validate the Trotter circuit, compare the time-evolved magnetization against the exact result from matrix exponentiation:

from scipy.linalg import expm
from qiskit.quantum_info import Operator

def exact_time_evolution_mz(n, J, h, times):
    """Exact time evolution via matrix exponentiation."""
    from scipy.sparse import csr_matrix, kron, eye
    sx = np.array([[0, 1], [1, 0]]) / 2
    sz = np.array([[1, 0], [0, -1]]) / 2
    I2 = np.eye(2)

    def full_op(op, i, n):
        ops = [I2] * n
        ops[i] = op
        result = ops[0]
        for o in ops[1:]:
            result = np.kron(result, o)
        return result

    H = np.zeros((2**n, 2**n), dtype=complex)
    for i in range(n - 1):
        H -= J * 4 * full_op(sz, i, n) @ full_op(sz, i + 1, n)
    for i in range(n):
        H -= h * 2 * full_op(sx, i, n)

    psi0 = np.ones(2**n) / np.sqrt(2**n)  # |+>^n
    Mz = sum(full_op(sz, i, n) for i in range(n)) * (2 / n)

    mz_vals = []
    for t in times:
        U = expm(-1j * H * t)
        psi_t = U @ psi0
        mz_vals.append(np.real(psi_t.conj() @ Mz @ psi_t))
    return np.array(mz_vals)

h_sim = 0.5
times_check = np.linspace(0.1, 3.0, 30)
mz_trotter_fine, _ = zip(*[
    (measure_magnetization(n_qubits, J, h_sim, t, max(1, int(t / 0.1)))[1][-1], t)
    for t in times_check
])
mz_trotter_fine = np.array(mz_trotter_fine)
mz_exact_evol = exact_time_evolution_mz(n_qubits, J, h_sim, times_check)

plt.figure(figsize=(7, 4))
plt.plot(times_check, mz_exact_evol, label='Exact', linewidth=2)
plt.plot(times_check, mz_trotter_fine, 'o', markersize=4, label='Trotter (dt=0.1)')
plt.xlabel('Time')
plt.ylabel('<Mz>')
plt.title('Trotter vs Exact for h/J=0.5')
plt.legend()
plt.tight_layout()
plt.savefig('tfim_trotter_vs_exact.png', dpi=150)

Key Takeaways

The 1D TFIM demonstrates several foundational ideas in quantum simulation:

Trotterization converts a continuous Hamiltonian evolution into a sequence of native two-qubit and single-qubit gates. The gate count grows as $O(n \cdot r)$, making it tractable on near-term devices for small $r$.

Phase transitions are visible even at small system sizes. The magnetization order parameter drops near $h/J = 1$, with finite-size effects rounding the sharp transition expected in the thermodynamic limit.

Validation via exact diagonalization is essential. For $n \leq 20$ qubits, numpy and scipy can represent the full $2^n \times 2^n$ Hamiltonian, giving a reference to verify Trotter accuracy.

Circuit depth vs fidelity trade-off: finer Trotter steps (smaller $\Delta t$) reduce algorithmic error but increase circuit depth and therefore physical noise on hardware. For real device execution, you must balance these competing factors by choosing the fewest steps that still capture the physics of interest.

To extend this tutorial, try implementing second-order Trotter (the Suzuki-Trotter formula), adding periodic boundary conditions (connecting site $n-1$ back to site $0$), or mapping the model to a free-fermion chain via the Jordan-Wigner transformation and comparing the resulting spectrum to your Trotter simulation.