Two-point-correlation in the 3D ising model

Question

I am currently coding a 3D (Monte-Carlo) implementation of the Ising model, using the single spin-flip & Wolff algorithm.

So far, I was able to calculate all the interesting observables, like $M$ and $\chi$.

Now I want to implement the two-point-correlation function, but I am struggling here both in understanding the function itself and how to code it.

So in David Landau and Kurt Binder's book "A guide to Monte-Carlo simulations in statistical physics" (2021), they introduce the two-point correlation function as follows:

$$\Gamma(r) = \langle \rho(0) \cdot \rho(r) \rangle$$

where $r$ is the spatial distance and $\rho$ is the quantity whose correlation is being measured. So from my understanding, this should then look like this in the Ising model:

$$\Gamma(r) = \langle \sigma_{0} \cdot \sigma_{i} \rangle$$

Where $\sigma_{0}$ would be an arbitrarily chosen spin, called "origin site".

Now I am confused because I don't see how one should think of this in terms of Monte-Carlo simulations. Does that mean I choose (in each iteration) one (random?) $\sigma_{0}$ spin and calculate its product with another random spin ($\sigma_{i}$) and save that value to calculate the mean over all iterations later? Is it always the same spin, in a sense the two-point function is always calculated for a "specific radius $r$", where $r$ is set by $\sigma_i$?

Or do I take the sum of that product over all possible spins, in that sense:

$$\langle \sigma_{0} \cdot \sigma_{i} \rangle = \frac{1}{N}\sum_{i}^N \sigma_{0} \cdot \sigma_{i}$$ for each iteration?

Then, on the other hand I see this formulation in that thread here: $$\sum_{i,j} \left[ \langle \sigma_i \cdot \sigma_j \rangle - \langle \sigma_i \rangle \cdot \langle \sigma_j \rangle \right]$$

which confuses me even more, because now where is the "origin site" $\sigma_0$? Where does that formula come from?

Thank you.

Answer 1

Does that mean i choose (in each iteration) one (random?) $\sigma_0$ spin and calculate it's product with another random spin ($\sigma_i$) and save that value to calculate the mean over all iterations later?

If you have sampled $S$ "independent" configurations $\sigma^1,\dots,\sigma^S$ from the Gibbs distribution, then the 2-point function associated to the spins at fixed vertices $i$ and $j$ can be computed using $$ \langle \sigma_i \sigma_j \rangle \approx \frac{1}{S}\sum_{k=1}^S \sigma^k_i\sigma^k_j. $$ By the law of large numbers, this will converge to the correct value when $S\to\infty$.

In principle, you could also average over the box, that is, compute $$ \frac{1}{S}\sum_{k=1}^S \frac{1}{N}\sum_{\ell}\sigma^k_{\ell}\sigma^k_{\ell+(j-i)}, $$ where the sum is over all the $N$ vertices of your system (assuming periodic boundary conditions). This would still converge to the same quantity, but you'd need more than the law of large numbers to prove this.

Then on the other hand i see this formulation in that thread here [...] which confuses me even more, because now where is the "origin site" σ0 ? Where does that formula come from?

The formula you quote is what you obtain when you differentiate twice the log of the partition function w.r.t. the magnetic field. This gives you (up to some constants) the magnetic susceptibility of the model. Note that what appears in this formula is not the sum of the 2-point functions, but the sum of the truncated 2-point functions $\langle\sigma_i\sigma_j\rangle-\langle\sigma_i\rangle\langle\sigma_j\rangle$. These are actually the natural statistical quantities (they are just the covariances between the spin random variables at $i$ and $j$).

Finally, concerning some other of your remarks in the comments:

• Yes, the sum over all 2-point functions $$ \sum_i \langle\sigma_0\sigma_i\rangle $$ diverges at the critical point (for an infinite system!). This is also true for the truncated 2-point functions.
• Beware that the 2-point function is not rotationally invariant, even at macroscopic distances, except at the critical point. So you cannot just fix some distance $r$ and consider any pair of spins $i$ and $j$ at distance $r$ from each other as being representative of all such pairs. Even the rate of exponential decay is anisotropic.