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

\begin{document} 
\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}) &\eqiv \frac{E(\mathbf{s})}{N} \ ,
    \label{eq:energy_spin} \\
    m (\mathbf{s}) &\eqiv \frac{M(\mathbf{s})}{N} \ ,
    \label{eq:magnetization_spin}
\end{align}
for a $N \times N$ lattice.

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}


\subsection{Sampling approach}\label{subsec:sampling}
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

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

\end{document}

