Methods of solving the system of equations for the energy gap in the revisited BCS theory of superconductivity

The Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity has been revisited in a series of papers [1], [2], [3] and in this context the equation for the energy gap was generalized to a system of integral equations. The physical consequences of this change are major, leading not only to the change of the critical temperature and of energy gap, but even to a change of the order of the phase transition and to multiple solutions for the energy gap. Nevertheless, finding the solutions of the proposed system of equations is much more complicated than solving the typical BCS gap equation and requires a careful analysis. This analysis is done here and consists of the following steps:• writing the system of equations at finite temperature (kBT comparable with the energy gap Δ) and in the low temperature limit (kBT≪Δ);• separate analysis of the equations and of their solutions in the two temperature ranges, (first) kBT≪Δ and (second) kBT comparable with Δ;• presenting the methods to consistenlty searching the solutions.


Method details
Let us detail here the methods used in Ref. [3] . The starting point is the BCS Hamiltonian [4] written in terms of the creation c † k,s and annihilation c k,s operators on the states determined by the quantum numbers ( k, s ) (concretely, k is the wavevector and s ≡↑ , ↓ is the spin projection).
The free particle energy is ∈ (0) k and the number operator is denoted by ˆ n ks ≡ c † k,s c k,s . We assume and V kl ≡ 0 otherwise. Then, μˆ N is subtracted from the Hamiltonian ( ˆ N ≡ k,s c † ks c ks ) and H M = H − μ N was liniarized and diagonalized by the Bogoliubov-Valatin transformations [5 , 6] , where ξ k ≡∈ (0) k −μ, ∈ k ≡ ξ 2 k + 2 , is the energy gap , γ k 0 = u k c k ↑ − v k c † −k ↓ , and γ k 1 = v k c † k ↑ + u k c −k ↓ ( γ † and γ are quasiparticle creation and annihilation operators , respectively). The coefficients u k and v k satisfy In these notations, the energy gap satisfies [4 , 1] , Eq. (4) ensures the consistency of the formalism. In equilibrium, the populations n ki are functions of temperature, as we shall see below.

Standard formalism
In the standard formalism, μ is the chemical potential and the quasiparticle populations are [4] where β ≡ 1 / ( k B T ) and k B is the Boltzmann's constant. Plugging Eq. (5) into Eq. (4) , one obtains an equation in the unknown energy gap . Since the wavevector k influences the populations only through the single-particle energies ξ k and ∈ k , in the following we shall drop the subscript k and we shall replace the summations over k by (formal) summations over ξ or ∈ (preserving the number of states), whenever this does not lead to confusion. Furthermore, since we analyze large (macroscopic) systems, we introduce the density of states σ 0 along the ξ (or ∈ (0) ) axis and replace the summations by integrals. In the case of Eq. (4) , this leads to where we used the symmetries ∈ ξ = ∈ −ξ and n ξ (T ) = n −ξ (T ) . At T = 0 , Eq. (6) is particularly simple to solve, since n ξ ( T = 0 ) = 0 for any ξ , and one obtains which leads to σ 0 V 1 [4] .
At finite temperatures, the right hand side of Eq. (6) is a monotonically decreasing function of (for any given T below the critical temperature T c ), so Eq. (6) admits a single solution. At T = T c , the solution is = 0 and this leads to Euler's constant [4] .

The new set of equations
In Refs. [1][2][3] a new procedure was proposed, which eventually imposes some consistency restriction over the standard formalism. The Hamiltonian ˆ H M may be constructed using any parameter μ-not necessary the chemical potential-and then it can be diagonalized, to bring it in the form (2). Denoting the chemical potential by μ R , we can write the partition function as which, after the maximization, leads to the quasiparticle populations [1] is a correction to the quasiparticle energy (in the previous papers we denoted it by ˜ μ, omitting the subscript k ). The energy gap and the quasiparticle populations are obtained by solving the coupled Eqs. (4) and (9) [1 , 2] . In the quasi-continuous limit, we assume constant density of states σ 0 . Then, Eqs. (4) and (9) may be written as [1 , 2 , 3] Therefore, instead of Eq. (6) (with the populations given by Eq. 5 ), that we solve in the standard formalism, we have to solve the system (10), which is more difficult. In order to do that, we write Eq.
(10) in dimensionless quantities, defining x F ≡ βF , x = β ∈ , y ≡ β , and y R ≡ β( μ R − μ) , and introduced the notations χ F ( x F , y, y R ) and I ( x F , y, y R ) for the expressions in the middle of Eqs. (11a) and ( 11d ). We observe that n ξ and n −ξ (or, equivalently, n x and n −x ) may be different, so we wrote them explicitly.

Low temperature limit
The system of Eq. (11) has multiple solutions and the integrals that appear in Eqs. (11a) and ( 11d ) make it difficult to consistently search for them. For this reason, we start with the study of the equations in the low temperature limit, where we can analytically calculate the integrals and find the solutions easier. As in Ref. [2] , we write The advantage of using the forms (12) is that when β → ∞ (i.e., when T → 0 ), n ξx and n −ξx converge to either 0 or 1, depending on the signs of m ξx and m −ξx .
The method for finding the quasiparticle population at zero temperature. By defining t ≡ ξ / = ± r 2 − 1 , we may write m ξ and m −ξ in a single, more convenient form, where E( t; a, b ) is a function of t, which depends on the parameters a and b. The sign of E determines the sign of m t and the discriminant of this second order polynomial is , then m (t) > 0 and n t = 0 , plus n t 1 = In the following we shall discuss only the case a > 0 and the situation a < 0 may be obtained from this one, by changing b → −b and t → −t. Then, if D > 0 , the solutions (20) exist and t 2 > 0 , whereas t 1 > 0 if and only if ab + 1 > 0 ( t 1 < 0 if ab + 1 < 0 ). Using these results, we can calculate the quantities in Eq. (11) : where we introduced the notation A ≡h ω c / , which, in general, is much bigger than 1, and B ( A, a, b ) , which is the expression on the right hand side (r.h.s.) From Eq. (15a) we find the minimum value A 0 of A , given by the standard zero temperature BCS equation Furthermore, noticing that Eq. (15a) may be written We introduce the notation M ≡ A/ A 2 + 1 and M 0 ≡ A 0 / A 2 0 + 1 . Since A 1 in general, then M 1 .
If we chose σ 0 V = 0 . 2 and = 2 × 10 −4 eV (the value of energy gap in an Al superconductor at zero temperature), we obtain A 0 ≈ 74 . 2 and 1 − M ≤ 1 − M 0 10 −4 1 . Nevertheless, in the following we shall work with finite (but large) values of A and we will show that in general the relative differences in the physical quantities calculated with M 1 and M = 1 are of the order of 1 − M. Therefore, the approximation M = 1 is good enough for all the physically relevant cases. This will make (15b) an equation for b(a ) , which, if plugged into (15a) or (18), determines (a ) [2] .
To find the solutions of the system (15) we proceed in two steps: in Step 1 ( S1 ) we find the solution of Eq. (15b) and, using this, in Step 2 ( S2 ) we solve Eq. (15a) . S1 is further divided into smaller sub-steps. In Step 1.1 ( S1.1 ) we take a ≤ 2 . In this case, we notice that D ≥ 0 if and only if b ≤ 0 , so we loo k for negative solutions of Eq. (15b) . First, we notice that there is always a solution with b = 0 , which corresponds to m (t) > 0 and n t = 0 for any t. this solution leads to the typical BCS solutions for at zero temperature. Therefore, we have only to look for the second solution. For this, we notice that the numerator B num ≡ 1 Step 1.2 ( S1.2 ) corresponds to a > 2 . In this case, b 0 ≡ ( a 2 − 4 ) / 4 a > 0 , so, in principle, a solution b > 0 of Eq. (15b) may exist. But, in such a case, the solution b sol should be smaller than b in v [so that First, we notice that for any a > 2 and Then, we calculate the series expansion which implies 1 /a inv 1 − M.
In the general case, when M 1 , if a ≥ a in v , then b in v ≥ 0 and, at b = 0 , we obtain (since a and A 1 ). To calculate b in v , we use the fact that a > a in v 1 and b in v ≥ 0 Using also the assumption that b in v /a 1 , we obtain the expansion Neglecting the term proportional to 1 /a , we obtain which implies b inv 1 − M and, therefore, 0 < b in v 1 . Calculating the derivative, we obtain Fig. 2. ∂ B ( A , a , b ) /∂ a vs. a (for a > 2 ) and  1 , then b → b in v , and both, B ( A , a , b → b in v ) and ∂ B ( A , a , b →  To prove that there are also no negative solutions of Eq. (15b) for a > 2 , we calculate the derivative ∂ B ( A, a, b ) /∂ a and numerically show that it is negative (see Fig. 2 ). Since we can also notice that A, a, b ) for any b < 0 and a > 2 , so Eq. (15b) has no solution for a > 2 . Physical quantities at zero temperature and special cases Once we determine the solution b of Eq. (15b) , we can proceed to Step 2 and calculate directly from Eq. (18) . Noticing that the relative error in the determination of / 0 is of the order 10 −4 or smaller, we can ignore the effect of finite A and A 0 and calculate for M = 1 (i.e., A, A 0 → ∞ ). We shall analyze the situation when a ∈ [ 0 , 2 ] , since for a > 2 there are no solutions, whereas for a < 0 the solutions are similar to the ones for a > 0 . As mentioned above, for any a ∈ [ 0 , 2 ] there are two solutions, the first of which corresponds to b = 0 and = 0 (the typical BCS solution at zero temperature), whereas the second solution corresponds to b < 0 and < 0 , as presented in Fig. 3 . We observe that li m a → 0 b(a ) = −∞ , such that li m a → 0 ab(a ) ≡ l ab is finite (and negative). To calculate the limit l ab , we first observe that which, for M = 1 , gives [2] where 2 is the second solution of the system (15) In Fig. 4 we show the variation of t 1 , 2 and r 1 , 2 with the scaled asymmetry ( μ R − μ) / 0 . The quantities r 1 , 2 represent the scaled BCS quasiparticle energies corresponding to t 1 , 2 , namely r 2 = t 2 2 + 1 and r 1 = s gn ( t 1 ) t 2 In (b) we also plot the energy gap ( ± ) and the relative chemical potential ( μ R − μ) / , for exemplification.

Conservation of the number of particles.
For the first solution of Eq. (15), the average total number of particle N is equal to N μ , which is the number of states up to the free single-particle energy μ. For the second solution, the total number of particles in the low temperature limit is [1 , 2] where we used the notation N to show explicitly that the total number of particles N represents an average value. In Fig. 5 we plot the function ( N − N μ ) / ( 2 σ 0 ) . Therefore, in the case when N is fixed and is different from N μ , only the second solution may be realized and μ R becomes a function of μ.

Finite temperatures
To find the solutions of Eq. (15) , we start from the solutions at T = 0 , obtained in the Section 1.2.1. In this case, analytical calculations are more difficult, so our investigations are mostly numerical.
The function χ F ( x F , y, Y R ) (11a) is the finite temperature correspondent of the function B ( A, a, b ) Fig. 6 . We use the logarithmic scale on the vertical axis to observe the variation of χ F,den with x F for any value of y . For x F ≥ x F, 0 , χ F,den is negative, so we do not plot the function in this region.  We notice that for positive y R , χ F,num ( x F , y, y R ) is negative, whereas χ F,den ( x F , y, y R ) changes sign as a function of x F , going through zero at some point x F, 0 ( y, y R ) . Being a monotonically increasing function of x F (see Fig. 7 , for the positive part of χ F,den ( x F , y, y R ) ), the zero of χ F,den ( x F , y, y R ) is easily found numerically for any fixed y and y R . In Fig. 6 we plot x F, 0 ( y, y R ) for some relevant values of the parameters μ R − μ and T ( y R = μ R −μ k B T ≡ 0 k B T c μ R −μ 0 T c T , where 0 / ( k B T c ) = 2 /A and A was defined aboves [4] ). In Fig. 7 we exemplify the dependence of χ F,den ( x F , y, y R ) on x F / x F, 0 and y , for the values of μ R − μ and T used in Fig. 6 -we use the variable x F / x F, 0 in order to re-scale the intervals ( x F, 0 ( y, y R ) , 0 ) into the interval ( 0 , 1 ) for all the values of y and y R . For x F < x F, 0 the function is negative.
At x F = x F, 0 ( y, y R ) , χ F ( x F , y, y R ) is divergent. Since χ F ( x F , y, y R ) is also monotonically decreasing with | x F | , Eq. (11a) has a single solution in the interval ( x F, 0 ( y, y R ) , 0 ) , which can be found using a simple numerical algorithm. In Fig. 8 we plot χ F as a function of x F / | x F, 0 | and y (as we did in Fig. 7 ), for the values of μ R − μ and T used in Figs. 6 and 7 , to emphasize the general behavior of the function, especially its divergence at x F = x F, 0 ( y, y R ) .
Once the solution of Eq. (11a) is obtained, which is a function x F ( y, y R ) , we can plug it into Eq. (11d) and solve for y , for given y R . In Fig. 9 we plot 1 / ( σ 0 V ) − I ( x F , y, y R ) and we see that the function may have two solutions in y for some values of x F and y R . Because of this property, Eq. (11d) , with x F replaced by the function x F ( y, y R ) , obtained from Eq. (11a) , has two solutions, as shown in the accompanying paper [3] .