Understanding what the Bose-Einstein distribution I'm currently studying Kittel's Solid State Physics and in his chapter on Phonon heat capacity, we need to first calculate the total energy $U$. Phonons have energy $E_n = (n+1/2)\hbar\omega$ and he first calculated the average energy $\langle E\rangle$ and using the Boltzmann factor, he showed:
$$\langle E\rangle = \dfrac{1}{2}\hbar\omega+\hbar\omega\dfrac{1}{e^{\hbar\omega/k_BT}-1}=\dfrac{1}{2}\hbar\omega+\hbar\omega\langle n \rangle$$ so then we must have
$$\langle n \rangle=\dfrac{1}{e^{\hbar\omega/k_BT}-1}$$
I recognise this as the Bose-Einstein but I'm surprised to see this as being interpreted as an average of the number of states. I always thought this was a probabilistic distribution and in fact, Kittel does seem to use this as a probability since he later writes:
$$U=\int d\omega\ \hbar\omega D(\omega)\langle n \rangle$$ where $D(\omega)$ is the density of state. In this expression $D$ already accounts for the number of photons so $\langle n \rangle$ must be some probability weight? I'm sure something is flawed in my understanding so any help is appreciated!
 A: Distribution does not automatically means a probabilistic distribution - rather, in case of of Fermi-Dirac and Bose-Einstein statistics we speak about distributions over energy/frequency. E.g., one can easily see that they are not normalizes, if integrated by energy.
Fermi-Dirac distribution can be interpreted as a probability that the state of given energy is occupied or empty, but this interpretation is problematic with Bose-Einstein distribution, which may take values greater than 1.
A: I wouldn't say that $D\left(\omega\right)$ accounts for the number of phonons; I'd say it accounts for the number of vibrational (phonon) modes at a given energy. A phonon is a unit of energy in a vibrational mode. Then, $\left<n\right>$ is the average number of phonons in a given vibrational mode (what Samuel Weir said).
So the integrand is the energy of a phonon $\hbar \omega$ times the average number of those phonons in a mode $\left<n\right>$ times the density of modes at that energy $D\left(\omega\right)$.
A: For a given oscillation frequency $\omega$, the average energy of this oscillation mode is:
$$
E_\omega = \frac{1}{Z} \sum_{n=0}^\infty \left\{ \left(n+\frac{1}{2}\right)\hbar \omega \, e^{-\beta \left(n+\frac{1}{2}\right) \hbar \omega} \right\} \tag{1}
$$
where $Z$ is the probability normalization (also known as the partition function)
$$
  Z = \sum_{n=0}^\infty \left\{  e^{-\beta \left(n+\frac{1}{2}\right) \hbar \omega} \right\} = e^{-\beta\frac{1}{2}\hbar \omega}\sum_{n=0}^\infty \left\{  e^{-\beta n \hbar \omega} \right\}=e^{-\beta\frac{1}{2}\hbar \omega} \frac{1}{1-e^{-\beta\hbar \omega}}=\frac{1}{e^{\beta\hbar \omega/2}-e^{-\beta\hbar \omega/2}} \tag{2}
$$
Using Eq.(2) to express Eq.(1) as
\begin{align*}
E_\omega &= -\frac{1}{Z} \frac{\partial}{\partial \beta} \sum_{n=0}^\infty \left\{  e^{-\beta \left(n+\frac{1}{2}\right) \hbar \omega} \right\} = -\frac{1}{Z} \frac{\partial Z}{\partial \beta} =  -\frac{\partial \ln Z}{\partial \beta}\\
&=  -\frac{\partial }{\partial \beta} \left\{ -\ln \left(e^{\beta\hbar \omega/2}-e^{-\beta\hbar \omega/2} \right)\right\} \\
&=\frac{e^{\beta\hbar \omega/2}+e^{-\beta\hbar \omega/2}}{e^{\beta\hbar \omega/2}-e^{-\beta\hbar \omega/2}}\frac{1}{2}\hbar \omega =  \coth\left( \frac{\beta \hbar \omega}{2}  \right)\frac{1}{2}\hbar \omega
\end{align*}
Therefore, for a fixed $\omega$, the average oscillation mode of phonon $\langle n \rangle$:
$$
\langle n \rangle = \frac{E_\omega}{\hbar\omega} = \frac{1}{2} \coth\left( \frac{\beta \hbar \omega}{2}  \right).
$$
Then we consider all different oscillation frequencies
$$
   U = \int_0^{\omega_D} E_\omega D(\omega) d\omega \tag{3}
$$
where $D(\omega)$ is the density of state as the number of oscillation modes between $\omega$ and $\omega + d\omega$.
Seemingly, that the Debye's model will be  adopted, known as the acoustic mode.$^1$ It assumes a linear dispersion
$$
 \omega = v_s  k.
$$
where $v_s$ is the sound speed.
Therefore, the 1-d density if state
$$
D(\omega) = \frac{L}{2\pi}\frac{dk}{d\omega}= \frac{L}{2\pi v_s}
$$
$L$ is the length of the system.
In Eq. (3), there is one more parameter, $\omega_D$, the Debye frequency, defined as the upper limit of the cut-off frequency for rendering the total number of mode being equals to the total number of oscillator (atoms or number of cells).
$$
N = \int_0^{\omega_D} D(\omega) d\omega =  \int_0^{\omega_D} \frac{L}{2\pi v_s} d\omega =  \frac{L}{2\pi v_s} \omega_D.
$$
Therefore
$$
\omega_D = 2\pi v_s \frac{N}{L}.
$$
*1. The other well known model is the Einstein model, which is easier
$$
D(\omega) = \delta(\omega - \omega_o).
$$
also know as the optical mode of phonon.
