text:
# Electrostatic analogy for integrable pairing force Hamiltonians

## I Introduction

Pairing force interactions have been successfully employed to explain phenomena in different contexts such as superconductivity TINKHAM , nuclear physics IACHELLO94 , QCD RISCHKE  and astrophysicsASTRO . The original idea traces back to Bardeen, Cooper, and Schrieffer (BCS) who proposed pairing of time-reversed electrons as the crucial mechanism for superconductivityBCS . Physical implications of the corresponding Hamiltonian have been extracted resorting to a variety of analytical or numerical techniques. The mean field ground state was shown to be exact in the limit of large number of electrons where fluctuations can be neglected TINKHAM ; BCS ; BOGOLIUBOV . This constituted a great success of the BCS variational ansatz. However there are relevant physical situations where approximations are not reliable and exact treatments are highly desirable. Fortunately, by properly choosing the pairing couplings, the model admits an exact solutions.
A simplified, but still non-trivial model is the “reduced BCS model” (BCS model, for brevity) which assumes a uniform pairing $`g`$ among all the electrons within the Debye shell. This model is integrable CAMBIAGGIO  and was diagonalized long ago by Richardson through Bethe Ansatz (BA) RICHARDSON  (see also Ref. GAUDIN, ). Only very recently this generated a lot of interest both in nuclear and condensed matter physics whose communities benefited from the simple algorithm of Richardson’s BA solution to tackle the BCS model in the canonical ensemble BRT ; MASTELLONE ; VONDELFT ; NUCLEAR . Much work has been done to merge the model in the schemes of the Quantum Inverse Scattering (QIS) and Conformal Field Theory (CFT), which are modern arenas where quantum integrability and exact solutions can be treated on equal footing. These studies allowed significant steps forward. QIS studies identified the Richardson BA solution as the quasi classical limit of the exact solution of (twisted) disordered six vertex models of the XXX-type SKLYANIN-GAUDIN ; AMICO . This paved the way towards the exact evaluation of correlation functions of the BCS model AO ; AUSTRALIANS . The field theoretical study of the BCS model is due to Sierra and coworkers. The first step was to relate the Richardson BA solution with WZNW-$`su(2)_k`$ models SIERRA-CFT ; SIERRA-REV . This stimulated the discovery that the BA solution is the quasi classical limit of the Babujian’s off shell BA of the (untwisted) disordered six-vertex model OSBA ; AMICO . The field theoretical origin of the integrability of the BCS model has been clarified definitely SIERRA-CS  to be a twisted Chern-Simons (CS) theory on a torus. The Richardson BA solution together with the underlying integrability of the theory arises from the Knizhnik-Zamolodchikov-Bernard equations (which are the Knizhnik-Zamolodchikov equations on the torus). The emergence of the CS theory is quite interesting and relates at a formal level the BCS theory of superconductivity with the Fractional Quantum Hall Effect (FQHE) SIERRA-CFT  (also pointing out important differences). In fact, the exact BCS wave function admits a “Coulomb gas” representation, which corresponds to the “Coulomb gas” representation of the Laughlin wave functions (plasma analogy) PLASMA ; SIERRA-CFT . Accordingly, also for the BCS model an electrostatic analogy does exist. It was illustrated first by Gaudin GAUDIN-LARGE : the Cooper pair energies are obtained as equilibrium positions of $`N`$ mobile charges in the background of $`\mathrm{\Omega }`$ fixed charges and a uniform electric field of strength $`1/g`$. This analogy was used by Richardson RICH77  and Gaudin GAUDIN-LARGE  to obtain the thermodynamic limit of the Richardson BA equations in two different approaches. Both reproduced the BCS mean-field gap equation, hence confirming the statement of Bogoliubov BOGOLIUBOV  that at $`T=0`$ the mean-field results become exact. Recently, Sierra and coworkers revived the attention on the thermodynamic limit of the BCS model, comparing Gaudin’s long forgotten results with the numerical solution of the Richardson BA equations for a finite number of particles SIERRA-LARGE .

By generalizing the integrals of motion of the BCS model, the class of known integrable pairing-force Hamiltonians could be enlarged considerably towards models with non-uniform interactions ADO ; RICHARDSON-NEW  (for bosonic versions see the Ref.DUKELSKY, ). Here, the term “non-uniform” indicates that the interactions depend on the energy levels occupied by the interacting electrons. As the uniform BCS model, also these models emerge from the quasi-classical limit of the QIS approach for twisted disordered six vertex models of the XXZ-typeAMICO ; POGHOSSIAN . The untwisted case, studied in the Ref. HIKAMI , corresponds to large effective pairing interaction.

In the present work we generalize the $`2d`$-electrostatic analogy to the class of integrable Hamiltonians obtained in Ref. ADO, . The thermodynamic limit of these Hamiltonians is obtained. We compare the obtained analytical results with numerical solutions of the BA equations. The electrostatic analogy of the Richardson-Sherman equations is a screening condition for a total electric field produced by charges distributed in two-dimensional space. Accordingly, the electric field obeys a Riccati-type differential equations.

The paper is laid out as follows. In the next section we review the integrability of pairing force Hamiltonians that arises from the underlying infinite dimensional Gaudin algebra $`𝒢[sl(2)]`$. In section III we present the electrostatic analogy, review basic facts of $`2d`$-dimensional electrostatics and apply them to obtain the thermodynamic limit of integrable Hamiltonians with non-uniform pairing couplings. At the end of the section, we obtain Gaudin and Richardson’s results for the BCS model as a limiting case of our equations. Section IV is devoted to conclusions. In Appendix A we sketch some mathematical aspects connected with the integrability of the models. In Appendix B we present the connection with the Riccati equations by reviewing the work of Richardson RICH77 . In Appendix C we collect details of the calculations. In Appendix D we discuss some features of the ground state and possible routes towards the study of the excitations.

## II Exactly solvable pairing models

The Hamiltonian for $`N`$ charged particles interacting through pairing force

$$H=\underset{i\sigma }{}\epsilon _{i\sigma }n_{i\sigma }+\underset{ij}{}U_{ij}n_in_j\underset{ij}{}g_{ij}c_i^{}c_i^{}c_jc_j,$$
(1)
has been proved ADO  to be exactly solvable if the Coulomb $`U_{ij}`$ and the pairing strength $`g_{ij}`$ are restricted to have the following form

$$\begin{array}{ccc}\hfill g_{ij}& =& p\stackrel{~}{g}\frac{\epsilon _i\epsilon _j}{\mathrm{sinh}[p(u_iu_j)]},\hfill \\ & & \\ \hfill 4U_{ij}& =& Ap\stackrel{~}{g}(\epsilon _i\epsilon _j)\mathrm{coth}[p(u_iu_j)],\hfill \end{array}\}\text{for }ij$$
(2)
where $`\epsilon _i`$ are the single-particle levels, while for convenience 111In formulas (2) instead of the single particle energies the quantities $`\eta _j=\epsilon _jg_{jj}/2+2_iU_{ij}`$ enter. We fixed $`g_{jj}`$ and $`U_{jj}`$ such that $`\eta _j=\epsilon _j`$ we fix $`g_{ii}=\varphi `$, $`U_{ii}=A\varphi `$, where $`2\varphi =p\stackrel{~}{g}_{ij}(\epsilon _i\epsilon _j)\mathrm{coth}[p(u_iu_j)]`$. The quantities $`A`$, $`u`$, and $`u_i`$ are arbitrary real parameters while $`p`$ can be real or pure imaginary. To study the integrability of the Hamiltonian (1) it is convenient to reformulate it as a spin chain model. In fact it can be expressed (up to a constant term) as

$`H`$ $`={\displaystyle \underset{i}{}}2\epsilon _iK_i^z+4{\displaystyle \underset{i,j=1}{\overset{\mathrm{\Omega }}{}}}U_{ij}K_i^zK_j^z`$ (3)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j=1}{\overset{\mathrm{\Omega }}{}}}g_{ij}(K_i^+K_j^{}+K_i^{}K_j^+),`$
where $`K_j^+=c_j^{}c_j^{},K_j^{}=\left(K_j^+\right)^{},K_j^z=\frac{1}{2}(n_j+n_j1),`$ are $`su(2)`$ operators. The Hamiltonian (1) acts non trivially only on the Hilbert space $``$ of doubly occupied and empty level pairs; the singly occupied levels are excluded from the dynamics associated to (1). This phenomenon is commonly called the “blocking of singly occupied levels”. It reduces the problem of finding the spectrum of $`H`$ in the restricted Hilbert space $``$ and it is an important ingredient for the integrability of the Hamiltonian (1). Ultimately the latter resides in the algebraic connection of the pairing Hamiltonians with the infinite dimensional Gaudin algebra (see Appendix A) $`𝒢[sl(2)]=span\{K^\pm (u),K^z(u):u\text{ }\mathrm{C}\}`$: the diagonalization of the Hamiltonian is equivalent to the diagonalization of the correspondent of Casimir operator in $`𝒢[sl(2)]`$, $`K(u)`$ . This is achieved either through algebraic or coordinate-wise BA, and corresponds to the quasi-classical limit of the inhomogeneous $`XXZ`$ model. The set of operators in involution are the residues of $`K(u)`$ and the Hamiltonian is a polynomial of these operators which are, by construction, the integrals of motion of the model.
In the following we will see these statements at work for the model (1), (2). Then we will recover the BCS model as a limit of the pairing couplings.

#### II.0.1 Non uniform couplings.

We review the integrability of the pairing models (1). The key observation is that the integrals of motion $`\stackrel{~}{\tau }_i`$ of the model are residues of $`\stackrel{~}{\tau }(u)`$, defined in Appendix A, in $`u=u_i`$. In fact, $`H`$ can be written as

$`H`$ $`=`$ $`{\displaystyle \underset{i}{\overset{\mathrm{\Omega }}{}}}2\epsilon _i\stackrel{~}{\tau }_i+A{\displaystyle \underset{i,j}{\overset{\mathrm{\Omega }}{}}}\stackrel{~}{\tau }_i\stackrel{~}{\tau }_j{\displaystyle \underset{i}{}}\varphi _i𝐊_i^2,`$ (4)
$`\stackrel{~}{\tau }_i=K_i^z+\stackrel{~}{\mathrm{\Xi }}_i,`$ (5)
$`\stackrel{~}{\mathrm{\Xi }}_i`$ $`:=`$ $`p\stackrel{~}{g}{\displaystyle \underset{ji}{}}{\displaystyle \frac{K_i^+K_j^{}+K_i^{}K_j^+}{2\mathrm{sinh}[p(u_iu_j)]}}`$
$`+K_i^zK_j^z\mathrm{coth}[p(u_iu_j)]`$
The $`\stackrel{~}{g}\mathrm{}`$–limit of this Hamiltonian was found in Ref. HIKAMI, . From Eq.(4) it is evident that $`[H,\stackrel{~}{\tau }_j]=0`$. By these integrals of motion, the model is connected with the anisotropic Gaudin Hamiltonians $`\stackrel{~}{\mathrm{\Xi }}_l`$. The property ADO  $`[\stackrel{~}{\tau }_j,\stackrel{~}{\tau }_l]=0`$ for $`j,l=1,\mathrm{},\mathrm{\Omega }`$ comes from $`[\tau (u),\tau (v)]=0`$ SKLYANIN-GAUDIN . The common eigenstates of $`H`$ and $`\stackrel{~}{\tau }_i`$ are

$`|\mathrm{\Psi }`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{N}{}}}K^+(e_\alpha )|0`$ (7)
where $`|0`$ is the Fock vacuum and $`N`$ is the number of pairs. The eigenvalues of the Hamiltonian $`H`$ and of the constants of the motion are

$$=p\stackrel{~}{g}\underset{j=1}{\overset{\mathrm{\Omega }}{}}\underset{\alpha =1}{\overset{N}{}}\epsilon _j\mathrm{coth}[p(e_\alpha u_j)]+AN^2$$
(8)
$$\tau _i=\frac{\stackrel{~}{g}}{2}\underset{\alpha =1}{\overset{N}{}}p\mathrm{coth}[p(e_\alpha u_i)]+\frac{\stackrel{~}{g}}{4}\underset{\genfrac{}{}{0.0pt}{}{j=1}{ji}}{\overset{\mathrm{\Omega }}{}}p\mathrm{coth}[p(u_ju_i)]\frac{1}{2}.$$
(9)
The quantities $`e_\alpha `$ are solutions of the following set of equations

$$\frac{2}{\stackrel{~}{g}}+2\underset{\genfrac{}{}{0pt}{}{\beta =1}{\beta \alpha }}{\overset{N}{}}p\mathrm{coth}[p(e_\beta e_\alpha )]\underset{j=1}{\overset{\mathrm{\Omega }}{}}p\mathrm{coth}[p(u_je_\alpha )]=0.$$
(10)
which are the BA equations of the model (1).

#### II.0.2  Uniform couplings: the BCS model

We now review the exact solution together with the integrability of the pairing Hamiltonian with uniform coupling constants. In this case the Hamiltonian(1) reads

$`H_{BCS}={\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}2\epsilon _iK_i^z{\displaystyle \frac{g}{2}}{\displaystyle \underset{i,j=1}{\overset{\mathrm{\Omega }}{}}}\left(K_i^+K_j^{}+K_i^{}K_j^+\right)`$ (11)
$`H`$ can be directly diagonalized through coordinate-wise RICHARDSON  or algebraic BRAUN  BA. Nevertheless its diagonalization can be achieved along the lines depicted above. The integrals of motion $`\tau _i`$ of the BCS model are the residue of $`\tau (u)`$ in $`u=2\epsilon _i`$. By these integrals of motion, the model becomes connected with isotropic Gaudin Hamiltonians $`\mathrm{\Xi }_l`$. In fact, $`H`$ can be written as

$`H_{BCS}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}2\epsilon _i\tau _i+g{\displaystyle \underset{i,j=1}{\overset{\mathrm{\Omega }}{}}}\tau _i\tau _j+const.`$ (12)
$`\tau _i`$ $`=`$ $`K_i^z+\mathrm{\Xi }_i,`$ (13)
$`\mathrm{\Xi }_i`$ $`:=`$ $`g{\displaystyle \underset{ji}{}}{\displaystyle \frac{\stackrel{}{K}_i\stackrel{}{K}_j}{\epsilon _i\epsilon _j}}`$ (14)
where the spin vectors are: $`𝑲_j:=(K_j^x,K_j^y,K_j^z)`$; $`K_j^\pm =K_j^x\pm iK_j^y`$. The eigenstates of both the BCS Hamiltonian RICHARDSON ; GAUDIN  and its constants of motion $`\tau _j`$SKLYANIN-GAUDIN  are

$$|\mathrm{\Psi }_{BCS}=\underset{\alpha =1}{\overset{N}{}}K^+(e_\alpha )|0$$
(15)
The eigenvalues of the Hamiltonian and of the integrals of the motion are respectively

$`_{BCS}`$ $`=`$ $`2{\displaystyle \underset{\alpha =1}{\overset{N}{}}}e_\alpha ,`$ (16)
$`\tau _i`$ $`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle \underset{\alpha =1}{\overset{N}{}}}{\displaystyle \frac{1}{e_\alpha \epsilon _i}}+{\displaystyle \frac{g}{4}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{j=1}{ji}}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1}{\epsilon _j\epsilon _i}}{\displaystyle \frac{1}{2}}.`$ (17)
where $`e_\alpha `$ are solutions of the Richardson-Sherman (RS) equations

$$\frac{1}{g}+\underset{j=1}{\overset{\mathrm{\Omega }}{}}\frac{1/2}{e_\alpha \epsilon _j}\underset{\genfrac{}{}{0.0pt}{}{\beta =1}{\beta \alpha }}{\overset{N}{}}\frac{1}{e_\alpha e_\beta }=0$$
(18)
(the parameters $`e_\alpha `$ are one half of the parameters $`E_\alpha `$ originally defined by Richardson RICHARDSON ). The eigenstates, eigenvalues and the RS equations can be obtained in the $`p0`$ limit of the BA equations (10) with the identification $`u_i=\epsilon _i`$, $`\stackrel{~}{g}=g`$. By the same limit the anisotropic (XXZ) Gaudin model reduces to its isotropic version (XXX).

## III Electrostatic analogy for exactly solvable pairing models

The RS equations (18) are stating that the total force acting on the unit charges located in $`e_\alpha `$, due to the $`2d`$ electric field222We recall that in two dimensions the electric field generated by a point charge $`q`$ having complex coordinate $`z_0`$ can be represented by the holomorphic function $`E=q/(zz_0)`$. generated by the $`\mathrm{\Omega }`$ charges $`q_{fixed}=1/2`$ fixed in $`\epsilon _i`$, by the remaining $`N1`$ mobile unit charges $`q_{mobile}=+1`$ and by the external constant electric field of strength $`1/g`$ is zero. In other words, the solutions of the RS equations correspond to the (unstable) equilibrium configurations of $`N`$ charges in the complex plane under the influence of a given electric field. This electrostatic analogy was first pointed out by GaudinGAUDIN . It allows the exact access to the thermodynamic limit of the BCS modelGAUDIN-LARGE  by exploiting the formalism of complex analysis. The analytic structure of the ansatz electric field is prescribed by the positions of the charges.

We will now present an electrostatic analogy for the generalized BCS models (1), (2).

In the variables

$$z_\alpha :=\frac{\mathrm{tanh}(pe_\alpha )}{p},x_i:=\frac{\mathrm{tanh}(pu_i)}{p},$$
(19)
the Eqs. (10) are algebraic

$$\frac{q_+^{(0)}}{z_\alpha 1/p}\frac{q_{}^{(0)}}{z_\alpha +1/p}\underset{i=1}{\overset{\mathrm{\Omega }}{}}\frac{1/2}{z_\alpha x_i}+\underset{\genfrac{}{}{0.0pt}{}{\beta =1}{\alpha \beta }}{\overset{N}{}}\frac{1}{z_\alpha z_\beta }=0,$$
(20)
where $`2q_\pm ^{(0)}=1/p\stackrel{~}{g}[\mathrm{\Omega }2(N1)]/2`$. We shall assume that the parameters $`u_i`$ are given by $`u_i=f(\epsilon _i)`$ where $`f`$ is a monotonic function. For the case $`p=\mathrm{i}`$ note that due to the $`\pi `$-periodicity of the transformation (19) the $`u_i`$ can be restricted to lie in $`(\pi /2,\pi /2)`$. Hence, the Bethe equations for the generalized BCS models (10) can be recast in the same form as the original RS equations Eq. (18), except that the constant external electric field of strength $`1/g`$ in the homogeneous case is replaced by the presence of two isolated charges $`q_\pm ^{(0)}`$ fixed at points $`\pm 1/p`$.333The analogy is by no means unique. In fact, a general Moebius transformation $`\mu (z)=re^{i\varphi }(z\omega _1)/(z\omega _2)\omega _1,\omega _2\text{ }\mathrm{C}`$ preserves the algebraic structure of the Bethe equations (20), varying the locations of the charges. Note that for $`p=\mathrm{i}`$ the charges are complex. The interpretation of this is discussed later on.

### III.1 Basics of $`2d`$ electrostatics

We sketch briefly the main ingredients of two-dimensional electrostatics in terms of holomorphic functions in $`\text{ }\mathrm{C}`$ (i.e. real harmonic functions in $`\text{ }\mathrm{R}^2`$).

Let the electric field $`𝐄=(E_x,E_y)`$ be associated to the complex number $`E=E_x\mathrm{i}E_y`$. The Maxwell equations can then be summarized as

$$_{\overline{z}}E=\frac{1}{2}(\mathrm{div}\stackrel{}{E}+\mathrm{i}(\mathrm{curl}\stackrel{}{E})_z)=\pi (\rho +\mathrm{i}(\stackrel{}{j}_{mag})_z)$$
(21)
with real charge density $`\rho `$ and real magnetic monopole current density $`\stackrel{}{j}_{mag}`$ in $`z`$-direction; the derivatives are $`_z=(_x\mathrm{i}_y)/2,_{\overline{z}}=(_x+\mathrm{i}_y)/2`$. Any integral over a closed curve $`\mathrm{\Gamma }`$ gives

$`{\displaystyle \frac{1}{2\pi \mathrm{i}}}{\displaystyle _\mathrm{\Gamma }}𝑑zE(z)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _\mathrm{\Gamma }}(E_xdyE_ydx)\mathrm{i}(E_xdx+E_ydy)`$
$`\widehat{=}`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _\mathrm{\Gamma }}𝐄\mathrm{𝐝𝐬}{\displaystyle \frac{\mathrm{i}}{2\pi }}{\displaystyle _\mathrm{\Gamma }}𝐄\mathrm{𝐝𝐥}`$
$`=`$ $`Q_\mathrm{\Gamma }\mathrm{i}I_\mathrm{\Gamma },`$
where $`Q_\mathrm{\Gamma }`$ and $`I_\mathrm{\Gamma }`$ are the total electric charge and magnetic monopole current enclosed by $`\mathrm{\Gamma }`$ (in the following we refer to $`Q_\mathrm{\Gamma }\mathrm{i}I_\mathrm{\Gamma }`$ as charge). The surface element $`\mathrm{𝐝𝐬}`$ is a vector perpendicularly pointing outwards of $`\mathrm{\Gamma }`$ with the same length as the line element $`\mathrm{𝐝𝐥}`$.

A charge $`q`$ contributes to the electric field with a simple pole with residue $`q`$. A line of charges gives a holomorphic function on a Riemann surface with a branch cut along the charge line. The discontinuity of the field in crossing the cut gives the charge density $`\frac{1}{2\pi \mathrm{i}}\left(E_{}(z)E_+(z)\right)`$, where $`E_{}(z)`$ ($`E_+(z)`$) are the limiting values of the field when $`z`$ tends to the cut from the right (left) with respect to the orientation of the curve.

### III.2 Thermodynamic limit

The thermodynamic limit of the model (1) can be obtained in the following way: we first divide Eq.(20) by $`\mathrm{\Omega }`$ and define the (positive) charge densities

$`\rho (x_j)`$ $`:=`$ $`{\displaystyle \frac{1/2}{\mathrm{\Omega }(x_{j+1}x_j)}}`$
$`\sigma (z_\alpha )`$ $`:=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }|z_{\alpha +1}z_\alpha |}}`$
To obtain a sensible thermodynamic limit, we assume i) that the pairing strength scales as $`\stackrel{~}{g}=G/\mathrm{\Omega }`$, with fixed $`G`$; ii) that the Debye shell defined by the end points $`\epsilon _1`$, $`\epsilon _\mathrm{\Omega }`$ does not depend on $`\mathrm{\Omega }`$; iii) that the number of pairs increases with $`\mathrm{\Omega }`$ according to $`N=\nu \mathrm{\Omega }`$, where $`\nu `$ is the filling. In the limit $`\mathrm{\Omega }\mathrm{}`$, Equations (20) then become

$$\frac{Q_+^{(0)}}{z1/p}\frac{Q_{}^{(0)}}{z+1/p}_L𝑑x\frac{\rho (x)}{zx}+_\mathrm{\Gamma }|dz^{}|\frac{\sigma (z^{})}{zz^{}}=0.$$
(22)
where the integrals are meant in the sense of the principal value $`𝒫`$. After the transformations (19) the transformed Debye shell $`L`$ is still a segment of the real axis having end points $`a_0=\mathrm{tanh}(pu_1)/p`$, $`b_0=\mathrm{tanh}(pu_\mathrm{\Omega })/p`$; the isolated charges are $`Q_\pm ^{(0)}=\frac{1}{2}(\nu 1/2\pm 1/(pG))`$; the density $`\rho `$ is determined by the single-particle energy density $`\rho _\epsilon `$ as $`\rho (x)=\rho _\epsilon (\epsilon (x))/[f^{}(\epsilon (x))(1p^2x^2)]`$, where $`\epsilon (x)=f^1(\text{Arctanh}(px)/p)`$, and it fulfills

$$_{\epsilon _1}^{\epsilon _\mathrm{\Omega }}𝑑\epsilon \rho _\epsilon (\epsilon )=_L𝑑x\rho (x)=1/2.$$

The curve $`\mathrm{\Gamma }`$ and the density $`\sigma `$ have to be determined, with the constraint

$$_\mathrm{\Gamma }|dz|\sigma (z)=\nu .$$

The idea of Gaudin was to construct the total electric field $`E(z)`$ as a function with an analytic structure prescribed by the actual charge distribution. He further assumed that eventual solutions of the equation (22) are arranged in $`K`$ piece-wise differentiable arcs: $`\mathrm{\Gamma }=\mathrm{\Gamma }_\rho _{n=1}^K\mathrm{\Gamma }_n`$, where $`\mathrm{\Gamma }_\rho `$ is the common support of $`\rho (x)`$ and $`\sigma (x)`$ which is not contained in the arcs $`\mathrm{\Gamma }_{arc}_{n=1}^K\mathrm{\Gamma }_n`$. For example, a charge distribution along the line $`[a,b]`$ in the complex plane with total charge $`q`$ and constant charge density would lead to the total electric field $`E(z)=(q/|ab|)\mathrm{ln}[(zb)/(za)]`$ which diverges at the end points of $`[a,b]`$. Since in the case under consideration the electric field has to vanish at the end points of $`\mathrm{\Gamma }`$, we need the charge density to vanish “sufficiently fast” there. An example of such a function could be any $`[(za)(zb)]^\alpha \mathrm{ln}\frac{zb}{za}`$ with any $`\alpha \frac{1}{2}\mathrm{IN}`$ or just $`\left[(za)(zb)\right]^{\frac{\alpha }{2}}`$ with odd $`\alpha `$. Both functions indeed vanish in $`a`$ and $`b`$ and have a branch cut along $`[a,b]`$. The choice of the admissible functions is finally restricted by imposing that the free charges are distributed along a finite set $`\mathrm{\Gamma }`$ GAUDIN-LARGE ; SIERRA-LARGE . In fact, for $`\alpha \frac{1}{2}\mathrm{IN}`$ we either have infinite branch cuts or divergences at the end points. We want to stress that the run of the branch cut can be chosen in an arbitrary way as long as it joins continuously the branch points $`a`$ and $`b`$. The run is fixed requiring that it is an equipotential curve of $`E(z)`$. The ansatz function for the total electric field is

$$\begin{array}{cc}\hfill E(z)=& S(z)\left[\frac{Q_+}{z1/p}+\frac{Q_{}}{z+1/p}+_L𝑑x\frac{\phi (x)}{zx}\right],\hfill \\ \hfill S(z)=& \underset{n=1}{\overset{K}{}}\sqrt{(za_n)(zb_n)}.\hfill \end{array}$$
(23)
Note that the sign of the complex square root is uniquely given, once $`\mathrm{\Gamma }_{arc}`$ is fixed, by $`\sqrt{(za_n)(zb_n)}=\sqrt{|za_n||zb_n|}\mathrm{exp}\frac{\mathrm{i}}{2}\mathrm{arg}_{\mathrm{\Gamma }_n}(za_n)(zb_n)`$. The argument function $`\mathrm{arg}_{\mathrm{\Gamma }_n}(za_n)(zb_n)`$ is the sum of the angles between $`z\text{ }\mathrm{C}`$ and the two end points of $`\mathrm{\Gamma }_n`$, which has a $`2\pi `$–discontinuity crossing $`\mathrm{\Gamma }_n`$. The unknown quantities are the end points $`(a_n,b_n)`$ of $`\mathrm{\Gamma }_{arc}`$, $`Q_\pm `$ and $`\phi (x)`$.
To give a physical interpretation of such an ansatz function, look at $`S(z)Q_\pm `$ and $`S(z)\phi (x)`$ as screened charges and charge density, respectively. From the symmetry of the distribution of the fixed charges on $`L`$ and in $`\pm 1/p`$ it follows that if $`a_n`$ has non-zero imaginary part, then $`b_n=\overline{a}_n`$, and we choose the orientation of $`\mathrm{\Gamma }_n`$ from $`a_n`$ to $`b_n`$, with $`\mathrm{Im}(a_n)<0`$. If $`a_n`$ is real, so must be $`b_n`$. We argue however that this occurs only for distinct “critical” values of the pairing coupling at which the imaginary part of complex conjugate $`a_n`$ and $`b_n`$ vanishes and hence $`a_n=b_n`$. The corresponding arc is then describing a closed curve. Therefore we consider $`S(z)`$ to don’t have cuts along the transformed Debye shell $`L`$ (see Appendix D). We emphasize that this does not mean that $`\sigma (z)`$ has no support on $`L`$. A non-zero solution charge density on $`L`$ is accounted for by the screening density $`S(x)\phi (x)`$.

After having determined the $`a_n,b_n`$, the curves $`\mathrm{\Gamma }_n`$ are found as equipotential curves of the total electric field Eq.(23). The density $`\sigma `$ is then determined by the discontinuity of $`E`$ in crossing $`\mathrm{\Gamma }`$.

In the following, all the unknown quantities are obtained following a procedure à la Gaudin on which we report in the Appendix C. We will calculate the complex integrals involved in this procedure by exploiting basic knowledge of electrostatics.

Since the screened charges must tend to their bare values when the point $`z`$ is sufficiently close to the source, we have that

$`S(\pm 1/p)Q_\pm `$ $`=Q_\pm ^{(0)},`$ (24a)
In a similar way, using $`\sigma (z)\rho (z)=(E_{}(z)E_+(z))/(2\pi \mathrm{i})`$ , we find
$`S(x)\phi (x)`$ $`=\sigma (x)\rho (x).`$ (24b)
Next, imposing that the electric field asymptotically goes as $`E(z)1/z^2`$ (this is the leading term in the multi-pole expansion, since the total charge is zero), we obtain the following equations

$`{\displaystyle _L}𝑑xx^n\phi (x)={\displaystyle \frac{1}{p^n}}\left(Q_++(1)^nQ_{}\right),`$ $`0nK.`$ (25)
We present a detailed derivation of these equations in Appendix C. There we also show explicitly that the condition $`_\mathrm{\Gamma }\sigma (z)|dz|=\nu `$ is automatically fulfilled with Eqs.(25).

Using Eqs.(24) we can express $`Q_\pm `$ and $`\phi (x)`$ in terms of $`S(\pm 1/p)`$ and $`S(x)`$ respectively. The $`\mathrm{\Gamma }_n`$ are finally obtained as equipotential curves of the total electric field $`E(z)`$GAUDIN-LARGE ; SIERRA-LARGE

$$\mathrm{Re}_{a_n}^z𝑑z^{}E_{}(z^{})=E_x𝑑x+E_ydy=0,$$
(26)
The density $`\sigma `$ is determined by the discontinuity of $`E`$ along $`\mathrm{\Gamma }`$, according to:

$$\sigma (z)=\rho (z)+\frac{1}{\pi }|E_{}(z)|.$$
(27)
It is worth noting that Eqs.(25) are a set of $`K+1`$ real equations for the $`2K`$ real parameters determining the end points of $`\mathrm{\Gamma }_n`$. Gaudin stated that seemingly in any case there is a finite number of solutions. From a physical perspective, though, we expect that the $`K1`$ free parameters span a family of curves corresponding to a band of excitations of the system. This conjecture will be verified in forthcoming workFUTURE .

We summarize all the conditions found for the unknowns $`a_n`$, $`b_n`$, $`\phi `$, $`Q`$

$$\begin{array}{ccc}\hfill S(z)& =& \underset{n=1}{\overset{K}{}}\sqrt{(za_n)(zb_n)},\hfill \\ & & \\ \hfill E(z)& =& S(z)\left[\frac{Q_+}{z1/p}+\frac{Q_{}}{z+1/p}+\underset{L}{}𝑑x\frac{\phi (x)}{zx}\right],\hfill \\ & & \\ \hfill Q_\pm & =& \frac{\nu 1/2\pm 1/pG}{2S(\pm 1/p)}\hfill \\ & & \\ \hfill \phi (x)& =& \frac{\sigma (x)\rho (x)}{S(x)}\hfill \\ & & \\ \hfill \sigma (z)& =& \rho (z)+\frac{1}{\pi }|E_{}(z)|\hfill \end{array}$$
The arcs are then determined by
$$\underset{L}{}𝑑xx^n\frac{\rho (x)\sigma (x)}{S(x)}=\frac{Q_++()^nQ_{}}{p^n};n=0,\mathrm{},K$$

$$\mathrm{\Gamma }_n:\mathrm{Re}\underset{a_n}{\overset{z}{}}𝑑z^{}E_{}(z^{})=0.$$

As we stated above, the density $`\rho `$ of the transformed variables $`x=\mathrm{tanh}(pf(\epsilon ))/p`$ is connected with the given density $`\rho _\epsilon `$ of the single-particle levels via

$$\begin{array}{cc}\hfill \rho (x)& =\rho _\epsilon (\epsilon (x))/[f^{}(\epsilon (x))(1p^2x^2)],\hfill \\ \hfill \epsilon (x)& =f^1(\text{Arctanh}(px)/p).\hfill \end{array}$$

The energy of the state is finally given by

$`=`$ $`\alpha \nu ^2G{\displaystyle _L}𝑑x\rho (x)\epsilon (x){\displaystyle _\mathrm{\Gamma }}|dz|\sigma (z){\displaystyle \frac{1p^2xz}{zx}},`$ (28)
$`=`$ $`\alpha \nu ^2+p^2G\nu {\displaystyle _{\epsilon _1}^{\epsilon _\mathrm{\Omega }}}𝑑\epsilon \rho _\epsilon (\epsilon )\epsilon x(\epsilon )+`$
$`\text{ }G{\displaystyle _L}𝑑x{\displaystyle \frac{\rho _\epsilon (\epsilon (x))\epsilon (x)}{f^{}(\epsilon (x))}}{\displaystyle _\mathrm{\Gamma }}|dz|{\displaystyle \frac{\sigma (z)}{zx}},`$ (29)
where $`\alpha =A\mathrm{\Omega }`$.

#### III.2.1 Comparison with the thermodynamic limit of the BCS model

In the limit $`p0`$, the two isolated charges $`Q_\pm ^{(0)}`$ take an infinite value and they are displaced to infinity in such a way to give a uniform electric field $`1/G`$. In this limit, the charges $`Q_\pm `$ behave as

$$\begin{array}{cc}& Q_+\pm \left[\frac{1}{2G}\frac{1}{2}\left(\nu \frac{1}{2}+\frac{1}{2G}\underset{n=1}{\overset{K}{}}(a_n+b_n)\right)p\right]p^{K1},\hfill \\ & Q_{}\left[\frac{1}{2G}\frac{1}{2}\left(\nu \frac{1}{2}+\frac{1}{2G}\underset{n=1}{\overset{K}{}}(a_n+b_n)\right)p\right](p)^{K1},\hfill \end{array}$$

where we kept into account the dependence of the relative determination of $`S(\pm 1/p)`$ upon the number of arcs. Thus, the ansatz field reduces to

$$E_{p0}(z)E_0(z)=S(z)_L𝑑x\frac{\phi (x)}{zx},$$
(30)
Eq.(24b) still hold, while Eqs.(25) simplify to:

$`{\displaystyle _L}𝑑xx^n\phi (x)=0,0nK2,`$ (31a)
$`{\displaystyle _L}𝑑xx^{K1}\phi (x)={\displaystyle \frac{1}{G}},`$ (31b)
$`{\displaystyle _L}𝑑xx^K\phi (x)=\nu +{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2G}}{\displaystyle \underset{n=1}{\overset{K}{}}}(a_n+b_n),`$ (31c)
The last two equations come from the fact that the field should go asymptotically as $`E_0(z)1/G+Q_T/z`$ in this limit, since the total charge is now $`Q_T=\nu 1/2`$. In the approach by Gaudin the last condition was obtained by imposing $`|dz|\sigma (z)=\nu `$. This extracts the residue at infinity, which has to be such that the total charge is zero. In our approach, we already started from a globally neutral system and hence this normalization condition for $`\sigma (z)`$ was automatically fulfilled (see Appendix C)

We compare our findings for integrable pairing models with the results for the BCS model obtained in Refs. GAUDIN-LARGE, ; RICH77,  and reconsidered recently by Sierra and coworkers SIERRA-LARGE

$$\begin{array}{ccc}\hfill S(z)& =& \underset{n=1}{\overset{K}{}}\sqrt{(za_n)(zb_n)},\hfill \\ & & \\ \hfill E_0(z)& =& S(z)\underset{L}{}𝑑x\frac{\phi (x)}{zx},\hfill \\ & & \\ \hfill \phi (x)& =& \frac{\sigma (x)\rho (x)}{S(x)};\sigma (z)=\rho (z)+\frac{1}{\pi }|E_{}(z)|\hfill \end{array}$$

The arcs are then determined by

$$\begin{array}{c}\underset{L}{}𝑑xx^n\frac{\sigma (x)\rho (x)}{S(x)}=0,n=0,\mathrm{},K2\hfill \\ \\ \underset{L}{}𝑑xx^{K1}\frac{\rho (x)\sigma (x)}{S(x)}=\frac{1}{G}\hfill \\ \\ \underset{L}{}𝑑xx^K\frac{\rho (x)\sigma (x)}{S(x)}\frac{1}{2G}\underset{n=1}{\overset{K}{}}(a_n+b_n)=\nu \frac{1}{2}\hfill \end{array}$$

$$\mathrm{\Gamma }_n:\mathrm{Re}\underset{a_n}{\overset{z}{}}𝑑z^{}E_{0}^{}{}_{}{}^{}(z^{})=0.$$

#### III.2.2 Numerics and discussion of the solutions

To solve the BA equations (10), (20) we have chosen the parameters $`u_i=\epsilon _i`$ at half filling $`\mathrm{\Omega }=2N`$. The corresponding interactions $`U_{ij}`$ and $`g_{ij}`$ are shown in Fig.1. In the trigonometric case the Debye shell ranges in $`(1,1)`$; in the hyperbolic case it ranges in $`(0,2)`$ for numerical convenience (see the caption of Fig.2). The ground state of the system (which is the only state we consider in the numerics) is obtained by evolving the $`G=0`$ state

$$\underset{u0}{lim}z_\alpha =x_\alpha \alpha (1,\mathrm{},N).$$
(32)
to a finite value of $`G`$. At certain value of $`G`$ the equations are singular since a couple of pairing parameters $`z_{2\lambda 1}`$ and $`z_{2\lambda }`$ coincide with the energy level $`x_{2\lambda 1}`$. Such a singular behaviour can be smoothed out using the procedure developed by Richardson in Ref. RICH-NUM, . In the following we briefly summarize it. The corresponding divergences are removed by the following transformations

$$z_{2\lambda 1}=A_\lambda iB_\lambda ,z_{2\lambda }=A_\lambda +iB_\lambda ,$$
(33)
and then

$$X_\lambda =A_\lambda \frac{x_{2\lambda 1}+x_{2\lambda }}{2},Y_\lambda =\frac{B_\lambda ^2}{\delta _\lambda ^2X_\lambda ^2}.$$
(34)
where $`\delta _\lambda =(x_{2\lambda }x_{2\lambda 1})/2`$ is half the energy spacing between the corresponding energy values (with whom they coincide if the pairing coupling strength $`u`$ is zero).

By the first transformation Eqs. (20) can be written in a form which is manifestly real, whereas the second transformation (34) removes the divergences.

Now we discuss some features of the arcs shown in Figs.2 and 3. For $`K=1`$ (we recall that straight lines on the real axis are not counted by the index $`n`$), the end points are uniquely determined (if they exist) and the curve $`\mathrm{\Gamma }`$ corresponds to the ground state444With ground state we mean that state connected with the ground state for $`G=0`$. Up to some “critical” value of $`G`$ this is the ground state of the system. of the system. In this case the end points of the arc, $`\lambda \pm i\mathrm{\Delta }`$, are determined by the two coupled equations (obtained as linear combinations of Eqs. (25))

$`{\displaystyle _L}𝑑x{\displaystyle \frac{(1+px)\rho (x)}{\sqrt{(x\lambda )^2+\mathrm{\Delta }^2}}}`$ $`=2Q_+={\displaystyle \frac{1/pG+\nu 1/2}{\sqrt{(\lambda 1/p)^2+\mathrm{\Delta }^2}}}`$
$`{\displaystyle _L}𝑑x{\displaystyle \frac{(1px)\rho (x)}{\sqrt{(x\lambda )^2+\mathrm{\Delta }^2}}}`$ $`=2Q_{}={\displaystyle \frac{1/pG+\nu 1/2}{\sqrt{(\lambda +1/p)^2+\mathrm{\Delta }^2}}}`$
Why this set of equation corresponds to the ground state is discussed in some detail in Appendix D.

In the trigonometric case ($`p=\mathrm{i}`$), the equations are both complex, and they are conjugate. Thus, in order to find $`\lambda `$ and $`\mathrm{\Delta }`$ it suffices to solve separately real and imaginary part of one of them. In the hyperbolic case ($`p=1`$), assuming a uniform single-particle energy density $`\rho _\epsilon (\epsilon (x))=\rho _0`$ (i.e. $`\rho (x)=\rho _0/(1x^2)`$) and a Debye shell spanning the interval $`[0,2\omega _D]`$, we have that the system above admits a real solution (and hence two coinciding end-points)

$$\begin{array}{cc}\hfill \mathrm{\Delta }=& 0\hfill \\ \hfill \lambda =& \mathrm{sinh}(2\omega _D)/\left(\mathrm{cosh}2\omega _D\mathrm{exp}\left[(\nu 1/2)/\rho _0\right]\right),\hfill \end{array}$$

if $`1/(\rho _0G)=2\omega _D`$. This corresponds to the external curve in Fig. (2).

Situations corresponding to $`K>1`$ lead to excited states. We point out that the independent equations (25) are in number of $`K+1`$ leaving undetermined $`K1`$ parameters. This might be an evidence of “bands” of excitations above the ground stateFUTURE .

Looking at the behaviour of the quasi–momenta solutions of the BA equations, we can note that the more intense is the pairing constant $`G`$ the more evident is the tendency of the quasi–energies to be complex. In Fig. 4 we compare this tendency for the hyperbolic, trigonometric, and uniform BCS models. We found that the pairing tendency increases from the trigonometric over the uniform to the hyperbolic model. Looking at the explicit couplings plotted in Fig. 1, we should conclude that the pairing and Coulomb interaction in Eqs. (2) are competitive couplings.

## IV Summary and Conclusions

We have discussed the electrostatic analogy for the exactly solvable pairing force Hamiltonians found in Ref. ADO, . These models generalize the BCS model in that the coupling constants are not uniform (in the space of quantum labels $`i,j`$, see Eqs. (1)). The ordinary BCS model is obtainable as a limit ADO  $`p0`$ of the models (1), (2). The electrostatic analogy exists because the BA equations for the “rapidities” $`e_\alpha `$ can be recast in a special algebraic form. For the class of models we deal with such an algebraic form is obtained by the change of variables: $`z_\alpha =\mathrm{tanh}(pe_\alpha )`$, $`x_i=\mathrm{tanh}(p\epsilon _i)`$. After this transformation the Bethe equations express the condition for an equilibrium alignment of mobile charges $`q_{mobile}=+1`$ in the $`2d`$–real or complex plane in a neutralizing background of $`q_{fixed}=1/2`$ fixed charges placed in the positions $`x_i`$ and two further charges in the positions $`\pm 1/p`$. Remarkably, this analogy is very effective to obtain the exact thermodynamic limit (i.e. the large $`N`$ limit) of the models. In this limit the mobile charges (which are the solutions of the Bethe equations) arrange themselves along equipotential curves $`\mathrm{\Gamma }=_n^K\mathrm{\Gamma }_n`$ of the total electric field with vanishing charge density at its extremities. The electrostatic analogy presented in this work is by no means unique; the charge distribution can be transformed by the whole group of conformal Moebius transformations. We have obtained explicitly the thermodynamic limit for the ground state configuration at half–filling (see Fig. (2), (3)). We note, however, that the approach is also valid for excited states of the system, though the explicit calculation can involve technical subtleties (discussed in Appendix D).
The electrostatic analogy together with the constructive equations for the BCS model found by Gaudin and Richardson is demonstrated to be obtained from the limit $`p0`$ of the results presented here. A comparison between the hyperbolic, trigonometric, and uniform cases is seen in Fig. (4). Since this represent an inversion of the behaviour one would have expected from Fig. (1) we conclude that the pairing and Coulomb interactions in Eqs. (2) are competitive couplings. The exact thermodynamics of integrable pairing models is one of the future goals.

###### Acknowledgements.

We thank G. Sierra for helpful discussions. We further acknowledge constant support by G. Falci, R. Fazio, and G. Giaquinta.

## Appendix A The Gaudin algebra

The Gaudin algebra $`𝒢[sl(2)]`$ is constructed from $`sl(2)`$. The $`sl(2)`$ “lowest” weight module is generated by the vacuum vector $`|0_j`$, $`K_j^{}|0_j=0,K_j^z|0_j=k_j|0_j`$ where $`k_j`$ is the “lowest” weight ($`k_j=1/2`$ for spin $`1/2`$ which is the case considered for electrons). The infinite dimensional $`𝒢[sl(2)]`$ is generated by

$$K^\pm (\xi ):=\underset{j=1}{\overset{\mathrm{\Omega }}{}}\varphi (\xi u_j)K_j^\pm ,K^z(\xi ):=\underset{j=1}{\overset{\mathrm{\Omega }}{}}\psi (\xi u_j)K_j^z$$
(35)
where $`\xi \text{ }\mathrm{C}`$ and $`u_j\text{ }\mathrm{R}`$. The module of $`𝒢[sl(2)]`$ is characterized by the vacuum $`|0_{j=1}^\mathrm{\Omega }|0_j`$: $`K^{}(\xi )|0=0,K^z(\xi )|0=\kappa (\xi )|0,`$ where $`\kappa (\xi ):=_{j=1}^\mathrm{\Omega }k_j\psi (\xi u_j)`$ is the lowest weight of $`𝒢[sl(2)]`$. The element of $`𝒢[sl(2)]`$ corresponding to the $`su(2)`$ Casimir operator is $`K(\xi ):=K^z(\xi )K^z(\xi )+\frac{1}{2}\left(K^+(\xi )K^{}(\xi )+K^{}(\xi )K^+(\xi )\right)`$. The generating function of the integral of the motion of the pairing models is related to $`K(\xi )`$ through

$$\tau (\xi )=2\mathrm{\Lambda }K(\xi )+K^z(\xi )+c(\xi )1\mathrm{l}$$
(36)
where $`c(\xi )`$ is a $`\text{ }\mathrm{C}`$-function and $`\mathrm{\Lambda }`$ is a real parameter. The property $`[\tau (v),\tau (w)]=0`$ arises from the quasi-classical limit of the $`sl(2)`$ QIS theory; this is the ultimate reason for the integrability of the BCS model.

In the present paper $`\varphi `$ and $`\psi `$ are either rational or trigonometric/hyperbolic functions.

### A.1 Rational $`𝒢[sl(2)]`$

In this case: $`\varphi (\xi )=\psi (\xi ):=1/\xi `$, $`\kappa (\xi )k_0(\xi )=_{j=1}^\mathrm{\Omega }k_j/(\xi u_j)`$, and $`\mathrm{\Lambda }g`$. The operators (35) obey

$`[K^z(v),K^\pm (w)]`$ $`=`$ $`{\displaystyle \frac{K^\pm (v)K^\pm (w)}{vw}},`$
$`[K^{}(v),K^+(w)]`$ $`=`$ $`2{\displaystyle \frac{K^z(v)K^z(w)}{vw}},`$
where $`vw\text{ }\mathrm{C}`$. The Bethe equations are the Richardson-Sherman (RS) equations

$$k_0(e_\alpha )=\frac{1}{g}+\underset{\genfrac{}{}{0.0pt}{}{\beta =1}{\beta \alpha }}{\overset{N}{}}\frac{1}{e_\beta e_\alpha },\alpha =1,\mathrm{},N.$$
(37)
We note that RS equations (37) are intimately related to the algebraic structure of $`𝒢[sl(2)]`$ since they act as constraints on the lowest weight $`k(E_\alpha )`$. The difference between the BCS and Gaudin model results in a different constraint imposed on the lowest weight vector of $`𝒢[sl(2)]`$ which leads to different sets $``$, $`^{}`$ of solutions of the BA equations ($`^{}`$ is spanned by the solutions of (37) when $`g\mathrm{}`$). This fact has been used to extend the Sklyanin theorem for the Gaudin models to the BCS modelAO .

### A.2 Hyperbolic/Trigonometric $`𝒢[sl(2)]`$

In this case: $`\varphi (\xi ):=p/\mathrm{sinh}[p\xi ]`$, $`\psi (\xi ):=p\mathrm{coth}[p\xi ]`$ , $`\kappa (\xi )k(\xi )=_{j=1}^\mathrm{\Omega }pk_j\mathrm{coth}[p(\xi u_j)]`$and $`\mathrm{\Lambda }\stackrel{~}{g}`$. The operators (35) obey

$$\begin{array}{c}[K^z(v),K^\pm (w)]=p\frac{K^\pm (v)\mathrm{cosh}[p(vw)]K^\pm (w)}{\mathrm{sinh}[p(vw)]},\hfill \\ \\ [K^{}(v),K^+(w)]=2p\frac{K^z(v)K^z(w)}{\mathrm{sinh}[p(vw)]},\hfill \end{array}$$

where $`vw\text{ }\mathrm{C}`$.

The BA equations for the corresponding Hamiltonian (1) are

$$k(e_\alpha )=\frac{1}{\stackrel{~}{g}}+\underset{\genfrac{}{}{0.0pt}{}{\beta =1}{\beta \alpha }}{\overset{N}{}}p\mathrm{coth}[p(e_\beta e_\alpha )],\alpha =1,\mathrm{},N.$$
(38)
Also in this case the BA equations act as constraints on the lowest weight $`k(e_\alpha )`$ of the trigonometric $`𝒢[sl(2)]`$.

## Appendix B The Richardson route and the Riccati equation

In this appendix we first emphasize the key role played by the non linear differential equation of Riccati type in the electrostatic analogy of exactly solvable pairing models.RICCATI  Then we sketch the procedure originally employed by Richardson RICH77  to obtain the thermodynamic limit of the BCS model.

In the non-uniform case the electric field is

$$E(z):=\frac{1}{2}k(z)+\frac{1}{2\stackrel{~}{g}}+\frac{p}{2}\underset{\genfrac{}{}{0.0pt}{}{\beta =1}{\beta \alpha }}{\overset{N}{}}\mathrm{coth}[p(e_\beta z)]$$
(39)
with $`k(z)`$ as defined in the previous section. In order to have the field obeying a Riccati-type equation, it is crucial to map one of the isolated charges to infinity. This is done by the transformation

$$q_i:=\mathrm{exp}(2p\epsilon _i);\zeta _\alpha :=\mathrm{exp}(2pe_\alpha ).$$
(40)
In the variables $`q_i`$ and $`\zeta _\alpha `$, the Riccati equation is

$$\frac{dE(z)}{dz}+E^2(z)=\stackrel{~}{\tau }(z)+\stackrel{~}{f}(z)$$
(41)
where $`\stackrel{~}{\tau }_0(z):=1/\stackrel{~}{g}_{j=1}^\mathrm{\Omega }\stackrel{~}{\tau }_j/[q_j(q_jz)]`$ is a generating function of the eigenvalues of integrals of the motion of the model (1) (see Eqs. (9)) and

$`\stackrel{~}{f}(z)=`$ $`{\displaystyle \frac{Q_0(Q_0+1)}{\stackrel{~}{g}^2}}`$
$`+{\displaystyle \underset{j=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1}{q_jz}}\left[{\displaystyle \frac{k_j(k_j+1)}{q_jz}}2{\displaystyle \frac{Q_0}{z}}2{\displaystyle \frac{P_0}{q_j}}\right]`$
where $`P_0:=2Q_0+\mathrm{\Omega }2N+1/2`$.

For the isotropic limit $`p0`$ the BCS rational case is recovered. The electrostatic field is here given by (see RS Eqs. (37))

$`E_0(z)`$ $`=`$ $`k_0(z)+{\displaystyle \frac{1}{g}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{\beta =1}{\beta \alpha }}{\overset{N}{}}}{\displaystyle \frac{1}{ze_\beta }}`$
$`=`$ $`{\displaystyle \frac{1}{g}}+{\displaystyle \underset{j=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{z\epsilon _j}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{\beta =1}{\beta \alpha }}{\overset{N}{}}}{\displaystyle \frac{1}{ze_\beta }}.`$
In this case the Riccati equation reads

$$\frac{dE_0(z)}{dz}+E_0^2(z)=\tau (z)+f(z)$$
(44)
where

$$f(z)=\frac{N^2}{g^2}+\underset{j=1}{\overset{\mathrm{\Omega }}{}}\frac{k_j(k_j+1)}{(z\epsilon _j)^2}$$
(45)
where $`\tau (z)`$ is the generating function of the eigenvalues of integrals of the motion of the BCS model Eqs. (17). The Richardson BA equations are the zeros of the solutions of (45). The role of the Riccati equations in connection with the Gaudin models was also investigated in Ref.RICCATI, .

In both the rational or non-uniform cases it is evident that the Riccati equation plays an important role. This fact has been used by Richardson to derive the BCS theory from the expansion of the solution of the Riccati equation (in the rational case and $`k_j=1/2j`$) in powers of $`1/N`$. In the following we briefly summarize the Richardson’s physical arguments. The Richardson equations, via the electrostatic analogy, can be seen as a self-consistent evaluation of the field of the fixed charges in the presence of the screening due to the mobile charges. The strategy followed by Richardson is to attempt to eliminate any reference to the mobile charges in the right hand side of the above equation. In fact, using the Eqs. (37) the Riccati equation for the field can be recast in the following form:

$`{\displaystyle \frac{dE_0(z)}{dz}}+E_0^2(z)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{(zϵ_i)^2}}+\left[{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{zϵ_i}}+{\displaystyle \frac{1}{g}}\right]^2`$ (46)
$`2{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{zϵ_i}}{\displaystyle \underset{\beta =1}{\overset{N}{}}}{\displaystyle \frac{1}{ϵ_ie_\beta }}.`$
In the last term on the right hand side there appears the field due to the mobile charges at the location of the fixed charges. The self-consistency condition, and the effective screening, enters here. Richardson notices that the field of the mobile charges may be written in the following way

$$\frac{1}{ϵ_ie_\beta }=\frac{1}{2\pi \mathrm{i}}_Cdz\frac{E_0(z)}{ϵ_iz}$$
(47)
with the curve $`C`$ going around the singularities of $`E(z)`$ due to the mobile charges. By using this, one gets a integro-differential equation of the Riccati type

$`{\displaystyle \frac{dE_0(z)}{dz}}+E_0^2(z)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{(zϵ_i)^2}}+\left[{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{zϵ_i}}+{\displaystyle \frac{1}{g}}\right]^2`$
$`+`$ $`{\displaystyle \frac{1}{2\pi \mathrm{i}}}`$ $`{\displaystyle _C}dz^{}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{E_0(z^{})}{(zϵ_i)(z^{}ϵ_i)}}.`$ (48)
The nice feature of the above equation is that any reference to the mobile charges has disappeared. One may object that in the above equation there is still a reference to the mobile charges via the curve $`C`$ entering the contour integral. In fact, by knowing the property that the field has to satisfy at infinity one may observe that the integral over a closed curve may be written as the sum of two curves: the first is $`C`$ enclosing the singularities of the mobile charges, while the second enclosing the singularities of the fixed charges. Furthermore to evaluate the contour integral along the curve extending at infinity, it is enough to know the first term in the multi-pole expansion of the field. The final result is that the integral over $`C`$ may be written in terms of the curve enclosing the singularities of the fixed charges and the behavior of the field at infinity only. The elimination of the mobile charges from the effective equation for the field is closely in the spirit of the Thomas-Fermi approach to the evaluation of the effective potential due to the combined effect of fixed charges (nucleus in atoms, ions in metals) and mobile charges (electrons). The integro-differential Riccati equation plays in this context the role of the Poisson equation, which becomes the Thomas-Fermi equation once the charge of the mobile electron is written in terms of the potential.

The next trick used by Richardson is to expand the complex electric field $`E_0(z)`$ near the point at infinity, in powers of $`1/z`$. This amounts to a multi-pole expansion. The constant electric field gives the zeroth order term. The other terms may be expanded upon using the geometric series. One easily gets

$`E_0^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{g}}`$ (49)
$`E_0^{(1)}`$ $`=`$ $`N{\displaystyle \frac{1}{2}}\mathrm{\Omega }`$ (50)
$`E_0^{(2)}`$ $`=`$ $`{\displaystyle \underset{\beta =1}{\overset{N}{}}}e_\beta {\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}ϵ_i{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}ϵ_i.`$
The monopole and dipole terms give then the number of pairs and the energy. To perform the thermodynamic limit, we rescale the coupling constant as $`gG/\mathrm{\Omega }`$ in order that the energy remains an extensive quantity in this limit (see section III B). The field is expanded as $`E_0=E_N+E_1+E_{1/N}+\mathrm{}`$ where the subscript indicated terms of order $`N`$, $`1`$, and $`1/N`$, respectively. We are assuming that $`N`$ and $`\mathrm{\Omega }`$ keep a fixed ratio. One may rewrite the Riccati equation order by order

$`E_N^2`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{G}}+{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{zϵ_i}}\right)^2{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{H_{Ni}}{zϵ_i}}`$ (51)
$`2E_NE_1`$ $`=`$ $`{\displaystyle \frac{dE_N}{dz}}+{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1/2}{(zϵ_i)^2}}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{H_{1i}}{zϵ_i}}`$
where at each order in $`N`$

$$H_i=\frac{1}{2\pi \mathrm{i}}_Cdz\frac{E_0(z)}{ϵ_iz}.$$
(52)
Following the analogous lines presented in section III, Richardson wrote down the solution for $`E_N`$ as

$$E_N=\frac{1}{2}Z(z)\underset{i=1}{\overset{\mathrm{\Omega }}{}}\frac{1}{Z(ϵ_i)(zϵ_i)}$$
(53)
with $`Z(z)=\sqrt{(za)(zb)}`$ having a branch cut between the points $`a`$ and $`b`$. To see how the BCS limit is recovered, one has to expand $`E_N`$ in multiple’s and insert the $`E_N^{(m)}`$ back into (49). A direct calculation shows that

$`E_N^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1}{Z(ϵ_i)}}`$ (54)
$`E_N^{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1}{Z(ϵ_i)}}\left(ϵ_i{\displaystyle \frac{a+b}{2}}\right)`$ (55)
$`E_N^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}{\displaystyle \frac{1}{Z(ϵ_i)}}\left(ϵ_i^2ϵ_i{\displaystyle \frac{a+b}{2}}{\displaystyle \frac{1}{8}}(ab)^2\right).`$
By setting $`a=\lambda +\mathrm{i}\mathrm{\Delta }`$ and identifying the real and imaginary parts of $`a`$ with the chemical potential and the energy gap, the first two equations of (54) reproduce the gap equation and the condition for the chemical potential, while the third gives the ground state energy. This completes the evaluation of the field in the leading order in $`1/N`$. To proceed further, one has to solve the second of (51), which is complicated by the way the field $`E_1`$ enters via the quantities $`H_{1i}`$. Richardson showed, however, that there is an elegant way to circumvent the tackling of the integral equation. An important point to notice is that the singularities of $`E_N`$ exhaust all the charges, so that the correction $`E_1`$ may have no poles at the positions of the fixed charges. Inspection of the second of Eq.(51) shows that $`E_1`$ has poles at the zeros of $`E_N`$. However, $`E_1`$ cannot have poles because this would imply further charges. To avoid this one has to impose that the right hand side of the second equation of (51) has to vanish at the zeros of $`E_N`$. There are $`\mathrm{\Omega }1`$ zeros of $`E_N`$ and one gets $`\mathrm{\Omega }1`$ linear equations for the $`\mathrm{\Omega }`$ unknown $`H_{1i}`$. One more equation may be obtained by performing a multi-pole expansion in the Riccati equation. The Riccati equation for $`1/z`$ terms reads

$$2E_0^{(0)}E_0^{(1)}=\frac{\mathrm{\Omega }^2}{2G}\underset{1=1}{\overset{\mathrm{\Omega }}{}}H_i$$
(56)
which may be rewritten as

$$N=\frac{G}{\mathrm{\Omega }}\underset{i=1}{\overset{\mathrm{\Omega }}{}}H_i.$$
(57)
By expanding in powers of $`1/N`$ one gets

$`N`$ $`=`$ $`{\displaystyle \frac{G}{\mathrm{\Omega }}}{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}H_{Ni}`$ (58)
$`0`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{\mathrm{\Omega }}{}}}H_{1i}.`$ (59)
The second of the above equation provides the $`\mathrm{\Omega }`$th linear equation for the determination of the unknown $`H_{1i}`$. These latter, once determined, may be inserted back into the second of (51) to obtain the correction $`E_1`$.

## Appendix C Electric field from contour integration

In this appendix we present the derivation of Eqs.(24) and (25) exploiting standard results of complex analysis. Then we calculate the residue at infinity. Finally we discuss some subtleties associated with the normalization condition for the charge density $`\sigma (z)`$.

1. The contour integral $`{\displaystyle 𝑑zE(z)/(zz^{})}`$

Our first aim is to write the last term in Eq.(22) in terms of the ansatz field $`E`$ defined by Eq.(23). To do this, we consider the contour $`C`$ shown in Fig.(5).

From the generalized Cauchy theorem, we know that the sum of the contour integrations of $`E(z^{})/(zz^{})`$ inside $`C`$ is equal to the negative residue of $`E(z^{})`$ at infinity (i.e. the charge at infinity), which is zero in our case.

The contours surrounding $`\mathrm{\Gamma }`$ can be contracted to an integration along the right-hand side minus an integration along the left-hand side of $`\mathrm{\Gamma }`$. Since the discontinuity of $`E(z)`$ at $`z\mathrm{\Gamma }`$ equals the charge density at $`z`$ times $`2\pi \mathrm{i}`$, this gives

$$\frac{1}{2\pi \mathrm{i}}_{C_\mathrm{\Gamma }}dz\frac{E(z^{})}{zz^{}}=_\mathrm{\Gamma }|\mathrm{d}z|\frac{\sigma (z^{})\rho (z^{})}{zz^{}}.$$
(60)
The integrals along the small half-circles at the end-points of $`\mathrm{\Gamma }`$ tend to zero. If we don’t have an arc of $`\mathrm{\Gamma }`$ on $`L`$ defined within the function $`S(z)`$, we obtain

$$\begin{array}{cc}\hfill 0=& \frac{1}{2\pi \mathrm{i}}_C\frac{E(z^{})}{zz^{}}𝑑z^{}=\frac{S(1/p)Q_+}{z1/p}+\frac{S(1/p)Q_{}}{z+1/p}\hfill \\ & +_\mathrm{\Gamma }|\mathrm{d}z|\frac{\sigma (z^{})\rho (z^{})}{zz^{}}+_L\frac{S(x)\phi (x)}{zx}𝑑x\stackrel{!}{=}0.\hfill \end{array}$$
(61)
This finally leads to

$$\begin{array}{cc}\hfill _\mathrm{\Gamma }\frac{\sigma (z^{})}{zz^{}}|dz^{}|& =\frac{S(1/p)Q_+}{z1/p}\frac{S(1/p)Q_{}}{z+1/p}+\hfill \\ & _{L\mathrm{\Gamma }}\frac{S(x)\phi (x)}{zx}𝑑x+_{L\mathrm{\Gamma }}\frac{\rho (x)}{zx}𝑑x.\hfill \end{array}$$
(62)
By substituting Eq. (62) in the Bethe equations Eq.(22), we find the conditions in Eqs.(24).

Note that $`\sigma (z)`$ and $`\rho (z)`$ can have common support even though the screening function $`S(z)`$ has no cut along $`L`$. $`L\mathrm{\Gamma }`$ here means the common support of $`\sigma (z)`$ and $`\rho (z)`$.

2. Asymptotic behavior of $`E(z)`$ at $`z=\mathrm{}`$

In order to find the remaining equations (25), we expand the electric field in inverse powers of $`z`$, and impose that the leading order is $`1/z^2`$.

$$E(z)=z^{K1}\underset{j=0}{\overset{\mathrm{}}{}}c_j\left(\frac{1}{z}\right)^j\underset{k=0}{\overset{\mathrm{}}{}}I_k\left(\frac{1}{z}\right)^k,$$

where $`c_j`$ are the coefficients of the expansion of $`_{n=0}^K\sqrt{(1a_n/z)(1b_n/z)}`$.

$$I_j=Q_++(1)^jQ_{}+_L𝑑xx^j\phi (x).$$

Equaling all powers $`z^{i1}`$ for $`i=0,\mathrm{},K`$ to zero, we obtain the triangular set of linear homogeneous equations for $`I_j`$:

$`c_0I_0=0`$
$`c_0I_1+c_1I_0=0`$
$`\mathrm{}`$
$`c_0I_K+\mathrm{}+c_KI_0=0,`$
which has the only solution (all $`c_i0`$) $`I_0=\mathrm{}=I_K=0`$, which are the Eqs.(25).

3. The normalization of $`\sigma (z)`$ from Eqs.(25)

We finally prove that the condition $`|dz|\sigma (z)=\nu `$ is contained within Eqs.(25).

Contour integration of $`E(z)`$ yields

$`{\displaystyle |dz|\sigma (z)}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \mathrm{i}}}{\displaystyle _{C_\mathrm{\Gamma }}}𝑑zE(z)+{\displaystyle _{L\mathrm{\Gamma }}}𝑑z\rho (z)`$
$`=`$ $`S(1/p)Q_+S(1/p)Q_{}`$
$`{\displaystyle _{L\mathrm{\Gamma }}}𝑑xS(x)\phi (x)+{\displaystyle _{L\mathrm{\Gamma }}}𝑑z\rho (z)`$
$`=`$ $`Q_+^{(0)}+Q_{}^{(0)}+1/2=\nu ,`$
where we exploited the fact that $`\text{Res}\{E(z),\mathrm{}\}=0`$, which is the equation in (25) for $`n=K`$ if all the equations for $`0n<K`$ are fulfilled.

## Appendix D Ground state and Excited states

We argue that one should not include arcs with end-points on $`L`$. This is because the charges on such an ansatz arc have to be mobile. The charges on the real axis however find themselves arrested in the “cells” flanked by two adjacent fixed charges. Therefore the only possibility for them to become mobile (in the thermodynamic limit) is to “escape” into the complex plane. This however can happen only for solution charges in neighbored “cells”RICH-NUM . Zones of separated “cells” $`L_{sep}`$ are characterized by $`\stackrel{~}{\rho }(x):=\sigma (x)\rho (x)0`$; $`xL_{sep}`$. We can contribute for them, just taking $`\stackrel{~}{\rho }`$ as the density of fixed charges. The complementary region $`L_{pairing}:=LL_{sep}`$ is thus characterized by $`\sigma (x)\rho (x)>0`$ at $`G=0`$ and we let $`\stackrel{~}{\rho }(x):=\rho (x)`$. In the simplest case we have $`\stackrel{~}{\rho }(x)=\rho (x)`$ everywhere and also the modulus of the total charge is $`\rho (x)`$ in each point. This means $`\sigma (x)=0`$ for $`xL_{sep}`$ and distinct connected regions of neighbored “cells” with $`\sigma (x)=2\rho (x)`$ for $`xL_{pairing}`$ at $`G=0`$ (see Fig.6). This most simple situation applies in particular to the ground state, which we have considered in the numerics. The general situation can be attacked replacing $`\rho (x)`$ with $`\stackrel{~}{\rho }(x)`$ and the filling $`\nu `$ by

$`\stackrel{~}{\nu }`$ $`:=`$ $`{\displaystyle _\mathrm{\Gamma }}|\mathrm{d}z|\stackrel{~}{\sigma }(z)\stackrel{G=0}{=}{\displaystyle _{L_{pairing}}}dx\stackrel{~}{\rho }(x)`$
$`=`$ $`\nu {\displaystyle _{L_{sep}}}dx\left[\stackrel{~}{\rho }(x)+\rho (x)\right].`$
That the above argument applies also to the generalized class of BCS models discussed here, follows from the same structure of the equation after the transformation. But it can equally be seen from the electrostatic analogy, since no two charges can penetrate each other due to the diverging forces when sitting upon each other, which must be zero instead due to the Richardson equations. These infinite forces can only be overcome if either at least one of the external charges diverges (i.e. the case G=0) or if two solution charges together approach a fixed charge (in the presence of a degeneracy $`d`$, more than two solution charges are needed because the charge ratio $`|q_{mobile}/q_{fixed}|=2/d`$ is diminished by $`d`$). We come back to this discussion after having presented the solution of the problem.

Now we explain why the solution of Eq.(20) corresponds to the ground state of the system, which implies that up to the crossing point of $`\mathrm{\Gamma }`$ with the Debye shell $`L`$ all “cells” between fixed charges are occupied. This implies that there we have $`\sigma =2\rho `$ and hence the total charge density is $`\rho `$. This sign change of the charge density below the crossing point is already included in Eq.(20), since we apparently replaced

$$\frac{\rho (x)}{S(x)}\frac{\rho (x)}{\sqrt{(x\lambda )^2+\mathrm{\Delta }^2}}.$$
(63)
Notice that the complex square root changes sign at the crossing point $`x_c`$: indeed we have

$$S(x)=\{\begin{array}{ccc}\sqrt{(x\lambda )^2+\mathrm{\Delta }^2}\hfill & & x>x_c\hfill \\ & \text{for}& \\ \sqrt{(x\lambda )^2+\mathrm{\Delta }^2}\hfill & & x<x_c\hfill \end{array}.$$

This means that in Eq.(20) the charge density is $`\rho (x)`$ for $`x<x_c`$ and $`\rho (x)`$ for $`x>x_c`$, which correspond to the ground state. For general excited states, which might also include more than one arc, we must furnish the real square roots with the proper signs as indicated in Fig.6. As long as all the arcs cross the real axis in the same interval of $`L_{pairing}`$ where they are generated from, this procedure can be implemented in the formalism described in this work. Otherwise the crossing points could be determined numerically.