Statistical Mechanics of Ideal Fermi Systems

from Statistical Physics using Mathematica
© James J. Kelly, 1996-2002

Introduction

According to the Fermi-Dirac distribution function, which is derived within the notebook occupy.nb, the mean occupancy for a single-particle orbital with energy ε is

\(n\&_[ε]\  = \ 
      1\/\(Exp[\(ε - μ\)\/\(\(k\_B\) T\)]\  + \ 1\)\)

where the chemical potential μ is a function of density and temperature.  At T = 0, the argument of the exponential is -∞ when ε < μ or   +∞ when ε > μ; hence, the occupancy is unity for all states with ε below μ and is zero for all states with ε above μ.  Therefore, at absolute zero a Fermi gas is described as completely degenerate and is characterized by a frozen distribution in which all orbitals with ε < (Subscript[ε, F]) are occupied and all orbitals with ε > (Subscript[ε, F]) are vacant, where the Fermi energy and Fermi temperature are defined in terms of the chemical potential as

\(T = \(0\ \ \ \ ⟹\ \ \ \ μ = \(ε\_F\  = \ \(k\_B\) T\_F\)\)\)

and are functions of the density of the system.  Thus, it is useful to express the thermodynamic properties of Fermi systems at finite temperature, which are no longer completely degenerate, in terms of reduced energy and temperature variables scaled to (Subscript[ε, F]) or (Subscript[T, F]), respectively.

The dependence of the chemical potential upon density and temperature is determined by integrating the occupation numbers to obtain the total particle number

\(N\  = \ ∫\_0\%∞ n\&_[ε]\ 𝒟[ε] dε\)

where for large systems the density of states, 𝒟[ε], is approximated by a continuous function of energy.  For an ordinary nonrelativistic gas in three dimensions, the single-particle energy is simply \(TraditionalForm\`ε\  = \ \(\(ℏ\^2\) k\^2\)\/\(2\ m\)\) and the density of momentum states is 𝒟[(Overscript[k, ⇀])] (d^3) k =  
        g V  (((d^3) k)/(((2  π))^3)) where g = 2  s + 1 is the intrinsic degeneracy for a particle with spin s and V is the volume of the system.  Generalizations to one- or two-dimensional systems or to relativistic kinematics are developed in the problems at the end of the notebook.  The chemical potential enters the above integral through  \(TraditionalForm\`n\&_[ε]\) in a rather nontrivial fashion that requires numerical evaluation.  Similarly, other thermodynamic properties can be expressed in terms of integrals of this type.  Therefore, it is useful to define the fermi function of order ν as

\(fermi[ν, 
        z]\  = \ \(1\/Γ[ν]\) \(∫\_0\%∞\( x\^\(ν - 1\)\/\(\(z\^\(-1\)\) Exp[x] + 
                  1\)\) dx\)\)

where z  =  Exp[(μ/((Subscript[k, B]) T))] is known as the fugacity and where x = ε/(Subscript[ε, F]) is the energy in units of the Fermi energy.  The gamma function is included to ensure a convenient normalization, namely fermi[ν, z] -> z for small z.  Using these fermi functions, we will show below that some of the common thermodynamic functions for an ideal nonrelativistic Fermi gas in three dimensions become

(Subscript[n, Q]) =  ((N (λ^3))/(g V)) = 
        fermi[(3/2), z]

𝒢 = -p V  =  -N (Subscript[k, B]) 
          T fermi[(5/2), z]/fermi[(3/2), z]

U  =  (3/2) N (Subscript[k, B]) T 
         fermi[(5/2), z]/fermi[(3/2), z]

where (Subscript[n, Q]) is the quantum concentration and λ=(Sqrt[((2 π (ℏ^2))/(m (Subscript[k, B]) T))]) is the thermal wavelength.  The properties of fermi functions are investigated below in considerable detail.

Glossary

𝒩	= number of particles
	V	= volume
	T	= temperature
	(Subscript[k, B])	= Boltzmann constant
	g	= energy-level degeneracy
	ε	= single-particle energy
	k	= single-particle wave number
	𝒟[ε]	= density of states wrt energy
	𝒟[k]	= density of states wrt momentum
	λ	= thermal wavelength
	(Subscript[n, Q])	= quantum concentration
	μ	= chemical potential
	z	= fugacity
	(Subscript[𝒵, 1]^FD)	= single-particle grand-partition function
	\(TraditionalForm\`n\&_\)	= mean occupation number
	(Subscript[ε, F])	= Fermi energy
	(Subscript[k, F])	= Fermi momentum
	(Subscript[T, F])	= Fermi temperature ((Subscript[ε, F])/(Subscript[k, B]))
	τ	= reduced temperature (T/(Subscript[T, F]))
	𝒢	= grand potential
	U	= internal energy
	p	= pressure
	S	= entropy
	F	= free energy
	H	= enthalpy
	(Subscript[C, V])	= isochoric heat capacity
	(Subscript[C, p])	= isochoric heat capacity
	γ  =  ((Subscript[C, p])/(Subscript[C, V]))	= ratio of principal heat capacities
	fermi[n, z]	= Fermi-Dirac function of order n
	α  =  
        (1/V) (Subscript[((∂(V/∂)T)), p]) 	= isobaric expansivity
	\(TraditionalForm\`κ\_T = \ \ \(-\(1\/V\)\)\ \((∂V\/∂p)\)\_T\)	= isothermal compressibility
	\(TraditionalForm\`κ\_S = \ \ \(-\(1\/V\)\)\ \((∂V\/∂p)\)\_S\)	= adiabatic compressibility

Initialization

Defaults, packages, and symbols

ClearAll["Global`*"];
Off[General::spell, General::spell1];

$DefaultFont = {"Times", 12};
$TextStyle = {FontFamily -> "Times", FontSize -> 12, 
          FontSlant -> "Italic"};

Needs["Utilities`Notation`"];
Needs["Miscellaneous`PhysicalConstants`"];
Needs["Miscellaneous`Units`"];
Needs["Graphics`Master`"]

Symbolize[(Subscript[𝒵, 1]^FD)];Symbolize[(Subscript[k, B])];Symbolize[(Subscript[k, F])];Symbolize[(Subscript[ε, F])];Symbolize[(Subscript[κ, x])_]

Symbolize[(Subscript[λ, x])_];Symbolize[(Subscript[T, x])_];Symbolize[(Subscript[C, x])_];Symbolize[(Subscript[v, x])_];Symbolize[(Subscript[V, x])_]

Symbolize[(Subscript[n, x])_];Symbolize[(Subscript[N, x])_];Symbolize[\(n\&_\)];Symbolize[\(\(n\_x_\)\&_\)];Symbolize[(Subscript[k, max])]

Symbolize[\(B\&~\)];Symbolize[\(p\&~\)];Symbolize[\(U\&~\)];Symbolize[\(S\&~\)];Symbolize[\(F\&~\)];Symbolize[\(G\&~\)];Symbolize[\(C\&~\_x_\)];Symbolize[\(γ\&~\)]

SetAttributes[{(Subscript[k, B]),  ℏ,  g,  m},  Constant];

FundamentalConstants  =  {(Subscript[k, B]) -> 
            (BoltzmannConstant/((Joule/Kelvin))), ℏ -> 
            (PlanckConstantReduced/(Joule Second))};

Memory management

This notebook consumes enough memory that it becomes advantageous to run the memory conservation utility.  However, we are not interested in seeing the messages reporting the actions of that process.

Needs["Utilities`MemoryConserve`"];
Off[MemoryConserve::start,  MemoryConserve::end];

Error messages

Several annoying error messages can be suppressed, if desired, but should be enabled when developing or debugging code.

Off[General::"ovfl"];
Off[Solve::"ifun"];
Off[NIntegrate::"inum"];
Off[NIntegrate::"ncvb"];
Off[NIntegrate::"slwcon"];
Off[NIntegrate::"precw"];

Rules for changing variables

TtoLambda = {T -> ((2 π (ℏ^2))/(m (Subscript[k, B]) (λ^2))),  β -> (1/((Subscript[k, B]) T))};
lambdaToT = {λ -> (Sqrt[((2 π (ℏ^2))/(m (Subscript[k, B]) T))])};
lambdaToZ = {λ -> (((((g V fermi[(3/2), z])/𝒩)))^(1/3))};

zToTau = z -> fugacity[τ];
muToZ = μ  -> (Subscript[k, B]) T Log[z];
zToT  =  z -> Exp[(μ/((Subscript[k, B]) T))];

Density of states

Evaluate the density of states for nonrelativistic gas in 3 dimensions in terms of both momentum and energy.

𝒟[
          k] = ((4 π g V (k^2))/(((2 π))^3));

(Subscript[ε, k]) = ((((ℏ k))^2)/(2 m));

𝒟[ε] = 𝒟[k]/(∂_k) (Subscript[ε, k]) //. {k -> ((Sqrt[2 m ε])/ℏ)}  // PowerExpand

Occupation numbers

The statistics of occupation numbers have already been explored in the notebook occupy.nb.  Here we repeat the derivation of the mean occupation number, briefly, but then investigate its properties for nearly degenerate Fermi systems in more detail.  The goal in this section is to motivate the definition of the Fermi functions and to explain how expansions with respect to reduced temperature can be developed for various thermodynamic functions.

Mean occupation number

The mean occupation number, as a function of energy and chemical potential, is determined from the grand partition function for a single orbital.

(Subscript[𝒵, 1]^FD) = (∑_(n = 0)^1) 
          Exp[-n β (ε -  μ)]

\(\(n\_FD\)\&_ = 
      k\_B\ T\ ∂\_\(\(\ \)\(μ\)\)\((Log[𝒵\_1\%FD] /. {β -> 
                    1\/\(k\_B\ T\)})\)\  // FullSimplify\)

\(\(n\_FD\)\&_ /. muToZ\  // FullSimplify\)

Evaluate Fermi energy and momentum

The Fermi momentum, (Subscript[k, F]), is defined to be the highest occupied momentum when all lower  states are fully occupied.  The corresponding energy is the Fermi energy, (Subscript[ε, F]).

(Subscript[n, Q]) = ((𝒩 (λ^3))/(g V));

kFrule=Solve[(∫_0^((Subscript[k, F])))𝒟[
                  k] \[DifferentialD]k == 𝒩, 
          (Subscript[k, F])][[3]]

eFrule  =  {(Subscript[ε, F]) -> ((((ℏ (Subscript[k, F])))^2)/(2 m)) } /. kFrule

It is useful to recognize that the density of states, for either momentum or energy, can be expressed in terms of simple functions of the Fermi momentum or energy which display the scaling relations in a transparent manner.   Thus, if we define κ = k/(Subscript[k, F]) and ϵ=ε/(Subscript[ε, F]) we can produce the simpler density of state functions below.

𝒟[κ] = (𝒟[k]  ((Subscript[k, F])^3)  /. 
                kFrule // PowerExpand // Simplify) /. 
        k -> κ

𝒟[ϵ] = (𝒟[ε]  ((Subscript[ε, F])^(3/2))  /. eFrule // PowerExpand // 
            Simplify) /. ε -> ϵ

Thus, we can evaluate the mean energy or momentum, in units of the (Subscript[ε, F]) or (Subscript[k, F]), for a completely degenerate nonrelativistic Fermi gas in three dimensions using simple dimensional arguments.

\(ϵ\&_\  = \ \(∫\_0\%1 ϵ\ 𝒟[ϵ] \[DifferentialD]ϵ\)\/\(∫\_0\%1\ 𝒟[ϵ] \[DifferentialD]ϵ\)\)

\(κ\&_\  = \ \(∫\_0\%1 κ\ 𝒟[κ] \[DifferentialD]κ\)\/\(∫\_0\%1\ 𝒟[κ] \[DifferentialD]κ\)\)

Therefore, we conclude that the internal energy for a completely degenerate Fermi gas is simply \(TraditionalForm\`T = \(0\ \ ⟹\ \ U = \ \(3\/5\) 𝒩\ ε\_F\)\), which will provide an important consistency check on the more detailed numerical evaluations to follow for finite temperature.  Furthermore, recognizing that the pressure is simply two-thirds of the energy density for an ideal nonrelativistic gas, regardless of its permutation symmetry, we conclude that \(TraditionalForm\`T = \(0\ \ ⟹\ \ p\ V = \ \(2\/5\) 𝒩\ ε\_F\)\) for a Fermi gas.  The failure of the energy density and pressure to approach zero as the temperature approaches zero, as expected for a classical gas, represents the most dramatic consequence of the Pauli exclusion principle for fermions.  For dense systems the degeneracy energy and pressure, namely their values at zero temperature, can be extremely large.

Occupation numbers for low temperature

Permutation symmetry has its greatest effect on the thermodynamics of Fermi systems when the reduced temperature is small.  The mean occupation number then approaches a step function of temperature, near unity below the Fermi energy or near zero above the Fermi energy.  Particles occupying states far below the Fermi level cannot change orbitals without absorbing a relatively large amount of energy and hence do not participate in thermodynamic processes.  Therefore, many properties of the system are governed primarily by the states that are in the immediate vicinity of the Fermi level.  We illustrate by comparing the low temperature occupation number distribution with that for a completely degenerate system.  For this purpose it is convenient to express both the energy and temperature in units of the chemical potential, but it also will be important later to remember that the chemical potential is temperature dependent.

step[x_]  :=  If[x < 1, 1, 0]

nfermi[x_, τ_]  :=  
      (1/((E^((x - 1)/τ)) + 1))

\(\(FilledPlot[{step[x], nfermi[\ x, 0.05]}, {x, 0, 1.5}, 
        Frame -> True, 
        FrameLabel -> {"ε/μ", \*""\!\(n\&_\)""}, 
        PlotLabel -> \*""FD Occupancy for T=0.05 (Subscript[T, F])""];\)\)

Particles in the shaded or colored region can be considered active, while the remainder can be considered dormant.  Most of the action occurs near the Fermi surface.  From this plot it is clear that the chemical potential is equal to the single-particle energy for which the mean occupation number is (1/2).

Thermodynamic properties of nearly degenerate systems

In this section we express various thermodynamic functions as expansions with respect to reduced temperature, τ, using the properties of the occupation number distribution.  This procedure provides considerable insight into the behavior of the important class of systems that are nearly degenerate.  Some of these results will be rederived later as low-temperature limits of more general expressions based upon Fermi-Dirac functions.

Expansion with respect to temperature

Generally a thermodynamic function is constructed by integrating some function of the single-particle energy against the density of states and the mean occupation number.  Deviations with respect to the properties of a completely degenerate system arise at finite temperature because the mean occupation number spreads out a little near the Fermi level.  The effect of this change in occupation distributions can be studied using generic integrals of the form

A[τ]  = (∫_0^∞) 
            n[ϵ, τ] (∂_ϵ) 
              f[ϵ] \[DifferentialD]ϵ

where (∂_ϵ) f[ϵ] includes both the density of states and the single-particle function of interest.  Integration by parts is accomplished using the following rule.

\(partsRule\  = \ ∫\_a_\%b_ 
            g_\ ∂\_ϵ 
                f_[ϵ] \[DifferentialD]ϵ\  :> \((g\ f[ϵ] /. ϵ -> 
                b)\)\  - \ \((g\ f[ϵ] /. ϵ -> 
                a)\) - ∫\_a\%b 
                f[ϵ]\ ∂\_ϵ 
                g\ \[DifferentialD]ϵ\)

Hence, assuming that f[ϵ] includes the density of states and vanishes at ϵ=0 while the mean occupation number vanishes at ϵ=∞, we can discard the surface terms and obtain

A[τ] /. partsRule /.{f[0] -> 0, n[∞,τ]->0}

Examining the derivative of the mean occupation number, we recognize that when τ is small it is strongly peaked at the Fermi surface.

dnfermi[x_, τ_] := (∂_x) nfermi[x, τ]

Plot[Evaluate[dnfermi[x, 0.1]],  {x, 0, 2}];

Therefore, if f[ϵ] varies slowly with energy near the Fermi energy, we can use a Taylor series approximation about ϵ = 1.  Furthermore, we can extend the range of integration to (-∞, +∞) and use a change of variables to simplify the integrand, recognizing that dnfermi is symmetric about its peak.

integrand = -f[ϵ] Derivative[1, 0][n][ϵ, τ]/.{f[ϵ] -> 
                  Normal[Series[f[ϵ], {ϵ, 1, 4}]],Derivative[1, 0][n][ϵ, τ]->dnfermi[ϵ, τ]}/.{ϵ -> x τ  +  1}//Expand

expansion  = (τ (∫_(-∞)^∞) 
              integrand \[DifferentialD]x ) // Simplify

Chemical potential

We can apply this expansion to the relationship between density and chemical potential by identifying ((∂f)/(∂ϵ)) with 𝒟[ϵ] and making the necessary substitutions into the expansion formula.

(Subscript[f, N])  = ∫𝒟[ϵ] \[DifferentialD]ϵ

Nexpansion  = ( expansion  /. 
              f[x_] -> (Subscript[f, N])  /. Derivative[n_][f][1] -> 
              D[(Subscript[f, N]), {ϵ, 
                  n}]) /. {ϵ -> μ, τ -> (Subscript[k, B]) T}

Realizing that μ is near (Subscript[ε, F]) for low temperatures, we use the substitution μ->(Subscript[ε, F])+δ to develop an equation for the correction term valid to second order in δ.  After solving this equation, we expand with respect to temperature, being careful to limit that expansion to the same order used for Nexpansion, and select the physically meaningful solution.

stuff  =  ((Nexpansion )/(𝒩 ((Subscript[ε, F])^(3/2))))  /. μ -> (Subscript[ε, F]) + δ // Series[#, {δ, 0, 2}] &  // Normal

Solve[
        stuff == 
          1, δ][[1]]

δsol = (δ /. Solve[stuff == 1, δ]  // 
                Series[#, {T, 0, 4}] &  // Normal) // PowerExpand // 
        Simplify

μExpansion  =  (Subscript[ε, F]) + δsol[[2]]  /. {T -> τ (Subscript[ε, F]) / (Subscript[k, B])} // Collect[#, {(Subscript[ε, F])}] &

Internal energy

To evaluate the temperature dependence of the total energy, in units of  μ, we identify ((∂f)/(∂ϵ)) with  ϵ ρ[ϵ] and make the necessary substitutions into the expansion formula.

(Subscript[f, ϵ])  = ∫ϵ 𝒟[ϵ] \[DifferentialD]ϵ

Ustuff  = ( ((Subscript[ε, F])^(-3/2)) expansion  /. 
              f[x_] -> (Subscript[f, ϵ])  /. Derivative[n_][f][
                1] -> 
              D[(Subscript[f, ϵ]), {ϵ, 
                  n}]) /. {ϵ -> μ, τ -> (Subscript[k, B]) T}

Finally, we use our expansions for the temperature dependence of μ to obtain a simple formula for the temperature dependence of the internal energy.

Uexpansion  =  (Ustuff  /. μ -> μExpansion  /. {T -> τ (Subscript[ε, F]) / (Subscript[k, B])}) // 
            Series[#, {τ, 0, 4}] & // Normal  // 
        Collect[#, {𝒩, (Subscript[ε, F])}] &

Fermi-Dirac functions

Although the expansions developed in the preceding section are very useful at low temperatures, a more complete analysis over a broad range of temperatures requires that we develop more general techniques.  Thus, it is useful to define Fermi functions in terms of integrals of powers of the energy against the occupation number such that

\(fermi[ν, 
        z]\  = \ \(1\/Γ[ν]\) \(∫\_0\%∞\( x\^\(ν - 1\)\/\(\(z\^\(-1\)\) Exp[x] + 
                  1\)\) dx\)\)

where the fugacity is defined as z = Exp[ (μ/((Subscript[k, B]) T))] .  The gamma function is included to ensure a convenient normalization, namely  fermi[ν, z] -> z for small z.  Since the chemical potential for a dilute system is large and negative, the classical limit corresponds to small z.  Conversely, at low temperatures the chemical potential approaches the Fermi energy so that the fugacity for a Fermi system is large.  In later sections we will show that the thermodynamic properties of Fermi systems can be expressed in terms of these Fermi functions for arbitrary temperature and density.

In this section we define the Fermi function and develop rules to exploit some of its analytic properties.  We then display several fermi functions graphically.  Finally, we develop expansions that apply to both the classical and degenerate regimes.  The large z expansion, in particular, represents a generalization of the results of the preceding section.

Definitions

Fermi-Dirac functions are closely related to the polylogarithm function represented by Mathematica as PolyLog[ν,z]

fermi[ν_, z_] = 
      Integrate[((x^(ν - 1))/((z^(-1)) Exp[x] + 1)), {x, 0, ∞}, Assumptions -> {Re[ν] > 0}]/Gamma[ν]

The Fermi-Dirac functions for order {ν, (1/2), 3, (1/2)} are plotted as functions of z.  In addition, the limit ν -> ∞, which is simply the line f = z, is included.

\(\(Plot[
        Evaluate[{Table[fermi[ν, z], {ν, 1\/2, 3, 1\/2}], z}], {z, 0, 
          5}, PlotRange -> {{0, 5}, {0, 2.5}}, 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {"fugacity, z", \*""(Subscript[f, ν])[z]""}, 
        PlotLabel -> "Fermi-Dirac functions"];\)\)

Note that because z becomes large at small temperature, the divergence of Fermi-Dirac functions requires special handling via asymptotic expansions for low temperatures.

Small z expansion

Inspection of the Taylor series around z -> 0

Series[fermi[ν, z], {z, 0, 3}]

suggests that the general power series takes the form

fermiZexpansion[ν_, 
        z_, (Subscript[k, max]) : _] := -(∑_(k = 1)^((Subscript[k, max])))((((-z))^k)/(k^ν))

Below we verify that the infinite series does reproduce the desired function.

fermi[ν, z] == 
      fermiZexpansion[ν, z, ∞]

It will be useful to define a rule for small z as follows.

zSmallRule = -fermi[ν_, z] -> -fermiZexpansion[ν, z, 
            3]

Large z expansion

The degenerate Fermi gas is an important special case characterized by large z.  At low temperatures the occupation probability is nearly a square function, unity below and zero above the chemical potential, with deviations from the square function confined to a narrow region, of width approximately equal to (Subscript[k, B]) T, surrounding μ.  It is useful to separate the deviations from the completely degenerate behavior into two regions, one below and the other above μ.  The correction below μ is given by integrand1 and above μ by integrand2.  Since the integrands become very small for energies more than a few (Subscript[k, B]) T away from μ, we can extend the range of integration to -∞ for integrand1 or to ∞ for integrand2.  The exponentially small error we make in extending the energy range to-∞ can be treated in more detail, if desired.  With appropriate changes of variables, these two terms can be combined into a single integral over the range {x, 0, ∞}.   A series expansion of the integrand is then made so that the integral can be performed.

Clear[integrand1, integrand2, integrand]

integrand1[x_] := 
      (x^(ν - 1)) (1 - (1/(Exp[x - y] + 1))) /. {x -> 
            y - x}

integrand2[x_] := 
      ((x^(ν - 1))/(Exp[x - y] + 1)) /. {x -> y + x}

temp = integrand2[x] - integrand1[x]  // Simplify

integrand = 
      Normal[(Series[#1, {a, 0, 3}] &)[
            temp /. {(((x + y))^ν)_ -> 
                  (y^ν) (((((a x)/y) + 1))^ν), (((-x + y))^ν)_ -> (y^ν) (((-((a x)/y) + 1))^ν)}]] /. 
        a -> 1

correction = (∫_0^∞) 
          integrand \[DifferentialD]x

zLargeRule = -fermi[ν_, 
                  z] -> -(((((y^ν)/ν) + correction)/Gamma)[ν]) /. {Gamma[ν] -> Gamma[ν + 1]/ν} /. 
          y -> Log[z]  // Simplify

Therefore, we have obtained an asymptotic expansion of the Fermi function for large fugacity, known as Sommerfeld's lemma, that is very useful in developing low-temperature approximations.  With sufficient patience we could carry this expansion to high order, if desired.

Recursion relations for symbolic derivatives

It is also useful to note the following downward recursion relation for derivatives.

(∂_z) fermi[ν, z]

z (∂_z) fermi[ν, z] == 
      fermi[ν - 1, z]

Rules for manipulation of Fermi-Dirac integrals

Unfortunately, the Mathematica Integrate function still fails to recognize several simple variants of the Fermi-Dirac integral.  Furthermore, we would like to employ symbolic integrals without needing to specify the parameters.  Therefore, in order to guide Mathematica toward proper evaluation of these integrals, we need to temporarily unprotect the Integrate function so that we can install our own rules for the relevant types of integrals.  It is necessary to provide some rather redundant rules to ensure that all equivalent variations are recognized.

Unprotect[Integrate];

\(∫\_0\%∞\( x_\^n_\/\(Exp[β_. \ x_]\/z + 
                1\)\) \[DifferentialD]x_ := \(Gamma[n + 1]\ fermi[n + 1, z]\)\/β\^\(n + 1\)\)

\(∫\_0\%∞\( x_\^n_\/\(Exp[β_. \ x_] + 
                z\)\) \[DifferentialD]x_ := \(Gamma[n + 1]\ fermi[n + 1, z]\)\/\(z\ β\^\(n + 1\)\)\)

\(\(sumRule = ∫\_a_\%b_\((f_ + 
                  g_)\) \[DifferentialD]x_ :> ∫\_a\%b 
                f \[DifferentialD]x + ∫\_a\%b 
                g \[DifferentialD]x;\)\)

\(\(partsRule1 = ∫\_0\%∞\( x_\^m_\ y_\) \[DifferentialD]x_ :> ∫\_0\%∞\(\(-m\)\ x\^\(m - 1\)\ \(∫y \[DifferentialD]x\)\) \[DifferentialD]x;\)\)

\(\(partsRule2 = ∫\_0\%∞\( x_\^m_\ y_\) \[DifferentialD]x_ :> ∫\_0\%∞\(-\(\(x\^\(m + 1\)\ ∂\_x y\)\/\(m + 1\)\)\) \[DifferentialD]x;\)\)

Protect[Integrate];

gammaRule = n_ Gamma[n_] -> Gamma[n + 1];

Chemical potential and fugacity

Temperature dependence of chemical potential and occupancy

The chemical potential is related to the temperature and density by

Clear[eq1];
eq1[z_, τ_] = 
      PowerExpand[(Subscript[n, Q]) == fermi[(3/2), z] /. 
              lambdaToT /. {T -> ((τ (Subscript[ε, F]))/(Subscript[k, B]))} /. eFrule]

where we express energy and chemical potential in units of (Subscript[ε, F]) and temperature as τ = T/(Subscript[T, F]).  Although this equation usually must be solved numerically, it will be useful to develop symbolic relationship valid for either low or high temperature limits.  In both cases we must be careful to select the proper root and not to use the expansion too far outside its range of validity.

zSmallEq = eq1[z, τ] /. zSmallRule

zhigh  =  
        z /. Solve[zSmallEq, 
              z][[1]];

zLargeEq = (eq1[z, τ] /. zLargeRule // Simplify)

zlow = 
        Exp[Log[z] /. Solve[zLargeEq, 
                Log[z]][[2]]];

Interpolation function for chemical potential

In this section we construct an interpolation function from numerical solutions to the equation related chemical potential and temperature.  Then we compare with the limiting formulas derived in the preceding section.

zSolution[τ_] := 
      z /. FindRoot[Evaluate[eq1[z, τ]], {z, 0.6, 0.4}]

points = 
        Prepend[Table[{τ, τ Log[zSolution[τ]]}, {τ, 0.2, 
              5, 0.2}], {0. , 1. }];

Clear[chemicalPotential];
chem = Interpolation[points];
chemicalPotential[τ_] := chem[τ]

\(\(DisplayTogether[
        Plot[chemicalPotential[τ], {τ, 0, 3}, 
          PlotStyle -> RGBColor[1, 0, 0]], 
        Plot[τ\ Log[zlow], {τ, 0, 0.45}, 
          PlotStyle -> GrayLevel[0]], 
        Plot[τ\ Log[zhigh], {τ, 1, 3}], 
        ListPlot[points, PlotStyle -> AbsolutePointSize[4.0]], 
        PlotLabel -> "Chemical Potential", 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {\*""T/(Subscript[T, F])"", \*""μ/((Subscript[k, B])(Subscript[T, F]))""}, 
        PlotRange -> {{0, 3}, {\(-6\), 1}}];\)\)

The Fermi energy is defined to be the chemical potential when T = 0 for a specified density, such that \(TraditionalForm\`T = \(0\ ⟹μ = ε\_F\)\); similarly, T = (Subscript[T, F]) ⟹ μ = 
          0.  Thus, the chemical potential for a Fermi gas is positive for T < (Subscript[T, F]) and is negative for T >  (Subscript[T, F]).  A Fermi gas is described as completely degenerate for T ≪ (Subscript[T, F]).  For very large temperatures, T ≫ 
        (Subscript[T, F]) ⟹μ ≪ -(Subscript[k, B]) T, the chemical potential becomes large and negative such that the nondegenerate Fermi gas approaches classsical behavior.  For intermediate temperatures, T ∼ 
        (Subscript[T, F]) ⟹ μ  ∼  0, the partially degenerate Fermi gas can be expected to display thermodynamic properties intermediate between the degenerate and classical situations.

Next we compare the interpolation formula for fugacity with the limiting behaviors derived in the preceding section.

fugacity[τ_ /; τ > 0] := 
      Exp[chemicalPotential[τ]/τ]

DisplayTogether[
        Plot[fugacity[τ], {τ, 0.1, 5}, 
          PlotStyle -> RGBColor[1, 0, 0]], 
        Plot[zlow, {τ, 0.01, 0.45}, PlotStyle -> GrayLevel[0]], 
        Plot[zhigh, {τ, 0.5, 3}], PlotLabel -> "Fugacity", 
        PlotRange -> {Automatic, {0, 10}}, GridLines -> Automatic, 
        Frame -> True, 
        FrameLabel -> {"τ", "z"}];

The small τ formula appears to be accurate enough for τ ≲ 0.4 while the large τ formula is good for τ ≳ 1.  However, the divergence of the fugacity as τ -> 0 requires special handling of the low-temperature limits for many thermodynamic functions.  We will employ the large z expansion to provide well-behaved limiting formulas as needed.  We also find that small inaccuracies in the interpolation formula can sometimes produce small imaginary parts in some of the thermodynamic functions.  These spurious contributions are negligible for τ > 1 and easily suppressed using the Chop function, but the expansions are needed to obtain accurate results for smaller τ.

Mean occupation numbers for several temperatures

Clear[temp];
\(\(temp[τ_] := \(n\_FD\)\&_ //. {k\_B -> 1, 
            T -> τ, μ -> 
              chemicalPotential[τ]};\)\)
\(\(nFDplot[τ_] := 
        Plot[temp[τ], {ε, 0, 2}, 
          PlotLabel -> "Mean Occupation Number", 
          GridLines -> None, Frame -> True, 
          FrameLabel -> {\*""ε/(Subscript[ε, F])"", \*""\!\(n\&_\)""}, 
          DisplayFunction -> Identity];\)\)

Show[{nFDplot[0.01], nFDplot[0.5], nFDplot[1], nFDplot[1.5]}, 
        DisplayFunction -> $DisplayFunction];

\(\(Plot[chemicalPotential[τ], {τ, 0, 3}, 
        PlotLabel -> "Chemical Potential", 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {\*""T/(Subscript[T, F])"", \*""μ/((Subscript[k, B])(Subscript[T, F]))""}, 
        PlotRange -> {{0, 3}, {\(-6\), 1}}];\)\)

\(\(Plot3D[temp[τ], {ε, 0, 2}, {τ, 0.001, 1}, 
        PlotRange -> {{0, 2}, {0, 1}, {0, 1}}, PlotPoints -> 30, 
        PlotLabel -> "Occupation Number", 
        AxesLabel -> {\*""ε/((Subscript[k, B])(Subscript[T, F]))"", \*""T/(Subscript[T, F])"", "n"}];\)\)

Rules for derivatives with respect to temperature and volume

The following rules are provided to facilitate development of thermodynamic relationships.

dzdT = -((3 z fermi[(3/2), z])/(2 T fermi[(1/2), 
                  z]));

dLambdaDT = 
        Dt[λ, T, 
            Constants -> {𝒩, 
                V}] -> -(λ/(2 T));

eq2 = (Subscript[n, Q]) == fermi[(3/2), z]

dzRules = { 
            Solve[Dt[eq2, V, 
                Constants -> {𝒩, T, λ}], 
              Dt[z, V, Constants -> {𝒩, T, λ}]], 
            Solve[Dt[eq2, T, Constants -> {𝒩, V}], 
              Dt[z, T, Constants -> {𝒩, V}]]} /. 
          dLambdaDT  // Flatten

DtRules = Join[{dLambdaDT}, dzRules]

Thermodynamic Functions

The thermodynamic functions for an ideal Fermi system are obtained by integrating single-orbital functions against the density of energy states.  It is useful to express these integrals in terms of fermi functions, but we will also require forms suitable for numerical evaluation over a wide range of temperatures, including regions where the fugacity diverges.  Therefore, we use ordinary symbols, such as p for pressure, for symbolic expressions and include a tilde, such as \(TraditionalForm\`\(p\&~\)\), to represent functions of τ in reduced form suitable for plotting.

Grand potential

It is convenient to express the grand potential in terms of the so-called q-potential defined as \(TraditionalForm\`q\  = \ \(\(\(p\ V\)\/\(\(k\_B\) 
                  T\)\)\ ⟹\ 𝒢 = \(-k\_B\) T\ q\)\).

q  = (((∫_0^∞) 𝒟[ε] Log[1 + 
                            z Exp[-β ε]] \[DifferentialD]ε /. partsRule2 // Simplify) //. TtoLambda // 
            PowerExpand) /. lambdaToZ

𝒢 = -(Subscript[k, B]) T q

Mechanical equation of state

Degeneracy pressure

The mechanical equation of state is obtained directly from the grand potential, 𝒢  =  p V.  It is useful to express the mechanical equation of state in the form ((p V)/(N (Subscript[k, B]) T)) = B[τ], where B is sometimes called the compressibility factor.

p = -(𝒢/V)

B = ((p V)/(𝒩 (Subscript[k, B]) T))

Note that in these expressions the dependence upon density is implicit in z.  For low temperatures it is necessary to employ asymptotic expansions of the Fermi-Dirac functions with respect to z.

Bhigh[τ_] = B /. zToTau;
\(Blow[τ_]\  = \(B /. zLargeRule\) /. zToTau // Simplify; 
    B\&~[τ_] = 
      If[τ > 0.1, Chop[Bhigh[τ]], Blow[τ]];\)

\(\(compressibilityPlot = 
        Plot[B\&~[τ], {τ, 0.2, 5}, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {"τ", \*""pV/(𝒩(Subscript[k, B])T)""}, 
          PlotLabel -> "Mechanical Equation of State"];\)\)

The compressibility factor for an ideal Fermi gas approaches the classical limit from above at high temperatures.  The Pauli exclusion principle, which prevents two fermions from occupying exactly the same quantum state, produces a type of quantum repulsion which tends to keep fermions apart.  This effect increases the pressure for a given density and temperature, particularly at large density or low temperature.  The enhancement of the pressure over classical expectations, called the degeneracy pressure, is rather small for T > (Subscript[T, F]) but increases rapidly as the temperature falls below 0.5 (Subscript[T, F]).

Dimensionless form

The dependence is displayed in a dimensionless form.

pReduced = p /. 
          zToTau /. {g -> 1, 𝒩 -> 1, 
          (Subscript[k, B]) -> 1, τ -> T (V^(2/3))}

Plot[Evaluate[Table[pReduced, {T, 0.5, 2, 0.5}]], {V, 0.1, 2}, 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {"V", "p"}, 
        PlotLabel -> "Isotherms for ideal Fermi systems"];

Compressibility and expansivity

(Subscript[κ, T]) = (-(1/(V Dt[p, V, 
                        Constants -> {𝒩, 
                            T, λ}])) /. dzRules) /. 
          lambdaToZ  // Simplify

α = (Subscript[κ, T]) Dt[p, T, 
                  Constants -> {𝒩, V}]  /. dzRules /. 
            lambdaToZ  // Simplify  // Apart

Virial expansion

An expansion is developed for ((p V)/(N (Subscript[k, B]) T))  in terms of the quantum concentration, here designated x.

zToX = Solve[x == fermi[(3/2), z] /. zSmallRule /. 
              (z^3) -> 0, z][[1]] // 
        Simplify

virialExpansion = (Series[#1, {x, 0, 2}] &)[B /. lambdaToZ /. 
              zSmallRule /. zToX] // Simplify

virialExpansion  // N // Chop

Thermal equation of state

Degeneracy energy

The thermal equation of state for a nonrelativistic three-dimensional gas is obtained directly from the scaling relation U  =  (3/2) p V; other systems display similar scaling relations with different coefficients.

U = (3/2) p V

Recognizing that U can be expressed in terms of B, we define a function suitable for numerical evaluation of the entire range of τ as follows.

\(\(U\&~[τ_] = \(3\/2\) τ\ B\&~[τ];\)\)

\(\(plot[U] = 
        Plot[{\(3\ τ\)\/2, U\&~[τ]}, {τ, 0, 2}, 
          PlotStyle -> {Dashing[{0.02, 0.02}], {}}, 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, F])"", \*""U/((Subscript[Nk, B])(Subscript[T, F]))""}, PlotLabel -> "Internal Energy"];\)\)

The internal energy for an ideal Fermi gas also approaches the classical limit from above.  At low temperatures the energy of a completely degenerate Fermi gas approaches a constant value, (3/5) N (Subscript[ε, F]), whereas the energy of a classical gas is proportional to temperature and goes to zero.  If we imagine that the system is assembled by adding particles one at a time, then at low temperatures each particle will generally be added to the lowest available energy level.  However, because the Pauli exclusion principle prohibits more than one fermion from occupying each (fully specified) orbital, each particle must be added to the top of the energy distribution.  Therefore, even at absolute zero there will be a substantial degeneracy energy in a Fermi gas.   Recall the expansion we developed for a degenerate Fermi gas.

Uexpansion

Isochoric Heat Capacity

A symbolic expression for the isochoric heat capacity is obtained by differentiating the internal energy

(Subscript[C, V]) = Dt[U, T, Constants -> {𝒩, V}] /. 
              DtRules /. lambdaToZ  // 
          Collect[#1, {𝒩, (Subscript[k, B])}] &  // Simplify

and is expressed in a form suitable for numerical evaluation as follows.

CVhigh[τ_]  =  (Subscript[C, V]) /. zToTau;
CVlow[τ_]  =  (Subscript[C, V]) /. zLargeRule /. zToTau // 
          Simplify;
\(C\&~\_V[τ_] := 
      If[τ > 0.1, Chop[CVhigh[τ]], CVlow[τ]]\)

\(\(plot[C\_V] = 
        Plot[C\&~\_V[τ]\/\(𝒩\ k\_B\), {τ, 0, 5}, 
          PlotRange -> {{0, 3}, {0, 1.5}}, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, F])"", \*""(Subscript[C, V])/((Subscript[Nk, B]))""}, 
          PlotLabel -> "Isochoric Heat Capacity"];\)\)

The isochoric heat capacity is linear for small temperatures: it vanishes at absolute zero because with no vacant low-energy single-particle states the completely degenerate Fermi gas cannot absorb increments of energy less than (Subscript[ε, F]).  As vacancies develop at higher temperature the heat capacity rises until it approaches classical equipartition at high temperatures.  Thus, it is instructive to develop the following limiting cases.

(∂_τ) Uexpansion

Normal[
          Series[(Subscript[C, V]) /. zSmallRule /. 
              z -> zhigh, {τ, ∞, 2}]] // 
        Collect[#, {𝒩, (Subscript[k, B])}] & // N

Isobaric Heat Capacity

The isobaric heat capacity is defined as the derivative of enthalpy with respect to temperature for constant pressure.  The enthalpy, U  +  p V, for an ideal nonrelativistic gas is quite simple, reducing to ((5 p V)/2), which allows the isobaric heat capacity to be obtained easily.  Alternatively, it can also be obtained using a more general approach which relates the difference between (Subscript[C, p]) and (Subscript[C, V]) to the isobaric expansivity and the isothermal compressibility of the system.  The latter approach is more appropriate here.

(Subscript[C, p]) = 
      (Subscript[C, V]) + ((V T (α^2))/(Subscript[κ, T]))  // Simplify

Once again, we need a form suitable for numerical evaluation.

Cphigh[τ_]  =  (Subscript[C, p]) /. zToTau;
Cplow[τ_]  =  (Subscript[C, p]) /. zLargeRule /. zToTau // 
          Simplify;
\(C\&~\_p[τ_] := 
      If[τ > 0.1, Chop[Cphigh[τ]], Cplow[τ]]\)

\(\(plot[C\_p] = 
        Plot[C\&~\_p[τ]\/\(𝒩\ k\_B\), {τ, 0, 5}, 
          PlotRange -> {{0, 5}, {0, 3}}, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, F])"", \*""(Subscript[C, p])/((Subscript[Nk, B]))""}, 
          PlotLabel -> "Isobaric Heat Capacity"];\)\)

The isobaric heat capacity behaves much like the isochoric in that it approaches zero linearly at low temperature or the classical equipartition value asymptotically for high temperature.

Normal[
          Series[(Subscript[C, p]) /. zSmallRule /. 
              z -> zhigh, {τ, ∞, 2}]] // 
        Collect[#, {𝒩, (Subscript[k, B])}] & // N

Series[(Subscript[C, p]) /. zLargeRule /. 
                zToT /. {μ -> μExpansion, 
                T -> (Subscript[ε, F]) τ/(Subscript[k, B])} // 
            PowerExpand // Simplify, {τ, 0, 4}] // 
      Collect[#, {𝒩, (Subscript[k, B])}] &

Compare (Subscript[C, p]) and (Subscript[C, V])

For a classical ideal gas we would find (Subscript[C, V ])  =  (3/2) N (Subscript[k, B]), (Subscript[C, p])  =  (5/2) N (Subscript[k, B]), and ((Subscript[C, p])/(Subscript[C, V]))  =  (5/3)where p (V^γ) is constant for adiabats and where γ = (5/3).  For the Fermi gas we find a somewhat more complicated expression for γ in terms of Fermi functions, but simple expansions can be developed for the degenerate or the dilute limits.

γ = ((Subscript[C, p])/(Subscript[C, V]))  // Simplify

γhigh[τ_] = γ /. zToTau // Simplify;
γlow[τ_] = γ /. zLargeRule /. zToTau // 
          Simplify;
\(\(γ\&~[τ_] = 
        If[τ > 0.1, 
          Chop[γhigh[τ]], γlow[τ]];\)\)

\(\(plot[γ] = 
        Plot[γ\&~[τ], {τ, 0, 5}, 
          PlotRange -> {{0, 5}, {0, 2}}, GridLines -> Automatic, 
          Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, F])"", \*""(Subscript[C, p])/(Subscript[C, V])""}, 
          PlotLabel -> "Ratio of Heat Capacities"];\)\)

The ratio between the principal heat capacities approaches unity quadratically for low temperature, or the classical equipartition value asymptotically for high temperature.

Normal[
        Series[γ /. zSmallRule /. 
            z -> zhigh, {τ, ∞, 2}]] // N

Series[γ /. zLargeRule /. 
              zToT /. {μ -> μExpansion, 
              T -> (Subscript[ε, F]) τ/(Subscript[k, B])} // PowerExpand // 
        Simplify, {τ, 0, 4}]

Entropy and free energy

Symbolic forms

S = ((U + p V - 𝒩 μ)/T) /. muToZ  // 
        Collect[#1, {𝒩, (Subscript[k, B]), T}] &

F = U - T S // Collect[#1, {𝒩, (Subscript[k, B]), T}] &

Plot reduced forms

Shigh[τ_] = 
        (S/(𝒩 (Subscript[k, B]))) /. zToTau;
Slow[τ_] = (S/(𝒩 (Subscript[k, B]))) /. zLargeRule /. 
            zToTau // Simplify;
\(\(S\&~[τ_] = 
        If[τ > 0.1, Chop[Shigh[τ]], Slow[τ]];\)\)

Fhigh[τ_] = ((F τ)/(𝒩 (Subscript[k, B]) T)) /. 
          zToTau;
Flow[τ_] = ((F τ)/(𝒩 (Subscript[k, B]) T)) /. 
              zLargeRule /. zToTau // Simplify;
\(\(F\&~[τ_] = 
        If[τ > 0.1, Chop[Fhigh[τ]], Flow[τ]];\)\)

\(\(plot[S] = 
        Plot[S\&~[τ], {τ, 0, 5}, 
          PlotLabel -> "Reduced Entropy", 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, F])"", \*""S/((Subscript[Nk, B]))""}, DisplayFunction -> Identity];\)\)

\(\(plot[F] = 
        Plot[F\&~[τ], {τ, 0, 5}, 
          PlotLabel -> "Reduced Free Energy", 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {\*""T/(Subscript[T, F])"", \*""F/((Subscript[Nk, B])(Subscript[T, F]))""}, DisplayFunction -> Identity];\)\)

Show[ 
        GraphicsArray[{{plot[U], plot[S]},  {plot[(Subscript[C, V])], plot[F]}}, 
          GraphicsSpacing -> {-0.1, 0.05}], 
        DisplayFunction -> $DisplayFunction] ;

Degenerate Fermi gas

Expansions for low temperatures

Here we develop small τ expansions for the principal thermodynamic functions for a nearly degenerate Fermi gas.  We start by expressing temperature and chemical potential in reduced form (scaled to the Fermi energy) and put the dimensions back later.

temp = eq1[z, τ] /. 
              zLargeRule //. {z -> Exp[(μ/τ)]}  // 
          PowerExpand // Simplify

eq3 = 
      Simplify[Normal[(Thread[Series[#1, {μ, 1, 2}], Equal] &)[
            temp]]]

musol = Solve[eq3, μ]  // Simplify

muSmallTauRule = μ -> (Normal[(Series[#1, {τ, 0, 
                        4}] &)[μ (Subscript[ε, F]) /. 
                  musol[[2]]]] // 
            Collect[#, {(Subscript[ε, F]), τ}] &)

uDegenerate = U /. zLargeRule /. 
                  zToT  /. {(Subscript[k, B]) -> 1, 
                  T -> τ (Subscript[ε, F])} /. muSmallTauRule // 
            PowerExpand // Series[#1, {τ, 0, 4}] & // 
        Collect[#1, {𝒩, (Subscript[ε, F]), τ}] &

cvDegenerate = ((Subscript[k, B])/(Subscript[ε, F])) (∂_τ) 
            uDegenerate //Collect[#1, {𝒩, 
              (Subscript[k, B]), τ}] &

pDegenerate = ((2 uDegenerate)/(3 V)) // 
        Collect[#1, {𝒩, 
              V, (Subscript[ε, F]), τ}] &

The easiest way to evaluate the Helmholtz free energy is to remember that the Gibbs free enthalpy for a single-component system is proportional to the chemical potential and then to use the thermodynamic relationship between free energy and free enthalpy.

fDegenerate = 𝒩 μ  - 
            pDegenerate V /. {muSmallTauRule} // 
        Collect[#1, {𝒩, (Subscript[ε, F]), τ}] &

sDegenerate = (Subscript[k, B]) ((uDegenerate - fDegenerate)/(τ (Subscript[ε, F])))//Collect[#1, {𝒩, 
              (Subscript[k, B]), τ}] &

α /. zLargeRule /. zToT  /. {(Subscript[k, B]) -> 1, 
                T -> τ (Subscript[ε, F]) } /. muSmallTauRule // 
          PowerExpand // Series[#1, {τ, 0, 4}] & // 
      Collect[#1, {𝒩, (Subscript[ε, F]), τ}] &

(Subscript[κ, T]) /. zLargeRule /. zToT  /. {(Subscript[k, B]) -> 1, 
                T -> τ (Subscript[ε, F])} /. muSmallTauRule // 
          PowerExpand // Series[#1, {τ, 0, 4}] & // 
      Collect[#1, {𝒩, 
            V, (Subscript[ε, F]), τ}] &

Examples

Conduction electrons in copper

In metals conduction electrons are free to move throughout the material.  It is then useful to approximate the behavior of conduction electrons by a Fermi gas.   However, since the conduction electrons are not really free, the effect of the mean field can be represented by a modification of the effective mass associated with the electron.  The value below is based upon a fit to specific heat data for low temperatures.  It is also of interest to recognize that the enormous pressure associated with the degenerate electron gas must be balanced by the attractive mean field provided by the lattice ions.

copperValues = 
        Join[{g -> 2, 𝒩 -> V density, 
            density -> 8.50 (10^28), effectiveMass -> 1.39, 
            (Subscript[T, D]) -> 345, 
            m -> ((effectiveMass ElectronMass)/Kilogram)}, 
          FundamentalConstants];

(Subscript[ε, F])[copper]  =  (Subscript[ε, F])  /. eFrule //. 
        copperValues

(Subscript[ε, F])[copper]  Convert[Joule, ElectronVolt]

(Subscript[ε, F])[copper]/(Subscript[k, B]) /. copperValues

pDegenerate  Convert[Pascal, 
            Atmosphere] /. {(Subscript[ε, F]) -> (Subscript[ε, F])[
              copper], τ -> 0} //. copperValues

At room temperature, conduction electrons in a metal already constitute a highly degenerate Fermi gas; hence, the electronic contribution to the heat capacity is linear wrt temperature.  On the other hand, the lattice contribution, as given by the Debye model, is cubic wrt temperature.  Therefore, at temperatures well below the Debye temperature we can combine these two terms as follows.

cvCopperElectrons = 
      Expand[(cvDegenerate/(𝒩 (Subscript[k, B]))) //. {(Subscript[ε, F]) -> (Subscript[ε, F])[
                  copper], τ -> (((Subscript[k, B]) T)/(Subscript[ε, F]))} //. 
          copperValues]

cvCopperLattice = 
      (12/5) (π^4) ((((T/(Subscript[T, D]))))^3) /. copperValues  // N

cvCopper = (Collect[#1, {Joule, Kelvin}] &)[
        Expand[(1/Joule) ((cvCopperElectrons + 
                    cvCopperLattice) AvogadroConstant BoltzmannConstant Mole Kelvin)]]

Clearly, the electronic contribution is revealed by the coefficient of the linear term, whereas the lattice contribution dominates the cubic term.  These effects are easily separated by plotting (Subscript[C, V])/T versus (T^2).   The intercept then represents the electronic and the slope the lattice contribution.

cvTemp[x_] := Simplify[(cvCopper/T)] /. (T^2) -> x;
\(\(Plot[1000\ cvTemp[x], {x, 0, 10}, 
        PlotRange -> {{0, 10}, {0, 1.2}}, GridLines -> Automatic, 
        Frame -> True, 
        FrameLabel -> {\*""(T^2)  ((kelvin^2))"", \*""(Subscript[C, V])/T  (mJ/mole/(kelvin^2))""}, 
        PlotLabel -> "specific heat for copper"];\)\)

Nuclear matter

The atomic nucleus consists of protons and neutrons, which are both spin - (1/2) particles of almost equal mass.  It is useful to consider protons and neutrons to be states of the same particle, the nucleon, differing only in an internal quantum number called isospin.  Nuclear matter is a theoretical system consisting of equal numbers of protons and neutrons with the Coulomb interaction turned off.  Thus, the intrinsic degeneracy factor for momentum states in nuclear matter is g = 4.  The density of nuclear matter is based upon the central density of large nuclei, which is approximately constant.

nmValues = 
        Join[{g -> 4, 
            m -> Convert[AtomicMassUnit, Kilogram]/Kilogram, 
            density -> 
              0.16 SI[
                  (1/(((Femto Meter))^3))] (Meter^3), 𝒩 -> V density}, FundamentalConstants];

density Convert[(1/(Meter^3)), (1/(((Femto Meter))^3))] //. 
      nmValues

(Subscript[ε, F])[nm]  =  (Subscript[ε, F]) /. eFrule //. 
        nmValues

(Subscript[ε, F])[nm] Convert[Joule, ElectronVolt]

(Subscript[ε, F])[nm]/(Subscript[k, B]) //. nmValues  // PowerExpand

pDegenerate Convert[Pascal, Atmosphere] /. 
          eFrule  /. {τ -> 0} //. nmValues

Convert[((Mega ElectronVolt)/(((Femto Meter))^3)), 
      Atmosphere]

The density of nuclear matter is about ((1/6)) per cubic femtometer.  The Fermi energy of about 37 MeV corresponds to a temperature of about 4 ×(10^11) kelvin.  Clearly, atomic nuclei at room temperature are completely degenerate.  The enormous pressure within the degenerate Fermi gas must be balanced by the strong nuclear forces that bind the nucleus together.  The name "strong interaction" is obviously quite appropriate!

Problems

(((Subscript[n, c]) ((Subscript[N, D]) - (Subscript[n, D])))/(Subscript[n, D])) ≈ 
      A[T] Exp[-((Subscript[E, D])/((Subscript[k, B]) T))]

\(ε\  = \ \(\(ℏ\^2\) k\^2\)\/\(2\ m\)\  ∓ \ \(μ\_0\) B\)

\(χ\  ≈ \ \ μ\_0\^2\ \(∫\_0\%∞\( 𝒟\^′\)[ε]\ n\&_[ε] dε\)\)

\(R\  = \ \(\(4  π\ m\ \(k\_B\) 
              T\)\/\((2  π\ ℏ)\)\^3\) \(∫\_W\%∞ 
              Log[1\  + \ 
                Exp[\ \(μ - ε\_z\)\/\(\(k\_B\) T\)]] dε\_z\)\)

v  =  (Sqrt[(Subscript[((∂(p/∂)ρ)), S])])

\(v\^2 = \ \(\(\(5 \( k\_B\) T\)\/\(3  m\)\) 
          fermi[5\/2, z]\/fermi[3\/2, z]\  = \ \(5\/9\) \(u\^2\)\&_\)\)