May 14th 2024 Quantum Statistics Duong Vu

Mathematical tips and tricks

The Gamma function Z ∞
Γ(n) = xn−1 e−x dx
0
satisfies
Γ(n) = (n − 1)!, Γ(n + 1) = nΓ(n),
and has special values √
Γ(1) = 1, Γ(1/2) = π
A common trick of calculating the average value of some quantity is to transform
the integral into the Gamma function (for example, let’s say that the quantity of
interest Q is proportional to xk , and the energy is proportional to xn )
R ∞ k −βxn R∞ k 1 k+1
n · u n −1 e−βu du


x e dx u Γ
⟨Q⟩ = 0R ∞ −βxn = 0 R ∞ 1 −1 = β −k/n n

0 e dx 0 u
n e −βu du Γ n1

The Riemann zeta function

∞
X 1
ζ(s) =
n=1
ns

has nice values for even s

π2 π4
ζ(2) = , ζ(4) =
6 90
and has approximate values for odd s

ζ(3) ≈ 1.202, ζ(5) ≈ 1.037

This function appears in the following integral (photon gas, Debye solid, etc.)
Z ∞ n−1
x dx
= ζ(n)Γ(n)
0 ex − 1
The polylogarithm function
∞
X zk
Lis (z) =
ks
k=1
May 14th 2024 Quantum Statistics Duong Vu

this function can be defined as follows

Z z
Lis (t)
Li−1 (z) = − ln(1 − z), Lis+1 (z) = dt
0 t

and it appears in the following integral (quantum fluid)

Z ∞
xn−1 dx
−1 βx
= (kB T )n Γ(n)[∓Lin (∓z)]
0 z e ±1
Another thing to note is

Lin (1) = ζ(n), Lin (z) ≈ z for z ≪ 1

Surface and Volume of N-d sphere

2π n/2 n−1 π n/2

Sn =
R , Vn =
Rn
Γ n2 n
Γ 2 +1
May 14th 2024 Quantum Statistics Duong Vu

The Grand Canonical Ensemble

Consider a system connected to a reservoir that can exchange energy and particles.
The probability of the system being in a particular microstate j with energy ϵ and
N particles is proportional to the number of microstates of the reservoir with
energy U − ϵ and N − N particles

pj (ϵ, N ) ∝ Γ(U − ϵ, N − N ) = eS(U −ϵ,N −N )/kB

Since N ≪ N and ϵ ≪ U , we can Taylor expand

   
∂S ∂S
S(U − ϵ, N − N ) = S(U, N ) − ϵ −N
∂U N ,V ∂N U,V

Recall that    
∂S 1 ∂S µ
= , =−
∂U N ,V T ∂N U,V T
Substitute this in, we get
pj (ϵ, N ) ∝ e−β(ϵ−µN )
As a result, the probability of our system being in a microstate j is

e−β(ϵj −µNj )
pj =
Z
where Z is the grand partition function
X
Z= e−β(ϵj −µNj )
j

Similar to the case of the partition function, where we have Z = e−βF , we can
define the grand potential ΦG

Z = e−βΦG , ΦG = −kB T ln Z

Immediately from the grand partition function, we can get some of the many
thermodynamic quantities
   
∂ ln Z ∂ ln Z
N = kB T , U =− + µN
∂µ β ∂β µ
May 14th 2024 Quantum Statistics Duong Vu

We can also get the entropy

X
S = −kB pj ln pj
j
X
= −kB pj [−β(ϵj − µNj ) − ln Z]
j
U − µN + kB T ln Z
=
T
Rearranging this equation, we get

ΦG = −kB T ln Z = U − T S − µN

Therefore,
dΦG = −SdT − pdV − N dµ
and      
∂ΦG ∂ΦG ∂ΦG
S=− , p=− , N =−
∂T V,µ ∂V T,µ ∂µ T,V
We expect all extensive variables to scale by a factor of λ when we scale the system
by a factor λ
U → λU, S → λS, V → λV, N → λN
Writing the entropy S as a function of U, V, N we get

λS(U, V, N ) = S(λU, λV, λN )

We can consider λ as a variable and differentiate both sides with respect to λ

∂S ∂λU ∂S ∂λV ∂S ∂λN
S= + +
∂λU ∂λ ∂λV ∂λ ∂λN ∂λ
which gives
1 p µ
S= ·U + ·V − ·N
T T T
(note how all of the coefficients are intensive quantities!). From here, we get

U − T S + pV = µN

Our left-hand side is just the Gibbs free energy, so therefore

G = µN
May 14th 2024 Quantum Statistics Duong Vu

This gives us a new interpretation of the chemical potential: it is the Gibbs free
energy per particle!
G
µ=
N
Note that from here, we also get that

ΦG = U − T S − µN = G − µN − pV = −pV
May 14th 2024 Quantum Statistics Duong Vu

Statistics of identical particles

The grand canonical ensemble is the most straightforward way of dealing with
identical particles. To specify a particular microstate of a system with identical
particles, we specify the number of particles in each quantum state (character-
ized by a set of quantum numbers). The grand canonical ensemble saves us the
headache of counting the number of ways to divide N particles into all of the quan-
tum states. As we will see, all quantum states will “act” as their own independent
“thing.” Therefore, we only need to consider the grand partition function for a
specific quantum state α with energy Eα (note that α will be used to differentiate
different quantum states, while j will be used to differentiate different numbers of
particles) X X
−β(nj Eα −µnj )
Zα = e = e−nj β(Eα −µ)
j j
For Fermions, the nj can only take up values 0 and 1, while for Bosons, nj can
have any nonnegative integer value. We can calculate the grand partition function
for both.
For Fermions,
Zα = 1 + e−β(Eα −µ)
while for Bosons,
1
Zα =
1 − e−β(Eα −µ)
It’s easier to summarize the result if we take the natural logarithm of the expression
ln Zα = ± ln[1 ± e−β(Eα −µ) ]
where the upper one is for Fermions, and the bottom one is for Bosons.
Of course, the grand canonical ensemble wouldn’t be good if it didn’t easily gener-
alize to multiple energy states. To show this clearly, we have to remember that α
will be used to differentiate between different quantum states, and j will be used
to differentiate between different numbers of particles in one particular quantum
state. For each microstate of the system (characterized by a set of nα ), the factor
that enters the sum is
h i h i Y
−n1 β(E1 −µ) −n2 β(E2 −µ)
e × e × ··· = e−nα β(Eα −µ)
α

Therefore, our grand partition function is obtained by summing over all the dif-
ferent sets of {nα } X Y
Z= e−nα β(Eα −µ)
all {nα } α
May 14th 2024 Quantum Statistics Duong Vu

Since we are summing over all sets of {nα }, we can use the fact that the sum of
products is the product of sums. This is such an important point that it is worth
reviewing why this is true for the simplest case. Let’s say we have the product of
two sums ! !
X X
product of sums = ai × bj
i j

We can imagine drawing a table to expand this product

a1 a2 ··· am
b1 a1 b1 a2 b1 ··· am b 1
b2 a1 b2 a2 b2 ··· am b 2
.. .. .. ... ..
. . . .
bn a1 bn a2 bn ··· am b n

The product of sums can also be written as the sum of all the table values. As a
result, X
product of sums = ai bj = sum of products
i,j

We can imagine how this would generalize to an arbitrary number of sums. Back
to our problem, we can rewrite the sum of products as
YX
Z= e−nj β(Eα −µ)
α j

where now each of the sums goes over all of the possible numbers of particles that
can be in a particular quantum state α. Therefore, the grand partition function
of our system is just the product of all of the grand partition functions of each of
the individual quantum states.
If we take the natural logarithm of the grand partition function, we get
X X
ln Z = ln Zα = ± ln[1 ± e−β(Eα −µ) ]
α α

From here, we can find the number of particles in quantum state α

1 ± −βe−β(Eα −µ) )
   
1 ∂ ln Z 1
⟨nα ⟩ = − =∓ =
β ∂Eα β 1 ± e−β(Eα −µ) eβ(Eα −µ) ± 1
We can define the distribution function f (E) for fermions and bosons. It is the
mean occupation of a particular quantum state with energy E. The distribution
May 14th 2024 Quantum Statistics Duong Vu

function for fermions is known as the Fermi-Dirac distribution function

1
f (E) =
eβ(E−µ) + 1
and the distribution function for bosons is known as the Bose-Einstein distri-
bution function
1
f (E) = β(E−µ)
e −1
We can go ahead and plot this

e−β(E−µ) Bose-Einstein
f (E)

1
Fermi-Dirac

0
−4 −2 0 2 4
β(E − µ)

Both functions tend to the Boltzmann distribution e−β(E−µ) in the limit of high
energy β(E − µ) ≫ 1 (needs more explanation).
Quoted directly from Blundell:“This is because this limit corresponds to low den-
sity (µ small) and here there are many more states thermally accessible to the
particles than there are particles; thus double occupancy never occurs and the
requirements of exchange symmetry become irrelevant and both fermions and
bosons behave like classical particles. The differences, however, are particularly
felt at high density. In particular, note that the distribution function for bosons
diverges when µ = E. Thus for bosons, the chemical potential must always be
below, even if only slightly, the lowest-energy state. If it is not, then the lowest-
energy state would become occupied with an infinite number of particles, which
is unphysical.”
May 14th 2024 Quantum Statistics Duong Vu

Non-interacting quantum fluid

We can consider a fluid of non-interaction particles with spin S. Note that when
we consider a quantum state in the previous section, the spin is one of the quantum
numbers that help specify that quantum state. That’s why the number of fermions
in the last section could only be 0 or 1; we treat different spins as different quantum
states. The grand partition for each of the spin states is identical, however, so the
resulting grand partition function is
Y
Z= Zk2S+1
k

where
Zk = [1 ± e−β(Ek −µ) ]±1
and the k is used to denote a set of quantum numbers specifying that quantum
state (regardless of the spin). In our case of a particle in a box, this k is the wave
vector, and
ℏ2 k2
Ek =
2m
In the following section, we will assume that the temperature is high enough so
that Bose-Einstein Condensation does not occur. The grand potential is

ΦG = −kB T ln Z
X
= ∓kB T (2S + 1) ln[1 ± e−β(Ek −µ) ]
k
Z ∞
= ∓kB T ln[1 ± e−β(E−µ) ]g(E) dE
0

where we approximated the integral as a sum, and g(E) is the density of state.
We can calculate this in k-space
1 4πk 2 dk (2S + 1)V k 2 dk
g(k) dk = × × (2S + 1) =
8 (π/L)3 2π 2
Using E = ℏk 2 /2m we get
√
(2S + 1)V E 1/2 dE 2m 1 2m
g(E) dE = × 2 · ·
2π 2 ℏ 2 ℏ
1/2
 3/2
(2S + 1)V E dE 2m
=
(2π)2 ℏ2
May 14th 2024 Quantum Statistics Duong Vu

From here, we can evaluate the grand potential

 3/2 Z ∞
(2S + 1)V 2m
ΦG = ∓kB T 2 2
ln[1 ± e−β(E−µ) ]E 1/2 dE
(2π) ℏ 0
 3/2  ∞ Z ∞ −β(E−µ)

(2S + 1)V 2m 2 3/2 2 ∓βe
= ∓kB T 2 2
E ln[1 ± e−β(E−µ) ] − E 3/2 −β(E−µ)
dE
(2π) ℏ 3 0 0 3 1 ± e
 3/2 Z ∞ 3/2
2 (2S + 1)V 2m E dE
=− 2 2 β(E−µ)
3 (2π) ℏ 0 e ±1
We also have Z ∞
X g(E) dE
N= nk =
k 0 eβ(E−µ) ± 1
and Z ∞
X E g(E) dE
U= nk Ek =
k 0 eβ(E−µ) ± 1
We can write eβµ = z, where z is called the fugacity. With this, the expressions
for N and U are as follow
" 3/2 # Z ∞
E 1/2 dE

(2S + 1)V 2m
N=
(2π)2 ℏ2 0 z −1 eβE ± 1
and " 3/2 # Z ∞
E 3/2 dE

(2S + 1)V 2m
U=
(2π)2 ℏ2 0 z −1 eβE ± 1
The integrals can be written in terms of the polylogarithm function, where
Z ∞ n−1
E dE
−1 βE
= (kB T )n Γ(n)[∓Lin (∓z)]
0 z e ±1
√
Using this and rewriting things in terms of the thermal wavelength λth = h/ 2πmkB T
to clean up our factors, we get
(2S + 1)V
N= [∓Li3/2 (∓z)]
λ3th
and
3 (2S + 1)V
U = kB T [∓Li5/2 (∓z)]
2 λ3th
3 Li5/2 (∓z)
= N kB T
2 Li3/2 (∓z)
May 14th 2024 Quantum Statistics Duong Vu

Note that we also have

2
ΦG = − U
3
For very small z (eβµ ≪ 1), the polylogarithm can be approximated as

Lin (z) ≈ z

Plugging this into the equation for N , we get

 
(2S + 1)V z 1 N 3
N= →z= λ ≪1
3
λth 2S + 1 V th

We would then recover the ideal gas law

3
U = N kB T, −ΦG = pV = N kB T
2
May 14th 2024 Quantum Statistics Duong Vu

Fermi gas
Consider a Fermi gas at T = 0. At this temperature, fermions will first occupy
the lowest energy states (2S +1 states per energy level) until they reach the highest
energy level, usually denoted EF , known as the Fermi energy. It is defined as

EF = µ(T = 0)

This makes sense for two reasons. First, we have that

∂U U (N ) − U (N − 1)
µ(T = 0) = = = EF
∂N N − (N − 1)
and second, looking at our distribution function at T = 0,
1
f (E) = = Θ(µ − E) = Θ(EF − E)
eβ(E−µ) + 1
which says that all states with energy below EF will be occupied, and all states
with energy above EF will not.
The number of states occupied (and therefore the number of particles) is
Z kF
N= g(k) d3 k
0

where kF is the Fermi wave vector

ℏ2 kF2
EF =
2m
We then get
4 3
1 3 πkF (2S + 1)V kF3
N= × × (2S + 1) =
8 (π/L)3 2π 2 3
so that writing n = N/V , we have
1/3
6π 2 n

kF =
2S + 1

and hence 2/3

ℏ2 6π 2 n

EF =
2m 2S + 1
May 14th 2024 Quantum Statistics Duong Vu

We can calculate the total energy

Z EF
U= Eg(E) dE
0

but in this case, it’s slightly easier to calculate the average energy (over all
fermions)
R EF
Eg(E) dE
⟨E⟩ = R0 EF
0 g(E) dE
To evaluate this, all we need to know is that g(E) is proportional to E 1/2 . We
then get
3
⟨E⟩ = EF
5
and therefore
3
U = N E = N EF
5
At T = 0, the first law becomes

dU = T dS − pdV + µdN
= −pdV + µdN

Using this, the pressure can be written as

 
∂U 3 2 2U
p=− = N EF × =
∂V T,V 5 3V 3V

which is not zero, even at T = 0! It is called the degeneracy pressure.

The Sommerfeld expansion
It is useful to evaluate the integral of the following form for non-zero temperature
Z ∞
I= ϕ(E)f (E) dE
0

To do this, we first define the following function

Z E
ψ(E) = ϕ(u) du
0

so that we have
dψ
= ϕ(E), ψ(0) = 0
dE
May 14th 2024 Quantum Statistics Duong Vu

Next, we integrate by part our original integral

Z ∞ ∞ Z ∞
dψ df
I= f (E) dE = ψ(E)f (E) − ψ(E) dE
0 dE 0 0 dE
Z ∞
df
=− ψ(E) dE
0 dE

Then, we evaluate df /dE

df βeβ(E−µ)
= − β(E−µ)
dE (e + 1)2
Now, note that df /dE is the slope of the distribution function. It will be the
greatest around the point E = µ. This motivates us to change the integration
variable to x = β(E − µ). We then get
Z ∞
ex
I= ψ(x) x dx
−βµ (e + 1)2
There are two more important steps. The first one is to note that for small T
(large β), we can replace the lower bound of the integral by −∞. The second one
is to write ψ(x) as a Taylor expansion around x = 0
∞
xs ds ψ
X  
ψ(x) =
s=0
s! dxs x=0

Our integral becomes

∞  Z ∞
1 ds ψ xs ex
X 
I= s x + 1)2
dx
s=0
s! dx x=0 −∞ (e

Our next step is to evaluate this integral. Note that we can write
ex 1
=
(ex + 1)2 (ex/2 + e−x/2 )2
so this is an even function. Therefore, the integral would vanish for odd values of
s. For even s, this becomes
Z ∞ Z ∞
xs ex xs ex
x 2
dx = 2 dx
−∞ (e + 1) 0 (ex + 1)2
May 14th 2024 Quantum Statistics Duong Vu

To evaluate this, we try to write ex /(ex + 1)2 as an infinite sum. We know the
infinite sum of 1/(ex + 1)
∞
1 e−x −x
X
x
= −x
=e (−1)n e−nx
e + 1 1 − (−1) · e n=0

where we factored out e−x to make sure that the common ratio of the geometric
sum is smaller than or equal to 1. The infinite sum of ex /(ex + 1)2 is just the
negative of the derivative of this sum
∞ ∞
ex d X n −(n+1)x
X
x 2
=− (−1) e = (−1)n (n + 1)e−(n+1)x
(e + 1) dx n=0 n=0

Substituting this back into our integral, we get

Z ∞ Z ∞ X ∞
xs ex
x 2
dx = 2 x s
(−1)n (n + 1)e−(n+1)x dx
−∞ (e + 1) 0 n=0
X∞ Z ∞
=2 (−1)n (n + 1) xs e−(n+1)x dx
n=0 0
∞ Z ∞
X 1
=2 (−1) n−1
n· us e−u du
n=1
ns+1 0
∞
X (−1)n+1
= 2(s!)
n=1
ns
∞ ∞
" ! !#
X 1 X 1
= 2(s!) −2
n=1
ns n=1
(2n)s
1−s
= 2(s!)(1 − 2 )ζ(s)
As a result, our original integral evaluates to
∞  s 
X dψ
I= 2 s
(1 − 21−s )ζ(s)
s=0,2,...
dx x=0
π 2 d2 ψ 7π 4 d4 ψ
   
= ψ x=0 + + + ···
6 dx2 x=0 360 dx4 x=0
Z µ
π2 4
   3 
dϕ 7π dϕ
= ϕ(E) dE + (kB T )2 + (kB T )4 + ···
0 6 dE E=µ 360 dE 3 E=µ
May 14th 2024 Quantum Statistics Duong Vu

We can apply this to find the chemical potential of an electron gas (S = 1/2) at
non-zero temperature
Z ∞ "  3/2 # Z ∞
(2S + 1)V 2m
N= g(E)f (E) dE = 2 2
E 1/2 f (E) dE
0 (2π) ℏ 0
" 3/2 # " #
 2
  2 4
 4
2 (2S + 1)V 2m 3/2 π kB T 7π kB T
= µ 1 + + + ···
3 (2π)2 ℏ2 8 µ 640 µ

Since the change in µ is small, we can put µ3/2 to the other side and replace the µ
in the denominators with µ(0) to find the first order term. Finding higher-order
terms requires more work; according to Wikipedia,
" 2 4 #
2 4
 
π kB T π kB T
µ(T ) = µ(0) 1 − − + ···
12 µ(0) 80 µ(0)

Next, we can find the internal energy and, from there, the heat capacity.
 3/2 Z ∞
(2S + 1)V 2m
U= E 3/2 f (E) dE
(2π)2 ℏ2
0
3/2 " 2 #
2
 
2 (2S + 1)V 2m 5π kB T
= 2 2
µ(T )5/2 1 + + ···
5 (2π) ℏ 8 µ
" 2 #
π 2 kB T

3
= N µ(T ) 1 + + ···
5 2 µ(0)
" 2 #
2

3 5π kB T
= N µ(0) 1 + + ···
5 12 µ(0)

and hence
π 2 kB T
 
3
+ O T3


CV = N k B
2 3 µ(0)
May 14th 2024 Quantum Statistics Duong Vu

Bose gas
For a Bose gas, recall that we have
(2S + 1)V
N= Li3/2 (z)
λ3th

and
3 Li5/2 (z)
U = N kB T
2 Li3/2 (z)
We also have that
1
f (E) =
eβ(E−µ) − 1
Therefore, for a Bose gas, the chemical potential must always be smaller than
or equal to the smallest energy level; otherwise, the occupation number will be
negative at low energy. As a result, the chemical potential must be ≤ 0 for a
gapless dispersion relation like E = ℏ2 k 2 /2m, and hence the fugacity z = eβµ
must lie between 0 and 1. We can rearrange the equation for N to get

nλ3th
= Li3/2 (z)
2S + 1
Here, we see a problem. The right-hand side is maximum when z = 1, at which it
is equal to ζ(3/2). The left-hand side, however, can be made arbitrarily large by
adding more particles (and therefore increasing n) or decreasing the temperature.
Therefore, for some range of values of N and T , we can solve for z. But, if we
have
nλ3th
 
3
>ζ = 2.612
2S + 1 2
there is no solution. What has happened?
May 14th 2024 Quantum Statistics Duong Vu

Bose - Einstein condensation

The problem is that our approximation of going from a sum to an integral breaks
down at this point. At the critical temperature given by
2/3
2πℏ2

n
kB Tc =
m 2.612(2S + 1)
the ground state has become macroscopically occupied; thus, our approximation
becomes invalid. To treat this, we separate N into two terms
N = N0 + N1
where N0 is the number of particles in the ground state, and N1 is our original
integral representing the expected number of bosons in all the other states. Above
the critical temperature, we have
(2S + 1)V
N = N1 = Li3/2 (z)
λ3th
Below the critical temperature, N1 saturates to the calculated value
 
(2S + 1)V 3
N1 = ζ
λ3th 2
and therefore, the fraction of particles in the ground state is
   3/2
N0 N − N1 (2S + 1)V 3 T
= =1− ζ = 1 −
N N N λ3th 2 Tc
We can plot this as follows

1
N0 /N

0
0 1
T /Tc
May 14th 2024 Quantum Statistics Duong Vu

With this, we can calculate all thermodynamic quantities for the Bose gas gener-
ally. Starting with the internal energy, for T ≤ Tc ,

ζ 52


3
U = N1 kB T 3

2 ζ 2

 3/2
3 ζ 52 T
= N kB T 3

2 ζ 2 Tc
 5/2
T
= 0.77N kB Tc
Tc

and of course, for T > Tc we have

3 Li5/2 (z)
U = N kB T
2 Li3/2 (z)

Taking the derivative of this gives us the heat capacity. For T ≤ Tc , we have
5

 3/2
ζ
15 2 T
CV = 3

N kB
4 ζ 2 Tc
 3/2
T
= 1.93N kB
Tc

and for T > Tc , we get (after some tedious algebra, which I will do later)
 
3 5 Li5/2 (z) 3 Li3/2 (z)
CV = N k B −
2 2 Li3/2 (z) 2 Li1/2 (z)