Ground-state energies of the open and closed p+ip -pairing models from the Bethe Ansatz
GGround-state energies of the open and closed p + ip -pairing models from the Bethe Ansatz Yibing Shen, Phillip S. Isaac, Jon Links ∗ School of Mathematics and Physics,The University of Queensland, Brisbane, QLD 4072, Australia ∗ email: [email protected] Abstract
Using the exact Bethe Ansatz solution, we investigate methods for calculatingthe ground-state energy for the p + ip -pairing Hamiltonian. We first considerthe Hamiltonian isolated from its environment (closed model) through twoforms of Bethe Ansatz solutions, which generally have complex-valued Betheroots. A continuum limit approximation, leading to an integral equation,is applied to compute the ground-state energy. We discuss the evolution ofthe root distribution curve with respect to a range of parameters, and thelimitations of this method. We then consider an alternative approach thattransforms the Bethe Ansatz equations to an equivalent form, but in termsof the real-valued conserved operator eigenvalues. An integral equation isestablished for the transformed solution. This equation is shown to admitan exact solution associated with the ground state. Next we discuss resultsfor a recently derived Bethe Ansatz solution of the open model. With theaforementioned alternative approach based on real-valued roots, combinedwith mean-field analysis, we are able to establish an integral equation withan exact solution that corresponds to the ground-state for this case. Keywords:
Integrable systems, BCS model, Bethe Ansatz, Bethe rootdistributions
Preprint submitted to Elsevier July 3, 2018 a r X i v : . [ n li n . S I] J u l . Introduction The p + ip -pairing Hamiltonian is an example of a Bardeen-Cooper-Schrieffer (BCS) model which admits an exact Bethe Ansatz solution. Thisresult was initially established for the closed system which conserves particlenumber [1, 2, 3, 4, 5, 6]. For the open model there is no conservation oftotal particle number, due to interaction terms which accommodate particleexchange with the system’s environment. Consequently a u (1) symmetry isbroken, which generally renders the analysis of the system to be more com-plicated. Integrability of the open model was established in [7] through use ofthe Boundary Quantum Inverse Scattering Method. An alternative deriva-tion, which is less technical, was later provided in [8]. Topological propertiesof the open model, relating to zero-energy excitations, have been studied in[9]. Like its better known ancestor the Richardson model [10], which is asso-ciated 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 througha continuum limit approximation. Numerical studies suggest that, in thelimit of infinite particle number, the ground-state roots become dense andlie on a connected curve in the complex plane. For Richardson’s solution athalf-filling this approach was first investigated by Gaudin [11], and subse-quently re-examined by Rom´an et al. [12]. One way to view this problem isto use the language of a two-dimensional electrostatic analogy. It was foundin this manner that the solution for the ground-state energy coincides withthe prediction coming from mean-field calculations. It is worthy of mentionthat an extended discussion and application of the electrostatic analogy formore general Richardson-Gaudin systems can be found in [13].Following this approach, the continuum limit approximation has alsobeen adopted in [3, 4] for the closed p + ip -pairing model. Despite hav-ing much more complicated patterns of Bethe root distribution compared tothe Richardson model, the results were again found to give agreement withmean-field analysis. However, there are some technical issues concerningthe assumptions made regarding the Bethe root distributions which warrantcloser scrutiny. This is the first of the primary objectives of the currentstudy. In particular a specific example will be provided which establishesthat, for certain model parameters, the continuum limit approach fails to pro-vide fully consistent equations describing the Bethe root behaviour. However,a surprising outcome is that although the method of calculation is flawed in2ome instances, the results remain valid. That is, in the continuum limitthe 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 acompletely different approach which does not use the Bethe root distributionat all.The origins of the new approach that will be followed trace back to thework of Babelon and Talalaev [14] who showed that, through a change ofvariables, the Bethe Ansatz equations for Richardson-Gaudin type systemscould be recast into a set of coupled polynomial equations. The roots ofthese equations are related to the eigenvalues of the self-adjoint conservedoperators, and as such that are necessarily real-valued. The same form ofpolynomial equations were adopted in [15] as a means of efficient numericalsolution of the conserved operator spectrum, and in [16] to compute wave-function overlaps. In these instances the equations are quadratic. Extensionswere given in [17, 18] to a setting suitable for the p + ip Hamiltonian, forwhich the polynomial equations are also quadratic. Here it will be shown thatin this form, the continuum limit approach can be formulated and solved insuch a way that it does not require an Ansatz for the distribution of the rootsof the equations.The second primary objective is to apply this methodology for the calcu-lation 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 theground-state energy in the closed case easily extends to the open case. It willalso be shown that the result is again in complete agreement with mean-fieldcalculations.The general form of the integrable Hamiltonian is introduced in Sect.2. The closed model, for which the coupling constant of the environmentinteraction is set to zero, is then described in detail in Sect. 3. Two forms ofBethe Ansatz solution are presented, and the continuum limit approximationfor calculating the ground-state energy is formulated. Following from this ananalysis exposing the limitations of the continuum limit approximation isconducted. In Sect. 4, attention turns towards formulating an alternativeapproach, based on the set of transformed Bethe Ansatz equations whichresult in coupled quadratic equations. Then in Sect. 5 the process is extendedto accommodate the open model. Concluding remarks and discussion areoffered in Sect. 6. 3 . The Hamiltonian
The annihilation and creation operators for two-dimensional fermions ofmass m with momentum k = k x + i k y are denoted by c k , c † k , satisfying { c k , c k (cid:48) } = { c † k , c † k (cid:48) } = 0 , { c k , c † k (cid:48) } = δ kk (cid:48) I. We consider the following Hamiltonian of the pairing model interacting withits environment [9] H = (cid:88) k | k | m c † k c k − G m (cid:88) k (cid:54) = ± k (cid:48) k ¯ k (cid:48) c † k c †− k c − k (cid:48) c k (cid:48) + Γ2 (cid:88) k ( k c † k c †− k + ¯ k c − k c k ) , (1)where G , Γ are positive real constants. The sum of momenta is taken overan 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 2 L . The following equality is satisfied on this Hilbert subspace2 c † k c k c †− k c − k = c † k c k + c †− k c − k . Let k x + i k y = | k | exp(i φ k ), we then introduce the following notation: S + k = exp ( iφ k ) c † k c †− k , S − k = exp ( − iφ k ) c − k c k , S z k = c † k c †− k c − k c k − I . These operators satisfy the su (2) algebra commutation relations[ S z k , S ± k ] = ± S ± k , [ S + k , S − k ] = 2 S z k . 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 H = L (cid:88) k =1 z k S zk − G L (cid:88) k =1 L (cid:88) j (cid:54) = k z k z j S + k S − j + Γ L (cid:88) k =1 z k ( S + k + S − k ) + (cid:32) L (cid:88) k =1 z k (cid:33) I. (2)Defining the following operators T j = L (cid:88) k (cid:54) = j (cid:18) z k z k − z j (4 S zj S zk − I ) + 2 z j z k z k − z j ( S + j S − k + S − j S + k ) (cid:19) + 2 G S zj + 2Γ G z − j ( S + j + S − j ) , (3)4t can be shown that H = G L (cid:88) j =1 z j T j + (cid:32) L (cid:88) j =1 z j (cid:33) I, (4)and {T j } is a set of mutually commuting conserved operators. These opera-tors have been shown [9] to satisfy the following quadratic identities T j = 1 G + 4Γ z − j G + 2 L (cid:88) k (cid:54) = j z k ( T j − T k ) z j − z k . (5)Hence the eigenvalues { t j } corresponding to {T j } give the energy expression E = G L (cid:88) j =1 z j t j + 12 L (cid:88) j =1 z j , (6)and due to (5) { t j } satisfy t j = 1 G + 4Γ z − j G + 2 L (cid:88) k (cid:54) = j z k ( t j − t k ) z j − z k . (7)
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 H = L (cid:88) k =1 z k S zk − G L (cid:88) k =1 L (cid:88) j (cid:54) = k z k z j S + k S − j + (cid:32) L (cid:88) k =1 z k (cid:33) I (8)= G L (cid:88) j =1 z j T j + (cid:32) L (cid:88) j =1 z j (cid:33) I, where T j = L (cid:88) k (cid:54) = j (cid:18) z k z k − z j (4 S zj S zk − I ) + 2 z j z k z k − z j ( S + j S − k + S − j S + k ) (cid:19) + 2 G S zj . (9)5he quadratic identity (5) in this case becomes T j = 1 G + 2 L (cid:88) k (cid:54) = j z k ( T j − T k ) z j − z k . (10)The eigenvalues { t j } of { T j } in (9) then according to (10) satisfy t j = 1 G + 2 L (cid:88) k (cid:54) = j z k ( t j − t k ) z j − z k , j = 1 , , . . . , L. (11)Also note that in this case we have L (cid:88) j =1 T j = 2 (cid:32) L (cid:88) k =1 S zk (cid:33) + 2 G L (cid:88) k =1 S zk − L I. The Bethe Ansatz solution for (8) was obtained in [1]. We state thesolution and its connection to the { t j } in (11): of the coupled Bethe Ansatzequations G − + 2 M − L − y k + L (cid:88) l =1 y k − z l = M (cid:88) j (cid:54) = k y k − y j , k = 1 , . . . , M, (12)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 t j = − G − − M + 2 z j M (cid:88) k =1 z j − y k . (13)The techniques involved in achieving this connection (13) first appeared in[14] and are also adopted in [15, 16, 17, 18] later. The corresponding energyfor (8) is given by E = G L (cid:88) j =1 z j t j + 12 L (cid:88) j =1 z j = (1 + G ) M (cid:88) k =1 y k . | Φ (cid:105) = M (cid:89) k =1 C ( y k ) | (cid:105) , where | (cid:105) denotes the vacuum state and C ( y ) = L (cid:88) j =1 z j y − z j S + k . Alternatively, a second form of Bethe Ansatz solution can be derived fromthe hole-pair perspective [6] for the isolated case. Let P = L − M denotethe quantum number of hole-pairs. For each solution { v k } of the coupledequations − G − + 2 P − L − v k + L (cid:88) l =1 v k − z l = P (cid:88) j (cid:54) = k v k − v j , k = 1 , . . . , P, (14)there exists a correspondence between { v k } and { t j } given by the following t j = G − − P + 2 z j P (cid:88) k =1 z j − v k . The corresponding energy is given by E = G L (cid:88) j =1 z j t j + 12 L (cid:88) j =1 z j = L (cid:88) l =1 z l + ( G − P (cid:88) k =1 v k . Each eigenstate has the form | Ψ (cid:105) = P (cid:89) k =1 B ( y k ) | χ (cid:105) , where | χ (cid:105) denotes the completely filled state of L particle-pairs and B ( y ) = L (cid:88) j =1 z j y − z j S − k . − x . − x = 1 − g − IIIIIIIVVV I
Figure 1: Region I with x > − g − is known as the weak coupling BCS phase. Theboundary between I and II , i.e. x = 1 − g − , is known as the Moore-Read line. Region II with (1 − g − ) / < x < − g − and g − > II and III , i.e. x = (1 − g − ) / III with x < (1 − g − ) / g − > g < Introduce the following parameters, x = M/L, g = GL, such that x takes values in [0 ,
1] and g in ( −∞ , ∞ ) as shown in Fig. 1. Thereexists a rotational symmetry around the point ( g − = 0 , x = 0 .
5) betweenthe two forms of Bethe Ansatz solutions. For instance, if we have a solution { y k } to (12) with x = 0 . g − = 0 .
4, then this solution corresponds to asolution { v k = y k } to (14) with x = 0 . g − = − . { z l } being fixed.Apart from this correspondence, there exists another type of relationwhich we call inversion. Inversion establishes an invertible mapping, givenby a skewed reflection against the line x = 0 . − g − , between solution setsin 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 BetheAnsatz equations (14). By setting v k = u − k we can derive the followingexpression, G − − u k + L (cid:88) l =1 u k − z − l = P (cid:88) j (cid:54) = k u k − u j , k = 1 , . . . , P. { z l } is no longer preserved. In other words, know-ing a solution { v k } to (14) with a certain set of parameters { x, g − , z l } ,we also obtain a solution { u k = v − k } to (14) with transformed parameters { x, ¯ g − = − x − g − + 1 , ¯ z l = z − l } .When we combine the rotational symmetry and inversion, we achievea correspondence between solution sets to (12) and (14). For instance, asolution { y k } to (12) with x = 0 . g − = 0 . { z l } is also asolution { v k = y k } to (14) with x = 0 . g − = − . { z l } . Then the inversion of this solution { u k = v − k } also solves (14) with x = 0 . g − = 0 . { z − l } . Moreover, since G ,Γ and z k are all positive real constants, the Bethe Ansatz equations (12)and (14) are invariant under complex conjugation. This implies that, in theabsence of degeneracy in the spectrum of the set of conserved operators,every solution set of Bethe roots consists of complex-conjugate pairs or realnumbers. Numerical solution for the ground-state Bethe roots { y j } in (12) and itspeculiar behaviour under certain choices of parameters are discussed in [1]and [5]. The distribution of the Bethe roots suggests that they lie on curvesin the complex plane, allowing an integral approximation to be applied. Thecontinuum 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, andsimilarly G becomes small while g is finite.First we formally define the discrete density for each single particle energylevel z k . Since all the z k are real and positive it is natural to relabel themas z k with z k < z j whenever k < j such that z is the smallest and z L is thelargest. The discrete root density ˜ ρ is defined as˜ ρ ( z j ) = L ( L − z j +1 − z j ) , ˜ ρ ( z L ) = 0 , (15)such that L (cid:88) j =1 ˜ ρ ( z j ) · ( z j +1 − z j ) = L. In the continuum limit, we introduce ρ to be the continuum density for z k with connected support being a subset of (0 , ω ) and replace all z k with a9ontinuous variable (cid:15) . the continuum approximation for summation over anygiven function f is undertaken by replacing sum with integral according to L (cid:88) j =1 f ( z j ) −→ (cid:90) ω d (cid:15) ρ ( (cid:15) ) f ( (cid:15) ) . (16)Setting f = 1 in (16) gives the normalization condition for the density, (cid:90) ω d (cid:15) ρ ( (cid:15) ) = L. (17)We also introduce a continuous curve Ω, which is invariant under complexconjugation, 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 withsupport on Ω. Since on the RHS of (12), the expression f ( y j ) = 2 y k − y j gives rise to a singularity, we adopt the Cauchy principal value to approxi-mate the summation M (cid:88) j (cid:54) = k f ( y j ) −→ − (cid:90) Ω | d y | r ( y ) f ( y ) , where − (cid:82) denotes the Cauchy principal value of an integral. The continuumlimit approximation for the Bethe Ansatz solution { y k } for (12) reads as (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) − y − G − + 2 M − Ly = − (cid:90) Ω | d y (cid:48) | r ( y (cid:48) ) y (cid:48) − y , (18) (cid:90) Ω | d y | r ( y ) = M, (19)where equation (19) is the normalization condition for the density r ( y ). Theground-state energy is given by E = (cid:90) Ω | d y | yr ( y ) . The solution curve Ω and r ( y ) have been solved in [3] and are dependenton the choice of the parameters x , g − as shown in Fig. 1. The solution10 e y Im y (cid:15) A ω (a) Weak coupling BCS Re y Im y (cid:15) A ωb (b) Moore-Read line Re y Im y (cid:15) A ωba (c) Weak pairing Re y Im y ωba (d) Strong pairing Figure 2: The solution curve Ω evolves from (a) to (d) as we send g from 0 to + ∞ while x is fixed at a value between 0 and 0 . curve Ω consists of two parts, a complex part Ω C depicted by a solid lineand a real part Ω A depicted by one or more dashed lines, see Fig. 2. InFig. 2(a), the complex arc Ω C intersects the real part Ω A = (0 , (cid:15) A ) at (cid:15) A .In Fig. 2(b), the arc Ω C closes with an endpoint b on the negative real lineas ( g − , x ) arrives at the Moore-Read line from right. In Fig. 2(c), a linesegment ( 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 untilit vanishes when b = (cid:15) A = 0 and ( g − , x ) arrives at the Read-Green line fromright. In Fig. 2(d), as ( g − , x ) departs from the Read-Green line to the left, b becomes negative and Ω consists of one real part Ω A = ( a, b ). Following the approach as discussed in previous Sect. 3.4, the continuumlimit approximation for the ground-state Bethe Ansatz solution { v k } of (14)reads as (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) − v − − G − + 2 P − Lv = − (cid:90) Ω | d v (cid:48) | r ( v (cid:48) ) v (cid:48) − v , (20)11 e y Im y (cid:15) A ωab (a) Weak coupling BCS Re y Im y (cid:15) A ωb (b) Moore-Read line Re y Im y (cid:15) A ωb a (c) Weak pairing Re y Im y (cid:15) A ωb a (cid:15) B (d) Strong pairing Figure 3: The solution curve Ω proposed for the integral approximation for the secondform Bethe Ansatz solution consists of a real part Ω A depicted by one or more dashedlines and a complex part Ω C depicted by one or more solid lines. Again we fix x ∈ (0 , . g from 0 to + ∞ . where Ω is a continuous curve introduced to approximate the distribution ofthe Bethe roots { v k } , and r ( v ) is the density for { v k } in the continuum limitsatisfying (cid:90) Ω | d v | r ( v ) = P. The ground-state energy is given by E = (cid:90) ω d (cid:15) (cid:15)ρ ( (cid:15) ) − (cid:90) Ω | d v | vr ( v ) . Now we need to solve for Ω and density r ( v ) for the ground state. Thesolution curve Ω is proposed to be classified under the four phases as of theapproximation for the first form Bethe Ansatz solution (18), see Fig. 3. Herewe have modified shapes deduced from the rotational symmetry combinedwith inversion discussed in Sect. 3.3. Fig. 3(a)(b)(c) are topological inver-sion of Fig. 2(a)(b)(c). In Fig. 3(c), the point a approaches zero as ( g − , x )approaches the Read-Green line from right. As ( g − , x ) departs from theRead-Green line and continues to move left, a third real line segment (0 , (cid:15) B )12nd a second complex loop intersecting the real line at a and (cid:15) B appear, seeFig. 3(d). Further calculation gives rise to the same energy expression asfrom the integral approximation (18) and mean-field analysis [3]. The Moore-Read line is an example of a ground-state phase boundary lineassociated with changes in the topology of the root distribution. In previousSect. 3.4 and 3.5, as we send parameters ( g − , x ) from the weak couplingBCS 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 reachthe Moore-Read line. This discrepancy between the integral approximationand the discrete case is not present in models such as the Richardson s -wave pairing model and the d + id -wave pairing Hamiltonian [19, 20] where asimilar approach of integral approximation is adopted. In the case of the two-level Richardson s -wave pairing model [19], the solution curve of the integralapproximation evolves until it closes and forms a loop as some governingparameters approach a ground-state phase boundary line, meanwhile thediscrete Bethe roots do not condense and their distribution is predicted bythe solution curve within small error in the limit. In the case of the d + id -wavepairing Hamiltonian [19, 20], the solution curve of the integral approximationcontracts to a point at the origin and the discrete Bethe roots also condenseat the origin as the governing parameters approach a ground-state phaseboundary line.Due to the aforementioned discrepancy in the closed model, we performa closer inspection of the integral approximation (18) in the Moore-Read linecase with solution curve Ω depicted in Fig. 2(b). The general form of thedensity r ( y ) for this case is proposed in [3] to be the following, r ( y ) | d y | = (cid:40) s ( y ) d y if y ∈ Ω C ρ ( y ) d y if y ∈ Ω A , (21) s ( y ) = 12 iπ (cid:20)(cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) − y + My (cid:21) ,s ( a ) = 0 . C is obtained by solving the following integral equationIm (cid:20)(cid:90) yy d y (cid:48) s ( y (cid:48) ) (cid:21) = 0 , y ∈ Ω C , (22)where y is any point of Ω C .We first consider the limiting case where G →
0. In this case the Betheroots { y k } are all real and lie within the interval (0 , ω f ] where ω f is theupper bound for { y k } and ω f < ω . Hence in the integral approximation thesolution curve Ω reduces to Ω A = (0 , ω f ) and (cid:15) A = ω f . According to (21),the normalization condition (19) then reads as (cid:90) Ω | d y | r ( y ) = (cid:90) ω f d y ρ ( y ) = M, (23)which determines the value of ω f . As G increases, the complex part Ω C startsto form and (cid:15) A decreases away from ω f . Consequently the allowable boundfor (cid:15) A is (0 , ω f ).Now we consider the special case of an inverse-square density ρ ( (cid:15) ) = ω L(cid:15) , (cid:15) ∈ ( ω , + ∞ ) . This is a limiting case for the density ρ ( (cid:15) ) as we let it vanish on (0 , ω ) whilesending ω → + ∞ . When G →
0, substituting the inverse-square densityinto (23) we have ω f = ω / (1 − x ). Hence the constraint for (cid:15) A is (cid:15) A ∈ (cid:18) ω , ω − x (cid:19) . (24)The equation for Ω C is then derived from (22),Re (cid:20) ω L (cid:90) ∞ ω d (cid:15) (cid:15) log (cid:18) (cid:15) − y(cid:15) − y (cid:19) − M log yy (cid:21) = 0 . Since ω (cid:90) ∞ ω d (cid:15) (cid:15) log (cid:18) (cid:15) − y(cid:15) − y (cid:19) = (cid:18) − ω y (cid:19) log (cid:18) − yω (cid:19) − (cid:18) − ω y (cid:19) log (cid:18) − y ω (cid:19) , C is determined by the following equationRe (cid:20)(cid:18) − ω y (cid:19) log (cid:18) − yω (cid:19) − (cid:18) − ω y (cid:19) log (cid:18) − y ω (cid:19) − x log yy (cid:21) = 0 . Choosing y to be a , Ω C is given byRe (cid:20)(cid:18) − ω y (cid:19) log (cid:18) − yω (cid:19) − (cid:16) − ω a (cid:17) log (cid:18) − aω (cid:19) − x log ya (cid:21) = 0 . (25)The remaining constraint s ( a ) = 0 from (21) implies that ω L (cid:90) ∞ ω (cid:15) − a ) (cid:15) + Ma = 0 , hence we have exp (cid:18) (1 − x ) − aω (cid:19) = 1 + − aω . (26)Setting x = 0 .
4, from (26) we numerically determine that − a/ω ≈ . y/ω is plotted as in Fig. (4). (a) (b) Figure 4: In (a), equation (25) is plotted with z = y/ω and the filling fraction x = 0 . C . Note that in our setting ω < (cid:15) A < ω f = 5 ω /
3, hence 1 < (cid:15) A /ω < /
3. Numerical calculation shows that (cid:15) A /ω ≈ .
71 or 0 . However, with this choice of parameters and inverse-square density, nu-merical results show that the intersection (cid:15) A is outside the allowable bound15 ω , ω /
3) from (24). This inconsistency is verified when we study the solu-tion curve in the integral approximation for the second form of Bethe Ansatzsolution (20) in the Moore-Read line case with a constant density, which cor-responds to an inverse-square density for (18) in the Moore-Read line caseas is discussed in Sect. 3.3. The resulting solution curve again intersects thereal line at positions outside the allowable bound.Alternatively, if we choose a constant density ρ ( (cid:15) ) = L/ ( ω − ω ) with sup-port on ( ω , ω ) for (18), it can be shown that for the integral approximationto be valid ω /ω cannot exceed the numerically determined value 0 . ρ ( (cid:15) ).
4. 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 eliminatesthese requirements and accommodates to an arbitrary form of density ρ ( (cid:15) )subject to (17). In the following discussion, we continue with the approachthat first appeared in [14]. DefiningΛ j = M (cid:88) k =1 z j − y k , from (13) we have t j = − G − − M + 2 z j Λ j . (27)Since each solution { y k } to (12) consists of complex-conjugate pairs or realnumbers, Λ j and t j are all real. Substituting into (11) we obtain the followingquadratic equations [17, 18],Λ j − qz j Λ j = L (cid:88) k (cid:54) = j Λ j − Λ k z j − z k + 1 z j L (cid:88) l =1 Λ l , (28)where 2 q = G − + 2 M − L −
1. A continuum limit approximation applied to(28) leads to the following integral equation,Λ( (cid:15) ) − q(cid:15) Λ( (cid:15) ) = (cid:90) ω d δ ρ ( δ ) Λ( (cid:15) ) − Λ( δ ) (cid:15) − δ + 1 (cid:15) (cid:90) ω d δ ρ ( δ )Λ( δ ) . (29)16he objective now is to derive a solution to (29). We will establish severaluseful results before the solution is stated. Let R ( (cid:15) ) = (cid:112) ( (cid:15) − a )( (cid:15) − b ) , (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) = 2 | q |√ ab , (cid:90) ω d (cid:15) (cid:15)ρ ( (cid:15) ) R ( (cid:15) ) = 1 G . (30)Note that a, b are determined by G and q . We require them to be a complex-conjugate pair, or both negative real numbers, i.e. a = ¯ b ∈ C \ R or a ≤ b ≤ . (31)As a result a + b ∈ R and ab ≥
0. The function R ( (cid:15) ) is the elementaryreal-valued square-root with (cid:15) ∈ (0 , ω ). In the limit where q = 0, we require b to vanish also. In this case (30) becomes R ( (cid:15) ) = (cid:112) (cid:15) ( (cid:15) − a ) , (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) = C, (cid:90) ω d (cid:15) (cid:15)ρ ( (cid:15) ) R ( (cid:15) ) = 1 G , (32)where a is real and determined by G . All the following results and proofsleading to the expression for the ground-state energy for the closed modelassume q (cid:54) = 0 (cid:54) = 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 similarcalculation and is omitted. Lemma 1.
For all γ, δ, (cid:15) ∈ (0 , ω ) , the following identity holds R ( γ ) − R ( (cid:15) )( γ − (cid:15) )( δ − (cid:15) ) − R ( γ ) − R ( δ )( γ − δ )( δ − (cid:15) ) = R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) . (33) Proof: Let θ L ( γ, δ, (cid:15) ) = R ( γ ) − R ( (cid:15) )( γ − (cid:15) )( δ − (cid:15) ) − R ( γ ) − R ( δ )( γ − δ )( δ − (cid:15) ) ,θ R ( γ, δ, (cid:15) ) = R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) . For δ (cid:54) = γ (cid:54) = (cid:15) (cid:54) = δ , θ L ( γ, δ, (cid:15) ) − θ R ( γ, δ, (cid:15) ) = 0 , the proof of which is straightforward and omitted here. It can be shown that θ L ( γ, δ, (cid:15) ) = 1 (cid:0) R ( γ ) + R ( (cid:15) ) (cid:1)(cid:0) R ( γ ) + R ( δ ) (cid:1) × (cid:18) − R ( γ ) + (cid:0) γ − ( a + b ) (cid:1) δ + (cid:15) − ( a + b ) R ( δ ) + R ( (cid:15) ) + ( a + b ) (cid:15)δ − ab ( δ + (cid:15) ) (cid:15)R ( δ ) + δR ( (cid:15) ) (cid:19) . his representation of θ L ( γ, δ, (cid:15) ) has no singularities for all γ, δ, (cid:15) ∈ (0 , ω ) .The same expression can be found for θ R ( γ, δ, (cid:15) ) . Therefore θ L ( γ, δ, (cid:15) ) − θ R ( γ, δ, (cid:15) ) continuously extends to vanish for all γ, δ, (cid:15) ∈ (0 , ω ) . (cid:3) Corollary 2.
Let I = 2 (cid:90) ω (cid:90) ω d δ d γ ρ ( γ ) ρ ( δ ) R ( δ ) (cid:18) R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) (cid:19) , J = (cid:18)(cid:90) ω d γ ρ ( γ ) R ( γ ) R ( γ ) − R ( (cid:15) ) γ − (cid:15) (cid:19) , then I = J − q ab . Proof: Exchanging the variables δ and γ in the integrand, and then perform-ing a change of order of integration yields I = 2 (cid:90) ω (cid:90) ω d δ d γ ρ ( γ ) ρ ( δ ) R ( γ ) (cid:18) R ( γ ) − R ( (cid:15) )( γ − (cid:15) )( δ − (cid:15) ) − R ( γ ) − R ( δ )( γ − δ )( δ − (cid:15) ) (cid:19) = 2 (cid:90) ω (cid:90) ω d δ d γ ρ ( γ ) ρ ( δ ) R ( γ ) (cid:18) R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) (cid:19) , where the second step is due to Lemma 1 . Now we add up two distinctexpressions for I / , I = (cid:90) ω (cid:90) ω d δ d γ ρ ( γ ) ρ ( δ ) R ( δ ) (cid:18) R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) (cid:19) + (cid:90) ω (cid:90) ω d δ d γ ρ ( γ ) ρ ( δ ) R ( γ ) (cid:18) R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) (cid:19) = (cid:90) ω (cid:90) ω d δ d γ ρ ( γ ) ρ ( δ ) R ( δ ) + R ( γ ) R ( δ ) R ( γ ) (cid:18) R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) (cid:19) . Since (cid:0) R ( δ ) + R ( γ ) (cid:1) (cid:18) R ( δ ) − R ( (cid:15) )( δ − (cid:15) )( γ − (cid:15) ) − R ( δ ) − R ( γ )( δ − γ )( γ − (cid:15) ) (cid:19) = − (cid:0) R ( γ ) − R ( (cid:15) ) (cid:1)(cid:0) R ( δ ) − R ( (cid:15) ) (cid:1) ( γ − (cid:15) )( δ − (cid:15) ) , e conclude that I = (cid:90) ω d δ ρ ( δ ) R ( δ ) R ( δ ) − R ( (cid:15) ) δ − (cid:15) (cid:90) ω d γ ρ ( γ ) R ( γ ) R ( γ ) − R ( (cid:15) ) γ − (cid:15) − (cid:90) ω (cid:90) ω d δ d γ ρ ( δ ) ρ ( γ ) R ( δ ) R ( γ )= J − q ab . (cid:3) Lemma 3.
Let K = 2 (cid:15) (cid:90) ω (cid:90) ω d γ d δ ρ ( γ ) ρ ( δ ) R ( γ ) R ( γ ) − R ( δ ) γ − δ , O = (cid:90) ω d δ ρ ( δ ) R ( δ ) − R ( (cid:15) )( δ − (cid:15) ) (cid:15) + (cid:90) ω d δ ρ ( δ ) R ( (cid:15) ) R ( δ ) R ( δ ) − R ( (cid:15) )( δ − (cid:15) ) (cid:15) , P = (cid:90) ω (cid:90) ω d (cid:15) d δ ρ ( (cid:15) ) ρ ( δ ) (cid:15)δR ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ , then K = 4 | q | G √ ab (cid:15) − q ab a + b(cid:15) , O = 1 G(cid:15) − | q |√ ab a + b(cid:15) + 2 | q |√ ab , P = a + b G − | q |√ abG + 1 G (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) . Proof: Since K = 2 (cid:15) (cid:90) ω (cid:90) ω d γ d δ ρ ( γ ) ρ ( δ ) R ( γ ) R ( γ ) − R ( δ ) γ − δ = 2 (cid:15) (cid:90) ω (cid:90) ω d γ d δ ρ ( δ ) ρ ( γ ) R ( δ ) R ( δ ) − R ( γ ) δ − γ , again by adding up two distinct expressions for K / , we have K = 1 (cid:15) (cid:90) ω (cid:90) ω d γ d δ ρ ( γ ) ρ ( δ ) R ( γ ) R ( δ ) (cid:16) R ( γ ) + R ( δ ) (cid:17) R ( γ ) − R ( δ ) γ − δ = 4 | q | G √ ab (cid:15) − q ab a + b(cid:15) . ext O = (cid:90) ω d δ ρ ( δ ) R ( δ ) 1( δ − (cid:15) ) (cid:15) (cid:16) R ( δ ) − R ( δ ) R ( (cid:15) ) + R ( (cid:15) ) R ( δ ) − R ( (cid:15) ) (cid:17) = 1 G(cid:15) + 2 | q |√ ab − | q |√ ab a + b(cid:15) . Finally, following a similar calculation as that for K , P = 12 (cid:90) ω (cid:90) ω d (cid:15) d δ ρ ( (cid:15) ) ρ ( δ ) (cid:15)δ (cid:18) R ( δ ) + 1 R ( (cid:15) ) (cid:19) R ( (cid:15) ) − R ( δ ) (cid:15) − δ = a + b G − | q |√ abG + 1 G (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) . (cid:3) Now we prove the following:
Proposition 4.
The following function Λ( (cid:15) ) = | q |√ ab R ( (cid:15) ) (cid:15) + q(cid:15) − (cid:90) ω d γ ρ ( γ ) R ( γ ) · R ( γ ) − R ( (cid:15) ) γ − (cid:15) (34) is a solution to the integral equation (29) with a, b subject to (30) and (31) .Proof: Λ( (cid:15) ) = J q ab R ( (cid:15) ) (cid:15) + q (cid:15) − | q |√ ab R ( (cid:15) ) (cid:15) (cid:90) ω d γ ρ ( γ ) R ( γ ) R ( γ ) − R ( z ) γ − z − q(cid:15) (cid:90) ω d γ ρ ( γ ) R ( γ ) R ( γ ) − R ( (cid:15) ) γ − (cid:15) + 2 q | q |√ ab R ( (cid:15) ) (cid:15) , and − q(cid:15) Λ( (cid:15) ) = q(cid:15) (cid:90) ω d γ ρ ( γ ) R ( γ ) R ( γ ) − R ( (cid:15) ) γ − (cid:15) − q | q |√ ab R ( (cid:15) ) (cid:15) − q (cid:15) , hence Λ( (cid:15) ) − q(cid:15) Λ( (cid:15) ) = J − | q |√ ab (cid:90) ω d γ ρ ( γ ) R ( (cid:15) ) R ( γ ) R ( γ ) − R ( (cid:15) ) (cid:15) ( γ − (cid:15) )+ q ab − q ab ( a + b ) (cid:15) . n the RHS of (29) , we first use Lemma 1 and simplify the following term (cid:90) ω d δ ρ ( δ ) Λ( δ ) − Λ( (cid:15) ) δ − (cid:15) = I | q |√ ab (cid:90) ω d δ ρ ( δ ) R ( δ ) − R ( (cid:15) )( δ − (cid:15) ) (cid:15) − | q |√ ab (cid:90) ω d δ ρ ( δ ) R ( δ ) 1 δ(cid:15) − q (cid:90) ω d δ ρ ( δ ) 1 δ(cid:15) . Since (cid:15) (cid:90) ω d δ ρ ( δ )Λ( δ ) = − K | q |√ ab (cid:90) ω d δ ρ ( δ ) R ( δ ) δ(cid:15) + q (cid:90) ω d δ ρ ( δ ) δ(cid:15) , we have (cid:90) ω d δ ρ ( δ ) Λ( (cid:15) ) − Λ( δ ) (cid:15) − δ + 1 (cid:15) (cid:90) ω d δ ρ ( δ )Λ( δ )= I − K | q |√ ab (cid:90) ω d δ ρ ( δ ) R ( δ ) − R ( (cid:15) )( δ − (cid:15) ) (cid:15) . Finally, Λ( (cid:15) ) − q(cid:15) Λ( (cid:15) ) − (cid:90) ω d δ ρ ( δ ) Λ( (cid:15) ) − Λ( δ ) (cid:15) − δ − (cid:15) (cid:90) ω d δ ρ ( δ )Λ( δ )= J − I K − | q |√ ab O + q ab − q ab ( a + b ) (cid:15) = 0 . (cid:3) The integral approximation for (27) corresponding to the solution Λ( (cid:15) ) givenin (34) is t ( (cid:15) ) = − G − − M + 2 (cid:15) Λ( (cid:15) )= − (cid:90) ω d δ ρ ( δ ) δR ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ . (35)As a direct consequence of Proposition 4 , we have the following result,
Corollary 5.
The function t ( (cid:15) ) defined in (35) is a solution to the integralapproximation of (11) , t ( (cid:15) ) = 1 G + 2 (cid:90) ω d δ ρ ( δ ) (cid:15) ( t ( (cid:15) ) − t ( δ )) (cid:15) − δ . (36)21his particular solution t ( (cid:15) ) corresponds to the ground state since from (6)in the continuum limit, E = G (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15)t ( (cid:15) ) + 12 (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) = − G (cid:90) ω (cid:90) ω d (cid:15) d δ ρ ( (cid:15) ) ρ ( δ ) (cid:15)δR ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ + 12 (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) = − a + b G + | q |√ ab − (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) + 12 (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15), (37)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] andintegral approximation of both forms of Bethe Ansatz solutions discussed inSect. 3.4 and 3.5 despite its limitations.
5. The p + ip model interacting with its environment Consider the open model (2) with extra terms governed by parameter Γ, H = H + Γ L (cid:88) k =1 ( S + k + S − k ) . The Bethe Ansatz solution for (2) was derived in [7]: for each solution { v j } of the coupled equations α + L (cid:88) k (cid:54) = j v k v k − v j + L (cid:88) l =1 z l v j − z l = − β v j (cid:81) Ll =1 (1 − v j z − l ) (cid:81) Lk (cid:54) = j (1 − v j v − k ) , (38)where α = 1 + G − , β = Γ / G and j = 1 , . . . , L , there is a correspondencebetween { t j } in (7) and { v j } via a change of variables t j = −G − − L + 2 z j L (cid:88) k =1 z j − v k . (39)This translates to the fact that (7) and (38) are equivalent. It is worthmentioning that the difference between the closed and open model is that inthe closed model, it is possible for the number of Bethe roots M to be less22han L , while in the open model we must have exactly L Bethe roots [8, 21].The energy is given by (6), E = G L (cid:88) j =1 z j t j + 12 L (cid:88) k =1 z k = G α L (cid:88) k =1 v k + G β L (cid:88) j =1 (cid:81) Ll =1 (1 − v j z − l ) (cid:81) Lk (cid:54) = j (1 − v j v − k ) . (40)The corresponding eigenstate reads as [9] L (cid:89) j =1 (cid:32) Γ v j + G L (cid:88) k =1 z k z k − v j S † k (cid:33) | (cid:105) . We perform some numerical analysis to study the behaviour of the Betheroots { v j } . We consider a small-sized case where L = 5. As we send Γ from0 to a large number, numerical results (see Tab. 1) show that one of { v j } in(38) diverges around parameters ( G = 1 , Γ = 10 . G , Γ), we have two real roots meeting to form a complex-conjugate pairand vice versa, see Fig. 5, 6 and 7.
The numerical analysis of the distribution of the Bethe roots gives no clearindications of a solution curve for a large particle number. However, since(38) is invariant under complex conjugation, again assuming the absence ofdegeneracy in the spectrum of the set of conserved operators, each solutionset { 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 considerthe integral approximation for (7), which reads as t ( (cid:15) ) = 1 G + F (cid:15) − G + 2 (cid:90) ω d δ ρ ( δ ) δ(cid:15) − δ (cid:0) t ( (cid:15) ) − t ( δ ) (cid:1) , (41)where F = − able 1: Numerical solutions for { v k } in (38) for various values of Γ while L = 5, G = 1and z l = 1 , , . . . ,
5, are fixed. The energy is calculated according to (40). The results aregrouped corresponding to Fig. 5, 6 and 7 respectively. Γ Bethe roots Energy10 − ( − . ± . × − , − × − , − . − . − . . − . ± . − . − . − . − . . − . ± . − . − . − . − . . − . ± . − . ± . − . − . . − . ± . − . ± . − . − . − . ± . − . ± . − . − . . . ± . − . ± . − . − . . . ± . − . ± . − . − . .
15 0 . ± . − . ± . − . − . .
212 0 . ± . − . ± . − − . .
213 0 . ± . − . ± . − − . . . ± . − . ± . − . . . ± . − . ± . . − . . ± . . ± . . − . . ± . . ± . . − . . ± . . ± . . − . . . ± . . ± . . − . .
455 1 . ± . . ± . . − . . . . . ± . . − . . . . . ± . . − . . . . ± . . − . .
013 1 . . . . . − . . . . . . . − . .
34 1 . . . . . − . . . . ± . . . − . .
925 1 . . ± . . . − . .
96 1 . . ± . . ± . − . . . ± . . ± . − . .
59 1 . . ± . . ± . − . .
64 1 . . ± . . . − . .
565 1 . . ± . . . − . . . . . . − . . . . . . − . . . . . . − . . . . . . − . a) (b)(c) (d) Figure 5: The Bethe roots { v k } are depicted by solid dots. The crosses mark positions z l = 1 , , , ,
5. In (a), three roots are coincident at 0. As Γ increases, the threecoincident roots depart from 0. From (b) to (c), 2 real roots meet for a value of Γ between2 . .
42. After the two roots meet, they separate and form a complex-conjugate pair.In (d), the negative real root diverges as Γ continues to increase. (a) (b) Figure 6: The inverted roots { v − k ∈ C } are depicted by stars to visualize the divergingreal root. The triangles are the { z − l } . From (a) to (b), the inverted real root traverses 0and becomes positive. a) (b)(c) (d)(e) (f) Figure 7: In (a), the real root decreases from + ∞ and approaches 5 as Γ continues toincrease. From (b) to (c), two complex-conjugate roots meet at position 2 and becomereal. From (c) to (d), two complex-conjugate roots meet at position 4 and become real.From (d) to (e), two real roots meet and separate as a complex-conjugate pair at position3, and then two other real roots behave similarly at position 5. In (f) as Γ becomes large,all roots approach { z l } from right. T j defined in (3).From (55) and (56), we derive the mean-field expression for T j as (cid:104)T j (cid:105) = − z j + F (∆ + F ) GR ( z j ) + L (cid:88) k (cid:54) = j z k − z j (cid:18) z j R ( z k ) R ( z j ) − z k (cid:19) , j = 1 , . . . , L, (42)where R ( z ) = (cid:112) z + z (∆ + F ) , ∆ ∈ R . The integral approximation for (42) is T ( (cid:15) ) = − (cid:15) + F (∆ + F ) GR ( (cid:15) ) + (cid:90) ω d δ ρ ( δ ) 1 δ − (cid:15) (cid:18) (cid:15) R ( δ ) R ( (cid:15) ) − δ (cid:19) . (43)We will establish that (43) is in fact a solution to (41), i.e. t ( (cid:15) ) = T ( (cid:15) ) , with the “gap” equation determining the value of ∆,∆ G (∆ + F ) = (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) R ( (cid:15) ) . (44)In addition, this solution t ( (cid:15) ) corresponds to the ground state. We againestablish some useful results to assist our proof. Lemma 6.
Given equation (44) , let T ( (cid:15) ) = − (cid:90) ω d δ ρ ( δ ) δ R ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ , C = (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) , then T ( (cid:15) ) = ∆ G (∆ + F ) + 2 (cid:90) ω d δ ρ ( δ ) δ ( T ( (cid:15) ) − T ( δ )) (cid:15) − δ , (45) T ( (cid:15) ) = CR ( (cid:15) ) − L − (cid:15) (cid:90) ω d δ ρ ( δ ) 1 R ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ , (46) t ( (cid:15) ) = − F G (∆ + F ) R ( (cid:15) ) (cid:15) + T ( (cid:15) ) . (47)27 roof: In (32) , set G = ∆ G (∆ + F ) , q = b = 0 , a = − (∆ + F ) , then we have R ( (cid:15) ) = R ( (cid:15) ) and C = C . Hence from (35) and by Corollary5 , T ( (cid:15) ) = t ( (cid:15) ) satisfying (36) , T ( (cid:15) ) = ∆ G (∆ + F ) + 2 (cid:90) ω d δ ρ ( δ ) δ ( T ( (cid:15) ) − T ( δ )) (cid:15) − δ . Next T ( (cid:15) ) = − (cid:90) ω d δ ρ ( δ ) δ − (cid:15) + (cid:15) R ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ = CR ( (cid:15) ) − L − (cid:15) (cid:90) ω d δ ρ ( δ ) 1 R ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ . Finally, t ( (cid:15) ) = − (cid:15) + F (∆ + F ) GR ( (cid:15) ) + (cid:90) ω d δ ρ ( δ ) 1 δ − (cid:15) (cid:18) (cid:15) R ( δ ) R ( (cid:15) ) − (cid:15) − δ + (cid:15) (cid:19) = − F G (∆ + F ) R ( (cid:15) ) (cid:15) + CR ( (cid:15) ) − L − (cid:15) (cid:90) ω d δ ρ ( δ ) 1 R ( δ ) R ( δ ) − R ( (cid:15) ) δ − (cid:15) . (cid:3) Now we show the following:
Proposition 7.
The following function t ( (cid:15) ) = − F G (∆ + F ) R ( (cid:15) ) (cid:15) − (cid:90) ω d δ ρ ( δ ) δ R ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ (48) is a solution to the integral equation (41) with ∆ subject to (44) .Proof: We adopt the notation introduced in Lemma 6 , rewriting t ( (cid:15) ) = − F G (∆ + F ) R ( (cid:15) ) (cid:15) + T ( (cid:15) ) . Since t ( (cid:15) ) = F G (∆ + F ) + F (cid:15) − G − F G (∆ + F ) R ( (cid:15) ) (cid:15) T ( (cid:15) ) + T ( (cid:15) ) , nd (cid:90) ω d δ ρ ( δ ) δ(cid:15) − δ (cid:0) t ( (cid:15) ) − t ( δ ) (cid:1) = − F G (∆ + F ) (cid:90) ω d δ ρ ( δ ) δ(cid:15) − δ (cid:18) R ( (cid:15) ) (cid:15) − R ( δ ) δ (cid:19) + 2 (cid:90) ω d δ ρ ( δ ) δ(cid:15) − δ ( T ( (cid:15) ) − T ( δ )) , then from (45) we are able to simplify equation (41) as − F + ∆ G (∆ + F ) + 1 G + 2 F G (∆ + F ) (cid:18) R ( (cid:15) ) (cid:15) T ( (cid:15) ) − (cid:90) ω d (cid:15) ρ ( δ ) δ(cid:15) − δ (cid:18) R ( (cid:15) ) (cid:15) − R ( δ ) δ (cid:19)(cid:19) = 0 . (49) Since (cid:90) ω d δ ρ ( δ ) δ(cid:15) − δ (cid:18) R ( (cid:15) ) (cid:15) − R ( δ ) δ (cid:19) = (cid:90) ω d δ ρ ( δ ) (cid:18) R ( (cid:15) ) − R ( δ ) (cid:15) − δ − R ( (cid:15) ) (cid:15) (cid:19) , R ( (cid:15) ) (cid:15) T ( (cid:15) ) = C (cid:15) + (cid:15) (∆ + F ) (cid:15) − L R ( (cid:15) ) (cid:15) − R ( (cid:15) ) (cid:90) ω d δ ρ ( δ ) R ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ , and we use the special case of q = b = 0 for expression O in Lemma 3 , where (2 | q | / √ ab ) , ( a + b ) and G − are replaced by C , − (∆+ F ) and ∆ G − (∆+ F ) − respectively, consequently we have − R ( (cid:15) ) (cid:90) ω d δ ρ ( δ ) R ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ − (cid:90) ω d δ ρ ( δ ) R ( (cid:15) ) − R ( δ ) (cid:15) − δ = − ∆ G (∆ + F ) − C (∆ + F ) − C (cid:15). Hence R ( (cid:15) ) (cid:15) T ( (cid:15) ) − (cid:90) ω d δ ρ ( δ ) δ(cid:15) − δ (cid:18) R ( (cid:15) ) (cid:15) − R ( δ ) δ (cid:19) = − ∆ G (∆ + F ) . Now the LHS of (49) is reduced to the following − F + ∆ G (∆ + F ) + 1 G − F ∆ G (∆ + F ) = 0 . (cid:3)
29e derive the energy expression using t ( (cid:15) ) from (6), E = G (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) t ( (cid:15) ) + 12 (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) = − F F ) (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) + 12 (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) − G (cid:90) ω (cid:90) ω d (cid:15) d δ ρ ( (cid:15) ) ρ ( δ ) (cid:15)δ (cid:0) R ( (cid:15) ) − R ( δ ) (cid:1) R ( δ )( (cid:15) − δ ) . Again from the special case for expression P in Lemma 3 , where b , a and G − are replaced by 0, − (∆ + F ) and ∆ G − (∆ + F ) − respectively, we have − G (cid:90) ω (cid:90) ω d (cid:15) d δ ρ ( (cid:15) ) ρ ( δ ) (cid:15)δ (cid:0) R ( (cid:15) ) − R ( δ ) (cid:1) R ( δ )( (cid:15) − δ )= ∆ G − ∆2(∆ + F ) (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) . Hence E = ∆ G + 12 (cid:90) ω d (cid:15) ρ ( (cid:15) ) (cid:15) − (cid:90) ω d (cid:15) ρ ( (cid:15) ) R ( (cid:15) ) . This energy expression is consistent with (50) (see Appendix A) for theground state, which is derived from mean-field analysis.
6. Conclusion
The continuum limit approximation for calculating the ground-state en-ergy of the p + ip model, via the Bethe Ansatz solution, was studied. Startingwith the closed model, we revisited the formulations of [3, 4, 13] which un-dertake calculations by assuming a form of density function for the Betheroot distribution. It was found that this approach does not provide a con-sistent 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 collapseat the origin [1, 3, 4, 5]. An alternative approach, which avoids the need topostulate a form for the Bethe root density, was proposed in terms of thecoupled equations satisfied by the conserved operator eigenvalues [17, 18].In this case a solution corresponding to the ground state in the continuum30imit was found. Curiously, the expression obtained for the ground-state en-ergy coincides with that obtained by the Bethe root distribution. That is tosay that although the Bethe root approach involves a flawed methodology,it nonetheless produces the correct answer! In both cases the ground-stateenergy per particle is exactly the same as the mean-field prediction in thelimit of infinite particle number.The conserved operator eigenvalue approach was then extended to accom-modate the open case based on results from [8, 9, 21]. Again, it was foundthat the result for the ground-state energy is in agreement with mean-fieldcalculations. The open case can be considered as a model allowing for the ex-change of particle between the system and the environment. It is importantto note that in this case there is no signature of any quantum phase tran-sition, no matter how weak the environment coupling is, which is in starkcontrast to the closed system. A similar scenario was considered in [22] wherethe environment was modelled by a single bosonic degree of freedom in sucha way that integrability was preserved. There too, arbitrarily small couplingto the environment was found to annihilate the existence of any quantumphase transition.For future work it would be natural to extend this analysis to calculatethe leading order finite-size correction to the ground-state energy, for bothopen and closed models. For the closed s -wave pairing systems there areresults obtained by series expansions [23], which can be extrapolated andshown to be valid more generally [24]. Obtaining analogous expressions forthe open and closed p + ip -pairing systems appears to be entirely feasible. Acknowledgements
This work was supported by the Australian Research Council through Dis-covery Project DP150101294.
Appendix A Mean-field analysis for the open model
We introduce the following notation adopted in [6]: Q † = L (cid:88) k =1 z k S + k , Q = L (cid:88) k =1 z k S − k . G(cid:104) Q (cid:105) , ˆ∆ † = 2 G(cid:104) Q † (cid:105) and ∆ = | ˆ∆ | . The extended Hamiltonian(2) can be rewritten as H ≈ L (cid:88) k =1 z k S zk − G Q (cid:104) Q † (cid:105) − G Q † (cid:104) Q (cid:105) + G(cid:104) Q † (cid:105)(cid:104) Q (cid:105) + µ (cid:32) L (cid:88) k =1 (cid:104) S zk (cid:105) − L (cid:88) k =1 S zk (cid:33) − F Q † − F Q + 12 L (cid:88) k =1 z k = ∆ G + 12 L (cid:88) k =1 z k + L (cid:88) k =1 (cid:18) z k / − ( ˆ∆ + F ) z k / − ( ˆ∆ † + F ) z k / − z k / (cid:19) . Note that in the mean-field approximation for this extended model, the La-grange multiplier is µ = 0. The matrix form is derived from the representa-tion of H acting on ( C ) ⊗ L . Consider the following eigenvalue problem, (cid:18) z k / − ( ˆ∆ + F ) z k / − ( ˆ∆ † + F ) z k / − z k / (cid:19) (cid:18) v k u k (cid:19) = λ k (cid:18) v k u k (cid:19) , k = 1 , . . . , L, by minimizing each eigenvalue, we derive the ground-state energy E = ∆ G + 12 L (cid:88) k =1 z k − L (cid:88) k =1 (cid:113) z k + χ z k , (50)where χ = | ˆ χ | , ˆ χ = ˆ∆ + F . Then we calculate the following (cid:104) S zk (cid:105) = 12 ( | v k | − | u k | ) = − z k (cid:112) z k + χ z k , (cid:104) S + k (cid:105) = u k v ∗ k = ˆ χ † z k (cid:112) z k + χ z k , (cid:104) S + k (cid:105) = u ∗ k v k = ˆ χz k (cid:112) z k + χ z k . Apply Hellmann-Feynman theorem, since (cid:28) ∂ H ∂ G (cid:29) = − ∆ G , ∂ E ∂ G = ∆ G − L (cid:88) k =1 χz k (cid:112) z k + χ z k ∂χ∂ G , χ = ∆ + F ( ˆ∆ + ˆ∆ † ) + F , (51) ∂χ∂ G = ∆ G χ + F G χ Re ˆ∆ , we have L (cid:88) k =1 z k (cid:112) z k + χ z k = 1 G ∆ ∆ + F Re ˆ∆ . (52)Furthermore from (51), ∂χ∂F = Re ˆ∆ χ + Fχ , and (cid:28) ∂ H ∂F (cid:29) = − (cid:32) ˆ∆ ∗ G + ˆ∆2 G (cid:33) = − Re ˆ∆2 G ,∂ E ∂F = − L (cid:88) k =1 z k (cid:112) z k + χ z k (Re ˆ∆ + F ) , hence L (cid:88) k =1 z k (cid:112) z k + χ z k = 1 G Re ˆ∆Re ˆ∆ +
F . (53)Comparing (52) and (53), and assuming F (cid:54) = 0, we conclude thatIm ˆ∆ = 0 , (54)i.e. ˆ∆ = ˆ∆ † = ∆. It follows that χ = ˆ χ = ˆ χ † = ∆+ F . Hence we immediatelyhave the following results (cid:104) S zk (cid:105) = − z k (cid:112) z k + (∆ + F ) z k , (55) (cid:104) S + k (cid:105) = (cid:104) S − k (cid:105) = 12 (∆ + F ) z k (cid:112) z k + (∆ + F ) z k , (56)and the “gap” equation∆ G (∆ + F ) = L (cid:88) k =1 z k (cid:112) z k + (∆ + F ) z k . (57)33 eferences [1] M. Ibanez, J. Links, G. Sierra, S.-Y. Zhao, Exactly solvable pairingmodel for superconductors with p x + ip y -wave symmetry, Phys. Rev. B (2009) 180501(R).[2] T. Skrypnyk, Non-skew-symmetric classical r -matrices and integrablecases of the reduced BCS model, J. Phys. A, Math. Theor. 42 (2009)472004.[3] C. Dunning, M. Ibanez, J. Links, G. Sierra, S.-Y. Zhao, Exact solutionof the p + ip pairing Hamiltonian and a hierarchy of integrable models,J. Stat. Mech.: Theory Exp. (2010) P08025.[4] S. Rombouts, J. Dukelsky, G. Ortiz, Quantum phase diagram of theintegrable p x + ip y fermionic superfluid, Phys. Rev. B (2010) 224510.[5] M. Van Raemdonck, S. De Baerdemacker, D. Van Neck, Exact solutionof the p x + ip y pairing Hamiltonian by deforming the pairing algebra,Phys. Rev. B (2014) 155136.[6] J. Links, I. Marquette, A. Moghaddam, Exact solution of the p + ip model revisited: duality relations in the hole-pair picture, J. Phys. A:Math. Theor. (2015) 374001 (2015).[7] I. Lukyanenko, P.S. Isaac, J. Links, An integrable case of the p + ip pairing Hamiltonian interacting with its environment, J. Phys. A: Math.Theor. (2016) 084001.[8] J. Links, Solution of the classical YangBaxter equation with an exoticsymmetry, and integrability of a multi-species boson tunnelling model,Nucl. Phys. B (2017) 117.[9] P.W. Claeys, S. De Baerdemacker, D. Van Neck, Read-Green resonancesin a topological superconductor coupled to a bath, Phys. Rev. B (2016) 220503(R).[10] R.W. Richardson, A restricted class of exact eigenstates of the pairing-force Hamiltonian, Phys. Lett. (1963) 277.3411] M. Gaudin, ´Etats propres et valeurs propres de lHamiltonien dap-pariement, unpublished Saclay preprint, 1968, included in: Travaux deMichel Gaudin, Modles Exactament R´esolus, Les ´Editions de Physique,France, 1995.[12] J.M. Rom´an, G. Sierra, J. Dukelsky, Large- N limit of the exactly solv-able BCS model: analytics versus numerics, Nucl. Phys. B (2002)483.[13] L. Amico, A. Di Lorenzo, A. Mastellone, A. Osterloh, and R. Raimondi,Ann. Phys. (2002) 228.[14] O. Babelon, D. Talalaev, On the Bethe ansatz for the Jaynes-Cummings-Gaudin model, J. Stat. Mech.: Theor. Exp. (2007) P06013.[15] A. Faribault, O. El Araby, C. Str¨ater, V. Gritsev, Gaudin models solverbased on the correspondence between Bethe ansatz and ordinary differ-ential equations, Phys. Rev. B (2011) 235124.[16] A. Faribault, D. Schuricht, On the determinant representations ofGaudin models’ scalar products and form factors, J. Phys. A: Math.Theor. (2012) 485202.[17] P.W. Claeys, S. De Baerdemacker, M. Van Raemdonck, D. Van Neck, Aneigenvalue-based method and determinant representations for generalintegrable XXZ
Richardson–Gaudin models, Phys. Rev. B (2015)155102.[18] P.W. Claeys, D. Van Neck, S. De Baerdemacker, Inner products in in-tegrable Richardson–Gaudin models, SciPost Phys. (2017) 028.[19] I. Marquette, J. Links, Generalized Heine–Stieltjes and Van Vleck poly-nomials associated with two-level, integrable BCS models, J. Stat.Mech.: Theor. Exp. (2012) P08019.[20] I. Marquette, J. Links, Integrability of an extended d + id -wave pairingHamiltonian, Nucl. Phys. B (2013) 378.[21] J. Links, On completeness of Bethe Ansatz solutions for sl (2)Richardson–Gaudin systems, in “Physical and Mathematical Aspectsof Symmetries”, Proceedings of the 31st International Colloquium in35roup Theoretical Methods in Physics. Edited by S. Duarte, J.-P.Gazeau, S. Faci, T. Micklitz, R. Scherer, F. Toppan, (Springer, 2018).[22] C. Dunning, P.S. Isaac, J. Links, S.-Y. Zhao, BEC-BCS crossover in a( p + ip )-wave pairing Hamiltonian coupled to bosonic molecular pairs,Nucl.Phys. B (2011) 372.[23] R.W. Richardson, Pairing in the limit of a large number of particles, J.Math. Phys. (1977) 1802.[24] W.V. Pogosov, M. Combescot, From one to N Cooper pairs, step bystep, Physica C471