Two-Point correlation function of Ising data is not giving expected results

Question

Due to the nature of my problem, and readers perhaps not having a physics background, I need to give some context to my problem which is intrinsic to it. The problem however is not a physics one, but a coding one.

Coming off a previous question of mine posted here, I managed to write a code which calculates the two-point auto correlation function of a 2D grid (composed of $\pm 1$ as elements). To test the code, I ran it on some 2D Ising data.

Before I go into more detail, here is my problem: For an Ising model at the critical temperature, the expected behaviour of the two-point autocorrelation function looks like :

$$f(r) \propto \frac{1}{r^\alpha}$$

where for an Ising model, one knows that $\alpha = \frac{1}{4}$. When I run my code however, I never get anything similar to $\frac{1}{4}$. So there are a few potential issues here:

1. My code is not working properly
2. The data sets I am using do not accurately correspond to the 2D Ising data at the critical temperature.
3. The non-linear models fit is not accurate enough.

Point 3 is rather unlikely, because I am using a simple power law model. With regards to point 2, I have used many different data sets. What I have been doing is looking at plots of 2D Ising data at the critical temperature, such as this:

and then by getting the binary image data per pixel, I can create a 2D array, of the image given by $\pm 1$:

TempDat = ImageData[Binarize[Import["IsingExample.png"]]];
ImData = Replace[TempDat, n_?(# == 0 &) -> -1, {2}];

From here I have the data in the necessary format to calculate my correlation functions.

Note, here you can find a website that generates a whole bunch of these 2D Ising data image plots based on Monte Carlo simulations. Just set temperature at 2.27, and let the simulation run a while to generate data.

Based on this I have little reason to believe that these plots do not correspond to what they claim to be: a 2D plot representing spin states(up or down) of an Ising model at the critical temperature. Thus my final conclusion is that something is not right in my code, but for the life of me, I can't figure out what.

Since on a 2D lattice, the spin-spin correlation is rotationally symmetric at the critical point, we only need to consider vectors along the lattice directions. The way I calculate the two point correlation function is as follows:

Let $\varphi(p_i)$ return the value of the grid at point $p_i$, this will either be a $1$ or $-1$, then $\varphi(p_i)\varphi(p_j)$ will be 1 if $\varphi(p_i)=\varphi(p_j)$ or $-1$ if $\varphi(p_i)\neq\varphi(p_j)$. Taking the average over all points of distance $r$ between points $p_i,pj$ we have that:

$$f(r) = <\varphi(p_i)\varphi(p_j)> - <\varphi(p_i)><\varphi(p_j)>,$$

where $<\varphi(p_i)>=<\varphi(p_j)>$ is just the average value of everypoint on the grid, code wise

ave = N@Mean[Mean[array]]^2

Now that I have covered the set up, here is my code: First we input the data from the ".csv" file generated by the image(code above):

SetDirectory@NotebookDirectory[];
arrayA = Import["Ising1Data.csv"];
fig = 
  ArrayPlot[arrayA, 
    ColorFunction -> (Blend[{White, Black}, #] &), ImageSize -> Large]

Next I create my modules to sample every possible pair of points of a distance $r$, along the $x$- directions and $y$-directions:

PPx = 
  Compile[{{r, _Integer}, {ave, _Real}, {grid, _Real, 2}},
    Module[{dimx, dimy, e1, e2},
      {dimx, dimy} = Dimensions[grid];
      Mean @ Flatten @
        Table[ 
          e1 = grid[[i, j]];
          e2 = grid[[Mod[i + r, dimx, 1], j]];     
          (e1*e2) - ave, 
          {i, 1, dimx}, {j, 1, dimy}]],
    CompilationTarget -> "C", RuntimeOptions -> "Speed"];

PPy = 
  Compile[{{r, _Integer}, {ave, _Real}, {grid, _Real, 2}}, 
    Module[{dimx, dimy, e1, e2},
      {dimx, dimy} = Dimensions[grid];
      Mean @ Flatten @ 
        Table[   
          e1 = grid[[i, j]];
          e2 = grid[[i, Mod[j + r, dimy, 1]]];
          (e1*e2) - ave, 
          {i, 1, dimx}, {j, 1, dimy}]],
    CompilationTarget -> "C", RuntimeOptions -> "Speed"];

Since the data is generated along a torus, we need to implement boundary conditions as we move the line segment across the data — hence the Mod functions.

Next we need to change the length of $r$ — we want to plot the average correlation as a function of $r$ - from 1 to $\frac{m}{2}$ where $m$ is the dimension of the grid. Again, this is due to the periodic boundary conditions.

array = Developer`ToPackedArray@arrayA;
ave = N @ Mean[Mean[array]]^2;
g = Dimensions[array];
gx = IntegerPart[g[[1]]/2];
gy = IntegerPart[g[[2]]/2];

Corrx = 
  Monitor[
    Table[Flatten[{i, PPx[i, ave, array]}, 1], {i,1,gx}], 
    {N[i/gx]*100}]; // AbsoluteTiming

Corry = 
  Monitor[
    Table[Flatten[{i, PPy[i, ave, array]}, 1], {i,1,gy}], 
    {N[i/gy]*100}]; // AbsoluteTiming

We then take the average of the correlations for each $r$ in the $x$ and $y$ directions and plot the results.

CorrAll = 
  Mean /@ GatherBy[Flatten[{Corrx, Corry}, 1], First];

pltCorrAll = ListPlot[PphiAll, PlotRange -> All, ImageSize -> Large]

Lastly, I do a non-linear model fit, note I only do a model fit for about $\sim 10$% of the of the dimensions of the grid, this is because one then gets too much influence across the periodic boundaries. Thus in this case we do a fit over the first $100$ values of $r$.

nlmCorr = 
   NonlinearModelFit[CorrAll[[1 ;; 100]], a*x^(-b), {a, b}, x, MaxIterations -> 10000];
Normal[PnlmPhi]
RegPlotCorr = 
  Plot[nlmCorr[x], {x, 0, 100}, PlotStyle -> {Red, Thick}, PlotRange -> All]l
pltCorrAll = 
  ListPlot[CorrAll[[1;;100]], PlotRange -> All, ImageSize ->Large];
Show[pltCorrAll, RegPlotCorr]

This is the result of all of the above

With the following regression model:

$$f(r) \propto \frac{1}{r^{0.3674}}.$$

My given exponent here is $0.3674$, which is off by a good $47$%. This is way too far off from the expected value of $0.25$. I have tested my code on pure random data, and it works as expected(showing an average correlation of $0$).

No matter what Ising data set I look at, I cannot reproduce the expected behavior of the two point correlation function. I am hoping that by posting on this forum someone can spot an error somewhere in either my code or coding approach to the problem. I would really like to get this working; that way I can begin to apply the code to various problems of mine.

I appologise for the length of the question, I just want to give as thorough an explanation of my problem as possible to make it easier for readers to understand what I am doing.

P.S If the moderator(s) of the site feel that this question is better posted on a different sector, I will gladly accept its migration.

Answer 1

After some further investigations and discussion, I discovered the solution to this problem. Since it has been un answered it might be instructive to future readers.

The problem here, is not so much the coding, but instead a minor overlook on the physics side. For a system that is exactly at the critical point, the correlation functions looks like

$ f(r) = \frac{A}{r^\alpha} $

The problem is that, it is computationally enormously expensive to produce an Ising grid exactly at the critical point. So what is done instead; is we say we are close to the critical point. This was indeed already pointed out in the question. However, what must then be introduced is the correlation length, $\xi$. This then changes the correlation function to the more general form:

$ f(r) = \frac{A}{r^\alpha}\exp\left({\frac{-r}{\xi}}\right) $

At the critical point $\xi\rightarrow \infty$, hence $\exp\left({\frac{-r}{\xi}}\right) \rightarrow 1$.

If we now use this function to fit the above, we obtain

This results in the correct and expected value of $\alpha \simeq 0.25$. As expected from the physics.