Bose-Einstein statistics: meaning of the chemical potential

The average number of particles at the $$\varepsilon_r$$ energy level in the Bose-Einstein statistics is given by

$$\bar n_r=\frac{1}{e^{\beta(\varepsilon_r-\mu)}-1}$$

I wonder what the significance that the chemical potential $$\mu$$ has in these statistics. I see that $$\varepsilon_r$$ must be bigger than $$\mu$$, since $$n_r$$ would diverge at $$\varepsilon_r=\mu$$ and would take negative values for $$\varepsilon_r<\mu$$. This doesn't happen in the Fermi-Dirac nor Maxwell-Boltzman statistics (plotted in red and green in the picture, next to next).

I understand that, as bosons don't follow the Pauli exclusion principle, they tend to fill first the lowest energy levels, but what is the meaning of the divergence and negative values of $$n_r$$ at $$\varepsilon_r\leq\mu$$? What is the role of the chemical potential?

From a statistical-mechanical point of view, the chemical potential $$\mu$$ is just a parameter of a probability distribution. Together with other parameters, such as the inverse temperature $$\beta$$, this parameter determines some expectation values of the probability distribution, for example the expectations of the total energy and total number of quanta. These expectations reflect the values that we macroscopically observe on average.

From this point of view the constraint $$\mu < \epsilon_r$$ is simply necessary for such expectations to have sensible values – positive energy and positive number of quanta. The mathematics behind this is explained in this answer to another question of yours.

From the point of view of equilibrium thermodynamics, the chemical potential tells us the difference in the total energy between two equilibrium states having infinitesimally different number of quanta but otherwise identical extensive quantities such as entropy, volume, and number of other kinds of quanta: $$\mu(S,V,N,\dotsc) = \frac{\partial E(S,V,N,\dotsc)}{\partial N}$$

(though usually the free enthalpy is used instead of the energy; see Astarita's reference below). A negative chemical potential (at particular values of $$(S,V,N,\dotsc)$$ or equivalently $$(T,V,N,\dotsc)$$ or any other variables by which we parametrize the equilibrium states) means that an equilibrium state $$(S,V,N+\Delta N,\dotsc)$$ with a slightly larger number of quanta $$\Delta N >0$$ has lower energy than the state $$(S,V,N,\dotsc)$$. No physical law forbids such behaviour. The full functional dependence of $$\mu$$ on $$(S,V,N,\dotsc)$$ is of course such that we cannot decrease the energy indefinitely by adding quanta. Examples of the numerical dependence on the temperature are given in a paper by Cowan:

A microscopic explanation of how negative values can come about in Bosonic systems is given in a paper by Baierlein, which I warmly recommend:

see especially section VI. I can quote the relevant passage (with additional notes in brackets):

Adding an atom [or a quantum], which may be placed virtually anywhere, surely increases the spatial part of the multiplicity [exponential of the entropy] and hence tends to increase the entropy. To maintain the entropy constant [...] requires that the momentum part of the multiplicity decrease. In turn, that means less kinetic energy, and so the inequality $$\Delta E < 0$$ holds, which implies that the chemical potential is negative for an ideal gas in the semi-classical domain.

Let me add a slightly different and more general explanation of why a decrease in energy with an increase in quanta is physically possible.

The expected macroscopic energy $$\bar{E}$$ is, statistically, an average (expectation) over several possible energy values, some higher and some lower than $$\bar{E}$$. The contribution of these possible values is determined by the probabilities $$p(n_1,n_2,\dotsc )$$ of the various microstates $$(n_1,n_2,\dotsc)$$, which differ by the number of quanta in the various modes $$r$$. These probabilities depend parametrically on the macroscopic expected values $$\bar{E}, \bar{N}$$, etc. (see again this answer).

When we add a given, small number of quanta $$\Delta N$$ to the system (or more precisely: when on average we add this number of quanta, which is what "$$\Delta N$$" means), requiring that the final entropy be the same as the initial one, the final equilibrium distribution $$p^*(n_1,n_2,\dotsc )$$ must have an entropy $$-\sum_{(n_r)} p*(n_1,n_2,\dotsc )\ln p*(n_1,n_2,\dotsc )$$ equal to the one of the initial distribution. You can calculate that the change in entropy is $$-\sum_{(n_r)} \Delta p(n_1,n_2,\dotsc )\ln p(n_1,n_2,\dotsc )$$, with the normalization constraint $$\sum_{(n_r)} \Delta p(n_1,n_2,\dotsc )=0$$.

It may happen that the new probabilities give higher probability to states having energy lower than $$\bar{E}$$, and lower probability to states with energy higher than $$\bar{E}$$. The new expected value of the energy will therefore be lower than before.

Finally it's good to emphasize that any process leading from the first equilibrium state to the second usually happens outside of equilibrium (even if the two equilibria are very close). While the process occurs there is no probability distribution $$p(n_1,n_2,\dotsc)$$, and the chemical potential is still defined but may have a completely different functional dependence, depending also on non-equilibrium quantities (such as rate of change in temperature $$\mathrm{d}T/\mathrm{d}t$$). In some cases it even becomes a tensorial quantity (Eshelby tensor).

A brilliant text to understand the thermodynamic meaning of the chemical potential, also outside of equilibrium, is

I think the question may about what chemical potential is in general, so I will provide some pointers on that.

Basically, chemical potential is a bit like pressure or temperature, only instead of governing the direction of energy movements or volume changes, it governs the direction of movement of particles. The basic thermodynamic relation here is $$dU = T dS - p dV + \mu dN$$ for a simple mechanical system such as a gas. The main implication of this general formula, as it concerns chemical potential, is that particles will move, if they can, from regions of higher chemical potential to regions of low chemical potential, because such movement will give positive entropy changes. This is similar to heat movements from regions of higher temperature to regions of lower temperature.

In the case of a gas, the particles (or field exitations if you prefer) get distributed in such a way as to maximise the entropy subject to constraints. For a gas in a box (so that it is isolated) the constraints are fixed total energy and volume and number of particles (with the exception that if one treats thermal radiation as a gas (which one does not have to do) then there is no constraint on the number of particles).

So the particles (or exitations) distribute themselves over states of energy so as to maximise the entropy, and the BE distribution is saying that one gets, in the highest entropy configurations, a distribution over energy given by the distribution function $$n_r(T,\mu,\epsilon_r)$$. The parameter $$T$$ says how quickly the occupation of energy levels falls as you go to higher energy, and has such a value that the total energy is the given amount $$U$$. The parameter $$\mu$$ meanwhile takes such a value that the sum of all the $$n_r$$ will be equal to the given number of particles $$N$$.

In the classical limit (i.e. low number of particles per quantum state) $$\mu$$ is large and negative. At low temperatures it approaches the ground state energy (for bosons) or the energy of the lowest unoccupied state (for fermions).

Now to the question itself. The divergence as $$\mu$$ approaches the ground state energy $$\epsilon_g$$ is saying that $$\mu$$ never quite reaches the ground state, so the divergence does not happen in a finite system. I think it is easier to understand this the other way around, and simply note that at low temperature, when most of the particles occupy the ground state (for bosons), then the $$\mu$$ value must be close to but not quite equal to $$\epsilon_g$$.

In summary, then, $$\mu$$ can be seen as a parameter that gives information about the distribution over states on the one hand, and it is also that thermodynamic property which has to be equal, in equilibrium, between places or systems between which particles can move freely. The fact that it plays both these roles at once is not self-evident from any one formula, but it is readily deduced from standard techniques in statistical mechanics, such as by obtaining the free energy and its derivatives. Such methods show that the $$\mu$$ in the formula for $$n_r$$ is the same $$\mu$$ as the one in the formula for $$dU$$ in terms of $$dS$$, $$dV$$ and $$dN$$.