# 2D Ising model, understanding autocorrelations, Monte Carlo

I've been struggling quite a bit with implementing an autocorrelation code into my current project. The autocorrelation as it is now, is increasing exponentially from 1 at the start of my MC run, and hitting 2 halfway through the MC simulation regardless of how many sweeps I do through the lattice.

The system

10x10 square lattice with no external magnetic field and ferromagnetic coupling. The reason for 10x10 is for fast execution of the code in order to build it.

Here is what I've done so far:

1. Letting the Metropolis Monte Carlo work until the system is in equilibrium (checking this by running two different initial states with different random seeds).
2. Then I start sweeping through the lattice, updating the energies and magnetization for each attempt at flipping one spin. When I've done one sweep over the lattice, the last value of the energy and magnetization gets stored. Then it continues to the next sweep and updating the values further. Thus, the energies and magnetization only gets stored once per sweep.
3. When I've done, say 2000 sweeps, I calculate the autocorrelation for the system according to Newman & Barkema (Eq 3.21 in http://itf.fys.kuleuven.be/~fpspXIII/material/Barkema_FPSPXIII.pdf). The formula reads: $$\chi(t)=\frac{1}{t_{max}-t}\sum_{t'=0}^{t_{max}-t}m(t')m(t'+t)-\frac{1}{t_{max}-t}\sum_{t'=0}^{t_{max}-t}m(t')\times\frac{1}{t_{max}-t}\sum_{t'=0}^{t_{max}-t}m(t'+t)$$ where t defines number of sweeps of the lattice, i.e the displacement/lag from some value.

My problem From my simulation at $k_{b}T/J$, where $k_{b}$ is set to 1, $J=1$ is the ferromagnetic coupling, and $T=1$ the autocorrelation function grows. I've tried to normalize it by dividing by the first value resulting in a start at 1, but it acts strange as stated. Thus, I started to calculate by hand trying to see if I had coded something wrong. I worked with a system where each spin in the lattice had spin 1 for every $t\Rightarrow m(t')=m(t'+t)=\langle m \rangle=1$. I then cutoff the autocorrelation at $t$=1000 when running 2000 sweeps. The formula then reduces down to $$\chi(t)=\frac{1}{2000-t}\sum_{t'=0}^{2000-t}1-\frac{1}{2000-t}\sum_{t'=0}^{2000-t}1\times\frac{1}{2000-t}\sum_{t'=0}^{2000-t}1$$. Then for some values: $$\chi(0)=\frac{1}{2000}\sum_{t'=0}^{2000}1-\frac{1}{2000}\sum_{t'=0}^{2000}1\times\frac{1}{2000}\sum_{t'=0}^{2000}1$$ $$=\frac{2001}{2000}-\left(\frac{2001}{2000}\right)^{2}\approx-5\cdot10^{-4}$$ $$\chi(500)=\frac{1501}{1500}-\left(\frac{1501}{1500}\right)^{2}\approx-6.6\cdot10^{-4}$$ $$\chi(1000)=\frac{1001}{1000}-\left(\frac{1001}{1000}\right)^{2}\approx-1\cdot10^{-3}$$ As we see, the autocorrelation value has doubled when check for half of the number of sweeps through the lattice.

Questions:

1. I'd expect that the correlation function would behave as an exponentially decaying function like $e^{-t/\tau}$ where $\tau$ is the correlation time, but rather, the plots show exponential growth with values as calculated above.
2. The expression for the autocorrelation as stated in this post is a discretization of $$\chi(t)=\int dt'(m(t')-\langle m\rangle)(m(t'+t)-\langle m \rangle)$$ which implies that if all values are 1, the integral should yield some constant, not an increasing function. What went wrong?
• The autocorrelation function is defined between -1 and 1. I think you have forgotton to normalize by the variance of m(t). en.wikipedia.org/wiki/Autocorrelation – Greg Petersen Feb 16 '17 at 18:52
• Thanks for answering, I could try that, but the texts I'm consulting on the subject simply divide everything by the first value in order to normalize to 1, which should be the greatest value. It won't fix the growth I'm afraid. The equation I'm using is the one on page 61 in the link I've provided, and it is from the first 4 chapters of a released book: amazon.com/Monte-Carlo-Methods-Statistical-Physics/dp/… – AlX Feb 16 '17 at 19:42
• The first value is the variance so that's good. Also, i think your sum needs to go from 1 to 2000. If you go from 0 to 2000, that is actually 2001 measurements. But you are dividing by 2000 instead of 2001 when you take the expectation value. – Greg Petersen Feb 16 '17 at 20:17

I agree with Greg Petersen: the normalization factor should be the number of terms in the sum, i.e. $t_{\rm max}-t+1$ and not $t_{\rm max}-t$ if your sum goes from $t=0$ to $t_{\rm max}$.
Your choice of a system where each spin in the lattice had spin 1 for every $t$ is pathological: the spins are infinitely correlated at anytime! It corresponds to a temperature $T=0$.
At finite temperature, autocorrelations are expected to decay exponentially fast (except at the critical point) $\chi(t)\sim \chi(0)e^{-t/\tau}$ so it is useless to choose a too large value of $t_{\rm max}$. In practise, $t_{\rm max}$ should be equal to several times $\tau$. By choosing a too large value of $t_{\rm max}$, most of the terms in your sum are just noise and for a finite system, you cannot expect them to simply cancel. For a $10\times 10$ Ising model, $\tau$ is certainly much smaller than 2000.
Your sum needs to go from $1$ to $2000$ instead of $0$ to $2000$. Essentially you are taking an average so it doesn't make sense to divide by fewer values than you are summing together (i.e. $2001/2000$). If you make this correction, all your calculations above work out to $0$. This leads to an indeterminate auto-correlation function as you are left with $0/0$ when normalizing by the variance.