On the Chemical Potential of Ideal Fermi and Bose Gases

Knowledge of the chemical potential is essential in application of the Fermi–Dirac and the Bose–Einstein distribution functions for the calculation of properties of quantum gases. We give expressions for the chemical potential of ideal Fermi and Bose gases in 1, 2 and 3 dimensions in terms of inverse polylogarithm functions. We provide Mathematica functions for these chemical potentials together with low- and high-temperature series expansions. In the 3d Bose case we give also expansions about TB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{{{{\mathrm {B}}}}}$$\end{document}. The Mathematica routines for the series allow calculation to arbitrary order.


Introduction
The properties of ideal Fermi and Bose gases are the starting points for the understanding of the low-temperature behaviour of a broad range of physical systems, including electrons in metals [1,2], the helium liquids [3,4] and systems of trapped gases [5]. Associated with these is interest in systems of lower dimensionality including graphene [6], helium films [7] and ultra-cold atoms in quasi-1d and quasi-2d traps [8].
The properties of ideal quantum gases are expressed, conveniently and succinctly in terms of the so-called Fermi-Dirac and Bose-Einstein functions. And essential to this is knowledge of the chemical potential [9].
Traditionally Fermi-Dirac and Bose-Einstein functions have been found from tables and power series expansions. The 1938 paper by McDougall and Stoner [10] gave extensive tables for fermions and there was a discussion of the corresponding functions for bosons by London [11], with series expansions given by Robinson [12] and generalized by Clunie [13]. A consolidated treatment of these functions is found in Pathria [14].
The problem with the McDougal-Stoner tables and the above treatments was that the functions were given for different values of μ/kT which made for difficulties finding the chemical potential μ. In 1974 Betts [15] published tables of (reduced) chemical potential as a function of (reduced) temperature for 3d fermions and bosons; at that time, it was a heroic achievement of minicomputer programming. Ebner and Fu [16] gave useful series expressions for 3d fermions and analytic expressions for 2d fermions. Also they produced extensive tables [17], sadly now unavailable.
A further provision of series expansions for fermions was by Hore and Frankel [18] who gave expressions for the intermediate quantum region. These series were expanded about z = 1 where z is the fugacity. Expansions in the intermediate region for both bosons and fermions have been given by Sotnikov et al. [19].
The advent of symbolic mathematics software such as Mathematica has revolutionized the way Fermi-Dirac and Bose-Einstein calculations may be carried out. Such software can perform series expansions, reversion of series, symbolic integration, numerical integration to arbitrary precision and general symbolic manipulation. Moreover Mathematica "knows of" a large number of the "special functions" and their properties. Also it has inverse function support.
The Fermi-Dirac and the Bose-Einstein integrals may both be expressed in terms of the polylogarithm functions, once the chemical potential is known. In this paper we provide Mathematica functions to obtain the chemical potential. From these it is then straightforward to evaluate properties of Bose and Fermi gases. For reference we include the classical (Maxwell) gas. Of course at low temperatures such a gas is un-physical; it violates the third law of thermodynamics. But at high temperatures the way that the Bose and Fermi gases deviate from the classical is instructive.

The Polylogarithm Functions
The expectation value of an extensive function of energy Q(ε) is given bȳ where the sum is taken over the single particle energy states and α is the factor which accounts for the degeneracy of the particles' spin states. This is 2 for electrons (spin S = 1 2 ); more generally, it will be 2S + 1. Here n(ε) is the Fermi-Dirac or the Bose-Einstein distribution function where a is understood to be +1 for fermions and −1 for bosons. Classical particles also can be accommodated by taking a = 0. We may refer to such particles as "maxwellons".
The distribution function involves the chemical potential μ, and so, this must be known before any later calculations. A key result of this paper is the provision of Mathematica functions for evaluation of the chemical potential for free fermion and boson gases in one, two and three dimensions.
In the thermodynamic limit Eq. (1) usually may be converted to an integral where g(ε) is the energy density of states. Since Q(ε) and g(ε) are often proportional to powers of ε we require to evaluate thermodynamic integrals of the form where z = e μ/kT is the fugacity. These integrals are related to the polylogarithm functions Li s (z) [20,21]: The Maxwell case has been added for completeness. We see that the thermodynamic integrals in the Bose and the Fermi case are thus given in terms of the polylogarithm functions: It then follows that the calculation of thermodynamic functions of ideal Fermi and Bose gases boils down to the evaluation of polylogarithm functions-once the chemical potential is known.

Common Energy Scale
When discussing fermions it is customary to use the Fermi energy as a convenient energy scale. This is the energy of the highest filled state at T = 0. Such a definition is clearly not applicable to maxwellons and to bosons. But it would be convenient to have a common energy scale applicable to particles whatever their statistics. In 3d the Fermi energy ε F and the Fermi wave vector k F = (2mε F ) 1/2 / are given by The Fermi wave vector is a measure of the inverse particle spacing. And this measure is appropriate in the discussion of maxwellons and bosons. In this spirit we shall introduce a characteristic wave vector, which we call the quantum wave vector k q , and relate this to a quantum energy ε q = 2 k 2 q /2m. This would be the Fermi energy in the case of fermions. In 1, 2 and 3d this is We note that Betts [15] uses different definitions for Bose and Fermi characteristic energies; this can be confusing.

Density of States
The (energy) density of states for free particles in one, two and three dimensions is [22] In terms of the quantum characteristic energy these may be expressed (10) or, generally, in d dimensions We note that for atoms trapped in a harmonic trap the density of states is proportional to ε d−2 ; then, polylogarithm functions of different orders (usually integer) will be needed in the thermodynamic limit.

Chemical Potential and Fugacity
The number of particles in the system is given by However, the density of states g(ε) is proportional to N so that N cancels and Eq. (12) leads to an expression for the quantum energy This gives the quantum energy in terms of the chemical potential and the temperature. By inverting this relation we may find the chemical potential as a function of the quantum energy and the temperature. For the Maxwell, Fermi and Bose cases Eq. (13) gives In terms of the quantum energy we shall introduce the reduced chemical potential μ * , defined as μ * = μ/ε q ; (15) and this will be expressed as a function of τ , the reduced temperature 1 The fugacity is expressed in terms of τ as z(τ ) = e μ * (τ )/τ . Using our reduced variables, Eq. (14) become Then, from inversion of these equations, z(τ ) is given by and then μ * (τ ) by Mathematica routines for these functions (Fermi and Bose) are given in Appendix A.

Functions and Their Series Expressions
Details of the calculation of the low-temperature and high-temperature series are given in the appendices. Associated Mathematica Notebooks are included in Supplementary Information, allowing evaluation of these series to arbitrary order.
A key point about the Mathematica calculations is that while Mathematica provides the InverseFunction command used in the chemical potential and fugacity expressions, the InverseFunction provision does not extend to symbolic series calculations. For this reason we have to perform the series expansion first and then obtain the inverse by reversing the power series.
In this section below we summarize the results of these calculations. We shall use the Mathematica functions of Appendix A to create plots of the chemical potential and the series approximations given. A Mathematica notebook creating the plots below is given in Supplementary Information MOESM1_ESM.nb.

Fermions in 1d
In Eqs. (22) and (19) we put d = 1 giving Chemical potential-low-temperature series Chemical potential-high-temperature series The first term is the Maxwell chemical potential. The 1d Fermi chemical potential, together with low-T and high-T approximations, is shown in Fig. 1. Fugacity-low-temperature series. Since z F1 (τ ) diverges as τ → 0 there is no simple low-temperature power series. But the low-temperature behaviour can be expressed as with the limiting low-temperature behaviour Fugacity-high-temperature series The first term is the Maxwell fugacity.
Chemical potential-high-temperature series The first term is the Maxwell chemical potential. The 1d Bose chemical potential, together with low-T and high-T approximations, is shown in Fig. 2.
The first term is the Maxwell fugacity.
Chemical potential-low-temperature series There is no series in ascending powers of τ , but the following series of exponentials follows from Eq. (43) at low temperatures.
Chemical potential-high-temperature series corresponding to Eq. (2) of Ebner and Fu. The first term of Eq. (46) is the Maxwell chemical potential. The 2d Fermi chemical potential, together with low-T and high-T approximations, is shown in Fig. 3.
Fugacity-low temperatures. Since z F2 (τ ) diverges as τ → 0 there is no simple low-temperature power series. But in the spirit of Eq. (45) and writing z F2 (τ ) as we may regard the two terms in the brackets as a terminating low-temperature expansion, with the divergent giving the limiting low-temperature behaviour.
Fugacity-high-temperature series The first term is the Maxwell fugacity.

Bosons in 2d
In Eqs. (23) and (20) we put d = 2 giving But since Li 1 (z) = − ln(1 − z), it follows that in 2d explicit expressions may be obtained. Thus: In 2d we have the special results manifestations of May's theorem on Fermi-Bose correspondence in 2d [23]. We note that May's theorem holds only in the thermodynamic limit; it breaks down for finite systems [24].
Chemical potential-low-temperature series There is no series in ascending powers of τ , but the following series of exponentials follows from Eq. (52) at low temperatures.
Chemical potential-high-temperature series (57) The first term is the Maxwell chemical potential. The 2d Bose chemical potential, together with low-T and high-T approximations, is shown in Fig. 4.

Fugacity-low temperatures
There is no low-temperature series for the fugacity, but we may regard the expression for z B2 (τ ), Eq. (53), as a terminating low-temperature expansion: with low-temperature limit z B2 (0) = 1.
Fugacity-high-temperature series The first term is the Maxwell fugacity.

Fermions in 3d
In Eqs. (22) and (19) we put d = 3 giving Chemical potential-low-temperature series Chemical potential-high-temperature series This equation corresponds to that given by Ebner and Fu after their Eq. (3). But their equation has a typo. The first term of Eq. (65) is the Maxwell chemical potential.
The 3d Fermi chemical potential, together with low-T and high-T approximations, is shown in Fig. 5.
Fugacity-low-temperature series Since z F3 diverges as τ → 0 there is no simple low-temperature power series. But the low-temperature behaviour can be expressed as with the limiting low-temperature behaviour Fugacity-high-temperature series The first term is the Maxwell fugacity.

Bosons in 3d
In three dimensions bosons can undergo Bose-Einstein condensation. When this happens the chemical potential will be zero and the fugacity will be unity. We denote the reduced Bose temperature by τ B . In Eqs. (23) and (20) we put d = 3 to give μ * and z when τ > τ B . Then the reduced chemical potential and the fugacity are given by The Bose temperature τ B is the zero of μ * (τ ) of Eq. (69), that is, Chemical potential-low-temperature series An expression equivalent to the first term of the series is given by London [11]; the author thanks Bill Mullin for drawing his attention to this.
Chemical potential-high-temperature series The first term is the Maxwell chemical potential.
The 3d Bose chemical potential, together with low-T and high-T approximations, is shown in Fig. 6.

Fugacity-low-temperature series
Fugacity-high-temperature series The first term is the Maxwell fugacity.

Chemical Potential Plots in 1, 2 and 3d
We now use the Mathematica functions to create plots of the chemical potential in one-, two-and three dimensions. These will contrast the differences and similarities of the Fermi, Bose and Maxwell cases.

Three Dimensions
The temperature dependence of the 3d chemical potential is shown in Fig. 7. At low temperatures the fermion μ * (τ ) → 1 as τ → 0 (μ(T ) → ε F as T → 0). Upon cooling, bosons in 3d undergo BEC at the Bose temperature τ B where the chemical potential goes to zero. Just above the Bose temperature μ * (τ ) increases quadratically in (τ − τ B ), while below the Bose temperature μ * (τ ) is identically zero. There is a discontinuity in the second derivative of μ * (τ ) at τ = τ B . We shall in the following sections see that there is no BEC for d < 3; this is an example of the Mermin-Wagner theorem [25]. Then there is macroscopic occupation only at T = 0.
The maxwellon chemical potential increases from zero as the temperature increases from zero, with a nonzero slope; this is a violation of the third law of thermodynamics, which requires ∂μ/∂T → 0 as T → 0.

Two Dimensions
The temperature dependence of the 2d chemical potential is shown in Fig. 8.
but they never actually get there. In 2d we have the May's theorem result and at high temperatures, the Maxwell μ is midway between those of the Fermi and the Bose cases. In other words, in the high-temperature limit the Fermi chemical potential is ε q /2 higher than the Maxwell value and the Bose chemical potential is ε q /2 below. At low temperatures the Bose chemical potential goes to zero. There is no BEC, so μ is zero only at zero temperature. But the low-temperature μ is "very flat" (see Fig. 8 right plot); it is as if the 2d bosons are "trying to" condense. We shall see in Sect. 5.3 that there is certainly no BEC in 1d. But here, in 2d, we might say that BEC is "marginal".
The low-temperature Fermi chemical potential is similar to the Bose case, just shifted up by ε q .
The Maxwell chemical potential increases from zero as the temperature increases from zero with nonzero slope, in violation of the third law of thermodynamics.

One Dimension
The temperature dependence of the 1d chemical potential is shown in Fig. 9.
At high temperatures the Fermi and Bose chemical potentials, Eq. (91/92), are going in the same direction as the Maxwell case: the Fermi above and the Bose below but as the temperature increases the Fermi and Bose curves move further away from the Maxwell. The fractional deviation (μ * F,B (τ ) − μ * M (τ ))/μ * M (τ ) does go to zero as τ → ∞.
The low-temperature Bose chemical potential goes to zero as τ 2 . There is not even a hint of BEC; the macroscopic occupation of the ground state occurs only at T = 0.
The low-temperature Fermi chemical potential has an interesting form. As the temperature increases from zero μ increases a little before turning over and decreasing. This is not a violation of the third law since ∂μ/∂T → 0 as T → 0.
The Maxwell chemical potential increases from zero with nonzero slope as the temperature increases from zero, in violation of the third law of thermodynamics.
with Mathematica implementation

A.1.9 Fermi Chemical Potential in 3d
This is given by Eq. (62): with Mathematica implementation

A.2 Implementation Notes
The calculations for d = 1 and d = 3 can be relatively time-consuming as these involve the evaluation of inverse functions. However, in the d = 2 case the evaluations will be much quicker. In order to make comparisons I timed the evaluation and plotting of each function over the range 0 ≤ τ ≤ 2. The timings, in seconds, are given in the For fermions in 1 and 3d we see it is approximately ten times faster to evaluate the chemical potential than the fugacity. In these cases, if the fugacity is required, it is better to evaluate the chemical potential and obtain the fugacity from z = e −μ * /τ . For bosons in 1 and 3d it is faster to evaluate the fugacity. Indeed numerical evaluation of μ * B (τ ) can give (erroneous) small imaginary contributions because of the accumulation of round-off errors. The number in the brackets for the timing of the 1d Bose chemical potential corresponds to the time for evaluating and plotting the "real" part of μ * .

B High-Temperature Expansions
The high-temperature expansions are able to accommodate the Fermi, Bose and Maxwell cases together, indicating how they tend to the same value in the T → ∞ limit. Moreover, in what follows we shall treat the dimensionality d as a variable. It is instructive to keep d general at this stage and then specify it at the end. Alternatively d could be set now, whereupon the following expressions are much simplified.

B.1 Quantum Energy
We start from Eq. (13) for the quantum energy: in the high-temperature case, it is preferable to work, initially, in terms of the fugacity z.
Upon changing the variable of integration to the dimensionless x = ε/kT we have Here τ is the reduced temperature. This is the starting point for the Mathematica series calculation.

B.2 Series Expansion
At high temperatures n(ε) is small, and so, z is small. We expand the bracket factor in powers of z in terms which Eq. (81) becomes the series The integrals are straightforward so that

B.3 Fugacity
Equation (85) is a series for 1/τ d/2 in powers of z. We invert this to give a series for z in powers of 1/τ d/2 : In the Maxwell case (a = 0) we find

B.4 Chemical Potential
The Maxwell (reduced) chemical potential μ * corresponding to Eq. (21). In expressing the high-temperature power series for Fermi and Bose gases it is convenient to write the chemical potential as the Maxwell chemical potential plus the Fermi/Bose correction. Thus we shall examine the high-temperature expansion of μ * (τ ) − μ * M (τ ). Now The The logarithm of this series multiplied by τ then gives the μ series The d = 1, d = 2 and d = 3 cases of this equation, for a = 1 and a = −1 (Fermi and Bose) correspond to the high-temperature chemical potential series, Eqs. (29), (36), (46), (57), (65), (73).

B.5 Mathematica Calculation
This Mathematica calculation is given in the Notebook MOESM2_ESM.nb in Supplementary Information. These are the steps.
To perform the Mathematica calculation of the high-temperature expansions we start from inverse temperature integral, Eq. (81): 1. Expand the denominator in the integrand in powers of the small z.
2. The integration is then done term by term on the expansion. 3. This gives τ −d/2 as a power series in z. 4. This series is "reversed" to give z as a series in (decreasing) powers of τ d/2 .

C Fermi Low-Temperature (Sommerfeld) Expansions
The low-temperature Fermi expansions procedure was pioneered by Sommerfeld [26]. We shall follow the intuitive treatment of Reif [27], based on the observation that the derivative of the Fermi-Dirac distribution is sharply peaked at low temperatures. As in the high-temperature case we shall keep d general at this stage and then specify it at the end. Alternatively d could be set now, whereupon the following expressions are much simplified.

C.1 Fermi Energy
We start from Eq. (13) for the quantum energy: Since here we are restricted to fermions only we shall revert to the terminology ε F . We integrate by parts: Here n (ε) is the derivative of the Fermi-Dirac distribution function with respect to energy and we note that the first term vanishes so long as d > 0. Now n (ε) is sharply peaked about ε = μ. We transform the variable of integration to x = (ε − μ)/kT : and then, only small values of x will contribute to the integral. This means that "negligible error" will be introduced by extending the lower limit of the integral to −∞.
Thus we are making an approximation for the Fermi energy integral: We call this a Sommerfeld integral.

Upon dividing both sides of this equation by (kT ) d/2 we have an equation in reduced variables
This is an expression for τ in terms of τ/μ * . We shall invert this to obtain τ/μ * in terms of τ , and from this, obtain μ * in terms of τ . Equation (97) is the starting point for the Mathematica series calculation.

C.2 Sommerfeld Expansion
The Sommerfeld expansion then follows by expanding the rightmost bracket of Eq. (97) in powers of x and integrating term by term.
These integrals are pure numbers. Denote The d = 1 and d = 3 cases of this correspond to Eqs. (28) and (64). We note immediately a problem for d = 2; in this case, there is no power series in τ (see Sect. C.5).

C.5 Errors in the Sommerfeld Expansion
Landau and Lifshitz [28], in treating the 3d Fermi gas, note that the error introduced in the Sommerfeld procedure, whereby the lower limit of the integral is extended to −∞ involves neglecting "exponentially small" terms. These terms are negligibly small so long as the Sommerfeld series does not terminate. However, for even d the series does terminate, and then, these exponential terms can dominate. Indeed when d = 2 we have the Sommerfeld approximation μ * = 1; there are no powers of τ . Fortunately in 2d we have the exact result, expressed as the series of Eq. (45). This shows the "neglected" exponential terms. There is a further discussion of the errors in Sachs [2], p. 165.

C.6 Mathematica Calculation
This Mathematica calculation is given in the Notebook MOESM3_ESM.nb in Supplementary Information. These are the steps.
The (reduced) Bose temperature τ B is the z = 1 value of this expression We note that for all 0 < z < 1 (the allowed range for z), Eq. (115) will give values of τ only greater or equal to τ B . In other words the inversion of this equation will apply only when τ ≥ τ B .

E.1 Series Expansion
We are interested in the behaviour at small τ . For τ < τ B we have z = 1 (macroscopic number of particles in the ground state) and then as τ increases τ B , z decreases from 1 towards zero. So at low temperatures 1 − z is a small quantity so we shall (or at least then, doing the multiplications by τ − τ B and by τ B , adding and collecting the τ − τ B powers give (122) corresponding to Eq. (72).
The Mathematica calculation of these results is given in the Notebook MOESM5_ESM.nb in Supplementary Information. It follows the steps outlined above.