## Quantum Gas Integral

To explicitly include indistinguishability of the particles we use the grad canonical ensemble and treat each single particle state as a system in equilibrium with a reservoir of all the other single particle states at a temperature $T$ and chemical potential $\mu$. The number of particles and energy are given by $$N = \int_0^{\infty} g(\epsilon) \overline{n}(\epsilon) d\epsilon, $$ $$E = \int_0^{\infty} \epsilon g(\epsilon) \overline{n}(\epsilon) d\epsilon, $$ where $g(\epsilon)$ is the single particle density of states, which depends on the volume $V$ and $\overline{n}(\epsilon)$ is the mean number of particles in a single particle state of energy $\epsilon$, which depends on $T$, $\mu$ and whether the particles are fermions or bosons.

Our goal is to find thermodynamic quantities as a function of $T$, $V$, and $N$. The program numerically finds by trial and error the value of $\mu$ that solves Eq. (1) for a given $T$. Once we have $\mu$ we find the energy as a function of $T$. Reduced units are used by scaling the temperature by the Fermi temperature for fermions and the critical temperature for bosons.

**Problem: Evaluation of the chemical potential of an ideal Fermi gas**

To find $\mu(T > 0)$, we need to find the
value of $\mu$ that yields the desired mean number of particles. We have
\begin{equation}
\label{eq:noninteract/munum}
N = \!\int_0^\infty \overline{n}(\epsilon)g(\epsilon)d\epsilon = \frac{V(2m)^{3/2}}{2 \pi^2
\hbar^3}\!\!\int_0^\infty\!\frac{\epsilon^{1/2}d\epsilon}{e^{\beta(\epsilon-\mu)}+1},
\end{equation}
where
\begin{equation}
g(\epsilon)\,d\epsilon = {n_{\rm spin} V \over {4\pi^2\hbar^3}} (2m)^{3/2} \,
\epsilon^{1/2} \, d\epsilon. \label{eq:noninteract/pe1}
\end{equation}
We let $\epsilon
= x \epsilon_f $, $\mu = \mu^* \epsilon_f$, and $T^*=kT/\epsilon_f$, and rewrite
\eqref{eq:noninteract/munum} as
\begin{equation}
\rho = \frac{N}{V} = \frac{(2m)^{3/2}}{2 \pi^2 \hbar^3}
\epsilon_f^{3/2}\!\int_0^\infty\!\frac{x^{1/2}\,dx}{e^{(x-\mu^*)/T^*}+1}, \label{eq:6/rhoFermi}
\end{equation}
or
\begin{equation}
\label{eq:noninteract/mufermicalc}
1 = \frac32\!\int_0^\infty\!\frac{x^{1/2}\,dx}{e^{(x-\mu^*)/T^*}+1},
\end{equation}
where we have substituted
\begin{equation}
\epsilon_f = {\hbar^2 \over 2m} (3 \pi^2 \rho)^{2/3} \qquad
\label{eq:noninteract/ef}
\end{equation}
for $\epsilon_f$. To find the dependence
of $\mu^*$ on $T^*$ use Program `QuantumGasIntegral` to evaluate the integral on the right-hand side of
\eqref{eq:noninteract/mufermicalc} numerically.

- Verify that $\mu^* = 1$ at $T^*=0$.
- Start with $T^* = 0.05$ and find $\mu^*$ such that \eqref{eq:noninteract/mufermicalc} is satisfied. Does $\mu^*$ initially increase or decrease as $T^*$ is increased from zero? Is the change of $\mu^*$ linear or quadratic in $T^*$ as $T^*$ is increased?
- Estimate the value of $T^*$ at which $\mu^* \approx 0$.
- What is the sign of $\mu^*$ for $T^* \gg 1$?
- Given the values of $\mu^*(T^*)$, the program computes the values of $E(T)$. Describe its qualitative $T$ dependence and the $T$ dependence of $C_V$.

**Problem: Evaluation of $\mu$ for an ideal Bose gas**

Program `QuantumGasIntegral` numerically evaluates the
integral on the right-hand side of
\begin{equation}
\rho = \frac{\overline{N}}{V} = \frac{(2m)^{3/2}}{4 \pi^2 \hbar^3}\!\int_0^\infty \! \!
{\epsilon^{1/2} \,d\epsilon \over e^{\beta(\epsilon-\mu)} - 1}. \label{eq:noninteract/nbose2}
\end{equation}
The goal is to
find the value of $\mu$ for a given value of $\beta$ that yields the desired
value of $\rho$.
To put \eqref{eq:noninteract/nbose2} in a convenient form we introduce dimensionless variables and
let $\epsilon = kT_c y$, $T=T_cT^*$, and $\mu = kT_c \mu^*$ and rewrite \eqref{eq:noninteract/nbose2} as
\begin{equation}
1 = \frac{2}{2.612\sqrt{\pi}}\!\int_0^\infty \frac{y^{1/2} dy}{e^{(y - \mu^*)/T^*} - 1}, \label{eq:6/numereval.a}
\end{equation}
\begin{equation}
1 =0.432\!\int_0^\infty \frac{y^{1/2} dy}{e^{(y - \mu^*)/T^*} - 1},\label{eq:6/numereval.b}
\end{equation}
where we have used
\begin{equation}
kT_c= \frac{1}{2.612^{2/3}}
\frac{2\pi \hbar^2}{m}\rho^{2/3}.\label{eq:noninteract/tcbose}
\end{equation}

- Fill in the missing steps and derive \eqref{eq:6/numereval.b}.
- The program evaluates the right-hand side of \eqref{eq:6/numereval.b}. The idea is to find $\mu^*$ for a given value of $T^*$ such that the right-hand side of \eqref{eq:6/numereval.b} equals 1. Begin with $T^*=10$. First choose $\mu^* = -10$ and find the value of the integral. Do you have to increase or decrease the value of $\mu^*$ to make the numerical value of the right-hand side of \eqref{eq:6/numereval.b} closer to 1? Change $\mu^*$ by trial and error until you find the desired result. You should find that $\mu^* \approx -25.2$ for $T^*=10$.
- Let $T^*=5$ and find the value of $\mu^*$ so that the right-hand side of \eqref{eq:6/numereval.b} equals 1.
Does $\mu$ increase or decrease in magnitude? Continue this procedure for lower values of $T^*$.
Click on the
`plot`button each time an approximately correct value of $\mu$ is found and generate a plot of $\mu^*$ versus $T^*$. - Discuss the qualitative temperature dependence of $\mu$ for fixed density.

## Resources

Problems 6.36 for fermions and 6.40 for bosons in*Statistical and Thermal Physics: With Computer Applications*, 2nd ed., Harvey Gould and Jan Tobochnik, Princeton University Press (2021).