ICFP M2 Advanced Quantum Mechanics:

Homework Problem:
Fermi gas in the BCS regime
Solutions to selected questions
D.J. Papoular
December 10, 2019

1 Question 3: second–quantised expression for |ΨN i

We consider the wavefunction |ΨN i representing a well–defined number N of atoms which are
paired into N/2 Cooper pairs:
ΨN (r1 , . . . , rN ) = A [φ(r1 , r2 ) |χ12 i , . . . φ(rN −1 , rN ) |χN −1,N i] , (1)
where the antisymmetriser A acts on the fermionic positions r1 , . . . , rN , φ(ri , rj ) is the spatial
wavefunction of the pair involving
√ the fermions i and j. We choose the two–fermion internal
state |χij i = (|↑↓i − |↓↑i)/ 2 to be the singlet state so that it is antisymmetric, which ensures
that the spatial wavefunction φ(ri , rj ) is symmetric and, hence, allows for s–wave pairing. The
wavefunction |ΨN i is antisymmetric and normalised to unity:
Z
hΨN |ΨN i = d3 r1 . . . d3 rN |ΨN (r1 , . . . , rN )|2 = 1 . (2)

In order to write |ΨN i in second–quantised form, we introduce the pair creation operator b† :
Z
b† = d3 r1 d3 r2 φ(r1 , r2 ) Ψ̂†↑ (r1 )Ψ̂†↓ (r2 ) , (3)

where the field operator Ψ̂σ (r) creates a particle with the spin σ =↑ or ↓ at the point r.

1.1 A single Cooper pair

We first consider the action of b† on the vacuum state |vaci:
Z Z
† † †
b |vaci = d r1 d r2 φ(r1 , r2 ) Ψ̂↑ (r1 )Ψ̂↓ (r2 ) |vaci = d3 r1 d3 r2 φ(r1 , r2 ) |r1 ↑, r2 ↓i .
3 3
(4)

The ket |r1 ↑, r2 ↓i appearing in the right–hand side of Eq. (4) is written in second–quantised
notation. It is a two–fermion state which is properly antisymmetrised, and whose first–quantised
expression is:
1
|r1 ↑, r2 ↓i = √ (|1 : r1 ↑, 2 : r2 ↓i − |1 : r2 ↓, 2 : r1 ↑i) . (5)
2
Replacing Eq. (5) into Eq. (4), we obtain:
Z
† 1
b |vaci = d3 r1 d3 r2 φ(r1 , r2 ) √ (|1 : r1 ↑, 2 : r2 ↓i − |1 : r2 ↓, 2 : r1 ↑i) . (6)
2

1
We exchange the role of the integration variables r1 and r2 in the term involving |1 : r2 ↓, 2 : r1 ↑i
Recalling that φ(r1 , r2 ) = φ(r2 , r1 ), we obtain:
Z
1
b† |vaci = d3 r1 d3 r2 φ(r1 , r2 ) √ (|1 : r1 ↑, 2 : r2 ↓i − |1 : r1 ↓, 2 : r2 ↑i) (7)
2
Z 
= d3 r1 d3 r2 φ(r1 , r2 ) |1 : r1 , 2 : r2 i |χ12 i , (8)

where we have factorised the spatial and internal–state parts of the wavefunction. Hence, b† |vaci
is the properly normalised and antisymmetric wavefunction for a single Cooper pair.

1.2 N/2 Cooper pairs

We now wish to write the wavefunction |ΨN i in terms of b† . This wavefunction consists of N
atoms and, hence, of N/2 Cooper pairs. Therefore, we consider:
Z
†N/2
b |vaci = d3 r1 . . . d3 rN φ(r1 , r2 ) . . . φ(rN −1 , rN ) |r1 ↑, r2 ↓, . . . , rN −1 ↑, rN ↓i , (9)

where the ket |r1 ↑, r2 ↓, . . . , rN −1 ↑, rN ↓i is written in second–quantised notation, so that the

indices on the variables r1 , . . . , rN have nothing to do with particle ordering. The wavefunction
b†N/2 |vaci is a superposition of antisymmetrised wavefunctions, hence, it is itself antisymmetrised.
Therefore, it is equal to |ΨN i up to a normalisation factor.
We now calculate the norm of b†N/2 |vaci:
Z Z
hvac|bN/2 b†N/2 |vaci = d3 r1 . . . d3 rN d3 r10 . . . d3 rN
0

φ∗ (r10 , r20 ) . . . φ∗ (rN

0 0
−1 , rN ) φ(r1 , r2 ) . . . φ(rN −1 , rN )
hr10 ↑, r20 ↓, . . . , rN
0 0
−1 ↑, rN ↓ | r1 ↑, r2 ↓, . . . , rN −1 ↑, rN ↓i . (10)

The scalar product between the two N –particle states on the last line of Eq. (10) is non–zero if
the pairs appearing in the bra h(r10 , r20 ), . . . , (rN0 0
−1 , rN )| are the same as those appearing in the
ket |(r1 , r2 ), . . . , (rN −1 , rN )i up to a permutation of the pairs (and not of the particles). There
are (N/2)! such transpositions. For each of them, the scalar product is equal to the product of N
Dirac peaks δ(r1 − r10 ) · · · δ(rN − rN 0
). Therefore, Eq. (10) reduces to:
Z Z
N/2 †N/2
hvac|b b |vaci = (N/2)! d r1 d r2 |φ(r1 , r2 )| · · · d3 rN −1 d3 rN |φ(rN −1 , rN )|2 = (N/2)!
3 3 2

p (11)
Thus, the norm of b†N/2 |vaci is (N/2)! , and |ΨN i is expressed in terms of b† as follows:

b†N/2
|ΨN i = p |vaci , (12)
(N/2)!
where the numerical prefactor is the same as for a true bosonic creation operator.

2 Questions 7–9: the BCS wavefunction

2.1 The many–particle Hilbert space
We consider the BCS wavefunction |ΨBCS i, defined by Eq. 6 in the problem set:
1 p 
|ΨBCS i = exp Np b† |vaci , (13)
N
where Np = N/2 is the average number of pairs, which is half the average number of particles
† †
N = hN̂ i, b† =
P
k φk ck↑ c−k↓ is the creation operator for a Cooper pair, and the (real and

2
positive) number N is a normalisation factor which ensures that hΨBCS |ΨBCS i = 1. The state
|ΨBCS i does not have a well–defined particle number, hence, the Hilbert space H it belongs to is
best described in the language of second quantisation.
The Hilbert space H is usually decomposed into subspaces Hn with well–defined particle num-
bers n: M
H = H0 ⊕ H 1 ⊕ H 2 ⊕ · · · = Hn . (14)
n

In Eq. (14), H is written in terms of the subspaces Hn as a direct sum, because if |ψn1 i is a valid
n1 –particle quantum state (i.e. an element of Hn1 ) and |ψn2 i is a valid n2 –particle quantum state
(i.e. an element of Hn2√), then both |ψn1 i and |ψn2 i are valid elements of H, and the superposition
state (|ψn1 i + |ψn2 i)/ 2 is also a valid quantum state.
In the context of BCS theory, the decomposition of H in terms of particle numbers is not the
most adequate one: it is more convenient to decompose H in terms of spatial modes. We first
consider the plane–wave state |ki. Because of Pauli’s exclusion principle, there are only four licit
fermionic states that may be constructed from this spatial state, which may be non–populated
(|vack i), or populated by a single particle (|k ↑i or |k ↓i), or populated by two particles (|k ↑, k ↓i).
We call Hk0 the four–dimensional Hilbert space spanned by these four states. Then, the complete
many–body Hilbert space H is the tensor product of the spaces Hk0 over all plane–wave indices k:
O
H= Hk0 . (15)
k

This new decomposition of H involves tensor products and not direct sums, because it involves
independent spatial modes which may be independently populated. For instance, let us consider
two different plane–wave states |k1 i and |k2 i. Then, |k1 ↑i is a single–mode state, which is an
element of Hk0 1 , and |k2 ↑, k2 ↓i is another single–mode state, which is an element of Hk0 2 ; these two
states may be combined to obtain an element of H which is |k1 ↑i⊗|k2 ↑, k2 ↓i = |k1 ↑, k2 ↑, k2 ↓i,
and which contains two populated modes1 .
The most convenient decomposition of H in the context of BCS theory is a variant of Eq. (15).
We wish the elementary Hilbert spaces Hk to encode the presence or the absence of Cooper pairs,
each of which consists of one atom in the state |k ↑i and another in the state |−k ↓i. Hence,
instead of using the spaces Hk0 , we introduce the four–dimensional spaces Hk spanned by the
states |vack↑,−k↓ i, |k ↑i, |−k ↓i, and |k ↑, −k ↓i. Here, the vacuum state |vack↑,−k↓ i should be
understood as the absence of any particle in the states |k ↑i and |−k ↓i. We may finally write:
O
H= Hk . (16)
k

2.2 Factorisation of the BCS wavefunction

The commutator [c†k1 ↑ c†−k1 ↓ , c†k2 ↑ c†−k2 ↓ ] = 0 for any two plane waves |k1 i and |k2 i. Therefore, the
exponential of the sum of operators in Eq. (13) may be written as a product of exponentials:
p 
Np φk c†k↑ c†−k↓ |vaci .
Y
N |ΨBCS i = exp (17)
k

Now, we recall that c†k↑ and c†−k↓ are fermionic operators, so that c†2 †2
k↑ = 0 and c−k↓ = 0. The
exponential appearing in Eq. (17) is a series involving the integer powers (c†k↑ c†−k↓ )n , which are all
zero for n ≥ 2. Therefore, one may replace each exponential in Eq. (17) by the only two surviving
terms, which leads to:
Y 
1 + Np φk c†k↑ c†−k↓ |vaci .
p
N |ΨBCS i = (18)
k

1 One
√ √
may also consider the superposition (|k1 ↑i+|k2 ↑, k2 ↓i)/ 2 = (|k1 ↑i⊗|vack2 i+|vack1 i⊗|k2 ↑, k2 ↓i)/ 2.

3
 Now, we shall factorise |ΨBCS i into a form which is compatible with the decomposition of
Eq. (16). The vacuum state |vaci appearing in Eq. (18) is the product of all vacuum states
|vack↑,−k↓ i: Y
|vaci = |vack↑,−k↓ i . (19)
k

Furthermore, for two different wavevectors k1 and k2 and any two complex numbers αk1 and αk2 ,
  
1 + αk1 c†k1 ↑ c†−k1 ↓ 1 + αk2 c†k2 ↑ c†−k2 ↓ |vack1 ↑,−k1 ↓ i |vack2 ↑,−k2 ↓ i (20)
  
= 1 + αk1 c†k1 ↑ c†−k1 ↓ |vack1 ↑,−k1 ↓ i |vack2 ↑,−k2 ↓ i + αk2 |vack1 ↑,−k1 ↓ i c†k2 ↑ c†−k2 ↓ |vack2 ↑,−k2 ↓ i
(21)
h  i h  i
= 1 + αk1 c†k1 ↑ c†−k1 ↓ |vack1 ↑,−k1 ↓ i 1 + αk2 c†k2 ↑ c†−k2 ↓ |vack2 ↑,−k2 ↓ i . (22)

Therefore, Eq. (18) may be rewritten as:

Y h  i
1 + Np φk c†k↑ c†−k↓ |vack↑,−k↓ i .
p
N |ΨBCS i = (23)
k
Q
Equation (23) is of the form |ΨBCS i = k |Ψk i, where each |Ψk i is an element of Hk given by:

1  p  Y
|Ψk i = |vack↑,−k↓ i + Np φk |k ↑, −k ↓i , with N = Nk . (24)
Nk
k

p to normalise |ΨBCS i, it is sufficient to impose hΨk |Ψk i = 1 for all k. This leads to
In order
Nk = 1 + Np |φk |2 , so that we may finally write:
Y
|ΨBCS i = (uk |vack↑,−k↓ i + vk |k ↑, −k ↓i) , (25)
k

where uk and vk read:

p
1 Np φ k
uk = p and vk = p . (26)
1 + Np |φk |2 1 + Np |φk |2

2.3 Interpretation of the BCS wavefunction

Q
The factorisation |ΨBCS i = k |Ψk i provided by Eq. (25) allows for a simple physical interpre-
tation of the BCS wavefunction, which may be understood as follows. It describes a collection of
independent subsystems labelled by the wavevectors k, and whose wavefunctions are |Ψk i. Each
of these independent systems is in a superposition state involving the state |vack↑,−k↓ i (mean-
ing that the Cooper pair |k ↑, −k ↓i is absent) with the probability amplitude uk , and the state
|k ↑, −k ↓i (the Cooper pair is present) with the probability amplitude vk .
The factorisation described above is only possible with the BCS state |ΨBCS i. It does not hold
for the state |ΨN i with a fixed total number of particles. This is why |ΨBCS i is more amenable
to calculations and interpretations than |ΨN i.

The coefficients uk and vk may be determined exactly in the case of vanishing interactions
(a = 0). In this limit, the N –particle system consists of a Fermi sea defined by the Fermi
wavevector kF = (3π 2 n)1/3 , where n is the spatial density (cf. Problem 2, the different numerical
prefactor reflects the fact that all single–particle states may now be populated by two particles
with the internal states |↑i and |↓i). Then, below the Fermi surface, all single–particle states are
occupied, so that uk = 0 and vk = 1; above the Fermi surface, all single–particle states are empty,
meaning that uk = 1 and vk = 0.

4
2.4 Statistics of the number of particles in the BCS state
In order to determine the statistics of the total atom number N̂ in this state, we follow a procedure
similar to the one we employed in Problem 1 when calculating the photon or electron number
statistics in two–particle interference experiments. More precisely, we first calculate the properties
of the atom number n̂k = n̂k↑ + n̂−k↓ for a given wavevector k, and then use the independence of
P
the random variables nk to conclude as to the properties of N̂ = k n̂k .
The random variable nk takes the value 0 with probability |uk |2 , and the value 2 with the
probability |vk |2 (the value is 2 and not 1 because each Cooper pair contains two particles).
Therefore, its mean value is hn̂k i = 2|vk |2 . Similarly, the random value n2k takes the value 0 with
the probability |uk |2 , and the value 22 = 4 with the probability |vk |2 , so that its mean value is
hn̂2k i = 4|vk |2 . Therefore, the variance of the random variable nk is:
2
∆n2k = hn2k i − hnk i = 4(|vk |2 − |vk |4 ) = 4|uk |2 |vk |2 , (27)
2 2
P the normalisation condition |uk | + |vk | = 1.
where the last step follows from
Now, we recall that N̂ = k n̂k . The random variables nk are independent, so that the mean
2
value N = hN̂ i of the particle number and its variance ∆N 2 = hN̂ 2 i − hN̂ i satisfy:
X X X X
hN̂ i = hnk i = 2 |vk |2 and ∆N 2 = ∆n2k = 4 |uk |2 |vk |2 . (28)
k k k k

Finally, we show that for large particle numbers N , the fluctuation ∆N of the particle number is
small. For each wavevector k, the normalisation |uk |2 +|vk |2 = 1 entails that |uk |2 ≤ 1. Combining
this inequality with the expression for ∆N 2 in Eq. (28), we obtain:
s
2
X
2 ∆N 2
∆N ≤ 4 |vk | = 2 hN̂ i , so that ≤ . (29)
k
hN̂ i hN̂ i

The second inequality in Eq. (29) shows that ∆N/N becomes arbitrarily small for sufficiently large
N . Therefore, for large mean atom numbers, the fluctuations of the atom number in the BCS
state are negligible, meaning that in the large–N limit the state |ΨBCS i should capture the physics
of a system containing a fixed number of particles (i.e. described by the wavefunction |ΨN i).

3 Quest. 16: BCS ground state through a variational approach

Section 2.2 of the problem presents an analysis of the BCS Hamiltonian HBCS using the Bogoliubov
approach, which involves three steps. First, one exploits off–diagonal long–range order to replace
the quartic Hamiltonian HBCS by an approximate quadratic Hamiltonian. Second, Bogoliubov
rotations, parametrised by the coefficients uk and vk , are applied to the operators ck↑ and c†−k↓ to
diagonalise this quadratic Hamiltonian. Third, the coefficients uk and vk are interpreted as those
appearing in Eq. (25), i.e. they define the Fourier coefficient φk of the pair wavefunction φ(r1 , r2 ).
In the case of a uniform Fermi gas, a more direct approach is possible, which relies on a varia-
tional argument. In this approach, one chooses the coefficients uk and vk appearing in Eq. (25) so
as to minimise the grand–canonical energy hΨBCS |HBCS − µN̂ |ΨBCS i. This variational approach
yields the same values for uk and vk as the Bogoliubov approach presented in the problem.

We start from the expression for HBCS which appears as Eq. 10 in the problem. Substracting
from it the operator µN̂ , we obtain the grand–canonical Hamiltonian:
  g X
(k − µ) c†k↑ ck↑ + c†−k↑ c−k↑ + c† c† ck0 ↓ c−k0 ↑ .
X
HBCS − µN̂ = (30)
Ω 0 k↑ −k↓
k k,k

In the spirit of the Hilbert space decomposition of Eq. (16), we replace the sum over k0 by a sum
over −k0 . We also split the double sum over k and k0 according to whether the two wavevectors

5
are equal or different, and exploit the anticommutation rules in the case of equal wavevectors:
  g X g X † †
(k −µ) c†k↑ ck↑ + c†−k↑ c−k↑ + c†k↑ ck↑ c†−k↓ c−k↓ +
X
HBCS −µN̂ = ck↑ c−k↓ c−k0 ↓ ck0 ↑ .
Ω Ω 0
k k k6=k
(31)
The first two sums in Eq. (31) do not mix the subspaces Hk of Eq. (16). In order to evaluate their
average, we split the wavefunction |ΨBCS i into two terms |Ψk0 i and |Ψk1 i, corresponding to the
cases where the pair (k ↑, −k ↓) is absent or present, respectively. Thanks to Eq. (25):

|ΨBCS i = uk |Ψk0 i + vk |Ψk1 i . (32)

Using this decomposition, one finds that the average of the kinetic energy term related to the
wavevector k is 2(k − µ)vk2 , and the average of the interaction term for k = k0 is (g/Ω)vk2 .
The third sum in Eq. (31) mixes two subspaces Hk and Hk0 with k 6= k0 . In order to evaluate its
average, we split the wavefunction |ΨBCS i into four terms |Ψk0k0 0 i, |Ψk0k0 1 i, |Ψk1k0 0 i, |Ψk1k0 1 i,
describing the absence or the presence of the two pairs (k ↑, −k ↓) and (k0 ↑, −k0 ↓):

|ΨBCS i = uk uk0 |Ψk0k0 0 i + uk vk0 |Ψk0k0 1 i + vk uk0 |Ψk1k0 0 i + vk vk0 |Ψk1k0 1 i . (33)

Using this decomposition, one finds that the average of the interaction term corresponding to
k 6= k0 is (g/Ω)uk vk uk0 vk0 . Combining the averages for all three sums in Eq. (31), we obtain:
X g 2 g X
hΨBCS |HBCS − µN̂ |ΨBCS i = k − µ + vk + u k u k 0 vk vk 0 . (34)
Ω Ω 0
k k6=k

The interaction terms for k = k0 shift the chemical potential by |g|/Ω. In the BCS limit, this shift
is negligible compared to µ ∼ EF , so that the average grand–canonical energy finally reduces to:
X g X
hΨBCS |HBCS − µN̂ |ΨBCS i = (k − µ) vk2 + u k u k 0 vk vk 0 . (35)
Ω 0
k k6=k

In order to minimise Eq. (35), we recall that u2k + vk2 = 1. Hence, we let uk = cos θk and
vk = sin θk , and write Eq. (35) in terms of the angles θk :
1X 1g X
hΨBCS |HBCS − µN̂ |ΨBCS i = (k − µ) [1 − cos(2θk )] + sin(2θk ) sin(2θk0 ) . (36)
2 4Ω 0
k k6=k

Derivating Eq. (36) with respect to θk , we obtain the minimisation conditions:

1g X
(k − µ) tan(2θk ) = − sin(2θk0 ) . (37)
2Ω 0
k 6=k

We introduce the gap parameter ∆, defined as an average in the state |ΨBCS i:

g X g X 1 g X
∆= hΨBCS |ck↑ c−k↓ |ΨBCS i =− uk vk =− sin(2θk ) , (38)
Ω Ω 2 Ω
k k k

where the second equality follows from Eq. (32). Combining Eqs. (37) and (38):
 
1g
k − µ − cos(2θk ) tan(2θk ) = ∆ . (39)
2Ω

We again neglect |g|/Ω compared to µ on the left–hand side, and obtain:

∆
tan(2θk ) = , (40)
k − µ

6
from which we extract the absolute values cos(2θk ) and | sin(2θk )|:

|k − µ| |∆| p
| cos(2θk )| = and | sin(2θk )| = , where Ek = (k − µ)2 + ∆2 . (41)
Ek Ek
We choose the sign of cos(2θk ) to ensure the convergence of the average total number of particles:

hΨBCS |c†k↑ ck↑ + c†−k↓ c−k↓ |ΨBCS i =

X X
N= 2vk2 , (42)
k k

where the last step follows from the decomposition of Eq. (32). In order for the sum of Eq. (42)
to converge, the coefficients vk must go to zero for large k = |k|. Therefore,

For large k, cos(2θk ) = cos2 θk − sin2 θk = u2k − vk2 > 0 . (43)

We expect ∆ > 0; additionally, for large k, k − µ > 0. Hence, cos(2θk ) and sin(2θk ) satisfy:

k − µ ∆
cos(2θk ) = + and sin(2θk ) = + . (44)
Ek Ek

Using cos2 θk = [1 + cos(2θk )]/2, sin2 θk = [1 − cos(2θk )]/2, and cos θk sin θk = sin(2θk )/2:
 
2 2 1 k − µ ∆
u k = 1 − vk = 1+ and uk vk = , (45)
2 Ek 2Ek

in accordance with Eq. 18 in the problem.

The present variational method is more direct than the Bogoliubov approach outlined in the
problem. However, it is less general than the Bogoliubov approach, for two reasons. First, not
only does the Bogoliubov approach yield the ground state, it also shows that the energies Ek
correspond to the excitation spectrum. Second, the Bogoliubov approach may be applied in an
inhomogeneous system, i.e. in the presence of a trapping potential V (r) (see e.g. Ref. [1, chap. 5]).

4 Quest. 19–27: chemical potential and gap in terms of a and N

The BCS model is parametrised by two key quantities: the chemical potential µ and the gap ∆.
Our last step is to relate µ and ∆ to two parameters which may be tuned experimentally: the
scattering length a and the particle number N . This is done by simultaneously solving the gap
equation and the number equation, which appear together as Eq. 22 in the problem set.
In the BCS limit, the interaction between the particles is attractive and arbitrarily weak. This
leads to a gap parameter ∆ which is positive and very small. We shall make this assumption
throughout this section and check its validity at the end of the calculation. In particular, we shall
evaluate the integrals appearing in the gap and number equations to leading order in ∆/EF .

4.1 Chemical potential

We start from Eq. (42) above, which gives the total number of particles. Replacing vk2 by its value
given by Eq. (45):
X d3 k
  
k − µ k − µ
Z
N= 1− = 1− . (46)
Ek (2π)3 /Ω Ek
k

The gas is homogeneous, therefore its spatial density is n = N/Ω:

d3 k
 
k − µ
Z
n= 1 − . (47)
(2π)3 Ek

7
In Eq. (47), the integrand f (k ) is a function of k = ~2 k 2 /(2m) only. Introducing the Fermi
energy EF and the Fermi wavevector kF defined by EF = ~2 kF2 /(2m), we change variables from
d3 k to d in the integral using the general formula:
1/2
d3 k kF3
Z Z 
d 
f (k ) = f () . (48)
(2π)3 4π 2 EF EF
Applying Eq. (48) to Eq. (47), we find:
1/2 
kF3
 
k − µ
Z
d 
n= 1− . (49)
4π 2 EF EF Ek

Replacing kF by its expression for a Fermi gas with two accessible internal states, kF = (3π 2 n)1/3 ,
the spatial density cancels out and we obtain:
Z ∞  1/2  
4 d  −µ
= 1− . (50)
3 0 EF EF Ek

We consider the integral appearing in Eq. (50) in the limit where ∆ → 0. The excitation energy
 1/2
Ek = (k − µ)2 + ∆2 reduces to Ek = |k − µ|. The content of the second pair of parentheses
in Eq. (50) is equal to 2 for  < µ, and equal to 0 for  > µ. Therefore, Eq. (50) reduces to:
Z µ  1/2  3/2
4 d  4 µ
= 2 = . (51)
3 0 EF EF 3 EF
We conclude from Eq. (51) that µ = EF : In the BCS limit, the attractive interaction is so weak
that it does not affect the chemical potential, which retains its non–interacting value.

4.2 Gap parameter

4.2.1 The gap equation
We now turn to the gap parameter ∆, which we express in terms of the coefficients uk and vk
thanks to Eq. (38) above. We replace uk vk by its value given by Eq. (45):
g X g X ∆
∆=− uk vk = − . (52)
Ω Ω 2Ek
k k

The gap ∆ cancels out from Eq. (52) (although it is still present implicitly through Ek ):
1g X 1
1=− . (53)
2Ω Ek
k

We would now like to replace the discrete sum in Eq. (53) by an integral. This must be done with
care, because the integrand k 2 /k is not integrable for large values of k. Therefore, we introduce a
cutoff function η(k), which is equal to 1 for small values of k and goes to zero for large k sufficiently
rapidly to ensure the convergence of the integral:
d3 k η(k)
Z
1 1
=− . (54)
g 2 (2π)3 Ek
We compare this equation with the one linking the interaction strength g to the scattering length
a, which is established in the “Complement” to the problem (in particular, cf. question 33):

d3 k η(k)
Z
1 m 1
= − . (55)
g 4π~2 a 2 (2π)3 k

8
Substracting Eq. (54) from Eq. (55), we find:
d3 k
Z  
m 1 1
= η(k) − . (56)
4π~2 a (2π)3 2k 2Ek

In Eq. (56), the integrand involves k 2 (1/k − 1/Ek ), which behaves like k 2 µ/(22k ) ∼ 1/k 2 for
large values of k. Therefore, the integral converges even in the absence of the cutoff function, and
its limiting value does not depend on the shape of η(k), which can be dropped:
d3 k
Z  
m 1 1
= − . (57)
4π~2 a (2π)3 2k 2Ek

Finally, we change variables from d3 k to d in the integral thanks to Eq. (48), and use EF =
~2 kF2 /(2m) to find the adimensional gap equation:
Z ∞  1/2  
d  EF EF π
− =− . (58)
0 EF EF Ek  kF a

4.3 Shifting the zero–energy singularity to the Fermi surface

The extraction from Eq. (58) of the gap ∆ in terms of kF a requires a few technical manipulations.
These are guided by the intuition that the Fermi surface plays a key role, an idea which is confirmed
by the result above that EF = µ. Therefore, we expect the main contribution to the integral in
Eq. (58) to come from energies of the order of µ.
In order to make this √contribution
√ as explicit as possible, we wish to split the first factor in the
√ √
integral√into two parts:  = (  − µ) + µ. However, this is not directly possible, because the
factor  compensates for the divergence of the factor (1/Ek − 1/) for small .
Therefore, our first step is to shift the small– divergence to the Fermi surface. In Eq. (58), we
would like to replace the second factor (1/Ek − 1/) by (1/Ek − 1/( − µ)). For that purpose, we
wish to give a meaning to the following undefined integral:
Z ∞  
1 1  µ
dx x1/2 − , where x = and α = . (59)
0 x−α x EF EF

This integral is reminiscent of those which involve Cauchy’s principal value2 . This resemblance
motivates its calculation using the formalism of complex analysis. Instead of integrating from
0 to +∞ along the positive real axis, we integrate along the closed contour C in the complex
plane shown on Fig. 1. This contour avoids both the pole at x = α and the branch cut for the
square–root, which coincides with the negative real axis. Hence, we calculate:
Z  
1 1
IC = dz z 1/2 − = 0, (60)
C z−α z
where the last step follows from the fact that the closed contour C circles no poles of the integrand.

We now split the integral of Eq. (60) in three pieces: 0 = IC = I1 + I2 + I3 , corresponding to

the three parts of the contour C:
• The contribution I1 coming from the (bumpy) horizontal part of the contour is:
Z ∞  
1/2 1 1
I1 = dx x − , (61)
0 x − α + i0+ x
where the small imaginary part i0+ shifts the pole at x = a below the horizontal axis.
2 Cauchy’s
R∞
principal value is used to give a meaning to undefined integrals of the form: −∞ dx f (x)/(x − x0 ). It is
frequently encountered in Physics: for instance, it plays a key role in establishing the Kramers–Krönig relations,
which link the real and imaginary parts of the susceptibility in linear response theory. An accessible review of
the definition and physical applications of Cauchy’s principal value may be found in [2, § 8.1].

9
 y

large radius R
Figure 1 Closed integration contour C used in
Eq. (60). The contour avoids both the pole at x = a
and the branch cut for the square root, which coin-
cides with the negative real axis.
x
branch cut 0 α

• The integrand behaves like 1/z 3/2 for √

large |z|, so that the contribution of the large–radius
quarter–circle is bounded: I2 ≤ π/(2 R). Hence, I2 goes to zero for large values of the
radius R.
• The contribution coming from the imaginary
√ axis
p may be calculated
p exactly using the prim-
itive dt t1/2 (1/(t − u) − 1/t) = − u ln[(1 + t/u)/(1 − t/u)]:
R

Z ∞
√
 
1/2 1 1
I3 = − idy (iy) − = +iπ α . (62)
0 iy − α iy

Bringing all three pieces together, we obtain I1 = −I3 , that is:

Z ∞
√
 
1/2 1 1
dx x +
− = −iπ α . (63)
0 x − α + i0 x

Combining Eq. (63) with the gap Eq. (58), we finally obtain:
Z ∞  
π 1/2 1 1
− = Re dx x − , where δ = ∆/EF . (64)
kF a 0 [(x − α)2 + δ 2 ]1/2 x − α + i0+

4.4 Evaluation of the integral in the gap equation

√ √
Now that the zero—energy singularity has been eliminated, it is possible to split  = (  −
√ √
µ) + µ in order to make the contribution of the Fermi surface explicit, as suggested at the
beginning of Sec. 4.3. Hence, we write −π/(kF a) = Re[J1 + J2 ], where J1 and J2 are defined as:
Z ∞  
1 1
J1 = dx (x1/2 − α1/2 ) − , (65)
0 [(x − α)2 + δ 2 ]1/2 x − α + i0+
Z ∞  
1/2 1 1
J2 = dx α − . (66)
0 [(x − α)2 + δ 2 ]1/2 x − α + i0+

First, we calculate J1 . The integral appearing in Eq. (65) is well–behaved in the limit ∆ → 0
(that is, δ → 0). We find:
∞ a
 a1/2 − x1/2
Z  Z
 1 1
J1 = dx x1/2 − a1/2 − = −2 dx . (67)
0 |x − a| x − a 0 a−x

The integral in Eq. (67) may be evaluated exactly:

Z a
1 √
J1 = −2 dx √ √ = −4 a (1 − ln 2) . (68)
0 a+ x

10
 Now, we turn to the calculation of J2 . Both terms appearing in the integral of Eq. (66) have
explicit primitives, and we find:

∞
     
1 x−a 1 x−a +
J2 = ln 1 + − ln 1 − − ln x − a + i0 .
2 [(x − a)2 + δ 2 ]1/2 2 [(x − a)2 + δ 2 ]1/2 0
(69)

Let us examine the two bounds:

• For x → +∞, the first term inside the straight brackets goes to (ln 2)/2. The second
and third term diverge when considered separately, but their sum is finite and equal to
− ln(δ) + (ln 2)/2.
• For x → 0, the third term inside the straight brackets goes to − ln(a) + iπ. The second term
also goes to a finite value, which is −(ln 2)/2 if δ  a. The first term goes to −(ln 2)/2 −
ln a + ln δ.
We recall that µ = EF (see Sec. 4.1 above). Therefore, a = µ/EF = 1, and J2 = −2 ln(δ/2) − iπ.

Adding J1 and J2 together, we obtain:

 
−π ∆
= Re[−4(1 − ln 2) − 2 ln(δ/2) − iπ] = −4(1 − ln 2) − 2 ln , (70)
kF a 2EF

which finally leads to:  

∆ 8 π
= 2 exp . (71)
EF e 2kF a
Equation (71) is the famous result of BCS theory which expresses the gap parameter ∆ in terms
of the scattering length a, and more precisely in terms of the adimensional parameter kF a. This
equation validates the two assumptions we had made on ∆ throughout the calculation: First,
∆ > 0. Second, in the deep BCS regime where kF a is negative and small, the exponential is
extremely small, so that ∆  EF .
The prefactor in the exponential differs by a factor 1/2 from the the solution of the Cooper
problem (cf. Problem #4, Eq. 14). This is because the approach of the Cooper problem accounts
for the formation of a single Cooper pair on top of a Fermi sea wich remains unpaired, whereas
the full BCS treatment relies on the assumption that all N atoms are in a quantum superposition
between being paired and being unpaired (see Sec. 2.3 above). One of the conclusions of the
BCS treatment is that the many–body state defined by Eq. (25) is indeed energetically favoured
compared to the non–interacting ground state (cf. question 17 and Eq. 19 in the problem).

5 Question 29: Validity and experimental considerations

5.1 BEC–BCS crossover
This problem has explored the BCS model applied in the BCS limit, i.e. kF a negative and small.
However, this same model is also applicable in the deep BEC regime, i.e. kF a positive and small.
The BCS model is a mean–field theory3 , therefore it cannot be expected to provide an accurate
solution in the strongly–correlated regime where kF |a| & 1. Nevertheless, it provides useful orders
of magnitude, and in the context of this model all thermodynamic functions vary smoothly when
1/(kF a) is varied along the crossover from the BCS to the BEC regimes [3, §16.7].
Among these thermodynamic quantites, the critical temperature Tc below which the system
is superfluid is an important example. The order of magnitude for Tc is dictated by the gap:
3 Inthe Bogoliubov treatment outlined in the problem, the mean–field approximation intervenes when one replaces
the BCS Hamiltonian of Eq. (30), which is of order 4 in the operators ckσ and c†kσ , by the quadratic Hamiltonian
Heff (cf. question 12).

11
 Figure 2 Critical temperature for superfluidity along
the BEC–BCS crossover. The smooth solid line la-
belled “condensation” results from the BCS approach
discussed in the problem. The diamond indicates the
value of Tc /TF in the unitary regime obtained using
numerical quantum Monte Carlo calculations. Re-
produced from Ref. [3].

kB Tc ∼ ∆, with kB being Boltzmann’s constant. An improved approximation may be obtained

by solving the BCS equations at non–zero temperature [4, §4.8]. The resulting adimensional ratio
Tc /TF is represented as a function of 1/(kF a) on Fig. 2 (solid line labelled ‘condensation’).
The two weakly–interacting limits for the critical temperature Tc are well–understood. First,
for kF a small and negative, the BCS theory yields:
eγ
kB TcBCS = ∆(T = 0) , (72)
π
where γ ≈ 0.577 is Euler’s constant, and ∆(T = 0) is the zero–temperature value of the gap, which
is given by Eq. (71) and is exponentially small. Second, for kF a small and positive, the pairs are
tightly–bound bosonic molecules, and the critical temperature is given by the standard formula
for a uniform Bose gas, with the mass of a single boson being mmol = 2m and the bosonic density
being nmol = n/2:
nmol Λ3Tc = ζ(3/2) . (73)
In Eq. (73), ΛT = [2π~2 /(mmol kB T )]1/2 is the thermal de Broglie wavelength and ζ(3/2) =
1/n3/2 ≈ 2.61. Solving Eq. (73) for Tc , one obtains the critical temperature in the BEC limit:
P

1/3
kB TcBEC

2 1
= ≈ 0.218 . (74)
EF π [3ζ(3/2)]2

The investigation of the strongly–interacting regime is much more challenging. The prediction of
the BCS model is smooth but non–monotonic near the unitary limit where 1/(kF a) = 0 (see the
solid curve on Fig. 2). In the regime where |kF a| & 1, mean–field theories are expected to be
inaccurate, and at present no exact solution for the many–body problem is available. Therefore,
one turns to more involved numerical Quantum Monte Carlo calculations, which yield the following
universal4 result for 1/(kF a) = 0 (see e.g. Ref. [3, §16.6]):

kB Tcunitarity
≈ 0.167 , (75)
EF
which is in excellent agreement with recent experimental results [5], but is slightly below the
mean–field prediction, thus confirming that the latter is not sufficiently accurate.

4 The unitary Fermi gas is defined by two conditions: (i) kF r0  1, where r0 is the range of the interaction
potential, meaning that the system is dilute; and (ii) |kF a|  1, meaning that the system is strongly–interacting.
In this limit, the scattering length a is no longer a relevant lengthscale: the only remaining dimensional parameter
is the Fermi energy EF . Therefore, the adimensional ratio kB Tcunitarity /EF is a universal number.

12
 Figure 3 Experimental proof of superfluidity in an
interacting ultracold gas of fermionic 6 Li: Abrikosov
lattices of vortices have been nucleated and observed
for 1/(kF a) = 1.6 (BEC regime), 1/(kF a) = 0 (uni-
tary regime), and 1/(kF a) = −0.7 (BCS regime). Re-
produced from Ref. [4].

Figure 4 Abrikosov lattice of vortices observed by

scanning tunnelling microscopy in the superconduc-
tor NbSe2 , which exhibits BCS superconductivity for
T < Tc = 7.2 K. Reproduced from Ref. [6].

5.2 Observation of Abrikosov lattices

The observation of vortex lattices in interacting Fermi gases (see e.g. Ref. [4, Sec. 6.6]) has provided
spectacular experimental proof of their superfluid character. These lattices have been nucleated
for various values of a in the BEC regime (1/(kF a) = 1.6 on Fig. 3), at unitarity (1/(kF a) = 0),
and in the BCS regime (1/(kF a) = −0.7 on Fig. 3). For these values of the interaction parameter
1/(kF a), Eqs. (74) and (75) show that the critical temperature is of the order of a few tenths
of the Fermi energy, so that Tc ∼ 200 nK, a temperature which is nowadays routinely accessed
experimentally. However, bringing a Fermi gas deep into the BCS regime is an experimental
challenge because of the exponentially small critical temperature in the BCS limit (see Eq. (72)).
To the best of my knowledge, this challenge remains to be met. For now, the equivalent of the very
deep BCS regime on the far left of Fig. 2 is more easily accessed in condensed–matter systems, i.e.
metallic superconductors, where TcBCS is of the order of a few Kelvin. For example, Fig. 4 shows
a vortex lattice observed by scanning tunnelling microscopy in niobium diselenide, whose critical
temperature is Tc = 7.2 K, in the presence of a magnetic field.

References
[1] P. de Gennes, Superconductivity of metals and alloys, Addison-Wesley (1966).
[2] W. Appel, Mathematics for Physics and Physicists, Princeton University Press (2007).
[3] L. P. Pitaevskii, S. Stringari, Bose-Einstein condensation and superfluidity, Oxford University Press, 2nd ed.
(2016).
[4] W. Ketterle, M. Zwierlein, in Proceedings of the International School of Physics Enrico Fermi, course CLXIV
(Varenna, June 2006) (2008).
[5] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, M. W. Zwierlein, Science 335, 563 (2012).
[6] A. A. Abrikosov, Rev. Mod. Phys. 76, 975 (2004).

13