Reproducible Many-Body Notebooks

A small calculation is easiest to understand when its Hamiltonian, code, and expected answer fit together. The examples below can be pasted into notebook cells or saved as a Python script. They calculate an open Ising chain and a finite-cutoff BCS gap using only Python and NumPy. No separate Many-Body notebook download is required.

Required background. The computational overview explains finite-system calculations, and benchmark problems supplies the reference answers. Reproducibility status explains how to interpret downloadable notebook results elsewhere in the Library.

An open Ising chain in a few lines

Use Pauli matrices, $J=1$, and open boundaries:

$$H=-J\sum_{i=0}^{L-2}Z_iZ_{i+1}-h\sum_{i=0}^{L-1}X_i.$$

The integer $s$ labels the product state with bit $i$ equal to zero for $Z_i=+1$ and one for $Z_i=-1$. Flipping bit $i$ implements $X_i$. This makes the diagonal interaction and off-diagonal field visible separately in the code. The dense matrix grows as $4^L$ entries, so keep these examples small. See exact diagonalization for larger sparse calculations.

Run these cells in order in a fresh Python process or restarted notebook kernel. The examples were tested with Python 3.13 and NumPy 2.1.2; the first lines print the versions actually in use.

import sys
import numpy as np

print(sys.version.split()[0], np.__version__)

def ising_open(L, J=1.0, h=1.0):
    if L < 2:
        raise ValueError("Use at least two sites")
    H = np.zeros((2**L, 2**L))
    for s in range(2**L):
        z = [1 - 2*((s >> i) & 1) for i in range(L)]
        H[s, s] = -J * sum(z[i]*z[i+1] for i in range(L-1))
        for i in range(L):
            H[s ^ (1 << i), s] -= h
    return H

H2 = ising_open(2, h=0.7)
E2 = np.linalg.eigvalsh(H2)
target = np.array([-np.sqrt(2.96), -1.0, 1.0, np.sqrt(2.96)])
np.testing.assert_allclose(H2, H2.T, atol=1e-14)
np.testing.assert_allclose(E2, target, atol=1e-12, rtol=0)
print("Two-site energies:", E2)

The four energies are approximately $(-1.72046505,-1,1,1.72046505)$. This comparison tests the matrix construction as well as diagonalization. Merely verifying $H\psi=E\psi$ would only test an eigenpair of whatever matrix was assembled.

Recovering a finite-size gap

At $h=J=1$, the exact global gap of the open chain is

$$\Delta_L=4\sin\!\left(\frac{\pi}{4L+2}\right).$$

Diagonalize the complete Hilbert space here; keeping only one parity sector would change which gap you calculate. The following cell implements MB-B009.

for L in range(2, 8):
    E = np.linalg.eigvalsh(ising_open(L))
    gap = E[1] - E[0]
    exact = 4*np.sin(np.pi/(4*L + 2))
    np.testing.assert_allclose(gap, exact, atol=2e-12, rtol=0)
    scaled = (L + 0.5)*gap/np.pi
    print(f"L={L}: gap={gap:.10f}, scaled={scaled:.10f}")

The scaled quantity approaches one. Its remaining correction is $-\pi^2/[96(L+1/2)^2]$ to leading order. Thus an imperfect straight line in a log–log plot need not indicate an incorrect eigensolver: the finite-size correction is physical.

Thermal averages and quench dynamics from one spectrum

For a temperature-independent Hamiltonian, the Gibbs weights give

$$U=\sum_n p_nE_n,\qquad \frac{C_V}{k_{\mathrm B}}=\beta^2\sum_n p_n(E_n-U)^2.$$

Subtracting $E_0$ before exponentiating avoids overflow without changing normalized probabilities. If the partition function itself is needed, restore the shift through $\ln Z=-\beta E_0+\ln\sum_n e^{-\beta(E_n-E_0)}$.

H = ising_open(4, h=0.7)
E, V = np.linalg.eigh(H)
beta = 1.3  # inverse energy, in units J=1
weights = np.exp(-beta*(E-E[0]))
p = weights/weights.sum()
U = p @ E
Cv_over_kB = beta**2 * (p @ ((E-U)**2))
rho = (V*p) @ V.conj().T
np.testing.assert_allclose(np.trace(rho), 1.0, atol=1e-13)
np.testing.assert_allclose(np.trace(rho @ H), U, atol=1e-12)
print("U, Cv/kB:", U, Cv_over_kB)

For a quench, prepare the ground state at $h_i=1.4$ and evolve with $h_f=0.7$. The columns of $V$ are eigenvectors of the final Hamiltonian. In units $\hbar=1$, each energy component acquires the phase $e^{-iE_nt}$.

_, Vi = np.linalg.eigh(ising_open(4, h=1.4))
psi0 = Vi[:, 0].astype(complex)
coeff = V.conj().T @ psi0
z0 = np.array([1 - 2*(s & 1) for s in range(16)])
z1 = np.array([1 - 2*((s >> 1) & 1) for s in range(16)])
energy0 = np.vdot(psi0, H @ psi0).real
for t in [0.0, 0.5, 1.0, 2.0]:
    psi = V @ (np.exp(-1j*E*t)*coeff)
    np.testing.assert_allclose(np.vdot(psi, psi), 1.0, atol=1e-12)
    np.testing.assert_allclose(np.vdot(psi, H @ psi), energy0, atol=1e-12)
    if t == 0:
        np.testing.assert_allclose(psi, psi0, atol=1e-12)
    zz = np.vdot(psi, z0*z1*psi).real
    print(f"t={t:.1f}: <Z0 Z1>={zz:.8f}")

The bond correlation changes while energy and norm remain constant. A parity-odd one-point magnetization would vanish for this symmetry-preserving initial state and Hamiltonian, making it a poor choice for demonstrating the dynamics. This four-site evolution can show recurrences; it cannot establish irreversible relaxation in an infinite chain.

Solving the finite-cutoff BCS equation

Take the dimensionless pairing strength $\lambda=0.25$, cutoff energy $E_c=1$, and $k_{\mathrm B}=1$. Here the density of states entering $\lambda$ is per spin. At zero temperature,

$$\frac1\lambda=\int_0^1\frac{d\xi}{\sqrt{\xi^2+\Delta_0^2}},\qquad \Delta_0=\frac1{\sinh(1/\lambda)}.$$

At the transition, the gap tends to zero and

$$\frac1\lambda=\int_0^1\frac{d\xi}{\xi}\tanh\!\left(\frac{\xi}{2T_c}\right).$$

The integrand has the finite limit $1/(2T_c)$ at the origin. Gauss–Legendre nodes lie inside the integration interval, so the following code never divides by zero. Bisection keeps the positive temperature root bracketed; increasing the quadrature order checks integration accuracy independently of the root tolerance.

from numpy.polynomial.legendre import leggauss

lam = 0.25

def bisect_root(f, lo, hi, steps=80):
    flo, fhi = f(lo), f(hi)
    if flo*fhi >= 0:
        raise ValueError("Endpoints must bracket the root")
    for _ in range(steps):
        mid = (lo + hi)/2
        fm = f(mid)
        if fm == 0:
            return mid
        if flo*fm > 0:
            lo, flo = mid, fm
        else:
            hi = mid
    return (lo + hi)/2

Delta0 = 1/np.sinh(1/lam)
for order in [64, 128, 256]:
    nodes, weights = leggauss(order)
    xi, weights = (nodes + 1)/2, weights/2
    def residual(T):
        return weights @ (np.tanh(xi/(2*T))/xi) - 1/lam
    Tc = bisect_root(residual, 0.001, 0.2)
    assert abs(residual(Tc)) < 1e-12
    np.testing.assert_allclose(Tc, 0.0207674786897134,
                               atol=2e-12, rtol=0)
    print(order, Delta0, Tc, 2*Delta0/Tc)

The answer is $\Delta_0\simeq0.0366435703$, $T_c\simeq0.0207674787$, and $2\Delta_0/T_c\simeq3.52893780$. The small departure from the weak-coupling ratio $2\pi/e^\gamma$ is primarily a finite-cutoff effect. It is not removed by tightening a numerical tolerance. See BCS mean-field theory for the approximation and its physical limitations.

Extending the examples

Change one physical or numerical parameter at a time and retain the printed values alongside any plots. For the chain, compare more lengths before fitting a bulk exponent. For the quench, compare another observable and a second method such as a converged short-time product formula. For the BCS integral, change quadrature order separately from the root bracket and stopping tolerance.

When adapting these cells to a notebook, restart and run all cells after each substantial change. Put parameters before their uses; record the interpreter and package versions as above. This makes another reader’s first run follow the same calculation as yours. For stochastic extensions, a seed alone is insufficient: also estimate autocorrelation and uncertainty, as explained in quantum Monte Carlo.

Exercises

Why shift the Gibbs energies?

Show that replacing $E_n$ by $E_n-E_0$ leaves the normalized probabilities unchanged. Would the same replacement leave the internal energy unchanged if used everywhere in the formula for $U$?

Solution

The common factor $e^{\beta E_0}$ cancels between every weight and their sum. However, averaging the shifted energies gives $U-E_0$. Use the original energies for $U$, or add $E_0$ back. The heat capacity is unchanged for a temperature-independent shift.

A phase error that conserves norm

Replace every evolution factor by $e^{-iE_nt/2}$ in the quench cell. Which displayed checks still pass, and what additional comparison detects the error?

Solution

Norm, final-Hamiltonian energy, and the initial state at $t=0$ still pass. The code now evolves to time $t/2$. Comparing the bond correlation with an independently converged propagator at the stated time detects the wrong phase scale. Conservation laws alone do not fix the speed of evolution.

Exact cutoff or asymptotic gap?

Derive the ratio between the exact zero-temperature gap and $2e^{-1/\lambda}$, and explain the sign of the difference.

Solution

Using $\sinh x=(e^x-e^{-x})/2$ gives

$$\frac{1/\sinh(1/\lambda)}{2e^{-1/\lambda}} =\frac1{1-e^{-2/\lambda}}>1.$$

The weak-coupling expression drops the positive finite-cutoff correction. Its relative error measured against the exact gap is $e^{-2/\lambda}$, approximately $3.35\times10^{-4}$ at $\lambda=0.25$.

References

• C. R. Harris et al., “Array programming with NumPy,” Nature 585, 357–362 (2020).
• A. Rule et al., “Ten simple rules for writing and sharing computational analyses in Jupyter Notebooks,” PLOS Computational Biology 15, e1007007 (2019).
• J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of superconductivity,” Physical Review 108, 1175–1204 (1957).
• P. Pfeuty, “The one-dimensional Ising model with a transverse field,” Annals of Physics 57, 79–90 (1970).