In Understanding Molecular Simulation, the following Maxwell Boltzmann distribution for momentum is given:
$$\mathcal{P}(p) = \left(\frac{\beta}{2\pi m}\right)^{3/2}\text{exp}\left[\frac{-\beta p^2}{2m}\right]$$
From this the relative variance of the kinetic energy can be calculated as:
$$\frac{\sigma_{p^2}^2}{\left< \ p^2\right>^2} \equiv \frac{\left< \ p^4\right> - \left< \ p^2\right>^2}{\left< \ p^2\right>^2} = \frac{2}{3} $$
This is used to motivate the fact that the instantaneous temperature of a system in the canonical ensemble will fluctuate also, and its variance can be given as:
$$\begin{split}\frac{\sigma_{T_k}^2}{\left< \ T_k \right>^2_{NVT}} & \equiv \frac{\left< \ T_k^2\right>_{NVT} - \left< \ T_k\right>_{NVT}^2}{\left< \ T_k\right>_{NVT}^2} \\ & = \frac{N\left< \ p^4\right> + N(N-1)\left< \ p^2\right>\left< \ p^2\right>- N^2\left< \ p^2\right>^2}{N^2\left< \ p^2\right>^2} \\ & = \frac{1}{N} \frac{\left< \ p^4\right> - \left< \ p^2\right>^2}{\left< \ p^2\right>^2} = \frac{2}{3N}\end{split}$$
This seems to indicate that:
$$\left< \ T_k^2\right>_{NVT} = N\left< \ p^4\right> + N(N-1)\left< \ p^2\right>\left< \ p^2\right>$$
The instantaneous temperature is:
$$k_B T = T_k = \sum^N_i{\frac{p_i^2}{2mN_f}}$$
Where $N_f$ represents the degrees of freedom and is given as $N_f = 3N-3$. Given this expression, $T_k^2$ is the following:
$$\begin{split}T_k^2 & = \left(\sum_i^N{\frac{p_i^2}{2mN_f}}\right)^2 \\ & = \sum_i^N\sum_j^N{\frac{p_i^2p_j^2}{4m^2N_f^2}} \\ & = \frac{1}{4m^2N_f^2}\left[(p_1^4+p_2^4+... +p_N^4) + (p_1^2p_2^2 + p_1^2p_3^2 + ... + p_N^2p_{N-1}^2)\right]\end{split}$$
When taking the average does this actually simplify down to what is given above? If it does why is it okay to simply say each particles momentum can be represented by the average? If this is incorrect, what method should be used to reach the desired equation for $\left< \ T_k^2 \right>$?
