Numerical Ising Model - Wolff algorithm and correlations I'm doing some numerical Monte Carlo analysis on the 2 dimensional Ising model at the critical point. I was using the Metropolis 'single flip' evolution at first with success, though it suffers from critical slow down and makes studying large lattices unlikely possible. I'm now looking at cluster flip algorithms, specifically the Wolff algorithm.
I managed to implement it, and it looks to be working as it should (flipping a unique spin at $T = +\infty$, the whole lattice at $T=0$, matches the right energy density in the thermodynamic limit...) but I don't get the right behaviour for the two point $<\sigma_i\sigma_j>$ correlation function.
According to CFT it should behave like:
$$<\sigma_i\sigma_j> \;\propto \;\frac{1}{|i-j|^{\frac{1}{4}}}$$
I'm more and more convinced that it has to do with boundary conditions, I use non periodic free boundaries. The literature on the subject doesn't say much on this point.
Am I missing a subtlety (or an evidence) in this procedure, or in the use of this algorithm?
 A: *

*It is only exactly at the critical temperature that this CFT result works. You haven't mentioned if you have used the critical temperature when you did the monte-carlo.

*At/near critical point, autocorrelation time becomes huge. (If I am not mistaken, autocorrelation time must blow up exactly at critical temperature, however it is cut-off due to finiteness of system). So, its best to record measurements once every 3-5 autocorrelation times. Notice that, for a given random seed, autocorrelation effects can give a systematic error, while statistical error due to autocorrelation effect can be estimated by using different random seeds. Look into this book: http://www.amazon.com/Monte-Carlo-Methods-Statistical-Physics/dp/0198517971 

*My guess is that, as long as u are 10-20 spins away from the boundary when u do the measurement of spin-spin corr., you wont see boundary effects. This is only a guess. Why dont you put in periodic boundary conditions and see? The neighbour of i,j is (i+1)%L, (j+1)%L and 3 others like this.
A: I have tried to reproduce your problem using my own code, but I couldn't: I got the correct value for the exponent. I can't tell you what went wrong in your calculation, but I can tell you exactly what I did!
You can inspect my code on GitHub. It was the first thing I ever wrote in C++ and today I would do many things differently, so please don't judge the code too harshly. There is a readme file in the repository that explains pretty much everything about it.
I uploaded my runscript, the relevant part of the corresponding results and the small Pyxplot script that does the fit into this GitHub Gist. My value for the exponent is:
$$\beta = 0.243 \pm 0.001$$
You can see the fit plotted together with the Monte Carlo data in the sscorr_fit.pdf file.
Here is a list of things that I think might be relevant.


*

*I'm sure you know that, but just for the record: In the equation in the question one is not supposed to take the absolute value of the difference of the indices $i$ and $j$ of the two sites. The indices you assign to the sites are completely arbitrary. What matters is the distance between the two sites, so I would prefer to write it like this:
$$\langle s_i s_j \rangle \propto \frac{1}{|\vec r_j - \vec r_i|^\frac{1}{4}}$$

*I only measured the correlation along the direction of the lattice vectors and not along any diagonal. That should not matter though since there is a proof that the spin-spin correlation is rotationally symmetric at the critical point.

*I only fitted your formular in the range $1 \leq |\vec r_j - \vec r_i| \leq 50$ for a 512x512 model. I think if you go beyond that you get too much influence across the periodic boundary conditions. You can definitely not fit the entire $1 \leq |\vec r_j - \vec r_i| \leq 256$ range, as the correlation must have a minimum at $|\vec r_j - \vec r_i| = 256$.

*The correlation function changes quite quickly with temperature around the critical point, so make sure you have the right temperature and everthing is equilibrated. In my experience it is best to start the Wolff algorithm with a model where all spins point in the same direction. If you take random spins and then cool down from $T=\infty$ to $T_c$ the Wolff algorithm is extremely inefficient in the beginning as it can not build large clusters in this noise of random spins.
Hope that helps!
