Benchmarks for Many-Body Calculations

Nine small problems give numerical answers against which to compare a many-body calculation: complete spin and particle spectra, ideal-gas integrals, a BCS gap, and an exactly known finite-size sequence. Matching several quantities helps locate an error: an energy alone can conceal an incorrect basis, while degeneracies, trace moments, and observables test different parts of the calculation.

Required background. The computational overview explains what each numerical method approximates. Validation tests explains residuals and convergence checks.

The examples use their stated units, sectors, and boundaries. Passing them establishes agreement for these particular problems; a new geometry, interaction, or observable may require a new comparison. The runnable examples implement the open Ising chain and finite-cutoff BCS equation.

Comparing finite-cluster spectra and gas integrals

For a scalar observable $O$, record an absolute error and a scale-aware error,

$$\begin{aligned} \epsilon_{\mathrm{abs}} &= \left| O_{\mathrm{num}}-O_{\mathrm{ref}} \right|, \\ \epsilon_{\mathrm{scale}} &= \frac{ \epsilon_{\mathrm{abs}} }{ \max\!\left( |O_{\mathrm{ref}}|, O_{\mathrm{scale}} \right) }. \end{aligned}$$

The explicit $O_{\mathrm{scale}}$ avoids meaningless relative errors when the reference value is zero. For a normalized approximate eigenpair, report a residual such as

$$r_E = \frac{ \left\| H|\widetilde\psi\rangle - \widetilde E|\widetilde\psi\rangle \right\|_2 }{E_{\mathrm{scale}}}.$$

No one tolerance is appropriate for every method. The following reference profile is suitable for the tiny dense matrices below when they are assembled and diagonalized in IEEE double precision:

Quantity	Illustrative pass threshold
relative Hermiticity defect $\lVert H-H^\dagger\rVert_F/\max(\lVert H\rVert_F,1)$	$10^{-13}$
listed eigenvalue error in units of the stated coupling	$10^{-11}$
normalized trace-moment error	$10^{-12}$
listed ground-state observable error	$10^{-10}$
normalized eigenpair residual	$10^{-11}$

Sparse iterative solvers, stochastic estimates, arbitrary precision, and thermodynamic extrapolations require their own tolerances. A looser threshold can be legitimate, but the reason must be part of the calculation rather than chosen to hide a discrepancy.

Suite at a Glance

ID	Problem	Primary implementation layer	Trusted anchor
MB-B001	transverse-field Ising chain	spin-bit assembly and boundary bonds	exact $L=2$ spectrum and $L=4$ trace moments
MB-B002	spin-$1/2$ Heisenberg ring	spin normalization and multiplets	complete $L=4$ spectrum
MB-B003	Hubbard dimer	fermionic signs and spin sectors	analytic six-state spectrum
MB-B004	four-site Hubbard chain	many-fermion basis and observables	independently reproducible 36-state block
MB-B005	two-site Bose–Hubbard model	bosonic ladder factors	analytic three-state spectrum
MB-B006	ideal Bose gas	Bose functions and condensation branch	polylogarithm root and condensate fraction
MB-B007	ideal Fermi gas	Fermi integrals and Sommerfeld limit	zero-temperature values and low-$T$ coefficients
MB-B008	BCS gap equation	nonlinear self-consistency	exact zero-temperature gap
MB-B009	critical Ising gap scaling	size sequence and extrapolation	exact open-chain gap for every $L$

The order is useful. Run algebraic finite-cluster tests before interpreting a large simulation. A sophisticated extrapolation cannot repair a basis or sign error.

MB-B001: Transverse-Field Ising Chain

Hamiltonian and parameters

Use Pauli matrices with eigenvalues $\pm1$ and

$$H_{\mathrm{OBC}} = -J\sum_{i=1}^{L-1} \sigma_i^z\sigma_{i+1}^z -h\sum_{i=1}^{L}\sigma_i^x.$$

The computational basis obeys

$$\sigma^z|0\rangle=|0\rangle, \qquad \sigma^z|1\rangle=-|1\rangle.$$

For the periodic benchmarks take $L\geq3$ and add the bond $-J\sigma_L^z\sigma_1^z$ once. The two-site example below is open: a periodic two-edge ring and a single undirected dimer have different bond multiplicities. No constant shift is applied.

The primary analytic case is $L=2$ with open boundaries. Its four eigenvalues are

$$\left\{ -\sqrt{J^2+4h^2}, -J, +J, +\sqrt{J^2+4h^2} \right\}.$$

At $J=1$ and $h=0.7$, the sorted spectrum is

$$E/J = \{-1.720465053409, -1, 1, 1.720465053409\}.$$

Stronger assembly check

Use $L=4$, periodic boundaries, $J=1$, and $h=0.7$. The distinct energy levels and degeneracies are

$E/J$	degeneracy
$-4.563856065203$	1
$-4.441311123147$	1
$-1.735286090565$	1
$-1.400000000000$	2
$-0.441311123147$	1
$0$	4
$0.441311123147$	1
$1.400000000000$	2
$1.735286090565$	1
$4.441311123147$	1
$4.563856065203$	1

The full Hilbert-space dimension is $16$, and two useful basis-independent moments are

$$\begin{aligned} \operatorname{Tr}H &=0, \\ \frac{1}{16}\operatorname{Tr}H^2 &= 4(J^2+h^2) \\ &=5.96. \end{aligned}$$

Calculations to compare

1. Confirm the dimension $2^L$ and Hermiticity.
2. Reproduce the $L=2$ spectrum before adding periodic boundaries.
3. For $L=4$, reproduce the degeneracy sum, ground energy, trace, and second moment.
4. Verify the global spin-flip symmetry $P=\prod_i\sigma_i^x$ through $[H,P]=0$.
5. If symmetry blocks are used, reconstruct the complete spectrum from all blocks.

Interpreting a discrepancy

• Substituting $S^\alpha=\sigma^\alpha/2$ without changing the couplings reduces the interaction by four and the field by two; it does not rescale this whole Hamiltonian by one common factor.
• Mixing a two-edge periodic ring with a one-bond dimer changes the physical Hamiltonian. Use the open two-site and periodic four-site examples as specified.
• A correct ground energy with a wrong trace moment points to incomplete basis enumeration or a hidden energy shift.
• Missing degeneracies can indicate accidental symmetry breaking in the matrix assembly.

Passing MB-B001 validates these finite-chain conventions. It does not validate the thermodynamic critical point, long-distance correlations, or a time-evolution routine.

MB-B002: Spin-1/2 Heisenberg Ring

Hamiltonian and parameters

Use dimensionless spin operators $\mathbf S_i=\boldsymbol\sigma_i/2$ and the four-site periodic Hamiltonian

$$H = J\sum_{i=1}^{4} \mathbf S_i\mathbin{\cdot}\mathbf S_{i+1}, \qquad \mathbf S_5\equiv\mathbf S_1.$$

The four undirected bonds are $(1,2)$, $(2,3)$, $(3,4)$, and $(4,1)$. Set $J=1$ and apply no constant offset.

The complete $16$-state spectrum is

$E/J$	degeneracy	diagnostic content
$-2$	1	singlet ground state
$-1$	3	spin-one multiplet
$0$	7	remaining singlet and triplet states
$+1$	5	fully symmetric spin-two multiplet

The trace moments are

$$\operatorname{Tr}H=0, \qquad \frac{1}{16}\operatorname{Tr}H^2 = \frac{3}{4}J^2.$$

The fully polarized state supplies a one-line bond check: every bond has $\mathbf S_i\cdot\mathbf S_{i+1}=1/4$, so its total energy is $J$.

Calculations to compare

1. Reproduce all four distinct energies and their degeneracies.
2. Check $[H,S^z_{\mathrm{tot}}]=0$ and $[H,\mathbf S_{\mathrm{tot}}^2]=0$.
3. Verify that states grouped by total spin are degenerate to the declared tolerance.
4. Reconstruct the full trace moments from any $S^z_{\mathrm{tot}}$ blocks used in the calculation.
5. Compare the lowest energies obtained in the full basis and in the $S^z_{\mathrm{tot}}=0$ block.

Interpreting a discrepancy

A spectrum larger by a factor of four almost always signals the use of Pauli matrices where spin operators were intended. An isolated error in the polarized energy suggests a bond-counting problem. Correct eigenvalues but incorrect multiplet labels point to a faulty total-spin operator or an inconsistent tensor-factor ordering.

The thermodynamic Bethe-ansatz energy density,

$$\lim_{L\to\infty} \frac{E_0(L)}{LJ} = \frac14-\ln2,$$

is a valuable later test, but it is not part of the tiny-ring acceptance condition. The canonical Heisenberg page owns its interpretation.

MB-B003: Hubbard Dimer

The Hubbard dimer explains the singlet and triplet states behind this spectrum.

Hamiltonian and parameters

Use two sites, two spin species, and

$$H ={} -t\sum_{\sigma} \left( c_{1\sigma}^\dagger c_{2\sigma} + c_{2\sigma}^\dagger c_{1\sigma} \right) + U\sum_{i=1}^{2} n_{i\uparrow}n_{i\downarrow}.$$

Fix the fermionic mode order to

$$(1\uparrow,1\downarrow,2\uparrow,2\downarrow)$$

and benchmark the total-particle sector $N=2$, whose dimension is $\binom42=6$. The $N_\uparrow=N_\downarrow=1$ block has dimension $4$ and contains the $S^z=0$ member of the triplet.

With suitable phases for the singlet and even/odd doublon states, the singlet block is

$$H_{S=0} = \begin{pmatrix} 0 & -2t & 0\\ -2t & U & 0\\ 0 & 0 & U \end{pmatrix}.$$

The full $N=2$ spectrum is

$$\begin{gathered} E_- = \frac{U-\sqrt{U^2+16t^2}}{2}, \\ E_T=0 \quad\text{with degeneracy }3, \\ E_D=U, \qquad E_+ = \frac{U+\sqrt{U^2+16t^2}}{2}. \end{gathered}$$

At $t=1$ and $U=4$,

$$\begin{aligned} E_0/t &=-0.828427124746190, \\ \Delta_{S-T}/t &=0.828427124746190, \\ \langle D\rangle_0 &=0.146446609406726, \end{aligned}$$

where

$$D = \sum_i n_{i\uparrow}n_{i\downarrow}.$$

The last target follows independently from Feynman–Hellmann differentiation,

$$\langle D\rangle_0 = \frac{\partial E_-}{\partial U} = \frac12 \left( 1- \frac{U}{\sqrt{U^2+16t^2}} \right).$$

Calculations to compare

1. Confirm the dimensions $6$ for fixed $N=2$ and $4$ for fixed $N_\uparrow=N_\downarrow=1$.
2. Reproduce the complete analytic spectrum in the full sector.
3. Verify the three triplet states are degenerate at zero energy.
4. Compare direct evaluation of $\langle D\rangle_0$ with $\partial E_-/\partial U$.
5. Repeat after changing the mode ordering and confirm that observables and eigenvalues remain invariant when all fermionic signs are transformed consistently.

Interpreting a discrepancy

The dimer is especially sensitive to fermionic sign conventions. Incorrect singlet mixing, a split triplet, or a missing factor of two in the square root usually identifies an inconsistent hopping sign or an incomplete spin sum. Agreement only after taking absolute values of matrix entries is not a pass.

MB-B004: Four-Site Hubbard Chain

The Hubbard chain connects this finite open chain with the bulk model.

Hamiltonian and parameters

Use the same Hamiltonian and site-major mode ordering as MB-B003, now on an open chain of $L=4$ sites. The bonds are $(1,2)$, $(2,3)$, and $(3,4)$. Fix

$$t=1, \qquad U=4, \qquad N_\uparrow=N_\downarrow=2.$$

The symmetry-block dimension is

$$\dim\mathcal H_{2,2} = \binom42\binom42 =36.$$

This benchmark is deliberately larger than the dimer but still small enough for complete dense diagonalization. The trusted ground-state data are

$$\begin{aligned} E_0 &=-1.953145308684556, \\ E_1 &=-1.412898695919094, \\ E_1-E_0 &=0.540246612765462, \\ \langle D\rangle_0 &=0.339585650339361. \end{aligned}$$

Here $D=\sum_i n_{i\uparrow}n_{i\downarrow}$. The energy decomposition provides two additional checks:

$$\begin{aligned} E_{\mathrm{int}} &=U\langle D\rangle_0, \\ &=1.358342601357445, \\ E_{\mathrm{kin}} &=E_0-E_{\mathrm{int}}, \\ &=-3.311487910042001. \end{aligned}$$

Over the complete $36$-state block,

$$\begin{aligned} \operatorname{Tr}H &=144, \\ \frac{1}{36}\operatorname{Tr}H^2 &=25.3333333333333. \end{aligned}$$

Calculations to compare

1. Verify the combinatorial dimension before constructing the Hamiltonian.
2. Check Hermiticity and closure of every hopping transition inside the fixed particle-number block.
3. Reproduce $E_0$, $E_1$, the trace, and the second moment.
4. Compute $\langle D\rangle_0$ directly and verify $E_0=E_{\mathrm{kin}}+U\langle D\rangle_0$.
5. Cross-check the result with two independent representations when practical, such as a generic occupation-bit implementation and factorized up/down bit strings.

The quoted $E_1-E_0$ is the first excitation gap within the fixed $(N_\uparrow,N_\downarrow)=(2,2)$ block. It is not, by itself, the thermodynamic charge gap, spin gap, or unrestricted global gap. Those quantities require specified neighboring sectors and a size sequence.

Interpreting a discrepancy

• Dimension or closure failures indicate sector enumeration errors.
• A correct trace but wrong low spectrum often indicates hopping signs or bond orientation.
• A correct energy with wrong double occupancy suggests an observable-basis mismatch.
• Agreement between two solvers using the same matrix does not independently test matrix assembly; an independent construction also tests the matrix elements.

MB-B005: Two-Site Bose–Hubbard Model

Hamiltonian and parameters

Use

$$H ={} -J \left( b_1^\dagger b_2 + b_2^\dagger b_1 \right) + \frac{U}{2} \sum_{i=1}^{2} n_i(n_i-1).$$

in the fixed-$N=2$ basis

$$\bigl\{ |2,0\rangle, |1,1\rangle, |0,2\rangle \bigr\}.$$

The matrix is

$$H = \begin{pmatrix} U & -\sqrt2J & 0\\ -\sqrt2J & 0 & -\sqrt2J\\ 0 & -\sqrt2J & U \end{pmatrix},$$

with spectrum

$$\begin{aligned} E_- &= \frac{U-\sqrt{U^2+16J^2}}{2}, \\ E_D &=U, \\ E_+ &= \frac{U+\sqrt{U^2+16J^2}}{2}. \end{aligned}$$

At $J=1$ and $U=2$,

$$\begin{aligned} E_-/J &=1-\sqrt5, \\ E_D/J &=2, \\ E_+/J &=1+\sqrt5. \end{aligned}$$

or numerically

$$\begin{aligned} E_-/J &=-1.236067977499790, \\ E_D/J &=2, \\ E_+/J &=3.236067977499790. \end{aligned}$$

Define the total on-site pair indicator

$$D_B = \frac12\sum_i n_i(n_i-1).$$

Its ground-state expectation is

$$\begin{aligned} \langle D_B\rangle_0 &= \frac12 \left( 1- \frac{U}{\sqrt{U^2+16J^2}} \right) \\ &=0.276393202250021. \end{aligned}$$

at the benchmark point.

Calculations to compare

1. Confirm that the fixed-number basis has dimension $3$.
2. Generate hopping amplitudes from ladder operators rather than hard-coding them.
3. Reproduce the full spectrum and the reflection parity of all three eigenstates.
4. Verify $\langle D_B\rangle_0=\partial E_0/\partial U$.
5. Check the noninteracting limit $U=0$, where the two bosons occupy the bonding orbital and $E_0=-2J$.

Interpreting a discrepancy

Off-diagonal entries of magnitude $J$ instead of $\sqrt2J$ reveal missing bosonic enhancement factors. An interaction diagonal of $2U$ on $|2,0\rangle$ indicates that the factor $1/2$ in $n(n-1)$ was omitted. The Bose–Hubbard Dimer dossier owns the exact finite-model interpretation and additional invariant checks; the general sparse construction lives in Occupation-Number Representation.

MB-B006: Ideal Bose Gas

Hamiltonian and parameters

Use a uniform, spinless, three-dimensional ideal Bose gas in the thermodynamic limit. Define

$$\lambda_T = \sqrt{ \frac{2\pi\hbar^2}{mk_{\mathrm B}T} }, \qquad z=e^{\beta\mu},$$

and, above the condensation point,

$$n\lambda_T^3 = \operatorname{Li}_{3/2}(z), \qquad 0<z<1.$$

The implementation should reproduce

$$\begin{aligned} \zeta(3/2) &=2.612375348685488\ldots, \\ \zeta(5/2) &=1.341487257250917\ldots. \end{aligned}$$

For the normal-state target

$$n\lambda_T^3=1,$$

the physical root is

$$z=0.698614359135065.$$

For a second target below the transition, set $T/T_c=1/2$. The thermodynamic condensate fraction is

$$\begin{aligned} \frac{N_0}{N} &= 1- \left( \frac{T}{T_c} \right)^{3/2} \\ &=0.646446609406726. \end{aligned}$$

Calculations to compare

1. Reproduce the two zeta values independently of the root solver.
2. Solve the normal-state number equation on the physical interval $0<z<1$ and report its residual.
3. Refine quadrature order or arithmetic precision until the quoted digits stabilize.
4. Below $T_c$, separate the ground mode and pin the excited-state fugacity to its thermodynamic limiting value $z=1$.
5. Verify that the excited fraction scales as $(T/T_c)^{3/2}$ for several temperatures below $T_c$.

Interpreting a discrepancy

A root $z>1$ is unphysical for the continuum ideal Bose gas and usually reflects unconstrained root finding. Failure below $T_c$ often occurs when the code keeps trying to satisfy the entire density with excited states. In a finite box the chemical potential remains below the ground-state energy; the exact $z=1$ prescription is a thermodynamic-limit problem, not a finite-volume identity.

MB-B007: Ideal Fermi Gas

Hamiltonian and parameters

Use a uniform three-dimensional gas with spin degeneracy $g=2$ and units

$$\frac{\hbar^2}{2m}=1.$$

Choose

$$n=\frac{1}{3\pi^2},$$

so that $k_F=1$ and $\epsilon_F=1$. With

$$f(\epsilon) = \frac{1}{e^{(\epsilon-\mu)/(k_{\mathrm B}T)}+1},$$

the number and energy densities can be evaluated directly as

$$\begin{aligned} n &= \frac{1}{\pi^2} \int_0^\infty k^2 f(k^2)\,dk, \\ u &= \frac{1}{\pi^2} \int_0^\infty k^4 f(k^2)\,dk. \end{aligned}$$

At zero temperature the exact targets are

$$\begin{aligned} \mu &=\epsilon_F=1, \\ \frac{U}{N} &=\frac35\epsilon_F=0.6, \\ \frac{P}{n\epsilon_F} &=\frac25. \end{aligned}$$

For $\theta=T/T_F\to0$, the Sommerfeld coefficients provide a refinement benchmark:

$$\begin{aligned} \frac{\mu}{\epsilon_F} &= 1- \frac{\pi^2}{12}\theta^2 +O(\theta^4), \\ \frac{U}{N\epsilon_F} &= \frac35 \left[ 1+ \frac{5\pi^2}{12}\theta^2 +O(\theta^4) \right], \\ \frac{C_V}{Nk_{\mathrm B}} &= \frac{\pi^2}{2}\theta +O(\theta^3). \end{aligned}$$

Calculations to compare

1. Reproduce the $T=0$ density, energy per particle, and pressure ratio.
2. At finite $T$, solve the number equation for $\mu$ rather than holding $\mu=\epsilon_F$.
3. Use a decreasing sequence such as $\theta=0.20,0.10,0.05,0.025$ and monitor the scaled corrections.
4. Verify the limits

$$\begin{aligned} \frac{1-\mu/\epsilon_F}{\theta^2} &\longrightarrow \frac{\pi^2}{12}, \\ \frac{ U/(N\epsilon_F)-3/5 }{\theta^2} &\longrightarrow \frac{\pi^2}{4}, \\ \frac{C_V}{Nk_{\mathrm B}\theta} &\longrightarrow \frac{\pi^2}{2}. \end{aligned}$$

5. Record quadrature truncation and root residual separately from the asymptotic $O(\theta^4)$ or $O(\theta^3)$ error.

The low-temperature formulas are convergence laws, not exact finite-$\theta$ targets. Demanding twelve-digit agreement at $\theta=0.1$ would confuse truncation of the Sommerfeld series with numerical error.

Interpreting a discrepancy

An incorrect zero-temperature density commonly signals a missing spin degeneracy or phase-space factor. A chemical potential that stays fixed at finite temperature violates the fixed-density problem. Unstable heat capacity estimates may come from differentiating noisy energy data; analytic derivatives, higher precision, or a fitted low-temperature expansion can separate differentiation error from quadrature error.

MB-B008: Mean-Field BCS Gap Equation

Hamiltonian and parameters

Use a constant density of states within a symmetric energy cutoff $|\xi|\leq\omega_c$. Let $\nu(0)$ be the density of states per spin in this pairing convention and $g>0$ the attractive coupling. The cutoff $\omega_c$ here has energy units, despite its frequency-like symbol. Choose

$$\lambda=g\nu(0)=0.25, \qquad \omega_c=1.$$

The nonzero mean-field gap satisfies

$$\frac{1}{\lambda} = \int_0^{\omega_c} \frac{d\xi}{\sqrt{\xi^2+\Delta^2}} \tanh\!\left( \frac{\sqrt{\xi^2+\Delta^2}}{2k_{\mathrm B}T} \right).$$

At $T=0$, the integral is elementary and gives the exact finite-cutoff target

$$\frac{\Delta_0}{\omega_c} = \frac{1}{\sinh(1/\lambda)} = 0.0366435703258656.$$

The weak-coupling approximation is

$$\frac{\Delta_0}{\omega_c} \approx 2e^{-1/\lambda} = 0.0366312777774684,$$

whose relative error at this parameter point is approximately $3.35\times10^{-4}$. A solver should reproduce the exact finite-cutoff result, not be judged against the asymptotic expression at more precision than the approximation warrants.

The transition temperature follows from the linearized equation

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

For the benchmark parameters,

$$\frac{k_{\mathrm B}T_c}{\omega_c} = 0.0207674786897134,$$

and therefore

$$\frac{2\Delta_0}{k_{\mathrm B}T_c} = 3.52893780447368.$$

The slight difference from the asymptotic weak-coupling value $2\pi/e^\gamma\approx3.52775$ is primarily the finite-cutoff correction retained in $\Delta_0$.

Calculations to compare

1. Evaluate the $T=0$ integral and reproduce the exact nonzero root.
2. Report both the gap error and the residual of the integral equation.
3. Bracket the physical root on $\Delta>0$; do not let the trivial normal solution masquerade as convergence below $T_c$.
4. Solve the linearized equation for $T_c$ and compare the finite-cutoff ratio with the weak-coupling limit.
5. Check that the superconducting stationary point has lower mean-field free energy than the normal state for $T<T_c$.
6. If the calculation is performed at fixed density rather than fixed chemical potential, solve and report the number equation as an additional problem.

Interpreting a discrepancy

False convergence to $\Delta=0$ usually reflects an unbracketed nonlinear solve or use of the undivided stationarity equation without phase selection. A correct root with a large integral residual suggests cancellation or quadrature error. Agreement with $2e^{-1/\lambda}$ but not $1/\sinh(1/\lambda)$ can simply mean the implementation silently replaced the finite cutoff by the weak-coupling approximation.

MB-B009: Finite-Size Scaling of the Ising Gap

Hamiltonian and parameters

Return to the open transverse-field Ising chain of MB-B001 at its critical point,

$$h=J>0.$$

The gap between the two lowest states of the complete finite Hilbert space is

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

Useful exact values are

$L$	$\Delta_L/J$
2	1.236067977499790
3	0.890083735825258
4	0.694592710667721
5	0.569259353093141
6	0.482146721021292
7	0.418113853070614
8	0.369073437853208

The boundary-shifted scaling variable exposes the known correction structure. With $\ell=L+1/2$,

$$\begin{aligned} \frac{\ell\Delta_L}{\pi J} &= \frac{\sin x}{x}, \\ x &= \frac{\pi}{4\ell}, \end{aligned}$$

so that

$$\frac{\ell\Delta_L}{\pi J} = 1- \frac{\pi^2}{96\ell^2} +O(\ell^{-4}).$$

This single problem tests both spectral extraction and the scaling analysis. The leading critical law is

$$\Delta_L\sim\frac{\pi J}{L},$$

corresponding to dynamical exponent $z=1$, while the exact formula reveals how open-boundary corrections approach that limit.

Calculations to compare

1. Reproduce the exact gaps for every size used in the fit.
2. State whether the gap is global or restricted to a symmetry sector. This problem uses the global gap.
3. Verify that $\Delta_L$ decreases monotonically over the benchmark sequence.
4. Fit at least two lower-size cutoffs $L_{\min}$ and report the drift of the inferred exponent and amplitude.
5. Compare an unshifted fit in $L$ with the boundary-aware variable $\ell=L+1/2$.
6. Test the normalized residual

$$R_L = \frac{\Delta_L}{ 4J\sin[\pi/(4L+2)] }-1$$

before interpreting any extrapolation.

Interpreting a discrepancy

An exact-formula residual at small $L$ is a spectrum, boundary, or sector error rather than a scaling uncertainty. A good log–log line with drifting exponent can result from omitted boundary corrections. Using a parity-restricted gap without saying so changes the physical excitation and invalidates comparison with this problem.

Passing the fit only confirms that the analysis can recover a known $z=1$ sequence over the tested window. It does not validate an ansatz for a different universality class.

Exercises

Exercise 1: The factor-of-four failure

A four-site Heisenberg implementation returns the levels $-8J$, $-4J$, $0$, and $4J$ with the correct degeneracies. Identify the likely error and give a direct one-state test.

Solution

The spectrum is four times the target, so the code probably used Pauli matrices $\boldsymbol\sigma$ in

$$J\sum_i \boldsymbol\sigma_i\mathbin{\cdot}\boldsymbol\sigma_{i+1}$$

instead of spin operators $\mathbf S_i=\boldsymbol\sigma_i/2$. Test the fully polarized state. Each intended bond contributes $J/4$, and the four-site ring must therefore have energy $J$. A result $4J$ confirms the normalization error.

Exercise 3: Why the sector gap is not the charge gap

Explain why $E_1-E_0$ in the four-site $(N_\uparrow,N_\downarrow)=(2,2)$ Hubbard block is not a charge gap. Write a finite-system expression that probes particle addition and removal instead.

Solution

Both states in the quoted difference have the same particle numbers, so the excitation does not test the energy cost of changing the total charge. A common finite-system charge-gap estimator at total particle number $N$ is

$$\Delta_c(N) ={} E_0(N+1) + E_0(N-1) -2E_0(N).$$

where $E_0(N)$ is minimized over all allowed spin sectors at that total particle number. Restricting the spin sectors defines a different addition/removal energy. Other conventions use pair addition and removal to preserve spin balance. In either case, neighboring particle-number sectors are essential.

Exercise 5: Boundary-aware Ising scaling

Starting from the exact critical gap, show that the shifted length $\ell=L+1/2$ removes the leading $1/L$ boundary correction from the normalized amplitude.

Solution

Write

$$\Delta_L = 4J\sin\!\left( \frac{\pi}{4\ell} \right), \qquad \ell=L+\frac12.$$

With $x=\pi/(4\ell)$,

$$\frac{\ell\Delta_L}{\pi J} = \frac{\sin x}{x}.$$

Expanding gives

$$\begin{aligned} \frac{\sin x}{x} &= 1- \frac{x^2}{6} +O(x^4) \\ &= 1- \frac{\pi^2}{96\ell^2} +O(\ell^{-4}). \end{aligned}$$

There is no term proportional to $1/\ell$. If one normalizes with $L$ instead, expanding $L/(L+1/2)$ reintroduces a leading $1/L$ correction.

Exercise 6: Design a fair low-temperature test

Why should the ideal-Fermi-gas value at one temperature not be compared directly with the truncated Sommerfeld formula at machine precision? Propose a better acceptance test.

Solution

The displayed Sommerfeld formulas omit terms of order $\theta^4$ in $\mu$ and $U/N$, and order $\theta^3$ in $C_V$. Even exact quadrature should therefore differ from the truncated formula at finite $\theta$.

A better test uses a decreasing temperature sequence, solves the number equation at each point, and examines scaled corrections such as

$$a_\mu(\theta) = \frac{1-\mu(\theta)/\epsilon_F}{\theta^2}.$$

The sequence should approach $\pi^2/12$ with drift consistent with $O(\theta^2)$. Quadrature and root residuals should be much smaller than that expected asymptotic drift and reported separately.

References

• J. P. F. LeBlanc et al., “Solutions of the two-dimensional Hubbard model: benchmarks and results from a wide range of numerical algorithms,” Physical Review X 5, 041041 (2015).
• J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of superconductivity,” Physical Review 108, 1175–1204 (1957).
• R. K. Pathria and P. D. Beale, Statistical Mechanics, 3rd ed., Elsevier (2011).
• H. Bethe, “Zur Theorie der Metalle. I. Eigenwerte und Eigenfunktionen der linearen Atomkette,” Zeitschrift für Physik 71, 205–226 (1931).
• E. Dagotto, “Correlated electrons in high-temperature superconductors,” Reviews of Modern Physics 66, 763–840 (1994).
• M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, “Boson localization and the superfluid–insulator transition,” Physical Review B 40, 546–570 (1989).
• K. Huang, Statistical Mechanics, 2nd ed., Wiley (1987).
• J. Hubbard, “Electron correlations in narrow energy bands,” Proceedings of the Royal Society A 276, 238–257 (1963).
• L. N. Trefethen and D. Bau III, Numerical Linear Algebra, SIAM (1997).
• E. Lieb, T. Schultz, and D. Mattis, “Two soluble models of an antiferromagnetic chain,” Annals of Physics 16, 407–466 (1961).
• P. Pfeuty, “The one-dimensional Ising model with a transverse field,” Annals of Physics 57, 79–90 (1970).
• A. W. Sandvik, “Computational studies of quantum spin systems,” AIP Conference Proceedings 1297, 135–338 (2010).
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press (2007).
• E. H. Lieb and F. Y. Wu, “Absence of Mott transition in an exact solution of the short-range, one-band model in one dimension,” Physical Review Letters 20, 1445–1448 (1968).