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:

2D Ising data on a Torus(periodic boundary conditions)

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):

arrayA = Import["Ising1Data.csv"];
fig = 
    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 @
          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 @ 
          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 = 
    Table[Flatten[{i, PPx[i, ave, array]}, 1], {i,1,gx}], 
    {N[i/gx]*100}]; // AbsoluteTiming

Corry = 
    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];
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

Non-linear model fit of the data

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.

  • $\begingroup$ I did not go through all of your post, but here's an idea/observation: your system is finite. Phase transitions happen in the thermodynamic limit. For finite systems, not everything works the same way. E.g., the magnetic susceptibility has a maximum (instead of diverging), and the location of this maximum will depend on the size of the system. In other words, the best approximation of the critical behaviour may not be exactly at 2.27. $\endgroup$
    – Szabolcs
    Mar 2, 2017 at 18:05
  • 1
    $\begingroup$ There is something called "finite size effects" within the topic of Ising models - which as you stated is introduced when you consider the discrete version of the theory. I will be honest, if this is the case, I'm not sure how it's affecting the data. Since however, I have seen other people claim that they have recovered the expected two-point correlation behaviour for their simulations; I am inclined to believe that, even for a discrete approximation, the critical behaviour still occurs at T=2.27, and thus the correlations should reflect this. $\endgroup$ Mar 3, 2017 at 9:19
  • $\begingroup$ See this link for someone who claims they get the correct correlation behaviour. $\endgroup$ Mar 3, 2017 at 9:21
  • $\begingroup$ Instead of going to NonlinearModelFit, try plotting your function on a log-log scale and see how the slope (representing the exponent) changes. $\endgroup$
    – Szabolcs
    Mar 3, 2017 at 11:50
  • $\begingroup$ This is the susceptibility and magnetization, both indicators of critical behaviour, versus $\beta = 1/T$ in 2D Ising models of increasing lattice sizes. The vertical line shows the critical point in the thermodynamic limit. I am not saying that the effect is non-negligible for a $1000\times 1000$ lattice as well. $\endgroup$
    – Szabolcs
    Mar 3, 2017 at 13:08

1 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

enter image description here

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.