Peeter Joot [email protected]

Huang 9.3

Exercise 1.1

Relativisitic Fermi gas ([1], pr 9.3)

Consider a relativisitic gas of N particles of spin 1/2 obeying Fermi statistics, enclosed in volume V, at absolute zero. The energy-momentum relation is q e = (pc)2 + e02 , (1.1) where e0 = mc2 , and m is the rest mass. a. Find the Fermi energy at density n. b. With the pressure P defined as the average force per unit area exerted on a perfectly-reflecting wall of the container. Set up expressions for this in the form of an integral. c. Define the internal energy U as the average e − e0 . Set up expressions for this in the form of an integral. d. Show that PV = 2U /3 at low densities, and PV = U /3 at high densities. State the criteria for low and high densities. e. There may exist a gas of neutrinos (and/or antineutrinos) in the cosmos. (Neutrinos are massless Fermions of spin 1/2.) Calculate the Fermi energy (in eV) of such a gas, assuming a density of one particle per cm3 . f. Attempt exact evaluation of the various integrals. Answer for Exercise 1.1 Part a.

We’ve found [3] that the density of states associated with a 3D relativisitic system is q 4πV e e2 − e02 , D (e) = (ch)3

(1.2)

For a given density n, we can find the Fermi energy in the same way as we did for the nonrelativisitic energies, with the exception that we have to integrate from a lowest energy of e0 instead of 0 (the energy at p = 0). That is

1

N V   Z eF q 1 4π = 2 +1 e2 − e02 dee 2 (ch)3 e0
eF 8π 1 2 2 3/2 = x − e 0 (ch)3 3 e0
8π 3 / 2 = e2 − e02 . 3(ch)3 F

n=

(1.3)

Solving for eF /e0 we have v u  2/3 eF u 3(ch)3 n t = + 1. e0 8πe03

(1.4)

We’ll see the constant factor above a number of times below and designate it n0 =

8π  e0 3 , 3 ch

(1.5)

so that the Fermi energy is s eF = e0



n n0

 2/3 + 1.

(1.6)

Part b. For the pressure calculation, let’s suppose that we have a configuration with a plane in the x, y orientation as in fig. 1.1.

Figure 1.1: Pressure against x, y oriented plane

2

It’s argued in [4] §6.4 that the pressure for such a configuration is P=n

Z

pz uz f (u)d3 u,

(1.7)

where n is the number density and f (u) is a normalized distribution function for the velocities. The velocity and momentum components are related by the Hamiltonian equations. From the Hamiltonian eq. (1.1) we find 1 (for the x-component which is representative) ∂e ∂p x q ∂ = (pc)2 + e02 ∂p x p x c2 =q . (pc)2 + e02

ux =

(1.8)

For α ∈ {1, 2, 3} we can summarize these velocity-momentum relationships as uα cpα = . c e

(1.9)

Should we attempt to calculate the pressure with this parameterization of the velocity space we end up with convergence problems, and can’t express the results in terms of f ν+ (z). Let’s try instead with a distribution over momentum space P=n

Z

(cpz )2 f (cp)d3 (cp). e

(1.10)

Here the momenta have been scaled to have units of energy since we want to express this integral in terms of energy in the end. Our normalized distribution function is f (cp) ∝ R

1 z−1 e βe +1 , 1 d3 (cp) z−1 e βe +1

(1.11)

but before evaluating anything, we first want to change our integration variable from momentum to energy. In spherical coordinates our volume element takes the form d3 (cp) = 2π(cp)2 d(cp) sin θdθ d(cp) = 2π(cp)2 de sin θdθ. de

(1.12)

c2 p2 = e2 − e02 ,

(1.13)

Implicit derivatives of

1

Observe that by squaring and summing one can show that this is equivalent to the standard relativisitic momentum p x = √ mvx2 2 . 1− u / c

3

gives us d(cp) e e = =q . de cp e2 − e02

(1.14)

Our momentum volume element becomes e

d3 (cp) = 2π(cp)2 q

de sin θdθ e2 − e02
e = 2π e2 − e02 q de sin θdθ e2 − e02 q = 2πe e2 − e02 de sin θdθ.

For our distribution function, we can now write q e e2 − e02 de 2π sin θdθ f (cp)d3 (cp) = C −1 βe , z e +1 4πe03 R where C is determined by the requirement f (cp)d3 (cp) = 1 p Z ∞ (y + 1) (y + 1)2 − 1dy −1 C = . z−1 e βe0 (y+1) + 1 0

(1.15)

(1.16)

(1.17)

The z component of our momentum can be written in spherical coordinates as
(cpz )2 = (cp)2 cos2 θ = e2 − e02 cos2 θ,

(1.18)

Noting that Z π 0

cos2 θ sin θdθ = −

Z π 0

cos2 θd(cos θ) =

2 , 3

(1.19)

all the bits come together as Z
3/2 Cn ∞ 2 1 e − e02 de 3 − 1 z e βe + 1 3e0 e0 Z
3/2 1 ne0 ∞ 2 x −1 = dx. βe − 1 3 1 z e 0x + 1

P=

(1.20)

Letting y = x − 1, this is Cne0 P= 3


3/2 Z ∞ (y + 1)2 − 1 0

z−1 e βe0 (y+1) + 1

dy.

(1.21)

We could conceivable expand the numerators of each of these integrals in power series, which could then be evaluated as a sum of f ν+ (ze− βe0 ) terms.

4

Note that above the Fermi energy n also has an integral representation    Z ∞ 1 1 2 +1 deD (e) −1 βe 2 z e +1 e q0 Z ∞ e e2 − e02 8π = de −1 βe (ch)3 e0 z e +1 p Z ∞ 3 8πe0 (y + 1) (y + 1)2 − 1 = dy , (ch)3 0 z−1 e βe0 (y+1) + 1

n=

(1.22)

or 3n0 (1.23) . C Observe that we can use this result to remove the dependence of pressure on this constant C n=

Z ∞

P = n 0 e0 Part c.

0


3/2 (y + 1)2 − 1 dy −1 βe (y+1) . z e 0 +1

(1.24)

Now for the average energy difference from the rest energy e0 U = h e − e0 i Z ∞

deD (e) f (e)(e − e0 ) p Z 8πV ∞ e(e − e0 ) e2 − e0 de = (ch)3 e0 z−1 e βe + 1 p Z 8πVe04 ∞ y(y − 1) (y + 1)2 − 1 = . dy (ch)3 0 z−1 e βe + 1

=

e0

So the average energy density difference from the rest energy, relative to the rest energy, is p Z ∞ y(y + 1) (y + 1)2 − 1 h e − e0 i = 3n0 dy . Ve0 z−1 e βe0 (y+1) + 1 0 Part d.

(1.25)

(1.26)

From eq. (1.24) and eq. (1.26) we have 1 Ve0 =3 n0 h e − e0 i =

Z e0 ∞

P

0

p Z ∞ y(y + 1) (y + 1)2 − 1dy

z−1 e βe0 (y+1) + 1
3/2 (y + 1)2 − 1 dy, z−1 e βe0 (y+1) + 1 0

(1.27)

or

R ∞ ((y+1)2 −1)3/2 U 0 z−1 eβe0 (y+1) +1 dy √ PV = R . 3 ∞ y(y+1) (y+1)2 −1dy z−1 e βe0 (y+1) +1

0

5

(1.28)

This ratio of integrals is supposed to resolve to 1 and 2 in the low and high density limits. To consider this let’s perform one final non-dimensionalization, writing

x = βe0 y 1 kB T θ= = βe0 e0

(1.29)

µ = µ − e0 z = e βµ . The density, pressure, and energy take the form Z ∞

p (θx + 1) (θx + 1)2 − 1 dx z −1 e x + 1 0
3/2 Z ∞ (θx + 1)2 − 1 P =θ dx n 0 e0 z −1 e x + 1 0 p Z ∞ x(θx + 1) (θx + 1)2 − 1 h e − e0 i 2 = 3θ dx . Ve0 n0 z −1 e x + 1 0 n = 3θ n0

(1.30a) (1.30b) (1.30c)

We can rewrite the square roots in the number density and energy density expressions by expanding out the completion of the square (1 + θx)

p

√ √ (1 + θx)2 − 1 = (1 + θx) 1 + θx + 1 1 + θx − 1 r √ θx 1/2 , = 2θx (1 + θx) 1 + 2

(1.31)

¯ − x = 0, we have Expanding the distribution about ze ∞
1 ze− x −x s −x s ( − 1) = ze = ze , ∑ z−1 e x + 1 1 + ze− x s=0

(1.32)

allowing us to write, in the low density limit with respect to z ∞ √ n = 3 2θ 3/2 ∑(−1)s zs+1 n0 s=0 ∞ P = θ ∑(−1)s zs+1 n 0 e0 s=0

Z ∞ 0

Z ∞ 0

∞ √ h e − e0 i = 3 2θ 5/2 ∑(−1)s zs+1 Ve0 n0 s=0

r dxx1/2 (1 + θx)

dx (θx + 1)2 − 1 Z ∞ 0

6

1+
3/2

θx − x(1+s) e 2

e− x(1+s)

r dxx

3/2

(1 + θx)

1+

θx − x(1+s) e . 2

(1.33a) (1.33b)

(1.33c)

Low density result An exact integration of the various integrals above is possible in terms of special functions. However, that attempt (included below) introduced an erroneous extra factor of θ. Given that this end result was obtained by tossing all but the lowest order terms in θ and z, let’s try that right from the get go. For the pressure we have an integrand containing a factor 2


3/2

 3/2

√

(θx + 1) − 1 θx + 2

3/2

3/2

= (θx + 1 − 1)

= 2 2θ 3/2 x3/2



θx 1+ 2

(θx + 1 + 1)  3/2

=θ

3/2 3/2 3/2

x

2

 1 (1.34)

√

≈ 2 2θ 3/2 x3/2

Our pressure, to lowest order in θ and z is then

√ 5/2 Z ∞ 3/2 − x P x e dx = 2 2θ z e0 n 0 √ 5/2 0 = 2 2θ zΓ(5/2).

(1.35)

Our energy density to lowest order in θ and z from eq. (1.33c) is Z ∞ √ U = 3 2θ 5/2 z dxx3/2 e− x Ve0 n0 0 √ = 3 2θ 5/2 zΓ(5/2).

(1.36)

Comparing these, we have

e0 n 0

√

1 2θ 5/2 zΓ(5/2)

=3

V 2 = , U P

(1.37)

or in this low density limit 2 PV = U. 3 High density limit

(1.38)

For the high density limit write z = ey , so that the distribution takes the form

1 1 = . z −1 e x + 1 e x − y + 1 This can be approximated by a step function, so that f (z) =

P ≈ n 0 e0 U ≈3 Ve0 n0

Z y 0

Z ∞ 0

θdx (θx + 1)2 − 1

θdxθx(θx + 1)

7

p


3/2

(θx + 1)2 − 1

(1.39)

(1.40a) (1.40b)

With a change of variables u = θx + 1, we have P ≈ n 0 e0

Z θy+1x 1

du u2 − 1


3/2

  p p 1 = (2θy(θy + 2) − 3) θy(θy + 2)(θy + 1) + 3 ln θy + θy(θy + 2) + 1 8 1 ≈ (θ ln z)4 4 Z θy+1x p U ≈3 (u2 − u) u2 − 1 Ve0 n0 1    p 3 p = θy(θy + 2)(θy(2θy(3θy + 5) − 1) + 3) − 3 ln θy + θy(θy + 2) + 1 24 3 ≈ (θ ln z)4 4 Comparing both we have

or

(1.41a)

(1.41b)

1 3V 4 = = , e0 n0 (θ ln z) P U

(1.42)

1 PV = U. 3

(1.43)

Part e. eF |n =1/(0.01)3 = 6.12402 × 10−35 J × 6.24150934 × 1018

eV = 3.82231 × 10−16 eV J

(1.44)

Wow. That’s pretty low! Part f. Pressure integral Z ∞ 0

Of these the pressure integral is yields directly to Mathematica


3/2 − x(1+s) dx (θx + 1)2 − 1 e   3θe(s+1)/θ s+1 = K2 (s + 1)2 θ p π 3/2 p π 5/2 p p p 45 2 θ 315 π2 θ 7/2 945 π2 θ 9/2 31185 π2 θ 11/2 3 2θ = + + − + +··· (s + 1)5/2 8(s + 1)7/2 128(s + 1)9/2 1024(s + 1)11/2 32768(s + 1)13/2

where K2 (z) is a modified Bessel function [5] of the second kind as plotted in fig. 1.2.

8

(1.45)

Figure 1.2: Modified Bessel function of the second kind Plugging this into the series for the pressure, we have P =3 n 0 e0



kB T e0

2



∞

∑(−1)s

s=0

zee0 /kB T

s+1

(s + 1)2

K2 ((s + 1)e0 /kB T ) .

(1.46)

s+1  θ2 1/ θ Plotting the summands 3(−1)s (s+1) ze K2 ((s + 1)/θ ) for z = 1 in fig. 1.4 shows that this mix 2 of exponential Bessel and quadratic terms decreases with s. Plotting this sum in fig. 1.3 numerically to 10 terms, shows that we have a function that appears roughly polynomial in z and θ.

1 ´ 106

500 000 4 0 0 2

z

2 Θ 4 0

Figure 1.3: Pressure to ten terms in z and θ

9

Figure 1.4: Pressure summands For small z it can be seen graphically that there is very little contribution from anything but the s = 0 term of this sum. An expansion in series for a few terms in z and θ gives us   √ 3z P 3z2 z3 3z4 3z5 = πθ 5/2 √ − + √ − √ + √ e0 n 0 8 3 6 32 2 25 10 2   √ 7/2 45z 45z2 5z3 45z4 9z5 (1.47) √ √ + √ − + √ − + πθ 128 24 6 1024 2 200 10 8 2   √ 9/2 315z 315z4 315z2 35z3 63z5 √ − √ − √ + √ + πθ + . 4096 1152 6 65536 2 16000 10 128 2 This allows a kB T  mc2 and z  1 approximation of the pressure P 3√ 2πzθ 5/2 . = e0 n 0 2

(1.48)

Number density integral For the number density, it appears that we can evaluate the integral using integration from parts applied to eq. (1.30) Z ∞

p 3(θx + 1) (θx + 1)2 − 1 dx z −1 e x + 1 0   Z ∞
3/2 d 1 2 =θ dx (θx + 1) − 1 − 1 dx z ex + 1 0 ∞ Z ∞
3/2
3/2 1 2 2 = θ (θx + 1) − 1 − θ dx (θx + 1) − 1 z −1 e x + 1 0

n =θ n0

0

=θ

Z ∞ 0

dx (θx + 1)2 − 1


3/2

ze− x (1 + ze− x )2

10

.

(1.49)

− z −1 e x z −1 e x + 1


2

Expanding in series, gives us  Z ∞ 
3/2 − x(s+1) n −2 s+1 ∞ =θ∑ z dx (θx + 1)2 − 1 e n0 s 0 s=0     ze1/θ s+1  ∞  s+1 −2 K . = 3θ 2 ∑ 2 (s + 1)2 θ s s=0

(1.50)

Here the binomial coefficient has the meaning given in the definitions of [2], where for negative integral values of b we have     −b + s b s −b ≡ (−1) . (1.51) s −b + s −b Expanding in series to a couple of orders in θ and z we have

√ √   √  √  n 2π 1/2  √ 5 2π 3/2  √ θ 2 3z − 9/ 2 z + 18 z + θ 4 3z − 27/ 2 z + 108 z + · · · = n0 36 576 To first order in θ and z this is

(1.52)

1√ n = 2πzθ 1/2 , n0 2

(1.53)

PV = 3N(kB T)2 /e0 .

(1.54)

which allows a relation to pressure

It’s kind of odd seeming that this is quadratic in temperature. Is there an error? Energy integral

Starting from eq. (1.30c) and integrating by parts we have p Z ∞ x(θx + 1) (θx + 1)2 − 1 h e − e0 i 2 = 3θ dx Ve0 n0 z −1 e x + 1 0 Z ∞ 
3/2 d  x = −θ 2 dx (θx + 1)2 − 1 dx z−1 e x + 1 0 ! Z ∞
3/2 xz−1 e x 1 2 2 − = −θ dx (θx + 1) − 1
2 z −1 e x + 1 0 z −1 e x + 1 = θ2 = θ2

Z ∞ 0

Z ∞

dx (θx + 1)2 − 1


3/2 (x − 1)z−1 e x − 1
2 z −1 e x + 1
3/2 (x − 1)ze− x − z2 e−2x

dx (θx + 1)2 − 1 0 (1 + ze− x )2   Z ∞ ∞
3/2
−2 = θ2 ∑ dx (θx + 1)2 − 1 (x − 1)ze− x − z2 e−2x (ze− x )s s 0 s=0  Z  ∞ 
3/2  −2 s+1 ∞ = θ2 ∑ z dx (θx + 1)2 − 1 (x − 1)e− x(s+1) − ze− x(s+2) . s 0 s=0

11

(1.55)

The integral with the factor of x doesn’t have a nice closed form as before (if you consider the K2 a nice closed form), but instead evaluates to a confluent hypergeometric function [6]. That integral is  √ 3  3 2(s+1) Z ∞ 15 , − 4, πθ U −
3/2 − x(1+s) 2 θ x (θx + 1)2 − 1 e dx = , (1.56) 5 8(s + 1) 0 and looks like fig. 1.5. Series expansion shows that this hypergeometricU function has a θ 3/2 singularity at the origin

Figure 1.5: Plot of HypergeometricU, and with θ 5 scaling 

3 2(s + 1) U − , −4, 2 θ



√ √ √ √ √ 2 2 s + 1s + 2 2 s + 1 21 s + 1 = + √ √ +··· θ 3/2 2 2 θ

(1.57)

so our multiplication by θ 5 brings us to zero as seen in the plot. Evaluating the complete integral yields the unholy mess ∞



  √ 105 πθ 3 U − 21 , −4, 2(s+1) θ

h e − e0 i = ∑ θ 2 (−1)s (s + 1)zs+1  Ve0 n0 s=0  √ 2  1 3 πθ zU − 2 , −2, 2(s+2) θ

16(s + 1)5

−

  √ 3 πθ 2 U − 12 , −2, 2(s+1) θ


s+1 (θ − 2)(−3θ + 2s + 2)e θ K2 s+1 θ + 2(s + 2)3 θ(s + 1)2

s+2 s+2 s+1 2(θ − 2)e θ K1 s+1 z(−3θ + 2s + 4)e θ K2 s+2 2ze θ K1 θ θ − + − θ(s + 1) (s + 2)2 s+2

2(s + 1)3

− (1.58)

s+2 θ


 ,

to first order in z and θ this is

h e − e0 i 9 √ = 2πzθ 7/2 . Ve0 n0 4

12

(1.59)

Comparing pressure and energy we have for low densities (where z ≈ 0)

e0 n 0

√

1 2πzθ 5/2

=

or

31 9 V = θ , 2P 4 U

(1.60)

2 θPV = U. (1.61) 3 It appears that I’ve picked up an extra factor of θ somewhere, but at least I’ve got the 2/3 low density expression. Given that I’ve Taylor expanded everything anyways around z and θ this could likely have been done right from the get go, instead of dragging along the messy geometric integrals. Reworking this part of this problem like that was done above.

13

Bibliography [1] Kerson Huang. Introduction to statistical physics. CRC Press, 2001. 1.1 [2] Peeter Joot. Basic statistical mechanics., chapter Non integral binomial coefficient. . URL http: //sites.google.com/site/peeterjoot2/math2013/phy452.pdf. 1 [3] Peeter Joot. Basic statistical mechanics., chapter Relativisitic density of states. . URL http://sites. google.com/site/peeterjoot2/math2013/phy452.pdf. 1 [4] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996. 1 [5] Wolfram. BesselK, . URL http://reference.wolfram.com/mathematica/ref/BesselK.html. [Online; accessed 11-April-2013]. 1 [6] Wolfram. HyperGeometricU, . URL http://reference.wolfram.com/mathematica/ref/ HypergeometricU.html. [Online; accessed 17-April-2013]. 1

14