Finite-temperature Lieb–Liniger equilibrium¶

The library implements the repulsive Yang–Yang equation in the grand-canonical ensemble: give c>0, T>0, and chemical potential mu; obtain the density rather than prescribing it. A fixed-density wrapper instead accepts c,T,n and solves for mu. bethe-lieb-liniger-thermal exposes both ensembles, including temperature scans. Chemical-potential scans are not yet implemented.

The Hamiltonian is H=-sum d_i²+2c sum delta(x_i-x_j), with \(\hbar =2m =k_{B} =1\). Pressure and energy are per length, density is particles per length, and entropy is also per length. The grand potential per length is minus the pressure. The uniform thermodynamic limit is not a finite-ring spectrum or a trapped gas.

Command-line examples¶

# Grand-canonical equilibrium, with density determined by mu.
bethe-lieb-liniger-thermal --c 4 --temperature 1 --mu -1

# Fixed-density temperature scan; mu is solved independently at each point.
bethe-lieb-liniger-thermal --c 4 --temperature 0.5 \
  --temperature-end 2 --points 9 --density 1 --max-nodes 512 --threads 4 --csv thermal.csv

# Extended precision and simultaneous structured exports.
bethe-lieb-liniger-thermal --c 4 --temperature 1 --density 0.1 \
  --precision fp128 --json thermal.json --tsv thermal.tsv

# Explicit 256-bit arithmetic; quadrature still has its own accuracy budget.
bethe-lieb-liniger-thermal --c 4 --temperature 1 --mu -1 \
  --precision bits-256 --tolerance 1e-55 --max-nodes 512

bethe-lieb-liniger-thermal --references

Choose exactly one of --mu and --density. Without --temperature-end, the tool evaluates one point. With it, --points defaults to 17 samples, includes both endpoints, and permits ascending or descending positive temperatures. The specified c and ensemble parameter stay fixed throughout the scan.

--threads N uses Uni20's scheduler for independent temperatures (default 1). Each point owns its equilibrium solve and, for fixed density, its chemical-potential search; those internal iterations remain sequential. The shared bounded observable grid collects rows on the calling thread in the requested temperature order, including failures. It admits at most 128 active or completed-but-uncollected temperatures. No continuation seed is passed between temperatures, so concurrent execution preserves the serial calculation. Set OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 to avoid nested threading. Single-temperature requests have no grid parallelism. Explicit MPFR scans use the same scheduler; working precision and convergence controls stay local to each request.

The named table thermodynamics reports temperature, chemical potential, density, pressure, energy/length, entropy/length, quadrature nodes and cutoff, error estimates, iteration/evaluation counts, and status. Mesh/cutoff diagnostics are from the final inner solve, not totals; Newton updates and equilibrium evaluations are totals for that row. The density mismatch applies only to fixed-density runs. After a failed canonical solve, inner diagnostics and mu are unavailable; after a failed grand-canonical solve, the supplied mu remains known. Neither publishes failed physical observables.

The numerical controls mirror the C++ options below: --tolerance, --initial-nodes, --max-nodes, --initial-cutoff, --max-cutoffs, and --max-iterations. Fixed-density runs additionally accept --density-tolerance and --max-evaluations. Defaults depend on the selected working precision.

The shared output options provide presentation-layer screen output, CSV/TSV/JSON files, provenance metadata and CPU timing. Rows stream to attached exports during the scan; --no-retain avoids keeping a second copy in memory. CPU timing covers the complete temperature grid, including ordered row output; the exported CPU timing scope states this explicitly. Failure at a point does not stop subsequent temperatures: its physical values are null in JSON or empty in CSV/TSV, and the final exit code is 2. Invalid input returns 1, before opening output files. References appear only with --references.

Equations¶

With C(x)=c/(pi*(c²+x²)), the Yang–Yang equation is

\[ \begin{aligned} \epsilon(k)&=k^2-\mu-\int C(k-q)T\log\!\left(1+e^{-\epsilon(q)/T}\right)\,dq,\\{} \mathrm{filling}(k)&=\frac{1}{1+e^{\epsilon(k)/T}},\\{} \rho_{\mathrm{total}}(k)&=\frac1{2\pi} +\int C(k-q)\mathrm{filling}(q)\rho_{\mathrm{total}}(q)\,dq. \end{aligned} \]

All integrals cover the real axis. The filling has a Fermi-like form because it counts occupied Bethe quantum numbers; the physical particles remain bosons. Pressure integrates \(T \log (1+\exp (-\epsilon /T))/(2\,\pi)\). Particle and energy densities integrate filling*rho_total and k²*filling*rho_total respectively. Entropy integrates the binary entropy of the filling times rho_total.

Sources: Yang–Yang (1969), Guan–Batchelor (2011), equations (2)–(5), and Franchini's notes, section 2.11. We solve the nonlinear equation numerically; no finite-c polylog approximation is substituted. The ideal-Fermi fugacity series supplies an independent test in the Tonks limit.

C++ interface¶

#include <bethe/lieb_liniger_thermal.hpp>

namespace yy = bethe::lieb_liniger::thermal;
auto state = yy::equilibrium(4.0, 1.0, -1.0); // c, T, mu
if (state.converged) {
    auto pressure = *state.pressure;
    auto density = *state.density;
    auto energy_per_length = *state.energy_per_length;
    auto entropy_per_length = *state.entropy_per_length;
}

Use long double, configured binary128, or uni20::mpreal for extended arithmetic. MPFR equilibrium solves use normal Uni20 dispatch, including CPU LU. Runtime precision is explicit in both the inputs and controls:

auto p = uni20::Precision::bits(128);
using R = uni20::mpreal;
yy::Options<R> options(p);
options.tolerance = R("1e-20", p);
auto state = yy::equilibrium(R(4, p), R(1, p), R(-1, p), options);
yy::DensityOptions<R> canonical(p);
canonical.equilibrium.tolerance = R("1e-18", p);
canonical.tolerance = R("1e-16", p);
auto fixed = yy::at_density(R(4, p), R(1, p), R("0.1", p), canonical);

These explicit tolerances keep the introductory examples inexpensive. Increase arithmetic precision and numerical accuracy independently; a high-accuracy fixed-density request performs several complete equilibrium solves.

Inputs, tolerances and an optional cutoff must have the controls' finite precision. Mixed, exact-only and unset runtime values are rejected before solving. Constants, quadrature, compensated reductions, Newton corrections, thermal factors and failure diagnostics retain that precision. The solver changes neither ambient MPFR precision nor its exponent range. There is no narrowing to binary64 inside the solver or its quadrature rules.

The CLI accepts --precision bits-N and digits-N when Uni20 MPFR scalars are configured. More arithmetic bits do not establish quadrature or tail convergence: insufficient mesh/cutoff budgets still return a failure. Very small observables are rejected if their requested absolute accuracy cannot be resolved on the scalar's finite grid; MPFR uses its current minimum positive value rather than a native floating-point underflow threshold. The default Options<Real> has relative tolerance 16384*epsilon, 16 initial nodes, at most 256 nodes, 512 accepted Newton updates across the whole solve, and at most four rapidity cutoffs. Node budgets may be increased to 512; cutoff budgets to 32. Zero iteration/cutoff budgets are allowed for diagnostics.

initial_cutoff, if specified, is a positive physical rapidity. Otherwise the initial domain is chosen from T, mu, and the requested tolerance. Successive cutoffs grow by 5/4. Mesh sizes double at each fixed cutoff, starting again at initial_nodes; a non-doubling max_nodes need not itself be sampled.

Three distinct convergence checks¶

A damped Newton solve converges the discretized pseudoenergy equation. Its Jacobian is also the dressed linear operator for rho_total.
Successive Gauss–Legendre meshes must agree for all four thermodynamic observables at a fixed cutoff.
Two separately mesh-converged cutoffs must agree for all four observables.

nodes counts the positive-half nodes: evenness supplies the negative half through \(C (k -q)+C (k +q)\). cutoff is physical; cutoffs counts attempted domains. mesh_error and cutoff_error are eight times the largest relative observable change. They are estimates, not rigorous bounds on omitted tails or quadrature. The library reports no physical observable until both checks pass. Small c, low T, or a large user-selected initial cutoff may need more nodes than allowed.

The mesh comparison is conservative: its difference can be dominated by the coarser rule even when the finer result is much more accurate. In a 256-bit Tonks-limit check at T=1, mu=-1 and c=1e120, tolerances from 1e-55 to 1e-65 converge at 512 nodes and agree with the independent fugacity-series values to about 1e-76 relative error. A 1e-70 request instead reports mesh_limit. Arithmetic accuracy in this example is established by that independent reference; the convergence diagnostic alone does not certify those extra digits.

Calculations scale energies by \(S =\max (T,\lvert \mu \rvert)\) and rapidities by sqrt(S). nonlinear_residual is the largest scaled-equation residual divided by max(1,abs(epsilon),abs(k²-mu)) in these dimensionless variables. Its target is tolerance*(T/S)/64, so the equation's precision requirement tracks temperature.

Stable shared thermal factors avoid exponentially overflowing expressions. Entropy is calculated directly, not by subtracting \(p +e -\mu \,n\) at low T. Precision-generic Gauss–Legendre rules, compensated sums, rational scattering kernels, and recoverable pivoting solves are shared with existing continuum/Hubbard solvers.

iteration_limit, mesh_limit, cutoff_limit, and precision_limit distinguish failure paths; all four optional observables remain empty on failure. A failed coarse mesh is retried at higher order; if the last mesh cannot resolve the arithmetic, the result is precision_limit. Invalid parameters throw. No asymptotic zero-density or T=0 answer is silently substituted on underflow. Use the separate zero-temperature API for T=0.

Validation includes the original integral equations, scale covariance, \(p +e -\mu \,n =T \,s\), pressure derivatives at fixed c, a native-precision fugacity series in the Tonks limit, decreasing-temperature comparison with the ground state, and independently exhausted nonlinear, mesh, and cutoff budgets.

Fixed density¶

auto canonical = yy::at_density(4.0, 1.0, 1.0); // c, T, n
if (canonical.converged) {
    auto mu = canonical.state->chemical_potential;
    auto free_energy_per_length = mu * *canonical.state->density
                               - *canonical.state->pressure;
}

DensityOptions<Real> contains the inner equilibrium options, an outer relative density tolerance (default \(65536\,\epsilon\)), and max_evaluations (default 128). Each evaluation is a complete grand-canonical solve with its own inner budgets; the inner tolerance is tightened to at most one eighth of the density tolerance. evaluations counts those solves and iterations sums their accepted Newton updates.

The chemical potential starts at the classical-gas estimate, but every density used is calculated from the full Yang–Yang equation. Exponentially expanding steps bracket the target; Illinois regula falsi then interpolates inside that bracket, falling back to its midpoint when rounding places a trial at an end. The stopping test requires a relative density mismatch at most half the outer tolerance, leaving room for the inner numerical error estimate. Small bracket width alone is not accepted as convergence.

density_error is the last successfully evaluated relative density mismatch. density_limit means the outer evaluation budget was exhausted. Inner failures retain their original status; no failed density is treated as zero and no partially converged thermodynamics are published. The optional state is set only on success, and its density is the calculated value, not overwritten by the requested density. Extreme parameters can still encounter the inner mesh or precision limits described above. Tests cover native-precision round trips, scale covariance, dilute and degenerate gases, and independent outer/inner budget exhaustion.