Test finite-chain Bethe solutions

Does a proposed Bethe solution actually produce an eigenstate of the finite spin chain? In this project you will reconstruct and normalize a two-magnon wavefunction, apply two independently built Hamiltonians, and compare its energy with direct diagonalization. The six-site example gives a concrete success; a reversed scattering ratio, a missing closing bond, and a zero-vector solution show what can go wrong.

Required background. Use the coordinate wavefunction and ring equations from Quantize magnons on a ring. You should be able to run Python and form a complex inner product; the entry repair reviews normalization and residuals.

Helpful background. The XXX course gives the sequence, the model defines its regime, and the two-magnon derivation explains the contact equation behind the formulas tested here.

The periodic XXX eigenstate test

Take $N\geq3$ spin-$1/2$ sites on a ring, with $S_i^\alpha=\sigma_i^\alpha/2$, $\hbar=1$, lattice spacing one, and ferromagnetic coupling $J\gt0$. The Hamiltonian and periodic identification are

$$\begin{aligned} H&=J\sum_{i=0}^{N-1}\left(\frac14-\mathbf S_i\mathbin{\cdot}\mathbf S_{i+1}\right) =\frac J2\sum_{i=0}^{N-1}(1-P_{i,i+1}),\\ \mathbf S_N&=\mathbf S_0. \end{aligned}$$

Here $P_{i,i+1}$ swaps two neighboring spins. The all-up state has energy zero. We index sites by $0,\ldots,N-1$ and let $M$ count down spins; the $M=2$ sector has dimension $\binom N2$. The two-site ring is excluded because it needs a separate convention for counting bonds.

For $x\lt y$, set $A_{12}=1$ and use

$$\begin{aligned} \psi_{\mathrm{raw}}(x,y)&=z_1^xz_2^y+S_{12}z_2^xz_1^y,\qquad z_j=e^{ik_j},\\ S_{12}&=-\frac{1+z_1z_2-2z_2}{1+z_1z_2-2z_1},\\ z_1^N&=S_{12}^{-1},\qquad z_2^N=S_{12}. \end{aligned}$$

This amplitude orientation matters. In Karbach and Müller’s notation the ratio $A/A'$ is the inverse of our $S_{12}=A_{21}/A_{12}$; their energy is shifted by the vacuum value $E_0=-JN/4$. See Karbach and Müller 1997, arXiv v1 (submitted 1998), pp. 1–3, equations (1), (9)–(18), PDF, and the site’s spin-chain conventions.

We test the explicit real-root branch

$$k_1=-k_2=k=\frac{2\pi m}{N-1},\qquad 1\leq m\lt\frac{N-1}{2}.$$

Here $S_{12}=e^{-ik}$, and the ring equations reduce to $e^{ik(N-1)}=1$. The energy prediction is

$$E=2J(1-\cos k).$$

The code evaluates these analytic roots; it does not run a root finder or enumerate all solutions.

Entry check and repair

For $v=(1,i)^{\mathsf T}$, compute $v^\dagger v$ and normalize the vector. Let $A=\operatorname{diag}(1,3)$. Its normalized expectation value is $2$; does that make $v$ an eigenvector with eigenvalue $2$?

Repair. The inner product uses complex conjugation: $v^\dagger v=1+(-i)i=2$. Thus $\widehat v=v/\sqrt2$. The residual is

$$(A-2\mathrm{Id})\widehat v=\frac1{\sqrt2}(-1,i)^{\mathsf T}, \qquad \left\|(A-2\mathrm{Id})\widehat v\right\|_2=1.$$

An expectation value is not an eigenvalue equation. In NumPy, use np.vdot(v, v) for the complex inner product and np.linalg.norm(H @ v - E * v) for the residual. Normalize first: scaling a wrong vector toward zero would otherwise make its unnormalized residual small.

Reproduce the calculation

Download and extract the complete finite-chain experiment (ZIP), then open its xxx folder. Individual files are also available: experiment.py, inputs.json, results.json, requirements.txt, and the README.

The recorded run used Python 3.9.6 and NumPy 2.0.2. For a fresh environment with Python 3.9–3.12, run from the extracted xxx folder:

python3 -m venv .venv
.venv/bin/python -m pip install -r requirements.txt
.venv/bin/python experiment.py --check

On Windows use .venv\Scripts\python.exe. With NumPy already installed, use python3 experiment.py --check in the extracted folder. In a repository checkout, run python3 public/computations/xxx/experiment.py --check from the repository root.

The check recomputes the experiment and compares it with the saved values without changing any files. With no option the script prints fresh JSON; --write deliberately replaces results.json. Make extensions in a copied folder so that the original experiment remains available for comparison.

Build the Hamiltonian in two ways

The full construction uses explicit Pauli tensor products in the $2^N$-dimensional space. The local basis is up $(1,0)^{\mathsf T}$ and down $(0,1)^{\mathsf T}$, with site zero the leftmost tensor factor.

The independent sector construction uses increasing tuples of down-spin positions. On one bond,

$$\begin{aligned} h\lvert\uparrow\uparrow\rangle&=0, &h\lvert\downarrow\downarrow\rangle&=0,\\ h\lvert\uparrow\downarrow\rangle &=\frac J2(\lvert\uparrow\downarrow\rangle-\lvert\downarrow\uparrow\rangle),\\ h\lvert\downarrow\uparrow\rangle &=\frac J2(\lvert\downarrow\uparrow\rangle-\lvert\uparrow\downarrow\rangle). \end{aligned}$$

These rules produce a diagonal contribution $J/2$ and a swapped-state contribution $-J/2$ whenever neighboring spins differ. The sector routine does not call the Pauli construction or share its bond iterator. Comparing the resulting $M=1,2$ matrices catches implementation errors that a single representation can conceal.

Before testing Bethe states, the script checks Hermiticity, total-$S^z$ conservation, the zero-energy vacuum, and the entire one-magnon spectrum

$$E_m=J\left(1-\cos\frac{2\pi m}{N}\right),\qquad m=0,\ldots,N-1.$$

It also applies the one-magnon matrix to every corresponding plane wave. These checks fix the energy scale, sign, and ring boundary condition before introducing the two-magnon scattering phase.

Work through six sites

Set $N=6$, $J=1$, and $k=2\pi/5$. Writing $z=e^{ik}$ and $d=y-x$ gives

$$\psi_{\mathrm{raw}}(x,y)=z^{-d}+z^{d-1} =2e^{-ik/2}\cos\left[k\left(d-\frac12\right)\right].$$

Its squared norm is $30$, so the normalized coefficients are $\psi=\psi_{\mathrm{raw}}/\sqrt{30}$. One way to verify that normalization is to group the 15 basis states by separation:

$$\sum_{x\lt y}\lvert\psi_{\mathrm{raw}}(x,y)\rvert^2 =\sum_{d=1}^{5}(6-d)\left\lvert z^{-d}+z^{d-1}\right\rvert^2=30.$$

The prediction $E/J=(5-\sqrt5)/2$ is checked against the Hamiltonian and its spectrum. The second selected six-site state has $k=4\pi/5$ and $E/J=(5+\sqrt5)/2$.

Define the dimensionless unit-vector residual and the independent energy discrepancy by

$$r=\frac{\lVert H\psi-E\psi\rVert_2}{J}, \qquad \delta_E=\min_{\ell}\frac{\lvert E-E_\ell^{\mathrm{diag}}\rvert}{J}.$$

The saved run gives the following values. For $r$, the table reports the larger residual from the full tensor and sector constructions.

$k_1=-k_2$	$E/J$	$r$	$\delta_E$
$2\pi/5$	$1.381966011250105$	$1.91\times10^{-16}$	$4.44\times10^{-16}$
$4\pi/5$	$3.618033988749895$	$7.34\times10^{-16}$	$1.33\times10^{-15}$

The script also checks both algebraic ring equations and the cyclic seam $\psi(x,y)=\psi(y,x+N)$ for every ordered pair. The saved file retains the complex normalized coefficients. When comparing to diagonalization, it projects onto the whole matching eigenspace, so arbitrary eigenvector phases or a different basis in a degenerate eigenspace do not create false failures.

What the wider checks and failures show

Across $N=4,5,6,7,8$, the script tests nine selected two-magnon states. The full and sector Hamiltonian blocks agree exactly in the recorded arithmetic. The largest one-magnon spectral discrepancy divided by $J$ is $8.88\times10^{-16}$, the largest two-magnon eigenstate residual $r$ is $8.23\times10^{-16}$, and the largest algebraic Bethe residual is $1.13\times10^{-15}$.

These values are consistent with binary64 rounding. The equation checks use tolerance $10^{-11}$. Saved-data comparison allows absolute tolerance $5\times10^{-12}$ and relative tolerance $10^{-10}$ for platform-dependent last digits. There is no time step or spatial discretization to refine: the matrices are the finite system being tested.

Two deliberate mistakes give much larger residuals. For the first six-site state:

Calculation	$\lVert H\psi-E\psi\rVert_2/J$
Correct periodic state	$1.91\times10^{-16}$
Replace $S_{12}$ by $S_{12}^{-1}$ only	$1.0454$
Delete the closing bond only	$0.47992$

The wrong-ratio control keeps the same roots and energy; its predicted energy still occurs in the spectrum. Thus a small nearest-energy discrepancy alone would miss this error. Reversing the ratio while also consistently relabeling the roots would be a different, equivalent convention change; this control deliberately changes only the ratio.

A formal solution can give a zero vector

At $N=5$, take exactly $z_1=z_2=-1$. The displayed scattering formula gives $S_{12}=-1$, and both algebraic equations hold: $(-1)^5=-1$. But

$$\psi_{\mathrm{raw}}(x,y)=(-1)^{x+y}-(-1)^{x+y}=0$$

for every basis state. The code rejects its zero norm instead of pretending to normalize it. It also rejects the undefined $0/0$ ratio at $z_1=z_2=1$. Other singular limits may require separate limiting constructions; these tests do not classify them.

There are also physical states outside the selected branch. The script independently constructs $(S^-_{\mathrm{tot}})^2\lvert\uparrow\cdots\uparrow\rangle$ in the tensor basis. After normalization it is the uniform two-down-spin state, and its energy is zero. Every bond annihilates this uniform vector because swapping neighboring spins leaves it unchanged. This is a descendant of the ferromagnetic vacuum, not one of the selected nonzero-momentum root pairs.

Two tested six-site states cannot span a 15-dimensional sector. Neither the selected real roots, the descendant check, nor small residuals establish completeness. Other total momenta, complex roots, and singular constructions remain outside this project’s enumeration.

Exercises

Guided practice: normalize before testing

For $N=6$ and $k=2\pi/5$, form the coefficients for all $0\leq x\lt y\leq5$. Verify the raw squared norm by the separation sum above, normalize, and calculate the Bethe residual and $r$. Then retain the same energy and roots but replace only $S_{12}$ by its inverse. Explain why the nearest diagonalized energy does not diagnose this mistake while $r$ does.

Hint · Full solution

Independent practice: another state and another size

Without copying the first state’s energy, repeat the reconstruction for $N=6$, $k=4\pi/5$. Then use $N=7$, $k=2\pi/6$ and predict its energy and raw squared norm before reading the saved results. Explain which of the following you have established: an eigenstate at each tested size, every state in those sectors, or a thermodynamic-limit result.

Hint · Full solution

Changed setting: cut the ring

For the first six-site state, remove only the bond between sites $5$ and $0$, keeping the same normalized wavefunction and predicted energy. Identify the new operator $H_{\mathrm{open}}$, relate its residual to the removed bond, and reproduce the negative-control value. Does that residual mean that the open Hamiltonian is non-Hermitian or that every open-chain model is nonintegrable? State what would need to change in a valid open-boundary Bethe calculation.

Hint · Full solution

Hints

Normalization and orientation. There are $6-d$ ordered pairs with separation $d$. Use absolute squares for complex coefficients. The wrong ratio changes the vector but not the energy formula evaluated at the original roots.

Other cases. Evaluate $2(1-\cos k)$ separately for each angle. For the selected regular branch, the finite separation sum gives raw squared norm $N(N-1)$.

Open boundary. Write $H_{\mathrm{open}}=H-h_{5,0}$. Start from the exact periodic equation $(H-E)\psi=0$ before considering floating-point error.

Solutions and checks

Normalization and orientation

The separation sum is $30$ and the normalized state uses $1/\sqrt{30}$. The energy is $E/J=(5-\sqrt5)/2$. The saved Bethe-equation residual is $2.78\times10^{-16}$, and the larger of the two Hamiltonian residuals is $1.91\times10^{-16}$.

With only the ratio inverted, the dimensionless Hamiltonian residual $r$ becomes about $1.0454$ and the algebraic Bethe residual about $1.9021$. The energy number is unchanged, so it is still within roundoff of an eigenvalue. The reconstructed vector, however, fails the contact and ring conditions in the retained convention. An energy match checks a scalar; an eigenstate residual checks every basis component.

Other cases

For the second six-site state,

$$\frac EJ=2\left(1-\cos\frac{4\pi}{5}\right) =\frac{5+\sqrt5}{2},\qquad \lVert\psi_{\mathrm{raw}}\rVert_2^2=30.$$

Its larger dimensionless Hamiltonian residual is $7.34\times10^{-16}$. For $N=7$ and $k=\pi/3$, $E/J=1$ and the raw squared norm is $7\cdot6=42$. The code verifies the corresponding state with the same independent constructions and residual threshold.

These are checks of particular finite-size eigenstates. The two-magnon sectors have dimensions $15$ and $21$, respectively, and the chosen branch supplies only a subset. Varying a few small chain lengths neither enumerates those sectors nor controls an $N\to\infty$ limit.

Open boundary

The missing bond gives

$$H_{\mathrm{open}}=H-\frac J2(1-P_{5,0}).$$

For the exact periodic eigenvector,

$$(H_{\mathrm{open}}-E)\psi=-\frac J2(1-P_{5,0})\psi.$$

Its dimensionless residual $\lVert(H_{\mathrm{open}}-E)\psi\rVert_2/J$ is about $0.47992465$, agreeing with the implemented missing-bond control. The changed Hamiltonian remains Hermitian and preserves the number of down spins. The failed test says that this periodic state and energy do not solve that changed boundary problem.

A valid open-boundary construction must satisfy the endpoint equations and account for reflected waves with the appropriate quantization conditions. The periodic seam equations cannot simply be retained after deleting a bond. No general nonintegrability conclusion follows from this mismatch.

Interpret the result

A successful project result includes the fixed convention, a nonzero normalized state, residual definitions, agreement between independent Hamiltonians, and controls that detect meaningful mistakes. Here those checks support the stated finite eigenpairs. They do not substitute for the analytic contact derivation, prove Bethe completeness, or establish the full commuting-transfer-matrix structure.

Return to the XXX sequence to connect the calculation with the model and derivation. For extensions, preserve the original inputs and state exactly which boundary, root family, or observable has changed.

The singular-state reproduction develops a different finite-chain test: a root pair can solve cleared Bethe equations while its unregularized vector vanishes. The Research benchmark separates a physical pair with a delicate limit from a formally similar pair rejected by translation symmetry.

References

• Karbach, Michael, and Gerhard Müller. “Introduction to the Bethe ansatz I.” Computers in Physics 11 (1997), pp. 36–43. DOI. Author version arXiv:cond-mat/9809162v1, submitted 1998. Open PDF, with printed pages 1–8.