D.58 The particle energy distributions

This note derives the Maxwell-Boltzmann, Fermi-Dirac, and Bose-Einstein energy distributions of weakly inter­acting particles for a system for which the net energy is precisely known.

The objective is to find the shelf numbers $\vec{I}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $(I_1,I_2,I_3,\ldots)$ for which the number of eigen­functions $Q_{\vec{I}}$ is maximal. Actually, it is mathemati­cally easier to find the maximum of $\ln(Q_{\vec{I}})$, and that is the same thing: if $Q_{\vec{I}}$ is as big as it can be, then so is $\ln(Q_{\vec{I}})$. The advantage of working with $\ln(Q_{\vec{I}})$ is that it simplifies all the products in the expressions for the $Q_{\vec{I}}$ derived in derivation {D.57} into sums: mathematics says that $\ln(ab)$ equals $\ln(a)$ plus $\ln(b)$ for any (positive) $a$ and $b$.

It will be assumed, following derivation {N.24}, that if the maximum value is found among all shelf occupation numbers, whole numbers or not, it suffices. More daringly, errors less than a particle are not going to be taken seriously.

In finding the maximum of $\ln(Q_{\vec{I}})$, the shelf numbers cannot be completely arbitrary; they are constrained by the conditions that the sum of the shelf numbers must equal the total number of particles $I$, and that the particle energies must sum together to the given total energy $E$:

$$\sum_s I_s = I \qquad \sum_s I_s {\vphantom' E}^{\rm p}_s = E.$$

Mathema­ticians call this a constrained maxim­ization problem.

According to calculus, without the constraints, you can just put the derivatives of $\ln(Q_{\vec{I}})$ with respect to all the shelf numbers $I_s$ to zero to find the maximum. With the constraints, you have to add “penalty terms” that correct for any going out of bounds, {D.48}, and the correct function whose derivatives must be zero is

$$F = \ln(Q_{\vec I}) - \epsilon_1\bigg(\sum_s I_s - I\big... ...epsilon_2\bigg(\sum_s I_s {\vphantom' E}^{\rm p}_s - E\bigg)$$

where the constants $\epsilon_1$ and $\epsilon_2$ are unknown penalty factors called the Lagrangian multi­pliers.

At the shelf numbers for which the number of eigen­functions is largest, the derivatives $\partial{F}$$\raisebox{.5pt}{$/$}$$\partial{I}_s$ must be zero. However, that condition is difficult to apply exactly, because the expressions for $Q_{\vec{I}}$ as given in the text involve the factorial function, or rather, the gamma function. The gamma function does not have a simple derivative. Here typical textbooks will flip out the Stirling approxi­mation of the factorial, but this approxi­mation is simply incorrect in parts of the range of interest, and where it applies, the error is unknown.

It is a much better idea to approximate the differen­tial quotient by a difference quotient, as in

$$\begin{eqnarray*} \lefteqn{0 =\frac{\partial F}{\partial I_s} \approx \frac{... ... F(I_1,I_2,\ldots,I_{s-1},I_s,I_{s+1},\ldots)}{I_s + 1 - I_s}. \end{eqnarray*}$$

This approxi­mation is very minor, since according to the so-called mean value theorem of mathematics, the location where $\Delta{F}$$\raisebox{.5pt}{$/$}$$\Delta{I}_s$ is zero is at most one particle away from the desired location where $\partial{F}$$\raisebox{.5pt}{$/$}$$\partial{I}_s$ is zero. Better still, $I_s+\frac12$ $\vphantom0\raisebox{1.5pt}{$\equiv$}$ $I_{s,{\rm {best}}}$ will be no more that half a particle off, and the analysis already had to commit itself to ignoring fractional parts of particles anyway. The difference quotient leads to simple formulae because the gamma function satisfies the condition $(n+1)!$ $\vphantom0\raisebox{1.5pt}{$=$}$ $(n+1)\,n!$ for any value of $n$, compare the notations section under “!”.

Now consider first distin­guishable particles. The function $F$ to differen­tiate is defined above, and plugging in the expression for $Q^{\rm {d}}_{\vec{I}}$ as found in derivation {D.57} produces

$$F = \ln(I!) + \sum_s\left[I_s\ln(N_s)-\ln(I_s!)\right] -... ...psilon_2 \bigg(\sum_s I_s {\vphantom' E}^{\rm p}_s - E\bigg)$$

For any value of the shelf number $s$, in the limit $I_s\downarrow-1$, $F$ tends to negative infinity because $I_s!$ tends to positive infinity in that limit and its logarithm appears with a minus sign. In the limit $I_s\uparrow+\infty$, $F$ tends once more to negative infinity, since $\ln(I_s!)$ for large values of $I_s$ is according to the so-called Stirling formula approximately equal to $I_s\ln(I_s)-I_s$, so the $-\ln(I_s!)$ term in $F$ goes to minus infinity more strongly than the terms propor­tional to $I_s$ might go to plus infinity. If $F$ tends to minus infinity at both ends of the range $\vphantom0\raisebox{1.5pt}{$-$}$1 $\raisebox{.3pt}{$<$}$ $I_s$ $\raisebox{.3pt}{$<$}$ $\infty$, there must be a maximum value of $F$ somewhere within that range where the derivative with respect to $I_s$ is zero. More specifi­cally, working out the difference quotient:

$$\frac{\Delta F}{\Delta I_s} = \ln(N_s) - \ln(I_s+1) - \epsilon_1 -\epsilon_2 {\vphantom' E}^{\rm p}_s = 0$$

and $-\ln(I_s+1)$ is infinity at $I_s$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\vphantom0\raisebox{1.5pt}{$-$}$1 and minus infinity at $I_s$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\infty$. Somewhere in between, $\Delta{F}$$\raisebox{.5pt}{$/$}$$\Delta{I}_s$ will cross zero. In particular, combining the logarithms and then taking an exponential, the best estimate for the shelf occupation number is

$$I_{s,{\rm best}}= I_s + {\textstyle\frac{1}{2}} = \frac{... ...vphantom' E}^{\rm p}_s+\epsilon_1}} - {\textstyle\frac{1}{2}}$$

The correct­ness of the final half particle is clearly doubtful within the made approxi­mations. In fact, it is best ignored since it only makes a difference at high energies where the number of particles per shelf becomes small, and surely, the correct proba­bility of finding a particle must go to zero at infinite energies, not to minus half a particle! Therefore, the best estimate $\iota^{\rm {d}}$ $\vphantom0\raisebox{1.5pt}{$\equiv$}$ $I_{s,{\rm {best}}}$$\raisebox{.5pt}{$/$}$$N_s$ for the number of particles per single-particle energy state becomes the Maxwell-Boltzmann distribution. Note that the derivation might be off by a particle for the lower energy shelves. But there are a lot of particles in a macroscopic system, so it is no big deal.

The case of identical fermions is next. The function to differen­tiate is now

$$\begin{eqnarray*} F & = & \sum_s\left[\ln(N_s!)-\ln(I_s!)-\ln((N_s-I_s)!)\righ... ... \epsilon_2\bigg(\sum_s I_s {\vphantom' E}^{\rm p}_s - E\bigg) \end{eqnarray*}$$

This time $F$ is minus infinity when a shelf number reaches $I_s$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\vphantom0\raisebox{1.5pt}{$-$}$1 or $I_s$ $\vphantom0\raisebox{1.5pt}{$=$}$ $N_s+1$. So there must be a maximum to $F$ when $I_s$ varies between those limits. The difference quotient approxi­mation produces

$$\frac{\Delta F}{\Delta I_s} = -\ln(I_s+1) + \ln(N_s-I_s) - \epsilon_1 -\epsilon_2 {\vphantom' E}^{\rm p}_s = 0$$

which can be solved to give

$$I_{s,{\rm best}}= I_s + {\textstyle\frac{1}{2}} = \frac{... ...ilon_1}}{1+e^{\epsilon_2{\vphantom' E}^{\rm p}_s+\epsilon_1}}$$

The final term, less than half a particle, is again best left away, to ensure that 0 $\raisebox{-.3pt}{$\leqslant$}$ $I_{s,{\rm {best}}}$ $\raisebox{-.3pt}{$\leqslant$}$ $N_s$ as it should. That gives the Fermi-Dirac distribution.

Finally, the case of identical bosons, is, once more, the tricky one. The function to differen­tiate is now

$$\begin{eqnarray*} F & = & \sum_s\left[\ln((I_s+N_s-1)!)-\ln(I_s!)-\ln((N_s-1)!... ... \epsilon_2\bigg(\sum_s I_s {\vphantom' E}^{\rm p}_s - E\bigg) \end{eqnarray*}$$

For now, assume that $N_s$ $\raisebox{.3pt}{$>$}$ 1 for all shelves. Then $F$ is again minus infinity for $I_s$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\vphantom0\raisebox{1.5pt}{$-$}$1. For $I_s\uparrow\infty$, however, $F$ will behave like $-(\epsilon_1+\epsilon_2{\vphantom' E}^{\rm p}_s)I_s$. This tends to minus infinity if $\epsilon_1+\epsilon_2{\vphantom' E}^{\rm p}_s$ is positive, so for now assume it is. Then the difference quotient approxi­mation produces

$$\frac{\Delta F}{\Delta I_s} = \ln(I_s+N_s) - \ln(I_s+1) - \epsilon_1 -\epsilon_2 {\vphantom' E}^{\rm p}_s = 0$$

which can be solved to give

$$I_{s,{\rm best}}= I_s + {\textstyle\frac{1}{2}} = \frac{... ...antom' E}^{\rm p}_s+\epsilon_1}-1} - {\textstyle\frac{1}{2}}.$$

The final half particle is again best ignored to get the number of particles to become zero at large energies. Then, if it is assumed that the number $N_s$ of single-particle states on the shelves is large, the Bose-Einstein distribution is obtained. If $N_s$ is not large, the number of particles could be less than the predicted one by up to a factor 2, and if $N_s$ is one, the entire story comes part. And so it does if $\epsilon_1+\epsilon_2{\vphantom' E}^{\rm p}_s$ is not positive.

Before addressing these nasty problems, first the physical meaning of the Lagrangian multi­plier $\epsilon_2$ needs to be established. It can be inferred from examining the case that two different systems, call them $A$ and $B$, are in thermal contact. Since the inter­actions are assumed weak, the eigen­functions of the combined system are the products of those of the separate systems. That means that the number of eigen­functions of the combined system $Q_{\vec{I}_A\vec{I}_B}$ is the product of those of the individual systems. Therefore the function to differen­tiate becomes

$$\begin{eqnarray*} F & = &\ln(Q_{\vec I_A}Q_{\vec I_B}) \\ && - \epsilon_{1,... ...} + \sum_{s_B} I_{s_B} {\vphantom' E}^{\rm p}_{s_B}- E \bigg) \end{eqnarray*}$$

Note the constraints: the number of particles in system $A$ must be the correct number $I_A$ of particles in that system, and similar for system $B$. However, since the systems are in thermal contact, they can exchange energy through the weak inter­actions and there is no longer a constraint on the energy of the individual systems. Only the combined energy must equal the given total. That means the two systems share the same Lagrangian variable $\epsilon_2$. For the rest, the equations for the two systems are just like if they were not in thermal contact, because the logarithm in $F$ separates, and then the differen­tiations with respect to the shelf numbers $I_{s_A}$ and $I_{s_B}$ give the same results as before.

It follows that two systems that have the same value of $\epsilon_2$ can be brought into thermal contact and nothing happens, macroscopi­cally. However, if two systems with different values of $\epsilon_2$ are brought into contact, the systems will adjust, and energy will transfer between them, until the two $\epsilon_2$ values have become equal. That means that $\epsilon_2$ is a temperature variable. From here on, the temperature will be defined as $T$ = 1/${\epsilon_2}k_{\rm B}$, so that $\epsilon_2$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1/${k_{\rm B}}T$, with $k_{\rm B}$ the Boltzmann constant. The same way, for now the chemical potential $\mu$ will simply be defined to be the constant $-\epsilon_1$$\raisebox{.5pt}{$/$}$$\epsilon_2$. Chapter 11.14.4 will eventually establish that the temperature defined here is the ideal gas temperature, while derivation {D.62} will establish that $\mu$ is the Gibbs free energy per atom that is normally defined as the chemical potential.

Returning now to the nasty problems of the distribution for bosons, first assume that every shelf has at least two states, and that $({\vphantom' E}^{\rm p}_s-\mu)$$\raisebox{.5pt}{$/$}$${k_{\rm B}}T$ is positive even for the ground state. In that case there is no problem with the derived solution. However, Bose-Einstein condens­ation will occur when either the number density is increased by putting more particles in the system, or the temperature is decreased. Increasing particle density is associated with increasing chemical potential $\mu$ because

$$I_s = \frac{N_s-1}{e^{({\vphantom' E}^{\rm p}_s-\mu)/{k_{\rm B}}T}-1}$$

implies that every shelf particle number increases when $\mu$ increases. Decreasing temperature by itself decreases the number of particles, and to compensate and keep the number of particles the same, $\mu$ must then once again increase. When $\mu$ gets very close to the ground state energy, the exponential in the expression for the number of particles on the ground state shelf $s$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1 becomes very close to one, making the total denominator very close to zero, so the number of particles $I_1$ in the ground state blows up. When it becomes a finite fraction of the total number of particles $I$ even when $I$ is macroscopi­cally large, Bose-Einstein condens­ation is said to have occurred.

Note that under reasonable assumptions, it will only be the ground state shelf that ever acquires a finite fraction of the particles. For, assume the contrary, that shelf 2 also holds a finite fraction of the particles. Using Taylor series expansion of the exponential for small values of its argument, the shelf occupation numbers are

$$\begin{eqnarray*} && I_1 = \frac{(N_1-1){k_{\rm B}}T}{{\vphantom' E}^{\rm p}_1... ...hantom' E}^{\rm p}_3-{\vphantom' E}^{\rm p}_2)} \\ && \vdots \end{eqnarray*}$$

For $I_2$ to also be a finite fraction of the total number of particles, ${\vphantom' E}^{\rm p}_2-{\vphantom' E}^{\rm p}_1$ must be similarly small as ${\vphantom' E}^{\rm p}_1-\mu$. But then, reasonably assuming that the energy levels are at least roughly equally spaced, and that the number of states will not decrease with energy, so must $I_3$ be a finite fraction of the total, and so on. You cannot have a large number of shelves each having a finite fraction of the particles, because there are not so many particles. More precisely, a sum roughly like $\sum_{s=2}^\infty{\rm {const}}/s\Delta{E}$, (or worse), sums to an amount that is much larger than the term for $s$ $\vphantom0\raisebox{1.5pt}{$=$}$ 2 alone. So if $I_2$ would be a finite fraction of $I$, then the sum would be much larger than $I$.

What happens during condens­ation is that $\mu$ becomes much closer to ${\vphantom' E}^{\rm p}_1$ than ${\vphantom' E}^{\rm p}_1$ is to the next energy level ${\vphantom' E}^{\rm p}_2$, and only the ground state shelf ends up with a finite fraction of the particles. The remainder is spread out so much that the shelf numbers immediately above the ground state only contain a negligible fraction of the particles. It also follows that for all shelves except the ground state one, $\mu$ may be approximated as being ${\vphantom' E}^{\rm p}_1$. (Specific data for particles in a box is given in chapter 11.14.1. The entire story may of course need to be modified in the presence of confinement, compare chapter 6.12.)

The other problem with the analysis of the occupation numbers for bosons is that the number of single-particle states on the shelves had to be at least two. There is no reason why a system of weakly-inter­acting spinless bosons could not have a unique single-particle ground state. And combining the ground state with the next one on a single shelf is surely not an acceptable approxi­mation in the presence of potential Bose-Einstein condens­ation. Fortunately, the mathematics still partly works:

$$\frac{\Delta F}{\Delta I_1} = \ln(I_1+1) - \ln(I_1+1) - \epsilon_1 -\epsilon_2 {\vphantom' E}^{\rm p}_1 = 0$$

implies that $\epsilon_1-\epsilon_2{\vphantom' E}^{\rm p}_1$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0. In other words, $\mu$ is equal to the ground state energy ${\vphantom' E}^{\rm p}_1$ exactly, rather than just extremely closely as above.

That then is the condensed state. Without a chemical potential that can be adjusted, for any given temperature the states above the ground state contain a number of particles that is completely unrelated to the actual number of particles that is present. Whatever is left can be dumped into the ground state, since there is no constraint on $I_1$.

Condens­ation stops when the number of particles in the states above the ground state wants to become larger than the actual number of particles present. Now the mathematics changes, because nature says “Wait a minute, there is no such thing as a negative number of particles in the ground state!” Nature now adds the constraint that $I_1$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0 rather than negative. That adds another penalty term, $\epsilon_3I_1$ to $F$ and $\epsilon_3$ takes care of satisfying the equation for the ground state shelf number. It is a sad story, really: below the condens­ation temperature, the ground state was awash in particles, above it, it has zero. None.

A system of weakly inter­acting helium atoms, spinless bosons, would have a unique single-particle ground state like this. Since below the condens­ation temperature, the elevated energy states have no clue about an impending lack of particles actually present, physical properties such as the specific heat stay analytical until condens­ation ends.

It may be noted that above the condens­ation temperature it is only the most probable set of the occupation numbers that have exactly zero particles in the unique ground state. The expec­tation value of the number in the ground state will include neighboring sets of occupation numbers to the most probable one, and the number has nowhere to go but up, compare {D.62}.