The Ising model

The model

The Euclidean phi-four theory in $d$ dimensions with $O(N)$ symmetry is given by the action
\begin{eqnarray}
S[\phi]=\int d^dx \bigg(\frac{1}{2}(\partial_{\mu}\phi^i)^2+\frac{\mu_0^2}{2}\phi^i\phi^i+\frac{\lambda}{4}(\phi^i\phi^i)^2\bigg).\label{contaction}
\end{eqnarray}
We will employ lattice regularization in which $x=an$, $\int d^dx =a^d\sum_n$, $\phi^i(x)=\phi_n^i$ and $\partial_{\mu}\phi^i=(\phi_{n+\hat{\mu}}^i-\phi_n^i)/a$. The lattice action reads then
\begin{eqnarray}
S[\phi]
&=&\sum_n\bigg(-2\kappa \sum_{\mu}\Phi_n^i\Phi_{n+\hat{\mu}}^i+{\Phi}_n^i\Phi_n^i+g(\Phi_n^i\Phi_n^i-1)^2 \bigg).
\end{eqnarray}
The mass parameter $\mu_0^2$ is replaced by the so-called hopping parameter $\kappa$ and the coupling constant $\lambda$ is replaced by the coupling constant $g$ where
\begin{eqnarray}
{\mu}_{0L}^2\equiv \mu_0^2a^2=\frac{1-2g}{\kappa}-2d~,~\lambda_{L}\equiv \frac{\lambda}{a^{d-4}}=\frac{g}{\kappa^2}.
\end{eqnarray}
The fields $\phi_n^i$ and $\Phi_n^i$ are related by an irrelevant rescaling.

In the simulations we will start by fixing the lattice quartic coupling $\lambda_L$ and the lattice mass parameter $\mu_{0L}^2$ which then allows us to fix $\kappa$ and $g$ as
\begin{eqnarray}
\kappa=\frac{\sqrt{8\lambda_L+(\mu_{0L}^2+4)^2}-(\mu_{0L}^2+4)}{4\lambda_L}.
\end{eqnarray}
\begin{eqnarray}
g=\kappa^2\lambda_L.
\end{eqnarray}
The phase diagram will be drawn originally in the $\mu_{0L}^2-\lambda_L$ plane which is the lattice phase diagram.

The limit  $g\longrightarrow \infty $ of the $O(1)$ model is precisely the Ising model in $d$ dimensions.

The mean field approximation

There are two phases in this model. A disordered (paramagnetic) phase characterized by $\langle \Phi_n^i\rangle=0$  and an ordered (ferromagnetic) phase characterized by $\langle \Phi_n^i\rangle=v_i\neq 0$.

This can be seen, for example, by using  the mean field approximation which is an exact prescription in dimension $d=4$. In this method, we approximate the spins $\Phi_n^i$ at the $2d$ nearest neighbors of each spin $\Phi_n^i$ by the average $v^i=\langle \Phi_n^i\rangle$, viz
\begin{eqnarray}
\frac{\sum_{\mu}(\Phi_{n+\hat{\mu}}^i+\Phi_{n-\hat{\mu}}^i)}{2d}=v^i.
\end{eqnarray}
However, from the definition of $\langle \Phi_n^i\rangle$  we have
 \begin{eqnarray}
\langle \Phi_n^i\rangle
&=&\frac{\int d\mu(\Phi)~\Phi_n^i e^{2\kappa\sum_n\sum_{\mu}\Phi_n^i\Phi_{n+\hat{\mu}}^i}}{\int d\mu(\Phi)~e^{2\kappa\sum_n\sum_{\mu}\Phi_n^i\Phi_{n+\hat{\mu}}^i}}
\end{eqnarray}
In the mean field approximation this definition becomes a condition on $v^i$ given explicitly by
\begin{eqnarray}
v^i
&=&\frac{\int d\mu(\Phi)~\Phi_n^i e^{4\kappa d \sum_n \Phi_n^i v^i}}{\int d\mu(\Phi)~e^{4\kappa d \sum_n\Phi_n^i v^i}}.
\end{eqnarray}
From the one hand, we get in the limit $g\longrightarrow 0$ the solution
\begin{eqnarray}
v^i=2\kappa d v^i+...\Rightarrow \kappa_c=\frac{1}{2d}.
\end{eqnarray}
On the other hand, in the limit $g\longrightarrow \infty$ we get the solution
\begin{eqnarray}
v^i=\frac{4\kappa d v^i}{N}+...\Rightarrow \kappa_c=\frac{N}{4d}.
\end{eqnarray}
The critical line $\kappa_c=\kappa_c(g)$ interpolates in the $\kappa-g$ plane between these two lines.

For example, for the $O(1)$ model in the limit $g\longrightarrow \infty$ the solution must satisfy the condition
\begin{eqnarray}
v=\tanh 4\kappa d v.
\end{eqnarray}
For $\kappa<\kappa_c$ (where $\kappa_c=1/4d$) there is only one solution at $v=0$ whereas for $\kappa>\kappa_c$ there are two solutions away from the zero, i.e. $v\neq 0$.  Clearly for $\kappa$ near $\kappa_c$ the solution $v$ is near $0$ and thus we can expand the above equation as

\begin{eqnarray}
\frac{1}{3}(4d)^2\kappa^3 v^2=\kappa-\kappa_c.
\end{eqnarray}
Thus only for $\kappa>\kappa_c$ there is a non-zero solution.

In summary we have the two phases
\begin{eqnarray}
\kappa>\kappa_c:~{\rm broken, ordered, ferromagnetic}
\end{eqnarray}
\begin{eqnarray}
\kappa<\kappa_c:~{\rm symmetric, disordered, paramagnetic}.
\end{eqnarray} 

The one-loop result in four dimensions

As we have already mentioned the mean field approximation becomes exact in dimension four.  This translates into the fact that in four dimensions the critical behavior is fully controlled by the Gaussian fixed point (since the mean field approximation  effectively reduces the theory to free fields).

Indeed, for $d=4$ the critical value at $g=0$ is $\kappa_c=1/8$ for all $N$ which can be derived from the one-loop order of the continuum phi-four  theory with $O(N)$ symmetry as follows. The renormalized mass at one-loop order in this theory  is given by the equation (with $a=1/\Lambda$)

\begin{eqnarray}
a^2m_R^2
&=&am^2+\frac{(N+2)\lambda}{16\pi^2}+\frac{(N+2)\lambda}{16\pi^2}a^2m^2\ln a^2m^2\nonumber\\
&+&\frac{(N+2)\lambda}{16\pi^2}a^2m^2{\bf C}+a^2\times {\rm finite}~{\rm terms}.
\end{eqnarray}
Thus we obtain in the continuum limit $a\longrightarrow 0$ the result
\begin{eqnarray}
a^2m^2&\longrightarrow &-\frac{(N+2)\lambda}{16\pi^2}-\frac{(N+2)\lambda}{16\pi^2}a^2m^2\ln a^2m^2\nonumber\\
&-&\frac{(N+2)\lambda}{16\pi^2}a^2m^2{\bf C}-a^2\times {\rm finite}~{\rm terms}.
\end{eqnarray}
In other words, we have (with $r_0=(N+2)/8\pi^2$)
\begin{eqnarray}
a^2m^2\longrightarrow a^2m_c^2=-\frac{r_0}{2}\lambda+O(\lambda^2).
\end{eqnarray}
This is the critical line for small values of the coupling constant which when expressed in terms of $\kappa$ and $g$ we obtain

\begin{eqnarray}
\kappa \longrightarrow \kappa_c=\frac{1}{8}+(\frac{r_0}{2}-\frac{1}{4})g+O(g^2).
\end{eqnarray}
The continuum limit $a\longrightarrow 0$ corresponds precisely to the limit in which the mass approaches its critical value. This happens for every value of the coupling constant and hence the continuum limit  $a\longrightarrow 0$ is the limit in which we approach the critical line. The continuum limit is therefore a second order phase transition.

The continuum limit in two dimensions: Monte Carlo and renormalization group equation methods

We are interested in $d=2$ and $N=1$ which is the Ising model in two dimensions. This shares with the more realistic case in $d=3$ and $N=1$ (which describes all second order phase transitions in nature) the same critical behavior. However, in two dimensions the model admits a quasi-exact solution in terms of quantum field theory (super-renormalizability) and the renormalization group equation (Wilson-Fisher fixed point) which in the limit $g\longrightarrow\infty$ should reduce to the remarkable Onsager's exact solution.

Monte Carlo method

In two dimensions the lattice points are labelled by  $n=(n_1,n_2)$ where $n_1=1,...,L_1$ and  $n_2=1,...,L_2$.  We will also impose periodic boundary conditions, i.e. $\Phi_{n_1+L_1,n_2}\equiv \Phi_{n_1,n_2}$ and $\Phi_{n_1,n_2+L_2}\equiv \Phi_{n_1,n_2}$.

The lattice action for $\Phi^4$ theory in two dimensions with $O(N)$ symmetry reads explicitly


\begin{eqnarray}
S[\phi]
&=&\sum_{n_1=1}^{L_1}\sum_{n_2=1}^{L_2}\bigg(-2\kappa \Phi_{n_1,n_2}^i(\Phi_{n_1+1,n_2}^i+\Phi_{n_1,n_2+1}^i)+{\Phi}_{n_1,n_2}^i\Phi_{n_1,n_2}^i+g(\Phi_{n_1,n_2}^i\Phi_{n_1,n_2}^i-1)^2 \bigg).\label{act}
\end{eqnarray}
This Euclidean quantum field theory is a statistical dynamical system with a Boltzmann weight given by
\begin{eqnarray}
{\cal P}[\Phi]=\frac{\exp(-S[\Phi])}{Z}.
\end{eqnarray}
$Z$ is the partition function of the system which is defined in an obvious way.

In Monte Carlo simulations we can use the Metropolis algorithm to generate a collection of thermalized configurations $\Phi_n^i$ according to this probability distribution.

Indeed, in the Metropolis algorithm we make small changes to the field configuration, for example  $\Phi_{p_1,p_2}^j\longrightarrow \Phi_{p_1,p_2}^j+h$, then accept/reject these changes using precisely the Boltzmann probability distribution
\begin{eqnarray}
W(\Phi_{p_1,p_2}^j\longrightarrow \Phi_{p_1,p_2}^j+h) ={\rm min}(1,\exp(-\Delta S_{p_1,p_2}^j(h)).\label{met}
\end{eqnarray}
The variation $\Delta S_{p_1,p_2}^j(h)$ due to the change $\Phi_{p_1,p_2}^j\longrightarrow \Phi_{p_1,p_2}^j+h$ is given explicitly by
\begin{eqnarray}
 \Delta S_{p_1,p_2}^j(h)&=&-2\kappa h\big(\Phi_{p_1+1,p_2}^j+\Phi_{p_1-1,p_2}^j+\Phi_{p_1,p_2+1}^j+\Phi_{p_1,p_2-1}^j\big)+(1-2g)\big(h^2+2h\Phi_{p_1,p_2}^j\big)\nonumber\\
&+&g\big(h^2+2h\Phi_{p_1,p_2}^j\big)\big(h^2+2h\Phi_{p_1,p_2}^j+2\Phi_{p_1,p_2}^i\Phi_{p_1,p_2}^i\big).\label{var}
\end{eqnarray}
In the actual code (written in Fortran which is our preferred programming language) we will go through the following steps:

- We write a subroutine which computes the energy given by the action (\ref{act}). We can also include other primary observables in this subroutine for example the magnetization
\begin{eqnarray}
M=\frac{1}{NL_1L_2}\langle m\rangle~,~m=|\sum_{i,n_1,n_2}\Phi_{n_1,n_2}^i|.
\end{eqnarray}

-We write a subroutine which computes the variation (\ref{var}). This depends on $j$, $p_1$, $p_2$ and $h$.

- In preparation for writing the Metropolis process we include a random number generator. We prefer Mersenne Twistor or ran2 of numerical recipes.

- We write a subroutine which implements the Metropolis algorithm  (\ref{met}). So we change all the fields at all lattice points successively and accept/reject according to the probability (\ref{met}). The value of the variation $h$ is drawn randomly and its range is tuned so that the acceptance rate is fixed around $1/3$.

- In preparation for writing thermalization and measurement processes we write a subroutine which computes the average, the error and the auto-correlation time for arbitrary observables.

- We write down the main program in which the random number generator, the physical parameters of the model and the field configuration are initialized followed by a  thermalization process which changes the system until thermalization (or equilibrium) is reached then followed by a measurement process in which a  set of thermalized configurations $\Phi_n^i$ is collected and  primary observables are measured.

- We measure the secondary observables such as second moments like the specific heat (associated with the energy) and susceptibility (associated with  the magnetization).  These are given by
\begin{eqnarray}
C_v=\langle S^2\rangle -\langle S\rangle^2.
\end{eqnarray}

\begin{eqnarray}
\chi=\langle m^2\rangle -\langle m\rangle^2.
\end{eqnarray}
Renormalization group equation

In the simulations we will start by fixing the lattice quartic coupling $\lambda_L$ and the lattice mass parameter $\mu_{0L}^2$ which then allows us to fix $\kappa$ and $g$. The phase diagram will be drawn originally in the $\mu_{0L}^2-\lambda_L$ plane which is the lattice phase diagram. This should be extrapolated to the infinite volume limit  defined by $l_1=L_1a\longrightarrow \infty$ and  $l_2=L_2a\longrightarrow \infty$ with $a$ fixed.

However, the Euclidean quantum field theory phase diagram should be drawn in terms of the renormalized parameters and is obtained from the lattice phase diagram by taking the limit $a\longrightarrow 0$.

In two dimensions the $\Phi^4$ theory requires only mass renormalization while the quartic coupling constant is finite (recall that in four dimensions three renormalizations are required: the mass, the quartic coupling and the wave function renormalizations).

Indeed, in two dimensions the bare mass $\mu_0^2$ diverges logarithmically when we remove the cutoff, i.e. in the limit $\Lambda\longrightarrow \infty$ where $\Lambda=1/a$ while $\lambda$ is independent of $a$. As a consequence, the lattice parameters will go to zero in the continuum limit $a\longrightarrow 0$.

We know that mass renormalization is due to the tadpole diagram which is the only divergent Feynman diagram in the theory and takes the form of a simple reparametrization given by
\begin{eqnarray}
\mu_0^2=\mu^2-\delta\mu^2.\label{ren}
\end{eqnarray}
$\mu^2$ is the renormalized mass parameter and $\delta\mu^2$ is the counter term which is fixed via an appropriate renormalization condition. The unltraviolet divergence $\ln\Lambda$  of $\mu_0^2$ is contained in $\delta\mu^2$ while the renormalization condition will split the finite part of $\mu_0^2$ between $\mu^2$ and $\delta\mu^2$. The choice of the renormalization condition can be quite arbitrary. A convenient choice suitable for Monte Carlo measurements and which distinguishes between the two phases of the theory is given by the usual normal ordering prescription \cite{Loinaz:1997az} .




Quantization at one-loop  (together with a self-consistent Hartree treatement) gives in the symmetric phase where $\mu^2>0$ the $2-$point function
\begin{eqnarray}
\Gamma^{(2)}(p)
&=&p^2+\mu^2+\Sigma(p)~,~\Sigma(p)=3\lambda A_{\mu^2}-\delta\mu^2+....
\end{eqnarray}
The dots $...$ denote two-loop effect. $A_{\mu^2}$ is precisely the value of the tadpole diagram given by
\begin{eqnarray}
A_{\mu^2}=\int\frac{d^2k}{(2\pi)^2}\frac{1}{k^2+\mu^2}.
\end{eqnarray}
The renormalization condition which is equivalent to normal ordering the interaction  is equivalent to the choice
 \begin{eqnarray}
\delta\mu^2=3\lambda A_{\mu^2}.\label{ren1}
\end{eqnarray}
A dimensionless coupling constant can the be defined by
\begin{eqnarray}
f=\frac{\lambda}{\mu^2}.
\end{eqnarray}
The action becomes
\begin{eqnarray}
S[\phi]=\int d^2x \bigg(\frac{1}{2}(\partial_{\mu}\phi)^2+\frac{1}{2}\mu^2(1-3fA_{\mu^2})\phi^2+\frac{f\mu^2}{4}\phi^4\bigg).
\end{eqnarray}
For sufficiently small $f$ the exact effective potential is well approximated by the classical potential with a single minimum at $\phi_{\rm cl}=0$. For larger $f$, the coefficient of the mass term in the above action can become negative and as a consequence a transition to the broken symmetry phase is possible, although in this regime the effective potential is no longer well approximated by the classical potential.  Indeed, a transition to the broken symmetry phase was shown to be present in \cite{Chang:1976ek}, where a duality between the strong coupling regime of the above action and a weakly coupled theory normal-ordered with respect to the broken phase was explicitly constructed.



On a finite volume lattice with periodic boundary conditions equation (\ref{ren}), together with equation (\ref{ren1}), becomes
\begin{eqnarray}
\mu_L^2-3\lambda_L A_{\mu_L^2}-\mu_{0L}^2=0.\label{reno}
\end{eqnarray}
The Feynman diagram $A_{\mu^2}$ takes explicitly the form
\begin{eqnarray}
A_{\mu^2}
&=&\frac{1}{N^2}\sum_{n_1=1}^N\sum_{n_2=1}^N\frac{1}{4\sin^2 {\pi n_1}/{N}+4\sin^2{\pi n_2}/{N}+\mu_L^2}.
\end{eqnarray}
Given the critical value of $\mu_{0L}^2$ for every value of $\lambda_L$ we need then to determine the corresponding critical value of $\mu_{L}^2$. This can be done numerically using the Newton-Raphson algorithm. The continuum limit $a\longrightarrow 0$ is then given by extrapolating the results into the origin, i.e. taking $\lambda_L=a^2\lambda\longrightarrow 0$, $\mu_L^2=a^2\mu^2\longrightarrow 0$ in order to determine the critical value
\begin{eqnarray}
f_*={\rm lim}_{\lambda_L,\mu_L^2\longrightarrow 0}\frac{\lambda_L}{\mu_{L*}^2}.
\end{eqnarray}

References

%\cite{Loinaz:1997az}
\bibitem{Loinaz:1997az}
  W.~Loinaz and R.~S.~Willey,
  'Monte Carlo simulation calculation of critical coupling constant for continuum phi**4 in two-dimensions,'
  Phys.\ Rev.\ D {\bf 58}, 076003 (1998)
  [hep-lat/9712008].
  %%CITATION = HEP-LAT/9712008;%%
  %31 citations counted in INSPIRE as of 23 May 2015

%\cite{Chang:1976ek}
\bibitem{Chang:1976ek}
  S.~J.~Chang,
  'The existence of a second order phase transition in the two-dimensional phi**4 field theory,'
  Phys.\ Rev.\ D {\bf 13}, 2778 (1976)
  [Phys.\ Rev.\ D {\bf 16}, 1979 (1977)].
  %%CITATION = PHRVA,D13,2778;%%
  %113 citations counted in INSPIRE as of 25 May 2015