Calculating the mean separation using a canonical ensemble I am having trouble with the following physics question.
Two atoms, each of mass $m$, interact with each other by a force that can be derived by the mutual potential energy equation given by
$$U = U_{0}\left[\left(\frac{a}{x}\right)^{12} - 2\left(\frac{a}{x}\right)^{6}\right],$$
where $a$ is a constant, and $x$ is the separation between two particles. The particles are in contact with a heat bath at temperature $T$, low enough so that $kT << U_{0}$, but the temperature is high enough for statistical mechanical techniques to be applied. I need help calculating the mean separation $x(T)$ of the particles.

I was able to recognize the above equation as a Lennard-Jones Potential equation. Furthermore, since the system is in contact with a heat bath, I am aware that we should use a canonical ensemble approach. I am having some trouble making progress on this question. Help would be very appreciated :)
I was thinking about expanding the potential energy function as a Taylor series, but I am not too sure about this approach.
 A: This is a very difficult problem in general, but the prompt says to take $kT\ll U_0$ which suggests we treat the potential like a quadratic well. This means that the system asks like a quantum harmonic oscillator which we know has $\langle x \rangle=0$. This means the atoms will spend equal times closer than equilibrium and farther than equilibrium so the average will just be the equilibrium separation.
We can find the equilibrium position $x_0$ by requiring $$\left.\frac{\partial U}{\partial x}\right|_{x_0}=0\to U_0 \left(\frac{12 a^6}{x_0^7}-\frac{12 a^{12}}{x_0^{13}}\right)=0\to  x_0=a$$
This says $\langle x \rangle=a$. However this didn't really make use of any statistical mechanics. From here you could get the partition function if you like by expanding $U(x)$ about $x=a$ gives (skipping lots of algebra) $$U(x)\approx U_0\left[\frac{36 x^2}{a^2}-1\right]\to k=\frac{72U_0}{a^2}$$
So we have 
$$H\approx \frac{p^2}{2m}+\frac{1}{2}\frac{72U_0}{a^2} x^2-U_0
$$
The extra $-U_0$ can just be absorbed into the eigenvalues and then we have a quantum harmonic oscillator as promised. This gives energy eigenvalues $\varepsilon_n=\hbar\omega(n+1/2)-U_0$ with $\omega=\sqrt{k/m}=\sqrt{72U_0/ma^2}$.
From here we can use the standard techniques to derive the partition function
$$
Z = \sum_{n=0}^\infty e^{-\beta \varepsilon_n}=\sum_{n=0}^\infty e^{-\beta \left(\hbar\omega(n+1/2)-U_0\right)}=\frac{e^{-\beta\left(\hbar\omega/2-U_0\right)}}{1-e^{-\beta\hbar\omega}}
$$
From which any quantity of interest can be derived.
If your professor really wants to see $x(T)$ where $x(T)\neq a$ i.e. having a real functional dependence on $T$, then you might need to do your Taylor series to higher order than 2 in which case you'd get anharmonic terms and be forced (I believe) to do some sort of perturbation to get the energy levels and from there the partition function. I won't go into that here but you can find resources online about how to potentially go about that.
Hope this helps!
