In most basic statistical physics/condensed matter discussion the density of states is used to convert a discrete sum to a continuous integral
$$\sum_{\alpha} \mapsto \int d\epsilon \ g(\epsilon).$$
For example, a common derivation of the density of states can be done by considering the volume of a sphere in $k$-space (corresponding to the total number of states with wavevectors $<k$) and then taking the derivative e.g. Wikipedia – k-space topologies.
However under what conditions is this accurate?
Indeed, the continuous density of states method breaks down for Bose-Einstein condensates because the density of states
$$
g(\epsilon) = \frac{\sigma}{4 \pi^2} \left( \frac{2m}{\hbar^2} \right)^{3/2} \sqrt{\epsilon}
$$
gives no weight to the states at $\epsilon = 0 $ since $g(\epsilon) \propto \sqrt{\epsilon}$.
For example, jacob1729's answer here describes that
One option for this is of course $g(\epsilon)=\sum_{\alpha}\delta(\epsilon-\epsilon_\alpha)$ but in practice this might be smeared out a little whilst not hurting the approximation (eg convolve with a Gaussian of width $\Delta\ll kT$).
This is an approximation. We expect this approximation to become exact as the level spacing goes to zero under the assumption that the states do not become macroscopically occupied.
My questions
- How can we justify/quantify the statement "under the assumption that the states do not become macroscopically occupied" and what are the mathematical conditions for the convergence of the integral?
- In the derivation using $k$-space topologies, where is the assumption that the density of states cannot be small made? In the case of the Bose-Einstein condensates, if this density of states is wrong at $\epsilon = 0$ why does it work at all for larger energies?
The question How to prove that sum converges to integral using density of states? appears relevant but I still don't understand why the density of states fails for low energy states (for example, as in Bose-Einstein condensates).
Best Answer
The energy levels of bound states in a quantum system are always discrete, and so expressions like $g(\epsilon) \propto \sqrt{\epsilon}$ are always an approximation. How well such an approximation will work depends on what you're trying to calculate, and what approximation you are using. It is hard to say anything more general than that, other than a general word of caution that extra care should be taken if the function you are trying to evaluate has singularities (like the Bose-Einstein distribution in the limit $\epsilon \rightarrow \mu$).
We are so accustomed to seeing these continuous approximations for the density of states that one easily forgets that they are (1) not exact and (2) not unique.
An example might be helpful. Let us take the system you described of a 3D gas of free non-interacting particles confined to a box of dimensions $L\times L \times L$. We know that we will have eigenenergies $$\epsilon_k = \frac{\hbar^2}{2m} (k_x^2 + k_y^2 + k_z^2),$$ and that the allowed $k$ values will be discretized by the boundary conditions $k_\alpha L = 2\pi n$, where $n\in\mathbb{Z}$ and $\alpha \in \{x,y,z\}$.
As an example, assume that we want to calculate the counting function $N(\epsilon)$ which counts the total number of states below an energy $\epsilon$. We can calculate this as $$N(\epsilon) = \sum_k \theta(\epsilon - \epsilon_k) = \int_0^\epsilon d\epsilon' g(\epsilon') \theta(\epsilon - \epsilon'),$$ where $g(\epsilon) = \sum_k \delta(\epsilon - \epsilon')$ is the (exact) density of states.
Let us do this numerically to test various continuous approximations to the DOS. I'll set $2m = 1$, $\hbar = 1$ and $L = \frac{1}{2\pi}$ to make all the annoying numerical constants drop out. Then we have $$\epsilon_k = k_x^2 + k_y^2 + k_z^2 = n_x^2 + n_y^2 + n_z^2.$$ The normal approach to approximating the density of states as a continuous function would result in $g(\epsilon) \approx 2\pi \sqrt{\epsilon}$ which would then result in $N(\epsilon) \approx \frac{4\pi}{3} \epsilon^{3/2}$. We could also make another approximation by replacing the delta-functions in the exact DOS by narrow Lorentzians $\delta(\epsilon) \approx \frac{\eta}{\pi(\epsilon^2 + \eta^2)}$ for some small $\eta$.
In the three figures below I've calculated $N(\epsilon)$ using the exact DOS and the two approximation schemes outlined above (for the smeared deltas, I set $\eta = 0.02$). The only difference between the figures is the range shown on the x- and y-axes.
From the figures it is quite clear that the ordinary continuous approximation becomes very good at large energies, when there are many contributing states, but is not so great for low energies when only a smaller number of states contributes to the sum. From this we can infer that for this system, we can safely use $g(\epsilon) \approx 2\pi \sqrt{\epsilon}$ if the main contribution to the sum comes from the multitude of states at high energies. However if a single state by itself contributes significantly to the sum, it is necessary to have a very good approximation to $g(\epsilon)$ for that particular state. If you have a macroscopically occupied state in your system, it is likely to contribute significantly by itself to any function that you are trying to calculate, and in that case the approximation you are using must take that into account.