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 gasTo 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$.
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.
