I have this dilemma with thermodynamics I can't wrap my head around.
Let's say we have a molecule of at least 3 atoms that well-defined, non-zero moments of inertia around each of its axes, $I_x, I_y, I_z$, let's assume for simplicity that $I_x = I_y = I_z$. Classically, the rotational energy of such molecule is $$ E_{\text{rot.}} = \frac{1}{2} I \omega_x^2 + \frac{1}{2} I \omega_y^2 + \frac{1}{2} I \omega_z^2 $$
Since each of the degrees of freedom is a quadratic expression of the dynamical variable ($\omega$), the mean energy of $N$ such molecules in a high-temperature limit is $$ \left\langle E_{\text{rot.}} \right\rangle = \frac{1}{2} N k_B T + \frac{1}{2} N k_B T + \frac{1}{2} N k_B T = \frac{3}{2} N k_B T $$ courtesy of the equipartition theorem.
Now, for quantum treatment, the quantum version of the energy is given by $$ E_{\text{rot.}} = \frac{L^2}{2 I} = \frac{\hbar^2}{2 I} \ell (\ell+1) $$
The mean energy can again be calculated, in high-temperature limit by obtaining the partition sum for a single molecule first (the factor $2 \ell + 1$ is for multiplicity in each of the states with total angular momentum $\ell$) $$ Z = \sum_{\ell = 0}^\infty (2 \ell + 1) \exp \left( - \frac{\hbar^2}{2 I k_B T} \ell (\ell + 1) \right) $$ we can approximate this by an integral, which is valid for $k_B T \gg \frac{\hbar^2}{I}$ $$ Z \approx \int \limits_0^\infty \mathrm{d} \ell \, (2 \ell + 1) \exp \left( - \frac{\hbar^2}{2 I k_B T} \ell (\ell + 1) \right) = \int \limits_0^\infty \mathrm{d} u \, \exp \left( - \frac{\hbar^2}{2 I k_B T} u \right) = \frac{\hbar^2}{2 I k_B T} $$
The mean energy then is (for $N$ such molecules) $$ \left\langle E_{\text{rot.}} \right\rangle = - \frac{\partial}{\partial \beta} \log Z^N = N k_B T $$ which is not $(3/2) N k_B T$ as we should get according to the equipartition theorem at high temperatures. To me, this seems like at some point, we ditched one degree of freedom accidentally, since the expression is a result of the equipartition theorem if the number of degrees of freedom is 2, and not 3.
Note, that I can repeat the same derivation using a one-dimensional quantum rotator, using the eigenvalues of $L_z$, $\hbar m$ $$ E_{\text{rot.}} = \frac{L_z^2}{2 I} = \frac{\hbar^2}{2 I} m^2 $$ which correctly gives me $\left\langle E_{\text{rot.}} \right\rangle = \frac{1}{2} N k_B T$.
If I wanted a 2D rotator, I'd take something like this: $$ E_{\text{rot.}} = \frac{L_x^2 + L_y^2}{2 I} = \frac{\hbar^2}{2 I} \left( \ell (\ell+1) - m^2 \right) $$ however, I don't know how to evaluate the partition sum, even in approximate form. I'd expect to get mean energy $(2/2) N k_B T$ in this case, though.
So, where did I make a mistake with the 3D rotator? Where has the one degree of freedom gone?