Ising model: How can I spot the critical point?

Question

Consider a zero-field Ising model with $N$ spins and periodic boundary conditions, with the Hamiltonian given by: $$H = -K \sum _{(ij)} s_i s_j\tag{1}$$ in 1D and 2D, where $K = \frac{J}{k_BT}$, where $J$ is the coupling constant energy. If I run a Monte Carlo simulation with $N$ particles with Metropolis sampling, how do I spot the critical $K_c$ value? Onsager's solution in 2D predicted a critical point at $K_c = 0.4406868$.

My question is, when I run a simulation with $N$ particles and I track the Hamiltonian per particle $(H/N)$ and the magnetization per particle $\left(\sum _i s_i /N\right)$, with $K$ values going from $0.1$ to $0.7$ in increments of $0.1$, how do I spot the region of the critical coupling constant?

There has to be a signature of the critical point that is stronger when $K=0.4$ and $K=0.5$ that should lead me to think that the critical value of $K$ is between these two numbers, right?

Answer 1

My question is, when I run a simulation with $N$ particles and I track the Hamiltonian per particle $(H/N)$ and the magnetization per particle $\left(\sum _i s_i /N\right)$, with $K$ values going from $0.1$ to $0.7$ in increments of $0.1$, how do I spot the region of the critical coupling constant?

There has to be a signature of the critical point that is stronger when $K=0.4$ and $K=0.5$ that should lead me to think that the critical value of $K$ is between these two numbers, right?

The thermal phase transition in the 2D Ising model is of continuous nature, meaning among other things that the magnetization per particle, $$ \frac{M}{N} = \left| \frac{\sum_i s_i}{N} \right| \propto \left( T_\mathrm{c} - T \right)^\beta,\quad T<T_\mathrm{c}, $$ where $\beta=1/8$ in 2D, and $M/N=0$ at $T>T_\mathrm{c}$. In other words, it decays continuously as the critical temperature is approached from below (which corresponds to reducing $K$). You can set a threshold (e.g. machine precision) and look for the phase transition this way (e.g. using a binary search), but it is an imprecise method. A much sharper signature is found in the specific heat, $C_V=\partial E/\partial T$, which (in 2D) diverges logarithmically as $T\rightarrow T_\mathrm{c}$. Generally, it is useful to look also at such derivatives (higher order derivatives of the free energy) when studying transitions. In the presence of a field, the magnetic susceptibility is useful, for example.

Also, what is the criticality condition for 1 dimension?

Actually, the 1D Ising chain in zero field does not have a finite-temperature phase transition. Instead it's disordered at all finite temperatures. However, there exists a critical point at $T=0$.

Answer 2

As Anyon correctly pointed out, there is no phase transition at finite temperature in 1D.

In 2D there are a number of different ways to identify the phase transition (I'm assuming you're using Monte Carlo). You could directly look at the magnetization, but a more reliable signature is the magnetic susceptibility, $\chi_m (K)$, which is strongly peaked around $T_c$ (the same is true for the specific heat). This is a much stronger signal.

$\chi$ can be described by a partial derivative w.r.t. the field, but a more useful way to measure it is in terms of the fluctuations in the magnetization: $$\chi = \frac{d \langle m \rangle}{dh} \Big|_{h\to 0} = \frac{N}{kT} \left( \langle m^2 \rangle - \langle | m |\rangle^2 \right) \tag{1}$$ Using this definition avoids the need for numerical derivatives (there is a similar trick for the specific heat). This method is described on p. 34 of this excellent review article by Sandvik.

But, any simulation is at finite size. To make an even more precise measurement, you could take advantage of finite-size scaling, which allows you use the finite-size effects to determine the $L\to\infty$ transition point. Using the critical exponents (which are known), you can collapse $\chi(T,L)$ curves onto a single rescaled $$L^{-\kappa/\nu} \chi(T,L) = t L^{1/\nu} \tag{2}$$ (using the reduced temperature $t \equiv \frac{T-T_c}{T_c}$). For a full description of this approach, see Sandvik arXiv:1101.3281, p. 35.

Answer 3

I think Anyon's and taciteloquence's answers are perfect. I just want to add an emphasis on the following fact that frequently leads to confusion for beginners.

The formal definition of the magnetization \begin{equation} m = \frac{\sum_i s_i}{N} \end{equation} has a symmetry that $\mathrm{Prob}[m=+m_0]=\mathrm{Prob}[m=-m_0]$, since the energy of a particular configuration and a totally flipped version will have exactly the same energy. Therefore, $\langle m\rangle =0$ for all temperature for a finite system. Only in an infinite system (thermodynamic limit) can this symmetry be broken, which is called spontaneous symmetry breaking.

In principle, if you manage to correctly calculate the expectation value $\langle m\rangle$, it will always be 0, but when you are running for only a short time with a large system, the Monte Carlo simulation can be trapped in one side of the two macro states (corresponding to $\langle m\rangle =\pm m_0$). In order to correctly spot when this starts to happen (i.e. where the transition is), one way is to calculate the expectation value of the absolute magnetization $\langle |m|\rangle$ as Anyon suggests. As you go to larger system sizes, $\langle |m|\rangle\rightarrow0$ in the disordered phase, and $\langle |m|\rangle\rightarrow +m_0$ in the ordered phase. Another way is to calculate the squared value $\langle m^2\rangle$ also suggested by taciteloquence. This also tends to either $0$ or $m_0^2$ depending on which phase you're in. $\langle m^2\rangle$ can be more common since it corresponds to the magnetic susceptibility when multiplied by the system volume, thus behaving more nicely.

There are also other ways of determining the critical point such as spotting the "crossing point" for the Binder cumulant $\langle m^4\rangle/\langle m^2\rangle ^2$, or the correlation ratio $\xi/L$.

Finally, looking at the histogram of the order parameter $m$ is always a choice. You measure the probability of the magnetization $m$ being a particular value $m^*$, and do that for all $m^*$, getting a histogram $P(m)$ telling you the distribution of possible $m$ for a particular temperature. If this distribution is a Gaussian (when the system is large enough) centered at $m=0$, then you're in the disordered phase. If it is a combination of two Gaussian distributions centered at $m=\pm m_0$, you're in the ordered phase. At the critical point, you have a very broad (non-Gaussian) distribution at the center, barely becoming two separate peaks. All of the common methods I mentioned above are essentially different ways of detecting this "splitting" process of the Gaussian histogram.