\documentclass[../ising_model.tex]{subfiles} 

\begin{document} 
\section{Methods}\label{sec:methods}
\subsection{The Ising model}\label{subsec:ising_model}
The Ising model consists of a lattice of spins, which can be thought of as atoms 
in a grid. In two dimensions, the length of a lattice is given by $L$, and the 
number of spins within a lattice is given by $N = L^{2}$. When we consider the 
entire lattice, the system spin configuration is represented as a matrix $L \times L$
\begin{align*}
    \mathbf{s} &= \begin{pmatrix}
        s_{1,1} & s_{1,2} & \dots & s_{1,L} \\
        s_{2,1} & s_{2,2} & \dots & s_{2,L} \\
        \vdots & \vdots & \ddots & \vdots \\
        s_{L,1} & s_{L,2} & \dots & s_{L,L}
    \end{pmatrix} \ .
\end{align*} % $\mathbf{s} = [s_{1}, s_{2}, \dots, s_{N}]$.
The total number of possible spin configurations, also called microstates, is 
$|\mathbf{s}| = 2^{N}$. 

A given spin $i$ can take one of two possible discrete values $s_{i} \in \{-1, +1\}$. 
The spins interact with its nearest neighbors, and in a two-dimensional lattice 
each spin has up to four nearest neighbors. In our experiments we will use periodic 
boundary conditions, resulting in all spins having exactly four nearest neighbors. 
To find the analytical expressions necessary for validating our model implementation, 
we assume a $2 \times 2$ lattice. 

The hamiltonian of the Ising model is given by
\begin{equation}
    E(\mathbf{s}) = -J \sum_{\langle k l \rangle}^{N} s_{k} s_{l} - B \sum_{k}^{N} s_{k}\ ,
    \label{eq:energy_hamiltonian}
\end{equation}
where $\langle k l \rangle$ denotes a spin pair. $J$ is the coupling constant, and 
$B$ is the external magnetic field. For simplicity, we consider the Ising model 
where $B = 0$, and find the total system energy given by 
\begin{equation}
    E(\mathbf{s}) = -J \sum_{\langle k l \rangle}^{N} s_{k} s_{l} \ .
    \label{eq:energy_total}
\end{equation}
We count the neighboring spins using the pattern visualized in Figure \ref{fig:tikz_boundary}, 
to avoid counting a given spin pair several times.
\begin{figure}
    \centering
    \begin{tikzpicture}
        \draw[help lines] (0, 0) grid (2, 2);
        % \node[inner] (s1) at (0.5, 0.5) {s};
        \node (s1) at (0.5, 1.5) {$s_1$};
        \node (s2) at (1.5, 1.5) {$s_2$};
        \node (s3) at (0.5, 0.5) {$s_3$};
        \node (s4) at (1.5, 0.5) {$s_4$};
        \node[gray] (s12) at (2.5, 1.5) {$s_1$};
        \node[gray] (s34) at (2.5, 0.5) {$s_3$};
        \node[gray] (s13) at (0.5, -0.5) {$s_1$};
        \node[gray] (s24) at (1.5, -0.5) {$s_2$};
    
        \draw[red, ->] (s1.east) -- (s2.west);
        \draw[red, ->] (s2.east) -- (s12.west);
        \draw[red, ->] (s1.south) -- (s3.north);
        \draw[red, ->] (s2.south) -- (s4.north);
        \draw[red, ->] (s3.east) -- (s4.west);
        \draw[red, ->] (s4.east) -- (s34.west);
        \draw[red, ->] (s3.south) -- (s13.north);
        \draw[red, ->] (s4.south) -- (s24.north);
      \end{tikzpicture}
      \caption{Visualization of counting rules for each pair of spins in a 
      $2 \times 2$ lattice, where periodic boundary conditions are applied.}
      \label{fig:tikz_boundary}
\end{figure} 
We also find the total magnetization of the system, which is given by
\begin{equation}
    M(\mathbf{s}) = \sum_{i}^{N} s_{i} \ .
    \label{eq:magnetization_total}
\end{equation}
In addition, we have to consider the state degeneracy, the number of different 
microstates sharing the same value of total magnetization. In the case where we 
have two spins oriented up, the total energy have two possible values, as shown 
in Appendix \ref{sec:energy_special}.
\begin{table}[H]
    \centering
    \begin{tabular}{cccc} 
        \hline
        Spins up & $E(\mathbf{s})$ & $M(\mathbf{s})$ & Degeneracy \\
        \hline 
        4 & -8J & 4 & 1 \\
        3 & 0 & 2 & 4 \\
        \multirow{2}*{2} & 0 & 0 & 4 \\
            & 8J & 0 & 2 \\ 
        1 & 0 & -2 & 4 \\
        0 & -8J & -4 & 1 \\
        \hline
    \end{tabular}
    \caption{Values of the total energy, total magnetization, and degeneracy for the 
    possible states of a system for a $2 \times 2$ lattice, with periodic boundary 
    conditions.}
    \label{tab:lattice_config}
\end{table}
We calculate the analytical values for a $2 \times 2$ lattice, found in Table \ref{tab:lattice_config}.
However, we scale the total energy and total magnetization by the number of spins, 
to compare these quantities for lattices where $L > 2$. Energy per spin is given by
\begin{equation}
    \epsilon (\mathbf{s}) = \frac{E(\mathbf{s})}{N} \ ,
    \label{eq:energy_spin}
\end{equation}
and magnetization per spin given by 
\begin{equation}
    m (\mathbf{s}) = \frac{M(\mathbf{s})}{N} \ .
    \label{eq:magnetization_spin}
\end{equation}


\subsection{Statistical mechanics}\label{subsec:statistical_mechanics}
When we study the behavior of the Ising model, we have to consider the probability 
of a microstate $\mathbf{s}$ at a fixed temperature $T$. The probability distribution 
function (pdf) is given by 
\begin{equation}
    p(\mathbf{s} \ | \ T) = \frac{1}{Z} \exp^{-\beta E(\mathbf{s})},
    \label{eq:boltzmann_distribution}
\end{equation}
known as the Boltzmann distribution. This is an exponential distribution, where 
$\beta$ is given by 
\begin{equation}
    \beta = \frac{1}{k_{B}} \ ,
    \label{eq:beta}
\end{equation}
where and $k_{B}$ is the Boltzmann constant. $Z$ is a normalizing factor of the 
pdf, given by
\begin{equation}
    Z = \sum_{\text{all possible } \mathbf{s}} e^{-\beta E(\mathbf{s})} \ , 
    \label{eq:partition}
\end{equation}
and is known as the partition function. We derive $Z$ in Appendix \ref{sec:partition_function},
which gives us
\begin{equation*}
    Z = 4 \cosh (8 \beta J) + 12 \ .
\end{equation*}
Using the partition function and Equation \eqref{eq:boltzmann_distribution}, the 
probability of a microstate at a fixed temperature is given by 
\begin{equation}
    p(\mathbf{s} \ | \ T) = \frac{1}{4 \cosh (8 \beta J) + 12} e^{-\beta E(\mathbf{s})} \ .
    \label{eq:pdf}
\end{equation}
% Add something about why we use the expectation values?
We find the expected total energy
\begin{equation}\label{eq:energy_total_result}
	\langle E \rangle = -\frac{8J \sinh(8 \beta J)}{\cosh(8 \beta J) + 3} \ ,
\end{equation}
and the expected energy per spin
\begin{equation}\label{eq:energy_spin_result}
	\langle \epsilon \rangle = \frac{-2J \sinh(8 \beta J)}{ \cosh(8 \beta J) + 3} \ .
\end{equation}
We find the expected absolute total magnetization
\begin{equation}\label{eq:magnetization_total_result}
	\langle |M| \rangle = \frac{2(e^{8 \beta J} + 2)}{\cosh(8 \beta J) + 3} \ ,
\end{equation}
and the expected magnetization per spin
\begin{equation}\label{eq:magnetization_spin_result}
	\langle |m| \rangle = \frac{e^{8 \beta J} + 1}{2( \cosh(8 \beta J) + 3)} \ .
\end{equation}
We derive the analytical expressions for expectation values in Appendix 
\ref{sec:expectation_values}. 

We also need to determine the heat capacity 
\begin{equation}
    C_{V} = \frac{1}{k_{B} T^{2}} (\mathbb{E}(E^{2}) - [\mathbb{E}(E)]^{2}) \ ,
    \label{eq:heat_capacity}
\end{equation}
and the magnetic susceptibility 
\begin{equation}
    \chi = \frac{1}{k_{\text{B} T}} (\mathbb{E}(M^{2}) - [\mathbb{E}(M)]^{2}) \ .
    \label{eq:susceptibility}
\end{equation}
In Appendix \ref{sec:heat_susceptibility} we derive expressions for the heat 
capacity and the susceptibility, and find the heat capacity 
\begin{align*}
    \frac{C_{V}}{N} &= \frac{64J^{2} }{N k_{B} T^{2}} \bigg( \frac{3\cosh(8 \beta J) + 1}{(\cosh(8 \beta J) + 3)^{2}} \bigg) \ ,
\end{align*}
and magnetic susceptibility
\begin{align*}
    \frac{\chi}{N} &= \frac{4}{N k_{B} T} \bigg( \frac{3e^{8 \beta J} + e^{-8 \beta J} + 3}{(\cosh(8 \beta J) + 3)^{2}} \bigg) \ .
\end{align*}

For scalability of our model, we want to work with unitless spins. To obtain this, 
we introduce the coupling constant $J$ as the base unit for energy. With the 
Boltzmann constant we derive the remaining units, which can be found in Table 
\ref{tab:units}.
\begin{table}[H]
    \centering
    \begin{tabular}{ccc} % @{\extracolsep{\fill}}
        \hline
        Value & Unit &  \\
        \hline 
        $[E]$ & $J$ &  \\
        $[M]$ & $1$ & (unitless) \\
        $[T]$ & $J / k_{B}$ &  \\
        $[C_{V}]$ & $k_{B}$ &  \\
        $[\chi]$ & $1 / J$ &  \\
        \hline
    \end{tabular}
    \caption{Values of the total energy, total magnetization, and degeneracy for the possible states of a system for a $2 \times 2$ lattice, with periodic boundary conditions.}
    \label{tab:units}
\end{table}


\subsection{Phase transition and critical temperature}\label{subsec:phase_critical}
We consider the Ising model in two dimensions, with no external external magnetic 
field. At temperatures below the critical temperature $T_{c}$, the Ising model will 
magnetize spontaneously. When increasing the temperature of the external field, 
the Ising model transitions from an ordered to an unordered phase. The spins become 
more correlated, and we can measure the discontinuous behavior as an increase in 
correlation length $\xi (T)$ \cite[p. 432]{hj:2015:comp_phys}. At $T_{c}$, the 
correlation length is proportional to the lattice size, resulting in the critical 
temperature scaling relation
\begin{equation}
    T_{c}(L) - T_{c}(L = \infty) = aL^{-1} \ ,
    \label{eq:critical_infinite}
\end{equation}
where $a$ is a constant. For the Ising model in two dimensions, with a lattice of 
infinite size, the analytical solution of the critical temperature is 
\begin{equation}
    T_{c}(L = \infty) = \frac{2}{\ln (1 + \sqrt{2})} J/k_{B} \approx 2.269 J/k_{B}
    \label{eq:critical_solution}
\end{equation}
This result was found analytically by Lars Onsager in 1944 \cite{onsager:1944}. We can estimate the 
critical temperature of an infinite lattice, using the critical temperature of 
finite lattices of different sizes and linear regression. 


\subsection{The Markov chain Monte Carlo method}\label{subsec:mcmc_method}
Markov chains consist of a sequence of samples, where the probability of the next 
sample depends on the probability of the current sample. Whereas the Monte Carlo 
method introduces randomness to the sampling, which allows us to approximate 
statistical quantities \cite{tds:2021:mcmc} without using the pdf directly. We 
generate samples of spin configuration, to approximate $\langle \epsilon \rangle$, 
$\langle |m| \rangle$, $C_{V}$, and $\chi$, using the Markov chain Monte Carlo 
(MCMC) method. 

The Ising model is initialized in a random state, and the Markov chain is evolved 
until the model reaches an equilibrium state. However, generating new random states 
requires ergodicity and detailed balance. A Markov chain is ergoditic when all system 
states can be reached at every current state, whereas detailed balance implies no 
net flux of probability. To satisfy these criteria we use the Metropolis-Hastings 
algorithm, found in Algorithm \ref{algo:metropolis}, to generate samples of microstates. 

One Monte Carlo cycle consist of changing the initial configuration of the lattice,
by randomly flipping a spin. When a spin is flipped, the change in energy is evaluated as
\begin{align*}
    \Delta E &= E_{\text{after}} - E_{\text{before}} \ ,
\end{align*}
and accepted if $\Delta E \leq 0$. However, if $\Delta E > 0$ we have to compute
the Boltzmann factor
\begin{equation}
    p(\mathbf{s} | T) = e^{-\beta \Delta E}
    \label{eq:boltzmann_factor}
\end{equation} 
Since the total system energy only takes three different values, the change in 
energy can take $3^{2}$ values. However, there are only five distinct values 
$\Delta E = \{-16J, -8J, 0, 8J, 16J\}$, we derive these values in Appendix \ref{sec:delta_energy}.
We can avoid computing the Boltzmann factor, by using a look up table (LUT). 
We use $\Delta E$ as an index to access the resulting value of the exponential function, 
in an array. 
\begin{algorithm}[H]
\caption{Metropolis-Hastings Algorithm}
\label{algo:metropolis}
    \begin{algorithmic}
        \Procedure{Metropolis}{$lattice, energy, magnetization$}
        \For{ Each spin in the lattice }
        \State $ri \leftarrow \text{random index of the lattice}$
        \State $rj \leftarrow \text{random index of the lattice}$
        \State $dE \leftarrow 2 \cdot lattice_{ri,rj} \cdot (lattice_{ri,rj-1} + lattice_{ri,rj+1} 
        + lattice_{ri-1,rj} + lattice_{ri+1,rj} )$

        \State $rn \leftarrow \text{random number} \in [0,1)$
        \If{$rn \leq e^{- \frac{dE}{T}}$}
            \State $lattice_{ri,rj} = -1 \cdot lattice_{ri,rj}$
            \State $magnetization = magnetization + 2 \cdot lattice_{ri,rj}$
            \State $energy = energy + dE$
        \EndIf
        \EndFor
        \EndProcedure
    \end{algorithmic}
\end{algorithm}
The Markov process reaches an equilibrium after a certain number of Monte Carlo cycles,  
called burn-in time. After the burn-in time, the system state reflects the state of a real system,
and we can start sampling microstates. The probability distribution of the samples 
will tend toward the actual probability distribution of the system. 


\subsection{Implementation and testing}\label{subsec:implementation_test}
We implemented a test suite, and compared the numerical estimates to the analytical 
results from Section \ref{subsec:statistical_mechanics}. In addition, we set a tolerance to 
verify convergence, and a maximum number of Monte Carlo cycles to avoid longer 
runtimes during execution.

We used a pattern to access all neighboring spins, where all indices of the neighboring 
spins are put in an $L \times 2$ matrix. The indices are accessed using pre-defined 
constants, where the first column contain the indices for neighbors to the left 
and up, and the second column right and down. This method avoids the use of if-tests,
and takes advantage of the parallel optimization. 

We parallelized our code using both message passing interface (OpenMPI) and 
multi-threading (OpenMP). First, we divided the temperatures into smaller sub-ranges, 
and each MPI process received a sub-range of temperatures. Every MPI initialize a 
parallel region with a set of threads, which then initializes an Ising model and 
performs the Metropolis-Hastings algorithm. We limited the number of times threads 
are spawned and joined, by using single parallel regions, reducing parallel overhead. 
We used Fox \footnote{Technical specifications for Fox can be found at \url{https://www.uio.no/english/services/it/research/platforms/edu-research/help/fox/system-overview.md}}, 
a high-performance computing cluster, to run our program for $1$ and $10$ million 
Monte Carlo cycles when analyzing phase transitions. 



\subsection{Tools}\label{subsec:tools}
The Ising model and MCMC methods are implemented in C++, and parallelized using 
both \verb|OpenMPI| \cite{gabriel:2004:open_mpi} and \verb|OpenMP| \cite{openmp:2018}. 
We use the Python library \verb|matplotlib| \cite{hunter:2007:matplotlib} to produce 
all the plots, \verb|seaborn| \cite{waskom:2021:seaborn} to set the theme in the 
figures. We also use a linear regression method and the normal distribution method 
from the Python library \verb|SciPy|. To optimize our implementation, we used a 
profiling tool \verb|Score-P| \cite{scorep}.


\end{document}