# Two-point-correlation in the 3D ising model

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.

• (assuming your system is uniform) take the average over i of sigma_i sigma_i+r . The second version with the means substracted is actually the most useful in terms of defining the correlation length around the critical point. They coincide on the Paramagnetic Side (where the single site expectation values will be zero), but differ on the Ferromagnetic side. Commented Jan 4 at 16:14
• yes i use periodic boundary conditions.okay, so that means whenever i want to calculate the two-point correlation function i need to specify a specific radius r? I mean it's obvious, when you look at the formula, but what confuses me, is that at $T_C$ the Correlation length should diverge. This is something i could only see if i would let the radius approach infinity, right? also do i have to calculate the sum with all spin-combinations at distance r, or is it enough if i use only one axis, to save computation time? Commented Jan 4 at 16:29

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$$).

• 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.
• The magnetic susceptibility (up to a factor $kT$). (The precise definition is given in the other answer of mine you refer to in your post.) Commented Jan 4 at 18:05