D.63 Chemical potential in the distributions

The following convoluted derivation of the distribution functions comes fairly straightly from Baierlein [4, pp. 170-]. Let it not deter you from reading the rest of this otherwise very clearly written and engaging little book. Even a nonengineering author should be allowed one mistake.

The derivations of the Maxwell-Boltzmann, Fermi-Dirac, and Bose-Einstein distributions given previously, {D.59} and {D.60}, were based on finding the most numerous or most probable distribution. That implicitly assumes that significant deviations from the most numerous/probable distributions will be so rare that they can be ignored. This note will bypass the need for such an assumption since it will directly derive the actual expectation values of the single-particle state occupation numbers $\iota$. In particular for fermions, the derivation will be solid as a rock.

The mission is to derive the expectation number $\iota_n$ of particles in an arbitrary single-particle state $\pp{n}////$. This expectation value, as any expectation value, is given by the possible values times their probability:

$$\iota_n = \sum_q i_n P_q$$

where $i_n$ is the number of particles that system energy eigenfunction $\psi^{\rm S}_q$ has in single-particle state $\pp{n}////$, and $P_q$ the probability of the eigenfunction. Since thermal equilibrium is assumed, the canonical probability value $e^{-{\vphantom' E}^{\rm S}_q/{k_{\rm B}}T}/Z$ can be substituted for $P_q$. Then, if the energy ${\vphantom' E}^{\rm S}_q$ is written as the sum of the ones of the single particle states times the number of particles in that state, it gives:

$$\iota_n = \frac{1}{Z} \sum_q i_n e^{-(i_1{\vphantom' E}^... ...n +i_{n+1}{\vphantom' E}^{\rm p}_{n+1}+\ldots)/k_{\rm B}T}.$$

Note that $i_n$ is the occupation number of single-particle state $\pp{n}////$, just like $I_s$ was the occupation number of shelf $s$. Dealing with single-particle state occupation numbers has an advantage over dealing with shelf numbers: you do not have to figure out how many system eigenfunctions there are. For a given set of single-particle state occupation numbers $\vec\imath=\vert i_1,i_2,\ldots\rangle$, there is exactly one system energy eigenfunction. Compare figures 11.2 and 11.3: if you know how many particles there are in each single-particle state, you know everything there is to know about the eigenfunction depicted. (This does not apply to distinguishable particles, figure 11.1, because for them the numbers on the particles can still vary for given occupation numbers, but as noted in chapter 11.11, there is no such thing as identical distinguishable particles anyway.)

It has the big consequence that the sum over the eigenfunctions can be replaced by sums over all sets of occupation numbers:

$$\iota_n = \frac{1}{Z} \underbrace{\sum_{i_1}\sum_{i_2}\l... ..._n +i_{n+1}{\vphantom' E}^{\rm p}_{n+1}+\ldots)/k_{\rm B}T}$$

Each set of single-particle state occupation numbers corresponds to exactly one eigenfunction, so each eigenfunction is still counted exactly once. Of course, the occupation numbers do have to add up to the correct number of particles in the system.

Now consider first the case of $I$ identical bosons. For them the occupation numbers may have values up to a maximum of I:

$$\begin{eqnarray*} \lefteqn{\iota_n = \frac{1}{Z} \underbrace{\sum_{i_1=0}^I\... ..._n +i_{n+1}{\vphantom' E}^{\rm p}_{n+1}+\ldots)/{k_{\rm B}}T} \end{eqnarray*}$$

One simplification that is immediately evident is that all the terms that have $i_n=0$ are zero and can be ignored. Now apply a trick that only a mathematician would think of: define a new summation index $i_n'$ by setting $i_n=1+i_n'$. Then the summation over $i_n'$ can start at 0 and will run up to $I-1$. Plugging $i_n=1+i_n'$ into the sum above gives

$$\begin{eqnarray*} \lefteqn{\iota_n = \frac{1}{Z} \underbrace{\sum_{i_1=0}^I\... ...p}_n +i_{n+1}{\vphantom' E}^{\rm p}_{n+1}+\ldots)/k_{\rm B}T} \end{eqnarray*}$$

This can be simplified by taking the constant part of the exponential out of the summation. Also, the constraint in the bottom shows that the occupation numbers can no longer be any larger than $I-1$ (since the original $i_n$ is at least one), so the upper limits can be reduced to $I-1$. Finally, the prime on $i_n'$ may as well be dropped, since it is just a summation index and it does not make a difference what name you give it. So, altogether,

$$\begin{eqnarray*} \lefteqn{\iota_n = \frac{1}{Z} e^{-{\vphantom' E}^{\rm p}_n/... ...p}_n +i_{n+1}{\vphantom' E}^{\rm p}_{n+1}+\ldots)/k_{\rm B}T} \end{eqnarray*}$$

The right hand side falls apart into two sums: one for the 1 in $1+i_n$ and one for the $i_n$ in $1+i_n$. The first sum is essentially the partition function $Z^-$ for a system with $I-1$ particles instead of $I$. The second sum is essentially $Z^-$ times the expectation value $\iota_n^-$ for such a system. To be precise

$$\iota_n = \frac{1}{Z} e^{-{\vphantom' E}^{\rm p}_n/{k_{\rm B}}T} Z^- \left[1 + \iota_n^-\right]$$

This equation is exact, no approximations have been made yet.

The system with $I-1$ particles is the same in all respects to the one for $I$ particles, except that it has one less particle. In particular, the single-particle energy eigenfunctions are the same, which means the volume of the box is the same, and the expression for the canonical probability is the same, meaning that the temperature is the same.

But when the system is macroscopic, the occupation counts for $I-1$ particles must be virtually identical to those for $I$ particles. Clearly the physics should not change noticeably depending on whether 10$\POW9,{20}$ or 10$\POW9,{20}$ + 1 particles are present. If $\iota_n^-=\iota_n$, then the above equation can be solved to give:

$$\iota_n = 1\left/\left[\frac{Z}{Z^-} e^{{\vphantom' E}^{\rm p}_n/{k_{\rm B}}T} - 1 \right]\right.$$

The final formula is the Bose-Einstein distribution with

$$e^{-\mu/{k_{\rm B}}T} = \frac{Z}{Z^-}$$

Solve for $\mu$:

$$\mu = - k_{\rm B}T \ln\left(\frac{Z}{Z^-}\right) = \frac{ - k_{\rm B}T\ln(Z)+k_{\rm B}T\ln(Z^-)}{I - (I-1)}$$

The final fraction is a difference quotient approximation for the derivative of the Helmholtz free energy with respect to the number of particles. Now a single particle change is an extremely small change in the number of particles, so the difference quotient will be to very great accuracy be equal to the derivative of the Helmholtz free energy with respect to the number of particles. And as noted earlier, in the obtained expressions, volume and temperature are held constant. So, $\mu=\left(\partial{F}/\partial{I}\right)_{T,V}$, and (11.39) identified that as the chemical potential. Do note that $\mu$ is on a single-particle basis, while $\bar\mu$ was taken to be on a molar basis. The Avogadro number $I_{\rm A}$ = 6.022,1 10$\POW9,{26}$ particles per kmol converts between the two.

Now consider the case of $I$ identical fermions. Then, according to the exclusion principle, there are only two allowed possibilities for the occupation numbers: they can be zero or one:

$$\iota_n = \frac{1}{Z} \underbrace{\sum_{i_1=0}^1\ldots\s... ... +i_{n+1}{\vphantom' E}^{\rm p}_{n+1}+\ldots)/{k_{\rm B}}T}$$

Again, all terms with $i_n=0$ are zero, so you can set $i_n=1+i_n'$ and get

$$\begin{eqnarray*} \lefteqn{\iota_n = \frac{1}{Z} \underbrace{\sum_{i_1=0}^1\... ..._n +i_{n+1}{\vphantom' E}^{\rm p}_{n+1}+\ldots)/{k_{\rm B}}T} \end{eqnarray*}$$

But now there is a difference: even for a system with $I-1$ particles $i_n'$ can still have the value 1 but the upper limit is zero. Fortunately, since the above sum only sums over the single value $i_n'=0$, the factor $(1+i_n')$ can be replaced by $(1-i_n')$ without changing the answer. And then the summation can include $i_n'=1$ again, because $(1-i_n')$ is zero when $i_n'=1$. This sign change produces the sign change in the Fermi-Dirac distribution compared to the Bose-Einstein one; the rest of the analysis is the same.

Here are some additional remarks about the only approximation made, that the systems with $I$ and $I-1$ particles have the same expectation occupation numbers. For fermions, this approximation is justified to the gills, because it can be easily be seen that the obtained value for the occupation number is in between those of the systems with $I-1$ and $I$ particles. Since nobody is going to count whether a macroscopic system has 10$\POW9,{20}$ particles or 10$\POW9,{20}$ + 1, this is truly as good as any theoretical prediction can possibly get.

But for bosons, it is a bit trickier because of the possibility of condensation. Assume, reasonably, that when a particle is added, the occupation numbers will not go down. Then the derived expression overestimates both expectation occupation numbers $\iota_n$ and $\iota_n^-$. However, it could at most be wrong, (i.e. have a finite relative error) for a finite number of states, and the number of single-particle states will be large. (In the earlier derivation using shelf numbers, the actual $\iota_n$ was found to be lower than the Bose-Einstein value by a factor $(N_s-1)/N_s$ with $N_s$ the number of states on the shelf.)

If the factor $Ze^{{\vphantom' E}^{\rm p}_1/{k_{\rm B}}T}/Z^-$ is one exactly, which definitely means Bose-Einstein condensation, then $i_1=1+i_1^-$. In that case, the additional particle that the system with $I$ particles has goes with certainty into the ground state. So the ground state better be unique then; the particle cannot go into two ground states.