2D Ising model, understanding autocorrelations, Monte Carlo

Question

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?

Answer 1

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.

Answer 2

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.