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