\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 Figure \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. In the case where we 
have two spins oriented up the total energy have two possible values, as shown 
in \ref{sec:energy_special}.
% 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$ 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 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})} \ .
    \label{eq:pdf}
\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_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 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 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}
$\boldsymbol{Draft}$
% P9 critical temperature
% Add something about critical temperature in ferromagnets
At temperatures below the critical temperature $T_{c}$, the Ising model will 
magnetize spontaneous. Based on a $2 \times 2$ lattice, we can show that the total 
energy is equal to the energy where all spins have the orientation up \cite[p. 426]{hj:2015:comp_phys}.

When a ferromagnetic material is heated, it will change at a macroscopic level.
When increasing the temperature of the external field, the Ising model move from 
an ordered to an unordered phase. At the critical temperature the heat capacity 
$C_{V}$, and the magnetic susceptibility $\chi$ diverge \cite[p. 431]{hj:2015:comp_phys}.

We can describe the behavior of the physical system close to the critical temperature 
using power laws and critical exponents. For an Ising model of infinite lattice size 
in two dimensions we have 
\begin{align}
    \langle |m| \rangle &\propto |T - T_{c}(L = \infty)|^{\beta} \\
    C_{V} &\propto |T - T_{c}(L = \infty)|^{-\alpha} \\ 
    \chi &\propto |T - T_{c}(L = \infty)|^{-\gamma}
\end{align}

\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 Figure \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\}$, we derive these values 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 and testing}\label{subsec:implementation_test}
$\boldsymbol{Draft}$
% 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.
% P4 testing, validation
% P7 parallelization
10 MPI processes
- 10 threads per process
= 100 threads total
First we divide the temperature range, to give each MPI process a set of temperatures 
to work with. 2.1-2-4 divided into 40 steps, which gives us a step size of 0.0075.

Not a lot of downtime for the threads

However, when the temperature is close to the critical point, we observe an increase 
in expected energy and a decrease in magnetization. Suggesting a higher energy and 
a loss of magnetization close to the critical temperature.

% We did not set the seed for the random number generator, which resulted in 
% different numerical estimates each time we ran the model. However, all expectation 
% values are calculated using the same data. The burn-in time varied each time. 
% We see a burn-in time t = 5000-10000 MC cycles. However, this changed between runs. 
We decided with a burn-in time parallelization trade-off. That is, we set the 
burn-in time lower in favor of sampling. To take advantage of the parallelization 
and not to waste computational resources. The argument to discard samples generated 
during the burn-in time is ... Increasing number of samples outweigh the ...

It is worth mentioning that the time (number of MC cycles) necessary to get a 
good numerical estimate, compared to the analytical result, foreshadowing the 
burn-in time.

Markov chain starting point can differ, resulting in different simulation. By 
discarding the first samples, the ones generated before system equilibrium we can 
get an estimate closer to the real solution. Since we want to estimate expectation 
values at a given temperature, the samples should represent the system at that 
temperature.

Depending on number of samples used in numerical estimates, using the samples 
generated during burn-in can in high bias and high variance if the ratio is skewed. 
However, if most samples are generated after burn-in the effect is not as visible. 
Can't remove randomness by starting around equilibrium, since samples are generated 
using several ising models we need to sample using the same conditions that is 
system state equilibrium.


\subsection{Tools}\label{subsec:tools}
The Ising model and MCMC methods are implemented in C++, and parallelized using 
both \verb|OpenMPI| \cite[text]{gabriel:2004:open_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, we used \verb|Scalasca| \cite{scalasca} and \verb|Score-P| \cite{scorep} 
to profile and optimize our implementation.
% Add profiler

\end{document}