Skip to content

Complex roots for the quantum-group-invariant open XXZ chain

XXZ family · Supported models · Biquadratic application

This is a separate solver for the same opposite-end-field reference Hamiltonian used by the TL layer, not the ordinary zero-field XXZ front end. Even N, M<=N/2 and finite Delta>1 are supported. It keeps finite complex-root deviations; no ideal-string hypothesis is imposed. All arithmetic stays in the selected fp64, long-double, binary128 or explicit MPFR precision. Explicit MPFR is available in the C++ interface and in the biquadratic frontend's --q-seed and --q-spectrum modes. Use --precision bits-N or digits-N; the same even-chain and bounded-search restrictions apply.

Why solve a polynomial?

Real-root logarithmic equations use ordered roots and integer phases. They cannot simply be continued through arbitrary complex-string configurations. Instead we encode the roots in a monic polynomial

\[ \begin{aligned} x&=\cosh(2u)=\cos\alpha,\\{} Q(x)&=\prod_j(x-x_j)=x^M+\sum_{k=0}^{M-1}c[k]x^k. \end{aligned} \]

The coefficients are real; zeros can be real or form conjugate pairs. Reflection u -> -u and periodicity u -> u+i*pi are absorbed by x. This does not assert support for arbitrary non-self-conjugate root configurations.

Bajnok et al., section 5 give a Q-system for precisely these boundaries. We use its Wronskian relation (5.17). With eta=acosh(Delta), s=sinh(eta), and x_±=Delta*x ± s*sqrt(x^2-1), our normalization is

\[ \frac{P(x_+)Q(x_-)-P(x_-)Q(x_+)}{2s\sqrt{x^2-1}}=(x-1)^N. \]

P has degree N-M+1; its leading coefficient is fixed by this identity, and we fix the P -> P+aQ freedom with P_M=0. Both sides are polynomials in x. The implementation evaluates them without square roots: the antisymmetric basis uses \(x _+\,x _-=x ^{2}+s ^{2}\) and the recurrence

\[ D_0=0,\quad D_1=1,\qquad D_{k+1}=2\Delta xD_k-(x^2+s^2)D_{k-1}. \]

For fixed Q, descending coefficient elimination determines P using only constant, nonzero q-number pivots. The remaining M coefficients (powers 0,...,M-2 and 2M-1) form the nonlinear system. The Jacobian differentiates that elimination analytically. A damped Newton solve uses the raw residual for descent, preventing normalization alone from rewarding diverging coefficients. The stopping residual is the full Wronskian's componentwise coefficient backward error, not the old logarithmic phase residual or an energy error estimate.

The reference energy is obtained directly from Q:

\[ E_{\mathrm{ref}}=\frac{(N-1)\Delta}{4}-(\Delta^2-1)\frac{Q'(\Delta)}{Q(\Delta)}. \]

The TL and physical biquadratic energy maps and module multiplicities are unchanged. Root extraction is for independent validation and display, not the energy formula.

Numerical acceptance and limitations

A small coefficient residual is not by itself an accepted physical state. The solver also:

  • recovers separated roots at the selected precision using the polynomial-root engine;
  • rejects roots numerically indistinguishable from x=1, -1 or Delta (the excluded u=0, i*pi/2 and eta/2 solutions);
  • checks the original multiplicative Bethe equations, including reflected scattering, using scaled products to avoid product overflow;
  • propagates the root engine's estimated uncertainty through the sinh factors. Near a string pole, this can be much larger than machine epsilon. The per-equation residual must fit its own estimated bound, and the propagated relative uncertainty must be below 1e-4. Otherwise the state remains unverified, even if its polynomial residual is tiny.

The reported bethe_residual and bethe_residual_bound are maxima over equations; acceptance is checked equation by equation. These estimates are first-order floating-point diagnostics, not interval certificates. An unresolved root or exceptional configuration is never rescued by silently switching precision or relaxing the user's coefficient tolerance. The caller can explicitly rerun in higher precision.

There is no site cutoff. The CLI warns beyond the N=8 spectrum-validation range: the cached polynomial basis takes O(N^3) storage, and conditioning can become limiting well before memory does. Attempt and iteration budgets bound the search, not its wall time or memory footprint. Allocation-size overflow is rejected before constructing the workspace. This method is not a scalable long-chain excitation solver; use the separate two-string singlet solver for that targeted branch. Continuation can reuse a converged Q at a nearby Delta; automated adaptive continuation remains future work.

State discovery versus solving a state

spectrum() samples deterministic coefficient seeds and polynomials built from real roots and conjugate pairs. It uses neither ED seeds nor stored eigenvalues. All attempts, including failures and duplicate solutions, consume max_attempts. A zero budget returns no discoveries. Different precisions can follow different basins.

Solutions are deduplicated by their Q coefficients, not by energy, and sorted only after discovery. The target dimension is choose(N,M)-choose(N,M-1). The count is optional: integer overflow leaves it unavailable, never wrapped or approximate, and disables a completeness claim. A zero attempt budget returns before allocating the polynomial basis. Matching the count is a numerical completeness check; close solutions may be conservatively merged and leave a search incomplete. An incomplete set does not establish the lowest k energies, even if sorted.

C++ interface

#include <bethe/biquadratic_qsystem.hpp>

std::vector<long double> seed{1.5L, -4.0L/3.0L};
auto level = bethe::biquadratic::qsystem::solve<long double>(4, seed);
if (level.reference.converged) {
  auto energy = *level.energy;
  auto copies = *level.multiplicity; // 1 for this ell=0 singlet
}
auto scan = bethe::biquadratic::qsystem::spectrum<long double>(6, 2);
// scan.complete(), expected_count, attempts, failed_attempts, states

The underlying bethe::xxz::quantum_group::qsystem::{System,solve,spectrum} API takes Delta; bethe::temperley_lieb::qsystem::solve takes loop weight. Only the biquadratic adapter fixes Delta=3/2 and attaches spin-1 multiplicities. Selected failures preserve their coefficients and available energy estimate; missing/pole energies are optional, not plausible-looking zeroes.

The library also accepts explicit MPFR precision, with MPC complex root recovery and Uni20's normal LU dispatch for Newton corrections. CPU MPFR works without MPLAPACK; an enabled provider is preferred by the dispatch:

using R = uni20::mpreal;
auto p = uni20::Precision::bits(160);
bethe::SolverOptions<R> controls(p);
std::vector<R> seed{R(3, p) / 2, -R(4, p) / 3};
auto level = bethe::biquadratic::qsystem::solve<R>(4, seed, controls);

Delta, seeds and solver controls must share that finite precision; exact, unset or mixed-precision inputs are rejected. Coefficients, recovered roots, energies and failure diagnostics retain the requested precision. The bounded multistart search keeps its deterministic seed sequence and attempt budget; more arithmetic bits do not establish completeness or replace root admissibility checks. Independent solves can use owning Async query grids, while Newton iteration, continuation and the multistart search stay sequential. The biquadratic frontend exposes this precision through --precision bits-N and --precision digits-N, including its other supported complex-state modes.

The regression suite checks the analytic Jacobian, the complex Wronskian identity, the original Bethe equations, branch continuation, native-precision four-site singlets, failure budgets, small TL-module ED spectra and the full physical N=2,4,6 Hilbert-space spectra with representation weights. Eight-site module spectra are checked against independent auxiliary XXZ ED: with the tested extended precision and sufficient attempts all module counts are recovered, while fp64 can leave nearly ideal strings unverified. MPFR checks at 80/160/256 bits compare the four-site complex branch with an independent TL eigenvalue beyond binary128 accuracy. They also check the shifted hyperbolic Wronskian, analytic one-root bands, bounded searches and complete mixed-precision success/failure results under Async execution.