Angular distribution of dipole moments? Please tell me if my reasoning is correct. The energy of a dipole in an electric field is given by $\vec{a}\cdot\vec{E}=-aE\cos(\theta)$. For a collection of non-interacting dipoles in a thermodynamic system, the angular distribution is just given by the Gibbs-Boltzmann factor:
$$f(\theta)=Ae^{\beta aE\cos\theta}$$
where $\beta$ is the inverse temperature and $A$ is the normalization coefficient:
$$A=\int_0^{2\pi} e^{\beta a E \cos\theta} d\theta=2\pi I_0(\beta a E)$$
The normalization constant seems too complicated/weird. Is this correct?
 A: The Hamiltonian of a rigid rotor electric dipole is given by
$$H=\underbrace{\frac{p^2}{2m}}_{\text{translation}}+\underbrace{\frac{L^2_{\theta}}{2mR^2}+\frac{L^2_{\phi}}{2mR^2\sin^2\theta}}_{\text{rotational}}-\underbrace{2qRE\cos\theta}_{\text{dipole-field interaction}}$$
We see that the $(\theta,L_{\theta})$ conjugate variables don't separate cleanly from the $(\phi,L_{\phi})$ conjugate variables in the Hamiltonian, i.e. we have to talk about the distribution along the $\theta$, $\phi$, $L_{\theta}$, and $L_{\phi}$ coordinates all together. This will actually work in our favor.
Ignoring the translational part of the partition function, we have:
$$\begin{align}
Z_{\text{rot+dip-field}}&=\frac{1}{(2\pi\hbar)^2}\int e^{-\frac{\beta L^2_{\theta}}{2mR^2}-\frac{\beta L^2_{\phi}}{2mR^2\sin^2\theta}+2\beta qRE\cos\theta}\,dL_{\theta}dL_{\phi}d\theta d\phi\\
&=\frac{1}{(2\pi\hbar)^2}2\pi\left(\int_{-\infty}^{\infty}e^{-\frac{\beta L^2_{\theta}}{2mR^2}}dL_{\theta}\right)\left(\int_{0}^{\pi}e^{2\beta qRE\cos\theta}\left(\int_{-\infty}^{\infty}e^{-\frac{\beta L^2_{\phi}}{2mR^2\sin^2\theta}}dL_{\phi}\right)d\theta\right)\\
&=\frac{1}{(2\pi\hbar)^2}2\pi\left(\sqrt{\frac{2\pi mR^2}{\beta}}\right)\left(\int_{0}^{\pi}e^{2\beta qRE\cos\theta}\left(\sqrt{\frac{2\pi mR^2\sin^2\theta}{\beta}}\right)d\theta\right)\\
&=\frac{1}{(2\pi\hbar)^2}2\pi\frac{2\pi mR^2}{\beta}\int_{-1}^{1}e^{\alpha\cos{\theta}}\,d(\cos\theta),\,\,\,(\alpha=2\beta qRE)\\
&=\frac{1}{(2\pi\hbar)^2}2\pi\frac{2\pi mR^2}{\beta}\frac{2}{\alpha}\sinh(\alpha )\\
&=\boxed{\frac{mRT^2}{\hbar^2qE}\sinh\left(\frac{2qRE}{T}\right)}
\end{align}$$
The true distribution in $\theta$ can be just read off the above calculation by stopping short of the final $\theta$ integration and dividing by the partition function.
$$\boxed{f(\theta )=\frac{\alpha}{2\sinh(\alpha )}e^{\alpha \cos(\theta )}\sin(\theta )}$$
A: No, it's not. You have to integrate in the spherical angle. Using spherical coordinates $(\theta, \varphi)$.
For example the classical canonical distribution  we have $p(\theta)= \frac{e^{\beta pE\cos(\theta)}}{Z}$, where $p=qa$ is the dipole moment.
This function $Z$ (partition function) is given by
\[
Z= \int_0^{\pi} e^{\eta\cos\theta}\sin(\theta)d\theta \int_0^{2\varphi}d\varphi= \frac{2\pi}{\eta}\int_{-\eta}^{\eta}e^udu=\frac{4\pi}{\eta}\sinh(\eta),
\]
where $\eta=\beta pE$ and we used the change of coordinates $u=\eta\cos\theta$.
For ex. we can get the average potential energy in this distribution
\[
\langle -\vec p\cdot E\rangle= -pE\langle\cos\theta\rangle,
\]
where 
\[
\langle\cos\theta\rangle=\frac{1}{Z}\int_0^{\pi}  e^{\eta\cos\theta}\cos\theta\sin\theta d\theta \int_0^{2\varphi}d\varphi=\frac{\partial}{\partial\eta}\ln Z=\coth(\eta)-1/\eta
\]
