\documentclass[english,notitlepage,reprint,nofootinbib]{revtex4-1}  % defines the basic parameters of the document
% 
% For preview: skriv i terminal: latexmk -pdf -pvc filnavn
% If you want a single-column, remove "reprint"

% Silence warning of revtex4-1
\usepackage{silence}
\WarningFilter{revtex4-1}{Repair the float}

% Allows special characters (including æøå)
\usepackage[utf8]{inputenc}
\usepackage{fontawesome}
\usepackage{multirow}
\usepackage[table]{xcolor}
% \rowcolors{2}{gray!15}{white}
% \usepackage{tabularx}
% \usepackage[english]{babel}

%% Note that you may need to download some of these packages manually, it depends on your setup.
%% I recommend downloading TeXMaker, because it includes a large library of the most common packages.

\usepackage{physics,amssymb}  % mathematical symbols (physics imports amsmath)

\usepackage{graphicx}         % include graphics such as plots
\graphicspath{../images/}
\usepackage{xcolor}           % set colors
\usepackage{hyperref}         % automagic cross-referencing
\usepackage{listings}         % display code
% \usepackage{subfigure}        % imports a lot of cool and useful figure commands
\usepackage{subcaption}
% \usepackage{float}
%\usepackage[section]{placeins}
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{tikz}
% \usetikzlibrary{quantikz}
% defines the color of hyperref objects
% Blending two colors:  blue!80!black  =  80% blue and 20% black
\hypersetup{ % this is just my personal choice, feel free to change things
    colorlinks,
    linkcolor={red!50!black},
    citecolor={blue!50!black},
    urlcolor={blue!80!black}}

% Biblio stuff
% \def\biblio{\bibliographystyle{plain}\bibliography{../references/references}}
\newcommand\numberthis{\addtocounter{equation}{1}\tag{\theequation}}

\usepackage{xr}
\usepackage{subfiles}
% \externaldocument[M-]{\subfix{main}}

\begin{document}

\title{The Ising model}  % self-explanatory
\author{Cory Alexander Balaton \& Janita Ovidie Sandtrøen Willumsen \\ \faGithub \, \url{https://github.uio.no/FYS3150-G2-2023/Project-4}} % self-explanatory
\date{\today}                             % self-explanatory
\noaffiliation                            % ignore this, but keep it.

% Abstract
\begin{abstract}
	We have studied phase transitions in ferromagnets, using the Ising model and ...
\end{abstract}

% \section{Introduction}\label{sec:introduction}
% Magnetic materials are used for their particular property, in a wide range of technological devices. It is used in compasses for navigation, in computer hard drives to store data \cite{yale:2023:hdd}, and in MRI machines to produce detailed anatomical images. Some magnetic materials, called ferromagnets, are easily magnetized. When exposed to a strong magnetic field they become saturated. And heating the ferromagnet to a critical point, will cause it to lose its ferromagnetic property. 

% The Ising model is a mathematical description of a physical system. It was developed to describe ferromagnetism, and to understand the phase transition of a system. However, it has later been used to describe other phenomena where pairwise interactions are of interest, such as networks of neurons \cite{bialek:2008:ising}. 

% We will use the Ising model to study ferromagnet and its temperature-dependent behavior. In addition, we will numerically estimate the critical point where the system experience phase transition. 

% In Section \ref{sec:methods}, we will present the theoretical background for this experiment, as well as the implementation and tools used. Continuing with Section \ref{sec:results}, where we present our results and discuss our findings. Lastly, we will conclude our findings in Section \ref{sec:conclusion}.

% \subsection{Draft}
% We can study individual atoms 
% At a microscopic level we can study individual atoms, using microscopic theory. When a system consist of several atoms, we can study the atoms physical properties of matter using statistical physics. To determine quantities such as heat capacity and susceptibility, we have to use ...

% Some uncharged materials have the ability to acquire attractive powers, that is the atomic spins align as well as the magnetic moment. Neighboring spins interact, and the spins magnetic moment can influence each other. This physical phenomenon is called ferromagnetism \cite{britannica:2023:ferromagnetism}. 

% The Ising model is a mathematical description of a physical system, of a ferromagnet. At a certain teperature, the critical point, the behavior of the ferromagnet changes. We want to study the phase transition of a ferromagnet, and find the critical temperature.

% However, as the lattice size increase, so does the number of calculations. Since computational resources are finite, we use parallelization to decrease runtime and required resources. In parallelizing our implementation, we can approximate the critical temperature of a material.

% In Section \ref{sec:methods}, we will state the theoretical background for this experiment, as well as the implementation and tools used. Continuing with the results in Section \ref{sec:results}, in addition to a discussion of our results. Lastly, we will conclude our findings in Section \ref{sec:conclusion}.

% Study physical processes, phase transition and their behavior near critical points. That is, study how a magnetic material changes properties when exposed to thermal energy. The Ising model was originally a mathematical description meant to simulate ferromagnets. Has found many uses like hemoglobin in absorbing oxygen.

% Ferromagnets exhibit long-range ordering phenomenon at atomic level, line up spins with each other by locking the magnetic moments. Thermal agitation randomizes. Critical temperature of iron is 1043 K. Ising model simply describes how a magnetic material responds to thermal energy and an external magnetic field. The direction of the spins influences the total potential energy.

% Def: Ensemble is a collection of microphysics systems from which we derive expectation values and thermodynamic properties related to experiment. 

% Statistical physics bridge the gap between microscopic and macroscopic world, using microscopic theory we can derive thermodynamic quantities.

\section{Methods}\label{sec:methods}
% \subsection{The Ising Model}\label{subsec:ising_model}
% The Ising model in two dimensions consist of a lattice of spins, and can be though of as atoms in a grid. 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 spin configuration, is given by $\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, resulting in a total system energy given by Eq. \eqref{eq:energy_total}
% \begin{equation}\label{eq:energy_total}
%     E(\mathbf{s}) = -J \sum_{\langle k l \rangle}^{N} s_{k} s_{l} \ .
% \end{equation}
% $\langle k l \rangle$ denote a pair of spins, and $J$ is the coupling constant. In a two-dimensional lattice, a spin can have 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, such that all spins have exactly four nearest neighbors. We count the neighboring spins as visualized in Fig. \ref{fig:tikz_boundary}, also for lattices where $L > 2$.
% \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 pairs of spins in a $2 \times 2$ lattice, where periodic boundary conditions are applied.}
%       \label{fig:tikz_boundary}
% \end{figure}

% In addition, we consider the total magnetization of the system, which is given by Eq. \eqref{eq:magnetization_total}
% \begin{equation}\label{eq:magnetization_total}
%     M(\mathbf{s}) = \sum_{i}^{N} s_{i} \  .
% \end{equation}
% We also consider the state degeneracy, that is the number of different microstates with the same value of total magnetization. To find the analytical expressions necessary for validating our model implementation, we will assume a $2 \times 2$ lattice. Using the possible system states we find the values for the lattice, which can be found in table \ref{tab:lattice_config}.
% \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 will use energy per spin given by Eq. \eqref{eq:energy_spin} and magnetization per spin given by Eq. \eqref{eq:magnetization_spin}, to compare the quantities for different lattice sizes. 
% \begin{align}
%     \epsilon (\mathbf{s}) &= \frac{E(\mathbf{s})}{N} \ ,
%     \label{eq:energy_spin} \\
%     m (\mathbf{s}) &= \frac{M(\mathbf{s})}{N} \ ,
%     \label{eq:magnetization_spin}
% \end{align}

% When we study the phase transitions of ferromagnets, we have to consider the probability of a microstate $\mathbf{s}$, at a fixed temperature $T$, 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$ is 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 also called the partition function, and is the normalizing factor of the pdf. We derive the analytical expression for $Z$ in Appendix \ref{sec:analytical_expressions}, 
% \begin{align*}
%     z &= 4 \cosh (8 \beta J) + 12 \ .
% \end{align*}
% Using this partition function, we can derive analytical expressions for expectation values. For a linear function of a stochastic random variable $X$, with a known probability distribution, the expected value of $x$ is given by
% \begin{align*}
%     \mathbb{E}(aX + b) &= a \cdot \mathbb{E}(X) + b & \text{\cite[p. 131]{springer:2012:modernstat}}
% \end{align*}
% In our case the discrete random variable is the spin configuration, and we want to find the expected value of the function $E(\mathbf{s})$. Continuing, we will use the notation $\langle E \rangle$ for the expectation value of a given variable, in this case $E$. 

% Both energy per spin and magnetization per spin are functions of $\mathbf{s}$. In addition, the number of spins is given as a constant for each lattice. We can use the expression for $\langle E \rangle$ and $\langle M \rangle$ to find the expectation values per spin. For energy per spin 
% \begin{align*}
%     \langle \epsilon \rangle &= \sum_{i=1}^{N} \epsilon(s_{i}) p(s_{i} \ | \ T) \\
%     &= \sum_{i=1}^{N} \frac{E(\mathbf{s})}{N} p(s_{i} \ | \ T) \\
%     &= \frac{1}{N} \sum_{i=1}^{N} E(\mathbf{s}) p(s_{i} \ | \ T)
% \end{align*}
% The same applies to magnetization per spin
% \begin{align*}
%     \langle |m| \rangle = \frac{1}{N} \sum_{i=1}^{N} |M(s_{i})| p(s_{i} \ | \ T)
% \end{align*}
% We find the expected total energy to be
% \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 to be
% \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 determine the heat capacity and the susceptibility of our system, which require an expression of the variance of the energy and magnetization. We can use the definition of variance, given by
% \begin{align*}
%     \mathbb{V}(X) &= \sum_{x \in D} [(x - \mathbb{E}(x))^2 \cdot p(x)] & \text{\cite[p. 132]{springer:2012:modernstat}} \\
%     &= \mathbb{E}(X^{2}) - [\mathbb{E}(X)]^{2} \ ,
% \end{align*}
% for a stochastic variable $X$. We use this to find the heat capacity
% \begin{equation*} %\label{eq:specific_heat_capacity}
% 	C_{V} = \frac{64J^{2} }{N k_{\text{B}} T^{2}} \frac{(3\cosh(8 \beta J) + 1)}{(\cosh(8 \beta J) + 3)^{2}} \ ,
% \end{equation*}
% and the susceptibility
% \begin{equation*} %\label{eq:susceptibility}
% 	\chi = \frac{4}{N k_{\text{B}} T} \frac{(3e^{8 \beta J} + e^{-8 \beta J} + 3)}{(\cosh(8 \beta J) + 3)^{2}}
% \end{equation*}
% The derivation of all analytical expressions can be found in Appendix \ref{sec:analytical_expressions}.

% For scalability, we want to work with unitless spins. To obtain this, we introduce 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}{cccc} % @{\extracolsep{\fill}}
%         \hline
%         Value & Unit \\
%         \hline 
%         $[E]$ & $J$ \\
%         $[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_temp}
At the critical 
temperature the heat capacity $C_{V}$, and the magnetic susceptibility $\chi$ diverge \cite[p. 431]{hj:2015:comp_phys}.
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.
We can describe 
the behavior of the physical system, close to $T_{c}$, 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}
% We will generate samples from the Ising model, to approximate the probability distribution of the sample set. We sample using the Markov chain Monte Carlo 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 to determine statistical quantities \cite{tds:2021:mcmc}.

% The Ising model is initialized in a random state, and the chain is evolved until it reaches an equilibrium. At this point we start sampling, as the distribution of the set of samples will tend toward the actual distribution. At each step of flipping a spin, the change in energy will be evaluated as
% \begin{align*}
%     \Delta E &= E_{\text{after}} - E_{\text{before}} \ .
% \end{align*}
% We have to evaluate 
% The samples generated are microstates, each of which occur with some probability. Since the dependency between samples only apply to the neighboring samples, and the Markov chain is ergodic 

% Since we are generating samples from a stationary probability distribution, with a  in  
% We will use use In addition, we  Since t

% \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}

% \subsection{Implementation}\label{subsec:algo_implementation}
  

% \subsection{Tools}\label{subsec:tools}
% The Ising model and MCMC methods are implemented in C++, and parallelized using \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. 

% \subsection{Draft}
% The Ising model correspond to a lattice with a number of spin sites, where 
% each spin has the possible state up $\uparrow$ or down $\downarrow$. We will study the 
% two-dimensional Ising model for ferromagnets, however, the model is not resticted 
% to this dimentionality~\cite[p. 3]{obermeyer:2020:ising}. We will assume a square lattice, where the length $L$ of 
% the lattice give the number of spin sites $N = L^{2}$. When we
% consider the entire 2D lattice we can study the system spin configuration, or the 
% microstate, given by $\mathbf{s} = (s_{1}, s_{2}, \dots, s_{N})$. The value of 
% each spin $s_{i}$, where $i \in [1, N]$, is given by the direction of the spin 
% $\uparrow = +1$ and $\downarrow = -1$. 

% We find the total energy of the system using
% \begin{equation}\label{eq:energy}
%     E(\mathbf{s}) = -J \sum_{\langle k l \rangle}^{N} s_{k} s_{l}
% \end{equation}
% where $\langle k l \rangle$ denotes the neighboring spins.

% The total magnetization of the system is given by 
% \begin{equation}\label{eq:magnetization}
%     M(\mathbf{s}) = \sum_{i}^{N} s_{i},
% \end{equation}
% with degeneracy denoting how many states have the same values. We will consider 
% a 2D lattice with $L = 2$, with periodic boundary conditions, following the 
% counting pattern shown in figure~\ref{fig:tikz_counting}. The resulting values for
% the $2 \times 2$ lattice can be found in table \ref{tab:lattice_config}.

% \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}
%     \label{tab:lattice_config}
% \end{table}

% Probability stuff, for a system with a given temperature $T$. We use the Boltzmann 
% distribution in the analytical expression, to compare with our results. This is 
% a probability distribution from the family of exponential
% \begin{equation}\label{eq:boltzmann_distribution}
%     p(\mathbf{s} \ | \ T) = \frac{1}{Z} \exp^{-\beta E(\mathbf{s})},
% \end{equation}
% where $Z$ is the partition function $\beta$ is the related to the inverse of the 
% temperature as 
% \begin{equation}\label{eq:beta}
%     \beta = \frac{1}{k_{\text{B}} T}.
% \end{equation} 
% We derive an analytical expression for $Z$ in appendix \ref{sec:analytical_expressions}

% We will assume a $2 \times 2$ lattice to study the possible system states, and 
% to find analytical expressions necessary in the Markov Chain Monte Carlo method.
% These results are used to test our code during implementation.

% We study the 2D lattice and find the total energy and the total magnetization of
% of the system, in addition to the degeneracy. In addition, we want to work with 
% unitless spins. To obtain this, we introduce the base unit for energy $J$, and 
% with the Boltzmann constant we derive other units necessary. The resulting units 
% can be found in table \ref{tab:units}.

% % Problem 1

% \begin{equation}\label{eq:partition_function}
% 	Z = 4 \cosh (8 \beta J) + 12,
% \end{equation}
% %
% \begin{equation}\label{eq:energy_spin_first}
% 	\langle \epsilon \rangle = \frac{-2J \sinh(8 \beta J)}{ \cosh(8 \beta J) + 3}
% \end{equation}
% %
% \begin{equation}\label{eq:energy_spin_second}
% 	\langle \epsilon^{2} \rangle = \frac{4 J^{2} \cosh(8 \beta J)}{\cosh(8 \beta J) + 3}
% \end{equation}
% %
% \begin{equation}\label{eq:magnetization_spin_first}
% 	\langle |m| \rangle = \frac{e^{8 \beta J} + 1}{2( \cosh(8 \beta J) + 3)}
% \end{equation}
% %
% \begin{equation}\label{eq:magnetization_spin_second}
% 	\langle |m|^{2} \rangle = \frac{e^{8 \beta J} + 1}{2( \cosh(8 \beta J) + 3)}
% \end{equation}
% %
% \begin{equation}\label{eq:energy_total_first}
% 	\langle E \rangle = -\frac{8J \sinh(8 \beta J)}{\cosh(8 \beta J) + 3} 
% \end{equation}
% %
% \begin{equation}\label{eq:energy_total_first_squared}
% 	\langle E \rangle^{2} = \frac{64J^{2} \sinh^{2}(8 \beta J)}{(\cosh(8 \beta J) + 3)^{2}}
% \end{equation}
% %
% \begin{equation}\label{eq:energy_total_second}
% 	\langle E^{2} \rangle = \frac{64J^{2} \cosh(8 \beta J)}{\cosh(8 \beta J) + 3}
% \end{equation}
% %
% \begin{equation}\label{eq:magnetization_total_first}
% 	\langle |M| \rangle = \frac{2(e^{8 \beta J} + 2)}{\cosh(8 \beta J) + 3}
% \end{equation}
% %
% \begin{equation}\label{eq:magnetization_total_first_squared}
% 	\langle |M| \rangle^{2} = \frac{4(e^{8 \beta J} + 2)^{2}}{(\cosh(8 \beta J) + 3)^{2}}
% \end{equation}
% %
% \begin{equation}\label{eq:magnetization_total_second}
% 	\langle M^{2} \rangle = \frac{8(e^{8 \beta J} + 1)}{\cosh(8 \beta J) + 3}
% \end{equation}
% %
% \begin{equation}\label{eq:specific_heat_capacity}
% 	C_{V} = \frac{64J^{2} }{N k_{\text{B}} T^{2}} \frac{(3\cosh(8 \beta J) + 1)}{(\cosh(8 \beta J) + 3)^{2}}
% \end{equation}
% %
% \begin{equation}\label{eq:sesceptibility}
% 	\chi = \frac{4}{N k_{\text{B}} T} \frac{(3e^{8 \beta J} + e^{-8 \beta J} + 3)}{(\cosh(8 \beta J) + 3)^{2}}
% \end{equation}
% The derivation of analytical expressions can be found in appendix \ref{sec:analytical_expressions}.


% The Ising model and MCMC methods are implemented in C++, and parallelized using 
% \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. 


% Def: Correlation length $\xi$ of phase transition, and is related to the spin-correlation function which defines the magnetic susceptibility. The spin-correlation function is the covariance of the magnetization, expressing the degree of correlation between spins. $\xi$ is normally finite at critical points. Proportional to the lattice size at the critical point, can be used to relate behavior for infinitly large lattice $T_{C}(L) - T_{C}(L = \infty) = \propto \alpha T^{-1}$

% Implementation: The metropolis algo only consider ratio, not necessary to calculate partition function.

% Metropolis algo: 
% \begin{algorithm}[H]
%     \caption{Metropolis algorithm p. 435 in lecture notes}
%     \label{algo:metropolis_algo}
%     \begin{algorithmic}
%         \Procedure{Metropolis algorithm}{$args$}
%         \State $E_{\text{before}} \leftarrow \text{Set initial random state}$
%         \State $ \text{Flip spin}$ 
%         \State $E_{\text{after}} \text{Calculate energy}$
%         \State $\Delta E \leftarrow E_{\text{after}} - E_{\text{before}}$
%         \If{$\Delta E \leq 0$} Accept new config if
%             \State $\text{Update expectation values}$
%         \ElsIf{$\Delta E > 0$} 
%             \State $\text{Calculate}$
%         \EndIf
%         \EndProcedure
%         \end{algorithmic}
% \end{algorithm}

% Probability: Using the expectation value require a stochastic random variable X (stoch being that it is tied to an event), a function from the "utfallsrommet" to the real numbers.

% Understanding a systems properties close to a critical point. For the 2D ising model, where L is large, we can observe discontinuity in $C_{V}$ and $\chi$. They diverge at the critical point in the thermosynamic limit, meaning the variance in energy and magnetization diverge. For lattice of finite size the variance scale as $1 / \sqrt{M}$ where $M = 2^{N}$ is the number of microstates.
% \begin{align*}
%     e^{-\beta \Delta E} &= \frac{p_{after}}{p_{before}} = \frac{e^{-\beta E_{after}}}{e^{-\beta E_{before}}}
% \end{align*}

\section{Results}\label{sec:results}
% We continue with a lattice where $L = 20$, to study the burn-in time, that is the number of Monte Carlo cycles necessary for the system to reach an equilibrium. We consider two different temperatures $T_{1} = 1.0 J/k_{B}$ and $T_{2} = 2.4 J/k_{B}$, where $T_{2}$ is close to the critical temperature. We can use the correlation time $\tau \approx L^{d + z}$ to determine time, where $d$ is the dimensionality of the system and $z = 2.1665 \pm 0.0012$ \footnote{This value was determined by Nightingale and Blöte for the Metropolis algorithm.}
% % Need to include a section of Onsager's analytical result.
% We show the numerical estimates for temperature $T_{1}$ of $\langle \epsilon \rangle$ in Figure \ref{fig:burn_in_energy_1_0} and $\langle |m| \rangle$ in Figure \ref{fig:burn_in_magnetization_1_0}. For temperature $T_{2}$, the numercal estimate of $\langle \epsilon \rangle$ is shown in Figure \ref{fig:burn_in_energy_2_4} and $\langle |m| \rangle$ in Figure \ref{fig:burn_in_magnetization_2_4}. The lattice is initialized in both an ordered and an unordered state. We observe that for $T_{1}$ there is no change in either expectation value with increasing number of Monte Carlo cycles, when we start with an ordered state. As for the unordered lattice, we observe a change for the first 5000 MC cycles, where it stabilizes. The approximated expected energy is $-2$ and expected magnetization is $1.0$, which is to be expected for temperature 0f $1.0$. T is below the critical and the pdf using $T = 1.0$ result in
% \begin{align*}
%     p(s|T=1.0) &= \frac{1}{e^{-\beta \sum E(s)}} e^{-\beta E(s)} \\
%     &= \frac{1}{e^{-(1/k_{B}) \sum E(s)}} e^{-(1/k_{B}) E(s)}
% \end{align*}
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/burn_in_time_magnetization_1_0.pdf}
%     \caption{$\langle |m| \rangle$ as a function of time, for $T = 1.0 J / k_{B}$}
%     \label{fig:burn_in_magnetization_1_0}
% \end{figure}
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/burn_in_time_energy_1_0.pdf}
%     \caption{$\langle \epsilon \rangle$ as a function of time, for $T = 1.0 J / k_{B}$}
%     \label{fig:burn_in_energy_1_0}
% \end{figure}
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/burn_in_time_magnetization_2_4.pdf}
%     \caption{$\langle |m| \rangle$ as a function of time, for $T = 2.4 J / k_{B}$}
%     \label{fig:burn_in_magnetization_2_4}
% \end{figure}
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/burn_in_time_energy_2_4.pdf}
%     \caption{$\langle \epsilon \rangle$ as a function of time, for $T = 2.4 J / k_{B}$}
%     \label{fig:burn_in_energy_2_4}
% \end{figure}


% 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.

% We use the estimated burn-in time to set starting time for sampling, then generate samples to plot in a histogram for $T_{1}$ in Figure \ref{fig:histogram_1_0} and $T_{2}$ in Figure \ref{fig:histogram_2_4}. For $T_{1}$ we can see that most samples have the expected value $-2$, we have a distribution with low variance.

% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/pd_estimate_1_0.pdf}
%     \caption{Histogram $T = 1.0 J / k_{B}$}
%     \label{fig:histogram_1_0}
% \end{figure}
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/pd_estimate_2_4.pdf}
%     \caption{Histogram $T = 2.4 J / k_{B}$}
%     \label{fig:histogram_2_4}
% \end{figure}

For $T_{2}$ the samples are distributed across a wider range of energy values, 
the expected value is at approx -1.2, corresponding to what we observed the value 
at equilibrium in Figure \ref{fig:burn_in_energy_2_4}. In addition, the distribution 
has a higher variance for temperature close to critical. We also observe that the 
shape of sample distribution closely resemble the normal distribution.

% For a $L \times L$ lattice with $N = L^2$ spins and $M = 2^{N}$ microstates, 
% the variance scale as Var$(L) = \frac{1}{\sqrt{2^{L \times L}}}$.

% The spin correlation function = covariance, degree of correlation between spins.

% We started parallel implementation using OpenMP, limiting number of threads to 
% ensure the program did not hog the entire cpu. We made sure the workload was 
% divided equally between workers (threads) to ensure balanced workload. In addition, 
% we created a look-up-table for values of energy shift to limit repeated calculations 
% and parallel overhead of if-statements. We limit the number of times threads are 
% spawned and joined, by using single parallel regions, reducing parallel overhead. 

% To improve the efficiancy we divided the range of temperatures between MPI processes, 
% parallelizing the simulation further. Each mpi process then spawn a number of threads 
% which initialize their own Ising model, running through burn-in time and generate 
% samples.

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 parallelize using MPI. We generated samples for the temperature range $T \in [2.1, 2.4]$. Using Fox we generated both 1 million samples and 10 million samples.

The physical behavior of physical quantities like magnetization, heat capasity etc. near the critical temperature can be characterized using power law (p.431 lecture).

We observe an ordered phase at low temperatures, meaning we have ferromagnetism. The Metropolis algorithm near critical temperature

The energy is given by the first derivative of the thermodynamic potental
\begin{align*}
    F &= \langle E \rangle - TS = -kT \log Z & \text{$log Z = -F/kT = -F\beta$}
\end{align*}
Specific heat is the second derivative

Timing the code using both serial, OpenMP, and MPI+OpenMP we observe a speed up in Table \ref{tab:timing}. The speed-up of runtime is found using strong-scaling 
\begin{align*}
    time &= burn-in + parallel/workers
\end{align*}
We are not dividing the work of simulating a single Ising model between threads, as there is dependency between the steps. We can use Amdahl's law \cite[p. 124]{hager:2011:hpc} to calculate application speed-up, which simplifies as
\begin{align*}
S_f &= \frac{s + p}{s + p/N} \\
    app speedup &= \frac{serial runtime}{parallel runtime}
\end{align*}

When the speed-up was satisfactory, we investigated the phase transition for lattices of size $L = 20, 40, 60, 80, 100$. The result of $\langle \epsilon \rangle$ can be found in Figure \ref{fig:phase_energy}, $\langle |m| \rangle$ in Figure \ref{fig:phase_magnetization}, $C_{V}$ in Figure \ref{fig:phase_heat}, and $\chi$ in Figure \ref{fig:phase_susceptibility}. 
% Add observation around critical temp
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/1M/energy.pdf}
%     \caption{$\langle \epsilon \rangle$ for $T \in [2.1, 2.4]$, 1000000 MC cycles.}
%     \label{fig:phase_energy}
% \end{figure}

% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/1M/magnetization.pdf}
%     \caption{$\langle |m| \rangle$ for $T \in [2.1, 2.4]$, 1000000 MC cycles.}
%     \label{fig:phase_magnetization}
% \end{figure}

% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/1M/heat_capacity.pdf}
%     \caption{$C_{V}$ for $T \in [2.1, 2.4]$, 1000000 MC cycles.}
%     \label{fig:phase_heat}
% \end{figure}

% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/1M/susceptibility.pdf}
%     \caption{$\chi$ for $T \in [2.1, 2.4]$, 1000000 MC cycles.}
%     \label{fig:phase_susceptibility}
% \end{figure}

% We include results using 10 million MC cycles in Appendix \ref{sec:extra}

% We use the critical temperatures found studying the phase transition, in addition to the scaling relation
% \begin{equation}
%     T_{c} - T_{c}(L = \infty) = \alpha L^{-1}
% \end{equation}
% to estimate the critical temperatur for a lattize of infinte size. We also compare the estimate with the analytical solution 
% \begin{equation}
%     T_{c}(L = \infty) = \frac{2}{\ln (1 + \sqrt{2})} J/k_{B} \approx 2.269 J/k_{B}
% \end{equation}
% using linear regression. In Figure \ref{fig:linreg} we show the critical temperatures as function of the inverse lattice size. As size goes toward infinity we reach inv L = 0 and find the intercept which is equal to the estimated critical of L=infty.
 
which is close to the analytical solution $T_{C}(L = \infty) \approx 2.269 J/k_{B}$ 
found by Lars Onsager. 
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/1M/linreg.pdf}
%     \caption{Linear regression, where $\beta_{0}$ is the intercept $T_{c}(L = \infty)$ and $\beta_{1}$ is the slope.}
%     \label{fig:linreg}
% \end{figure}


%\subsection{Draft}

% Test that the numerical stuff gets close to the analytical 
% - validate implementation a given number of times, find average number of cycles

% Initial state ordered vs random
% - start ordered both at up=1 and random 

% burn-in time
% - burn-in 5000 vs 10000 mc cycles, only mention how a longer burn-in time can in some cases result in better cv and x

% probability using histogram of samples
% - Comparing the histogram of $T_1 = 1.0$ and $T_2 = 2.4$. For $T_1 = 1.0$ we see a lower acceptance of flips, resulting in a low variance. The histogram looks more like an exponential function

% parallelization
% - temp mpi
% - each mpi process spawn thread that make own ising model
% - each ising model goes through burn in time (no benefit)
% - sampling done in parallel
% - mcmc sampling openmp

% phase transition
% critical temperature

\onecolumngrid

\bibliographystyle{unsrtnat}
\bibliography{ref}

\appendix
% \section{Ising model system states}\label{sec:system_states}
% Units 
% \begin{table}[H]
%     \centering
%     \begin{tabular}[c]{cc} % @{\extracolsep{\fill}}
%         \hline
%         Value & Unit \\
%         \hline 
%         $[ E ]$ & $J$ \\
%         $[ T ]$ & $J / k_{\text{B}}$ \\ 
%         $[ M ]$ & $\dots$ \\
%         $[ C_{\text{V}} ]$ & $k_{\text{B}}$ \\
%         $[ \chi ]$ & $1 / J$ \\
%         \hline
%     \end{tabular}
%     \caption{Units, given by the coupling constant $J$ and the Boltzmann constant $k_{\text{B}}$.}
%     \label{tab:units}
% \end{table}

% To avoid counting duplicates, we used 
% \begin{figure}\label{fig:tikz_counting}
%     \centering
%     \begin{tikzpicture}
%         \draw (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{Rules for multiplying spin pairs.}
% \end{figure}

% \begin{figure}\label{fig:tikz_neighbor}
%     \centering
%     \begin{subfigure}{0.4\linewidth}
%         \begin{tikzpicture}
%             \draw (0, 0) grid (2, 2);
%             \node (s1) at (0.5, 1.5) {$\uparrow$};
%             \node (s2) at (1.5, 1.5) {$\uparrow$};
%             \node (s3) at (0.5, 0.5) {$\downarrow$};
%             \node (s4) at (1.5, 0.5) {$\downarrow$};
%         \end{tikzpicture}
%         \caption{}
%         \label{fig:sub_tikz_neighbor_a}
%     \end{subfigure}
%     \
%     \begin{subfigure}{0.4\linewidth}
%         \begin{tikzpicture}
%             \draw (0, 0) grid (2, 2);
%             \node (s1) at (0.5, 1.5) {$\uparrow$};
%             \node (s2) at (1.5, 1.5) {$\downarrow$};
%             \node (s3) at (0.5, 0.5) {$\downarrow$};
%             \node (s4) at (1.5, 0.5) {$\uparrow$};
%         \end{tikzpicture}
%         \caption{}
%         \label{fig:sub_tikz_neighbor_b}
%     \end{subfigure}
%     \caption{Possible spin configurations for two spins up.}
% \end{figure}

%----------
% APPENDIX
%----------
% \section{Analytical expressions}\label{sec:analytical_expressions}
% The Boltzmann distribution is normalized using a partition function $Z$, given by
% \begin{equation}\label{eq:partition}
%     Z = \sum_{\text{all possible } \mathbf{s}}^{N}
% \end{equation}
% We use the values estimated for the $2 \times 2$ case, found in \ref{tab:lattice_config}, 
% and find the partition function
% \begin{align*}
% 	Z &= 1 \cdot e^{-\beta (-8J)} + 4 \cdot e^{-\beta (0)} + 4 \cdot e^{-\beta (0)} + 2 \cdot e^{-\beta (8J)} \\
%     & \quad + 4 \cdot e^{-\beta (0)} 1 \cdot e^{-\beta (-8J)} \\
%     &= 2e^{8 \beta J} + 2e^{-8 \beta J} + 12.
% \end{align*}
% We rewrite the expression using $\cosh(8 \beta J) = 1/2 \big( e^{8 \beta J} + e^{-8 \beta J})$, and get
% \begin{align*}
%     z &= 4 \cosh (8 \beta J) + 12 
% \end{align*}
% The Boltzmann distribution is given by 
% \begin{equation}\label{eq:boltzmann}
%     p(\mathbf{s} \ | \ T) = \frac{1}{Z} e^{-\beta E(\mathbf{s})},
% \end{equation}
% for a given temperature $T$. With our expression for the partition function, we 
% get the probability distribution
% \begin{align*}
%     p(\mathbf{s} \ | \ T) &= \frac{1}{4 \cosh (8 \beta J) + 12} e^{-\beta E(\mathbf{s})}
% \end{align*}

% For discrete random variables $X$, with a known probability distribution, the 
% expected value of $x$ is given by 
% \begin{align*}
%     \mathbb{E}(x) &= \sum_{x \in D} x \cdot p(x) & \text{\cite[p. 127]{springer:2012:modernstat}}.
% \end{align*}
% For a function of a stochastic random variable, the expected value of $x$ is
% \begin{align*}
%     \mathbb{E}(h(X)) &= \sum_{x \in D} h(x) \cdot p(x) 
% \end{align*}
% And in the case of a linear function we have
% \begin{align*}
%     \mathbb{E}(aX + b) &= a \cdot \mathbb{E}(X) + b & \text{\cite[p. 131]{springer:2012:modernstat}}
% \end{align*}
% In our case the discrete random variable is the spin configuration, and we want 
% to find the expected value of the function $E(\mathbf{s})$. In addition, we use 
% $\langle E \rangle$ as notation for expexted value of a given variable, since $\mathbf{s}$ is the stochastic
% random variable.

% The expression for total energy and total magnetization is given in eq. \eqref{eq:energy} and \eqref{eq:magnetization}. 
% The expected values for these is given by
% \begin{equation*}
%     \langle E(\mathbf{s}) \rangle = \sum_{i=1}^{N} E(s_{i}) p(s_{i} \ | \ T)
% \end{equation*}

% \begin{equation*}
%     \langle |M(\mathbf{s})| \rangle = \sum_{i=1}^{N} |M(s_{i})| p(s_{i} \ | \ T)
% \end{equation*}

% Since we want to compare expected values for different lattice sizes, we have to 
% find the expected values per spin. We normalize the total expressions for total 
% energy \eqref{eq:energy} and magnetizaation \eqref{eq:magnetization} by the 
% number of spins to get 
% \begin{equation}\label{eq:spin_energy}
%     \epsilon(\mathbf{s}) = \frac{E(\mathbf{s})}{N}
% \end{equation}

% \begin{equation}\label{eq:spin_magnetization}
%     m(\mathbf{s}) = \frac{M(\mathbf{s})}{N}
% \end{equation}

% Both energy per spin and magnetization per spin are functions of $\mathbf{s}$, 
% we will use the short hand notation $\langle \epsilon \rangle$ and $\langle |m| \rangle$. 
% In addition, the number of spins are given as a constant for each lattice. We can 
% rewrite and use this when we find the expectation values per spin, for energy per spin 
% \begin{align*}
%     \langle \epsilon \rangle &= \sum_{i=1}^{N} \epsilon(s_{i}) p(s_{i} \ | \ T) \\
%     &= \sum_{i=1}^{N} \frac{E(\mathbf{s})}{N} p(s_{i} \ | \ T) \\
%     &= \frac{1}{N} \sum_{i=1}^{N} E(\mathbf{s}) p(s_{i} \ | \ T)
% \end{align*}
% The same applies for magnetization per spin
% \begin{align*}
%     \langle |m| \rangle = \frac{1}{N} \sum_{i=1}^{N} |M(s_{i})| p(s_{i} \ | \ T)
% \end{align*}
% Continuing with the expectation values for a $2 \times 2$ lattice, excluding the terms which give zero we get
% \begin{align*}
%     \langle E \rangle &= (-8J) \cdot \frac{1}{Z} e^{8 \beta J} + 2 \cdot (8J) \cdot \frac{1}{Z} e^{-8 \beta J} + (-8J) \cdot \frac{1}{Z} e^{8 \beta J} \\
%     &= \frac{16J}{Z} \big(e^{-8 \beta J} - e^{8 \beta J}) \\
%     &= -\frac{32J \sinh(8 \beta J)}{4(\cosh(8 \beta J) + 3)} \\
%     &= -\frac{8J \sinh(8 \beta J)}{\cosh(8 \beta J) + 3},
% \end{align*}
% and
% \begin{align*}
%     \langle |M| \rangle &= 4 \cdot \frac{1}{Z} \cdot e^{8 \beta J} + 4 \cdot 2 \cdot \frac{1}{Z} \cdot e^{0} \\
%     & \quad + 4 \cdot | -2| \cdot \frac{1}{Z} \cdot e^{0} + | -4| \cdot e^{8 \beta J} \\
%     &= \frac{8 e^{8 \beta J} + 16}{Z} \\
%     &= \frac{4 (2e^{8 \beta J} + 4)}{4(\cosh(8 \beta J) + 3)} \\
%     &= \frac{2(e^{8 \beta J} + 2)}{\cosh(8 \beta J) + 3}
% \end{align*}
% The squared function 
% \begin{align*}
%     \langle E^{2} \rangle &= (-8J)^{2} \cdot \frac{1}{Z} e^{8 \beta J} + 2 \cdot (8J)^{2} \cdot \frac{1}{Z} e^{-8 \beta J} + (-8J)^{2} \cdot \frac{1}{Z} e^{8 \beta J} \\
%     &= \frac{128J^{2}}{Z} \big(e^{8 \beta J} + e^{-8 \beta J} \big) \\
%     &= \frac{128J^{2} \cosh(8 \beta J)}{4(\cosh(8 \beta J) + 3)} \\
%     &= \frac{64J^{2} \cosh(8 \beta J)}{\cosh(8 \beta J) + 3},
% \end{align*}
% and
% \begin{align*}
%     \langle M^{2} \rangle &= 4^{2} \cdot \frac{1}{Z} \cdot e^{8 \beta J} + 4 \cdot 2^{2} \cdot \frac{1}{Z} \cdot e^{0} \\
%     & \quad + 4 \cdot (-2)^{2} \cdot \frac{1}{Z} \cdot e^{0} + (-4)^{2}\cdot e^{8 \beta J} \\
%     &= \frac{32e^{8 \beta J} + 32}{Z} \\
%     &= \frac{4 (8e^{8 \beta J} + 8)}{4(\cosh(8 \beta J) + 3)} \\
%     &= \frac{8e^{8 \beta J} + 8}{\cosh(8 \beta J) + 3} 
% \end{align*}
% The squared expectation value is given by 
% \begin{align*}
%     \langle E \rangle^{2} &= \bigg(-\frac{8J \sinh(8 \beta J)}{\cosh(8 \beta J) + 3} \bigg)^{2} \\
%     &= \frac{64J^{2} \sinh^{2}(8 \beta J)}{(\cosh(8 \beta J) + 3)^{2}},
% \end{align*}
% and
% \begin{align*}
%     \langle |M| \rangle^{2} &= \Big( \frac{2(e^{8 \beta J} + 2)}{\cosh(8 \beta J) + 3} \Big)^{2} \\
%     &= \frac{4(e^{8 \beta J} + 2)^{2}}{(\cosh(8 \beta J) + 3)^{2}}
% \end{align*}

% Calculating the heat capacity and susceptibility, we need the variance of both total 
% energy and total magnetizaation. We obtain this using the definition 
% \begin{align*}
%     \mathbb{V}(X) &= \sum_{x \in D} [(x - \mathbb{E}(x))^2 \cdot p(x)] & \text{\cite[p. 132]{springer:2012:modernstat}}. \\
%     &= \mathbb{E}(X^{2}) - [\mathbb{E}(X)]^{2}
% \end{align*}

% The variance of total energy is given by
% \begin{align*}
%     \mathbb{V}(E) &= \mathbb{E}(E^{2}) - [\mathbb{E}(E)]^{2} \\
%     &= \frac{64J^{2} \cosh(8 \beta J)}{\cosh(8 \beta J) + 3} - \frac{64J^{2} \sinh^{2}(8 \beta J)}{(\cosh(8 \beta J) + 3)^{2}} \\
%     &= 64J^{2} \bigg( \frac{\cosh(8 \beta J)}{\cosh(8 \beta J) + 3} - \frac{\sinh^{2}(8 \beta J)}{(\cosh(8 \beta J) + 3)^{2}} \bigg) \\
%     &= 64J^{2} \bigg( \frac{(\cosh(8 \beta J)) \cdot (\cosh(8 \beta J) + 3)}{(\cosh(8 \beta J) + 3)^{2}} \\
%     & \quad - \frac{\sinh^{2}(8 \beta J)}{(\cosh(8 \beta J) + 3)^{2}} \bigg) \\
%     &= 64J^{2} \bigg( \frac{\cosh^{2}(8 \beta J) + 3\cosh(8 \beta J) - \sinh^{2}(8 \beta J)}{(\cosh(8 \beta J) + 3)^{2}} \bigg) \\
%     &= 64J^{2} \bigg( \frac{3\cosh(8 \beta J) + 1}{(\cosh(8 \beta J) + 3)^{2}} \bigg)
% \end{align*}
% \begin{align*}
%     \mathbb{V}(M) &= \mathbb{E}(M^{2}) - [\mathbb{E}(|M|)]^{2} \\
%     &= \frac{8e^{8 \beta J} + 8}{\cosh(8 \beta J) + 3} - \frac{4(e^{8 \beta J} + 2)^{2}}{(\cosh(8 \beta J) + 3)^{2}} \\
%     &= \frac{(8(e^{8 \beta J} + 1)) \cdot (\cosh(8 \beta J) + 3) - 4(e^{8 \beta J} + 2)^{2}}{(\cosh(8 \beta J) + 3)^{2}} \\
%     &= \frac{4(e^{8 \beta J} + 1) \cdot (e^{8 \beta J} + e^{-8 \beta J}) + 24(e^{8 \beta J} + 1) - 4(e^{8 \beta J} + 1)^{2}}{(\cosh(8 \beta J) + 3)^{2}} \\
%     &= \frac{4e^{2(8 \beta J)} + 4e^{8 \beta J} 4e^{0} + 4e^{-8 \beta J} 24e^{8 \beta J} + 24 - 4e^{2(8 \beta J)} - 16e^{8 \beta J} - 16}{(\cosh(8 \beta J) + 3)^{2}} \\
%     &= \frac{4(3e^{8 \beta J} + e^{-8 \beta J} + 3)}{(\cosh(8 \beta J) + 3)^{2}}
% \end{align*}

% \begin{align*}
%     C_{\text{V}} &= \frac{1}{N} \frac{1}{k_{\text{B} T^{2}}} (\mathbb{E}(E^{2}) - [\mathbb{E}(E)]^{2}) \\
%     &= \frac{1}{N k_{\text{B} T^{2}}} \mathbb{V}(E) \\
%     &= \frac{64J^{2} }{N k_{\text{B}} T^{2}} \bigg( \frac{3\cosh(8 \beta J) + 1}{(\cosh(8 \beta J) + 3)^{2}} \bigg)
% \end{align*}

% \begin{align*}
%     \chi &= \frac{1}{N} \frac{1}{k_{\text{B} T}} (\mathbb{E}(M^{2}) - [\mathbb{E}(M)]^{2}) \\
%     &= \frac{1}{N k_{\text{B} T}} \mathbb{V}(M) \\
%     &= \frac{4}{N k_{\text{B} T}} \bigg( \frac{3e^{8 \beta J} + e^{-8 \beta J} + 3}{(\cosh(8 \beta J) + 3)^{2}} \bigg)
% \end{align*}

% \begin{align*}
%     \langle \epsilon^{2} \rangle &= \frac{1}{N^{2}} \sum_{i=1}^{N} E(\mathbf{s})^{2} p(s_{i} \ | \ T) \\
%     &=
% \end{align*}
% The same applies for magnetization per spin
% \begin{align*}
%     \langle |m| \rangle = \frac{1}{N} \sum_{i=1}^{N} |M(s_{i})| p(s_{i} \ | \ T)
% \end{align*}

% \section{Extra}\label{sec:extra}
% We increased number of MC cycles to 10 million
% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/10M/energy.pdf}
%     \caption{$\langle \epsilon \rangle$ for $T \in [2.1, 2.4]$, 10000000 MC cycles.}
%     \label{fig:phase_energy_10M}
% \end{figure}

% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/10M/magnetization.pdf}
%     \caption{$\langle |m| \rangle$ for $T \in [2.1, 2.4]$, 10000000 MC cycles.}
%     \label{fig:phase_magnetization_10M}
% \end{figure}

% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/10M/heat_capacity.pdf}
%     \caption{$C_{V}$ for $T \in [2.1, 2.4]$, 10000000 MC cycles.}
%     \label{fig:phase_heat_10M}
% \end{figure}

% \begin{figure}
%     \centering
%     \includegraphics[width=\linewidth]{../images/phase_transition/fox/wide/10M/susceptibility.pdf}
%     \caption{$\chi$ for $T \in [2.1, 2.4]$, 10000000 MC cycles.}
%     \label{fig:phase_susceptibility_10M}
% \end{figure}


\end{document}