Ground-state energy of a Richardson-Gaudin integrable BCS model

We investigate the ground-state energy of a Richardson-Gaudin integrable BCS model, generalizing the closed and open p+ip models. The Hamiltonian supports a family of mutually commuting conserved operators satisfying quadratic relations. From the eigenvalues of the conserved operators we derive, in the continuum limit, an integral equation for which a solution corresponding to the ground state is established. The energy expression from this solution agrees with the BCS mean-field result.

better known ancestor the s-wave pairing model [2] which is of Richardson-Gaudin type [4][5][6], the p x +i p y BCS Hamiltonian also supports a family of mutually commuting operators [7]. The Bethe Ansatz solution of the p x + i p y -pairing model was initially established with conserved particle number and no interaction with the environment [1,8,9]. For this reason, we call the integrable p x + i p y BCS model the 'closed' p + i p model. An extension of the closed model was later studied where the model is coupled to its environment. This extended model no longer conserves particle number, and thus cannot be considered as describing a closed system. Therefore we refer to it as an open model. The interaction term corresponding to particle exchange with the system's environment leads to u(1) symmetry being broken, which generally renders the analysis of the system to be more complicated.
Integrability of the open model was established in [10] through use of the boundary Quantum Inverse Scattering Method. An alternative derivation, which is less technical, was later provided in [11].
The p x + i p y -wave interaction also gives rise to non-trivial topological phase transitions by breaking time-reversal symmetry [12,13]. Recently, there has been much attention to the topological properties of this model using various methods. In particular, the exact Bethe Ansatz solutions provide insight into the topological behaviour of the system through the peculiar behaviour of Bethe roots [7,9,14]. Topological properties of the open model have been studied in [15,16].
In the recent work by Skrypnyk [17], a generalization of the closed and open p + i p models was constructed by considering a non-skew-symmetric r-matrix satisfying a generalized classical Yang-Baxter equation. The underlying conserved operators of this generalized model were also independently constructed by Claeys et al [18], as a spin-1/2 Richardson-Gaudin XYZ solution where a non-zero magnetic field component is present. This generalized model opens up an avenue to investigate more complicated interaction with the environment.
Originating in the work of Babelon and Talalaev [19], it was shown, through a change of variables, that the Bethe Ansatz Equations (BAE) for certain Richardson-Gaudin type systems can be recast into a set of coupled polynomial equations. The roots of these equations can be expressed in the form of eigenvalues of the self-adjoint conserved operators, and consequently are real-valued. The same form of polynomial equations were adopted in [20] as means for efficient numerical solution of the conserved operator spectrum, and in [21] to compute wavefunction overlaps. This technique was shown to be generalizable in [22,23]. It was subsequently shown for some systems that the conserved operator eigenvalues (COE), satisfying quadratic relations, are inherited from the same relations at the operator level [15,24]. The generalized model considered below also shares this property of underlying conserved operators satisfying quadratic relations [17,18]. However the polynomial equations for the generalized model, leading to the BAE, remain unknown.
The conventional Bethe Ansatz analysis of Richardson-Gaudin models in the thermodynamic limit suffers complications in that one needs to solve for a density function on an unknown curve in the complex plane which captures the behaviors of the Bethe roots [25][26][27]. For the closed model, the complex curve and its density were studied extensively in [9] for different regions of the phase diagram, with the resulting ground-state energies shown to be consistent with mean-field analysis. However, limitations of this method remain and are discussed in [28]. For the open model, even more intricate behavior of the Bethe roots is revealed from numerical analysis which renders it a difficult task to capture the Bethe roots with a piecewise continuous complex curve. Alternatively, we proposed a new method [28] for deriving the continuum limit ground-state energy in terms of the COE, initially for the closed model. This method does not rely on any hypothesis on the distribution of Bethe roots, hence its validity extends to regions of the phase diagram where the Bethe root evolution is not well-described by a complex curve. This new method was also shown to accommodate the open model.
The objective of our study is to apply this new approach to derive the ground-state energy of the generalized model in the continuum limit. The main focus of the task is solving the integral equation arising from the quadratic relations of the COE in the continuum limit. The proof of the solution will be given in detail, which is a generalization of the similar approach for the open case. Here it will also be shown that the result is in agreement with mean-field calculations.
The generalized integrable BCS Richardson-Gaudin Hamiltonian is introduced in Sect. 2 with its underlying conserved operators satisfying quadratic relations. In Sect. 3, results from a BCS mean-field analysis of the generalized model, including the ground-state energy, are given. In Sect. 4, attention turns towards using the alternative approach to derive the groundstate energy, assisted by establishing a solution to the integral equation which corresponds to the ground state of the model. Concluding remarks and discussion are offered in Sect. 5.

The Hamiltonian
Consider a system of fermions and let c α j , c † α k , α, α ∈ {+, −}, j, k ∈ 1, 2, . . . , L denote the annihilation and creation operators, satisfying We consider the following integrable BCS Richardson-Gaudin Hamiltonian reported in [17], where g, γ, λ ∈ and with ε k > 0, ε k + β x > 0 and ε k + β y > 0 for all k = 1, . . . , L. Assuming β x ≥ β y and setting The Hamiltonian has the form of a standard BCS model, (the terms in the first line of (1)), with additional terms. The standard terms are associated with the single-particle energy spectrum and the scattering of Cooper pairs. The additional terms are responsible for the breaking of the u(1)-symmetry associated to the total particle number, and are interpreted as interaction with the system's environment. The Hamiltonian (1) is a generalisation of the open p + i p Hamiltonian [10]. By setting λ = 0 and we recover the conserved operators underlying the open p + i p model. Furthermore, if we also set γ = 0 these reduce to those of the closed p + i p model. Now we replace the fermion operators with the generators S + k , S − k and S z k of the sl(2) ⊗L algebra, through where k = 1, . . . , L. These operators satisfy the following commutation relations The Hamiltonian H BC S can be rewritten (omitting the constant term) in the following form Define the following operators with S x i , S y i given by It can be shown that where {Q i } is a set of mutually commuting conserved operators. These operators have been shown [17,18] to satisfy the following quadratic relations: The eigenvalues {q i } corresponding to {Q i } give the energy expression and due to (5), satisfy the relation

Mean-field analysis
The Hamiltonian H in (3) can be rewritten as Introducing the parameter G = −g in the following calculation, we proceed with the mean-field analysis of the model. The techniques below extend the original methods for the s-wave BCS model [2]. First let where ∆ x , ∆ y ∈ . Then linearizing (8) through The detailed calculations are given in Appendix A. The ground-state energy reads subject to the 'gap' equations That is, given G, γ and λ the two equations (10) and (11) The 'gap' equations for the open model read If we require a non-zero γ, it is straightforward from these two equations to conclude that ∆ y = 0, i.e. we have one 'gap' equation In this case, we recover the same mean-field results, i.e. energies and gap equations, as in [28] for the open model. Next we will introduce an integral approximation to show that the general mean-field ground-state energy (9) is consistent with the exact solution. The techniques used require an extension of those in [28], to account for the two gap equations.

Integral approximation of the quadratic identities
Recall that for an arbitrary function F (x), we may consider the following integral approximation (or continuum limit) of a summation, where a, b are the lower and upper bounds for the set of real numbers {a = x 1 < x 2 < · · · < x L = b} and ρ(x) is a density function introduced for the distribution of x i such that In order to proceed with the integral approximation of (7), we assume that approach finite limits as L → ∞. We will then require the parameter G ∼ 1/L, such that G L remains finite. The last term in (7) then vanishes as L → ∞ since G 2 β 2 L ∼ 1/L. We will adopt the notation G = G L in the following calculations. The Hamiltonian (1) is a function of many independent variables, including {z i }. It is physically plausible to interpret these particular variables as the momentum spectrum in the non-interacting limit with β = 0. Consider a density ρ for the distribution of the variables {z 2 i }. This density is associated with the kinetic energies in the non-interacting limit with β = 0. The continuum limit for the quadratic identity (7) then becomes where ω, ω 0 > β are the upper and lower bounds introduced for the continuous variable replacing {z 2 i } in the continuum limit, and Assuming that ∆ x , ∆ y ∼ L, the continuum limit for the gap equations (10) and (11) read where∆ x = ∆ x /L,∆ y = ∆ y /L, and we have also set

Solution of the integral equation
The main result of our study is the following: Proposition 1. The expression is a solution to the integral equation (13) with the 'gap' equations (14) satisfied. Furthermore the continuum limit of (6), corresponds to the ground-state energy of the model. Proof: We set The main challenge is to calculate the following term in (13), Rewriting the integral equation (13) as our task is to show that (17) holds. We adopt the following change of variables: That is, parameters χ x , χ y replace γ, λ in the following calculations. We begin with the square term in (17) We then calculate the following term in (17) Combining the gap equations (14), we set the following shorthand First, Likewise By exploiting symmetries, the first term in (24) can be reduced as where we used the gap equations (14) in the last step.
Using a similar approach, we compute the second term in (24), where the shorthand " " for " " is adopted and the last step is assisted by symbolic calculation in Mathematica.
Similarly, it can be shown that Before proceeding, we set Substituting (18) and (19) into (17), it remains to verify that From (22), we simplify C 2 to obtain Likewise, from (23) expression C 3 can be reduced to Substituting (25) and (26) into (24), and also from (27) expression C 4 can be simplified as Lastly, using (20) and (21), it can be shown that (28) holds, thus concluding our proof that q( ) given by (16) is a solution to the integral equation (13). Next, using the gap equations (14), the energy expression is calculated from (6) in the following we obtain the energy expression as Hence the integral approximation of the conserved operators Q i leads, through use of (15), to the same ground-state energy as that from the mean-field result (9). Q.E.D. This result is a generalized solution to a similar problem in the continuum limit for the open p + i p model [28].

Conclusion
In this study, we derived the ground-state energy of the generalized BCS Richardson-Gaudin model in the continuum limit through the COE. This was achieved by finding a solution to the integral equation derived from the quadratic relations of the COE. The solution corresponds to the ground state as its energy expression is consistent with a mean-field analysis. This method was originally established for the closed and open p + i p models [28] in order to circumvent the difficulties in finding a complex curve to approximate the distribution of the roots of the BAE. For future work, it would be natural to derive the BAE for the generalized model and establish the connection between the BAE and the quadratic relations of the COE. In order to find the Bethe Ansatz solution to the generalized model, the methods of [10,19,20] could be extended to derive the solution. Alternatively, an extension of the method using polynomial equations, adopted for the closed [24] and open [11] models, could be developed to derive the BAE from the quadratic equations of the COE.
Setting β = 0 in (5), the method in [11] will suffice to derive the BAE as the quadratic relations for the COE coincide with the open case. The challenge, in the case when β = 0, is to generalize the previous method using polynomials to accommodate the quadratic term in (5), and f + i f − i generalizing z 2 i . Once the BAE are derived, the Bethe roots can be studied numerically and analytically. Furthermore, insights into topological properties of this model can be derived from both the mean-field analysis following methods adopted in [13], and the Bethe Ansatz solution following [14,15].

Hence we have
Now we consider the following eigenvalue problem, leading to Minimising the energy of the Hamiltonian, we choose The ground-state energy is given as Gz 2 i 2 , and the ground-state reads Let Now we apply the Hellmann-Feynman theorem, leading to the equation Next, leading to (10). Equation (11) can be derived similarly. It is obvious to see that the expressions for ∆ x in (10) and ∆ y in (11) lead to (29). The mean-field results give the expectation value of the conserved operator (4) as This expression is consistent with that of equation (16) in Proposition 1.