Statistics of Occupation Numbers

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

Introduction

Consider an ideal quantum fluid composed of identical fermions or bosons under conditions for which their mutual interactions may be neglected.  Under those conditions, the total energy (Subscript[E, ν]) for state ν can be expressed as a summation

\(\(\(\ \ \)\(E\_ν = \ ∑\+α n\_α[ν]\ ε\_α\)\)\)

where (Subscript[ε, α]) is the single-particle energy for orbital α and (Subscript[n, α]) is the number of particles in that orbital when the system is found in state ν.  Rather than specifying the coordinates or other quantum numbers for each particle in the system, the occupation representation specifies how many, but not which, particles are found in each orbital.  This is the most natural representation of the states for a system of identical particles and avoids the complications of enforcing permutation symmetry and evaluating degeneracy factors.  Thus, the many-body wave function is fully determined by the set ν = {(Subscript[n, α])} of occupation numbers for each orbital.

The occupancy of single-particle states in a many-body system is a statistical quantity subject to thermodynamic fluctuations.  Since the number of particles occupying any state is indeterminate, each orbital can be visualized as a subsystem in equilibrium with a reservoir of energy and particles; of course, the reservoir consists of all particles occupying orbitals other than the one being analyzed.  Therefore, the statistics of occupancy is susceptible to analysis using the grand canonical ensemble.  For an ideal quantum fluid, the grand partition function takes the form

\(𝒵 = ∑\+\(N = 0\)\%∞\(\(∑\+ν\)\&~ 
            Exp[\(-β\)\ \((E\_ν - \ μ\ N\_ν)\)]\)\)

where μ is the chemical potential and β=(1/((Subscript[k, B])T)).  The tilde on the summation over states ν indicates that for each value of N the summation over occupation numbers {(Subscript[n, α])} is restricted to states with total particle number (Subscript[N, ν]) = N where

\(\(\(\ \)\(N\_ν = \ \(\(∑\+α\)\(n\_α[ν]\)\(\ \)\)\)\)\)

The temperature and chemical potential are determined by the constraints upon energy and density, respectively.  Fortunately, the separability of the total energy and particle number for an ideal quantum fluid allows the grand partition function to be expressed in the form

\(𝒵 = ∏\+α\((∑\+\(n\_α\)Exp[\(-β\)\ \(n\_α\) \((ε\_α - \ μ)\)])\)\)

where each sum over (Subscript[n, α]) is now independent, restricted only by the Pauli exclusion principle for fermions.  Thus, the summation over total particle number N effectively eliminates the constraint upon the occupation numbers.  Therefore, the grand partition function for the many-body system factorizes into products of single-orbital grand partition functions

\(\(\(\ \ \)\(E\_ν = \ \(∑\+α n\_α[ν]\ ε\_α\ \ \ ⟹\ \ \ 𝒵\  = ∏\+α 𝒵\_α\)\)\)\)

of the form

\(𝒵\_α = ∑\+\(n\_α\)Exp[\(-β\)\ \(n\_α\) \((ε\_α - \ μ)\)]\)

where μ is the chemical potential established by the reservoir and (Subscript[n, α]) is the number of particles in orbital α.  Finally, the thermodynamic functions are obtained using the grand potential

𝒢  =  (∑_α) (Subscript[𝒢, α])
(Subscript[𝒢, α])  =  -(Subscript[k, B]) 
        T Log[(Subscript[𝒵, α])]

which for an ideal quantum fluid reduces to a summation over the contributions made by each orbital.

The statistics of quantum systems are intimately related to the intrinsic angular momentum (spin) of their constituents.  Systems composed of identical particles with integer spin, generically classified as bosons, obey Bose-Einstein (BE) statistics for which the many-body wave function is symmetric with respect to permutation of particle coordinates and the summation over occupation numbers extends from 0 to ∞.  Conversely, systems composed of identical particles with half-integral spin, generally classified as fermions, obey Fermi-Dirac (FD) statistics for which the many-body wave function is antisymmetric with respect to particle exchange and the summation over occupation numbers is limited to 0 or 1 particles in each unique (fully specified) single-particle orbital.  The prohibition against more than one fermion occupying any orbital is known as the Pauli exclusion principle and has profound consequences for both the structure and the thermodynamics of Fermi-Dirac systems.  Finally, it is instructional to compare BE and FD statistics with classical or Maxwell-Boltzmann (MB) statistics for which permutation symmetry is restored using the Gibbs prescription even though MB statistics do not accurately describe real quantum systems.  The summation over occupation numbers extends from 0 to ∞ for the MB distribution also, but the Gibbs prescription must be applied during the factorization of the grand partition function, as discussed below.

For each distribution, the probability that a given orbital contains n particles has the form

\(P[n\_α]\  = \ 
      Exp[\(-β\)\ \(n\_α\) \((ε\_α - \ μ)\)]/𝒵\_α\)

The mean occupancy of that orbital is then

\(\(n\_α\)\&_\  = \ \(∑\+\(n\_α\)P[n\_α] 
            n\_α\  = \ \(-∂\_μ 𝒢\_α\)\)\)

and the variance in occupancy is

(Subscript[σ, α]^2)= \(\(\((n\_α - \(n\_α\)\&_\ )\)\^2\)\&_\)= (Subscript[k, B])T \(∂\_μ\( n\_α\)\&_\)

Note that P[(Subscript[n, α])] depends implicitly upon T and μ and upon any external variables, such as volume, which govern the values of (Subscript[ε, α]).  Hence, when comparing the various probability distributions, often it is useful to eliminate μ in favor of \(TraditionalForm\`\(\(\(n\_α\)\&_\)\(\ \)\)\).  Finally, it is useful to express the relative fluctuation in occupation numbers as \(TraditionalForm\`Ϟ\_α = \ σ\_α\/\(n\_α\)\&_\).

The chemical potential is determined by the mean density of the system.  We can obtain the total particle number from the occupation probabilities using the summation

\(N\&_ = \(∑\+α\ \(n\_α\)\&_\ \  = \ \(\(-∂\_μ 𝒢\)\  = \ ∫\_0\%∞\ \(\(n[ε]\)\&_\) 𝒟[ε] dε\)\)\)

For a continuous spectrum of single-particle energy levels, the summation over occupation numbers is accomplished through integration against the density of states, 𝒟[ε].  Although the mean particle number is a statistical variable within the grand canonical ensemble, its  distribution is nevertheless extremely sharp for macroscopic systems.  Hence, we can generally assume that N and V are determined by external constraints and thereby related the chemical potential directly to the density.

In this notebook we concentrate on the statistics of occupation numbers, comparing the BE and FD distributions with the classical MB model.  A summary of the symbolic results is provided at the end of this notebook.  The thermodynamics of ideal Fermi and Bose gases are developed in more detail in the notebooks fermi.nb and bose.nb, respectively.

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])	= grand partition function for single orbital
	(Subscript[𝒢, 1])	= grand potential for single orbital
	(Subscript[S, 1])	= entropy for single orbital
	P[n]	= probability for finding n particles in a particular state
	\(TraditionalForm\`n\&_\)	= mean occupation number
	σ	= standard deviation for occupancy
	Ϟ	= relative fluctuation in occupancy = \(TraditionalForm\`σ/n\&_\)
	(Subscript[T, c])	= critical temperature for Bose-Einstein condensation
	(Subscript[ε, F])	= Fermi energy
	(Subscript[k, F])	= Fermi momentum
	
Many quantities below also include either superscripts or subscripts identify the symmetry type as FD, BE, or MB.

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[k, B])];Symbolize[(Subscript[n, x])_];Symbolize[(Subscript[μ, x])_];Symbolize[(Subscript[λ, x])_];Symbolize[(Subscript[T, x])_];Symbolize[(Subscript[σ, x])_];Symbolize[(Subscript[𝒵, 1]^x)_];Symbolize[(Subscript[𝒢, 1]^x)_];Symbolize[(Subscript[S, 1]^x)_]

Symbolize[\(n\&_\)];Symbolize[\(\(n\_x_\)\&_\)]

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

Error messages

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

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

Establish  transformation rules

\(\(MyAssumptions = {m > 0, V > 0, g > 0, k\_B > 0, n_ > 0, n\&_ > 0, 
          T > 0, ℏ > 0, λ > 0, β > 0, 𝒩 > 
            0};\)\)

MySimplify = 
        Simplify[#, MyAssumptions] &;
MyFullSimplify = FullSimplify[#, MyAssumptions] &;

Rules for changing variables

\(\(TtoY = {μ -> 0, a_. \ k\_B\ T -> a\/β, 
          a_. \/\(k\_B\ T\) -> a\ β, 
          a_. \ β\ ε -> a\ y};\)\)
TtoLambda = {T -> ((2 π (ℏ^2))/(m (Subscript[k, B]) (λ^2)))};
lambdaToT = {λ -> (Sqrt[((2 π (ℏ^2))/(m (Subscript[k, B]) T))])};
lambdaToQC = {λ -> (((((g V (Subscript[n, Q]))/𝒩)))^(1/3))};

muToZ = μ -> (Subscript[k, B]) T Log[z];

Thermodynamic functions

Given a single-orbital grand partition function or grand potential, these thermodynamic functionals evaluate the grand potential, the mean energy, the single-orbital entropy, the mean occupancy, and occupation variance for that orbital.

grand[𝒵_] := -(Subscript[k, B]) T Log[𝒵 /. {β -> (1/((Subscript[k, B]) T))}];
energy[𝒵_] := 
        (Subscript[k, B]) (T^2) (∂_T) 
              Log[𝒵 /. {β -> 
                    (1/((Subscript[k, B]) T))}];
entropy[𝒢_] := -(∂_T) 𝒢;

nOccupy[𝒢_] := -(∂_μ) 𝒢;
nVariance[𝒢_] := 
        (Subscript[k, B]) T (∂_μ) nOccupy[𝒢];

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

Maxwell-Boltzmann distribution

Evaluate mean occupation number and chemical potential

Although the semi-classical Maxwell-Boltzmann distribution doesn't actually pertain to any real quantum system, it is nevertheless instructive to compare the FD and BE distributions to this classical model.  Therefore, we begin with that model.  Classically, the degeneracy for a state of N identical particles described by occupation numbers ν = {(Subscript[n, α])} would be \(TraditionalForm\`g\_ν\^′ = \(N!\)/\(∏\+α\( n\_α!\)\)\), but we must now divide by N! to resolve the Gibbs paradox so that (Subscript[g, ν]) = 
        1/(∏_α) (Subscript[n, α])! becomes the proper statistical weight for the MB distribution.  Then, as discussed in the introduction, the grand partition function can be factored, such that

(Subscript[𝒵, 1]^MB) = (∑_(n = 0)^∞) Exp[-n β (ε -  μ)]/(n!)

(Subscript[𝒢, 1]^MB) = 
      grand[(Subscript[𝒵, 1]^MB)] // PowerExpand

\(\(n\_MB\)\&_ = 
      nOccupy[𝒢\_1\%MB]\  // Simplify\)

\(muTonMB = {\(Solve[
            n\&_ == \(n\_MB\)\&_, μ]\)[[1, 
          1]], \ β -> 
          1\/\(k\_B\ T\)}\)

(Subscript[𝒢, 1]^MB)  /. muTonMB

Probability distribution for occupation numbers

The probability that an orbital with energy ε contains n particles is given by the grand canonical distribution where the statistical weight is determined by the Gibbs prescription.  The Boltzmann factor for the grand canonical ensemble is determined by the energy relative to the chemical potential, where the temperature and chemical potential are established by the reservoir (here composed of particles in other orbitals).  Hence, it is instructive to express that probability in terms of the mean occupancy, \(TraditionalForm\`n\&_\), instead of the chemical potential.

(Subscript[P, MB])[n_, ε, β, μ]  =  
      Exp[-β n (ε - μ)]/(n! (Subscript[𝒵, 1]^MB))

\(P\_MB[n_, n\&_]\  = \ 
      P\_MB[n, ε, β, μ] /. muTonMB // 
        Simplify\)

Thus, in classical statistics the occupation number is governed by the Poisson distribution based upon a mean occupation number \(TraditionalForm\`n\&_\) given by a Boltzmann factor with the energy expressed relative to the chemical potential.  The chemical potential depends upon density and temperature, as derived above.  The fluctuation in occupation numbers is given by the variance, which for the Poisson distribution is simply \(TraditionalForm\`n\&_\).  Thus, the relative fluctuation in occupancy is \(TraditionalForm\`1\/\@n\&_\) for the MB gas.

(Subscript[σ, MB]) = (Sqrt[nVariance])[(Subscript[𝒢, 1]^MB)] /. 
        muTonMB

\(Ϟ\_MB = σ\_MB\/\(n\_MB\)\&_ /. muTonMB\)

For later comparisons it is also useful to note that the ratio of probabilities for successive values of n is inversely proportional to n itself.

\(\(P\_MB[n + 1, n\&_]/P\_MB[n, n\&_]\  /. 
          muTonMB\)\  /. {\(\((n + 1)\)!\) -> n\ \(n!\)} // 
      Simplify\)

Relate chemical potential to density

The total particle number is obtained by integrating the mean occupation number as a function of energy against the density of states.

\(μ\_MB = \ \(\(μ\)\(\ \)\(/.\)\((\(Solve[\ 𝒩\  == \ Integrate[𝒟[ε] \(\(\ \)\(n\_MB\)\)\&_, {ε, 0, ∞}, 
                      Assumptions -> {Re[1\/\(k\_B\ T\)] > 
                            0}], \ μ]\)[[1]] /. TtoLambda // MySimplify)\)\(\ \)\)\)

MySimplify[(Subscript[μ, MB]) /. lambdaToQC]

Recognizing that the Maxwell-Boltzmann distribution can only be applied when the quantum concentration, (Subscript[n, Q]) is small, we conclude that μ should be large and negative when the MB distribution is applicable.

Examples

Here we evaluate the quantum concentrations and chemical potentials for several systems of interest.  Unless explicitly stated otherwise, SI units are employed.  Hence, when defining the physical parameters we convert to SI units and then remove the tags for the units when performing numerical calculations.

Helium gas at STP

heliumSTPValues = 
        Join[{g -> 1, 𝒩 -> Mole AvogadroConstant, 
            V -> (MolarVolume/(((Meter^3)/Mole))), 
            m -> ((4 Convert[AtomicMassUnit, Kilogram])/Kilogram)}, 
          FundamentalConstants];

λ[TT_, 
        values_] := λ /. lambdaToT /. T -> TT //. 
        values

λ[273, heliumSTPValues]

quantumConcentration = ((λ[273, heliumSTPValues]^3 𝒩)/V)  /. heliumSTPValues

((Subscript[μ, MB])/((Subscript[k, B]) T)) /. lambdaToT /. heliumSTPValues /. 
      T -> 273

Since the quantum concentration is quite small for this system, the classical approximation should be adequate.  The chemical potential is large and negative.

Liquid helium

liquidHeValues = 
        Join[{g -> 1, 𝒩 -> V density, 
            density -> (146/m), 
            m -> ((4 Convert[AtomicMassUnit, Kilogram])/Kilogram)}, 
          FundamentalConstants];

quantumConcentration = 
      N[λ[T, liquidHeValues]^3 density //. liquidHeValues  /. 
          T -> 2]

((Subscript[μ, MB])/((Subscript[k, B]) T)) /. lambdaToT //. liquidHeValues /. 
      T -> 2

Since the quantum concentration for liquid helium is quite large, the MB distribution is no longer applicable.

Bose-Einstein distribution

Systems composed of identical particles with integer spin, namely bosons, obey Bose-Einstein statistics.  In this section we derive the mean occupation number and its dispersion using the formalism of the grand canonical ensemble, as discussed in the Introduction.  We also evaluate the dependence of the chemical potential upon density and temperature and consider several representative systems.

Evaluate mean occupation number

The summation over occupation numbers extends from 0 to ∞ for bosons.

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

(Subscript[𝒢, 1]^BE) = grand[(Subscript[𝒵, 1]^BE)]

\(\(n\_BE\)\&_ = 
      nOccupy[𝒢\_1\%BE]\  // FullSimplify\)

Because \(TraditionalForm\`\(n\_BE\)\&_\) must be nonnegative, it is clear that ε - μ must also be positive.  Hence, if we choose the energy scale to place the ground state at zero energy, we conclude that the chemical potential for an ideal BE gas must be negative.  For a given density we shall find that the chemical potential is a monotonically decreasing function of temperature.  Thus, for very low temperatures the chemical potential may become so small that the occupancy of the ground state becomes so large that practically all particles are found in that one state.  This is the origin of Bose condensation.  However, under those conditions we must examine the excitation spectrum more carefully — the customary density of states formula is not sufficient when μ becomes similar to the spacing between energy levels.  The critical temperature where the chemical potential goes to zero depends upon density, and is determined below.

\(muTonBE = {\(Solve[
              n\&_ == \(n\_BE\)\&_, μ]\)[[1, 
            1]], β -> 1\/\(k\_B\ T\)}\  // 
        MySimplify\)

(Subscript[𝒢, 1]^BE)  /. muTonBE  // Simplify

Probability distribution for occupation numbers

The probability that an orbital with energy ε contains n particles is given by the grand canonical distribution as specified by the temperature and chemical potential established by the reservoir (here composed of particles in other orbitals).  It is instructive to express that probability in terms of the mean occupancy, \(TraditionalForm\`n\&_\), instead of the chemical potential.

(Subscript[P, BE])[n_, ε, β, μ]  =  
      Exp[-β n (ε - μ)]/( (Subscript[𝒵, 1]^BE))  // FullSimplify

\(P\_BE[n_, n\&_]\  = \ 
      P\_BE[n, ε, β, μ]\  /. muTonBE\  // 
        MySimplify\)

\(P\_BE[n + 1, n\&_]/P\_BE[n, n\&_] /. muTonBE\  // Simplify\)

Thus, in BE statistics the occupation number is governed by the geometric distribution for which  the ratio of probabilities for successive values of n is constant, whereas for the MB distribution that ratio is inversely proportional to n.  Hence, the permutation symmetry for bosons enhances occupancies relative to classical statistics for independent particles.  Similarly, the fluctuations are also enhanced by a factor of (Sqrt[1+ (OverBar[n])]) relative to classical (MB) statistics by the requirements of permutation symmetry for identical bosons.

(Subscript[σ, BE]) = (Sqrt[nVariance])[(Subscript[𝒢, 1]^BE)] /. muTonBE  // 
        Simplify

\(Ϟ\_BE = σ\_BE\/\(n\_BE\)\&_ /. muTonBE\  // 
        MySimplify\)

Determine critical temperature

Bose condensation occurs when the chemical potential reaches zero.  For lower temperatures it is then longer possible to accommodate all particles in excited states and it becomes necessary for a macroscopic number of particles to reside in the ground state.  The fraction of the total number of particles in the ground state increases as the temperature is reduced, reaching unity at absolute zero.  The total number of particles in excited states can be obtained by integrating the mean occupation number against the density of states.  Note that the ground state at zero energy is not counted when a continuous density of states is employed.

\(dNdx\  = \ \ \(k\_B\ T\ 𝒟[ε] \(\(\ \)\(n\_BE\)\(\ \)\)\&_ /. {μ -> 0, ε -> 
              x\ k\_B\ T}\) /. TtoLambda; 
    NexcCritical = ∫\_0\%∞ dNdx\ \[DifferentialD]x\  /. 
          lambdaToT\  // MySimplify\)

quantumConcentrationCritical = (((λ^3) NexcCritical)/(g V)) /. lambdaToT // MySimplify  // N

For a given density, the Bose condensation occurs for the critical temperature for which the quantum concentration reaches its critical value.

(Subscript[λ, c]) = ((((quantumConcentrationCritical/density)))^(1/3))

(Subscript[T, c]) = 
      N[PowerExpand[((2 π (ℏ^2))/(m (Subscript[k, B]) ((Subscript[λ, c])^2)))]]

Temperature dependence of chemical potential and occupancy

The chemical potential for fixed temperature and density must be determined by numerical solution of an integral equation.  An interpolation function is then used to represent the temperature dependence of the chemical potential and to plot the occupation number distribution for finite temperature.

Here we express energy and chemical potential in units of (Subscript[k, B]) (Subscript[T, c]) and temperature as τ = T/(Subscript[T, c]).

Nexc[
        z_] = (V/(λ^3)) ((2 Integrate[((Sqrt[ε])/(Exp[ε]/z - 1)), {ε, 0, ∞}])/(Sqrt[π]))

eq1[
        z_, τ_] := ((Nexc[z] (λ^3))/V) == ((𝒩 ((Subscript[λ, c])^3) (τ^(-3/2)))/V) /. 
        density -> (𝒩/V)

zBE[τ_] := 
      z /. FindRoot[Evaluate[eq1[z, τ]], {z, 0.7, 0.75}]

pointsNearTc = 
        Prepend[Table[{τ, τ Log[zBE[τ]]}, {τ, 1.2, 1.5, 
              0.1}], {1. , 0. }];

pointsAtLargerT = 
        Table[{τ, τ Log[zBE[τ]]}, {τ, 2, 5, 
            0.5}];

points = Join[pointsNearTc, pointsAtLargerT];

Clear[chemicalPotentialBE];
chem = Interpolation[points];
chemicalPotentialBE[τ_] := 
        chem[τ] /; τ > 1;
chemicalPotentialBE[τ_] := 
        0 /; τ ≤ 1;

\(\(DisplayTogether[
        ListPlot[points, PlotStyle -> AbsolutePointSize[4]], 
        Plot[chemicalPotentialBE[τ], {τ, 0, 5}], 
        PlotLabel -> "BE chemical potential", 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {\*""T/(Subscript[T, c])"", \*""μ/((Subscript[k, B])(Subscript[T, c]))""}];\)\)

Clear[temp];
\(\(temp[τ_] := \(n\_BE\)\&_ //. {k\_B -> 1, 
            T -> τ, μ -> 
              chemicalPotentialBE[τ]};\)\)
\(\(nBEplot[τ_] := 
        Plot[temp[τ], {ε, 0, 2}, 
          PlotRange -> {{0, 2}, {0, 5}}, 
          PlotLabel -> "BE occupation number", 
          GridLines -> None, Frame -> True, 
          FrameLabel -> {\*""ε/((Subscript[k, B])(Subscript[T, c]))"", \*""\!\(n\&_\)""}, 
          DisplayFunction -> Identity];\)\)

Show[{nBEplot[1.1],nBEplot[2],nBEplot[5]},
DisplayFunction->$DisplayFunction];

\(\(Plot3D[temp[τ], {ε, 0, 2}, {τ, 1.01, 2}, 
        PlotRange -> {{0, 2}, {1, 2}, {5, 1}}, 
        PlotLabel -> "BE Mean Occupancy", 
        AxesLabel -> {\*""ε/((Subscript[k, B])(Subscript[T, c]))"", \*""T/(Subscript[T, c])"", \*""\!\(n\&_\)""}, 
        DisplayFunction -> $DisplayFunction];\)\)

Examples

Here we evaluate the Bose critical temperature for several systems of interest.  Unless explicitly stated otherwise, SI units are employed.  Hence, when defining the physical parameters we convert to SI units and then remove the tags for the units when performing numerical calculations.

Helium-4

Liquid helium-4 undergoes a phase transition to a superfluid phase at a temperature of 2.17 kelvins at atmospheric pressure.  This phase transition is related to the Bose-Einstein condensation.

liqHe4Values = 
        Join[{g -> 1, 𝒩 -> V density, 
            density -> (146/m), 
            m -> ((4 Convert[AtomicMassUnit, Kilogram])/Kilogram)}, 
          FundamentalConstants];

(Subscript[λ, c]) //. liqHe4Values

(Subscript[T, c]) //. liqHe4Values

Although liquid helium is not really an ideal (noninteracting) gas, it is intriguing to observe that the Bose condensation temperature is close to the superfluid transition temperature.

Rubidium in ultracold atomic trap

Recently, there have been several reports of Bose-Einstein condensation in ultracold gases in low-density atomic traps.  The densities are so low that interactions are relatively unimportant and the gas is nearly ideal.  The first observation using rubidium was made by Anderson et al., Science 269 (95) 198.  Although it is difficult to determine the density of the system, the following estimate is probably reasonable.

rubidiumValues = 
        Join[{g -> 1, 𝒩 -> V density, 
            density -> (10^19), 
            m -> ((87 Convert[AtomicMassUnit, Kilogram])/Kilogram)}, 
          FundamentalConstants];

(Subscript[T, c]) //. rubidiumValues

(λ/(Subscript[λ, c])) /. lambdaToT /. T -> (50/(10^9)) //. 
      rubidiumValues

meanSpacing = 
      (density^(-1/3)) Convert[Meter, Angstrom] //. rubidiumValues // 
        N

Thus, at the density estimated for the trap the mean interatomic spacing is much larger than the typical interaction range and for a temperature of 50 nanokelvins the thermal wavelength is already larger than the critical value for Bose-Einstein condensation.  Therefore, a trap with these characteristics is capable of producing ideal Bose-Einstein condensation.

Electron pairs

Superconductivity occurs when the interaction between electrons and the lattice result in pairing of electrons with opposite spins and momenta.  These Cooper pairs are then bosons with twice the electron mass.  Although hardly an ideal Bose gas, it is nevertheless of some interest to determine the Bose critical temperature even though a proper theory of superconductivity is beyond the scope of this course.

cooperValues = 
        Join[{g -> 1, 𝒩 -> V density, 
            density -> 5 (10^28), 
            m -> ((2 ElectronMass)/Kilogram)}, 
          FundamentalConstants];

(Subscript[T, c]) //. cooperValues

A gas of Cooper pairs is clearly highly degenerate.  However, a proper treatment of the Cooper-pair model must be rather more sophisticated than this "frivolous" model.

Fermi-Dirac distribution

Systems composed of identical particles with half-integral spin, namely fermions, obey Fermi-Dirac statistics.  In this section we derive the mean occupation number and its dispersion using the formalism of the grand canonical ensemble, as discussed in the Introduction.  We also evaluate the dependence of the chemical potential upon density and temperature and consider several representative systems.

Evaluate mean occupation number

The Pauli exclusion principle for fermions limits the occupation number for each orbital to the values 0 or 1.

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

(Subscript[𝒢, 1]^FD) = grand[(Subscript[𝒵, 1]^FD)]

\(\(n\_FD\)\&_ = 
      nOccupy[𝒢\_1\%FD]\  // FullSimplify\)

\(muTonFD = \(\({\(Solve[
              n\&_ == \(n\_FD\)\&_, μ]\)[[1, 
            1]], β -> 
            1\/\(k\_B\ T\)}\)\(\ \)\(//\)\(MySimplify\)\(\ \)\)\)

(Subscript[𝒢, 1]^FD)  /. muTonFD  // MySimplify

Probability distribution for occupation numbers

The probability that a particular orbital contains n particles, where n can be only either 0 or 1 for fermions, is given by the grand canonical ensemble as follows.

(Subscript[P, FD])[n_, ε, β, μ]  =  
      Exp[-β n (ε - μ)]/( (Subscript[𝒵, 1]^FD))  // FullSimplify

\(P\_FD[n_, n\&_]\  = \ 
      P\_FD[n, ε, β, μ] /. muTonFD\  // 
        MySimplify\)

\(Table[\ P\_FD[n, n\&_], \ {n, 0, 1}\ ]\  /. muTonFD // 
      Simplify\)

Thus, in FD statistics the mean occupation number, \(TraditionalForm\`n\&_\), can be interpreted as the  probability that a state is occupied and \(TraditionalForm\`\((1 - n\&_)\)\) is then the probability that it is vacant, those being the only two possibilities for fermions.  The fluctuation in occupation numbers is suppressed by a factor of \(TraditionalForm\`\@\(1 - n\&_\)\) relative to the classical occupation probability due to the requirements of permutation symmetry.  Thus, at low temperatures these fluctuations are quite severely damped for states below the chemical potential for which the occupancy is near unity (see figure below).

(Subscript[σ, FD]) = 
      Simplify[(Sqrt[nVariance])[(Subscript[𝒢, 1]^FD)] /. muTonFD ]

\(Ϟ\_FD = σ\_FD\/\(n\_FD\)\&_ /. muTonFD\  // 
        FullSimplify\)

Evaluate Fermi energy and momentum

The Fermi momentum, (Subscript[k, F]), is defined to be the momentum of the highest occupied state when each available state is filled in order of increasing energy until all particles have been accommodated.  The corresponding energy is called the Fermi energy, (Subscript[ε, F]), and is equal to the chemical potential for T = 0.  For a completely degenerate Fermi gas, the average single-particle energy is then (3/5) (Subscript[ε, F]) and the total internal energy is simply (3/5) N (Subscript[ε, F]).  Similarly, the pressure for a nonrelativistic ideal gas is (2/3) of the energy density.  Thus, we define pDegenerate to be the pressure for a completely degenerate Fermi gas.

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

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

pDegenerate = ((2 𝒩 (Subscript[ε, F]))/(5 V))

Temperature dependence of chemical potential and occupancy

The chemical potential for fixed temperature and density must be obtained by numerical solution of an integral equation.  An interpolation function is then used to represent the temperature dependence of the chemical potential and to plot the occupation number distribution for finite temperature.

We express energy and chemical potential in units of (Subscript[k, B])(Subscript[T, F]) and temperature as τ=T/(Subscript[T, F]).

fdIntegral[μ_] = 
      1.5 Integrate[((Sqrt[ε])/(1 + Exp[ε - μ])), {ε, 0, ∞}]

muFD[τ_] := τ μ /. 
        FindRoot[
          fdIntegral[μ] == (τ^(-3/2)), {μ, 0.6, 
            0.4}]

points = 
        Prepend[Table[{τ, muFD[τ]}, {τ, 0.2, 2, 0.2}], {0. , 
            1. }];

chemicalPotentialFD = Interpolation[points];

\(\(DisplayTogether[
        ListPlot[points, PlotStyle -> AbsolutePointSize[4]], 
        Plot[chemicalPotentialFD[τ], {τ, 0.01, 1.99}], 
        PlotLabel -> "FD chemical potential", 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {\*""T/(Subscript[T, F])"", \*""μ/((Subscript[k, B])(Subscript[T, F]))""}];\)\)

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

Show[{nFDplot[0.01], nFDplot[0.1], nFDplot[0.5], nFDplot[1.0]}, 
        DisplayFunction -> $DisplayFunction];

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

Examples

Here we evaluate the Fermi temperature for several systems of interest.  Unless explicitly stated otherwise, SI units are employed.  Hence, when defining the physical parameters we convert to SI units and then remove the tags for the units when performing numerical calculations.

Helium-3

Helium-3 has two electrons paired to spin-0, but nuclear spin-(1/2).  Thus, even though the chemical properties of helium-3 are practically indistinguishable from those of the more abundant helium-4 isotope, the nuclear spin requires atomic helium-3 to behave as a fermion and helium-4 as a boson.  Even though the nucleus plays no direct role in the interatomic interactions, the effect of its spin upon the permutation symmetry of the many-body system has a profound effect upon the thermodynamics of the liquid phase at low temperatures.

liqHe3Values = 
        Join[{g -> 2, 𝒩 -> V density, 
            density -> (81/m), 
            m -> ((3 Convert[AtomicMassUnit, Kilogram])/Kilogram)}, 
          FundamentalConstants];

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

At temperatures of a few kelvins, liquid helium-3 is a moderately-degenerate Fermi system.  However, since interatomic interactions are important, its behavior is not entirely ideal.

\(\(gasHe3Values = 
        Join[
		{g -> 2, 𝒩 -> 
              Mole\ AvogadroConstant, 
            V -> MolarVolume\/\(Meter\^3\/Mole\), 
            m -> \(3\ Convert[AtomicMassUnit, Kilogram]\)\/Kilogram}, 
          FundamentalConstants];\)\)

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

As a gas near STP, even helium-3 remains quite nondegenerate.  Hence, the MB distribution is applicable under these conditions.

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, 
            Tdebye -> 345, 
            m -> \(effectiveMass\ ElectronMass\)\/Kilogram}, 
          FundamentalConstants];\)\)

(Subscript[ε, F]) Convert[Joule, ElectronVolt] //. 
      copperValues

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

pDegenerate Convert[Pascal, Atmosphere] /. τ -> 0 //. 
      copperValues

(((Sqrt[((2 (Subscript[ε, F]))/m)]) Meter)/(Second SpeedOfLight)) //. copperValues

At room temperature, conduction electrons in a metal constitute a highly degenerate Fermi gas.  Since the Fermi velocity is a small fraction of light speed, it is appropriate to treat the electron gas as nonrelativistic.  However, the pressure exerted by the electron gas is enormous.  This outward pressure must be balanced by the attractive electrostatic forces binding the electrons to the lattice.  The confinement of this gas is accomplished by a nearly uniform potential rather than by walls at the boundaries of its volume.  Hence, we can hardly call the electron gas free.  Nevertheless, since the interactions with the lattice can be represented as a smooth mean field in which all of the electrons move more or less independently of each other, the noninteracting Fermi gas model is still appropriate.  The primary effect of the mean field is to shift the energy scale.  The relatively small momentum dependence of the electron-lattice interaction is responsible for the effective mass correction to the energy-momentum relationship.

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]) Convert[Joule, ElectronVolt] //. 
      nmValues

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

pDegenerate Convert[Pascal, Atmosphere] /. τ -> 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!  In this case, the confinement of the system is accomplished by the mutual interactions between the particles themselves. Fortunately, it can be shown that these interactions can be represented fairly well by a mean field approximation that makes the Fermi gas a useful model of the system.  However, the nonrelativistic approximation may not be completely adequate because the Fermi velocity is an appreciable fraction of light speed.

(((Sqrt[((2 (Subscript[ε, F]))/m)]) Meter)/(Second SpeedOfLight)) //. nmValues

Compare the FD, BE, and MB distributions

Occupancy

To compare the mean occupancies, it is useful to employ a reduced energy variable that is referred to the chemical potential such that y =(( ε-μ)/((Subscript[k, B])T)).

Simplify[Solve[
          y == ((ε - μ)/((Subscript[k, B]) T)), μ][[1]]]

quantumConcentration = 
      MySimplify[(((λ^3) density)/g) /. Solve[μ == ((Subscript[μ, MB])[[1]] /. {N -> 
                          V density}), 
                density][[1]] /. 
          lambdaToT]

\(\(plotNvsY[n_, ymin_, ymax_, style_] := 
        Plot[Evaluate[n //. {TtoY}], {y, ymin, ymax}, 
          PlotStyle -> style, 
          PlotLabel -> "Mean Occupation Number", 
          PlotRange -> {{\(-ymax\), ymax}, {0, 5}}, 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {"y", \*""\!\(n\&_\)""}, 
          DisplayFunction -> Identity];\)\)

\(\(arrayNvsY = {plotNvsY[\(n\_FD\)\&_, \(-5.0\), 
            5.0, {RGBColor[1, 0, 0]}], 
          plotNvsY[\(n\_BE\)\&_, 0.01, 5.0, {RGBColor[0, 1, 0]}], 
          plotNvsY[\(n\_MB\)\&_, \(-5.0\), 5.0, {RGBColor[1, 0, 1]}], 
          Graphics[{RGBColor[1, 0, 0], Text["FD", {\(-2.5\), 1.2}]}], 
          Graphics[{RGBColor[0, 1, 0], Text["BE", {0.7, 3.5}]}], 
          Graphics[{RGBColor[1, 0, 1], 
              Text["MB", {\(-1.8\), 3.5}]}]};\)\)

Show[arrayNvsY, 
        DisplayFunction -> $DisplayFunction];

The classical limit pertains when y is large.  It is then useful to distinguish two cases.

1.  For moderate single-particle energies, large y corresponds to large negative chemical potential and hence to small quantum concentration.  
2. For large single-particle energies, the classical limit can remain valid even when μ is small or even positive provided that the energies are great enough.

In either case, the mean occupancy is small.  For smaller values of y, the two quantum distribution functions diverge from the classical MB model, with the occupancies being strongly enhanced for bosons or strongly suppressed for fermions relative to the classical Boltzmann factor.

Single-orbital entropy

The entropy for ideal gases can be expressed as summations with respect to the contributions made by each single-particle state available to the system.  These single-orbital entropies are obtained by differentiating the appropriate grand potential with respect to temperature and are expressed here as functions of the mean occupancy for that orbital, as determined by its energy, the chemical potential, and the distribution appropriate to the particle type (FD, BE, MB).

(Subscript[S, 1]^MB)  =  
      entropy[(Subscript[𝒢, 1]^MB)]  /. muTonMB  // Simplify

(Subscript[S, 1]^FD)  =  entropy[(Subscript[𝒢, 1]^FD)]  /. 
          muTonFD //MySimplify

(Subscript[S, 1]^BE)  =  entropy[(Subscript[𝒢, 1]^BE)]  /. 
          muTonBE //MySimplify

\(\(plotFD = 
        Plot[Evaluate[S\_1\%FD\/k\_B], {n\&_, 0.001, 0.999}, 
          PlotStyle -> {RGBColor[1, 0, 0]}, 
          DisplayFunction -> Identity];\)\)

\(\(plotBE = 
        Plot[Evaluate[S\_1\%BE\/k\_B], {n\&_, 0.001, 10}, 
          PlotStyle -> {RGBColor[0, 1, 0]}, 
          DisplayFunction -> Identity];\)\)

\(\(plotMB = 
        Plot[Evaluate[S\_1\%MB\/k\_B], {n\&_, 0.001, 10}, 
          PlotStyle -> {RGBColor[1, 0, 1]}, 
          DisplayFunction -> Identity];\)\)

\(\(Show[{plotFD, plotBE, plotMB, 
          Graphics[{RGBColor[1, 0, 0], Text["FD", {0.5, 0.2}]}], 
          Graphics[{RGBColor[0, 1, 0], Text["BE", {2.7, 2.4}]}], 
          Graphics[{RGBColor[1, 0, 1], Text["MB", {2.7, \(-0.5\)}]}]}, 
        PlotRange -> {{0, 4}, {\(-2\), 3}}, 
        PlotLabel -> "single-orbital entropy", Frame -> True, 
        GridLines -> Automatic, 
        FrameLabel -> {"mean occupation number", \*""entropy/(Subscript[k, B])""}, 
        DisplayFunction -> $DisplayFunction];\)\)

All three distributions give similar results for small occupancy where the classical limit  applies.  However, the MB distribution clearly becomes pathological for large occupancies and should not be used when (Subscript[n, MB]) exceeds unity (or even 0.5).

\(Series[S\_1\%MB, {n\&_, 0, 2}]\)

\(Series[S\_1\%FD, {n\&_, 0, 2}]\)

\(Series[S\_1\%BE, {n\&_, 0, 2}]\)

The single-orbital entropy for the BE distribution is a monotonically increasing function of \(TraditionalForm\`n\&_\), approaching \(TraditionalForm\`k\_B\ Log[n\&_]\) for large \(TraditionalForm\`n\&_\).  This behavior is consistent with one's expectation that the entropy should grow as more particles are available to share the energy.   However, some care must be exercised when applying this result below the critical temperature where the macroscopic occupation of the ground state might suggest that the ground state dominates the entropy of the system, which is not true.  Although the ground-state contribution is then of order (Subscript[k, B]) Log[N], that contribution is still quite small compared to N (Subscript[k, B]).  Provided that a macroscopic number of particles remain in the excited states, the integrated contribution from them will still be of order N (Subscript[k, B]), even if each individual contribution is small.  Therefore, the excited states continue to dominate the entropy below the condensation temperature.  Furthermore, when most of the particles are in the ground state, the grand canonical ensemble predicts extremely large fluctuations in the particle number, which may not be consistent with the external constraints.  It is then more appropriate to apply the canonical distribution.

On the other hand, the single-orbital entropy for the FD distribution, although nonnegative, is not a monotonic function of (OverBar[n]): it grows for small (OverBar[n]), reaches a  maximum for (OverBar[n])=0.5, and then decreases.  This behavior is clearly a consequence of the Pauli exclusion principle.  For sufficiently low temperatures, the  FD occupancy  approaches unity for states below the Fermi energy and approaches zero for states above the Fermi energy.  The exclusion principle inhibits occupancy fluctuations for low-lying states and hence reduces the entropy associated with those states.  The thermodynamics of the degenerate Fermi gases is clearly dominated by states in the immediate vicinity of the Fermi energy itself, where the occupancy is 0.5 and the contribution to the entropy is maximal.  All of the action occurs at the Fermi surface.

Summary

In this section we collect some of our previous results for more direct comparison of the symbolic expressions for the three distributions.  Many of these results can be represented in a more compact form using a signature function.

signature[FD]  =  -1;
signature[MB]  =  0;
signature[BE]  =  +1;

Occupation probabilities

Probability distribution

\({{"FD", P\_FD[n, n\&_]}, {"BE", P\_BE[n, n\&_]}, {"MB", 
          P\_MB[n, n\&_]}}\  // TableForm\)

Mean occupancy

\({{"FD", \(n\_FD\)\&_}, {"BE", \(n\_BE\)\&_}, {"MB", \(n\_MB\)\&_}}\  // TableForm\)

These results can be combined into a single formula using the signature function, as follows.

\(\(n\&_[
          type_]\  = \ \((Exp[\(ε - μ\)\/\(k\_B\ T\)] - signature[type]\ )\)\^\(-1\)\ ;\)\)

\(\(\(\(n\&_[#] &\)\  /@ {FD, BE, MB}\)\(\ \)\)\)

Fluctuations

\(\({{"FD", σ\_FD}, {"BE", σ\_BE}, {"MB", σ\_MB}}\  // Simplify\) // TableForm\)

{(Subscript[σ, FD]), (Subscript[σ, BE])}/(Subscript[σ, MB])  // MySimplify

These results can be combined into a single formula using the signature function, as follows.

\(σ[type_]\  = \ 
      Sqrt[n\&_\ \((1 + n\&_\ signature[type])\)]\)

σ[#] &  /@ {FD, BE, MB}

The width of the probability distribution for occupation numbers is "normal" for MB statistics, "subnormal" for FD statistics, and "supranormal" for BE statistics.

Grand potential

\(\({{"FD", 𝒢\_1\%FD}, {"BE", 𝒢\_1\%BE}, {"MB", 𝒢\_1\%MB}}\  // 
        MySimplify\) // TableForm\)

{{"FD", (Subscript[𝒢, 1]^FD) /. 
              muTonFD}, {"BE", (Subscript[𝒢, 1]^BE) /. 
              muTonBE}, {"MB", (Subscript[𝒢, 1]^MB) /. muTonMB}}  // 
        MySimplify // TableForm

The grand potential also can be combined into a single formula using the signature function, but a limiting process is needed for the MB distribution to evaluate a 0/0 condition.

(Subscript[𝒢, 1])[type_]  =  
      (Subscript[k, B]) T Log[1  -  signature[type] Exp[-((ε - μ)/((Subscript[k, B]) T))]]/signature[type]

(Subscript[𝒢, 1])[#] &  /@ {FD, BE}

Limit[(Subscript[𝒢, 1])[type],  
      signature[type] -> signature[MB]]

Chemical potential

{{"FD", μ /. muTonFD}, {"BE", μ /. 
            muTonBE}, {"MB", μ /. muTonMB}}  // TableForm

Clearly these formulas for the chemical can be combined using the signature function.

\(muTon[
        type_]\  = \ μ -> ε + 
          k\_B\ T\ Log[n\&_\/\(1 + n\&_\ signature[type]\)]\)

Entropy

{{"FD", (Subscript[S, 1]^FD)}, {"BE", (Subscript[S, 1]^BE)}, {"MB", (Subscript[S, 1]^MB)}}  // 
      TableForm

The formulas for single-orbital entropy can be combined using the signature function and can be expressed either in terms of energy and temperature or in terms of mean occupation numbers, as follows.

(Subscript[S, 1])[type_] = entropy[(Subscript[𝒢, 1])[type]]

{(Subscript[S, 1])[FD],  (Subscript[S, 1])[BE], 
          Limit[(Subscript[S, 1])[type], signature[type] -> signature[MB]]} // 
        MySimplify // TableForm

(Subscript[S, 1])[type]  /. muTon[type]  // Simplify

{(Subscript[S, 1])[FD] /. muTon[FD], (Subscript[S, 1])[BE]  /. muTon[BE] , 
          Limit[(Subscript[S, 1])[type] /. muTon[type], 
            signature[type] -> signature[MB]] } // MySimplify  // 
      TableForm

Problems

(Subscript[S, 1])  =  -(Subscript[k, B]) (∑_n) P[n] Log[P[n]]