---
title: Richardson–Gaudin Integrable Models
url: https://www.emergentmind.com/topics/richardson-gaudin-integrable-models
type: topic
---

# Richardson–Gaudin Integrable Models

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 [1501.05827] [1805.03427] [1812.06428].

## 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
\[
X_{ij}=-X_{ji}, \qquad Z_{ij}=-Z_{ji}, \qquad X_{ij}X_{jk}-X_{ik}(Z_{ij}+Z_{jk})=0
\]
[1501.05827]. In this framework, integrable Hamiltonians are linear combinations of the \(R_i\).

The standard algebraic classification is organized by the invariant
\[
(X_{ij})^2-(Z_{ij})^2=\Gamma.
\]
The rational or XXX family corresponds to \(\Gamma=0\), the trigonometric family to \(\Gamma>0\), and the hyperbolic family to \(\Gamma<0\). Standard parameterizations quoted in the literature include
\[
X_{ij}=Z_{ij}=\frac{1}{\epsilon_i-\epsilon_j}
\]
for the rational model,
\[
X_{ij}=\frac{1}{\sin(\epsilon_i-\epsilon_j)}, \qquad Z_{ij}=\cot(\epsilon_i-\epsilon_j)
\]
for the trigonometric model, and
\[
X_{ij}=\frac{1}{\sinh(\epsilon_i-\epsilon_j)}, \qquad Z_{ij}=\coth(\epsilon_i-\epsilon_j)
\]
for the hyperbolic model [1501.05827].

A broader spin-\(\tfrac12\) formulation replaces the conventional antisymmetric Gaudin couplings by the most general local-plus-bilinear ansatz
\[
Q_i = B_i^x S_i^x + B_i^y S_i^y + B_i^z S_i^z + \sum_{j\neq i}^L \left[ \Gamma_{ij}^x S_i^x S_j^x + \Gamma_{ij}^y S_i^y S_j^y + \Gamma_{ij}^z S_i^z S_j^z \right].
\]
For spin-\(\tfrac12\), the local identity
\[
S_i^\alpha S_i^\beta = \frac{i}{2}\sum_\gamma \epsilon_{\alpha\beta\gamma} S_i^\gamma +\frac14 \delta_{\alpha\beta}\,\mathbbm 1
\]
relaxes the usual Gaudin constraints and permits non-antisymmetric couplings, including a fully anisotropic XYZ class in arbitrary field [1810.06059]. In that setting, the conserved charges can be parametrized by
\[
B_i^x = \frac{\gamma}{\sqrt{\alpha_x \epsilon_i+\beta_x}}, \qquad
B_i^y = \frac{\lambda}{\sqrt{\alpha_y \epsilon_i+\beta_y}}, \qquad
B_i^z = 1,
\]
together with site-dependent \(xx\), \(yy\), and \(zz\) couplings dressed by the same affine functions of \(\epsilon_i\) [1810.06059].

## 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
\[
\ket{\psi_N}=\prod_{\alpha=1}^N\left(\sum_{i=1}^n X_{i\alpha}S^{\dagger}_i\right)\ket{\theta},
\]
with \(\ket{\theta}\) the lowest-weight vacuum, while the rapidities obey
\[
1+g\sum_{i=1}^n Z_{i\alpha}d_i-g\sum_{\beta \neq \alpha}^N Z_{\beta\alpha}=0
\]
[1501.05827]. In the rational spin-\(\tfrac12\) case this becomes the familiar Richardson system
\[
1+\frac{g}{2}\sum_{i=1}^nZ_{i\alpha} = g \sum_{\beta \neq \alpha}^NZ_{\beta\alpha},
\]
or, in the standard pairing parametrization,
\[
1+\frac{g}{2}\sum_{i=1}^L \frac{1}{\epsilon_i-\lambda_{\alpha}}-g \sum_{\beta \neq \alpha}^N\frac{1}{\lambda_{\beta}-\lambda_{\alpha}}=0
\]
[1707.06793].

A central development in the modern theory is the shift from rapidity variables to conserved-charge eigenvalues. For XXZ models one introduces
\[
\Lambda_i \equiv \sum_{\alpha=1}^N Z_{i\alpha},
\]
which are directly related to the conserved-charge eigenvalues, and for spin-\(\tfrac12\) they satisfy the closed equations
\[
\Lambda_i^2=N(n-N)\Gamma-\frac{2}{g}\Lambda_i+\sum_{j \neq i}^n Z_{ji}(\Lambda_j-\Lambda_i),
\]
supplemented by
\[
-\frac{g}{2}\sum_{i=1}^n \Lambda_i=N
\]
[1501.05827]. In the rational limit,
\[
\Lambda_i=\sum_{\alpha=1}^N\frac{1}{\epsilon_i-x_\alpha}, \qquad
\Lambda_i^2=-\frac{2}{g}\Lambda_i+\sum_{j \neq i}^n \frac{\Lambda_j-\Lambda_i}{\epsilon_j-\epsilon_i}.
\]

For the most general spin-\(\tfrac12\) RG models in field, the commuting charges satisfy quadratic operator identities
\[
R_i^2 = \sum_{j\neq i} C_{ij} R_j + K_i,
\]
which descend to quadratic equations for the eigenvalues \(r_i\) on any common eigenstate [1805.03427]. This result covers XYZ, XXZ, and XXX models, including cases without \(U(1)\) symmetry and without a properly defined pseudo-vacuum [1805.03427]. A further step is the “Bethe-Ansatz-free” construction of eigenstates: for an on-shell solution \((r_1,\dots,r_N)\), one defines
\[
\hat P(r_1,\dots,r_N) \equiv \det \hat J, \qquad
\hat J_{ii}=r_i+R_i, \qquad \hat J_{ij}=-\Gamma_{ij}\quad (i\neq j),
\]
and the normalized projector onto the corresponding eigenstate is
\[
|\psi_n\rangle\langle \psi_n| =
\frac{ \hat P\bigl(r_1^{(n)},\dots,r_N^{(n)}\bigr) }{ N\bigl(r_1^{(n)},\dots,r_N^{(n)}\bigr) }
\]
[1812.06428]. This removes the distinction between models with and without \(U(1)\) symmetry at the level of eigenstate construction [1812.06428].

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 \(Q_i\) satisfy quadratic relations, and the continuum limit is formulated directly in terms of their eigenvalues \(q_i\), leading to a nonlinear singular integral equation whose solution yields the ground-state energy [1912.05692].

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

The canonical physical realization is the fermionic pairing problem. In the rational Richardson model,
\[
H=\sum_{j=1}^N \epsilon_j S_j^z + g \sum_{j,k=1}^N S_j^- S_k^+,
\]
with quasi-spin operators
\[
S_l^z = \frac12 \sum_m c^\dagger_{lm} c_{lm} - \frac12,\qquad
S_l^+ = \frac12 \sum_m c^\dagger_{lm} c^\dagger_{l\bar m},\qquad
S_l^-=(S_l^+)^\dagger,
\]
the spin Hamiltonian maps to the pairing Hamiltonian
\[
H_P=\sum_l \epsilon_l n_l + \frac g2 \sum_{l,l'} A_l^\dagger A_{l'}
\]
[1402.4838]. 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 [1402.4838].

Bosonic contraction limits provide another major branch of the subject. By pseudo-deforming one \(su(2)\) copy into a bosonic degree of freedom, the Dicke model is obtained as the contraction limit of an \(su(2)\)-based trigonometric RG system. In that limit,
\[
H_{\mathrm{Dicke}}=\hbar\omega\, b^\dagger b+\sum_{k=1}^m\varepsilon_k S_k^0 +G\sum_{k=1}^m \left(b^\dagger S_k + S_k^\dagger b\right),
\]
and the Dicke Bethe state becomes
\[
|\psi\rangle\propto \prod_{\alpha=1}^N \left( b^\dagger -G\sum_{k=1}^m \frac{S_k^\dagger}{\varepsilon_k-x_\alpha} \right) |\theta\rangle,
\]
with Richardson–Gaudin equations
\[
(\hbar\omega-x_\alpha)-2G^2\sum_{k=1}^m\frac{s_k}{\varepsilon_k-x_\alpha} +2G^2\sum_{\beta\neq\alpha}^N\frac{1}{x_\beta-x_\alpha}=0
\]
[1411.5199].

A closely related bosonic realization is the Lipkin–Meshkov–Glick model. Via the Schwinger boson map
\[
S_z=\frac{ b^\dagger b-a^\dagger a}{2}, \qquad S_+=b^\dagger a, \qquad S_-= a^\dagger b,
\]
the LMG Hamiltonian becomes a two-level boson pairing problem and can be embedded into an \(SU(1,1)\otimes SU(1,1)\) RG model. In that description, exact states are determined by pairons \(e_\alpha\) satisfying reduced Richardson equations, and pairon trajectories provide a direct description of first-, second-, and third-order phase transitions and of spectral crossings [1212.3238].

Richardson–Gaudin models also admit spin-boson realizations relevant to quantum optics. In the inhomogeneous Tavis–Cummings model,
\[
H_{TC}=\omega b^\dagger b+\sum_{j=1}^N \epsilon_j S_j^z + V\sum_{j=1}^N \left(b^\dagger S_j^-+S_j^+ b\right),
\]
one introduces
\[
\mathrm{S}^+(u)= b^\dagger + \sum_{j=1}^N \frac{V}{u-\epsilon_j}S_j^+,\qquad
\mathrm{S}^-(u)= b + \sum_{j=1}^N \frac{V}{u-\epsilon_j}S_j^-,
\]
\[
\mathrm{S}^z(u)=\frac{\omega-u}{2V}-\sum_{j=1}^N \frac{V}{u-\epsilon_j}S_j^z,
\]
and the commuting generating function \(S^2(u)\) yields local conserved charges \(R_i\) and the conserved excitation number
\[
\hat M = b^\dagger b+\sum_{i=1}^N \left(S_i^z+\frac12\right)
\]
[1805.03479]. The diagonal ensemble for quenches then depends only on the shared eigenstates of \(H=\gamma \hat M+\sum_i \alpha_i R_i\), not on the detailed choice of coefficients \(\gamma,\alpha_i\) [1805.03479].

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 \(p+ip\) models, and has commuting conserved operators \(Q_i\) satisfying quadratic relations [1912.05692].

## 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 \(L\times L\) naturally expressed through the eigenvalue-based variables \(\Lambda_i\), and one of size \(2N\times 2N\) naturally expressed through rapidities; both arise from the domain-wall boundary partition function and the Cauchy-matrix structure underlying the model [1706.05511]. In the same framework, the usual Slavnov determinant is recovered as a reduction of the \(2N\times 2N\) formula [1706.05511].

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 \(\Lambda_i\), showing that many rational-model constructions are in fact parametrization-independent consequences of the Gaudin algebra itself [1501.05827]. 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 [1501.05827].

A mixed-spin extension was constructed for the rational model with one arbitrary spin \(S\) and \(N-1\) spins \(\tfrac12\). There the domain-wall partition function, originally a permanent of a Cauchy-like \(\Omega\times\Omega\) matrix with \(\Omega=2S+(N-1)\), is rewritten as an \(N\times N\) determinant in the variables
\[
\Gamma_1(\epsilon_1),\dots,\Gamma_{2S}(\epsilon_1),\Gamma_1(\epsilon_2),\dots,\Gamma_1(\epsilon_N),
\]
which are polynomial combinations of
\[
\Lambda(\epsilon_1),\Lambda'(\epsilon_1),\dots,\Lambda^{(2S)}(\epsilon_1),\Lambda(\epsilon_2),\dots,\Lambda(\epsilon_N)
\]
[1511.03127]. This is explicitly designed to use the quadratic Bethe-equation variables rather than the rapidities themselves [1511.03127].

The same program extends to integrable models with a bosonic mode. Starting from \(su(2)\) 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 \((p+ip)\)-wave pairing Hamiltonian [1506.03702].

These determinant and eigenvalue-based methods are not merely algebraic conveniences. They are repeatedly presented as numerically advantageous because the \(\Lambda\)-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 [1511.03127] [1805.03479].

## 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 \(sl(2)\)-invariant \(r\)-matrix
\[
r(\lambda,\mu)=\frac{C_2^\otimes}{\lambda-\mu}
\]
is replaced by
\[
r^{(J)}(\lambda,\mu)=\frac{C_2^\otimes}{\lambda-\mu}+\xi\,(h\otimes X^+ - X^+\otimes h),
\]
and preserving integrability requires the insertion of a nilpotent auxiliary term
\[
L(\lambda;\xi,c)= c\, X_0^+ + L(\lambda,\xi).
\]
The resulting transfer matrix still satisfies
\[
[t^{(J)}(\lambda),t^{(J)}(\mu)]=0,
\]
and exact eigenstates are created by the deformed operators
\[
B_M(\mu_1,\dots,\mu_M) = X^-(\mu_1)\bigl(X^-(\mu_2)+\xi\bigr)\cdots \bigl(X^-(\mu_M)+\xi(M-1)\bigr)
\]
[1402.4838]. The paper explicitly remarks that these exact eigenstates have a “highly complex entanglement structure” requiring further investigation [1402.4838].

A second line of development is the fully anisotropic spin-\(\tfrac12\) XYZ class in arbitrary magnetic field. Here the commuting charges have distinct \(xx\), \(yy\), and \(zz\) couplings and need not be antisymmetric under \(i\leftrightarrow j\). The central point is that spin-\(\tfrac12\) permits this generalization because the local quadratic identity collapses the terms that, for higher spin, would enforce the usual antisymmetry [1810.06059]. This led directly to the generic quadratic operator relations later used to formulate eigenvalue-based Bethe equations for arbitrary-field spin-\(\tfrac12\) RG models [1805.03427].

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,
\[
\cL = i\sum_{j=1}^L\left[\Omega+\omega_j\right]\s^z_j  - \gamma_{0}\sum_{j,k=1}^L \s^z_j \s^z_k
 - \gamma \sum_{j,k=1}^L \sqrt{\omega_j\omega_k}\left(\s^x_j \s^x_k+\s^y_j \s^y_k\right),
\]
with Bethe states, rapidity equations, pseudo-Hermiticity
\[
\cL^\dagger = P^{\dagger}\cL P,
\]
exceptional points, dissipative quantum phase transitions, a nontrivial steady state in the homogeneous limit, and logarithmic gap growth away from homogeneity [2108.01677]. 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 \(\epsilon_j=i\omega_j/2\), giving exact spectral equations for decay modes and long-lived correlations [1711.00828].

A more recent development is the explicitly \(\mathcal{PT}\)-symmetric spin-\(\tfrac12\) RG model in arbitrary magnetic field. There one defines parity and time-reversal transformations, constructs the metric operator
\[
\rho = \mathcal{P}\mathcal{C} = \eta^\dagger \eta,
\]
uses it to derive Hermitian counterparts of the \(\mathcal{PT}\)-symmetric conserved charges, and finds spectra containing both real eigenvalues and complex-conjugate pairs [2506.01110]. The same work reports that at weak coupling the system fails to reach a steady state, whereas at stronger coupling it eventually does so [2506.01110].

## 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 [1712.01673]. The paper’s central claim is that RG integrability supplies an optimized complete basis set for pairing correlations [1712.01673].

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 [1707.06793].

The thermodynamic limit can also be formulated directly in conserved-charge language. For the generalized integrable BCS model extending the open and closed \(p+ip\) systems, the continuum limit of the quadratic eigenvalue equations yields an integral equation whose ground-state solution produces
\[
\frac{E}{L} = \frac{\mathcal G \hat\Delta_x^2}{4} +\frac{\mathcal G \hat\Delta_y^2}{4} -\frac12 \int_{\omega_0}^{\omega} d\varepsilon\, \rho(\varepsilon)\, R(\varepsilon),
\]
in exact agreement with the BCS mean-field result [1912.05692].

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 \(s\)-\(d\) wave Richardson–Gaudin–Kitaev chain with pairing form factor
\[
\eta_k= 4 \sin^2(k/2)(t_1 + 4 t_2\cos^2(k/2)),
\]
critical coupling
\[
G_c^{-1}=-\sum_k \eta_k,
\]
and a topological phase transition for which the occupancy of the non-interacting mode serves as a topological order parameter [1712.09375].

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
\[
\overline{\langle S_i^z\rangle}\approx 0 \qquad \forall i,
\]
interpreted in the paper as a “superradiant-like” coherent steady state [1805.03479]. 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 [1805.03479].

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 [1712.01673] [1707.06793].

Source: https://www.emergentmind.com/topics/richardson-gaudin-integrable-models