format: text

text:
# Large 𝑁 limit of the exactly solvable BCS model: analytics versus numerics

## I) Introduction

One of the most studied models in Condensed Matter Physics is the BCS model of superconductivity, proposed by Bardeen, Cooper, and Schrieffer in 1957 BCS . The BCS model is usually formulated in the grand canonical ensemble, which is appropiated to systems with a macroscopic number of fermions. However for small systems, such as nuclei or ultrasmall metallic grains, one has to consider the canonical ensemble where the BCS wave function is not adequate.

Fortunately enough, there is a simplified version of the BCS model which is exactly solvable in the canonical ensemble. This occurs when the strengh of the pairing interaction between all the energy levels is the same. The exact solution of the reduced BCS model was obtained by Richardson in 1963 R1 ; R1bis ; RS , who also studied its consequences and properties in a series of papers in the 60’s and 70’s R-norm ; R-roots ; R-boson ; R-limit . The exactly solvable BCS model, in turn, is closely related to the rational spin model of Gaudin, defined in terms of a set of commuting Hamiltonians Gaudin ; G-book .

The integrability of the reduced BCS Hamiltonian was proved in 1997 by Cambiaggio, Rivas, and Saraceno (CRS) CRS , who constructed a set of conserved quantities in involution, commuting with the BCS Hamiltonian and closely related to the rational Gaudin’s Hamiltonians Gaudin . More recently Gaudin’s trigonometric model have been generalized, to include a $`g`$-term à la BCS, in references ALO ; DES . There are also bosonic pairing Hamiltonians satisfying the same type of equations as in the fermionic case R-boson ; DS ; DP . Generalizations of the $`SU(2)`$ BCS model to other Lie groups AFS  and supergroups KM  have been worked out. Both, the BCS and Gaudin’s models are intimately linked to conformal field theory B ; BF ; S , Chern-Simons theory AFS , and Integrable Models in Statistical Mechanics AFF ; ZLMG ; vDP  (see ref. \[Scomo, \] for a short review on these topics). The exact solution has also been used to study the effect of level statistics random  and finite temperature Finite-T ; FFM  on the ultrasmall metallic grains.

An important property of the exact solution is that the energy of the states and the occupation numbers agree, to leading order in the number of particles, with the BCS theory BCS . This result was obtained by Gaudin G-book  and Richardson R-limit , using an electrostatic analogy of the Richardson’s equations, which have to be solved for finding the eigenstates of the BCS Hamiltonian. More recent applications of this electrostatic model to nuclear pairing can be found in ref. \[DEP, \].

In the limit of large number of levels $`N`$, Richardson’s equations become integral equations, which can be solved with techniques of complex analysis G-book . This is sufficient to determine the location and density of roots, as well as the total energy of the ground state to leading order in $`N`$. To study higher order corrections of the ground state energy and the excitations, Richardson developed a multipole expansion of the fields appearing in the electrostatic analogue model R-limit . We shall follow Gaudin’s approach which is more geometrical, but use also some of the results obtained by Richardson.

The aim of this paper is to make a revision and precise comparison between the numerical solution and the analytic expressions for the location and density of roots of the Richardson’s equations, as proposed originally by Gaudin in his unpublished paper in 1968, which got finally published in his collected papers in 1995 G-book . The comparison will be made in the limit of large $`N`$ for the two-level model and the equally spaced model.

These models have been extensively used in Nuclear Physics since they were introduced by Hogaasen-Feldman Ho  (two-level) and Richardson himself R-roots  (equally spaced). The latter model has also been used to describe the physics of ultrasmall metallic grains grains ; Lanczos ; BvD ; DMRG ; Sch  (for a review see ref. \[vDR, \]). In the large $`N`$ limit the equally spaced model is superconducting for all values of atractive BCS couplings $`g`$, while the two-level model displays a quantum phase transition between a superconducting state and a normal state as a function of the BCS coupling constant.

The organization of the paper is as follows: in sections II and III we give a primer of the Richardson’s exact solution of the BCS model, together with its integrability. In section IV we introduce the electrostatic analogue model of BCS. In section V we review Gaudin’s continuum limit of the Richardson’s equations and their general solution, using complex analysis. In section VI we derive the gap and chemical potential BCS equations under the assumption that the roots of Richardson equations collapse into a single arc, and study in detail the two-level and equally spaced models. In section VII we present our numerical results and compare them with the analytic ones obtained in section VI. We also study the scaling properties of the roots and the ground state energy. In section VIII we present our conclusions and prospects. Finally, we describe in an appendix the conformal mappings associated to the equally spaced model.

## II) Richardson’s exact solution of the reduced BCS model

The reduced BCS model is defined by the Hamiltonian BCS ; vDR ; grains

$$H_{BCS}=\frac{1}{2}\underset{j,\sigma =\pm }{}\epsilon _{j\sigma }c_{j\sigma }^{}c_{j\sigma }G\underset{j,j^{}}{}c_{j+}^{}c_j^{}c_j^{}c_{j^{}+},$$
(1)
where $`c_{j,\pm }`$ (resp. $`c_{j,\pm }^{}`$) is an electron annihilation (resp. creation) operator in the time-reversed states $`|j,\pm `$ with energies $`\epsilon _j/2`$, and $`G`$ is the BCS dimensionful coupling constant. The sums in (1) run over a set of $`N`$ doubly degenerate energy levels $`\epsilon _j/2(j=1,\mathrm{},N)`$. We adopt Gaudin’s notation, according to which $`\epsilon _j`$ denotes the energy of a pair occupying the level $`j`$ G-book . We shall assume that the energy levels are all distinct, i.e. $`\epsilon _i\epsilon _j`$ for $`ij`$. The Hamiltonian (1) is a simplified version of the reduced BCS Hamiltonian where all couplings have been set equal to a single one, namely $`G`$. Hereafter we shall refer to (1) simply as the BCS Hamiltonian.

Richardson had long ago solved this model exactly for an arbitrary set of levels, $`\epsilon _j`$ R1 ; R1bis ; RS . To simplify matters, we shall assume that there are not singly occupied electronic levels. As can be seen from (1), these levels decouple from the rest of the system; they are said to be blocked, contributing only with their energy $`\epsilon _j/2`$ to the total energy $`E`$. The above simplification implies that every energy level $`j`$ is either empty (i.e. $`|\mathrm{vac}`$), or occupied by a pair of electrons (i.e. $`c_{j,+}^{}c_{j,}^{}|\mathrm{vac}`$). Denote the total number of electrons pairs by $`M`$. Then of course $`MN`$. The most studied case in the literature corresponds to the half-filled situation, where the number of electrons, $`N_e=2M`$, equals the number of levels $`N`$ vDR . In the absence of interaction (i.e. $`G=0`$), all the pairs occupy the lowest energy levels, forming a Fermi sea. The pairing interaction promotes the pairs to higher energies and even for small values of $`M`$ or $`G`$ the levels are pair correlated Lanczos ; BvD ; DMRG ; vDR .

In order to describe Richardson’s solution one defines the hard-core boson operators

$$b_j=c_{j,}c_{j,+},b_j^{}=c_{j,+}^{}c_{j,}^{},N_j=b_j^{}b_j,$$
(2)
which satisfy the commutation relations,

$$[b_j,b_j^{}^{}]=\delta _{j,j^{}}(12N_j).$$
(3)
The Hamiltonian (1) can then be written as

$$H_{BCS}=\underset{j}{}\epsilon _jb_j^{}b_jG\underset{j,j^{}}{}b_j^{}b_j^{}.$$
(4)
Richardson showed that the eigenstates of this Hamiltonian with $`M`$ pairs have the (unnormalized) product form R1 ; R1bis ; RS

$$|M=\underset{\nu =1}{\overset{M}{}}B_\nu |\mathrm{vac},B_\nu =\underset{j=1}{\overset{N}{}}\frac{1}{\epsilon _jE_\nu }b_j^{},$$
(5)
where the parameters $`E_\nu `$ ($`\nu =1,\mathrm{},M`$) are, in general, complex solutions of the $`M`$ coupled algebraic equations

$$\frac{1}{G}=\underset{j=1}{\overset{N}{}}\frac{1}{\epsilon _jE_\nu }\underset{\mu =1(\nu )}{\overset{M}{}}\frac{2}{E_\mu E_\nu },\nu =1,\mathrm{},M$$
(6)
which play the role of Bethe ansatz equations for this problem B ; BF ; S ; AFF ; ZLMG ; vDP ; Scomo . The energy of these states is given by the sum of the auxiliary parameters $`E_\nu `$, i.e.

$$E(M)=\underset{\nu =1}{\overset{M}{}}E_\nu .$$
(7)
The ground state of $`H_{BCS}`$ is given by the solution of eqs. (6) which gives the lowest value of $`E(M)`$. The (normalized) states (5) can also be written as R-norm

$$|M=\frac{C}{\sqrt{M!}}\underset{j_1,\mathrm{},j_M}{}\psi (j_1,\mathrm{},j_M)b_{j_1}^{}\mathrm{}b_{j_M}^{}|\mathrm{vac},$$
(8)
where the sum excludes double occupancy of pair states and the wave function $`\psi `$ takes the form

$$\psi (j_1,\mathrm{},j_M)=\underset{𝒫}{}\underset{k=1}{\overset{M}{}}\frac{1}{\epsilon _{j_k}E_{𝒫k}}.$$
(9)
The sum in (9) runs over all the permutations, $`𝒫`$, of $`1,\mathrm{},M`$. The constant $`C`$ in (8) guarantees the normalization of the state R-norm  (i.e. $`M|M=1`$).

## III) Integrability of the reduced BCS Hamiltonian

A well known fact about the BCS Hamiltonian is that it is equivalent to that of a XY spin model with long range couplings, and a “position dependent” magnetic field proportional to $`\epsilon _j`$. To see this, let us represent the hard-core boson operators (2) in terms of the Pauli matrices as follows,

$$b_j=\sigma _j^+,b_j^{}=\sigma _j^{},N_j=\frac{1}{2}(1\sigma _j^z),$$
(10)
in which case the Hamiltonian (4) becomes

$`H_{BCS}`$ $`=`$ $`H_{XY}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}\epsilon _j+G(N/2M),`$ (11)
$`H_{XY}`$ $`=`$ $`{\displaystyle \underset{j}{}}\epsilon _jt_j^0{\displaystyle \frac{G}{2}}(T^+T^{}+T^{}T^+),`$
where the matrices

$$T^a=\underset{j=1}{\overset{N}{}}t_j^a,(a=0,+,),\text{with}t_j^0=\frac{1}{2}\sigma _j^z,t_j^+=\sigma _j^+,t_j^{}=\sigma _j^{},$$
(12)
satisfy the $`SU(2)`$ algebra,

$$[T^a,T^b]=f_c^{ab}T^c,\text{with}f_+^{+0}=f_{}^0=1,f_0^+=2,$$
(13)
whose Casimir is given by

$$𝐓𝐓=T^0T^0+\frac{1}{2}(T^+T^{}+T^{}T^+).$$
(14)
This spin representation of the BCS model is the appropiate one to study its integrability, as it was shown by Cambiaggio, Rivas, and Saraceno (CRS), who constructed a set of operators, $`R_j(j=1,\mathrm{},N)`$ CRS ,

$$R_i=t_i^02g\underset{j(i)}{\overset{N}{}}\frac{𝐭_i𝐭_j}{\epsilon _i\epsilon _j},(i=1,\mathrm{},\mathrm{\Omega }),$$
(15)
satisfying the properties

$$[H_{BCS},R_i]=[R_i,R_j]=0,(i,j=1,\mathrm{},N).$$
(16)
Let us denote by $`\lambda _j`$ the eigenvalue of $`R_j`$ acting on the state (8), i.e.,

$$R_j|M_R=\lambda _j|M_R.$$
(17)
CRS, unaware of Richardson’s solution, did not give an expression of $`\lambda _j`$ in their work. This was found in reference \[S, \] using CFT techniques,

$$\lambda _i=\frac{1}{2}+G\left(\underset{\nu =1}{\overset{M}{}}\frac{1}{\epsilon _iE_\nu }\frac{1}{2}\underset{j=1(i)}{\overset{N}{}}\frac{1}{\epsilon _i\epsilon _j}\right),$$
(18)
which verifies $`_i\lambda _i=(MN/2)`$ upon using eqs. (6), so that it vanishes at half filling.

However CRS showed that $`H_{XY}`$ given in eq. (11) can be expressed in terms of the operators $`R_i`$ as

$$H_{XY}=\underset{j}{}\epsilon _jR_j+G\left(\underset{j}{}R_j\right)^2\frac{3}{4}GN.$$
(19)
Then using the eqs. (11), (18), and (6) one can show that $`E(M)`$ is indeed given by eq. (7). (For the two-level system the last term in (19) must be replaced by $`2G(N/4+1)N/4`$ ).

## IV) The Gaudin-Richardson’s electrostatic model of BCS

The equations (6) admit a 2 dimensional electrostatic analogy, which will be very useful in the study of the limit $`N\mathrm{}`$ G-book ; R-limit . Let us consider a set of $`N`$ charges $`1/2`$ fixed at the positions $`\epsilon _j`$ on the real axis, and a uniform field parallel to this axis with stregth $`1/2G`$. The problem is to find the equilibrium positions of $`M`$ charges $`+1`$ at positions $`E_\mu `$ subject to their mutual repulsion, the attraction with the $`1/2`$ charges, and the action of the uniform field. The 2D electrostatic potential is given by $`W+W^{}`$, where $`W`$ reads,

$`W(\{E_\mu \},\{\epsilon _i\})`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,\mu }{}}\mathrm{log}(\epsilon _iE_\mu ){\displaystyle \underset{\mu >\nu }{}}\mathrm{log}(E_\nu E_\mu )`$ (20)
$`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i>j}{}}\mathrm{log}(\epsilon _i\epsilon _j)+{\displaystyle \frac{1}{2G}}\left({\displaystyle \underset{\mu }{}}E_\mu {\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\epsilon _i\right).`$
It is easy to see that the Richardson’s eqs. (6) arise as the stationary conditions $`W/E_\mu =0`$, while the conserved quantities (18) are proportional to the forces exerted on the fixed charges, i.e.

$$\lambda _i=2G\frac{W}{\epsilon _i}.$$
(21)
The total energy (7) is the center of gravity of the charges $`E_\mu `$, which must be located symmetrically with respect to the real axis, i.e. if $`E_\mu `$ is a solution to the eqs. (6), then $`E_\mu ^{}`$ must also be a solution. This condition is fullfilled in two cases: i) either $`E_\mu `$ is real, or ii) $`E_\mu `$ and $`E_\mu ^{}`$ form a complex conjugate pair. The formation of these complex pairs is usually an indication of the superconducting properties of the ground state.

## V) Continuum limit of Richardson’s equations

In reference \[G-book, \] Gaudin proposed a continuum limit of the equations (6), in order to obtain the BCS formulas for the ground state energy and occupation number of levels. We shall follow closely this reference. This continuum limit is defined by taking the number of levels $`N`$ going to infinity, while keeping fixed the following quantities:

$$g=GN,M/N.$$
(22)
At the same time the pair energy levels $`\epsilon _i`$ will be equivalent to a negative charge density $`\rho (\epsilon )`$ located on an interval $`\mathrm{\Omega }`$ of the real axis. The total charge of this interval is given by

$$_\mathrm{\Omega }\rho (\epsilon )𝑑\epsilon =\frac{N}{2}.$$
(23)
We shall suppose that $`\mathrm{\Omega }`$ is contained inside an interval $`(\omega ,\omega )`$. For example, in the BCS model with equally spaced levels, $`\omega /2`$ is equal to the Debye energy.

The basic asumption made by Gaudin, which is supported by numerical results, is that the solutions of eqs. (6) organize themselves into arcs $`\mathrm{\Gamma }_k(k=1,\mathrm{},K)`$, which are piece-wise differentiable and symmetric under reflection on the real axis. We shall call $`\mathrm{\Gamma }`$ the union of all these arcs, and $`r(\xi )`$ the linear charge density of roots $`E_\mu `$ in the complex plane. Hence, the total number of pairs, $`M`$, and the total energy, $`E`$, are given by

$`{\displaystyle _\mathrm{\Gamma }}r(\xi )|d\xi |`$ $`=`$ $`M,`$ (24)
$`{\displaystyle _\mathrm{\Gamma }}\xi r(\xi )|d\xi |`$ $`=`$ $`E.`$ (25)
The continuum limit of eqs. (6) is

$$_\mathrm{\Omega }\frac{\rho (\epsilon )d\epsilon }{\epsilon \xi }P_\mathrm{\Gamma }\frac{r(\xi ^{})|d\xi ^{}|}{\xi ^{}\xi }\frac{1}{2G}=0,\xi \mathrm{\Gamma },$$
(26)
which implies that the total electric field on every point of the arcs $`\mathrm{\Gamma }_k`$ is null. On the other hand, the values of the conserved quantities (18) are given in the continuum by

$$\lambda (\epsilon )=\frac{1}{2}+G\left(_\mathrm{\Gamma }\frac{r(\xi )|d\xi |}{\epsilon \xi }P_\mathrm{\Omega }\frac{\rho (\epsilon ^{})d\epsilon ^{}}{\epsilon \epsilon ^{}}\right),\epsilon \mathrm{\Omega }.$$
(27)
The formal solution of eqs. (26) can be found as follows. First of all, let us orient each arc $`\mathrm{\Gamma }_k`$ from the point $`a_k`$ to the point $`b_k`$, and call $`L_k`$ an anticlockwise path encircling $`\mathrm{\Gamma }_k`$. We look for an analytic field $`h(\xi )`$ outside $`\mathrm{\Gamma }`$ and the set $`\mathrm{\Omega }`$, such that

$$r(\xi )|d\xi |=\frac{1}{2\pi i}(h_+(\xi )h_{}(\xi ))d\xi ,\xi \mathrm{\Gamma },$$
(28)
where $`h_+(\xi )`$ and $`h_{}(\xi )`$ denote the limit values of $`h(\xi )`$ to the right and left of $`\mathrm{\Gamma }`$. This can be understood using the electrostatic equivalence, considering that the electric field presents a discontinuity proportional to the superficial density of charge when crossing such surface. Next, we define a function $`R(\xi )`$, with cuts along the curves $`\mathrm{\Gamma }_k`$, by the equation

$$R(\xi )=\left[\underset{k=1}{\overset{K}{}}(\xi a_k)(\xi b_k)\right]^{1/2},$$
(29)
and look for a solution which vanishes at the boundary points of $`\mathrm{\Gamma }`$, in the form

$$h(\xi )=R(\xi )_\mathrm{\Omega }\frac{\phi (\epsilon )d\epsilon }{\epsilon \xi }.$$
(30)
This field has to be constant at infinity, hence the first $`K2`$ moments of $`\phi `$ must vanish, i.e.

$$_\mathrm{\Omega }\epsilon ^{\mathrm{}}\phi (\epsilon )=0,0\mathrm{}<K1.$$
(31)
This formula is empty if $`K=1`$. The contour integral surrounding the charge density $`r(\xi )`$ in (26) can be expressed as

$$_\mathrm{\Gamma }\frac{r(\xi ^{})|d\xi ^{}|}{\xi \xi ^{}}=\frac{1}{2\pi i}_\mathrm{\Gamma }\frac{h_+(\xi ^{})h_{}(\xi ^{})}{\xi \xi ^{}}𝑑\xi ^{}=\frac{1}{2\pi i}_L\frac{h(\xi ^{})}{\xi \xi ^{}}𝑑\xi ^{},\xi 𝐂\mathrm{\Gamma },$$
(32)
where $`𝐂\mathrm{\Gamma }`$ is the region outside the curves in $`\mathrm{\Gamma }`$. Using eqs. (29) and (30), one finds for the principal value of (32)

$$P_\mathrm{\Gamma }\frac{r(\xi ^{})|d\xi ^{}|}{\xi \xi ^{}}=\frac{1}{2\pi i}_L\frac{d\xi ^{}R(\xi ^{})}{\xi \xi ^{}}_\mathrm{\Omega }\frac{\phi (\epsilon )d\epsilon }{\epsilon \xi ^{}}=_\mathrm{\Omega }\frac{\phi (\epsilon )R(\epsilon )}{\epsilon \xi }𝑑\epsilon +_\mathrm{\Omega }\epsilon ^{K1}\phi (\epsilon )𝑑\epsilon ,\xi \mathrm{\Gamma },$$
(33)
where we have deformed the contour of integration $`L`$ into two contours, one encircling the interval $`\mathrm{\Omega }`$ (first term) and another one around the infinity (second term). We are assuming that $`\mathrm{\Gamma }`$, and consequently $`L`$, do not cut the interval $`\mathrm{\Omega }`$. The case where $`L`$ intersects $`\mathrm{\Omega }`$ actually arises in the equally spaced model, and will be discussed later on.

Plugging eq. (33) into (26), we see that a solution is obtained provided the following equations are satisfied:

$`\phi (\epsilon )={\displaystyle \frac{\rho (\epsilon )}{R(\epsilon )}},`$ (34)
$`{\displaystyle _\mathrm{\Omega }}\epsilon ^{K1}{\displaystyle \frac{\rho (\epsilon )}{R(\epsilon )}}𝑑\epsilon ={\displaystyle \frac{1}{2G}},`$ (35)
which have to be suplemented by eq. (31), which can be rewritten as

$$_\mathrm{\Omega }\epsilon ^{\mathrm{}}\frac{\rho (\epsilon )}{R(\epsilon )}𝑑\epsilon =0,0\mathrm{}<K1.$$
(36)
The field $`h(\xi )`$ gives also the density of charges $`r(\xi )`$,

$$r(\xi )=\frac{1}{\pi }|h(\xi )|,\xi \mathrm{\Gamma },$$
(37)
and its value is given by replacing (34) into (30), i.e.

$$h(\xi )=R(\xi )_\mathrm{\Omega }\frac{\rho (\epsilon )}{R(\epsilon )}\frac{d\epsilon }{\epsilon \xi },$$
(38)
whose value at infinity is $`1/2G`$. Moreover, using eq.(32) and performing the same contour deformations that lead to eq.(33), one can show that the conserved quantities (27) are given by

$$\lambda (\epsilon )=P\left[GR(\epsilon )_\mathrm{\Omega }\frac{\rho (\epsilon ^{})}{R(\epsilon ^{})}\frac{d\epsilon ^{}}{\epsilon ^{}\epsilon }\right]=\frac{G}{2}\left(h(\epsilon +i\mathrm{\hspace{0.33em}0})+h(\epsilon i\mathrm{\hspace{0.33em}0})\right),\epsilon \mathrm{\Omega },$$
(39)
so that $`\lambda (\epsilon )`$ is the principal part of $`h(\epsilon )`$ on the set $`\mathrm{\Omega }`$. Finally the equations fixing the arcs $`\mathrm{\Gamma }_k`$ are the equipotential curves of the total distribution, i.e.

$$_{a_k}^\xi h(\xi ^{})𝑑\xi ^{}=0,\xi \mathrm{\Gamma }_k.$$
(40)
## VI) The BCS equations

In this section we shall show how the BCS equations describing the ground state of the model can be derived from the formalism of section V G-book .

The basic assumption is that for the ground state all the roots $`E_\mu `$ form a single arc, i.e. $`K=1`$, leaving only two complex parameters $`aa_1`$ and $`bb_1`$, which shall be denoted as

$$a=\epsilon _0i\mathrm{\Delta },b=\epsilon _0+i\mathrm{\Delta },$$
(41)
where $`\epsilon _0`$ and $`\mathrm{\Delta }>0`$ are real parameters.

From eq. (38) the electrostatic field is given by

$$h(\xi )=[(\xi \epsilon _0)^2+\mathrm{\Delta }^2]^{\frac{1}{2}}_\mathrm{\Omega }\frac{\rho (\epsilon )}{[(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2]^{\frac{1}{2}}}\frac{d\epsilon }{\epsilon \xi }.$$
(42)
And the conserved quantities $`\lambda (\epsilon )`$ can be derived from the eq. (39). Eq. (35) yields the BCS gap equation

$$_\mathrm{\Omega }\frac{\rho (\epsilon )d\epsilon }{\sqrt{(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2}}=\frac{1}{2G}.$$
(43)
Similarly, eq. (24) becomes the chemical potential equation

$$M=\frac{1}{2\pi i}_Lh(\xi )𝑑\xi =_\mathrm{\Omega }\left(1\frac{\epsilon \epsilon _0}{\sqrt{(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2}}\right)\rho (\epsilon )𝑑\epsilon ,$$
(44)
while eq. (25) gives the BCS expression for the ground state energy,

$$E=\frac{1}{2\pi i}_L\xi h(\xi )𝑑\xi =\frac{\mathrm{\Delta }^2}{4G}+_\mathrm{\Omega }\left(1\frac{\epsilon \epsilon _0}{\sqrt{(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2}}\right)\rho (\epsilon )\epsilon 𝑑\epsilon .$$
(45)
Gaudin’s paper G-book  contains a misprint in this equation since the term $`\frac{\mathrm{\Delta }^2}{4G}`$ is quoted as $`\frac{\mathrm{\Delta }^2}{2G}`$.

Comparing these equations with the corresponding ones in the BCS theory, we deduce the following relations between $`\mathrm{\Delta }`$, $`\epsilon _0`$, and $`\rho (\epsilon )`$, and $`\mathrm{\Delta }_{BCS}`$ (BCS gap), $`\mu `$ (chemical potential), and $`n(ϵ)`$ (single particle energy density):

$$\mathrm{\Delta }=2\mathrm{\Delta }_{BCS},\epsilon _0=2\mu ,\rho (\epsilon )=\frac{1}{4}n\left(\frac{\epsilon }{2}\right).$$
(46)
The factor $`1/4`$, in the last term of (46), emerges from the normalization $`1/2`$ of the charge density $`\rho (\epsilon )`$, times another factor $`1/2`$, due to the fact that the separation between energy pairs is twice the separation between single particle energy levels.

To determine the equation of the curve $`\mathrm{\Gamma }`$ we define the complex variable

$$z=\sqrt{(\xi \epsilon _0)^2+\mathrm{\Delta }^2}=x+iy.$$
(47)
Using eqs. (40) and (42) we deduce for $`\mathrm{\Gamma }`$,

$$0=_\mathrm{\Omega }𝑑\epsilon \rho (\epsilon )\left\{\frac{x}{\sqrt{(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2}}+\frac{1}{4}\mathrm{log}\frac{\left(x\sqrt{(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2}\right)^2+y^2}{\left(x+\sqrt{(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2}\right)^2+y^2}\right\}.$$
(48)
Eqs. (43), (44), and (45) have been derived under the condition that $`\mathrm{\Gamma }`$ does not cut the set $`\mathrm{\Omega }`$, which requires that (48) has no solutions with $`\xi \mathrm{\Omega }`$. In the equally spaced model the curve $`\mathrm{\Gamma }`$ may intersect $`\mathrm{\Omega }`$, so this case has to be treated with some care (see below).

Finally, the mean occupation of the energy level $`\epsilon _j`$ can be computed from the formula $`N_j=E/\epsilon _j`$. In the large $`N`$ limit $`N_j`$ can be shown to be given by G-book ; R-limit

$$N(\epsilon )=\frac{1}{2}\left(1\frac{\epsilon \epsilon _0}{\sqrt{(\epsilon \epsilon _0)^2+\mathrm{\Delta }^2}}\right),\epsilon \mathrm{\Omega },$$
(49)
which agrees with the BCS result BCS . Next we shall apply the previous formalism to two models.

### The two-level model

Let us suppose that we have two energy levels $`\pm \epsilon _1`$, which can accomodate at most $`N/2`$ pairs each. The charge density is given by

$$\rho (\epsilon )=\frac{N}{4}[\delta (\epsilon +\epsilon _1)+\delta (\epsilon \epsilon _1)].$$
(50)
From eq. (44) we deduce that the chemical potential vanishes, i.e. $`\epsilon _0=0`$, while the gap eq. (43) gives

$$\epsilon _1^2+\mathrm{\Delta }^2=g^2,g=GN.$$
(51)
This equation implies that $`g`$ must be greater than $`\epsilon _1`$. The total energy of this solution follows from (45),

$$E=\frac{N}{4g}(\epsilon _1^2+g^2).$$
(52)
Integrating eq. (42) we obtain the electric field, and accordingly the value of the conserved quantities

$$h(\xi )=\frac{N\xi }{2g}\frac{\sqrt{\xi ^2+\mathrm{\Delta }^2}}{\epsilon _1^2\xi ^2}\lambda (\epsilon _1)=\frac{1}{8g\epsilon _1}(2\epsilon _1^2+g^2).$$
(53)
Similarly one finds $`\lambda (\epsilon _1)=\lambda (\epsilon _1)`$, in agreement with the condition $`_i\lambda _i=0`$ derived in sect. III. The equation of the curve (48) becomes in this case

$$\frac{4x}{g}+\mathrm{log}\frac{(xg)^2+y^2}{(x+g)^2+y^2}=0,$$
(54)
or equivalently

$$x^2+y^2+g^2=\frac{2xg}{\mathrm{th}(2x/g)},$$
(55)
where $`x`$ and $`y`$ are given by eq. (47). In Gaudin’s paper G-book  there is a misprint in eqs. (52) and (55).

From the numerical results shown in fig. 1, when $`g>\epsilon _1`$ $`\mathrm{\Gamma }`$ is an open arc between the imaginary points $`\pm i\mathrm{\Delta }`$ in the complex $`\xi `$-plane. However when $`g=\epsilon _1`$, in the limit when $`N\mathrm{}`$, the arc closes at the origin, and for $`g<\epsilon _1`$ it remains closed surrounding the point $`\epsilon _1`$.

To treat this case Gaudin considers a general situation where $`\mathrm{\Gamma }`$ is a closed curve surrounding a piece $`\omega `$ of $`\mathrm{\Omega }`$ and defines a function $`s(\xi )`$ by the equation \[G-book, \]

$$s(\xi )d\xi =r(\xi )|d\xi |.$$
(56)
Then eq. (26) is solved by

$$s(\xi )=\frac{i}{\pi }\left[\frac{1}{2G}+_\omega \frac{\rho (\epsilon )d\epsilon }{\epsilon \xi }_{𝐂\omega }\frac{\rho (\epsilon )d\epsilon }{\epsilon \xi }\right],$$
(57)
which can be verified using $`P_\mathrm{\Gamma }𝑑\xi ^{}/(\xi ^{}\xi )=i\pi `$, when $`\xi \mathrm{\Gamma }`$. For the two-level model $`s(\xi )`$ is given by

$$s(\xi )=\frac{i}{\pi }\left[\frac{1}{2G}\frac{N}{2}\frac{\epsilon _1}{\epsilon _1^2\xi ^2}\right],\xi \mathrm{\Gamma }.$$
(58)
and eqs. (24) and (25) yield

$`M`$ $`=`$ $`{\displaystyle _\mathrm{\Gamma }}𝑑\xi s(\xi )={\displaystyle _\omega }𝑑\epsilon \mathrm{\hspace{0.33em}2}\rho (\epsilon )={\displaystyle \frac{N}{2}},`$ (59)
$`E`$ $`=`$ $`{\displaystyle _\mathrm{\Gamma }}𝑑\xi \xi s(\xi )={\displaystyle _\omega }𝑑\epsilon \mathrm{\hspace{0.33em}2}\epsilon \rho (\epsilon )={\displaystyle \frac{N}{2}}\epsilon _1,`$ (60)
Hence the state for $`g<\epsilon _1`$ is a Fermi sea to leading order in $`N`$.

The conserved quantities are derived from an equation similar to (39), where $`h(\xi )`$ is replaced by $`s(\xi )`$, i.e.

$$\lambda (\epsilon )=i\frac{\pi G}{2}\left(s(\epsilon +i0)+s(\epsilon i0)\right),\epsilon \omega \lambda (\epsilon _1)=\frac{1}{2}\left(1\frac{g}{4\epsilon _1}\right),$$
(61)
verifing, as in the previous case, the symmetry $`\lambda (\epsilon _1)=\lambda (\epsilon _1)`$.

To find the curve $`\mathrm{\Gamma }`$ we integrate (58), and impose the result to be real accordingly to (56). Impossing the imaginary part of eq. (56) to vanish we obtain a diferential equation for the curve $`\mathrm{\Gamma }`$, which is a check of consistency.

$$\left(\frac{\xi x_0}{g}+\frac{1}{2}\mathrm{log}\frac{\xi \epsilon _1}{\xi +\epsilon _1}\right)=0,\xi \mathrm{\Gamma },$$
(62)
where $`x_0`$ is an integration constant, which has a non vanishing value, contrary to Gaudin’s assumption. Eq. (62) can be written as

$$x^2+y^2+\epsilon _1^2=\frac{2x\epsilon _1}{\mathrm{th}(2(xx_0)/g)},\xi =x+iy.$$
(63)
To find $`x_0`$ we have first to look for the point $`\epsilon _0`$ where $`\mathrm{\Gamma }`$ intersects the real axis. From fig. 1, the density of roots at this point vanishes, i.e. $`s(\epsilon _0)=0`$, which implies

$$\epsilon _0=\epsilon _1\sqrt{1g/\epsilon _1}.$$
(64)
The value of $`x_0`$ is found by imposing that $`\epsilon _0`$ belongs to the curve (63), namely

$`x_0=\epsilon _0{\displaystyle \frac{g}{2}}\mathrm{log}{\displaystyle \frac{\epsilon _1+\epsilon _0}{\epsilon _1\epsilon _0}}.`$ (65)
As $`g`$ approaches $`\epsilon _1`$ from below, both $`\epsilon _0`$ and $`x_0`$ go to zero, and the curve (63) coincides with the curve (55) at the “critical point” $`g=\epsilon _1`$.

### Equally spaced model

Let us consider now a model with $`N`$ pair energy levels that are uniformly distributed in the interval $`[\omega ,\omega ]`$, with a density $`\rho (\epsilon )=\rho _0`$

$$\rho _0=\frac{N}{4\omega }.$$
(66)
The distance between the single particle energy levels, $`d`$, is given by $`d=\omega /N=2\omega _D/N`$, where $`\omega _D`$ is the Debye energy. This is the model used in most of the studies on ultrasmall superconducting grains vDR . At half-filling the number of pairs is $`M=N/2`$, which implies that $`\epsilon _0=0`$ (see eq. (44)). On the other hand, the gap eq. (43) yields the standard result BCS

$$\mathrm{\Delta }=\frac{\omega }{\mathrm{sinh}(1/g)},g=\frac{GN}{\omega },$$
(67)
where $`g`$ is the usual dimensionless BCS coupling constant. Recall that $`\mathrm{\Delta }`$ is twice the BCS gap $`\mathrm{\Delta }_{BCS}`$. The total energy of the BCS ground state can be derived from (45),

$$E=\frac{\omega ^2}{4d}\sqrt{1+\left(\frac{\mathrm{\Delta }}{\omega }\right)^2},$$
(68)
while the energy of the Fermi sea state is given by

$$E_{FS}=\frac{\omega ^2}{4d}\frac{GN}{2}.$$
(69)
The condensation energy is defined as $`E_C=EE_{FS}`$, which in the weak coupling limit, where $`\mathrm{\Delta }\omega `$, behaves as

$$E_C\frac{\mathrm{\Delta }_{BCS}^2}{2d},$$
(70)
the well known result.

Next, we shall study the shape of the arcs in this model. The electrostatic field $`h(\xi )`$ is given by

$$h(\xi )=\rho _0\mathrm{log}\frac{(\xi \omega )(\mathrm{\Delta }^2\xi \omega +\sqrt{(\mathrm{\Delta }^2+\xi ^2)(\mathrm{\Delta }^2+\omega ^2)}}{(\xi +\omega )(\mathrm{\Delta }^2+\xi \omega +\sqrt{(\mathrm{\Delta }^2+\xi ^2)(\mathrm{\Delta }^2+\omega ^2)}}.$$
(71)
To simplify matters, Gaudin considered the limit $`\mathrm{\Delta }\omega `$, where eq. (71) becomes

$$h(\xi )=\rho _0\mathrm{log}\frac{\xi \sqrt{\mathrm{\Delta }^2+\xi ^2}}{\xi +\sqrt{\mathrm{\Delta }^2+\xi ^2}},|\xi |\omega .$$
(72)
The conserved quantities $`\lambda (\epsilon )`$ are given by

$$\lambda (\epsilon )=\frac{g}{4}\mathrm{log}\frac{|(\epsilon \omega )(\mathrm{\Delta }^2\epsilon \omega +\sqrt{(\mathrm{\Delta }^2+\epsilon ^2)(\mathrm{\Delta }^2+\omega ^2)}|}{|(\epsilon +\omega )(\mathrm{\Delta }^2+\epsilon \omega +\sqrt{(\mathrm{\Delta }^2+\epsilon ^2)(\mathrm{\Delta }^2+\omega ^2)}|}\stackrel{|\epsilon |\omega }{}\frac{g}{4}\mathrm{log}\frac{|\epsilon \sqrt{\mathrm{\Delta }^2+\epsilon ^2}|}{|\epsilon +\sqrt{\mathrm{\Delta }^2+\epsilon ^2}|}.$$
(73)
In the low energy approximation the equation of the equipotential curves passing through the points $`\epsilon _0\pm i\mathrm{\Delta }`$ is given by

$$\left[\xi \mathrm{log}\frac{\xi \sqrt{\mathrm{\Delta }^2+\xi ^2}}{\xi +\sqrt{\mathrm{\Delta }^2+\xi ^2}}+2\sqrt{\mathrm{\Delta }^2+\xi ^2}\right]=0.$$
(74)
This curve cuts the real axis inside the interval $`[\omega ,\omega ]`$, which contradicts the assumptions made so far. From numerical studies we can see that for weak couplings only a fraction of energies $`E_\mu `$ form complex conjugated pairs, which in the limit $`N\mathrm{}`$ organize themselves into arcs, while the other energies remain real and located between the lower energy levels (see the appendix for more details on this issue). Taking this into account and using contour arguments, it can be shown that the BCS eqs. (43), (44), and (45) also hold in cases where $`\mathrm{\Gamma }`$ cuts the energy interval $`\mathrm{\Omega }`$.

In the general case the curve $`\mathrm{\Gamma }`$ is given by the equation

$$\left[i\xi \mathrm{arcsin}\sqrt{\frac{1+\xi ^2/\mathrm{\Delta }^2}{1\xi ^2/\omega ^2}}+\omega \mathrm{Argsinh}\sqrt{\frac{\mathrm{\Delta }^2+\xi ^2}{\omega ^2\xi ^2}}\right]=0.$$
(75)
In the limit $`\mathrm{\Delta }\omega `$ eq. (75) turns into (74). We show in the appendix that for $`\mathrm{\Delta }>\omega `$ the arc $`\mathrm{\Gamma }`$ does not touch the interval $`[\omega ,\omega ]`$ and all the roots $`E_\mu `$ are complex. This happens for $`g>g_0=1.13459`$ $`(\mathrm{sinh}1/g_0=1)`$.

## VII) Numerical results

In this section we present the comparison between the analytic results derived in section VI, and the numerical ones obtained by solving directly the Richardson’s eqs. (6) for the models considered above.

### Two-level model

In this model there are two energy levels, $`\pm \epsilon _1`$, with a degeneracy $`N/2`$. Hence the Richardson eqs. (6) become

$$\frac{1}{G}=\frac{N}{2}\left(\frac{1}{\epsilon _1E_\nu }\frac{1}{\epsilon _1+E_\nu }\right)\underset{\mu \nu }{\overset{M}{}}\frac{2}{E_\mu E_\nu },$$
(76)
with $`M=N/2`$. The total number of solutions of this system of equations is given by the combinatorial number $`C_M^N`$. The solution which corresponds to the ground state of the BCS Hamiltonian is the one for which $`E_\mu (G)\epsilon _1,\mu `$ in the limit where $`G0`$. For any non zero value of $`G`$, all the roots $`E_\mu `$ form complex conjugate pairs surrounding the lower energy $`\epsilon _1`$.

The reduced coupling constant for this case is taken as $`g=GN`$. Our numerical results are presented in units of $`\epsilon _1`$, which is equivalent to set $`\epsilon _1=1`$.

In fig. 1 we plot the distribution of the $`M=100`$ roots $`E_\mu `$, for three different values of the coupling constant $`g=GN`$, together with the analytic curves derived in section VI. There is an optimal agreement between the numerical and analytic results in the three regimes: normal ($`g<\epsilon _1`$), critical ($`g=\epsilon _1`$), and superconducting ($`g>\epsilon _1`$). As we increase the value of $`M`$ up to 800 pairs, the fit of the numerical data to the analytical curves improves, and in the limit $`M\mathrm{}`$ the discrete roots collapse into the continuous curves. In particular, the complex roots $`E_{max}(M)`$ and $`E_{max}^{}(M)`$, whose real part lies closer to $`\epsilon _0`$, approach their continuum limit value, namely

$$\underset{M\mathrm{}}{lim}E_{max}(M)=\epsilon _0+i\mathrm{\Delta },$$
(77)
where

$`\epsilon _0`$ $`=`$ $`\{\begin{array}{ccc}0& & g\epsilon _1\hfill \\ \epsilon _1\sqrt{1g/\epsilon _1}& & g\epsilon _1\hfill \end{array},`$ (80)
$`\mathrm{\Delta }`$ $`=`$ $`\{\begin{array}{ccc}\sqrt{g^2\epsilon _1^2}& & g\epsilon _1\hfill \\ 0& & g\epsilon _1\hfill \end{array},`$ (83)
in agreement with eqs. (51) and (64).

Fig. 2 shows the following scaling behaviour $`E_{max}(M)`$:

$$E_{max}(M)(\epsilon _0+i\mathrm{\Delta })\frac{1}{M^\theta },$$
(84)
where the exponent $`\theta `$ depends on the regime of the system. The numerical and analytic results of $`\theta `$ are given in table 1:

| $`\theta `$ | $`g=1.5\epsilon _1`$ | $`g=1.0\epsilon _1`$ | $`g=0.5\epsilon _1`$ |
| --- | --- | --- | --- |
| analytical | 2/3 | 1/3 | 1/2 |
| numerical $`(\mathrm{\Delta })`$ | 0.635 | 0.329 | 0.446 |
| numerical $`(\epsilon _0`$) | 0.672 | 0.328 | 0.513 |

Table 1.- Numerical and analytical values of the exponent $`\theta `$ appearing in eq. (84). The numerical values are extracted from the real and imaginary parts of $`E_{max}(M)`$.

The analytical value of $`\theta `$ can be derived as follows. For $`M`$ sufficiently large we can write

$$E_{max}(M)=(\epsilon _0+i\mathrm{\Delta })+\delta \xi ,$$
(85)
where $`\delta \xi `$ is a small number. Using eqs. (53), (56) and (58), the density of roots $`r(E_{max})`$ behaves as

$$r(E_{max})AM(\delta \xi )^\nu ,\nu =\{\begin{array}{cc}1/2& g>\epsilon _1\\ 2& g=\epsilon _1\\ 1& g<\epsilon _1\end{array},$$
(86)
where $`A`$ is a parameter which only depends on $`g`$ and $`\epsilon _1`$. The magnitude of $`\delta \xi `$ can be fixed by imposing the existence of at least one root in the interval $`(\epsilon _0+i\mathrm{\Delta },\epsilon _0+i\mathrm{\Delta }+\delta \xi `$), namely

$$1=r(E_{max})|\delta \xi |AM|\delta \xi |^{\nu +1}\delta \xi \frac{1}{M^{\frac{1}{1+\nu }}},$$
(87)
which is the desired result (84) with $`\theta =1/(1+\nu )`$.

As we have seen in fig. 1, the roots $`E_\mu `$ lie on the analytic curves to a very good degree of approximation. A further check of the analytic solution is to compare the density of roots $`r(\xi )`$, which is given by $`r(E)=|h(E)|/\pi `$ for $`g\epsilon _1`$ (see eq. (37)), and by $`r(\xi )=|s(\xi )|`$ for $`g\epsilon _1`$ (see eq. (56)). In fig. 3 we plot the fraction of roots per unit length $`r(E)/M`$ as a function of the arc length $`s[1,1]`$, where the $`s=\pm 1`$ corresponds to the end points of the arc $`\mathrm{\Gamma }`$. The agreement is extremely accurate for $`M=800`$.

The value of the conserved quantities obtained in eq. (53) also agree with the numerical ones in the large $`N`$ limit, in the three different regimes.

Finally, it is interesting to study the behaviour of the energy as a funtion of the number of pairs $`M`$. The leading term is given by eq. (52) for $`g\epsilon _1`$ and by eq. (60) for $`g\epsilon _1`$. In fig. 4 we plot $`e_{\mathrm{}}e(M)`$ versus $`1/M`$, where $`e(M)=E(M)/M`$ is the energy per pair. The linear behaviour of the numerical data means that the next leading correction to $`E(M)`$ is a constant term, $`E_1`$, whose value is given by the analytic expression

$$E_1/M=e_{\mathrm{}}e(M)=\{\begin{array}{ccc}(g\frac{1}{2}\mathrm{\Delta })/M\hfill & & g\epsilon _1\hfill \\ (\epsilon _1|\epsilon _0|)/M\hfill & & g\epsilon _1\hfill \end{array},$$
(88)
where $`\mathrm{\Delta }`$ and $`\epsilon _0`$ are given by eqs. (51) and (64) respectively. The formula for $`g>\epsilon _1`$ has been derived by Richardson in ref. \[R-limit, \], studying the next to leading order of the electrostatic model. The formula for $`g<\epsilon _1`$ is an educated guess which needs a proof. At $`g=\epsilon _1`$, the numerical data show a large deviation, for small $`M`$, from the formula (88), namely $`\epsilon _1/M`$. There are theoretical reasons to believe that the scaling behaviour of $`E_1/M`$ at $`g=\epsilon _1`$ , is given by the following formula (see fig. 4):

$$E_1/M=\epsilon _1(\frac{1}{M}\frac{A}{M^{4/3}},),A=0.62258,$$
(89)
which would imply that the next-to-next leading correction of the total energy goes as $`M^{1/3}`$, rather than $`M^1`$. In fact, the formula given by Richardson R-limit  for the $`O(M^1)`$ energy correction, which is given by $`E_2=g(3g\mathrm{\Delta })(g\mathrm{\Delta })/(4N\mathrm{\Delta }^2)`$, breaks down at $`g=\epsilon _1`$ where the gap vanishes. The derivation of eq. (89) shall be given elsewhere.

### Equally spaced model

This model is defined by $`N`$ non degenerate energy levels $`\epsilon _j=d(2jN1),j=1,\mathrm{},N`$, where $`d=\omega /N`$ is the single particle level spacing. The first numerical study of this model, using the exact solution, was done by Richardson in reference R-roots , where he introduced new set of variables, related to the roots $`E_\mu `$, in order to handle the singularities arising at particular values of the coupling constant $`g`$ (see below). This work was revised recently in connection to ultrasmall metallic grains in \[vDB1, ; vDB2, \].

The ground state of the system, at half-filling $`M=N/2`$, is the one for which $`lim_{g0}E_\mu (g)\epsilon _\mu (\mu =1,\mathrm{},N/2)`$.

The reduced coupling constant is given by $`g=G/d=GN/\omega `$ in this section. And all our numerical results are presented in units of $`\omega `$, which is equivalent to set $`\omega =1`$.

In fig. 5 we plot the real and imaginary parts of $`E_\mu `$, for a system with $`M=N/2=8`$ pairs. For small values of $`g`$, all the pair energies are real and below the corresponding energy levels at $`g=0`$. At a certain value $`g_1`$, the two closest pairs to the Fermi energy, namely $`E_M`$ and $`E_{M1}`$, coincide at $`\epsilon _{M1}`$ and then, for larger values of $`g>g_1`$ they form a complex conjugate pair, while the rest of the pair energies remain real. At higher critical values of $`g`$ the same phenomena happens for the other energies, until all of them become complex. At the critical values of $`g`$ the Richardson eqs. (6) are singular, but the singularities can be resolved by the change of variables introduced in \[R-roots, \].

In fig. 6 we show the distribution of roots for $`M=100`$ pairs, together with the theoretical curves obtained in section VI, for three different values of $`g`$. For $`g=0.5`$ and $`1`$ there is a fraction of roots which become complex, falling into the arcs $`\mathrm{\Gamma }`$ described by eq. (75), while the real roots extend from the cutting point of the arc $`\mathrm{\Gamma }`$ with the real axis down to the energy $`\omega `$. For $`g=1.5`$ all the roots are complex and fall into the theoretical curve (75). The numerical data are consistent with the value $`g_0=1.13459`$ above which all the energies become complex.

As in the two-level case, the root $`E_{max}(M)`$ closest to $`\epsilon _0+i\mathrm{\Delta }`$ satisfy the scaling law (84) with an exponent $`\theta =2/3`$. This result can be proved along the same lines as for the two-level case. In table 2 we give the comparison between the theoretical and numerical values of $`\theta `$, and in fig. 10 we show the scaling behaviour for three values of $`g`$.

| $`\theta `$ | $`g=1.5`$ | $`g=1`$ | $`g=0.5`$ |
| --- | --- | --- | --- |
| analytical | 2/3 | 2/3 | 2/3 |
| numerical $`(\mathrm{\Delta })`$ | 0.646 | 0.645 | 0.642 |
| numerical $`(\epsilon _0`$) | 0.658 | 0.659 | 0.661 |

Table 2.- Numerical and analytical values of the exponent $`\theta `$ appearing in eq. (84) for the equally spaced model. The notation is the same as in table 1.

In fig. 8 we give the fraction of roots per unit length $`r(E)/M`$ for $`M=800`$ and three values of $`g`$. Again, the numerical data fit accurately the theoretical result $`r(E)/M=|h(E)|/\pi `$, where $`h(E)`$ is given by eq. (71).

A further check is presented in fig. 9, where the distribution of the conserved quantities $`\lambda (\epsilon _i)`$ for $`M=20`$ is presented, together with their continuous limit curves given in eq. (73). A very good agreement is achieved even for such a small number of pairs.

We have finally studied the condensation energy $`E_C/M`$ as a function of $`M`$, obtaining the behaviour $`E_C(M)=e_{\mathrm{}}M+E_1`$, where $`E_1(g)`$ is plotted in fig. 10. An analytic expression for $`E_1(g)`$ is not known R-limit . In order to understand the numerical results we have made two fits. The first one is based on the result obtained for the two-level case, namely

$$E_1(g)=g\omega \frac{\mathrm{\Delta }}{2},$$
(90)
where $`\mathrm{\Delta }`$ is given by eq. (67). Fig. 10 shows that (90) is a good fit for large values of $`g`$, where all the roots are complex, forming a single arc. For $`gg_0=1.13459`$ there is a fraction $`M_{complex}/M`$ of complex roots, and a fraction $`M_{real}/M`$ of real roots. This suggests the following mixed fit:

$`E_1(g)/\omega `$ $`=`$ $`\left(g{\displaystyle \frac{\mathrm{\Delta }}{2\omega }}\right){\displaystyle \frac{M_{complex}}{M}}+(A+Bg)g{\displaystyle \frac{M_{real}}{M}},`$
$`A`$ $`=`$ $`0.6735,B=0.1729,`$
which describes pretty well $`E_1(g)`$ in the whole range of $`g^{}s`$. In the weak coupling regime there have been proposed several formulas for the finite size corrections of the condensation energy $`E_C`$ using the DMRG method DMRG  and the exact solution Sch . In the latter reference it is shown that for couplings $`g>g=1/lnN`$, the condensation energy $`E_C`$ can be splitted in two contributions, one due to the complex energy pairs ( the BCS term) and another one due to the real roots ( perturbative term). It would be interesting to clarify the relation between formula (LABEL:101) and eq.(24) in reference Sch .

## VIII) Conclusions and prospects

In this paper we have reviewed and completed Gaudin’s continuum limit of the exactly solvable BCS model, and compared its predictions with the numerical solutions describing the ground state of the two-level and the equally spaced models. We have confirmed previously known results R-roots ; random ; vDB1 ; vDB2 , and obtained new ones, which we list below:

* Analytic and numerical determination of the curves $`\mathrm{\Gamma }`$, where the roots $`E_\mu `$ of the Richardson equations lie, together with their density along $`\mathrm{\Gamma }`$. For the equally spaced model the curves $`\mathrm{\Gamma }`$ are given for arbitrary values of the coupling constant $`g`$, and for the two-level model in the normal regime, i.e. $`g<\epsilon _1`$, the analytic curves $`\mathrm{\Gamma }`$ differ from those in \[G-book, \].
* Analytic expression of the eigenvalues $`\lambda _j`$ of the conserved quantities found by Cambiaggio, Rivas, and Saraceno CRS , which fit accurately the numerical data, even for a few levels.
* Study of the scaling properties of the ground state energy $`E(M)`$ with the number of pairs $`M`$ which, in the large $`M`$ limit, behaves generically as $`E(M)=e_{\mathrm{}}M+E_1+O(1/M)`$. For the two-level model, in the superconducting regime, we find an agreement with the analytic result for $`E_1`$ \[R-limit, \], while in the normal regime we have guessed an analytic expression for $`E_1`$, which fits perfectly the numerical data. In the critical case there are large deviations in $`E(M)`$, which can be explained by a term $`O(M^{1/3})`$. In the equally spaced model the condensation energy also receives a constant contribution, $`E_1`$, which is fitted with a “phenomenological” formula which combines the effects of the complex and real roots.
* Scaling behaviour of the roots $`E_{max}(M)`$ which lie closest to the end points of the curves $`\mathrm{\Gamma }`$. This scaling can be characterized by a critical exponent $`\theta `$, whose analytic values fit well the numerical ones. In the superconducting regimes (i.e. $`g>0`$ for the equally spaced-model and $`g>\epsilon _1`$ for the two-level model) we find $`\theta =2/3`$, while in the normal (critical) regimes of the two-level model we find $`\theta =1/2`$ ($`1/3`$). In this manner the exponent $`\theta `$ characterizes the nature of the ground state. An interesting question is wether $`\theta `$ shows up in physical observables.
* Generalization of Gaudin’s conformal mapping for the equally spaced model using elliptic functions, which degenerate into trigonometric ones in the limit $`\mathrm{\Delta }\omega `$. This suggest the existence of a rich mathematical structure.

Finally, we shall mention some problems which are worth to investigate along the lines suggested by the present work:

* Study of the excited states and the finite size effects in a more geometrical way, completing the work initiated by Richardson in \[R-limit, \].
* Quantum dots at finite temperature is another difficult problem, which has been treated combining several techniques depending on the temperature regime Finite-T ; FFM ; FF . Is it possible to find a thermodynamic Bethe ansatz for this system?
* Study of solutions of Richardson equations with several arcs, i.e. $`K>1`$ in the notation of section V. For the standard BCS model they must describe very high excited states formed by separate condensates in interaction. This case may be relevant to systems such as arrays of superconducting grains or quantum dots. From a mathematical point of view, the cases with $`K>1`$ seem to be related to the theory of hyperelliptic curves and higher genus Riemann surfaces, which may shed some light on this physical problem.
* Generalization of Gaudin’s equations to Richardson models based on arbitrary Lie groups $`𝒢`$ AFS . The standard case corresponds to the choice $`𝒢=SU(2)`$. The electrostatic model for generic $`𝒢`$ has charges labelled by vectors of dimension $`\mathrm{rank}𝒢`$. It is naturally to expect that the roots $`E_\mu ^a(a=1,\mathrm{},\mathrm{rank}𝒢)`$ should fall into several arcs labelled by the index $`a=1,\mathrm{},\mathrm{rank}𝒢`$.

## Acknowledgments

We thank João Lopes dos Santos, Rodolfo Cuerno, Miguel Angel Martín-Delgado, and Javier Rodríguez-Laguna for computational help, and Didina Serban for the reference \[G-book, \]. This work has been supported by the Spanish grants BFM2000-1320-C02-01 (JMR and GS), and BFM2000-1320-C02-02 (JD).

## Appendix: Conformal mapping for the equally spaced system

In this appendix we study the analytic structure underlying the shapes of the curves defined by eqs. (74) and (75).

### i) Case $`\mathrm{\Delta }\omega `$

In this limit it is convenient to make a conformal mapping G-book  from the upper half-plane $`\mathrm{Im}\xi >0`$, with a cut in the segment $`(0,i\mathrm{\Delta })`$, into the half-strip $`\lambda >0`$, $`\pi /2<\mu <\pi /2`$,

$$\xi =i\mathrm{\Delta }\mathrm{cosh}(\lambda +i\mu ).$$
(92)
In the plane $`u=\lambda +i\mu `$ the curve (74) becomes

$$\lambda \mathrm{tanh}\lambda +\mu \mathrm{cot}\mu =0.$$
(93)
Figures 11 and 12 show, in the $`u`$ and $`\xi `$ planes, the arcs defined by eqs. (93) and (74). The point $`𝐎`$ ($`\xi =\epsilon _0+i\mathrm{\Delta })`$ is mapped into the origin ($`u=0`$), where the field $`h`$ vanishes. In fact $`h(\xi )`$, as a function of $`u`$, is simply $`h=2\rho _0u`$. The point $`𝐎^{}`$ ($`\xi =\epsilon _0i\mathrm{\Delta })`$ is the mirror reflection of $`𝐎`$. The arc $`\mathrm{𝐎𝐀}`$ is the one that corresponds to the actual numerical solution for $`\mathrm{Im}\xi >0`$, while its mirror image, $`𝐎^{}𝐀`$, is the one for $`\mathrm{Im}\xi <0`$. We thus see that the whole arc $`\mathrm{𝐎𝐀𝐎}^{}`$ cuts the real axis at the point $`𝐀`$ with an energy

$$\epsilon _A=\epsilon _0\mathrm{\Delta }\mathrm{sinh}\lambda _0,\lambda _0\mathrm{tanh}\lambda _0=1.$$
(94)
These results are expected from numerical studies, which show that for weak couplings only a fraction of energies $`E_\mu `$ form complex conjugated pairs, which in the limit $`N\mathrm{}`$ organize themselves into arcs, while the other energies remain real and located between the lower energy levels ( see fig. 6). We can check this result computing the number of real and complex roots $`E_\mu `$. From eq. (94) we deduce that the real roots occupy the interval $`[\omega ,\mathrm{\Delta }\mathrm{sinh}\lambda _0]`$, with a density which is $`2\rho _0`$. Hence the total number, $`M_{\mathrm{real}}`$ of real roots is

$$M_{\mathrm{real}}=2\rho _0(\omega \mathrm{\Delta }\mathrm{sinh}\lambda _0).$$
(95)
Since the total number of roots is $`M=M_{\mathrm{real}}+M_{\mathrm{complex}}=N/2=2\rho _0\omega `$, we deduce that

$$M_{\mathrm{complex}}=2\mathrm{\Delta }\rho _0\mathrm{sinh}\lambda _0.$$
(96)
This result agrees with the one obtained integrating the field $`h(\xi )`$ along the curve $`L`$ surrounding $`\mathrm{\Gamma }`$ (recall eq. (44)). This integration is more easily done in the $`u`$-strip

$$M_{\mathrm{complex}}=_Lh(\xi )𝑑\xi =4\times _0^{\lambda _0+i\pi /2}(2\rho _0u)(i\mathrm{\Delta }\mathrm{sinh}udu),$$
(97)
where the factor 4 comes from the contribution of the up and down, and left and right pieces of the contour $`L`$. Using contour arguments, it can be shown that the BCS eqs. (43), (44), and (45) also hold in cases where $`\mathrm{\Gamma }`$ cuts the energy interval $`\mathrm{\Omega }`$. As a side comment, we observe that $`M_{\mathrm{complex}}=\frac{\mathrm{\Delta }}{2d}\mathrm{sinh}\lambda _0`$, agrees with the heuristic argument stating that there are roughly $`\mathrm{\Delta }_{BCS}/d`$ energy levels around the Fermi level that are strongly affected by the pairing interactions T . See also reference Sch  for an approximate equation giving the number of complex pairs ( eq. (49)).

### ii) Generic case

In order to compare with the numerical results presented in the section VII, it is convenient not to make the approximation $`\mathrm{\Delta }\omega `$. The appropiate change of variables that generalizes (92) is given by

$$\xi =\frac{i\mathrm{\Delta }\mathrm{cosh}u}{\sqrt{1\left({\displaystyle \frac{\mathrm{\Delta }}{\omega }}\mathrm{sinh}u\right)^2}},$$
(98)
in terms of which $`h=2\rho _0u`$. The integral $`h(\xi )𝑑\xi `$ can be more easily done in the $`u`$ plane, yielding for the equation of the curve $`\mathrm{\Gamma }`$

$$\left[\frac{iu\mathrm{cosh}u}{\sqrt{1\left({\displaystyle \frac{\mathrm{\Delta }}{\omega }}\mathrm{sinh}u\right)^2}}\frac{i\omega }{\mathrm{\Delta }}\mathrm{arcsin}\left(\frac{\mathrm{\Delta }}{\omega }\mathrm{sinh}u\right)\right]=0,$$
(99)
which in the $`\xi `$-plane gives eq. (75), namely

$$\left[i\xi \mathrm{arcsin}\sqrt{\frac{1+\xi ^2/\mathrm{\Delta }^2}{1\xi ^2/\omega ^2}}+\omega \mathrm{Argsinh}\sqrt{\frac{\mathrm{\Delta }^2+\xi ^2}{\omega ^2\xi ^2}}\right]=0.$$
(100)
In the limit $`\mathrm{\Delta }\omega `$ this equation turns into eq. (74). We can ask when this curve cuts the interval $`[\omega ,\omega ]`$. A convenient parametrization of the cutting point is (recall eq. (94))

$$\xi _A=\frac{\mathrm{\Delta }\mathrm{sinh}\lambda _0}{\sqrt{1+\left({\displaystyle \frac{\mathrm{\Delta }}{\omega }}\mathrm{cosh}\lambda _0\right)^2}}.$$
(101)
Plugging (101) into (100) we deduce the equation for $`\lambda _0`$,

$$\frac{\lambda _0\mathrm{sinh}\lambda _0}{\sqrt{1+\left({\displaystyle \frac{\mathrm{\Delta }}{\omega }}\mathrm{cosh}\lambda _0\right)^2}}=\frac{\omega }{\mathrm{\Delta }}\mathrm{Argsinh}\left(\frac{\mathrm{\Delta }}{\omega }\mathrm{cosh}\lambda _0\right),$$
(102)
which only has real solutions provided $`\mathrm{\Delta }\omega `$. When $`\mathrm{\Delta }>\omega `$ the arc $`\mathrm{\Gamma }`$ does not touch the interval $`[\omega ,\omega ]`$ and all the roots $`E_\mu `$ are complex. This happens for $`g>g_0=1.13459`$. This value of $`g_0`$ is obtained solving the equation $`\mathrm{sinh}1/g_0=1`$.

To complete this appendix we shall show that eq. (98) yields a conformal mapping that generalizes the one found by Gaudin for generic values of $`\mathrm{\Delta }/\omega `$. To do so, let us first perform the change of variables

$$u=i\left(\frac{\pi }{2}\varphi \right),$$
(103)
and introduce the parameters

$$k=\frac{\mathrm{\Delta }}{\sqrt{\mathrm{\Delta }^2+\omega ^2}},k^{}=\frac{1}{\sqrt{\mathrm{\Delta }^2+\omega ^2}}=\sqrt{1k^2},$$
(104)
in terms of which (98) becomes

$$\xi =\frac{i\mathrm{\Delta }k^{}\mathrm{sin}\varphi }{\sqrt{1k^2\mathrm{sin}^2\varphi }}.$$
(105)
Using the elliptic functions with modulus $`k`$, and conjugate modulus $`k^{}`$, we define a new variable $`z`$ through WW

$`\mathrm{sin}\varphi `$ $`=`$ $`\mathrm{sn}(z,k),`$
$`\mathrm{cos}\varphi `$ $`=`$ $`\mathrm{cn}(z,k),`$ (106)
$`\sqrt{1k^2\mathrm{sin}^2\varphi }`$ $`=`$ $`\mathrm{dn}(z,k).`$
The relation between $`z`$ and $`\varphi `$ is given by the equation

$$z(k,\varphi )=_0^\varphi \frac{d\psi }{\sqrt{1k^2\mathrm{sin}^2\psi }},$$
(107)
while the inverse function is denoted as $`\varphi =\mathrm{am}(z,k)`$. The relation between $`\xi `$ and $`z`$ reads

$$\xi =i\mathrm{\Delta }k^{}\frac{\mathrm{sn}(z,k)}{\mathrm{dn}(z,k)}.$$
(108)
Finally, using the relation $`\mathrm{sn}(K+v,k)/\mathrm{dn}(K+v,k)=\mathrm{cn}(v,k)/k^{}`$, where $`K`$ is the half-period of the elliptic integrals, eq. (108) turns into

$$\xi =i\mathrm{\Delta }\mathrm{cn}(v,k),$$
(109)
where $`z=K+v`$. The relation between the variables $`u`$ and $`v`$ is given by

$$iu=\frac{\pi }{2}+\mathrm{am}(K+v,k)=\frac{\pi v}{2K}+\underset{n=1}{\overset{\mathrm{}}{}}\frac{2q^n}{n(1+q^{2n}}\mathrm{sin}\frac{n\pi (K+v)}{2K},$$
(110)
where $`q=e^{\pi K^{}/K}`$. In the limit $`\mathrm{\Delta }/\omega 1`$ we get $`k0`$, $`K\pi /2`$, and $`q0`$, which implies that $`iuv`$ and eq. (109) reduces to (92).

Eq. (109) gives a conformal map of the rectangle $`K<\mathrm{Re}v<K,\mathrm{\hspace{0.33em}0}<\mathrm{Im}v<K^{}`$ onto the upper half-plane $`\mathrm{Im}\xi >0`$, with a cut along the linear segment $`0<\mathrm{Im}\xi \mathrm{\Delta }`$ (see reference \[N, \] and fig. 13). Some examples of this mapping are given by

$$v=0,\pm K,\pm K+iK^{}\xi =i\mathrm{\Delta },0,\pm \omega ,$$
(111)
which imply the following identifications between the sides of the rectangle in the $`v`$-plane and the relevant energy regions in the $`\xi `$-plane:

| cut | $`\begin{array}{c}\mathrm{Re}\xi =0\\ 0<\mathrm{Im}\xi <\mathrm{\Delta }\end{array}`$ | $`\begin{array}{c}K<\mathrm{Re}v<K\\ \mathrm{Im}v=0\end{array}`$ |
| --- | --- | --- |
| lower band | $`\begin{array}{c}\omega <\mathrm{Re}\xi <0\\ \mathrm{Im}\xi =0\end{array}`$ | $`\begin{array}{c}\mathrm{Re}v=K\\ 0<\mathrm{Im}v<K^{}\end{array}`$ |
| upper band | $`\begin{array}{c}0<\mathrm{Re}\xi <\omega \\ \mathrm{Im}\xi =0\end{array}`$ | $`\begin{array}{c}\mathrm{Re}v=K\\ 0<\mathrm{Im}v<K^{}\end{array}`$ |
| off band | $`\begin{array}{c}|\mathrm{Re}\xi |>\omega \\ \mathrm{Im}\xi =0\end{array}`$ | $`\begin{array}{c}K<\mathrm{Re}v<K\\ \mathrm{Im}v=K^{}\end{array}`$ |

Table 3.- Correspondences of the boundary regions of the conformal map (109).