Richardson–Gaudin Integrable Models

• Richardson–Gaudin integrable models are exactly solvable quantum many-body systems defined by a commuting family of conserved charges, enabling exact spectral descriptions via Bethe ansatz or eigenvalue equations.
• They utilize determinant formulas and eigenvalue-based methods to compute overlaps, norms, and form factors, providing numerically efficient frameworks for complex pairing, spin, and bosonic systems.
• These models have broad applications in superconductivity, quantum optics, and topological superconductors, and have been extended to include anisotropic, non-Hermitian, and PT-symmetric cases.

Richardson–Gaudin integrable models are exactly solvable quantum many-body systems defined by a commuting family of conserved charges and realized in pairing Hamiltonians, long-range spin models, spin-boson systems, and bosonic pairing problems. In their standard form, they are built from Gaudin algebras and admit exact spectral descriptions either through Bethe roots or through conserved-charge eigenvalues satisfying quadratic relations; in more recent developments, the class has been extended to fully anisotropic spin-$\tfrac12$ models in arbitrary magnetic fields, bosonic contraction limits such as the Dicke model, non-Hermitian and $\mathcal{PT}$-symmetric settings, and determinant-based computational frameworks for overlaps, norms, and form factors (Claeys et al., 2015, Dimo et al., 2018, Faribault et al., 2018).

1. Algebraic definition and canonical families

A standard Richardson–Gaudin construction starts from $n$ copies of $su(2)$ with generators $S_i^\dagger$, $S_i$, $S_i^0$, satisfying

$$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$

The conserved charges are

$$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$

and mutual commutativity $[R_i,R_j]=0$ follows when the coefficients satisfy the Gaudin equations

$\mathcal{PT}$0

(Claeys et al., 2015). In this framework, integrable Hamiltonians are linear combinations of the $\mathcal{PT}$1.

The standard algebraic classification is organized by the invariant

$\mathcal{PT}$2

The rational or XXX family corresponds to $\mathcal{PT}$3, the trigonometric family to $\mathcal{PT}$4, and the hyperbolic family to $\mathcal{PT}$5. Standard parameterizations quoted in the literature include

$\mathcal{PT}$6

for the rational model,

$\mathcal{PT}$7

for the trigonometric model, and

$\mathcal{PT}$8

for the hyperbolic model (Claeys et al., 2015).

A broader spin-$\mathcal{PT}$9 formulation replaces the conventional antisymmetric Gaudin couplings by the most general local-plus-bilinear ansatz

$n$0

For spin-$n$1, the local identity

$n$2

relaxes the usual Gaudin constraints and permits non-antisymmetric couplings, including a fully anisotropic XYZ class in arbitrary field (Claeys et al., 2018). In that setting, the conserved charges can be parametrized by

$n$3

together with site-dependent $n$4, $n$5, and $n$6 couplings dressed by the same affine functions of $n$7 (Claeys et al., 2018).

2. Bethe ansatz, conserved-charge eigenvalues, and eigenstate constructions

In the conventional RG description, eigenstates are Bethe states. For the general XXZ framework one writes

$n$8

with $n$9 the lowest-weight vacuum, while the rapidities obey

$su(2)$0

(Claeys et al., 2015). In the rational spin-$su(2)$1 case this becomes the familiar Richardson system

$su(2)$2

or, in the standard pairing parametrization,

$su(2)$3

(Claeys et al., 2017).

A central development in the modern theory is the shift from rapidity variables to conserved-charge eigenvalues. For XXZ models one introduces

$su(2)$4

which are directly related to the conserved-charge eigenvalues, and for spin-$su(2)$5 they satisfy the closed equations

$su(2)$6

supplemented by

$su(2)$7

(Claeys et al., 2015). In the rational limit,

$su(2)$8

For the most general spin-$su(2)$9 RG models in field, the commuting charges satisfy quadratic operator identities

$S_i^\dagger$0

which descend to quadratic equations for the eigenvalues $S_i^\dagger$1 on any common eigenstate (Dimo et al., 2018). This result covers XYZ, XXZ, and XXX models, including cases without $S_i^\dagger$2 symmetry and without a properly defined pseudo-vacuum (Dimo et al., 2018). A further step is the “Bethe-Ansatz-free” construction of eigenstates: for an on-shell solution $S_i^\dagger$3, one defines

$S_i^\dagger$4

and the normalized projector onto the corresponding eigenstate is

$S_i^\dagger$5

(Faribault et al., 2018). This removes the distinction between models with and without $S_i^\dagger$6 symmetry at the level of eigenstate construction (Faribault et al., 2018).

The same eigenvalue-based viewpoint is also useful when the explicit Bethe equations are unavailable. For the generalized integrable BCS model studied by Skrypnyk, the conserved operators $S_i^\dagger$7 satisfy quadratic relations, and the continuum limit is formulated directly in terms of their eigenvalues $S_i^\dagger$8, leading to a nonlinear singular integral equation whose solution yields the ground-state energy (Shen et al., 2019).

3. Realizations in pairing, spin, bosonic, and spin-boson systems

The canonical physical realization is the fermionic pairing problem. In the rational Richardson model,

$S_i^\dagger$9

with quasi-spin operators

$S_i$0

the spin Hamiltonian maps to the pairing Hamiltonian

$S_i$1

(Kulish et al., 2014). This reduced BCS model is described there as central in superconductivity and especially relevant for mesoscopic and nuclear systems where particle number is too small for mean-field BCS to be reliable (Kulish et al., 2014).

Bosonic contraction limits provide another major branch of the subject. By pseudo-deforming one $S_i$2 copy into a bosonic degree of freedom, the Dicke model is obtained as the contraction limit of an $S_i$3-based trigonometric RG system. In that limit,

$S_i$4

and the Dicke Bethe state becomes

$S_i$5

with Richardson–Gaudin equations

$S_i$6

(Claeys et al., 2014).

A closely related bosonic realization is the Lipkin–Meshkov–Glick model. Via the Schwinger boson map

$S_i$7

the LMG Hamiltonian becomes a two-level boson pairing problem and can be embedded into an $S_i$8 RG model. In that description, exact states are determined by pairons $S_i$9 satisfying reduced Richardson equations, and pairon trajectories provide a direct description of first-, second-, and third-order phase transitions and of spectral crossings (Lerma et al., 2012).

Richardson–Gaudin models also admit spin-boson realizations relevant to quantum optics. In the inhomogeneous Tavis–Cummings model,

$S_i^0$0

one introduces

$S_i^0$1

$S_i^0$2

and the commuting generating function $S_i^0$3 yields local conserved charges $S_i^0$4 and the conserved excitation number

$S_i^0$5

(Tschirhart et al., 2018). The diagonal ensemble for quenches then depends only on the shared eigenstates of $S_i^0$6, not on the detailed choice of coefficients $S_i^0$7 (Tschirhart et al., 2018).

The fermionic side has likewise been enlarged beyond the reduced BCS model. The generalized integrable BCS Hamiltonian introduced by Skrypnyk contains both number-conserving and non-number-conserving pairing terms, interpolates between the open and closed $S_i^0$8 models, and has commuting conserved operators $S_i^0$9 satisfying quadratic relations (Shen et al., 2019).

4. Determinant formulas and computational frameworks

A defining practical feature of modern RG theory is the existence of determinant representations for scalar products, partition functions, norms, and form factors. For rational models, the overlap between an on-shell state and an off-shell state can be written in two determinant forms, one of size $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$0 naturally expressed through the eigenvalue-based variables $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$1, and one of size $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$2 naturally expressed through rapidities; both arise from the domain-wall boundary partition function and the Cauchy-matrix structure underlying the model (Claeys et al., 2017). In the same framework, the usual Slavnov determinant is recovered as a reduction of the $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$3 formula (Claeys et al., 2017).

The general XXZ case also admits determinant expressions. Starting solely from the Gaudin algebra, one obtains overlap, normalization, and local-form-factor formulas directly in terms of the $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$4, showing that many rational-model constructions are in fact parametrization-independent consequences of the Gaudin algebra itself (Claeys et al., 2015). The paper emphasizes that these results generalize those for rational RG models and Dicke–Jaynes–Cummings–Gaudin models and expose a universality linked to the underlying algebraic structure (Claeys et al., 2015).

A mixed-spin extension was constructed for the rational model with one arbitrary spin $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$5 and $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$6 spins $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$7. There the domain-wall partition function, originally a permanent of a Cauchy-like $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$8 matrix with $$[S_i^0,S^{\dagger}_j]=\delta_{ij}S^{\dagger}_i, \qquad [S_i^0,S_j]=-\delta_{ij}S_i, \qquad [S^{\dagger}_i,S_j]=2\delta_{ij}S_i^0.$$9, is rewritten as an $$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$0 determinant in the variables

$$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$1

which are polynomial combinations of

$$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$2

(Faribault et al., 2015). This is explicitly designed to use the quadratic Bethe-equation variables rather than the rapidities themselves (Faribault et al., 2015).

The same program extends to integrable models with a bosonic mode. Starting from $$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$3 XXZ RG models and taking the pseudo-deformation contraction limit, determinant expressions for scalar products and form factors were extended to the Dicke–Jaynes–Cummings–Gaudin models and to the two-channel $$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$4-wave pairing Hamiltonian (Claeys et al., 2015).

These determinant and eigenvalue-based methods are not merely algebraic conveniences. They are repeatedly presented as numerically advantageous because the $$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$5-type variables satisfy quadratic systems that are much easier to solve than the original highly nonlinear Bethe equations, and because many observables can then be evaluated without reconstructing the full rapidity set (Faribault et al., 2015, Tschirhart et al., 2018).

5. Deformations, anisotropy, and non-Hermitian generalizations

The RG class has been enlarged in several distinct directions. One is quantum-group deformation. In the Jordanian deformation of the Richardson model, the rational $$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$6-invariant $$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$7-matrix

$$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$8

is replaced by

$$R_i=S_i^0+g \sum_{k \neq i}^n\left[\frac{1}{2}X_{ik}(S^{\dagger}_i S_k+S_iS^{\dagger}_k)+Z_{ik}S_i^0 S_k^0\right],$$9

and preserving integrability requires the insertion of a nilpotent auxiliary term

$[R_i,R_j]=0$0

The resulting transfer matrix still satisfies

$[R_i,R_j]=0$1

and exact eigenstates are created by the deformed operators

$[R_i,R_j]=0$2

(Kulish et al., 2014). The paper explicitly remarks that these exact eigenstates have a “highly complex entanglement structure” requiring further investigation (Kulish et al., 2014).

A second line of development is the fully anisotropic spin-$[R_i,R_j]=0$3 XYZ class in arbitrary magnetic field. Here the commuting charges have distinct $[R_i,R_j]=0$4, $[R_i,R_j]=0$5, and $[R_i,R_j]=0$6 couplings and need not be antisymmetric under $[R_i,R_j]=0$7. The central point is that spin-$[R_i,R_j]=0$8 permits this generalization because the local quadratic identity collapses the terms that, for higher spin, would enforce the usual antisymmetry (Claeys et al., 2018). This led directly to the generic quadratic operator relations later used to formulate eigenvalue-based Bethe equations for arbitrary-field spin-$[R_i,R_j]=0$9 RG models (Dimo et al., 2018).

Open and non-Hermitian extensions form another major branch. In one exactly solvable Lindblad problem with collective dissipation, the Liouvillian maps to a non-Hermitian XXZ RG model acting on spin-1 triplet degrees of freedom,

$\mathcal{PT}$00

with Bethe states, rapidity equations, pseudo-Hermiticity

$\mathcal{PT}$01

exceptional points, dissipative quantum phase transitions, a nontrivial steady state in the homogeneous limit, and logarithmic gap growth away from homogeneity (Claeys et al., 2021). A related high-temperature noisy-spin problem maps the Liouvillian blocks of correlation tensors to a non-Hermitian rational spin-1 RG model with complex inhomogeneities $\mathcal{PT}$02, giving exact spectral equations for decay modes and long-lived correlations (Rowlands et al., 2017).

A more recent development is the explicitly $\mathcal{PT}$03-symmetric spin-$\mathcal{PT}$04 RG model in arbitrary magnetic field. There one defines parity and time-reversal transformations, constructs the metric operator

$\mathcal{PT}$05

uses it to derive Hermitian counterparts of the $\mathcal{PT}$06-symmetric conserved charges, and finds spectra containing both real eigenvalues and complex-conjugate pairs (AlMasri, 1 Jun 2025). The same work reports that at weak coupling the system fails to reach a steady state, whereas at stronger coupling it eventually does so (AlMasri, 1 Jun 2025).

6. Applications, thermodynamic limits, and variational use

Richardson–Gaudin models are not used only as exactly solvable Hamiltonians. They are also used as computational frameworks for nearby non-integrable problems. In nuclear pairing, the Richardson–Gaudin Configuration Interaction method first variationally optimizes an RG ground state and then uses the complete set of excited RG states as an optimized CI basis for a realistic non-integrable pairing Hamiltonian. In the benchmark for the Sn region, the variational RG step already reaches accuracies around the 1% level of the correlation energies, and the RGCI truncation exhibits an additional improvement scaling exponentially with the size of the effective Hilbert space (Baerdemacker et al., 2017). The paper’s central claim is that RG integrability supplies an optimized complete basis set for pairing correlations (Baerdemacker et al., 2017).

A closely related variational method uses on-shell RG eigenstates to approximate the ground state of integrability-breaking spin models. Because the trial manifold consists of exact RG states, Slavnov determinants, Gaudin norms, and eigenvalue-based variables remain available for efficient energy minimization. The method is exact in the integrable limit, improves substantially on perturbation theory for models close to integrability, and shows that for large integrability-breaking perturbations the relevant variational state may need to be an excited RG state rather than the integrable ground state (Claeys et al., 2017).

The thermodynamic limit can also be formulated directly in conserved-charge language. For the generalized integrable BCS model extending the open and closed $\mathcal{PT}$07 systems, the continuum limit of the quadratic eigenvalue equations yields an integral equation whose ground-state solution produces

$\mathcal{PT}$08

in exact agreement with the BCS mean-field result (Shen et al., 2019).

Integrable RG models have also been used to construct number-conserving topological superconducting chains. A new two-parameter rational solution of the Yang–Baxter/RG functional equation generates an integrable $\mathcal{PT}$09-$\mathcal{PT}$10 wave Richardson–Gaudin–Kitaev chain with pairing form factor

$\mathcal{PT}$11

critical coupling

$\mathcal{PT}$12

and a topological phase transition for which the occupancy of the non-interacting mode serves as a topological order parameter (Stouten et al., 2017).

In quantum optics, exact finite-size diagonal-ensemble calculations for the inhomogeneous spin-boson XXX RG family reveal a strong-coupling tendency toward a common steady state in which

$\mathcal{PT}$13

interpreted in the paper as a “superradiant-like” coherent steady state (Tschirhart et al., 2018). This does not establish a thermodynamic phase transition, but it demonstrates how the RG conserved structure can control non-equilibrium observables in spin-boson models (Tschirhart et al., 2018).

Taken together, these developments show that Richardson–Gaudin integrable models function both as exactly solvable many-body systems and as adaptable algebraic infrastructures. They provide commuting conserved quantities, exact or projector-based eigenstate constructions, quadratic eigenvalue equations, determinant expressions for observables, controlled contraction and deformation limits, and optimized bases for non-integrable many-body calculations across pairing physics, bosonic models, spin chains, open quantum systems, and topological superconductivity (Baerdemacker et al., 2017, Claeys et al., 2017).