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

\begin{document} 
\section{Methods}\label{sec:methods}
\subsection{The Ising model}\label{subsec:ising_model}
The Ising model consist of a lattice of spins, which can be though of as atoms 
in a grid. In two dimensions, the length of the 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 denoted as a vector $\mathbf{s} 
= [s_{1}, s_{2}, \dots, s_{N}]$. The total number of spin configurations, or 
microstates, is $2^{N}$. 

A given spin $i$ can take one of two possible discrete values $s_{i} \in \{-1, +1\}$, 
either spin down or spin up. The spins interact with its nearest neighbors, and 
in a two-dimensional lattice each spin has up to four nearest neighbors. However, 
the model is not restricted to this dimentionality \cite[p. 3]{obermeyer:2020:ising}. 
In our experiment we will use periodic boundary conditions, meaning all spins 
have exactly four nearest neighbors. To find the analytical expressions necessary 
for validating our model implementation, we will assume a $2 \times 2$ lattice. 

We count the neighboring spins as visualized in Fig. \ref{fig:tikz_boundary}, 
\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}
and find the total system energy given by % Eq. \eqref{eq:energy_total}
\begin{equation}
    E(\mathbf{s}) = -J \sum_{\langle k l \rangle}^{N} s_{k} s_{l} \ .
    \label{eq:energy_total}
\end{equation}
$\langle k l \rangle$ denote a pair of spins, and $J$ is the coupling constant. 
We also find the total magnetization of the system, which is given by % Eq. \eqref{eq:magnetization_total}
\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. The derivation of the 
analytical values can be found in Appendix % \ref{}
\begin{table}[H]
    \centering
    \begin{tabular}{cccc} % @{\extracolsep{\fill}}
        \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 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 use the analytical values, found in Table for both for lattices where $L = 2$ and $L > 2$. 

However, to compare the quantities for lattices where $L > 2$, we find energy 
per spin 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 ferromagnetism, we have to consider the probability for 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$ and $Z$ are given by
\begin{align*}
    \beta =& \frac{1}{k_{B}} \ , &
    Z &= \sum_{\text{all possible } \mathbf{s}} e^{-\beta E(\mathbf{s})} \ , \\
\end{align*}
and $k_{B}$ is the Boltzmann constant. $Z$ is a normalizing factor of the pdf, 
known as the partition function, which we derive in Appendix \ref{sec:partition_function}
\begin{equation*}
    Z = 4 \cosh (8 \beta J) + 12 \ .
\end{equation*}
Using the partition function and Eq. \eqref{eq:boltzmann_distribution}, the pdf 
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})} \ .
\end{equation}
% Add something about why we use the expectation values?
We derive the analytical expressions for expectation values in Appendix. 
\ref{sec:expectation_values}. We find the expected total energy
\begin{equation*} %\label{eq:energy_total_first}
	\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_first}
	\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_first}
	\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_first}
	\langle |m| \rangle = \frac{e^{8 \beta J} + 1}{2( \cosh(8 \beta J) + 3)} \ .
\end{equation*}

We will also 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 susceptibility 
\begin{align*}
    \chi = \frac{1}{k_{\text{B} T}} (\mathbb{E}(M^{2}) - [\mathbb{E}(M)]^{2}) \ .
    \label{eq:susceptibility}
\end{align*}

In Appendix \ref{sec:heat_susceptibility} we derive expressions for 
the heat capacity and the susceptibility, given by
\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*}
\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 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}
% P9 critical temperature

\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 depend 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 
require ergodicity and detailed balance. A Markov chain is ergoditic when all system 
states can be reached at every current state, whereas detaild balance implies no 
net flux of probability. To satisfy these criterias we use the Metropolis-Hastings 
algorithm, found in Fig. \ref{algo:metropolis}, to implement the MCMC method.

The Markov process will reach an equilibrium, reflecting the state of a real system.
% Add something about burn-in time
At the point of equilibrium we start sampling, as the distribution of the set of 
samples will tend toward the actual distribution. 
% Add something about spin flipping
At each step of flipping a spin, the change in energy is evaluated as
\begin{align*}
    \Delta E &= E_{\text{after}} - E_{\text{before}} \ .
\end{align*}
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\}$, which we find in Appendix \ref{sec:delta_energy}.
\begin{figure}[H]
    \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}
\end{figure}
We can avoid computing the Boltzmann factor
\begin{equation}
    p(\mathbf{s} | T) = e^{-\beta \Delta E}
    \label{eq:boltzmann_factor}
\end{equation}
at every spin flip, by using a look up table (LUT) with the possible values. We use 
the change in energy $\Delta E$ as a key for the resulting value of the exponential 
function in a hash map.

\subsection{Implementation}\label{subsec:implementation}
% P3 boundary condition and if-tests
To avoid the overhead of if-tests, and take advantage of the parallelization, we 
define an index for every edge case. That is, for a spin at a given boundary index 
we use a pre-set index (...), to avoid if-tests and reduce overhead and runtime.
% P7 parallelization

\subsection{Tools}\label{subsec:tools}
The Ising model and MCMC methods are implemented in C++, and parallelized using 
both \verb|MPI| and \verb|OpenMP| \cite{openmp:2018}. We used the Python library 
\verb|matplotlib| \cite{hunter:2007:matplotlib} to produce all the plots, and 
\verb|seaborn| \cite{waskom:2021:seaborn} to set the theme in the figures. In 
addition, while optimizing our implementation we used the profiler 
% Add profiler

\end{document}