What is the classical partition function for a system of anharmonic oscillators?

Question

This is an exam question that I'm trying to solve. Here's the text:

A system is composed of N localised, but independent one-dimensional classical oscillators. Assume that the potential energy for an oscillator contains a small anharmonic term

$$ V(x) = \frac{k_0x^2}{2} + \alpha x^4 $$

where $\alpha < x4 << kT$. Write down an expression for the Canonical partition function for this system of oscillators. Then, to first order in the parameter $\alpha$ ($\alpha > 0$), derive an expression for the internal energy and the isochoric heat capacity for this system and show that the anharmonic correction tends to reduce the energy per oscillator compared to the equipartition result of a perfectly harmonic oscillator. Explain also why letting $T \rightarrow \infty$ in the expression for the heat capacity represents an unphysical result.

---

Attempted solution

The classical canonical partition function, $Q$, in this case will be

$$ Q = \iint \limits_{-\infty}^{\infty} e^{-\left( \frac{k_0}{2}x^2 + \alpha x^4 + \frac{p^2}{2m} \right)\frac{1}{kT}} \text{d}p\text{d}x $$

Now comes the difficult part. I have tried to find the Maclaurin series of $V(x)$ around $\alpha = 0$ to first order, but this just returns the original expression (weird, huh?). Then I tried to evaluate the integrals as they stand above, but this seems to not be possible:

\begin{align} Q &= \iint \limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} e^{-\frac{\alpha x^4}{kT}} e^{-\frac{p^2}{2mkT}} \text{d}p\text{d}x \\ &= \int \limits_{-\infty}^{\infty} e^{-\frac{p^2}{2mkT}} \text{d}p \int \limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} e^{-\frac{\alpha x^4}{kT}} \text{d}x \end{align}

The first integral I can do, but not the second. However, if I use that $a < x^4 << kT$, which logically leads to that $ax^4 << kT$, then I perhaps could simplify the second integral by using that $e^{-\frac{ax^4}{kT}} \approx e^{-0} = 1$. But this totally neglects the anharmonic term, which clearly I should draw some conclusions about later on.

Instead, the trick I think is to do the correct Maclaurin expansion. Should I expand the expression for $Q$? Using

$$ Q \approx Q \rvert_{\alpha = 0} + a\frac{\partial Q}{\partial a}\rvert_{\alpha = 0} $$

results in

\begin{align} Q &\approx \iint \limits_{-\infty}^{\infty} e^{-\left( \frac{k_0x^2}{2kT} + \frac{p^2}{2mkT} \right)} \text{d}p\text{d}x + a\left[ \frac{\partial}{\partial a} \left( \iint \limits_{-\infty}^{\infty} e^{-\left( \frac{k_0x^2}{2kT} + \frac{p^2}{2mkT} \right)} \text{d}p\text{d}x \right) \right] \\ &= \int \limits_{-\infty}^{\infty} e^{-\frac{p^2}{2mkT}} \text{d}p \int \limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} \text{d}x + 0 \end{align}

Since the integrals in the second term do not depend on $a$, the whole thing reduces to $0$, and the anharmonic contribution here is lost. I expected from the start that the anharmonic term would give a contribution to the final result, but it does not when I do it like this. Anyway, continuing, the integrals that are left are easily evaluated, and we get a first order approximation to the canonical partition function

\begin{align} Q &\approx 2kT \sqrt{\frac{m}{k_0 }} \end{align}

where the square root sort of resembles the angular frequency of the oscillator (assuming that $k_0$ is the force constant).

Now, finally, we can start to actually answer the question, and find an expression for the internal energy $U$

$$ U =kT^2 \left( \frac{\partial \ln Q}{\partial T} \right)_{N,V} $$

We need the logarithm of $Q$

$$ \ln Q = \ln T + \ln \left( 2k\sqrt{\frac{m}{k_0}}\right) $$

And then get the expression for $U$

\begin{align} U &= kT^2 \frac{1}{T} \\ &= kT \end{align}

Taking the derivative of $U$ with respect to $T$, should give the isochoric heat capacity:

$$ C_V = k $$

So in this case, the isochoric heat capacity does not depend on $T$ at all, so I can't really evaluate whether the high-temperature limit is meaningful or not. My approach may have been incorrect from the start.

Answer 1

Here follows a complete mathematical derivation of the expressions for the internal energy and isochoric heat capacity. I am not sure why taking the high-temperature limit is unphysical, but maybe someone else knows.

---

The potential is

$$ V(x) = \frac{k_0x^2}{2} + \alpha x^4 $$

for one dimensional, localized oscillators. The Hamiltonian then becomes

$$ H(p,x) = \frac{k_0x^2}{2} + \alpha x^4 + \frac{p^2}{2m} $$

for each particle. Since we have a system of $N$ particles, we sum over all particles to get the full Hamiltonian for the entire system (assuming all oscillators have the same mass)

$$ H(p,x) = \sum\limits_{j=1}^{N} \frac{k_0x_j^2}{2} + \alpha x_j^4 + \frac{p_j^2}{2m} $$

The general expression for the classical canonical partition function is

\begin{align} Q_{\text{N,V,T}} &= \frac{1}{N! h^{3N}} \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty}e^{-H(x,p)/kT} \text{d}x\text{d}p \\ \end{align}

The text says that the oscillators are localized, so we should take away the N! in the expression for Q, since we are dealing with distinguishable particles. Plugging in the Hamiltonian, we then get for the partition function {dropping the subscript from now on}

\begin{align} Q = \frac{1}{h^{3N}} \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty}e^{-\left( \sum\limits_{j=1}^{N} \frac{k_0 x_j^2}{2} + \alpha x_j^4 + \frac{p_j^2}{2m} \right)\frac{1}{kT}} \text{d}x\text{d}p \end{align}

We know that an exponential where the exponent is a sum we can rewrite into a product of exponentials $e^{a+b+c} = e^ae^be^c$

\begin{align} Q &= \frac{1}{h^{3N}} \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \prod\limits_{j=1}^{N} e^{-\frac{k_0x_j^2}{2kT}} e^{-\frac{\alpha x_j^4}{kT}} e^{-\frac{p_j^2}{2mkT}} \text{d}x\text{d}p \\ &= \frac{1}{h^{3N}} \prod\limits_{j=1}^{N} \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} e^{-\frac{k_0x_j^2}{2kT}} e^{-\frac{\alpha x_j^4}{kT}} e^{-\frac{p_j^2}{2mkT}} \text{d}x\text{d}p \\ &= \frac{1}{h^{3N}} \prod\limits_{j=1}^{N} \int\limits_{-\infty}^{\infty} e^{-\frac{p_j^2}{2mkT}} \text{d}p \int\limits_{-\infty}^{\infty} e^{-\frac{k_0x_j^2}{2kT}} e^{-\frac{\alpha x_j^4}{kT}} \text{d}x \\ &= \frac{1}{h^{3N}} \left( \int\limits_{-\infty}^{\infty}e^{-\frac{p^2}{2mkT}} \text{d}p \right)^N \left( \int\limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} e^{-\frac{\alpha x^4}{kT}} \text{d}x \right)^N \end{align}

where in the last step we saw that the integrand will be the same for all values of $j$. The integral over the momenta are easily evaluated; just look up any table of definite exponential integrals to find the expression. We then get

$$ Q = \left( \frac{2mkT\pi}{h^2}\right)^{3N/2} \left( \int\limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} e^{-\frac{\alpha x^4}{kT}} \text{d}x \right)^N $$

The remaining integral is too complicated to evaluate analytically, so we perform a Taylor series to first order around $\alpha = 0$ (first order Maclaurin series in $\alpha$) for the anharmonic exponential. This is justified since, as the text explains, $a > x^4 >> kT$. Doing the expansion:

\begin{align} e^{-\frac{\alpha x^4}{kT}} &\approx e^0 + a \left( \frac{\text{d}}{\text{d}a} e^{-\frac{\alpha x^4}{kT}} \right)_{a=0} \\ &\approx 1 + a\left( e^{-\frac{\alpha x^4}{kT}} \left( -\frac{x^4}{kT} \right) \right)_{a=0} \\ &\approx 1 - \frac{ax^4}{kT} \end{align}

We now use this expression instead of the exponential in the integral. Lets just evaluate the integral for now, and then plug in for $Q$ once we have an expression for the integral. Calling the integral $I$, we get

\begin{align} I &= \int\limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} e^{-\frac{\alpha x^4}{kT}} \text{d}x \\ &= \int\limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} \left( 1 - \frac{ax^4}{kT} \right) \text{d}x \\ &= \int\limits_{-\infty}^{\infty} e^{-\frac{k_0x^2}{2kT}} \text{d}x - \frac{a}{kT} \int\limits_{-\infty}^{\infty} x^4 e^{-\frac{k_0x^2}{2kT}} \text{d}x \\ &= \left( \frac{2kT\pi}{k_0^2} \right)^{1/2} - \frac{a}{kT} \int\limits_{-\infty}^{\infty} x^4 e^{-\frac{k_0x^2}{2kT}} \text{d}x \end{align}

While this integral looks intimidating, we can look it up in standard tables. It has the following solution

$$ \int\limits_{0}^{\infty} x^{2n} e^{-cx^2} \text{d}x = \left( \frac{\pi}{c} \right)^{1/2} \frac{1\cdot 3 \cdots (2n-1)}{2^{n+1}c^n} $$

Since our limits go from $-\infty$ instead of $0$, we multiply with a factor of 2 to get the correct results for the integral. Using this, we obtain

\begin{align} I = \left( \frac{2kT\pi}{k_0^2} \right)^{1/2} - \frac{2a}{kT} \left( \frac{2kT\pi}{k_0^2} \right)^{1/2} \frac{3}{8} \left( \frac{2kT}{k_0^2} \right)^2 \end{align}

We can now cancel the 8 in the denominator with the $2^2$ and 2 in front of $a$, cancel a factor of $kT$ and then factorize by moving the square root out

\begin{align} I &= \left( \frac{2kT\pi}{k_0^2} \right)^{1/2} \left(1 - \frac{3akT}{k_0^4} \right) \end{align}

We can now insert this into our expression for the partition function:

$$ Q = \left( \frac{2mkT\pi}{h^2}\right)^{3N/2} \left( \frac{2kT\pi}{k_0^2} \right)^{N/2} \left(1 - \frac{3akT}{k_0^4} \right)^N $$

In statistical thermodynamics, all thermodynamic properties are related to the natural logarithm of the partition function, so we now take the logarithm. This is the reason for factorizing the expressions; it just makes taking the logarithm much easier.

$$ \ln Q = \frac{3N}{2} \ln \left( \frac{2mkT\pi}{h^2} \right) + \frac{N}{2}\ln\left( \frac{2kT\pi}{k_0^2} \right) + N\ln \left( 1 - \frac{3akT}{k_0^4} \right) $$

From here we can start to get expressions for the internal energy and isochoric heat capacity. First, internal energy

\begin{align} U &= kT^2 \left( \frac{\partial \ln Q}{\partial T} \right)_{N,V} \\ &= kT^2 \left( \frac{3N}{2} \frac{1}{T} + \frac{N}{2} \frac{1}{T} - N \frac{3ak}{k_0^4 - 3akT} \right) \\ &= \frac{3NkT}{2} + \frac{NkT}{2} - \frac{3Nak}{k_0^4 - 3akT} \end{align}

The isochoric heat capacity is

\begin{align} C_V &= \left( \frac{\partial U}{\partial T} \right)_{N,V} \\ &= \frac{3Nk}{2} + \frac{Nk}{2} - N\left( \frac{3ak}{k_0^4 - 3akT} \right)^2 \end{align}

The first term in $C_V$ comes from the kinetic energy, and is exactly what the equipartition theorem predicts. The second term is related to the potential energy, while the last term clearly is related to the anharmonic pertubation.

The exam question now asks why taking the high-temperature limit for $C_V$ is unphysical. This I am unable to answer. Mathematically, the limit seems to diverge, but I'm not sure what kind of answer is expected here. But at least I think the mathematical derivation is correct.

Answer 2

The first answer has one problem, you have N oscilators not N particles. So we have $$(\frac{2mkTπ}{h^{2}})^{N}$$ instead of $$(\frac{2mkTπ}{h^{2}})^{\frac{3N}{2}}$$. With this correction your internal energy becomes: $$U=NK_{b}T+\frac{3\alpha NK_{b}T}{3\alpha-\frac{k^{2}}{K_{b}T}} $$. Note that if $\alpha=0$ the energy becomes: $$U=NK_{b}T$$ which is the internal energy of N harmonic oscillators as expected.

Answer 3

In high temperatures the oscillator tends to occupy larger values of x, therefore the approximation you performed to turn the original partition function into a Gaussian integral (assuming small x) becomes invalid. You can also see this in the final result for $C_v$ which unphysically turns negative for high enough T value