# Leveraging Neural Network on Detecting Phases and Phase Transitions in the Ising Model
Author: Hao-Yang Yen (顏浩洋)

## Abstract
The classical Ising model stands out as one of the simplest yet most significant models in statistical mechanics, primarily due to its phase transition phenomena. In dimensions $d ≥ 2$, the Ising model exhibits two distinct phases: the ordered phase and the disordered phase, with a phase transition that can be described by an order parameter dependent on temperature. However, obtaining an exact solution for the Ising model can be complex and challenging. This complexity underscores the importance of numerical methods in studying the model. Machine learning techniques have become essential and powerful tools in the study of many-body quantum mechanics and theoretical condensed matter physics. My goal is to leverage neural networks to tackle this challenge effectively. One effective numerical tool for this purpose is the Metropolis algorithm, a Monte Carlo method that facilitates the exploration of the system’s configurations. In this project, I aim to leverage the Metropolis algorithm to generate a comprehensive training dataset representing the Ising model across different phases. This dataset will then be used to train a neural network, enabling it to characterize the phase properties effectively. Through this approach, I hope to enhance our understanding of phase transitions within the Ising model.

## Introduction

### Classical Ising Model Overview

The *Ising model* is a fundamental model in statistical mechanics used to study phase transitions and critical behaviors, such as changes between magnetized and non-magnetized states as temperature varies. This model is built on a lattice, which is a grid of points (or "lattice points") where each point represents a particle or site with a magnetic property known as *spin*.

In the Ising model, the lattice (denoted as $\Lambda$) consists of points that are connected to neighboring points, forming a structure in one, two, or higher dimensions. The simplest setup of the Ising model is on a finite $d$-dimensional lattice arranged in a rectangular structure, e.g., $\Lambda = \{1, 2, \ldots, L_1\} \times \{1, 2, \ldots, L_2\} \times \cdots \times \{1, 2, \ldots, L_d\}$. Here, each point is directly connected to its nearest neighbors, a defining feature of the Ising model.

For example, see the two-dimensional lattice in the figure below.
$\Lambda = \{1, 2, \ldots, 5\} \times \{1, 2, \ldots, 6\}$ for lattice in the figure. This is simplest kind lattice is called *square lattice*, each site in the lattice connected with its four nearest four sites, the Ising model having this lattice structure is called square Ising model.

### Spin and Spin Configuration

In this model, each lattice point can exist in one of two states: *spin up* $(\uparrow)$ or *spin down* $(\downarrow)$, represented as $+1$ and $-1$, respectively. For each lattice point $k \in \Lambda$, the spin is represented by a variable $\sigma_k$, where $\sigma_k \in \{ +1, -1 \}$. The entire collection of these spins across the lattice forms a *spin configuration* (denoted $\sigma = (\sigma_k)_{k \in \Lambda}$), illustrated below.

### Statistical Mechanics Background

In a many-body system like this lattice, the number of possible spin configurations grows exponentially with the number of lattice points. For a finite $d$-dimensional lattice with $N = \prod_{i=1}^d L_i$ points, there are $2^N$ possible configurations. To analyze the system, we use *statistical mechanics*.

The *Hamiltonian* measures the system's total energy based on spin interactions. For the Ising model, we assume each point interacts only with its nearest neighbors, illustrated below.


### Hamiltonian and Energy Calculation

The *Hamiltonian* (without an external magnetic field) is:

$$
H(\sigma) = -J \sum_{\langle ij \rangle} \sigma_i \sigma_j \quad \text{where} \quad J > 0
$$

Here, $J$ is a coupling constant that determines the interaction strength between neighboring spins, and $\langle ij \rangle$ denotes the summation over all pairs of nearest neighbors.
This interaction energy $H(\sigma)$ reflects the tendency of neighboring spins to align (if $J > 0$), which is crucial for studying how ordered phases (such as magnetized regions) emerge.

#### Example: Two-spin System
Consider a simple two-spin system with spins $\sigma_1, \sigma_2$. The Hamiltonian is:

\begin{align}
H(\sigma) = -J \sigma_1 \sigma_2\;.
\end{align}

The possible configurations are $\{\sigma \in (+1, +1), (+1, -1), (-1, +1), (-1, -1)\}$. The corresponding energies are:

\begin{align}
H(+1, +1) & = -J, & H(+1, -1) & = +J\;, \\
H(-1, +1) & = +J, & H(-1, -1) & = -J\;.
\end{align}

### Probability of Spin Configurations

The Hamiltonian determines the probability of a spin configuration through the *Boltzmann distribution*:

\begin{align}
P_T(\sigma) = \frac{\exp(-H(\sigma) / (k_B T))}{\displaystyle\sum_{\text{all possible } \sigma} \exp(-H(\sigma) / (k_B T))}\;,
\end{align}

where $k_B$ is the Boltzmann constant and $T$ is the temperature. As temperature varies, this distribution governs the likelihood of different configurations and the overall behavior of the system.

This distribution can also be expressed in a simpler form:

\begin{align}
P_\beta(\sigma) = \frac{\exp(-\beta H(\sigma))}{Z}\;,
\end{align}

where $\beta = 1 / (k_B T)$ and $$Z = \sum_{\text{all possible } \sigma} \exp(-\beta H(\sigma))\;,$$ is the *partition function*, which encodes all thermodynamic information about the system. The partition function, $Z$, is essential for calculating macroscopic quantities like internal energy, specific heat, and magnetization, providing insight into phase transitions and the overall properties of the system.

#### Example: Two-spin System
Consider a simple two-spin system with spins $\sigma_1, \sigma_2$. The Hamiltonian is:

\begin{align}
H(\sigma) = -J \sigma_1 \sigma_2\;.
\end{align}

The possible configurations are $\sigma = (+1, +1), (+1, -1), (-1, +1), (-1, -1)$. The corresponding energies are:

\begin{align}
H(+1, +1) & = -J, & H(+1, -1) & = +J\;, \\
H(-1, +1) & = +J, & H(-1, -1) & = -J\;.
\end{align}

The partition function is:

\begin{align}
Z &= \sum_{\sigma} \exp(-\beta H(\sigma)) \\
&= \exp(\beta J) + 2 \exp(-\beta J) + \exp(\beta J)\nonumber\\
&= 2 e^{\beta J} + 2 e^{-\beta J}.
\end{align}

Using $\cosh(x) = (e^x+e^{-x}) / 2$, we find:

\begin{align}
Z = 4 \cosh(\beta J)\;.
\end{align}

The probability of a configuration, for instance $(+1, +1)$, is then:

\begin{align}
P_T(+1, +1) = \frac{e^{\beta J}}{4 \cosh(\beta J)}\;.
\end{align}

### Phase Transitions and the Renormalization Group

The arrangement of spins in the Ising model is incredibly complex, making it difficult to analyze the system directly. Therefore, statistical mechanics provides several approaches—both theoretical and numerical—to simplify and study this system. One of the most powerful theoretical tools for this purpose is the renormalization group (RG) method. The main idea behind the renormalization group is to simplify a system by focusing only on the parameters that significantly impact its macroscopic (large-scale) properties while ignoring others that have little effect. Imagine a cup of water, for instance: many interactions occur between the water molecules, such as orbital interactions, dipole interactions, and even quantum effects. However, not all these details matter when considering water’s large-scale properties, like temperature or phase (solid, liquid, gas). The renormalization group allows us to ignore irrelevant details and build a simpler model with only a few key variables—called order parameters—that determine the system’s phase and how it changes under different conditions. An order parameter is a simplified variable that captures a system’s essential, observable behaviors. As an example, temperature often serves as an order parameter because temperature changes can cause a clear shift in the system’s phase (for instance, water freezing into ice or evaporating into steam). For the Ising model, the individual spins are less relevant on their own, while the temperature (or its equivalent variable $\beta$) is the main order parameter that dictates the system’s large-scale behavior. The phase under lower temperature is called ordered phase and the phase under higher temperature is called diordered phase. In the next section, we’ll explore this numerically. Numerically, renormalization group methods can be compared to machine learning techniques. Machine learning also extracts relevant features (or parameters) from data to make predictions. This similarity motivates using machine learning to identify and characterize phases in the Ising model, with machine learning potentially serving as a powerful tool to understand complex physical systems.

#### Example: RG for a 1D Ising Model
For the 1D Ising model, the RG transformation reduces the lattice size while modifying the effective coupling constant. Starting with a coupling ($K$), the fixed point analysis reveals that no phase transition occurs in 1D ($K \to 0$ as $T \to \infty$). However, in 2D, the fixed point structure is more complex, leading to a critical temperature $T_c$.

### Triangular Ising Model
In its two-dimensional form, the Ising model exhibits a phase transition, making it a cornerstone of theoretical and computational physics. The triangular Ising model is a variant of this system where spins are arranged on a triangular lattice. This variation introduces additional complexities and insights due to the geometric frustration and higher coordination number of the triangular lattice.

#### Triangular Lattice and Spin Configuration

Consider a triangular lattice with $N$ sites, where each site $i$ contains a spin variable $\sigma_i$ that can take values $+1$ or $-1$. The Hamiltonian of the system, which describes the energy of a particular spin configuration, is given by:

\begin{align}
H = -J \sum_{\langle i, j \rangle} \sigma_i \sigma_j\;,
\end{align}

This formulation is almost identical to that of the square Ising model, except that here $\langle i, j \rangle$ refers to the summation over the 6 nearest neighbors for each spin, as illustrated in the following figure:


The triangular lattice is characterized by a coordination number of 6, meaning each spin interacts with six neighboring spins. This increased connectivity compared to the square lattice leads to more complex spin interactions and behavior, making the triangular Ising model particularly interesting for studying phenomena like frustration and criticality.

### Critical Temperature and Phase Transition
In the absence of an external magnetic field ($h = 0$), the system undergoes a phase transition at a critical temperature $T_c$. Below $T_c$, the system is in the ordered ferromagnetic phase, where spins align, leading to spontaneous magnetization. Above $T_c$, thermal fluctuations dominate, and the system enters the disordered paramagnetic phase, where spins are oriented randomly and magnetization vanishes.

This transition is a key feature of the model, marking a shift from ordered to disordered states and serving as a cornerstone for studying phase transitions in statistical mechanics.

Thus, the triangular Ising model provides a rich framework for studying the interplay between lattice geometry, spin interactions, and phase transitions. Its unique properties, such as geometric frustration and higher connectivity, make it an intriguing system with applications in various fields such as physics, materials science, and computational problems like optimization. Analytical and numerical techniques continue to shed light on the fascinating behavior of this model, particularly in the presence of competing interactions and external fields.

### Application in Neural Network Studies

The triangular Ising model, due to its complex interactions and frustration effects, provides an excellent test case for evaluating neural networks' effectiveness in statistical mechanics. Neural networks, particularly those designed to predict phase transitions, can be employed to analyze the phase diagram of this system. By training the network on spin configurations at various temperatures, we can assess its ability to predict the system's behavior, especially in the critical region near $T_c$, where phase transitions occur.

Therefore, the triangular Ising model can serve as a valuable tool for testing whether neural network-based approaches are capable of accurately capturing phase transitions in complex systems like this one. As demonstrated in earlier studies, neural networks have shown impressive performance in classifying phases and predicting the transition temperature, suggesting their potential as powerful tools in statistical mechanics. Since we do not know the exact solution of the triangular Ising model now, we will try to find the critical temperature by the neural network method.

## Monte Carlo Simulation
Analyzing the Ising model theoretically is challenging and varies significantly across different lattice structures. However, we can often use the Monte Carlo method to approximate solutions numerically, even though this approach can be computationally intensive.

### Algorithm

Imagine we have a grid of particles (or "spins") arranged in a two-dimensional grid, where each particle can point either "up" ($+1$) or "down" ($-1$). This setup represents the spins on a lattice, with each point in the grid which is connected to its immediate neighbors. Each spin interacts with its neighbors, and these interactions determine the system's energy.

The simulation uses a popular algorithm called the Metropolis algorithm, which helps determine the likelihood of each spin changing direction based on the system's energy and temperature.

The psuedo code for the algorithms is as follows:

1. **Initialize** an $L \times L$ lattice with random spins $\sigma_i = \pm 1$ at each lattice point.
2. Set the temperature $T$ and calculate $\beta = 1 / (k_B T)$, where $k_B$ is the Boltzmann constant.
3. For each Monte Carlo step:
1. For each lattice site $(x, y)$:
- Select a random site $(x, y)$ on the lattice.
- Calculate the change in energy $\Delta E$ if the spin $\sigma_{x, y}$ were flipped: $$\Delta E = 2 J \sigma_{x, y} \sum_{\langle i,j \rangle} \sigma_j$$ where the sum is over the nearest neighbors of $(x, y)$.
- Generate a random number $r$ uniformly in the range $[0, 1]$.
- If $r < \exp(-\beta \Delta E)$:
- Flip the spin at $(x, y)$, that is, change the sign of $\sigma_{x,y}$: $$\sigma_{x, y} \to -\sigma_{x, y}\;.$$
2. Compute observables (e.g., average energy, magnetization) if needed.

Here, we iterate over each lattice site during each Monte Carlo step, computing the change in energy that would result from flipping the spin. The spin flip is accepted with a probability based on the temperature $T$, which allows the system to explore different configurations and gradually reach equilibrium.

* $L$: The lattice size (number of spins along each dimension).
* $J$: The interaction strength between neighboring spins.
* $\beta = 1 / (k_B T)$: The inverse temperature, where $k_B$ is the Boltzmann constant.
* $\sigma_{x, y}$: The spin at lattice site $(x, y)$.
* $\Delta E$: The energy difference if the spin $\sigma_{x, y}$ is flipped.

The algorithms simulates how the spins in this system evolve and eventually settle into an equilibrium state at a given temperature $T$.

The following animation demonstrates a dynamic process of the Ising model on a $30\times 30$ lattice becoming equilibrium at $T=1.5$ (left), $T=2$ (middle), $T=2.5$ (right) and the two colors in the animation represent spin up and spin down:

  

The MatLab code used to do the simulations is as follows:

using Random, Plots

function initialize_lattice(L, IC)
    if IC == 1
        spin = ones(Int, L, L)
    elseif IC == 2
        spin = -ones(Int, L, L)
    else
        spin = sign.(0.5 .- rand(L, L))
    end
    return spin
end


function monte_carlo_step!(spin, L, T)
    row, col = rand(1:L), rand(1:L)
    dE = 2 * spin[row, col] * (
        spin[mod1(row-1, L), col] + spin[mod1(row+1, L), col] + 
        spin[row, mod1(col-1, L)] + spin[row, mod1(col+1, L)]
    )
    
    if dE >= 0 
        if rand() < exp(-dE / T)
            spin[row, col] *= -1
            return 0  
        end
    else
        spin[row, col] *= -1
            return 0 
    end

end

# Parameters
L = 200 # Lattice size
IC = 3 # 1 for all spin up, 2 for all spin down, 3 for random
T = 1 # Temperature
steps = 10^5 # Number of simulation steps

# Initialize lattice
spin = initialize_lattice(L, IC)

# Run simulation with animation
anim = Animation()

for t in 1:steps
    monte_carlo_step!(spin, L, T)
    if t % 1000 == 1 
    heatmap(spin, color = :thermal, title = "lattice size=$(L)×$(L), temperature=$(T), time=$(t)")
    frame(anim)
    end
end

gif(anim, "ising.gif", fps = 50)

Notice that for different lattice structures, square lattice, and triangular lattice for example, we only need to change the definition of $\braket{i,j}$ in the algorithm. Therefore, this algorithm can be used to simulate various kinds of classical Ising models with different lattice structures by modifyting the function Neighbors in the MatLab code. The following animation demonstrates a dynamic process of the Ising model on a $30\times 30$ triangular lattice becoming equilibrium at $T=1$ (left), $T=2$ (middle), $T=3$ (right) and the two colors in the animation represent spin up and spin downas the previous animation:

this is an example of the lattice structure can afftect the dynamics of the system.

### Simulations
Here, we focus on four main physical properties in a two-dimensional Ising model simulation, which gives insights into how materials become magnetized or transition between different energy states with temperature changes. I simulate the two-dimensional square and triangular Ising model as an example in this section. Due to the computational complexity, we compute the average of the target physical quantities obtained from ten Monte Carlo simulations, which is a affordable number of samples. We use the following Julia code to do the simulations:

using Random, Plots

function initialize_lattice(L, IC)
    if IC == 1
        spin = ones(Int, L, L)
    elseif IC == 2
        spin = -ones(Int, L, L)
    else
        spin = sign.(0.5 .- rand(L, L))
    end
    
    En = 0
    for x in 1:L, y in 1:L
        En -= spin[x, y] * (spin[mod1(x+1, L), y] + spin[x, mod1(y+1, L)])
    end
    
    Mag = sum(spin)
    return spin, En, Mag
end

function monte_carlo_step!(spin, L, T)
    row, col = rand(1:L), rand(1:L)
    dE = 2 * spin[row, col] * (
        spin[mod1(row-1, L), col] + spin[mod1(row+1, L), col] + 
        spin[row, mod1(col-1, L)] + spin[row, mod1(col+1, L)]
    )
    
    if dE < 0 || rand() < exp(-dE / T)
        spin[row, col] *= -1
        return dE, 2 * spin[row, col]  # ΔE, ΔM
    end
    return 0, 0
end

# Parameters
L = 30
IC = 3
steps = 10^7
thermalization_steps = steps / 10

# Initialize lattice
spin, En, Mag = initialize_lattice(L, IC)

Temp = 1.5:0.1:3
Energy = Float64[]
Magnetization = Float64[]
MagSus = Float64[]
SpesHeat = Float64[]

# Thermalization
for _ in 1:steps
    monte_carlo_step!(spin, L, Temp[1])
end

for T in Temp
    # Measurement
    EnSum = 0.0
    MagSum = 0.0
    En2Sum = 0.0
    Mag2Sum = 0.0

    for _ in 1:steps
        dE, dM = monte_carlo_step!(spin, L, T)
        En += dE
        Mag += dM
        EnSum += En
        MagSum += abs(Mag)  # Use absolute value for symmetry breaking
        En2Sum += En^2
        Mag2Sum += Mag^2
    end

    EnMean = EnSum / steps / L^2
    MagMean = MagSum / steps / L^2
    En2Mean = En2Sum / steps / L^4
    Mag2Mean = Mag2Sum / steps / L^4

    push!(Energy, EnMean)
    push!(Magnetization, MagMean)
    push!(MagSus, ((L^2)/T) * (Mag2Mean - MagMean^2))
    push!(SpesHeat, ((L/T)^2) * (En2Mean - EnMean^2))
end

plot(Temp, Energy, label="Mean Energy", seriestype=:line, marker=:circle)
xlabel!("Temperature")
ylabel!("Internal Energy")
title!("Energy-Temperature Plot for 2D Ising Model")

plot(Temp, Magnetization, label="Mean Magnetization", seriestype=:line, marker=:circle)
xlabel!("Temperature")
ylabel!("Magnetization")
title!("Magnetization-Temperature Plot for 2D Ising Model")

plot(Temp, SpesHeat, label="Specific Heat", seriestype=:line, marker=:circle)
xlabel!("Temperature")
ylabel!("Spesific Heat")
title!("Spesific Heat-Temperature Plot for 2D Ising Model")

plot(Temp, MagSus, label="Magnetic Susceptibility", seriestype=:line, marker=:circle)
xlabel!("Temperature")
ylabel!("Magnetic Susceptibility")
title!("Magnetic Susceptibility-Temperature Plot for 2D Ising Model")

#### Internal Energy

The *internal energy* represents the average energy in the system. It's calculated using a formula called the Hamiltonian, which adds up the total energy for a particular state of the system. Mathematically, we represent it as:

\begin{align}
E =\mathbb{E}(H)= \sum_\sigma P_\beta(\sigma) H(\sigma)\;,
\end{align}

where $P_\beta(\sigma)$ gives the probability of each state. To make the energy independent of the system's size, we scale it by dividing by the total number of points on the lattice, $L^2$:

\begin{align}
E \to \frac{E}{L^2}\;.
\end{align}

We simulated this energy density across different temperatures, with results shown in the following figure. There’s a clear change, known as a *phase transition*, between temperatures $T = 2$ and $T = 2.5$.


On the other hand, to estimate a temperature interval contains the critical temperature, we also calculate the internal energy in a broad temperature range. The internal energy for the triangular Ising model shown in the following figure does not have an explicit phase transition point in the plot of internal energy. We will see the explicit phase transition point later.


#### Mean Magnetization

The *mean magnetization* describes how strongly the system becomes magnetized as spins (tiny magnetic units) align. The formula for magnetization in the Ising model sums up the spin values for all points in the system:

\begin{align}
\mathcal{M}(\sigma) = \sum_{k \in \Lambda} \sigma_k\;.
\end{align}

The *mean magnetization* takes the expected (average) value of this magnetization, adjusted by the probability of each state. We also scale magnetization by dividing by the system size:

\begin{align}
M = \mathbb{E}(\mathcal{M}) = \sum_\sigma P_\beta(\sigma) \mathcal{M}(\sigma)
\quad \text{and} \quad M \to \frac{M}{L^2}\;.
\end{align}

The result is shown in the following figure. It again shows a phase transition between $T = 2$ and $T = 2.5$.


Notice that the Ising model has so-called $\mathbb{Z}_2$ symmetry, that is if we change the sign of all spin, this does not change the partition functiom. As a consequence, we can set any direction parallel to the spin direction as $+z$ axis.

Moreover, in the following figure, we can see the explicit phase transition point between $T=1.8$ and $T=2$ in the plot of mean magnetization for the square Ising model as in square one.


#### Specific Heat

The *specific heat* measures how the energy of the system changes as temperature changes, especially around phase transitions where this change is rapid. It's calculated by examining the fluctuations in internal energy:

\begin{align}
C = \frac{1}{k_B T^2} \left( \mathbb{E}(H^2) - \mathbb{E}(H)^2 \right) = \frac{\beta}{T} \, \mathbf{Var}(H)\;,
\end{align}

where $k_B$ is the Boltzmann constant. The result in the following figure shows a rapid increase in specific heat near $T = 2$ and $T = 2.5$, indicating a phase transition.


We can not see the explicit phase transition point in the plot of specific heat for the square Ising model as the plot in the following figure:

Thus, we move to another physical quantity.

#### Magnetic Susceptibility

*Magnetic susceptibility* ($\chi$) shows how much the system becomes magnetized when exposed to an external magnetic field. It measures the fluctuations in magnetization, giving insight into how responsive the system is to magnetic fields:

\begin{align}
\chi = \frac{1}{k_B T} \left( \mathbb{E}(\mathcal{M}^2 )- \mathbb{E}(\mathcal{M} )^2 \right) = \frac{\beta}{T} \, \mathbf{Var}(\mathcal{M})\;;
\end{align}

As shown in the following figure, the magnetic susceptibility sharply changes around the phase transition point between $T = 2$ and $T = 2.5$, suggesting the system will be highly responsive to external fields.


Again, we can see the explicit phase transition point in the plot of magnetic susceptibility for the square Ising model as the following figure shows:

Magnetic susceptibility thus provides valuable insight into the system's behavior across different temperature regimes, particularly near the critical point.

In summary, these simulations reveal how energy, magnetization, specific heat, and magnetic susceptibility in a two-dimensional Ising model behave as temperature changes, especially highlighting phase transitions.

## Challanges
I aim to extend the application of neural network techniques to characterize phase transitions in high-dimensional Ising models. The Ising model is a fundamental framework in statistical mechanics that helps us understand how systems behave at different temperatures and the nature of phase transitions. However, when employing the Monte Carlo algorithm to simulate these systems, I face a significant challenge: we need to allow time for the magnetization to reach equilibrium.

This equilibration process can be quite slow, especially when the number of lattice points is large, which results in substantial computational costs associated with Monte Carlo simulations. The following figure demonstrates the internal energy and magnetization evolution in the Monte Carlo simulation, we can see that the magnetization in the disordered phase becomes equilibrium slowly, especially for the larger size ones:

 

We can also see that the fluctuation of the magnetization in the disordered phase is larger. This is the reason why we can see the simulation results in the disordered phase are more messy. This means that we may need more samples and a number of Monte Carlo simulation steps to improve the accuracy.

Moreover, the system having larger size becomes equilibrium slower than that having smaller size. The following animation demonstrates a dynamic process of the Ising model on a $200\times 200$ lattice becoming equilibrium at $T=1.5$ (left), $T=2$ (middle), $T=2.5$ (right) and the two colors in the animation represent spin up and spin down::

  

Therefore, as the system size increases, the complexity and the time required for the system to explore its configuration space grow, leading to longer waiting times for stable results. For this reason, I limit the system size and the system dimension on $30\times 30$ lattice to make the computational costs affordable.

While the two-dimensional Ising model serves as a valuable and illustrative example for understanding these concepts, it also highlights the computational limitations we face.

By using the two-dimensional Ising model as a case study, we can gain insights into the challenges involved in simulating phase transitions and begin to explore how machine learning methods might mitigate these computational demands. Ultimately, our goal is to develop a more efficient framework that can effectively characterize phase transitions in higher-dimensional systems without the extensive computational burden typically required by traditional Monte Carlo methods.

## Utilizing Neural Networks on Prediction Phases
Now I utilize the neural networks to predict the phases of the Ising model under different temperatures. I first generated the flattened matrices of the spin configurations in different phases. Part of the data is demonstrated in the following Figur:


These images are $30\times 30$ matrices and would be used as the training for the neural network model later.

The structure of the neural network is as follows:


and can be constructed by the following R code:

library(keras)
library(mlbench)
library(dplyr)
library(magrittr)
library(neuralnet)

# Step 1: Load the CSV files and add the phase column
df_1 <- read.csv('DS_ordered.csv')
df_1$phase <- factor("ordered")

df_2 <- read.csv('DS_ordered_higher.csv')
df_2$phase <- factor("ordered")
names(df_2) <- colnames(df_1)

df_3 <- read.csv('DS_disordered.csv')
df_3$phase <- factor("disordered")
names(df_3) <- colnames(df_1)

df_4 <- read.csv('DS_disordered_higher.csv')
df_4$phase <- factor("disordered")
names(df_4) <- colnames(df_1)

df_5 <- read.csv('DS_critical_lower.csv')
df_5$phase <- factor("ordered")
names(df_5) <- colnames(df_1)

df_6 <- read.csv('DS_critical_higher.csv')
df_6$phase <- factor("disordered")
names(df_6) <- colnames(df_1)

df_7 <- read.csv('DS_critic _higher.csv')
df_7$phase <- factor("disordered")
names(df_7) <- colnames(df_1)

df_8 <- read.csv('DS_critical.csv')
df_8$phase <- factor("ordered")
names(df_8) <- colnames(df_1)

# Step 2: Combine the data frames into one
df <- rbind(df_1, df_2, df_3, df_4)
df_near<-rbind(df_5,df_6)
df_critic<-rbind(df_8,df_7)

# Step 3: Split the data set into training and test sets
index<- sample(1:nrow(df), size = nrow(df)*0.5)
train<- df[index,]
test<- df[-index,]

# Step 4: Training the model 
model = neuralnet(
  phase~.,
  data=train,
  hidden=c(5,2),
  linear.output = FALSE
)
# Step 5: Predict the test set and see the accuracy 
pred <- predict(model, test)
labels <- c("ordered", "disordered")
prediction_label <- data.frame(max.col(pred)) %>%     
  mutate(pred=labels[max.col.pred.]) %>%
  select(2) %>%
  unlist()

test$phase<-as.numeric(test$phase)
test$phase<-3-test$phase
result_1<-data.frame(test$phase)
result_2<-data.frame(max.col(pred))
check = as.numeric(test$phase == max.col(pred))
accuracy = (sum(check)/nrow(test))*100
print(accuracy)
table(test$phase, prediction_label)

# Step 6: Predict the data set near the critical point
df <- rbind(df_1, df_2, df_3, df_4)

index<- sample(1:nrow(df), size = nrow(df)*0.5)
train<- df[index,]
test<-df_near

model = neuralnet(
  phase~.,
  data=train,
  hidden=c(5,2),
  linear.output = FALSE
)

pred <- predict(model, test)
labels <- c("ordered", "disordered")
prediction_label <- data.frame(max.col(pred)) %>%     
  mutate(pred=labels[max.col.pred.]) %>%
  select(2) %>%
  unlist()

pred <- predict(model, test)
labels <- c("ordered", "disordered")
prediction_label <- data.frame(max.col(pred)) %>%     
  mutate(pred=labels[max.col.pred.]) %>%
  select(2) %>%
  unlist()

test$phase<-as.numeric(test$phase)
result_1<-data.frame(test$phase)
result_2<-data.frame(max.col(pred))
check = as.numeric(test$phase == max.col(pred))
accuracy = (sum(check)/nrow(test))*100
print(accuracy)
table(test$phase, prediction_label)

The number of input parameters is the same as the number of the lattice points in the Ising model and the two output parameters correspond to the ordered and disordered phase. Notice that the number of hidden layers and the number of nodes of the hidden layers can be modified.

### Square Ising Model
We first consider utilizing the neural network with five hidden layers, which have parameters of 500, 200, 100, 50, and 10 to predict the phases of the square Ising model. To evaluate the model's performance, I generated 400 samples using the Monte Carlo method, evenly distributed across four different temperatures $T=1.5,\;2,\;2.5,\;3$. Each temperature contributed 100 samples to the dataset. These samples were then randomly split evenly into a training set and a test set, ensuring equal representation from each temperature range.

For the classification task, the ordered phase samples were taken from the temperature range $1.5 \leq T \leq 2$, while the disordered phase samples were selected from $2.5 \leq T \leq 3$. The model demonstrated exceptional predictive accuracy, as shown in the confusion matrix below:

\begin{array}{c|c|c}
& \text{Ordered} & \text{Disordered} \\
\hline
\text{Ordered} & 113 & 0 \\
\hline
\text{Disordered} & 1 & 86 \\
\end{array}

The overall prediction accuracy was calculated as:

\begin{align}
\text{Accuracy} = \frac{113 + 86}{200} = 99.5\%\;,
\end{align}

indicating a highly reliable classification performance. Since the data set was randomly split, we can consider different training and test sets, no matter what the training and test sets are, the prediction results are all very good, the accuracy of ten prediction results with different training and test sets is demonstrated below:

\begin{array}{|c|c|}
\hline
\text{Index} & \text{Prediction Accuracy (%)} \\
\hline
1 & 99 \\
2 & 99.5 \\
3 & 99.5 \\
4 & 100 \\
5 & 99 \\
6 & 99.5 \\
7 & 98.5 \\
8 & 100 \\
9 & 99.5 \\
10 & 100 \\
\hline
\end{array}

We now focus on temperatures close to the critical temperature. For the two-dimensional classical Ising model on a square lattice without an external magnetic field, the critical temperature ($T_c$) is expressed as:

\begin{align}
T_c = \frac{2}{\ln(1+\sqrt{2})} \cdot \frac{J}{k_B}\;,
\end{align}

where $J$ is the interaction strength, and ($k_B$) is the Boltzmann constant. Assuming $J=1$ and $k_B=1$, this evaluates numerically to approximately:

\begin{align}
T_c \approx 2.26917...\;.
\end{align}

At this temperature, the system undergoes a phase transition from a ferromagnetic (ordered) phase ($T < T_c$) to a paramagnetic (disordered) phase ($T > T_c$). This landmark result was first derived by Lars Onsager in 1944 through his exact solution for the two-dimensional Ising model.

To study behavior near the critical temperature, we consider two temperatures close to $T_c$, specifically $T=2.26$ and $T=2.28$. Using the Monte Carlo method, we prepare 25 samples for each temperature. The results show that predictions for $T=2.26$ and $T=2.28$ yield slightly lower accuracies than other tested temperatures. The prediction accuracies for $T=2.26$ and $T=2.28$ are as follows:

\begin{array}{|c|c|}
\hline
\text{Index} & \text{Prediction Accuracy (%)} \\
\hline
1 & 92 \\
2 & 100 \\
3 & 96 \\
4 & 76 \\
5 & 72 \\
6 & 80 \\
7 & 100 \\
8 & 92 \\
9 & 92 \\
10 & 68 \\
\hline
\end{array}

Additionally, when considering $T=2.269$ and $T=2.2693$, slightly below and above the critical temperature, different training and test sets result in significantly lower prediction accuracies compared to other temperatures. The prediction accuracies for $T=2.269$ and $T=2.2693$ are as follows:

\begin{array}{|c|c|}
\hline
\text{Index} & \text{Prediction Accuracy (%)} \\
\hline
1 & 96 \\
2 & 64 \\
3 & 58 \\
4 & 96 \\
5 & 98 \\
6 & 76 \\
7 & 60 \\
8 & 84 \\
9 & 88 \\
10 & 92 \\
\hline
\end{array}

These observations highlight the challenges and variability in prediction accuracy near the critical temperature, suggesting sensitivity to temperature variations and the importance of careful dataset preparation.

### Triangular Ising Model

The triangular Ising model has the lattice structure is a two-dimensional triangular grid instead of the traditional square lattice. This modification introduces distinct features that affect the behavior of the system, especially concerning phase transitions and the nature of interactions. The triangular Ising model exhibits phase transitions similar to those observed in the square Ising model, but with different critical temperatures and lattice geometries that impact the system’s dynamics.

As in the square Ising model, the triangular Ising model undergoes a phase transition between two primary states: an ordered ferromagnetic phase at low temperatures and a disordered paramagnetic phase at high temperatures. The ordered phase is characterized by aligned spins (i.e., all spins are either $+1$ or $-1$), while in the disordered phase, spins are randomly oriented.

In following figures, we demonstrate three different spin configurations corresponding to different temperatures:


These images clearly highlight the transition from order to disorder as the temperature increases, a hallmark of phase transition phenomena.

For a triangular lattice, when there is no external magnetic field ($h=0$), the system undergoes a phase transition at a critical temperature $T_c$. The critical temperature can be approximated by the Monte Carlo simulation we have performed:

\begin{align}
T_c \approx 1.9\;
\end{align}

I will try to make the value more precise later.

At $T_c$, the system transitions from a ferromagnetic phase ($T < T_c$) to a paramagnetic phase ($T > T_c$), just like the square Ising model. Below $T_c$, the interactions between spins are strong enough to align them, leading to a stable ordered phase. Above $T_c$, thermal fluctuations dominate, and the spins become disordered, resulting in a paramagnetic phase.

To study the phase transition in the triangular Ising model, we apply a neural network-based approach, similar to the method used for the square Ising model. We prepare a set of spin configurations at two different temperatures: $T = 1$ (well below $T_c$) and $T = 3$ (well above $T_c$), and train the neural network to classify these phases. The neural network prediction results are very promising, yielding high accuracy in distinguishing between the ordered and disordered phases:

\begin{array}{|c|c|}
\hline
\text{Index} & \text{Prediction Accuracy (%)} \\
\hline
1 & 99 \\
2 & 99.5 \\
3 & 99.5 \\
4 & 100 \\
5 & 99 \\
6 & 99.5 \\
7 & 98.5 \\
8 & 100 \\
9 & 99.5 \\
10 & 100 \\
\hline
\end{array}

These results show that the neural network is highly effective at capturing the phase transition and can accurately predict the system's state at different temperatures.

Next, we explore temperatures close to the critical temperature ($T \approx T_c$), where the system exhibits mixed behavior due to the competition between ordered and disordered interactions. We prepared the samples at temperatures $T=1.75$ and $T=2.05$. At these temperatures, the neural network’s performance slightly decreases, likely due to the increased complexity of the spin configurations in this critical region:

\begin{array}{|c|c|}
\hline
\text{Index} & \text{Prediction Accuracy (%)} \\
\hline
1 & 92 \\
2 & 92 \\
3 & 96 \\
4 & 96 \\
5 & 72 \\
6 & 84 \\
7 & 96 \\
8 & 92 \\
9 & 92 \\
10 & 72 \\
\hline
\end{array}

Although the accuracy drops somewhat near the critical temperature, the neural network still performs well, demonstrating its robustness in handling complex and nuanced phase transitions.

Now we consider the temperatures closer to each other, consider the temperatures $T=1.9$ and $T=1.92$, and make an ansatz $T=1.9$ is in ordered phase and $T=1.92$ is in disordered phase, the prediction results are as follows:

\begin{array}{|c|c|}
\hline
\text{Index} & \text{Prediction Accuracy (%)} \\
\hline
1 & 96 \\
2 & 60 \\
3 & 58 \\
4 & 96 \\
5 & 88 \\
6 & 76 \\
7 & 72 \\
8 & 84 \\
9 & 88 \\
10 & 72 \\
\hline
\end{array}

The accuracies of the prediction results significantly decrease. However, it is still good (all of them are greater than $50\%$). Thus, we can guess that the critical temperature is in the interval $T=1.9$ to $T=1.92$. We can consider a smaller interval to get a more precise estimation.

Based on these results, we can conclude that the neural network method is an effective tool for studying the phase transitions in the triangular Ising model. Despite some challenges near the critical temperature, the method shows strong performance in distinguishing between the ordered and disordered phases. This approach is promising for analyzing more complex systems and can be extended to other lattice types or even more intricate models involving frustration and higher-dimensional lattices.

The number of input parameters is the same as the number of the lattice points in the Ising model and the two output parameters correspond to the ordered and disordered phase. Notice that the number of hidden layers and the number of nodes of the hidden layers can be modified. I employ the samples obtained from the Monte Carlo simulation to train the model. Then we can utilize the trained model to predict the phases of the Ising model having different spin configurations. The prediction of the samples far from the phase transition point is very accurate; however, the accuracy will decrease as the system is close to the phase transition point. I am now trying to find the best number of layers and nodes of the model. After finding it, we can utilize the model to approximate the phase transition temperature as well since the prediction ability of the model will decrease significantly when the phase transition occurs. The computational complexity of the model is large as well, it will take me some time to finish this. If time allows, I will try to extend this approach to the higher-dimensional Ising model or other kinds of Ising models such as the triangular lattice Ising model.

## Challenges for Further Exploration

In this report, we have demonstrated that the computational complexity involved in applying both the Monte Carlo method and neural networks to the classical Ising model is already substantial. This challenge becomes even more pronounced when extending to the *quantum Ising model*, a crucial framework in modern condensed matter physics. The quantum Ising model incorporates quantum effects, making it indispensable for understanding phenomena in systems where quantum mechanical interactions dominate, such as those encountered in low-temperature physics and strongly correlated materials.

The inclusion of quantum effects renders the classical approach, which relies on binary spin variables $(\sigma_k)_{k \in \Lambda}$ and a simple Hamiltonian function $H(\sigma)$, inadequate. Instead, the quantum Ising model demands a more advanced mathematical formalism involving **tensor products of complex vectors and matrices** to represent the quantum states and Hamiltonians, significantly increasing the complexity and computational requirements of the model.

### Quantum State Representation

In the quantum Ising model, the state of the system is represented as a *quantum state vector* in a high-dimensional Hilbert space. For a system of $N$ spins, the state is given by the tensor product:

\begin{align}
|\psi\rangle = \begin{pmatrix}\phi_1 \\ \varphi_1\end{pmatrix} \otimes \begin{pmatrix}\phi_2 \\ \varphi_2\end{pmatrix} \otimes \cdots \otimes \begin{pmatrix}\phi_N \\ \varphi_N\end{pmatrix}\;,
\end{align}

where $\phi_i, \varphi_i \in \mathbb{C}$ represent the complex coefficients associated with the quantum states of spin $i$.

### Quantum Hamiltonian Representation

Similarly, the Hamiltonian of the system, which governs the dynamics and interactions, is described by a tensor product of complex matrices:

\begin{align}
\hat{H} = \begin{pmatrix} H^1_{11} & H^1_{12} \\ H^1_{21} & H^1_{22} \end{pmatrix} \otimes \begin{pmatrix} H^2_{11} & H^2_{12} \\ H^2_{21} & H^2_{22} \end{pmatrix} \otimes \cdots \otimes \begin{pmatrix} H^N_{11} & H^N_{12} \\ H^N_{21} & H^N_{22} \end{pmatrix}\;,
\end{align}

where $H^k_{ij} \in \mathbb{C}$ are the matrix elements representing the quantum mechanical interactions for each spin.

### Increased Computational Complexity

The inclusion of quantum effects significantly increases the **computational complexity** of the problem. For a system with $N$ spins, the state space grows exponentially as $2^N$. Similarly, the Hamiltonian becomes a $2^N \times 2^N$ matrix, making calculations for large systems computationally demanding.

For example, a 10-spin system already requires a $2^{10} = 1024$-dimensional state space and a 20-spin system requires a $2^{20} = 1,048,576$-dimensional state space, which is beyond the capabilities of standard computational resources.

Preparing data for neural network training in the quantum Ising model is resource-intensive, requiring simulations of large quantum systems, generating sample configurations, and training models on high-dimensional data. These challenges demand advanced computing infrastructure, such as high-performance computing (HPC) resources, specialized hardware like GPUs or TPUs, and potentially quantum simulators or quantum computers. Overcoming these limitations opens doors to studying larger systems and gaining deeper insights into the interplay between quantum mechanics and statistical physics, paving the way for innovative applications in physics, materials science, and computational research.

Thus, although the quantum Ising model is a powerful framework for exploring the interplay between quantum mechanics and statistical physics, its study requires significant computational resources due to the complexity of quantum state representation and dynamics. As computational techniques and hardware continue to improve, we can expect to gain deeper insights into this fascinating model and its applications in physics, materials science, and beyond.

## Conclusion
In the study of classical Ising models, two distinct phases are typically observed, representing fundamental aspects of their behavior. By leveraging the Monte Carlo method, researchers can generate samples that provide a detailed representation of these phases. Neural networks have proven to be highly effective in analyzing these samples, offering a reliable means of identifying and predicting the phases with remarkable accuracy.

However, when the scope shifts to more complex Ising models, the computational demands grow exponentially, creating significant challenges. These demands stem from the intricate interactions and increased dimensionality inherent in such models, which push the limits of existing computational methods. This challenge represents a major hurdle in fully utilizing neural networks for these advanced scenarios.

Nonetheless, if these computational barriers can be overcome—whether through algorithmic advancements, optimized hardware, or innovative approaches—neural networks have the potential to become an extraordinarily powerful tool. They could provide unprecedented insights into the phase transitions and critical phenomena associated with statistical mechanics models, advancing our understanding of these systems to new heights.

## References
[1] Edwin Bedolla, Luis Carlos Padierna, and Ram´on Castaneda-Priego. “Machine learning for condensed matter physics”. In: Journal of Physics: Condensed Matter 33.5 (2020), p. 053001.
[2] Juan Carrasquilla and Giacomo Torlai. “Neural networks in quantum many-body physics: a hands-on tutorial”. In: arXiv preprint arXiv:2101.11099 (2021).
[3] Nigel Goldenfeld. Lectures on phase transitions and the renormalization group. CRC Press, 2018.
[4] Roland H¨aggkvist et al. “A Monte Carlo sampling scheme for the Ising model”. In: Journal of statistical physics 114 (2004), pp. 455–480.