In a system with $N$ particles in some volume $V$ in contact with a reservoir of temperature $T$, we find that $$\bar{n_i}=\frac{g_i}{e^{\frac{{\epsilon}_i -\mu}{kT}} \pm 1}$$
depending on whether the particles are fermions or bosons. In the case of low $T$, we get $\bar{n_i} \approx g_{i} e^{\frac{\mu -{\epsilon}_{i}}{kT}}$ .
In some cosmology textbooks I've looked at, they state a formula for the number density of some species as $d_i=g_i e^{\frac{\mu}{kT}}\int{e^{\frac{{-\epsilon _i}}{kT}}d^3 p}$, which is essentially integrating over the approximation above. If I am not mistaken, this refers to the number of particles of species i (not energy level) per unit volume. There's two things I am confused about: where is any reference to volume made in the equation? The BE/FD distributions should give the total number of particles in a certain energy state over the total volume, correct?
The second thing im confused about: we have $\sum \bar{n_i}=N$ where we sum over all distinct energy states. If we assumed the energy spectrum was approximately continuous, shouldn't we have something like $N=\int{\bar{n_i} dE}$, not $N=\int{\bar{n_i} d^3 p}$? I dont get why the equation integrates over momentum instead of energy--its not like $\bar{n_i}$ is the average number of particles in some region of phase space, so why integrate over momentum?
Overall it seems like the equation for $d_i$ is treating the BE/FD distributions as the number density over phase space (i.e. the expected number of particles to be in a certain volume of phase space), while the equations dont seem to describe this. Am I wrong?