# Five-site XXX matrix elements and a finite correlation sum

This small, deterministic benchmark connects a regular algebraic Bethe vector to normalized matrix elements and a complete finite spectral sum. It concerns the periodic spin-1/2 XXX chain with

```text
H = (J/2) sum_x (I - P_(x,x+1)),  J > 0,  N = 5,
hbar = 1, lattice spacing = 1, site 6 = site 1.
```

There is no stochastic sampling, time-stepping solver or frequency broadening. Binary64 roundoff is tested against exact finite identities. This calculation is not a thermodynamic limit, a general Bethe-completeness result, a response commutator or a time-ordered correlation.

## Run the complete download

The [complete ZIP](https://integrable.org/computations/xxx-correlations.zip) preserves the required layout:

```text
integrable-xxx-correlations/
  xxx-correlations/
    experiment.py
    inputs.json
    results.json
    requirements.txt
    README.md
  xxx-algebra/
    experiment.py
```

Open the `xxx-correlations` folder, then run with Python 3.9–3.12:

```sh
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`. The recorded environment is Python 3.9.6 and NumPy 2.0.2. In the repository with NumPy installed:

```sh
python3 public/computations/xxx-correlations/experiment.py --check
```

The companion `xxx-algebra/experiment.py` supplies **only the actual upper-right monodromy-block Bethe-vector construction**. Importing it does not run its larger experiment. No companion inputs/results or checkout-specific absolute paths are needed. If downloading individual files, obtain the companion separately and preserve the sibling folders; do not rename it over this experiment. The ZIP is the simplest way to retain that layout. Routine site builds package the files without running either experiment.

`--check` recomputes the benchmark without writing any file, locks the complete inputs exactly, and compares required result fields, integer/Boolean values, finite floating-point values, and array lengths. Float tolerances are absolute `2e-12` and relative `1e-9`; runtime environment strings are reported but ignored in comparison. NaNs and infinities are rejected. Running without flags prints fresh JSON. `--write` deliberately replaces `results.json` after all scientific checks pass; use a separate copy for explorations.

## State and basis conventions

Single-site tensor order is `(up, down)`, site 1 is the leftmost factor, and `|0>` is all up. The companion uses `R(u)=uI+iP`, `L_x(lambda)=R_(a,x)(lambda-i/2)`, `T=L_5...L_1` and `B=T_12`. The state is

```text
Phi = B(1/2) B(-1/2) |0>.
```

In lexicographic occupied-pair order `(12,13,14,15,23,24,25,34,35,45)`, its coefficients are `-1/8` for cyclic adjacent pairs (including 15) and `+1/8` otherwise. Its exact norm squared is `5/32`. Normalizing gives coefficients `±1/sqrt(10)`, energy `2J` and active translation eigenvalue `1`. The script compares the actual B-block product with these explicit coefficients, verifies no amplitude lies outside the two-down-spin sector, and independently checks the physical bond Hamiltonian and translation. Fraction arithmetic verifies the dyadic raw coefficients' norm and bond action exactly; it does not turn the floating-point matrix multiplication into a symbolic calculation.

The final one-down-spin basis is `|x>`, `x=1,...,5`. The physical bond implementation uses zero for parallel spins and `(J/2)(|up,down>-|down,up>)` on an antiparallel input. It does not call the companion's Hamiltonian, swap or translation routines. A separately constructed `S_x^+` map removes a down spin with coefficient **one**. There is no extra spin-1/2 prefactor in a raising-operator matrix element.

## Fourier matrix elements and normalization

Use all five normalized states

```text
|k> = (1/sqrt(5)) sum_x exp(i k x) |x>,   k=2*pi*m/5,
U|x> = |x+1>,                         U|k>=exp(-ik)|k>,
O_q^+ = (1/sqrt(5)) sum_x exp(-i q x) S_x^+.
```

Since the initial normalized state `chi` has coordinate momentum zero, `O_q^+ chi` has coordinate momentum `-q` modulo `2*pi` and active translation eigenvalue `exp(+iq)`. Its only possible matrix element in this sector is

```text
M_q = <-q|O_q^+|chi> = 2[cos(2q)-cos(q)]/sqrt(10).
```

The script checks the complete complex amplitude, its orthogonal remainder and the active translation equation. The zero-momentum vector vanishes, so that case uses an absolute residual rather than division by its norm. The five squared weights are exactly `0,1/2,1/2,1/2,1/2`, summing to the two initial down spins. Independent rational polynomial reduction modulo `1+z+z^2+z^3+z^4` verifies these weights and the five-point Fourier completeness relation. This is completeness of a five-dimensional Fourier basis, not arbitrary Bethe-state completeness.

For a local site `y`, the amplitudes in that same basis have squared weights `0,1/10,1/10,1/10,1/10`, summing to `||S_y^+ chi||^2=2/5`. The code checks the local norm at all five sites and the individual Fourier amplitudes at the recorded site. Individual local amplitudes depend on `y` and phase conventions; weights do not.

## Correlation and discrete spectrum

The observable is the **unsymmetrized**, possibly complex correlation

```text
C_y(t) = <chi| S_y^-(t) S_y^+(0) |chi>
       = sum_k |<k|S_y^+|chi>|^2 exp[-i(E_k-2J)t],
E_k    = J(1-cos k).
```

The recorded local site is `y=1`. Inserting the complete Fourier basis gives

```text
C_y(t) = (2/5) exp(3i Jt/4) cos(sqrt(5) Jt/4).
```

Independently, the code diagonalizes the **physical bond-built 5x5 Hamiltonian**, applies its unitary exponential to `S_y^+ chi`, and contracts the resulting vector. The analytic Fourier basis and its known energies are not inputs to that calculation. All three answers are compared at 121 equally spaced dimensionless times `Jt` from 0 through 12. This is a finite matrix identity evaluated in binary64; it has no time-step convergence order or Monte Carlo confidence interval.

For `S_y(omega)=integral dt exp(i omega t) C_y(t)`, the spectrum is a sum of delta distributions. The table gives integrated weights under `domega/(2*pi)`, not heights of a broadened curve:

| One-magnon modes | Frequency divided by J | Integrated weight |
|---|---:|---:|
| 1 and 4 | `-(3+sqrt(5))/4 = -1.3090169943749475` | `1/5` |
| 2 and 3 | `-(3-sqrt(5))/4 = -0.19098300562505255` | `1/5` |

The mode at `k=0` has frequency `-2J` but zero weight. Both actual lines have negative frequency because `chi` is an **excited** state and the operator connects it to lower-energy states. This is not a ground-state zero-temperature structure factor. Degenerate partners are distinct orthonormal intermediate states, so their weights add. Numerical eigensolvers can rotate a degenerate eigenspace; the independently diagonalized check therefore compares its **total projector weight**, not an arbitrarily chosen eigenvector's individual weight.

Equal-time and first-moment identities supply two different checks:

```text
C_y(0) = sum_k w_k = 2/5,
integral domega/(2*pi) omega S_y(omega)
       = sum_k (E_k-2J) w_k
       = <chi|S_y^- [H,S_y^+]|chi> = -3J/10.
```

The first moment is computed both from the Fourier weights and from the physical operator `(H_1-2J) S_y^+ chi`, and checked exactly with Fractions on the raw state. Correspondingly, `C_y'(0)=+3iJ/10`.

## Recorded results and failure controls

The recorded maximum among the independent accepted numerical errors is `1.013e-15`. The maximum direct-evolution discrepancy from the closed formula is `1.013e-15`, and from the Fourier spectral sum is `6.510e-16`. Small differences in final digits across platforms are expected; the comparison tolerances are stated above. These numbers are roundoff-scale residuals, not error bars for an approximate physical model.

| Deliberate error | Computed consequence |
|---|---:|
| Use the raw B state without normalizing | Local equal-time value `1/16 = 0.0625` instead of `0.4` |
| Omit one intermediate mode with weight `0.1` | Equal-time value `0.3`; time-correlation maximum error `0.10000000000000016` |
| Keep one vector per distinct nonzero-weight energy, losing one partner at each energy | Equal-time value `0.2` instead of `0.4` |
| Omit the complete line at `omega/J=-(3+sqrt(5))/4` | Equal-time value `0.2` instead of `0.4` |
| Reverse the Fourier sign for mode 1 | Wrong-translation relative residual `1.902113032590307` |
| Reverse the Fourier sign for mode 2 | Wrong-translation relative residual `1.1755705045849463` |

The last error preserves the weights and energies in this symmetric example. Translation detects it: the wrong state has eigenvalue `exp(-iq)` instead of `exp(+iq)`, giving residual `2|sin(q)|`. Conversely, multiplying the initial raw ket by `3+4i` and then normalizing leaves all local weights unchanged. Rephasing the final Fourier kets by `exp(0.7i)` changes the complex amplitudes by `exp(-0.7i)` but leaves their squared moduli unchanged. These controls distinguish a physical normalization/direction error from a harmless phase convention.

`equation_errors` records absolute scalar/vector errors or Frobenius matrix errors as its names specify. Initial B-coefficient comparison is relative to the explicit nonzero raw vector. H-eigenvector errors use `H/J` and a unit state; selection-rule errors are absolute vector norms, including `q=0`. The wrong-direction control uses `||U v-exp(iq)v||_2/||v||_2`. Correlation comparisons use maximum absolute complex error over the stated time grid. Frequencies and moments in the JSON are divided by `J`. All participating norms are Euclidean unless a Frobenius matrix norm is named.

## Figure provenance

The original figure is `public/figures/xxx/correlation-spectrum.svg`; its editable source is `figures-src/xxx/correlation-spectrum.py`. It plots the closed finite correlation and the saved integrated line weights. The source verifies their signs, sum, first moment, and agreement with saved independent matrix evolution. It introduces no line broadening or interpolation-based frequency estimate. Regenerate from the repository with NumPy and Matplotlib (recorded version 3.9.4):

```sh
python3 figures-src/xxx/correlation-spectrum.py
python3 figures-src/xxx/correlation-spectrum.py --check
```

The second command regenerates the SVG in memory and compares it without writing a figure. Matplotlib is needed only to regenerate the figure, not to run the downloadable experiment.

## Sources and what is adapted

- Jean-Sébastien Caux, “Correlation functions of integrable models: a description of the ABACUS algorithm,” *Journal of Mathematical Physics* **50**, 095214 (2009), DOI [10.1063/1.3216474](https://doi.org/10.1063/1.3216474). [arXiv:0908.1660v1 PDF](https://arxiv.org/pdf/0908.1660v1), § II, printed pp. 3–4, equations (1)–(3). This supplies the correlation/Fourier/Lehmann framework. Its expectation is explicitly a ground-state expectation; the excited-state finite example and its negative frequencies here are derived independently. Our Fourier operator includes `1/sqrt(N)`, while the source places the corresponding `1/N` outside its squared matrix element.
- L. D. Faddeev, “How Algebraic Bethe Ansatz works for integrable model,” in *Quantum Symmetries / Symétries quantiques*, Les Houches Session LXIV, edited by A. Connes, K. Gawędzki and J. Zinn-Justin, North-Holland (1998), pp. 149–219. [hep-th/9605187v1 PDF](https://arxiv.org/pdf/hep-th/9605187v1), § 4, equations (82)–(97), supports the vacuum/creation-operator construction reused from the companion. The actual five-site vector is also checked directly against an independent finite spin Hamiltonian.

Neither source is cited as a proof that this five-site calculation computes generic thermodynamic correlation functions or implements the ABACUS algorithm.