Temperature Fluctuations in the Canonicial Ensemble

Before considering how to fiddle with particle velocities or forces to enforce constant temperature, it is worth considering what statistical thermodynamics has to say about temperature. When we have direct knowledge of instantaneous particle velocities, we know that the kinetic energy is

$$\mathscr{K} = \sum_{i=1}^{N} \sum_{\alpha\in\left\{x,y,z\right\}}\frac{p_{i,\alpha}^2}{2m} \equiv \sum_{j=1}^{3N}\frac{p_j^2}{2m}$$ (184)

where we recognize that all momentum components are independent variables. We also know that temperature (in reduced units) is directly proportional to kinetic energy:

$$\frac{3}{2}NT = \mathscr{K}$$ (185)

With this equivalence, we can consider the “instantaneous” temperature as

$$T = \frac{2}{3N}\sum_{j=1}^{3N}\frac{p_j^2}{2m}$$ (186)

Since momenta must fluctuate it is necessarily the case that instantaneous temperature also fluctuates. Let's see how much.

First, since all momenta are independent, it follows from the definition of the canonical partition function that a particle momentum component $p_j$ follows the Maxwell-Boltzmann distribution:

$$\rho(p_j) = \left(\frac{\beta}{2\pi m}\right)^{\frac32}\exp\left(-\frac{\beta p_j^2}{2m}\right)$$ (187)

We can characterize fluctuations in $p_j^2$ by dividing its variance $\sigma_{p_j^2}$ to the square of its average $\langle p_j^2\rangle^2$. The variance is defined

$$\sigma_{p_j^2} = \langle p_j^4\rangle-\langle p_j^2\rangle$$ (188)

Using the Maxwell-Boltzmann distribution:

$$\langle p_j^2\rangle = \int_{-\infty}^{\infty}dp_j p_j^2 \exp\left(-\frac{\beta p_j^2}{2m}\right) = \frac{3m}{\beta}$$ (189)

and

$$\langle p_j^4\rangle = \int_{-\infty}^{\infty}dp_j p_j^4 \exp\left(-\frac{\beta p_j^4}{2m}\right) = 15\left(\frac{m}{\beta}\right)^2$$ (190)

Thus,

$$\frac{\langle p_j^4\rangle-\langle p_j^2\rangle}{\langle p_j^2\rangle} = \frac{2}{3}.$$ (191)

Now, let's compute fluctuations in the instantaneous temperature. First, the average:

$$\langle T\rangle = \frac{2}{3N}\sum_{j=1}^{3N}\frac{\langle p_j^2... ... \frac{2}{3N}\frac{3N}{2m}\langle p_j^2\rangle = \frac{\langle p_j^2\rangle}{m}$$ (192)

Now the variance, $\langle T^2\rangle - \langle T\rangle$, starting with $\langle T^2\rangle$:

$$\langle T^2\rangle = \left(\frac{2}{3N}\right)^2\Bigg<\left(\sum_{j=1}^{3N}\frac{p_j^2}{2m}\right)\left(\sum_{k=1}^{3N}\frac{p_k^2}{2m}\right)\Bigg>$$ (193)

$$= \left(\frac{2}{3N}\right)^2 \sum_{jk}\frac{\langle p_j^2p_k^2\rangle}{(2m)^2}$$ (194)

$$= \left(\frac{2}{3N}\right)^2\left[9N\frac{\langle p_j^4\rangle}{(2m)^2}+9N(N-1)\frac{\langle p_j^2\rangle^2}{(2m)^2}\right]$$ (195)

$$= \frac{1}{(mN)^2}\left[N\langle p_j^4\rangle+N(N-1)\langle p_j^2\rangle^2\right)]$$ (196)

Note that in going from Eq. 195 to 196, we note the fact that

$$\langle p_j^2p_k^2 \rangle = \langle p_j^2\rangle \langle p_k^2\rangle = \langle p_j^2\rangle^2$$ (197)

since momenta are not correlated to each other. Putting these together:

$$\frac{\langle T^2\rangle - \langle T\rangle^2}{\langle T\rangle^2} = \frac{\frac{1}{(mN)^2}\left[N\langle p_j^4\rangle+N(N-1)\langle... ...angle p_j^2\rangle}{m}\right)^2}{\left(\frac{\langle p_j^2\rangle}{m}\right)^2}$$ (198)

$$= \frac{N\langle p_j^4\rangle+N(N-1)\langle p_j^2\rangle^2 - N^2\langle p_j^2\rangle^2}{N^2\langle p_j^2\rangle^2}$$ (199)

$$= \frac{1}{N}\frac{\langle p_j^4\rangle-\langle p_j^2\rangle^2}{\langle p_j^2\rangle^2} = \frac{2}{3N}$$ (200)

So clearly temperature fluctuates in the canonical ensemble. Of course, in the thermodynamic limit ( $N\rightarrow\infty$), these fluctuations vanish and we perceive a “constant” temperature, but in a simulation in which we resolve the momenta of a set of $N$ particles, we must observe that $T$ fluctuations as shown above. We can use this fact to decide whether or not a temperature-control scheme in MD is actually resulting in sampling the canonical ensemble.