Ground-state energies of the open and closed $p+ip$-pairing models from the Bethe Ansatz

Using the exact Bethe Ansatz solution, we investigate methods for calculating the ground-state energy for the $p + ip$-pairing Hamiltonian. We first consider the Hamiltonian isolated from its environment (closed model) through two forms of Bethe Ansatz solutions, which generally have complex-valued Bethe roots. A continuum limit approximation, leading to an integral equation, is applied to compute the ground-state energy. We discuss the evolution of the root distribution curve with respect to a range of parameters, and the limitations of this method. We then consider an alternative approach that transforms the Bethe Ansatz equations to an equivalent form, but in terms of the real-valued conserved operator eigenvalues. An integral equation is established for the transformed solution. This equation is shown to admit an exact solution associated with the ground state. Next we discuss results for a recently derived Bethe Ansatz solution of the open model. With the aforementioned alternative approach based on real-valued roots, combined with mean-field analysis, we are able to establish an integral equation with an exact solution that corresponds to the ground-state for this case.


Introduction
The p + ip-pairing Hamiltonian is an example of a Bardeen-Cooper-Schrieffer (BCS) model which admits an exact Bethe Ansatz solution. This result was initially established for the closed system which conserves particle number [1,2,3,4,5,6]. For the open model there is no conservation of total particle number, due to interaction terms which accommodate particle exchange with the system's environment. Consequently a u(1) symmetry is broken, which generally renders the analysis of the system to be more complicated.
Integrability of the open model was established in [7] through use of the Boundary Quantum Inverse Scattering Method. An alternative derivation, which is less technical, was later provided in [8]. Topological properties of the open model, relating to zero-energy excitations, have been studied in [9].
Like its better known ancestor the Richardson model [10], which is associated with s-wave pairing, the existence of the exact solution for the closed p + ip system provides a means to calculate the ground-state energy through a continuum limit approximation. Numerical studies suggest that, in the limit of infinite particle number, the ground-state roots become dense and lie on a connected curve in the complex plane. For Richardson's solution at half-filling this approach was first investigated by Gaudin [11], and subsequently re-examined by Román et al. [12]. One way to view this problem is to use the language of a two-dimensional electrostatic analogy. It was found in this manner that the solution for the ground-state energy coincides with the prediction coming from mean-field calculations. It is worthy of mention that an extended discussion and application of the electrostatic analogy for more general Richardson-Gaudin systems can be found in [13].
Following this approach, the continuum limit approximation has also been adopted in [3,4] for the closed p + ip-pairing model. Despite having much more complicated patterns of Bethe root distribution compared to the Richardson model, the results were again found to give agreement with mean-field analysis. However, there are some technical issues concerning the assumptions made regarding the Bethe root distributions which warrant closer scrutiny. This is the first of the primary objectives of the current study. In particular a specific example will be provided which establishes that, for certain model parameters, the continuum limit approach fails to provide fully consistent equations describing the Bethe root behaviour. However, a surprising outcome is that although the method of calculation is flawed in some instances, the results remain valid. That is, in the continuum limit the ground-state energy is the same as that predicted by mean-field theory, across all values of the model parameters. This will be proved by exploiting a completely different approach which does not use the Bethe root distribution at all.
The origins of the new approach that will be followed trace back to the work of Babelon and Talalaev [14] who showed that, through a change of variables, the Bethe Ansatz equations for Richardson-Gaudin type systems could be recast into a set of coupled polynomial equations. The roots of these equations are related to the eigenvalues of the self-adjoint conserved operators, and as such that are necessarily real-valued. The same form of polynomial equations were adopted in [15] as a means of efficient numerical solution of the conserved operator spectrum, and in [16] to compute wavefunction overlaps. In these instances the equations are quadratic. Extensions were given in [17,18] to a setting suitable for the p + ip Hamiltonian, for which the polynomial equations are also quadratic. Here it will be shown that in this form, the continuum limit approach can be formulated and solved in such a way that it does not require an Ansatz for the distribution of the roots of the equations.
The second primary objective is to apply this methodology for the calculation of the ground-state energy in the case of the open p + ip Hamiltonian. Here it will be shown how the alternative method developed to compute the ground-state energy in the closed case easily extends to the open case. It will also be shown that the result is again in complete agreement with mean-field calculations.
The general form of the integrable Hamiltonian is introduced in Sect. 2. The closed model, for which the coupling constant of the environment interaction is set to zero, is then described in detail in Sect. 3. Two forms of Bethe Ansatz solution are presented, and the continuum limit approximation for calculating the ground-state energy is formulated. Following from this an analysis exposing the limitations of the continuum limit approximation is conducted. In Sect. 4, attention turns towards formulating an alternative approach, based on the set of transformed Bethe Ansatz equations which result in coupled quadratic equations. Then in Sect. 5 the process is extended to accommodate the open model. Concluding remarks and discussion are offered in Sect. 6.

The Hamiltonian
The annihilation and creation operators for two-dimensional fermions of mass m with momentum k = k x + ik y are denoted by c k , c † k , satisfying We consider the following Hamiltonian of the pairing model interacting with its environment [9] where G, Γ are positive real constants. The sum of momenta is taken over an index set K with the properties (i) if k ∈ K, then −k ∈ K; (ii) for all k ∈ K we have |k| ≤ ω, where ω is called the cut-off. The cardinality of K is denoted as 2L. The following equality is satisfied on this Hilbert subspace Let k x + ik y = |k| exp(iφ k ), we then introduce the following notation: These operators satisfy the su(2) algebra commutation relations From now on, we use integers to enumerate the pairs of momentum states k and −k. Setting the mass to be m = 1 and z k = |k|, we rewrite (1) as Defining the following operators it can be shown that and {T j } is a set of mutually commuting conserved operators. These operators have been shown [9] to satisfy the following quadratic identities Hence the eigenvalues {t j } corresponding to {T j } give the energy expression and due to (5) {t j } satisfy 3. The p + ip model isolated from the environment The p + ip Hamiltonian is isolated from the environment (closed model) when Γ = 0. In this case, we adopt the letters G, E and T instead of G, E and T . The Hamiltonian reads The quadratic identity (5) in this case becomes The eigenvalues {t j } of {T j } in (9) then according to (10) satisfy Also note that in this case we have

First form of Bethe Ansatz solution
The Bethe Ansatz solution for (8) was obtained in [1]. We state the solution and its connection to the {t j } in (11): of the coupled Bethe Ansatz equations where M is the quantum number of particle-pairs, for each solution {y k } known as the Bethe roots, there exists a correspondence between {y k } and {t j } given by the following The techniques involved in achieving this connection (13) first appeared in [14] and are also adopted in [15,16,17,18] later. The corresponding energy for (8) is given by Each eigenstate has the form where |0 denotes the vacuum state and

Second form of Bethe Ansatz solution
Alternatively, a second form of Bethe Ansatz solution can be derived from the hole-pair perspective [6] for the isolated case. Let P = L − M denote the quantum number of hole-pairs. For each solution {v k } of the coupled equations there exists a correspondence between {v k } and {t j } given by the following The corresponding energy is given by Each eigenstate has the form where |χ denotes the completely filled state of L particle-pairs and Figure 1: Region I with x > 1 − g −1 is known as the weak coupling BCS phase. The boundary between I and II, i.e. x = 1 − g −1 , is known as the Moore-Read line. Region II with (1 − g −1 )/2 < x < 1 − g −1 and g −1 > 0 is known as the weak pairing phase. The boundary between II and III, i.e. x = (1 − g −1 )/2 is known as the Read-Green line. Region III with x < (1 − g −1 )/2 and g −1 > 0 is known as the strong pairing phase. The properties of the model for g < 0 have attracted little attention.

Symmetries of Bethe Ansatz solutions
Introduce the following parameters, Apart from this correspondence, there exists another type of relation which we call inversion. Inversion establishes an invertible mapping, given by a skewed reflection against the line x = 0.5 − g −1 , between solution sets in regions I and V I, and between solution sets in regions II and V . Regions IV and III are stable under inversion. Consider the second form of Bethe Ansatz equations (14). By setting v k = u −1 k we can derive the following expression, Note that under inversion, {z l } is no longer preserved. In other words, knowing a solution {v k } to (14) with a certain set of parameters {x, g −1 , z l }, we also obtain a solution When we combine the rotational symmetry and inversion, we achieve a correspondence between solution sets to (12) and (14). For instance, a solution {y k } to (12) (12) and (14) are invariant under complex conjugation. This implies that, in the absence of degeneracy in the spectrum of the set of conserved operators, every solution set of Bethe roots consists of complex-conjugate pairs or real numbers.

Integral approximation for the first form of Bethe Ansatz solution
Numerical solution for the ground-state Bethe roots {y j } in (12) and its peculiar behaviour under certain choices of parameters are discussed in [1] and [5]. The distribution of the Bethe roots suggests that they lie on curves in the complex plane, allowing an integral approximation to be applied. The continuum limit approximation for (12) where L is large is studied in [3,4]. In the limit, we require that M also becomes large while x is finite, and similarly G becomes small while g is finite.
First we formally define the discrete density for each single particle energy level z 2 k . Since all the z 2 k are real and positive it is natural to relabel them as z 2 k with z 2 k < z 2 j whenever k < j such that z 2 1 is the smallest and z 2 L is the largest. The discrete root densityρ is defined as In the continuum limit, we introduce ρ to be the continuum density for z 2 k with connected support being a subset of (0, ω) and replace all z 2 k with a continuous variable . the continuum approximation for summation over any given function f is undertaken by replacing sum with integral according to Setting f = 1 in (16) gives the normalization condition for the density, We also introduce a continuous curve Ω, which is invariant under complex conjugation, to approximate the distribution of the ground-state Bethe roots {y k } in (12). Let r(y) be the density for {y k } in the continuum limit with support on Ω. Since on the RHS of (12), the expression gives rise to a singularity, we adopt the Cauchy principal value to approximate the summation where − denotes the Cauchy principal value of an integral. The continuum limit approximation for the Bethe Ansatz solution {y k } for (12) reads as where equation (19) is the normalization condition for the density r(y). The ground-state energy is given by The solution curve Ω and r(y) have been solved in [3] and are dependent on the choice of the parameters x, g −1 as shown in Fig. 1. The solution  curve Ω consists of two parts, a complex part Ω C depicted by a solid line and a real part Ω A depicted by one or more dashed lines, see Fig. 2. In Fig. 2(a), the complex arc Ω C intersects the real part Ω A = (0, A ) at A . In Fig. 2(b), the arc Ω C closes with an endpoint b on the negative real line as (g −1 , x) arrives at the Moore-Read line from right. In Fig. 2(c), a line segment (a, b) forms on the negative real line adding a component to Ω A . Then the point b will approach 0 with the complex curve Ω C shrinking until it vanishes when b = A = 0 and (g −1 , x) arrives at the Read-Green line from right. In Fig. 2(d), as (g −1 , x) departs from the Read-Green line to the left, b becomes negative and Ω consists of one real part Ω A = (a, b).

Approximation for the second form of Bethe Ansatz solution
Following the approach as discussed in previous Sect. 3.4, the continuum limit approximation for the ground-state Bethe Ansatz solution {v k } of (14) reads as Re y where Ω is a continuous curve introduced to approximate the distribution of the Bethe roots {v k }, and r(v) is the density for {v k } in the continuum limit satisfying Ω |dv| r(v) = P.
The ground-state energy is given by Now we need to solve for Ω and density r(v) for the ground state. The solution curve Ω is proposed to be classified under the four phases as of the approximation for the first form Bethe Ansatz solution (18), see Fig. 3. Here we have modified shapes deduced from the rotational symmetry combined with inversion discussed in Sect. 3.3. Fig. 3(a)(b)(c) are topological inversion of Fig. 2(a)(b)(c). In Fig. 3(c), the point a approaches zero as (g −1 , x) approaches the Read-Green line from right. As (g −1 , x) departs from the Read-Green line and continues to move left, a third real line segment (0, B ) and a second complex loop intersecting the real line at a and B appear, see Fig. 3(d). Further calculation gives rise to the same energy expression as from the integral approximation (18) and mean-field analysis [3].
3.6. Limitations of the continuum limit approximation for the Bethe Ansatz solutions The Moore-Read line is an example of a ground-state phase boundary line associated with changes in the topology of the root distribution. In previous Sect. 3.4 and 3.5, as we send parameters (g −1 , x) from the weak coupling BCS phase to the Moore-Read line, the complex part of the solution curve Ω C evolves until it closes and forms a loop. However in the discrete case, all the Bethe roots condense at the origin [1] when the parameters reach the Moore-Read line. This discrepancy between the integral approximation and the discrete case is not present in models such as the Richardson swave pairing model and the d + id-wave pairing Hamiltonian [19,20] where a similar approach of integral approximation is adopted. In the case of the twolevel Richardson s-wave pairing model [19], the solution curve of the integral approximation evolves until it closes and forms a loop as some governing parameters approach a ground-state phase boundary line, meanwhile the discrete Bethe roots do not condense and their distribution is predicted by the solution curve within small error in the limit. In the case of the d+id-wave pairing Hamiltonian [19,20], the solution curve of the integral approximation contracts to a point at the origin and the discrete Bethe roots also condense at the origin as the governing parameters approach a ground-state phase boundary line.
Due to the aforementioned discrepancy in the closed model, we perform a closer inspection of the integral approximation (18) in the Moore-Read line case with solution curve Ω depicted in Fig. 2(b). The general form of the density r(y) for this case is proposed in [3] to be the following, The arc Ω C is obtained by solving the following integral equation where y 0 is any point of Ω C . We first consider the limiting case where G → 0. In this case the Bethe roots {y k } are all real and lie within the interval (0, ω f ] where ω f is the upper bound for {y k } and ω f < ω. Hence in the integral approximation the solution curve Ω reduces to Ω A = (0, ω f ) and A = ω f . According to (21), the normalization condition (19) then reads as which determines the value of ω f . As G increases, the complex part Ω C starts to form and A decreases away from ω f . Consequently the allowable bound for A is (0, ω f ). Now we consider the special case of an inverse-square density This is a limiting case for the density ρ( ) as we let it vanish on (0, ω 0 ) while sending ω → +∞. When G → 0, substituting the inverse-square density into (23) we have ω f = ω 0 /(1 − x). Hence the constraint for A is The equation for Ω C is then derived from (22), the arc Ω C is determined by the following equation Choosing y 0 to be a, Ω C is given by (25) The remaining constraint s(a) = 0 from (21) implies that Setting x = 0.4, from (26) we numerically determine that −a/ω 0 ≈ 1.58. From (25), the solution curve for y/ω 0 is plotted as in Fig. (4). However, with this choice of parameters and inverse-square density, numerical results show that the intersection A is outside the allowable bound (ω 0 , 5ω 0 /3) from (24). This inconsistency is verified when we study the solution curve in the integral approximation for the second form of Bethe Ansatz solution (20) in the Moore-Read line case with a constant density, which corresponds to an inverse-square density for (18) in the Moore-Read line case as is discussed in Sect. 3.3. The resulting solution curve again intersects the real line at positions outside the allowable bound.
Alternatively, if we choose a constant density ρ( ) = L/(ω −ω 0 ) with support on (ω 0 , ω) for (18), it can be shown that for the integral approximation to be valid ω 0 /ω cannot exceed the numerically determined value 0.04246. These results suggest that the validity of the continuum limit approximation is dependent on the choice of the density function ρ( ).

Conserved operator eigenvalue method
The difficulties in the integral approximation for the Bethe roots {y k } arise from the task to find a suitable solution curve Ω and solve for its density r(y) in (18). The next step is to adopt an alternative approach that eliminates these requirements and accommodates to an arbitrary form of density ρ( ) subject to (17). In the following discussion, we continue with the approach that first appeared in [14]. Defining Since each solution {y k } to (12) consists of complex-conjugate pairs or real numbers, Λ j and t j are all real. Substituting into (11) we obtain the following quadratic equations [17,18], where 2q = G −1 + 2M − L − 1. A continuum limit approximation applied to (28) leads to the following integral equation, The objective now is to derive a solution to (29). We will establish several useful results before the solution is stated. Let Note that a, b are determined by G and q. We require them to be a complexconjugate pair, or both negative real numbers, i.e.
a =b ∈ C \ R or a ≤ b ≤ 0.
As a result a + b ∈ R and ab ≥ 0. The function R( ) is the elementary real-valued square-root with ∈ (0, ω). In the limit where q = 0, we require b to vanish also. In this case (30) becomes where a is real and determined by G. All the following results and proofs leading to the expression for the ground-state energy for the closed model assume q = 0 = b. In the special case where q = b = 0, simply replace (2|q|/ √ ab) with the finite value C in (32), the proof of which follows a similar calculation and is omitted.

Corollary 2. Let
Proof: Exchanging the variables δ and γ in the integrand, and then performing a change of order of integration yields where the second step is due to Lemma 1. Now we add up two distinct expressions for I/2, Since we conclude that Proof: Since again by adding up two distinct expressions for K/2, we have Finally, following a similar calculation as that for K, Now we prove the following: Proposition 4. The following function is a solution to the integral equation (29) with a, b subject to (30) and (31). Proof: On the RHS of (29), we first use Lemma 1 and simplify the following term Finally, The integral approximation for (27) corresponding to the solution Λ( ) given in (34) is As a direct consequence of Proposition 4, we have the following result, Corollary 5. The function t( ) defined in (35) is a solution to the integral approximation of (11), This particular solution t( ) corresponds to the ground state since from (6) in the continuum limit, the last step is due to expression P in Lemma 3. The energy expression (37) coincides with the ground-state energy derived by mean-field analysis [3] and integral approximation of both forms of Bethe Ansatz solutions discussed in Sect. 3.4 and 3.5 despite its limitations.

The p + ip model interacting with its environment
Consider the open model (2) with extra terms governed by parameter Γ,

Bethe Ansatz equations and numerics
The Bethe Ansatz solution for (2) was derived in [7]: for each solution {v j } of the coupled equations where α = 1 + G −1 , β = Γ/G and j = 1, . . . , L, there is a correspondence between {t j } in (7) and {v j } via a change of variables This translates to the fact that (7) and (38) are equivalent. It is worth mentioning that the difference between the closed and open model is that in the closed model, it is possible for the number of Bethe roots M to be less than L, while in the open model we must have exactly L Bethe roots [8,21]. The energy is given by (6), The corresponding eigenstate reads as [9] L j=1 We perform some numerical analysis to study the behaviour of the Bethe roots {v j }. We consider a small-sized case where L = 5. As we send Γ from 0 to a large number, numerical results (see Tab. 1) show that one of {v j } in (38) diverges around parameters (G = 1, Γ 2 = 10.213). Also at certain values of (G, Γ), we have two real roots meeting to form a complex-conjugate pair and vice versa, see Fig. 5, 6 and 7.

Conserved operator eigenvalue method for the open model
The numerical analysis of the distribution of the Bethe roots gives no clear indications of a solution curve for a large particle number. However, since (38) is invariant under complex conjugation, again assuming the absence of degeneracy in the spectrum of the set of conserved operators, each solution set {v j } consists of complex-conjugate pairs or real numbers. Hence from (39), it is clear to see that all solution sets of {t j } are real. We now consider the integral approximation for (7), which reads as where F = −2Γ. Our task is to find a solution to (41) that corresponds to the ground state. We resort to mean-field analysis for suggestions of a possible solution. The mean-field analysis for the open model and its results  are included in Appendix A. We consider the operators T j defined in (3). From (55) and (56), we derive the mean-field expression for T j as where The integral approximation for (42) is We will establish that (43) is in fact a solution to (41), i.e.
with the "gap" equation determining the value of ∆, In addition, this solution t( ) corresponds to the ground state. We again establish some useful results to assist our proof.
Lemma 6. Given equation (44), let Proof: In (32), set then we have R( ) = R( ) and C = C. Hence from (35) and by Corollary 5, T ( ) = t( ) satisfying (36), Finally, Now we show the following: The following function is a solution to the integral equation (41) with ∆ subject to (44). Proof: We adopt the notation introduced in Lemma 6, rewriting then from (45) we are able to simplify equation (41) as and we use the special case of q = b = 0 for expression O in Lemma 3, where (2|q|/ √ ab), (a+b) and G −1 are replaced by C, −(∆+F ) 2 and ∆G −1 (∆+F ) −1 respectively, consequently we have Now the LHS of (49) is reduced to the following We derive the energy expression using t( ) from (6), Again from the special case for expression P in Lemma 3, where b, a and G −1 are replaced by 0, −(∆ + F ) 2 and ∆G −1 (∆ + F ) −1 respectively, we have This energy expression is consistent with (50) (see Appendix A) for the ground state, which is derived from mean-field analysis.

Conclusion
The continuum limit approximation for calculating the ground-state energy of the p+ip model, via the Bethe Ansatz solution, was studied. Starting with the closed model, we revisited the formulations of [3,4,13] which undertake calculations by assuming a form of density function for the Bethe root distribution. It was found that this approach does not provide a consistent solution for particular choices of the momentum density distribution. This was established by a close examination of the case known as the Moore-Read line, where it is known that all the ground-state Bethe roots collapse at the origin [1,3,4,5]. An alternative approach, which avoids the need to postulate a form for the Bethe root density, was proposed in terms of the coupled equations satisfied by the conserved operator eigenvalues [17,18]. In this case a solution corresponding to the ground state in the continuum limit was found. Curiously, the expression obtained for the ground-state energy coincides with that obtained by the Bethe root distribution. That is to say that although the Bethe root approach involves a flawed methodology, it nonetheless produces the correct answer! In both cases the ground-state energy per particle is exactly the same as the mean-field prediction in the limit of infinite particle number.
The conserved operator eigenvalue approach was then extended to accommodate the open case based on results from [8,9,21]. Again, it was found that the result for the ground-state energy is in agreement with mean-field calculations. The open case can be considered as a model allowing for the exchange of particle between the system and the environment. It is important to note that in this case there is no signature of any quantum phase transition, no matter how weak the environment coupling is, which is in stark contrast to the closed system. A similar scenario was considered in [22] where the environment was modelled by a single bosonic degree of freedom in such a way that integrability was preserved. There too, arbitrarily small coupling to the environment was found to annihilate the existence of any quantum phase transition.
For future work it would be natural to extend this analysis to calculate the leading order finite-size correction to the ground-state energy, for both open and closed models. For the closed s-wave pairing systems there are results obtained by series expansions [23], which can be extrapolated and shown to be valid more generally [24]. Obtaining analogous expressions for the open and closed p + ip-pairing systems appears to be entirely feasible.
Then let∆ = 2G Q ,∆ † = 2G Q † and ∆ = |∆|. The extended Hamiltonian (2) can be rewritten as Note that in the mean-field approximation for this extended model, the Lagrange multiplier is µ = 0. The matrix form is derived from the representation of H acting on (C 2 ) ⊗L . Consider the following eigenvalue problem, by minimizing each eigenvalue, we derive the ground-state energy where χ = |χ|,χ =∆ + F . Then we calculate the following Apply Hellmann-Feynman theorem, since and Furthermore from (51), Comparing (52) and (53), and assuming F = 0, we conclude that i.e.∆ =∆ † = ∆. It follows that χ =χ =χ † = ∆+F . Hence we immediately have the following results (57)