\subsection{Validation against Analytical Results}
In \cref{app:analytical} of the appendix, a summary of the expected behavior of a correctly sampled lattice of $L=2$ is provided, in terms of its expectation values. Subsequent to the generation of said results, a comparison is formulated against a set of 10 simulations, which were executed for a duration of \num{10000} Monte Carlo cycles at a temperature of $k_B T = \num{1.0}$. To assess the uncertainty associated with each ensemble mean, the standard deviation of the ensemble values is employed. The uncertainty in the relative error is derived through the application of Gaussian error propagation to the ensemble variances.  The results are presented in \cref{tab:comparison} and demonstrate a high degree of congruence with the anticipated analytical expectation values for the energy and magnetization of the system. Both quantities manifest relative deviations in the order of \qtyrange{e-3}{e-4}{\percent}. The comparison indicates discrepancies of \qty{3.4959}{\percent} for the specific heat capacity of the system, accompanied by a substantial uncertainty range of $>\qty{30}{\percent}$ . In a similar vein, the relative error of the susceptibility demonstrates a large deviation of $\qty{4.2046+-35.9493}{\percent}$. The significant uncertainty hinders the rise of any tension between the theoretical prediction and the simulated quantity. Consequently, the observed alignment is deemed to be highly satisfactory in accordance with theoretical expectations, thereby validating the implementation of the model.

A comparatively extended simulation time frame was selected for the lattice size. However, a convergence analysis reveals that a minimum of \num{1000} steps are necessary to achieve results that are equivalent in quality. In order to execute this analysis, the deviation between time averages $\hat \epsilon, |\hat m|$ and theoretical predictions $\expect{\epsilon}, \expect{|m|}$ as a function of the Monte Carlo cycle for each of the simulations. The results are displayed in \cref{fig:convergence_2x2}. Following a brief period of equilibration, the presence of artifacts becomes evident within the range of \numrange{e-2}{e-3}, with an upper boundary delineating the maximum deviation. This observation signifies the stability of the system within its equilibrium state.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{../../src/project4/python/convergence_analysis.pdf}
    \caption{Evolution of simulated 2x2 lattice structures over the duration of \num{10000} Monte Carlo cycles. Comparison of ensemble averages $\hat \epsilon, |\hat m|$ and theoretical predictions $\expect{\epsilon}, \expect{|m|}$. We use the values from \cref{tab:analytical_expectation_values} as the prediction values.}
    \label{fig:convergence_2x2}
\end{figure}


\begin{table}
    \caption{Comparison of ensemble properties with results obtained from analytical calculations in \cref{tab:analytical_expectation_values}. Ensemble properties sampled from 10 simulations with $L=2$, \num{10000} Monte Carlo cycles and random initialization.}
    \label{tab:comparison}
    \begin{tabular}{ccc}
        \toprule
        Quantity & Ensemble Average & Relative Error\\
        \midrule
        $\expect{\epsilon}$ & \num{-1.9958+-0.0014} & \qty{0.0071+-0.0692}{\percent} \\
        $\expect{|m|}$ & \num{0.9987+-.0005} & \qty{0.0004+-0.0463}{\percent} \\
        $\expect{\frac{C_V}{N}}$ & \num{0.033+-.011} & \qty{3.4959+-34.3238}{\percent} \\
        $\expect{\frac{\chi}{N}}$ & \num{0.0038+-.0014} & \qty{4.2046+-35.9493}{\percent} \\
        \bottomrule
    \end{tabular}
\end{table}

\subsection{Equilibration of the System} \label{subsec:equil}
In order to provide further validation of the working principle of our implementation, it is possible to evaluate the equilibration/thermalization process. The system evolution is recorded by initializing the system from either a fully ordered state, where all the spins are parallel, or a fully disordered state, where the spin directions are chosen from an entirely uncorrelated uniform distribution. At a constant temperature, a specific degree of ordering is anticipated to arise. Below the critical temperature of the system, the energetic advantage of ordering is sufficient to cause nearly complete sorting of the system. At temperatures above the critical temperature, a disruption in the microscopic ordering of spins becomes possible, resulting in the loss of magnetization at the macroscopic level of the system. The assessment of disruptions to the system's optimal order necessitates the utilization of the normalized energy of the system. An optimal ordering of the system will result in minimal energy, while any disruption to this ordered structure will result in an increase in energy. A steady normalized energy indicates successful equilibration, i.e., the energy is in a steady state relative to the system's temperature. In order to illustrate the discrepancies occurring on very small time scales in comparison to the behavior observed on the long time scale, the following quantities are graphed: $\epsilon$ is the current normalized energy of the system at a specific time step and the averaged normalized energy $\hat\epsilon$:
\begin{equation}
    \hat \epsilon(t) = \frac{1}{t} \int_0^t \mathrm{d} t' \epsilon(t').
\end{equation}

For a system of $L=20$, a total of \num{2000} Monte Carlo cycles are simulated and the normalized energies $\epsilon, \hat \epsilon$ are recorded. This procedure is repeated for two temperatures: $k_BT = \num{1.0}$ and $k_B T = \num{2.4}$ and the two extremes of ordering, i.e., fully ordered and fully randomized initial configurations.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{../../src/project4/python/energy_convergence.pdf}
    \caption{Evolution of the normalized current energy $\epsilon$ (slightly colored line) and the normalized time-averaged energy $\hat{\epsilon}$ (fully colored lines) for systems of lattice size $L=20$. }
    \label{fig:energy_evolution}
\end{figure}

The results of this study are shown in \cref{fig:energy_evolution}. The observed behavior is consistent with the anticipated outcomes. Specifically, it has been determined that the system approaches a state of equilibrium energy for both temperatures after approximately \num{1000} Markov Chain Monte Carlo cycles. When the fully unordered state is employed as the initial condition, a higher-than-equilibrium energy state is achieved, accompanied by a rapid decay towards the thermalized state. The temperature of $k_B T = 1$ results in complete ordering of the system, irrespective of the initial configuration. For lower temperatures, fewer stochastic deviations from equilibrium are observed, indicating that the system remains predominantly ordered even on short time scales. For the higher temperature point that was sampled, there are noticeable deviations from the time-averaged energy $\hat \epsilon$. On the long-time-scale average, the system nevertheless approaches an equilibrium state with respect to the energy. In summary, the model's behavior aligns with the theoretical framework.

\subsection{Phase Transition Behavior}
In order to observe the transition between the two phases of the two-dimensional Ising model, a simulation was conducted in which randomly initialized systems were observed at different temperature points for \num{50000} MCMC cycles. A range of temperatures, from $k_B T = \numrange{1.0}{2.4}$  is scanned at an offset between simulations of $k_B \Delta T = \num{1e-3}$. In order to accelerate the simulation, no output is recorded during the simulation, and the systems parameters are extracted solely at the end of the simulation. This procedure is carried out for four lattice sizes, namely for $L \in \{40, 60, 80, 100\}$. The resulting system parameters are then plotted as a function of the set temperature and the lattice size. A plot containing the observables can be found in  \cref{fig:observables}. The average energy of the system increases linearly with the temperature, irrespective of the system's size. For temperatures above the critical temperature, the system exhibits minimal spread, with larger systems demonstrating increased energies.  In the context of absolute magnetization, a decline is observed above the critical temperature threshold. This phenomenon can be attributed to the disruption of the global magnetization, which results in the manifestation of averaging effects within the lattice. The magnitude of the drop is more pronounced with increasing lattice sizes, as the averaging effect is statistically more significant. For temperatures below the critical temperature, a decrease in average absolute magnetization with an increase in temperature is observed, as the added dynamic in the system reduces perfect ordering. The specific heat capacity and the susceptibility demonstrate a diverging tendency around the critical temperature. A comparative analysis of the temperature-dependent specific heat capacity reveals a uniform trend across all systems, with the exception of a shift in the point of divergence. This variation in the critical temperature is attributed to the influence of the lattice size on the system's behavior. As the lattice size increases, the degree of divergence concomitantly rises.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{../../src/project4/python/observable_vs_temperature.pdf}
    \caption{Calculated normalized energy $\expect{\epsilon}$, magnetization $\expect{|m|}$, specific heat capacity $\expect{\frac{C_V}{N}}$ and susceptibility $\expect{\frac{\chi}{N}}$ as a function of system temperature for different lattice sizes. Systems were evolved for \num{50e3} MCMC cycles.}
    \label{fig:observables}
\end{figure}

\subsubsection{Limit of Infinite Lattice Size}
The diverging behavior is used in a second step to approximate the point of the critical temperature for each lattice size. For both the specific heat capacity and the susceptibility data a fit to the function
\begin{equation} \label{eq:fitfunc}
    f(T) = A e^{-\frac{(T-\mu)^2}{2\sigma^2}} + C
\end{equation}
is performed. As $\mu$ denotes the maximum point of the Gaussian function, it is employed as an estimator for the critical temperature. The results of all fits are summarized in \cref{fig:fit}. The fit is performed using the function \texttt{optimize.curve\_fit} from the \texttt{scipy} package, with the minimization of squared residuals. Subsequently, the values obtained for the critical temperature in dependence of the lattice size are employed to approximate the critical temperature of an infinitely sized lattice. To elaborate, the critical temperature is plotted against the lattice size, and a secondary least squares minimization is performed on the relation.
\begin{equation}
    T_C(L) = T_C(\infty) + a L^{-1}.
\end{equation}
It is evident that two distinct values are obtained for the critical temperature of infinite lattices. Specifically, the first value is determined to be \num{2.284} and is derived from the data points extracted from the susceptibility. The second value, \num{2.271}, is derived from the measurements on the specific heat capacity. The fit to obtain those results is shown in \cref{fig:critical_temp_fit}. The data points obtained via the specific heat capacity demonstrate a high degree of agreement with the anticipated result, $k_B T_C = 2.269$. The discrepancy between the measurements obtained from the susceptibility suggests the necessity for a more sophisticated function to ascertain the critical temperatures from the simulations, i.e., improving \cref{eq:fitfunc}. In order to achieve more precise measurements, it is necessary to increase the number of samples taken in the region of the critical temperature. Due to the scarcity of computing resources, this step has not been implemented in the present study.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{../../src/project4/python/phase_transition_analysis.pdf}
    \caption{Fit to the specific heat capacity $\expect{\frac{C_V}{N}}$ and susceptibility $\expect{\frac{\chi}{N}}$ as a function of system temperature for different lattice sizes. Fit to \cref{eq:fitfunc}.}
    \label{fig:fit}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{../../src/project4/python/critical_temperature_extrapolation.pdf}
    \caption{Fit of the values obtained from \cref{fig:fit} to the relation of the system size dependent critical temperature.}
    \label{fig:critical_temp_fit}
\end{figure}

\subsection{Performance Analysis}
The performance of the simulation was evaluated through benchmarking in both a single-threaded and a 12-thread multithreaded configuration using the system time utility. For each of the executions, a total of 500 identical systems of lattice size $L = \num{100}$ were simulated for \num{1000} Monte Carlo cycles. Large numbers of identical systems are deliberately selected to minimize the impact of a subset of threads that complete execution before others. This process enhances the overall run time without exhausting all available resources. In order to facilitate the uninterrupted execution of background processes, a configuration of 12 threads was selected on a system with a total of 16 threads. The single-threaded execution required a wall time of \qty{217.634}{\s}, while the multithreaded run completed in \qty{29.580}{\s}. This corresponds to a speedup of
\begin{equation}
\mathcal S = \frac{217.634}{29.580} \approx 7.36 ,
\end{equation}
which is substantial, though below the ideal linear scaling expected for perfectly parallel work on 12 threads. The resulting parallel efficiency is therefore $\mathcal E = \mathcal S/12 \approx 0.61$.

The reported CPU times provide additional insight into the system's performance. The single-threaded version accumulated \qty{215.592}{\s} of user time, consistent with full utilization of one core. In contrast, the multithreaded run consumed \qty{339.713}{\s} of user time over a wall interval of \qty{29.580}{\s}. This corresponds to an average of roughly $11.5$ concurrently active hardware threads. This finding suggests that the parallel implementation utilizes the majority of available cores during the majority of the runtime.

One potential explanation for the observed limited parallel efficiency is overhead due to concurrent processing, including task delegation, I/O, and an augmented thermal load on the CPU, which can result in thermal throttling. Given that the CPU utilized for the evaluation is part of a notebook computer, the final point is particularly salient, given the prevalence of inadequate cooling capacity in compact designs. Subsequent studies are required to ascertain the impact of all factors on the reduced parallel efficiency.

\subsection{Statistical Properties of the Ising Model}
The evaluation of the probability distribution of the normalized energy is achieved through a comprehensive analysis of the distribution of epsilon for systems of $L=20$ at the two temperatures previously specified in \cref{subsec:equil}. For each set of initial conditions, 25 identical simulations were executed, each with a different seed, and the resulting trajectories were collected. The results are then displayed in a normalized histogram, using the counts of each unique energy value present in the dataset. The approximate probability density function is displayed in \cref{fig:pdf}.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{../../src/project4/python/energy_histograms.pdf}
    \caption{Normalized counts per discrete energy value of  simulations using $L=20$. Total of 50 simulation of \num{2000} cycles.}
    \label{fig:pdf}
\end{figure}