I'm looking to find the thermodynamic (NVT) free energy of a classical coupled harmonic oscillator system such as the one below:

enter image description here

(image taken from http://openmetric.org/StatisticalPhysics/equilibrium/week3.html)

I would like a solution that allows an arbitrary $N$ number of masses, and ideally I would like to have a fully general expression with arbitrary (different) spring constants and arbitrary (different) masses.

I tried to compute this free energy by hand by computing the partition function:

$$ Z = \int_{-\infty}^{\infty}d\vec{p} \int_{-\infty}^{\infty}d\vec{x}\; e^{-\beta H}$$


$$ H = \sum_{i=1}^{N} \frac{p_i^2}{2m_i} + \sum_{i=0}^{N} m_i\omega_i^2(x_{i+1}-x_i)^2 $$

and $x_i$ denotes the displacement from equilibrium of the $i$th block, and $x_0=x_{N+1}=0$ represent the walls at the ends.

I was able to derive the expressions for the free energy $F = -\frac{1}{\beta}\ln Z$ for one, two, and three blocks with identical masses and identical springs (with the hope of seeing an extendible pattern) but sadly no obvious patterns emerged. The calculations are also rather tedious.

I don't doubt that this has already been done many times before -- does anyone have a reference to a solution?

  • $\begingroup$ I think you should be able to use the transfer matrix method to solve this. Models very similar to yours should be plentiful in the polymer physics literature. $\endgroup$
    – alarge
    Jul 11, 2015 at 3:49

2 Answers 2


Make the change of variable to $\delta_i = x_{i+1} - x_i,\; i=0\dots N-1.$ Then the system is uncoupled and $Z = \prod_i z_i$ with $$ z_i = \int e^{-\beta p^2/2m_i}dp \times \int e^{-\beta m_i\omega_i^2\delta^2}d\delta = \frac{\pi\sqrt{2}}{\beta\omega_i}.$$

  • $\begingroup$ The integrand doesn't decouple nicely as you claimed. On the exponential, there is a term involving $x_{N+1}-x_N = -x_N$, and note that $x_N = \sum_{i=0}^{N-1} \delta_i$. $\endgroup$
    – higgsss
    Aug 14, 2017 at 8:04

Let $k_i = m_i \omega_i^2$, and we have \begin{equation} H\big(\{p_i\},\{x_i\}\big) = \sum_{i=1}^{N} \frac{p_i^2}{2m_i} + \sum_{i=0}^{N} \frac{k_i}{2}(x_{i+1}-x_i)^2, \end{equation} where $x_{{\scriptscriptstyle N+1}} = x_{{\scriptscriptstyle 0}} = 0$. The partition function is given by \begin{equation} \begin{split} Z &= \int \left(\prod_{i=1}^N dp_i\right) \left(\prod_{i=1}^N dx_i\right) \,\exp\left[-\frac{H\big(\{p_i\},\{x_i\}\big)}{k_{B}T}\right]\\ &= \int \left(\prod_{i=1}^N dp_i\right)\, \exp\left(-\frac{1}{k_{B}T}\sum_{i=1}^{N} \frac{p_i^2}{2m_i}\right)\int \left(\prod_{i=1}^N dx_i\right)\, \exp\left[-\frac{1}{k_{B}T}\sum_{i=0}^{N}\frac{k_{i}}{2}(x_{i+1} - x_{i})^2\right]\\ &=I_p I_x\,, \end{split} \end{equation} where $I_p$ and $I_x$ respectively represent the $p$ and $x$ integrals in the second line.

The $p$ integral is a product of Gaussian integrals, which gives \begin{equation} I_p = (2\pi k_B T)^{N/2}\sqrt{\prod_{i=1}^{N}m_{i}} \end{equation}

To evaluate $I_x$, we first make the change of variables $x_{i}^\prime = x_{i} / \sqrt{k_B T} $. Then, $I_x = (k_B T)^{N/2} \tilde{I}_{x}$, where \begin{equation} \begin{split} \tilde{I}_{x} &= \int\left(\prod_{i=1}^N dx_i^\prime\right) \,\exp\left[-\sum_{i=0}^{N}\frac{k_{i}}{2}(x_{i+1}^\prime - x_{i}^\prime)^2\right]\\ &= \int\left(\prod_{i=1}^{N+1} dx_i^\prime\right) \,\delta(x_{N+1}^\prime)\,\exp\left[ -\sum_{i=0}^{N}\frac{k_{i}}{2}(x_{i+1}^\prime - x_{i}^\prime)^2\right]. \end{split} \end{equation}

In fact, for all practical purposes, it suffices to stop at this point. That is, \begin{equation} Z = I_p I_x = (2\pi)^{N/2} (k_B T)^N \tilde{I}_x, \end{equation} where $\tilde{I}_x$ is independent of any thermodynamic variable of the system (in this case $T$).

Still, it is possible to evaluate the integral $\tilde{I}_x$. To do this, we make another change of variables $y_{i} = x_{i}^\prime - x_{i - 1}^\prime \ \ (i = 1, 2, \ldots, N)$. The associated Jacobian determinant is equal to 1, and it follows that $x_{N+1}^\prime = \sum_{i=1}^{N+1} y_{i}$. Also, we can represent the delta function as \begin{equation} \delta(x_{N+1}^\prime) = \int \frac{dq}{2\pi}\, \exp(iqx_{N+1}^\prime) \end{equation} Then, \begin{equation} \begin{split} \tilde{I}_x &= \int\frac{dq}{2\pi}\left(\prod_{i=0}^{N} dy_{i+1}\right) \,\exp\big[iq({\scriptstyle \sum_{i=0}^{N} y_{i+1}})\big] \exp\left( -\sum_{i=0}^{N}\frac{k_{i}}{2}y_{i+1}^2\right)\\ &= \int \frac{dq}{2\pi} \prod_{i=0}^{N}\Bigg[\int dy_{i+1} \exp\bigg(-\frac{k_{i}}{2}y_{i+1}^2 + iqy_{i+1}\bigg)\Bigg]\\ &= \int \frac{dq}{2\pi} \prod_{i=0}^{N} \sqrt{\frac{2\pi}{k_i}}\exp\left(-\frac{q^2}{2k_i}\right)\\ &=\frac{(2\pi)^{N/2}}{\sqrt{\Big(\sum_{i=0}^{N}\frac{1}{k_{i}}\Big)\prod_{i=0}^{N} k_i}}. \end{split} \end{equation}

Collecting the results we have derived so far, we obtain \begin{equation} Z = I_p I_x = (2\pi k_B T)^{N} \sqrt{\frac{\prod_{i=1}^{N}m_{i}}{\Big(\sum_{i=0}^{N}\frac{1}{k_{i}}\Big)\prod_{i=0}^{N} k_i}}. \end{equation}


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.