New phenomenology from an old theory—The BCS theory of superconductivity revisited

Physica A 531 (2019) 121804

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

New phenomenology from an old theory—The BCS theory of superconductivity revisited Dragoş-Victor Anghel Institutul National de C&D pentru Fizica si Inginerie Nucleara – Horia Hulubei, Romania

highlights • • • • •

The BCS theory with attraction band (AB) asymmetric w.r.t. the chemical potential. We obtain additional solutions for the energy gap equation. Quasiparticle imbalance may appear even at zero temperature. When the AB is symmetric w.r.t. the chemical potential, we find two solutions. At zero temperature, the energy gaps for the two solutions are ∆0 and ∆0 /3.

article

info

Article history: Received 20 July 2017 Received in revised form 20 May 2019 Available online 11 June 2019

a b s t r a c t I analyze the low temperature limit of the BCS theory of s-wave single-band superconductors, when the attraction band (defined as the single-particle energy interval in which the pairing interaction is manifested) may be asymmetric with respect to the chemical potential. I discuss equilibrium systems, taking consistently into account the variation of the energy and of the total number of particles with the populations of the quasiparticle energy levels and I show that the equation for the energy gap may have several solutions. When the chemical potential is the center of the attraction band (the standard BCS assumption), one solution for the energy gap is ∆0 (the standard BCS gap, at zero temperature) and another one is ∆0 /3. If the absolute value of the difference between the center of the attraction band and the chemical potential is less than 2∆0 , then there are still two solutions, with energy gaps in the center of the attraction band and quasiparticle imbalance, which is present even at zero temperature. If the total number of electrons is conserved (not exactly, as in Richardson’s method, but on average, in the thermodynamic sense) and the attraction band is asymmetric, then one of the two solutions mentioned above cannot be realized. © 2019 Elsevier B.V. All rights reserved.

1. Introduction In Ref. [1], the Bardeen–Cooper–Schrieffer (BCS) theory of superconductivity [2,3] has been revisited and it has been shown that other solutions for the gap equation exist, beside the standard one. In the method employed in Ref. [1], the total number of particles, the total energy of the system, and the variation of the energy gap with the quasiparticle populations were explicitly taken into account in the calculation of the extremum of the thermodynamic potential and this led to a generalizations of the definition of the quasiparticle energies. It was observed that in the most simple case, i.e. when the density of states (DOS) is constant and the pairing interaction is acting in a single-particle energy interval which E-mail address: [email protected]. https://doi.org/10.1016/j.physa.2019.121804 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

D.-V. Anghel / Physica A 531 (2019) 121804

Fig. 1. The single-particle energy axis in the conduction band, with the attraction band marked with red. The chemical potential is µR , whereas µ is the center of the attraction band.

is symmetric with respect to the chemical potential, the standard BCS results are re-obtained. But if a small deviation from this simple case appears, the thermodynamics of the system may change dramatically: the energy gap and the superconductor-normal metal phase transition temperature change, the quasiparticle energies are redefined, and, most importantly, the normal metal-superconductor phase transition may become of the first order (see for example [4] for a discussion related to the order of phase transitions in different ensembles). In conclusion, the phenomenology predicted by the BCS model is much richer than it was thought before. The richness of the phenomenology observed in Ref. [1] entitles one to continue the investigation of the consequences of the method proposed there, since it sheds light on the consistency of other (well established) approaches of the BCS model. This is especially interesting, since most (if not all) high temperature superconducting compounds do not fall into the category described by the simplest BCS model. For example, if the density of particles is low, then the Fermi level may be close to the bottom of the conduction band at a distance which is comparable or smaller than the relevant phonons’ energies (see for example [5–7]). In such a case, the energy interval in which the pairing interaction is manifested (called the attraction band AB–see Fig. 1) is asymmetric with respect to the chemical potential (see also Ref. [8] and citations therein). In Ref. [1], the BCS model has been applied to systems in which the AB is asymmetric with respect to the chemical potential and it was shown that such systems may also be described by an energy gap which is centered in the AB, but has an imbalance in the quasiparticle population. This observation allowed us to find solutions which lead to the changes in the phenomenology of the BCS model specified above, which are valid solutions of the model. In this paper we present a consistent method to find the low temperature limit of the energy gap and the populations. While in Ref. [1] we found only one of the solutions, at arbitrary temperature (below the condensation temperature), here we show that there are two solutions, which we find for all the relevant values of the asymmetry of the attraction band, but only in the limit T → 0 K. In addition, in this manuscript we find the solutions for the case when the number of particles is conserved. 1.1. Basic formalism and notations Let us now briefly specify the basics of the formalism. The BCS effective Hamiltonian [2,3] is

ˆ BCS = H

∑

ϵk(0) cˆk†,s cˆk,s +

ks

∑

†

†

Vkl cˆk↑ cˆ−k↓ cˆ−l↓ cˆl↑ ,

(1)

kl

†

(0)

where cˆk,s and cˆk,s are the electron creation and annihilation operators on the single-particle state |k, s⟩, whereas ϵk is the corresponding single-particle energy. Introducing the notation bk ≡ ⟨ˆc−k↓ cˆk↑ ⟩, writing cˆ−l↓ cˆl↑ ≡ bk + (cˆ−k↓ cˆk↑ − bk ), ˆ BCS to the first order in the quantities cˆ−k↓ cˆk↑ − bk , one obtains a quadratic Hamiltonian, which can be and expanding H ˆ where µ diagonalized [9,10]. Nevertheless, the diagonalization is performed after subtracting a convenient quantity µN, is not necessary the chemical potential of the system. After the diagonalization, one obtains the simpler Hamiltonian

ˆ BCS ≈ Hˆ M ≡ µNˆ + H

∑ k

(ξk − ϵk + ∆k b∗k ) +

∑

† † ϵk (γˆk0 γˆk0 + γˆk1 γˆk1 ),

(2)

k

√ ξ 2 + ∆2k ∑k (0) are called the quasiparticle energies. In Eq. (2) I also introduced the standard notations ξk ≡ ϵk −µ and ∆k ≡ − k Vkl bl . † † The relations between the γˆ and cˆ operators are cˆk↑ = u∗k γˆk0 + vk γˆk1 and cˆ−k↓ = −vk γˆk0 + u∗k γˆk1 , where the coefficients 2 2 uk and vk satisfy the relation vk = 1 − uk = (1 − ξk /ϵk ) /2 (we consider vk and uk real). In these notations the expression †

where γˆki and γˆki (i = 0, 1) are the quasiparticle creation and annihilation operators, respectively, whereas ϵk ≡

D.-V. Anghel / Physica A 531 (2019) 121804

3

for the energy gap is [3]

∆k ≡ −

∑

†

†

Vkl u∗l vl ⟨1 − γˆl0 γˆl0 − γˆl1 γˆl1 ⟩.

(3)

l

†

Since nki ≡ γˆki γˆki are the number operators on the states of quasiparticle energy ϵk , if µ is chosen so that it is equal to the chemical potential of the system, then one obtains the standard result for the populations [1,3], namely

⟨nki ⟩ = f (βϵk ) ≡

1 eβϵk + 1

,

(4)

where β = 1/(kB T ) is the inverse temperature. (0) In the derivation of Eq. (3) there are no assumptions regarding the free-particle spectrum ϵk or the extent of the AB. Therefore, combining Eqs. (3) and (4) one obtains a general set of equations for the energy gap and the quasiparticle populations (see for example Ref. [8,11–15] and citations therein, for applications in superconductors with different types of energy bands), but only if µ is chosen equal to the chemical potential of the system. If µ is not equal to the chemical potential, then a different set of equations should be solved to find the equilibrium distribution [1]. If the formalism is consistent, the results of the new set of equations should be valid. The reason for which Eq. (4) is valid only in some special cases is that when the number of quasiparticles changes by 1 † – say, γˆki is applied on a state with nki = 0 – the total energy of the system does not change only by ϵk , but changes by a more complicated expression due to the simultaneous variation of ∆ and ϵk with nki [1]. This fact leads to a redefinition ˆ M ⟩/∂ nk , which should be used in the maximization of the of the quasiparticle energy – in the Landau sense – as ϵ˜k = ∂⟨H partition function and which is different from ϵk . Furthermore, in the maximization of the partition function one has to take into account explicitly also the variation of the total number of particles with the number of quasiparticles. Applying this recipe, in [1] it was obtained a self-consistent set of equations, different from the set (3)–(4), which should be solved in order to obtain the energy gap and the quasiparticle populations. 1.2. Consistency of the BCS model The simplifying assumption that leads to the basic BCS theory of superconductivity is that the DOS is constant and (0) (0) that Vkl ≡ −V is constant and different from zero only if both, ϵk and ϵl , belong to the AB, which is defined here as the interval IV ≡ [µ − h ω , µ + h ω ] , centered at the chemical potential [2]. These assumptions are justified only if small ¯ c ¯ c deviations from them do not lead to new phenomenology. But, as it was shown in Ref. [1], this is not the case. A slightly asymmetric AB leads to significant qualitative and quantitative changes in the phenomenology predicted by the model. If one denotes the chemical potential of the system by µR and keeps µ as the notation for the center of the AB, then, if µ = µR , one recovers the basic BCS model. But if µ ̸= µR (the chemical potential is not the center of the AB), then the energy gap and the populations of the quasiparticle energy levels are calculated by solving a system of integral equations. This system of equations was analyzed in [1] in the grandcanonical ensemble (i.e. keeping µR constant) and it was shown that the phase transition temperature decreases if |µR − µ| increases, a quasiparticle imbalance appears in equilibrium, and the phase transition is discontinuous. Furthermore, the system of equations for the energy gap may have more than one solution at fixed temperature and chemical potential. Quasiparticle imbalance in the context of the BCS theory have been reported before for non-equilibrium superconductors (see for example Refs. [16–20]). Such non-equilibrium situations can be described also by our approach, but here I focus only on equilibrium superconductivity. The equilibrium quasiparticle imbalance [21–23] appears in the model of hole superconductivity [24,25], but the concept and the predictions of this model are, in many respects, very different from the BCS theory and from our results. This study is important not only to check the consistency of the BCS model, but also, as mentioned above, because asymmetric attraction bands with respect to the chemical potential are known in the context of multi-band superconductors or whenever the Fermi energy is close to one extremum of the conduction band (see for example [5–8,11–15,26–30]), so that the asymmetry is imposed by the limits of the conduction band. In this paper I shall discuss only systems with one conduction band and constant density of states, under the assumption that the AB is asymmetric with respect to the chemical potential (µ ̸ = µR ). Using the method of calculation introduced in [1], I will analyze the energy gap and the populations in the low temperature limit in order to determine the number of solutions and their stability in the grandcanonical ensemble. Afterwards, I will impose the conservation of the number of particles, to obtain the results corresponding to the canonical ensemble. Some details of calculations are presented in the Appendix. 2. The self-consistent set of equations In Ref. [1], following the standard procedure (see for example Ref. [3]) and maximizing of the grandcanonical partition function of a superconductor, the populations of the quasiparticle states were obtained in the form †

nki ≡ ⟨γki γki ⟩ =

1 eβ (ϵk −µ˜ k ) + 1

,

i = 0, 1,

(5a)

4

D.-V. Anghel / Physica A 531 (2019) 121804

where i indicates the type of quasiparticle and

[ ] ∑ −3 µR − µ k (1 − nk0 − nk1 ) ξk ϵk µ ˜k ≡ ξk − ∑ . −3 ϵk k (1 − nk0 − nk1 ) ϵk

(5b)

is a correction to the quasiparticle energy, µ is the center of the attraction band IV and µR is the chemical potential. One can observe that both, µ ˜ k and ϵk are independent of i, whereas Eq. (5a) reduces to (4) when µR = µ. We consider that (0) (0) Vkl ≡ −V when both, ϵk and ϵl , belong to IV , so the energy gap equation (3) simplifies to 1=

V ∑ 1 − nk0 − nk1 2

ϵk

k

.

(6)

Eqs. (5) and (6) should be solved self-consistently to obtain the populations and the energy gap of the superconductor. In the quasicontinuous limit, the summation over k is transformed into an integral over ξk ≡ ξ . Since all the quantities of interest for us depend only on ξ , in the following we shall change the notations and instead of the variable or subscript k, we shall use the variable or subscript ξ . If the density of states (DOS) is constant, σ (ξ ) ≡ σ0 , then Eq. (6) becomes 2

σ0 V

h ¯ ωc

∫ =

1 − nξ 0 − nξ 1

√

−h¯ ωc

dξ ,

ξ 2 + ∆2

(7)

√

where the BCS quasiparticle energy is ϵ = ξ 2 + ∆2 . In the zero temperature limit, if we set nξ 0 = nξ 1 = 0 for any ξ (we shall see further that there is a metastable state in which this condition is not satisfied), one can obtain an analytical expression for the energy gap in the weak coupling limit (σ0 V ≪ 1): ∆0 = 2h ¯ ωc exp[−1/(σ0 V )]. Similarly, at the critical temperature Tc the energy gap should be zero, and form Eq. (7) one obtains kB Tc = Ah ¯ ωc e−1/(σ0 V ) , where A = 2eγ /π ≈ 1.13 and γ ≈ 0.577 is the Euler’s constant [3]. To be able to solve the self-consistent set of Eqs. (5) and (6), we first identify in Eq. (5b) the constant F , such that the effective chemical potential is written µ ˜ ≡ (µR − µ)(ξ − F )/ϵ . If the DOS is constant, Eqs. (5) get a simpler form and are equivalent to

∫ h¯ ωc

ξ (1 − nξ 0 − nξ 1 ) ϵ 3 dξ , ∫ h¯ ωc (1−nξ 0 −nξ 1 )dξ

−h¯ ωc

F ≡

−h¯ ωc

nξ i =

(8a)

ϵ3

1 eβ[ϵξ −(µR −µ)(ξ −F )/ϵξ ] + 1

.

(8b)

Introducing the dimensionless variables xF ≡ β F , x ≡ βϵ , y ≡ β ∆, and yR ≡ β (µR − µ), and assuming a constant density of states, the set of Eqs. (7) and (8) can be transformed into

∫ β √(h¯ ωc )2 +∆2 xF =

∫ β √(h¯ ωc )2 +∆2 y

nξ x = e n−ξx = 1

σ0 V

(n−ξx −nξx ) dx x2

y

x−yR

(√

(

y

,

)

(9b)

+1

1

,

)

x2 −y2 −xF /x

e ∫ β √(h¯ ωc )2 +∆2

(9a)

x2 −y2

1 x2 −y2 −xF /x

√

,

√

x2

x−yR −

=

(1−n−ξx −nξx ) dx

+1

1 − n−ξx − nξx

√

(9c)

x2 − y 2

dx,

(9d)

where we wrote explicitly the populations for the positive and negative branches in Eqs. (9b) and (9c). We observe that Eqs. (9) are symmetric under the exchange yR → −yR , xF → −xF , and ξ → −ξ . To underline the dependence of ξ on the parameter x = βϵ , I introduced the notation ξx . Solving self-consistently the set of Eqs. (9) we obtain the equilibrium populations and ∆. The system (9) depends on two parameters: yR and β = 1/(kB T ). In the next section we study the solutions of the system in the limit β → ∞ (or T → 0) for different values of the parameter µR − µ = yR /β . The details of calculations are presented in the Appendix. 3. Low temperature limit and constant DOS Since the solutions for yR < 0 can be obtained from the solutions with yR > 0, by the replacement xF → −xF and exchanging nξ with n−ξ , in the following we shall study the system (9) only in the case yR ≥ 0. For this, we analyze

D.-V. Anghel / Physica A 531 (2019) 121804

5

Fig. 2. The ratio ∆(T = 0)/∆0 vs a, for the solutions of b plotted in Fig. 8.

}−1

{

the argument of the exponential function in the denominator of nξx and n−ξx . If we write nξx ≡ exp[β mξx ] + 1

}−1

{

n−ξx ≡ exp[β m−ξx ] + 1 mξx ≡ m−ξx ≡

∆( r r

√

)

√

)

r 2 − a r 2 − 1 + ab ,

∆(

and

, then (10a)

r 2 + a r 2 − 1 + ab ,

(10b)

where r = ϵ/∆ = x/y ≥ 1, a = (µR − µ)/∆ = yR /y, and b = F /∆ = xF /y. When mξx > 0, then limT →0 β mξx = ∞ and limT →0 nξx = 0, whereas if mξx < 0, then limT →0 β mξx = −∞ and limT →0 nξx = 1. Similarly, in the limit T → 0, if m−ξx > 0, then n−ξx = 0 and if m−ξx < 0, then n−ξx = 1. In Appendix we calculate the asymptotic values for b and ab, for different values of a (Fig. 8). This allows us to find the populations, the energy gap, and the total number of particles. In Fig. 2 we plot ∆(T = 0)/∆0 for the two solutions of b plotted in Fig. 8(a). For the solution b ≡ 0, ∆(T = 0) = ∆0 for any x, whereas for the solution b < 0, ∆(T = 0) ≤ ∆0 . Using Eqs. (28) and (26c) we obtain lim ∆(b<0) (T = 0) = ∆0 /3,

(11)

a→0

for the solution with b < 0. Having the solutions for b and ∆, we can calculate the quasiparticle populations and the quasiparticle imbalance. For the solutions with b = 0, the situation is trivial: n[ξx = n−ξx = 0 for any r. ] For b < 0, if 1 < a < 2, then −1 < ab < 0 and

√

√

only the branch with ξ > 0 is populated for ξ ∈ ∆ r02 − 1, ∆ r22 − 1 . If 0 < a < 1, then −4/3 < ab < −1 and both

[

√

]

branches are populated: the branch ξ < 0 is populated in the interval ξ ∈ −∆ r02 − 1, 0 , whereas the branch ξ > 0

[

]

√

is populated in the interval ξ ∈ 0, ∆ r22 − 1 . I plot the branches populations in Fig. 3 and we observe that the branch imbalance is non-zero for any a > 0, whereas lima↗2 nξx /(σ0 ∆0 ) = 1. When a ↘ 0, the population imbalance disappears, although b → −∞ √ and ∆ ↘ ∆0 /3. This situation corresponds to µR = µ, ∆ = ∆0 /3, and the population for each branch equal to σ0 ∆0 /(3 3). The total number of particles is [1] N ≡ ⟨Nˆ ⟩ = N ′ +

∑

2vk2 +

∑ k, i

k

nk

ξk , ϵk

(12)

where N ′ is the contribution coming from outside the attraction band. Introducing the notation Nµ ≡ 2 −h¯ ωc ≤ξk <0 [

N = Nµ −

∑ k

1−

] |ξk | |ξk | + (nk0 + nk1 ) + ϵk ϵk

0≤ξk ≤h ¯ ωc

∑ k

[ 1−

] ξk ξk + (nk0 + nk1 ) , ϵk ϵk

∑k≤kF k

1 we get (13)

6

D.-V. Anghel / Physica A 531 (2019) 121804

Fig. 3. The total populations (integral over the ξ or ϵ , symbolically represented as a primitive) of the branches with ξ < 0 (blue dashed line) and ξ > 0 (red solid line), for the solutions with b < 0.

Fig. 4. The difference between the total number of particles N and number of free-particle states up to µ (Nµ ) and µR (NµR ) vs a = (µR − µ)/∆.

which, in the quasicontinuous limit (assuming constant DOS) becomes N = Nµ + 2σ0

∫

h ¯ ωc

−h¯ ωc

ξ nξ dξ ϵξ

(14)

For the situation when b ≡ 0 and ∆ = ∆0 , nξx = n−ξx = 0 and N = Nµ for any a. In the case of solutions with b < 0, from Eq. (14) we obtain N − Nµ 2σ0

∫ =

√ ∆ r22 −1

√ ∆ r02 −1

If we denote by NµR ≡ 2 N − NµR 2σ0

=

∑ϵ (0) ≤µR

N − Nµ 2σ0

ξ dξ √ = ∆(r2 − r0 ). ξ 2 + ∆2 k

(15a)

the number of free-particle states up to µR , then

− (µR − µ) = ∆(r2 − r0 − a).

(15b)

For a ∈ (0, 2), N − Nµ > 0, whereas N − NµR < 0. For a = 0, N − Nµ = N − NµR = 0, whereas for a = 2, N − Nµ = 0 and N − NµR = −4σ0 ∆0 (see Fig. 4).

D.-V. Anghel / Physica A 531 (2019) 121804

7

Let us now analyze the partition function ln(Z )βµ = −

∑ [(1 − nki ) ln(1 − nki ) + nki ln nki ] − β (E − µR N),

(16)

ki

and denote limT →0 Z ≡ Z0 . We observe that

−

kB T ln Z0 2σ0

E − µR N

=

2 σ0

=

h ¯ ωc

∫

E0

2σ0

[

+ 0

] ( ) ) E0 ξ ( ϵξ nξ + n−ξ − (µR − µ) nξ − n−ξ dξ ≡ + I1 + I2 , ϵ 2σ0

(17)

where [1] E0 = −

σ0 ∆2

[

(

1 + 2 ln

2

∆0 ∆

)]

.

(18)

The first and second integrals in (17) are I1 = ∆ 2

h ¯ ωc /∆

∫

√ 1

r2

(

r2 − 1

)

nξx + n−ξx dr

(19a)

and I2 = −∆2 a

∫

h ¯ ωc /∆

nξx − n−ξx dr ,

(

)

(19b)

1

respectively, where, as before, x = βϵ = ry. For the solution with b = 0, nξ = n−ξ = 0 for any ξ , so (b=0)

−

kB T ln Z0

=

2σ0

E0

2σ0

=−

∆20 4

.

(20)

For the solution with b < 0, if a ≤ 1, then ab ≤ −1 and the integral I1 becomes I1 =

∆2 2

[ √ r2

r22

− 1 + r0

√

r02

( − 1 + ln r2 +

√

r22

)

(

√

)]

,

(21a)

( )] √ − ln r0 + r02 − 1 .

(21b)

− 1 + ln r0 +

r02

−1

whereas if 1 < a < 2, then ab > −1 and I1 =

∆2 2

[ √ r2

r22 − 1 − r0

√

(

r02 − 1 + ln r2 +

√

r22 − 1

)

The integral I2 has the same expression for any a ∈ [0, 2], namely I2 = −∆2 a (r2 − r0 ) .

(21c)

Pugging Eqs. (21) into (17) we can calculate the logarithm of the partition function, which is plotted in Fig. 5. We observe that the solutions with b = 0 are the stable ones, since the partition function takes the bigger value. The solutions with b < 0 still satisfy the extremum condition and therefore represent metastable solutions. 4. Conservation of the number of particles We see from Fig. 4 that N is in general different from both Nµ and NµR . Since the DOS is constant, in the normal metal phase the number of particles is NµR . Therefore, upon condensation, the number of particles decreases for the situations depicted in Fig. 4 (i.e. if µR remains constant). If the system formed of neutral particles, this can be the case and one may describe the system in both, canonical and grandcanonical ensembles. But, if the system is formed of charged particles, like the electrons in a superconductor, we have to take into account the charges. In such a case, a change in the number of particles would lead to charging effects and therefore would not be physically possible. To solve this problem, we must work in the canonical ensemble, i.e. we must impose that – on average, not exactly, as in Richardson’s method [31–33] – the number of particles is the same before and after the condensation. If µ = µR , the problem is very simple, since Nµ = NµR ≡ N for both solutions of the system (9). On the other hand, if µR > µ, then for the solution with b ≡ 0 and ∆ = ∆0 , the number of particles in the normal phase is NµR , whereas in the superconducting phase is Nµ (̸ = NµR ). Therefore the solutions with b = 0 cannot conserve the number of particles when µR ̸ = µ, so they are not the physical solutions. The only physical solutions left for N ̸ = Nµ are the ones with b < 0, which I analyze below. We can see from Fig. 4 that for µ (and Nµ ) fixed, N ≤ Nmax ≈ 0.227(2σ0 ∆0 ) + Nµ ≡ N(amax , µ), where Nmax is the maximum value reached by N, when a = amax ≈ 1.468. Therefore, if the number of particles satisfies the relations Nµ < N < Nmax , then the superconducting phase may be reached for two values of a, such that N = NµR . These values are denoted by a1,2 , where 0 < a1 < amax and amax < a2 < 2. If N = Nmax , then a1 = a2 = amax , whereas if N = Nµ , then a = 0, but we have solutions for both, b = 0 and b → −∞ (the second solution corresponds to ∆ = ∆0 /3). These solutions are plotted in Fig. 6(a), whereas the values of the energy gap corresponding to them are plotted in Fig. 6(b).

8

D.-V. Anghel / Physica A 531 (2019) 121804

Fig. 5. The logarithm of the partition function in the limit T → 0 for the solutions with b = 0 (solid line) and b < 0 (dashed line).

Fig. 6. The conservation of the number of particles. For each N between Nµ and Nmax we can find two values of µR (a), such that N = NµR in the superconducting phase. In (b) we plot the values of the energy gap corresponding to the a values from (a).

5. Conclusions I analyzed the BCS formalism in the low temperature limit, under the assumption that the center of the AB (attraction ˆ BCS (Eq. (1)) and Nˆ are the Hamiltonian and band) µ does not coincide with the chemical potential of the system µR . If H the total number of particles in the system, respectively, then one can apply the Bogoliubov–Valatin transformation [9,10] ˆ which is then used in the grandcanonical formalism to calculate the equilibrium properties ˆ BCS − µR N, to diagonalize H ˆ BCS − µNˆ of the system. Nevertheless, the same statistical procedure may be applied if one diagonalizes the operator H ˆ before writing the grandcanonical partition function (with the chemical potential µR ). By doing ˆ BCS − µR N, instead of H so, we obtained some surprising results, which we briefly summarize below. First, if we denote by ∆0 the energy gap in the standard BCS theory (when µ = µR ) at zero temperature and if |µR − µ| > 2∆0 , the system of Eqs. (9), which gives the energy gap, has no solutions, so the superconducting state cannot exist even at zero temperature. If 0 < |µR − µ| < 2∆0 , then the system (9) has two solutions: one with ∆(T = 0) = ∆0 and another one, with ∆(T = 0) < ∆0 (see Fig. 2) and a quasiparticle imbalance in equilibrium (see Fig. 3). The solutions with ∆(T = 0) = ∆0 are stable, whereas the other ones are metastable, since in the first case the partition function takes bigger values than in the second case. However, it is interesting to note that in the limit µR → µ, when the standard BCS theory should be obtained, the system (9) still has two solutions, one with ∆(T = 0) = ∆0 (as expected) and another one, with ∆(T = 0) = ∆0 /3. When µR → µ, the quasiparticle imbalance converges to zero for both solutions, but in the case ∆(T = 0) = ∆0 /3 quasiparticles

D.-V. Anghel / Physica A 531 (2019) 121804

9

are present, symmetrically around µ (i.e. nξ = n−ξ , as one can see in Fig. 3 at a = (µR − µ)/∆ = 0). This means that even in the standard case there is a metastable solution with smaller energy gap and non-zero quasiparticle populations at zero temperature. In Section 4 it was shown that the conservation of the number of particles (the number of particles in the normal state is equal to the number of particles in the superconducting state) does not allow the formation of the stable states (that is, with ∆(T = 0) = ∆0 ), leaving only the metastable states (∆(T = 0) = ∆0 /3) as the valid solution—and therefore they become stable. This can be explained as follows. The chemical potential of the reservoir, µR , is determined by the total number of particles—in the normal metal state (above the phase transition temperature) we have N ≡ NµR , where NµR is the number of states below µR . If Nµ is the number of states below µ and µ = µR , then N = Nµ = NµR and we recover the BCS solution, as stated above. If N = NµR > Nµ , there is a maximum value Nmax ≈ 0.227(2σ0 ∆0 ) + Nµ (see Figs. 4 and 6), up to which solutions to the problem exist. For each N in the interval (Nµ , Nmax ), we can find two values of µR which satisfy the conservation of N. The formalism is symmetric under the simultaneous exchanges µR − µ → −(µR − µ), xF → −xF , and ξ → −ξ . Therefore the results presented can be easily extended to µR < µ (i.e. a < 0) and N < Nµ . One may, of course, choose another value µ′ (which may be different from both, µ and µR ), diagonalize the operator ˆ and apply the grandcanonical formalism to calculate the physical properties of the system, as we did above ˆ BCS − µ′ N, H ˆ Comparison between the results obtained for different choices of µ′ will be done elsewhere. ˆ BCS − µN. and in Ref. [1] for H Acknowledgments Discussions with Prof. N. Plakida, Dr. G. A. Nemnes, and Dr. S. Cojocaru are gratefully acknowledged. This work has been financially supported by the UEFISCDI (project PN 19060101/2019). Travel support from Romania-JINR Titeica–Markov program is gratefully acknowledged. Appendix. Derivation of the asymptotic equations and the low temperature limits of the energy gap and populations Let us find the values of r for which mξx < 0 or m−ξx < 0. In Eqs. (10) we denote t ≡ them as

√

r 2 − 1 ≥ 0 and we rewrite

∆

(t 2 − at + ab + 1), mξx ≡ √ t2 + 1

(22a)

m−ξx ≡ √ (t 2 + at + ab + 1), t2 + 1

(22b)

∆

We can now see that mξx and m−ξx may take negative values only if the discriminant of Eqs. (22), D ≡ a2 − 4ab − 4, is positive. Then, the solutions for Eq. (22a) are

√ t1 =

a−

√

a2 − 4ab − 4

and t2 =

2

a+

a2 − 4ab − 4 2

,

(23)

whereas the solutions for Eq. (22b) are t1′ = −t2 and t2′ = −t1 . Obviously, t2 > 0 and t1′ < 0. If we denote by Ir ≡ (r1 , r2 ), the interval on which mξx (r) < 0(10a), then r2 ≡

√

t22

+1=

√ ( a 2

a − 2b +

)

√

a2 − 4ab − 4 ≥ 1.

(24a)

Since t1 ≤ 0 if and only if ab ≤ −1, then r1 =

⎧ √ ( ⎨ a 2

⎩

a − 2b −

√

)

a2 − 4ab − 4 , if ab > −1,

(24b)

1, if ab ≤ −1.

Similarly, Ir′ ≡ (r1′ , r2′ ) is the interval on which m−ξx (r) < 0(10b). Then r1′ = 1 (since t1′ = −t2 < 0) and ′

r2 =

⎧ √ ( ⎨ a 2

⎩

a − 2b −

√

)

a2 − 4ab − 4 , if ab < −1,

(24c)

1, if ab ≥ −1.

For r ∈ Ir , limT →0 nξx (T ) = 1, whereas for r ∈ [1, ∞) \ [r1 , r2 ], limT →0 nξx (T ) = 0. Similarly, for r ∈ Ir′ , limT →0 n−ξx (T ) = 1, whereas for r ∈ [1, ∞) \ [r1′ , r2′ ], limT →0 n−ξx (T ) = 0. We also observe that Ir′ ⊂ Ir , because if t1 ≥ 0, then r1 ≥ 1 and Ir′ = ∅ ⊂ Ir , whereas if t1 ≤ 0, then r1 = r1′ = 1 and r2 > r2′ (Eqs. (24a) and (24c)), which implies again Ir′ ⊂ Ir . Using this observation we see that fn (x) ≡ n−ξx − nξx , which appears in the integrand in the numerator of Eq. (9a), is different from zero only if x/y = r ∈ Int(Ir \ Ir′ ) (where Int(·) denotes the interior of an interval), whereas fd (x) ≡ 1 − n−ξx − nξx , which appears in the integrand in the denominator, is zero in the same interval. Furthermore, fd (x) = −1, if x/y = r ∈ Ir′ , and fd (x) = 1, if x/y = r ∈ (r2 , ∞).

10

D.-V. Anghel / Physica A 531 (2019) 121804

Let us now calculate ∆ and xF . I introduce the notation r0 ≡

√

(

(a/2) a − 2b −

√

)

a2 − 4ab − 4

≥ 1. If ab ≥ −1,

n−ξ = 0 and from Eqs. (9), we obtain

][ ]⎫ ⎧[ √( ) √ h h ⎪ ⎪ ¯ ωc ¯ ωc 2 2 ⎪ ⎪ + 1 r + + r − 1 0 ⎨ ∆ ⎬ 0 ∆ 1 ] [ = ln , √ ⎪ ⎪ σ0 V ⎪ ⎪ ⎩ ⎭ r2 + r22 − 1 b = √

1 r2 1

1−

(h ¯ ωc

−

/∆)2 +1

− r1 0 √

1−

1 r22

√

1−

+

1 r02

(25a)

,

(25b)

and we observe that b < 0. If we define ∆0 by the equation

⎡ 1

σ0 V

≡ ln ⎣

√ (

h ¯ ωc

+

∆0

)2

h ¯ ωc

∆0

⎤ + 1⎦ ,

(25c)

(i.e. the value of ∆ when the gap is symmetric) then, from Eq. (25a) we can write h ¯ ωc

+

h ¯ ωc

+

∆0

∆

√(

h ¯ ωc

)2

h ¯ ωc

)2

∆0

√(

∆

+1

r0 +

= r2 +

+1

√

r02 − 1

√

r22 − 1

,

(25d)

If ab < −1, then,

⎡ 1

σ0 V

= log ⎣

b = √

h ¯ ωc

∆

√

(

+

h ¯ ωc

∆ 1 r2

1−

)2

1 (h ¯ ωc /∆)2 +1

−

⎤

) ( ) ( √ √ 2 2 ⎦ + 1 − log r2 + r2 − 1 − log r0 + r0 − 1 ,

1 r0

− √

1 − r2−2 −

√

1 − r0−2

,

(26a)

(26b)

Eliminating 1/(σ0 V ) from Eq. (26a) we write h ¯ ωc

∆0

h ¯ ωc

∆

+

√(

+

h ¯ ωc

)2

∆0

√(

+1

1

= (

) h ¯ ωc 2

+1

∆

r0 +

√

)(

r02 − 1

r2 +

√

r22 − 1

).

(26c)

If ab = −1, then r0 = 1 and Eqs. (25) and (26) give the same results, implying that the functions b(a) and ∆(a) are continuous. In the following we shall make the typical assumption that ∆, ∆0 ≪ h ¯ ωc , so that we can use the approximations h ¯ ωc

∆0

√ ( +

h ¯ ωc

)2

∆0

+1≈

2h ¯ ωc

∆0

and

h ¯ ωc

∆

√ +

(

h ¯ ωc

∆

)2 +1≈

2h ¯ ωc

∆

.

(27)

The approximations (27) simplify considerably the calculations, without significantly affecting the results. Numerical and analytical analysis of Eqs. (25) and (26) (in the approximation (27)) show that for 0 < a < 2 there are two solutions (see Fig. 7): one with b = 0, nξx (T = 0) = n−ξx (T = 0) = 0, and ∆(T = 0) = ∆0 , and another one, with b < 0 and ∆(T = 0) < ∆0 ; for the second solution, nξx (T = 0) and n−ξx (T = 0) take nonzero values for some values of x. For a = 2 (green curve in Fig. 7), only the solution with b = 0 remains, whereas for a > 2, Eqs. (25) and (26) have no solutions and the superconducting phase cannot be formed at T = 0. The solutions b(a), from Eqs. (25) and (26), are plotted in Fig. 8. We see that if a < 1, then ab < −1 (Fig. 8 b). In this case, nξx (T = 0) = 1, for r = x/y ∈ [1, r2 ), and n−ξx (T = 0) = 1, for r = x/y ∈ [1, r0 ). In the limit a ↘ 0, the product ab converges to a constant, which can be readily calculated from Eq. (26b), namely lim ab = −4/3.

a→0

For a ∈ (1, 2), we have ab ∈ (−1, 0) and nξx (T = 0) = 1, if r = x/y ∈ (r0 , r2 ), whereas n−ξx (T = 0) = 0 for any r.

(28)

D.-V. Anghel / Physica A 531 (2019) 121804

11

Fig. 7. The r.h.s. of Eq. (25b), for −1 ≤ ab ≤ 0, continued by (26b), for ab ≤ −1, plotted vs. b, for different values of a. The straight black line is b vs b.

Fig. 8. In (a) are the solutions b vs a, obtained from Eqs. (25) and (26). There are two solutions, b ≡ 0 (solid blue line) and b < 0 (dashed red line). In (b) we plot the product ab vs a. For the negative function (dashed red line), lima↘0 ab = 4/3.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

D.V. Anghel, G.A. Nemnes, Physica A 464 (2016) 74. J. Bardeen, L.N. Cooper, J.R. Schrieffer, Phys. Rev. 108 (1957) 1175. Michael Tinkham, Introduction to Superconductivity, second ed., McGraw Hill, Inc., 1996. A.S. Parvan, Nuclear Phys. A 887 (2012) 1. D.M. Eagles, Phys. Rev. 164 (1967) 489. D.M. Eagles, Phys. Rev. 178 (1969) 668. D.M. Eagles, Phys. Rev. 186 (1969) 456. D. Valentinis, D. van der Marel, C. Berthod, Phys. Rev. B 94 (2016) 024511. J.G. Valatin, Comments on the theory of superconductivity, Il Nuovo Cimento (1955-1965) 7 (1958) 843. N.N. Bogoljubov, Il Nuovo Cimento (1955-1965) 7 (1958) 794. V.A. Moskalenko, Fiz. Met. Metalloved. 8 (1959) 25. V.A. Moskalenko, Fiz. Met. Metalloved. 8 (1959) 503. V.A. Moskalenko, Sov. Phys. Usp. 17 (1974) 450. V.A. Moskalenko, M.E. Palistrant, V.M. Vakalyuk, Phys.-Usp. 34 (1991) 717. M. Palistrant, A. Surdu, V. Ursu, P. Petrenko, A. Sidorenko, Low Temp. Phys. (2011) 451. J. Clarke, Phys. Rev. Lett. 28 (1972) 1363. M. Tinkham, J. Clarke, Phys. Rev. Lett. 28 (1972) 1366. M. Tinkham, Phys. Rev. B 6 (1972) 1747. A.D. Smith, M. Tinkham, W.J. Skocpol, Phys. Rev. B 22 (1980) 4346.

12 [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

D.-V. Anghel / Physica A 531 (2019) 121804 A.D. Smith, W.J. Skocpol, M. Tinkham, Phys. Rev. B 21 (1980) 3879. J.E. Hirsch, Phys. Rev. Lett. 72 (1994) 558. J.E. Hirsch, Phys. Rev. B 58 (1998) 8727. J.E. Hirsch, Phys. Scr. 88 (2013) 035704. J.E. Hirsch, F. Marsiglio, Phys. Rev. B 39 (1989) 11515. F. Marsiglio, J.E. Hirsch, Hole superconductivity and the high-tc oxides, Phys. Rev. B 41 (1990) 6435. R.V. Parfeniev, V.I. Kozub, G.O. Andrianov, D.V. Shamshur, A.V. Chernyaev, N. Yu. Mikhailin, S.A. Nemov, Low Temp. Phys. 41 (2015) 112. F. Bouquet, R.A. Fisher, N.E. Phillips, D.G. Hinks, J.D. Jorgensen, Phys. Rev. Lett. 87 (2001) 047001. P. Szabó, P. Samuely, J. Kačmarčík, T. Klein, J. Marcus, D. Fruchart, S. Miraglia, C. Marcenat, A.G.M. Jansen, Phys. Rev. Lett. 87 (2001) 137005. K. Tanaka, W.S. Lee, D.H. Lu, A. Fujimori, T. Fujii, Risdiana, I. Terasaki, D.J. Scalapino, T.P. Devereaux, Z. Hussain, Z.-X. Shen, Science 314 (5807) (2006) 1910. T. i Kondo, T. Takeuchi, A. Kaminski, S. Tsuda, S. Shin, Phys. Rev. Lett. 98 (2007) 267004. R.W. Richardson, Phys. Lett. 3 (6) (1963) 277. R.W. Richardson, N. Sherman, Nucl. Phys. 52 (1964) 221. R.W. Richardson, J. Math. Phys. 9 (9) (1968) 1327.