Quantum rotator and equipartition theorem 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?
 A: There are some subtleties with commutation relations that you overlooked. In short, you are unwittingly using the Hamiltonian of the 2D rotator which is why the equipartition gives you the 2D result.
Quick warning on the application of the equipartition theorem in the classical case. You can add the contributions of different degrees of freedom only when they commute (with the Poisson bracket). While you could adapt the theorem to treat non-commuting variables, the easiest way is to revert to canonical variables, which have the added benefit of not having to adapt the canonical ensemble as well.
In the 3D case, the most natural way is using Euler angles. You get the Lagrangian:
$$
L = \frac{I}{2}\left(\dot\theta^2+2\cos\theta\dot\phi\dot\psi+\dot\phi^2+\dot\psi^2\right)
$$
from which you deduce:
$$
p_\theta = I\dot\theta\\
p_\phi = I(\dot\phi+\cos\theta\dot \psi)\\
p_\psi = I(\cos\theta\dot\phi+\dot\psi)\\
H = \frac{1}{2I}\left(p_\theta^2+\frac{1}{1-\cos^2\theta}(p_\phi^2-2\cos\theta p_\phi p_\psi+p_\psi^2)\right)
$$
As you can see, the application of the equipartition theorem is not so evident. You need the generalised version: $\langle x_i\frac{\partial H}{\partial x_j} \rangle = kT\delta_{ij}$ and since in our case the hamiltonian is quadratic in moment: $H = \frac{1}{2}\left(p_\theta\frac{\partial H}{\partial p_\theta}+p_\phi\frac{\partial H}{\partial p_\phi}+p_\psi\frac{\partial H}{\partial p_\psi}\right)$, so you do indeed get your result $\langle E_{rot}\rangle =\frac{3}{2}k_BT$. Note that you can explicitly calculate $Z$ (start by the momenta then the position) and you do get the same result.
To reconcile it with the QM result, you can directly quantize from the Euler angles, ie enforce canonical commutation relations for the three degrees of freedom. This is not equivalent to your original system. Indeed, you underestimate the degeneracy of the $L^2$ eigenspace as you are assuming that $L^2,L_z=p_\psi$ form a complete set of observables. Your system is 3D, so you have an additional degree of freedom, and in fact you also have $p_\phi$. Now you have a complete set of observables. The representation is therefore not quite the one of $SO(3)$ as you are used to. This is why your formula is wrong.
For the 1D version, it's simply the case of a protocole on a ring. Using polar coordinates, it's just as you'd think, no commutation subtleties, you can indeed apply the equipartition theorem since $H = \frac{p_\phi^2}{2I}$ and the quantum and classical partition functions coincide at high temperature:
$$
Z_{QM} = \sum_{m=-\infty}^{+\infty}e^{-\beta(\hbar m)^2/2I} \\
Z_c = \frac{1}{h}\int e^{-\beta H} d\phi dp_\phi \\
= \frac{2\pi}{h}\sqrt{2\pi k_BTI}
$$
with the $2\pi$ coming from the integration of $\phi\in[0,2\pi]$.
For the 2D version, you now have a particle on a sphere. You again have the commutation problem. This time, use spherical coordinates:
$$
H = \frac{1}{2I}\left(p_\theta^2+\frac{p_\phi^2}{\sin^2\theta}\right)
$$
The general equipartition theorem is still applicable and you do get what you expected $\langle E_{rot}\rangle = k_BT$.
To see the link with QM, you can also express it in terms of $L_x^2+L_y^2+L_z^2=L^2$:
$$
H = \frac{L^2}{2I}
$$
with $L_x,L_y,L_z$ having the PB $\{L_i,L_j\} = \epsilon_{ijk}$. The quantum analogue is therefore your first system which you applied in 3D. Technically, you need to be careful when deriving the canonical ensemble since the new variables don't commute. The easiest way is to go from the previously introduced canonical variables, using the change of variables
$$
p_\theta = L\cos\alpha\\
p_\phi = L\sin\theta\sin\alpha
$$
you get:
$$
Z_c = \frac{1}{h^2}\int  e^{-\beta H}dp_\theta dp_\phi d\theta d\phi \\
=\frac{1}{h^2}\int  e^{-\beta h}L\sin\theta dLd\alpha d\theta d\phi \\
=\frac{4\pi^2}{h^2}\int e^{-\beta L^2/2I}d(L^2)
$$
which is the $\hbar\to0$ limit of the discrete quantum sum:
$$
Z_{QM} = \sum_{l=0}^\infty (2l+1)e^{-\beta \hbar^2 l(l+1)/2I}
$$
Hope this helps.
