\subsection{Ising Model}
In a two-dimensional system of size $L \times L$, each lattice site $i$ contains a spin variable $s_i$ that can take on values of either +1 (up) or -1 (down). The energy of a given configuration of spins is described by the sum over nearest-neighbor interactions \cite{anderskvellestadProject4FYS3150}:
\begin{equation}
    E = -J \sum_{\langle i,j \rangle} s_i s_j.
\end{equation}

The coupling constant, $J$, is the parameter that quantifies the strength of the interaction between neighboring spins. The sum is taken over all pairs of nearest neighbors, defined as $\langle i,j \rangle$. In order to facilitate the process of summation without the occurrence of double counting, the total energy can be expressed in the following manner:
\begin{equation} \label{eq:total_energy}
    E = -\frac{J}{2} \sum_{i} \sum_{j \in \text{NN}(i)}  s_i s_j,
\end{equation}
where $\text{NN}(i)$ represents the set of nearest neighbors of site $i$. Depending on the configuration of the neighboring spins, the flipping of a single spin, $s_k$, will result in a change in the energy of the system by an amount
\begin{equation}
    \Delta E = \Delta s_k \cdot J \sum_{j \in \text{NN}(k)} s_j,
\end{equation}
where $\Delta s_k = s_k^{\text{new}} - s_k^{\text{old}} \in \{-2, 2\}$ is the change in the spin variable at site $k$.  To substantiate this claim, it is essential to recognize that the energy contribution from site $k$ prior to the flip is given by:
\begin{equation}
    E_{\text{old}} = -J s_k^{\text{old}} \sum_{j \in \text{NN}(k)} s_j.
\end{equation}
The discrepancy in the factor 2 in comparison to \autoref{eq:total_energy} emerges due to the fact that $s_k$ is a constituent of the sum over nearest neighbors for each of its neighbors, resulting in double counting. Subsequent to the flipping of the spin at site $k$, the contribution to the energy from that site is updated, and after the flip, it becomes
\begin{equation}
    E_{\text{new}} = -J s_k^{\text{new}} \sum_{j \in \text{NN}(k)} s_j.
\end{equation}
The alteration in energy resulting from the flip is then
\begin{equation}
    \Delta E = E_{\text{new}} - E_{\text{old}} = -J (s_k^{\text{new}} - s_k^{\text{old}}) \sum_{j \in \text{NN}(k)} s_j.
\end{equation}
This expression demonstrates how the energy change depends exclusively on the local configuration of spins surrounding site $k$.

The total magnetization of the system $M$ is derived from the summation of all spin variables \cite{anderskvellestadProject4FYS3150}
\begin{equation}
    M = \sum_{i} s_i.
\end{equation}
Both extensive quantities, defined as those that scale with the system size, can be normalized per spin to obtain intensive quantities
\begin{equation}
    \epsilon = \frac{E}{N}, \quad m = \frac{M}{N},
\end{equation}
where $N = L^2$ is the total number of spins in the lattice.

The Ising model displays a specific heat capacity, $C_V$, and magnetic susceptibility, $\chi$, both of which are defined as follows, once more normalized per spin \cite{anderskvellestadProject4FYS3150}:
\begin{equation} \label{eq:heat_capacity}
    \frac{C_V}{N} = \frac{N \left(\expect{\epsilon^2} - \expect{\epsilon}^2\right)}{k_B T^2},
\end{equation}
\begin{equation} \label{eq:susceptibility}
    \frac{\chi}{N} = \frac{N \left(\expect{m^2} - \expect{m}^2\right)}{k_B T},
\end{equation}
where $k_B$ is the Boltzmann constant and $T$ is the temperature of the system. The notation $\expect{\cdot}$ denotes the expectation value, which can be estimated through averaging over sampled configurations. The average over a sufficiently large number of Monte Carlo samples will converge to the true expectation value. Given that expectation values are computed from sampled configurations, it is imperative to ensure that the system's state space is sampled correctly. In the Ising model, it is assumed that the configurations follow the Boltzmann distribution. This distribution dictates the probability of a configuration with energy $E$ at temperature $T$ as follows \cite{anderskvellestadProject4FYS3150,kulskeIsingModelHighlights2025}:
\begin{equation}
    P(E) = \frac{e^{-E/(k_B T)}}{Z},
\end{equation}
with the partition function $Z$ defined as
\begin{equation}
    Z = \sum_{\text{all states}} e^{-E/(k_B T)}.
\end{equation}

\subsubsection{Critical Temperature}
In two dimensions, the Ising model undergoes a continuous phase transition at the critical temperature $T_C$ separating an ordered low-temperature phase from a disordered high-temperature phase. In the context of the model, spontaneous magnetization is exhibited when the temperature falls below the critical temperature, $T_C$. This phenomenon occurs because, even in the absence of an external magnetic field, the magnetization $m$ remains non-zero. As the temperature approaches $T_C$ from below, the magnetization decreases continuously, eventually vanishing according to a characteristic power law. Concurrently, fluctuations in energy and magnetization exhibit a marked increase, as evidenced by the divergent behavior of both the heat capacity \autoref{eq:heat_capacity}, and the susceptibility \autoref{eq:susceptibility}. Precisely at $T_C$, the system becomes scale-invariant, and configurations manifest correlations that extend across the entire lattice. In numerical simulations, such phenomena manifest as strong critical slowing down, necessitating sufficiently long sampling times to resolve the sharply peaked response functions. The exact solution for the square lattice with nearest-neighbor interactions yields the following result \cite{anderskvellestadProject4FYS3150,metraTemperaturedependentCriticalityRandom2021}:
$k_B T_C / J = 2 / \ln(1 + \sqrt{2}) \approx \num{2.269}$, providing a precise reference point for validating numerical results.

\subsubsection{Ergodic Property}
The implementation of Monte Carlo methods within the Ising model is contingent upon the assumption of ergodicity, which stipulates that the Markov chain generated by successive spin updates must possess the capacity to attain all physically relevant configurations within a sufficient time frame. In practice, this implies that single-spin flip dynamics must form a connected set over the configuration space, thereby avoiding traps in metastable regions. The Metropolis–Hastings update rule, which is employed for sampling purposes, ensures detailed balance by ensuring that the stationary distribution of the Markov chain corresponds to the Boltzmann distribution. Ergodicity guarantees the convergence of long-time averages of observables, such as  $\epsilon$ and $m$, to ensemble averages. As the critical point is approached, ergodicity remains valid; however, the time required to decorrelate successive configurations grows significantly due to the emergence of large correlated domains. This phenomenon, often referred to as critical slowing down, necessitates meticulous care in determining autocorrelation times to obtain statistically independent samples and reliable estimates of thermodynamic quantities.
When an ergodic Markov chain is considered, the equivalence between ensemble averages and long-time averages is guaranteed in the limit of a sufficient number of Monte Carlo steps. For any observable $A$, the ensemble expectation value $\expect{A}$ obtained from the Boltzmann distribution coincides with the time-averaged estimate $\hat{A}$ extracted from the simulation. It is imperative to note that the aforementioned equality is only valid once the chain has thoroughly explored configuration space. This underscores the significance of ergodicity and adequate sampling length as essential prerequisites for reliable numerical measurements. However, it does allow for the utilization of an adequately sampled time average in place of the expectation value. The time-averaged estimate is defined as follows \cite{ErgodicitySpringerLink}:
\begin{equation}
    \hat A(t) = \frac{1}{t} \int_0^t \mathrm{d} t' A(t').
\end{equation}


\subsection{Markov Chain Monte Carlo}
To ensure the proper sampling of the Boltzmann distribution and to facilitate efficient exploration of the configuration space, we employ Markov Chain Monte Carlo (MCMC) methods. MCMC methods generate a sequence of configurations (states) where each configuration depends only on the previous one, forming a Markov chain. The fundamental objective is to establish a transition probability between states that satisfies detailed balance, thereby ensuring the system's convergence to the desired equilibrium distribution. In each Monte Carlo cycle, the objective is to flip each spin in the lattice precisely once. The acceptance of a proposed spin flip is determined using the Metropolis-Hastings algorithm, which involves calculating the energy change, denoted by $\Delta E$, associated with the flip. In the event that $\Delta E \leq 0$ the flip is accepted without condition, as it results in a lower energy state. In the event that $\Delta E > 0$, the flip is accepted with a probability that is determined by the Boltzmann factors
\begin{equation} \label{eq:acceptance_prob}
    P_{\text{accept}} = e^{-\Delta E/(k_B T)}.
\end{equation}
This acceptance criterion ensures that the system can explore higher energy states, allowing it to escape local minima and sample the configuration space more effectively. The acceptance of the flip is contingent upon the random number $r$ being drawn uniformly from the interval $[0, 1)$. More specifically, the flip is accepted if $r < P_{accept}$.

\subsection{Parallelization}
Parallelization of MCMC methods is challenging due to the sequential nature of Markov chains. However, it is possible to exploit the independence of multiple simulations at different temperatures or initial conditions to achieve parallelism \cite{glatt-holtzParallelMCMCAlgorithms2024}. The implementation of multiple independent Markov chains operating concurrently, with each chain simulating the Ising model at an distinct temperature, has been demonstrated to result in a substantial acceleration of the overall sampling process. Each chain functions independently, facilitating straightforward parallelization across multiple processors. This approach is particularly beneficial when studying temperature-dependent properties of the Ising model, as it enables simultaneous exploration of the system's behavior across a range of temperatures.

\subsection{Implementation Details}

The simulations were executed using a custom C++ implementation of the two-dimensional Ising model with periodic boundary conditions. Each simulation instance constructs an $L\times L$ lattice with spins $s_{ij} \in \{\pm 1\}$. The lattice is initialized in one of two ways: either in a fully ordered configuration or with randomly assigned spins drawn from a uniform distribution. The generation of random numbers is achieved via a Mersenne Twister engine (\texttt{std::mt19937}), with each run being seeded individually. To ensure that periodic boundary conditions are met, the neighboring indices are computed using modulo arithmetic
\begin{equation}
    \begin{aligned}
        \text{NN}(i,j) = \{&((i+1) \mod L, j), \\&((i-1 + L) \mod L, j), \\&(i, (j+1) \mod L), \\&(i, (j-1 + L) \mod L)\}.
    \end{aligned}
\end{equation}

The spin dynamics follow the standard single-spin Metropolis algorithm. A Monte Carlo step involves the selection of a lattice site uniformly at random, the computation of the local energy change $\Delta E$ associated with the flipping of that spin, and the acceptance of the flip with probability $P_\text{accept}$. In order to minimize computational expenditure, all conceivable Boltzmann factors (\cref{eq:acceptance_prob}) associated with the four-neighbor interaction are initially calculated and subsequently stored in a lookup table.

During the simulation, the instantaneous energy and absolute magnetization are evaluated after each attempted update. The values and their squares are then accumulated to compute ensemble averages and variances. Subsequently, the code derives the heat capacity and magnetic susceptibility from these moments, as outlined in \cref{eq:heat_capacity,eq:susceptibility}. Optional time-resolved data-energy, magnetization, and other intermediate quantities, are stored at each Monte Carlo cycle when evolution recording is enabled.

Furthermore, to support large parameter sweeps, the code incorporates a scheduling layer that manages multiple independent simulations. The configuration of each run is delineated by a specific lattice size, a designated temperature, the number of Monte Carlo cycles employed, the initialization mode, and the RNG seed. A master seed generator ensures statistical independence across runs. The scheduler distributes the runs across a number of CPU threads specified by the user, executing them concurrently. Upon completion of a given run, an output file is generated, containing the averaged observables. When the relevant option is enabled, an additional file is also generated, which contains the full time evolution of the system.


\subsubsection{Tools and Usage of AI}
For concurrent execution, the implementation utilizes the C++ library \texttt{OpenMP}. The executables were compiled using \texttt{g++ (GCC) 15.2.1 20251022 (Red Hat 15.2.1-3)} with optimization flag \texttt{-O3} enabled. All executions were carried out on an AMD Ryzen 7 PRO 4750U Processor, equipped with 16 GB of RAM, operating on the Fedora Linux 42 distribution. A large impact of IO operations on performance is not anticipated, given the utilization of an NVMe SSD for all read and write operations.

For the purposes of analysis, the \texttt{Python 3.13} programming language was utilized in conjunction with the libraries \texttt{NumPy} \cite{harrisArrayProgrammingNumPy2020}, \texttt{Pandas} \cite{teamPandasdevPandasPandas2025}, \texttt{SciPy} \cite{2020SciPy-NMeth}, and \texttt{Matplotlib} \cite{hunterMatplotlib2DGraphics2007}. In the course of composing this report, the large language model \texttt{ChatGPT} was employed to proofread various sections and to propose enhancements in terms of phrasing and grammar. All figures were created manually, without the use of AI assistance. During the implementation of the code, the large language model provided by \texttt{GitHub Copilot} was used to improve coding speed, but not for generating large code blocks. A comprehensive review and testing of the code was conducted to ensure its correctness. 

