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?
Best Answer
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!