Abstract: This project revolves around simulating and exploring the Ising model, with a focus on studying its behavior and properties. We utilize the Metropolis-Hastings algorithm to obtain the stationary distribution by evaluating the Hamiltonian. By conducting MCMC simulations, We calculate the heat capacity, magnetization and susceptibility of the model to analyze the phase transition phenomena within the 2D Ising model. Additionally, we investigate the impact of various parameters like temperature and interaction coefficients on the patterns observed in a two-dimensional setting. Moving forward, our plans involve applying the Ising model to a fully connected graph, incorporating mean-field theory, and using this framework to simulate complex systems such as spin-glass and melting crystals.
Setting up the Ising Model in Wolfram Language
Introduction
Introduction
This project revolves around the exploration and application of the Ising model, a mathematical model of ferromagnetism in statistical mechanics. The primary aim of the project is to analyze the behavior of spins in two-dimensional systems under various parameters using the Ising model.It exhibits phenomena such as symmetry breaking in the low-temperature phase and possesses a distinctive “critical point” at a specific temperature. Calculating these quantities is typically challenging due to the exponentially increasing number of states in a system, making explicit computation unfeasible. The Ising model provides clever approaches to compute the partition function, which researchers greatly appreciate. Despite its simplicity, the Ising model finds applications in diverse systems, showcasing the concept of universality in critical phenomena. This universality, where the same theory applies to various phase transitions, fascinates physicists as it suggests an underlying order in our complex world.
Setting up the initial configuration
Setting up the initial configuration
The generalized Ising Model is a mathematical model that doesn’t correspond to an actual physical system. It’s a huge (square) lattice of sites, where each site can be in one of two states. We label each site with an index , and we call the two states -1 and +1. To say that the site is in the state -1, we write = -1.
i
i
th
σ
i
Create the initial Ising model configuration:
In[]:=
Wholespins[dimension_,scale_]:=RandomChoice[{-1,1},scale^dimension]Initialspins[dimension_,scale_]:=ArrayReshape[Wholespins[dimension,scale],ConstantArray[scale,dimension]](*dimension_meansthedimensionoftheconfiguration;scale_meansthesizeoftheconfiguration*)
Calculating the Hamiltonian of the model without the outside filed
Calculating the Hamiltonian of the model without the outside filed
The energy of a configuration σ is given by the Hamiltonian function the notation indicates that sites are nearest neighbors .
(<ij>)
iandj
H
σ
∑
<i,j>
J
ij
σ
i
σ
j
The sign of tells you whether neighbors prefer to be ferromagnetic or anti-ferromagnetic. Physically, whether a material is one or the other (or neither) depends on the exact quantum mechanical details of how the spins interact.Since the spin-spin interaction arises between two spins, we have to sum over pairs of sites to find its total interaction to the energy.
The configuration probability is given by the Boltzmann distribution with inverse temperature β ≥ 0,
J
The configuration probability is given by the Boltzmann distribution with inverse temperature β ≥ 0,
P
β
-βH(σ)
e
Z
β
Z
β
∑
σ
-βH(σ)
e
where β = , and is the normalization constant.
-1
(T)
K
B
Z
β
Using convolution to calculate the Hamiltonian of every spin:
In[]:=
kernel[initial_]:=Which[ArrayDepth[initial]==2,{{0,1,0},{1,0,1},{0,1,0}},ArrayDepth[initial]==1,{1,0,1}]nbresults[initial_]:=ListConvolve[kernel[initial],initial,2]interenergy[initial_,j_]:=-j*initial*nbresults[initial](*initial_istheinitialconfiguration,createdbyusingInitialspinsfunction*)
Markov chain Monte Carlo simulation by Metropolis-Hasting Algorithm
Markov chain Monte Carlo simulation by Metropolis-Hasting Algorithm
The Metropolis-Hastings algorithm is widely used as a Monte Carlo method for estimating Ising model properties. To perform the algorithm, selection probabilities denoted as g(μ, ν) are initially determined. These probabilities indicate the likelihood of selecting state ν from all possible states, given that the current state is μ. Acceptance probabilities, denoted as A(μ, ν), are then used to ensure detailed balance is satisfied. If the proposed state ν is accepted, the algorithm moves to that state and continues by selecting another state and deciding its acceptance.If ν is not accepted, the algorithm remains in state μ. This process iterates until a specified stopping criterion is met, often when the lattice becomes ferromagnetic, indicating that all sites point in the same direction.
During implementation, it is crucial to select g(μ, ν) in a manner that guarantees ergodicity. In thermal equilibrium, a system’s energy fluctuates within a narrow range. This concept motivates the adoption of single-spin-flip dynamics, which dictates that only one spin site on the lattice is altered in each transition. Furthermore, by employing single-spin-flip dynamics, it becomes possible to transition from any state to any other state by flipping each differing site between the two states, one at a time.
In the case of the Ising model and employing single-spin-flip dynamics, we can deduce the following relationship.
Considering that the lattice consists of L sites in total, and single-spin-flip is the sole mechanism for transitioning between states, we observe that there exist L new states ν from the current state μ. In this scenario, the algorithm assumes that the selection probabilities for these L states are equal, namely g(μ, ν) = 1/L. To satisfy the principle of detailed balance, the following equation must hold:
During implementation, it is crucial to select g(μ, ν) in a manner that guarantees ergodicity. In thermal equilibrium, a system’s energy fluctuates within a narrow range. This concept motivates the adoption of single-spin-flip dynamics, which dictates that only one spin site on the lattice is altered in each transition. Furthermore, by employing single-spin-flip dynamics, it becomes possible to transition from any state to any other state by flipping each differing site between the two states, one at a time.
In the case of the Ising model and employing single-spin-flip dynamics, we can deduce the following relationship.
Considering that the lattice consists of L sites in total, and single-spin-flip is the sole mechanism for transitioning between states, we observe that there exist L new states ν from the current state μ. In this scenario, the algorithm assumes that the selection probabilities for these L states are equal, namely g(μ, ν) = 1/L. To satisfy the principle of detailed balance, the following equation must hold:
P(μ,ν)
P(ν,μ)
g(μ,ν)A(μ,ν)
g(ν,μ)A(ν,μ)
A(μ,ν)
A(ν,μ)
P
β(ν)
P
β(μ)
1
Z
-β()
H
ν
e
1
Z
-β()
H
μ
e
-β(-)
H
ν
H
μ
e
Metropolis algorithm sets the acceptance probability to be the larger value between A(μ, ν) and A(ν, μ), ensuring one of them is always equal to 1. Therefore, the acceptance algorithm can be described as follows:
A(μ,ν)=
|
The implementation form of Markov Monte Carlo algorithm is based on the Boltzmann distribution.
In[]:=
deltaenergy[initial_,j_]:=-2interenergy[initial,j]Index[initial_]:=RandomInteger[{1,Length[initial]},ArrayDepth[initial]]boltzmann[initial_,T_,index_,j_]:=Exp[-Extract[deltaenergy[initial,j],index]/T]mcmcsimu[initial_,T_,index_,j_]:=If[RandomReal[]<boltzmann[initial,T,index,j],ReplacePart[initial,index->-1*Extract[initial,index]],ReplacePart[initial,index->Extract[initial,index]]](*intial_meanstheinitialconfiguration,createdbyusingInitialspinsfunction;jmeanstheinteractioncoefficientsindex_meanstheindexofspinintheconfigurationT_meanstemperature*)
Detection of phase transition
1D situation (No phase transition)
1D situation (No phase transition)
The energy of the ground state is given by -JL when all spins are in the same direction. In any other configuration, the additional energy is equal to 2J multiplied by the number of sign changes encountered while scanning the configuration from left to right.
Let’s denote the number of sign changes in a configuration as k. The energy difference from the ground state can then be expressed as 2k. Since energy is additive with respect to the number of flips, the probability p of a spin flip at each position is independent. The ratio of the probability of finding a flip to the probability of not finding one can be described by the Boltzmann factor:
Let’s denote the number of sign changes in a configuration as k. The energy difference from the ground state can then be expressed as 2k. Since energy is additive with respect to the number of flips, the probability p of a spin flip at each position is independent. The ratio of the probability of finding a flip to the probability of not finding one can be described by the Boltzmann factor:
p
1-p
-2βJ
e
The partition function can be understood as the total number of configurations, where each configuration is weighted by its corresponding Boltzmann weight. In this context, the configurations can be characterized by the number of sign changes they possess. Importantly, the partition function can be factorized, meaning it can be expressed as a product of individual factors associated with each sign-change configuration.
Z==(1+p)=
∑
configs
∑
k
S
k
e
∏
k
L
(1+p)
The logarithm divided by L is the free energy density:
βf=log(1+p)=log1++1
-2βJ
e
-2βJ
e
In the context of phase transitions, a key characteristic is the presence of a non-analytic free energy. Therefore, in the case of the one-dimensional model, it does not exhibit a phase transition since its free energy remains analytic throughout
2D situation -- Phase transition
2D situation -- Phase transition
In systems with two or more dimensions, the Ising model exhibits a significant change called a phase transition, characterized by a transition from an ordered state to a disordered state. To study this phenomenon in the 2D Ising model, we analyze two key order parameters: specific heat and susceptibility. By examining these observables, we can observe and understand the occurrence of the phase transition in the 2D Ising model.
Ising model simulation by applying algorithm
Ising model simulation by applying algorithm
To begin, we generate a 10x10 two-dimensional matrix with random values, which serves as our initial state. Subsequently, we execute 4000 iterations and retain the outcomes of the final 1000 iterations, guaranteeing that the system maintains an equilibrium state throughout the calculation.
Create the original configuration:
In[]:=
initialconfiguration=Initialspins[2,10]
Out[]=
{{-1,-1,-1,1,-1,-1,-1,-1,-1,1},{1,1,-1,-1,1,-1,-1,-1,-1,-1},{-1,-1,-1,1,-1,-1,1,-1,1,-1},{1,-1,-1,1,-1,-1,1,1,1,1},{-1,-1,1,-1,1,1,1,-1,-1,-1},{1,-1,1,-1,-1,-1,-1,-1,1,1},{1,-1,-1,1,-1,-1,1,-1,1,1},{-1,-1,-1,-1,-1,1,1,-1,-1,1},{1,-1,-1,1,-1,-1,1,-1,-1,-1},{-1,1,1,1,1,-1,1,1,-1,1}}
Record every step of simulation within a range of temperatures:
In[]:=
Tempstate[initial_,steps_,Tmin_,Tmax_,interval_]:=Table[NestList[mcmcsimu[#,t,Index[initial],1]&,initial,steps],{t,Tmin,Tmax,interval}](*steps_meansourwholesimulationstepsTmin_,Tmax_andinterval_meanstherangeofourtemperature*)
In[]:=
stateconfiguration=Tempstate[initialconfiguration,2500,0.1,10.1,0.1][[All,-1000;;-1,All,All]]
Out[]=
We take the last 1000steps results, to ensure that our calculation is based on the equilibrium state
Calculating Average Energy
Calculating Average Energy
Calculate the energy of interaction of every spin:
In[]:=
TempEnergystate=Map[interenergy[#,1]&,stateconfiguration,{-3}]
Out[]=
Calculate the whole energy of the configuration of every simulation steps:
In[]:=
wholeenergy=Total[Transpose[Map[Total[#]&,TempEnergystate,{-3,-2}]]]
Out[]=
{-400000,-400000,-400000,-400000,-400000,-400000,-400000,-400000,-299280,-400000,-391808,-398000,-400000,-283304,-389528,-380912,-361480,-227752,-290352,-277408,-219472,-320976,-187832,-210200,-275952,-180656,-189184,-164688,-180088,-156680,-152352,-172680,-148328,-150352,-121688,-127896,-134624,-128200,-117512,-113144,-102976,-105048,-115872,-90968,-85200,-96648,-78288,-86448,-86472,-76040,-101128,-76928,-61872,-85088,-71928,-90392,-60016,-77312,-73752,-78224,-76304,-65592,-57592,-60704,-64864,-59200,-52648,-53088,-71280,-55048,-59440,-56600,-62992,-56008,-73456,-46032,-45112,-65480,-55288,-55896,-58392,-45328,-47240,-51704,-56840,-53688,-48120,-52944,-47784,-46672,-33872,-51856,-36120,-50144,-50080,-20240,-48224,-50008,-43560,-34248,-41016}
Calculating heat capacity
Calculating heat capacity
so, we can get the average energy:
according to the definition of heat capacity and the average energy here, we can get the heat capacity formula:
Calculate the heat capacity by variants of Energy and T square,
here i means every temperature at which we do the calculation:
here i means every temperature at which we do the calculation:
We can clearly see the peak appears within a range of temperatures, which means the phase transition of the thermal properties, specifically heat capacity as order parameter of our system.
Despite the fluctuations, we can see the “jump” of heat capacity which represents as a peak in the plot. Since heat capacity is the second order derivatives of free energy, so we can infer from the “jump” here that a second-order phase transition occurred during the heating process. Considering the background of our Ising model, where spins are composed of neighboring spins, we are interested in observing the overall magnetization of the system, which refers to the overall alignment of spins in the system -- Magnetization.
Calculating Magnetization
Calculating Magnetization
Calculate the magnetization of the whole configuration
Calculate the absolute value of the average Magnetization
From the graph, we can observe that before and after the abrupt change, the magnetization undergoes a transition from 0 to a specific direction if we consider it from right to left. This indicates that the system transitioned from disorder to order, from a state of random and undefined direction to a state with a clear and specific direction. This also represents the breaking of a certain symmetry. In the current context, it implies the breaking of symmetry in the energy landscape with respect to spins being either +1 or -1. As is well known, this provides a physical explanation for the second-order phase transition.
After observing the system’s order parameter, “Magnetization,” we can further investigate the changes in the system’s properties. We can choose to examine the magnetic susceptibility as our observable. The magnetic susceptibility is defined as the sensitivity of the magnetization vector to an external field. We can calculate it using the partition function.
After observing the system’s order parameter, “Magnetization,” we can further investigate the changes in the system’s properties. We can choose to examine the magnetic susceptibility as our observable. The magnetic susceptibility is defined as the sensitivity of the magnetization vector to an external field. We can calculate it using the partition function.
Calculating Susceptibility
Calculating Susceptibility
We can also get the formula of susceptibility using the partition function above:
with:
Calculate the heat Susceptibility by variants of Magnetization and T
Unlike the magnetization vector, the magnetic susceptibility reflects the overall properties of the system. In theory, the magnetic susceptibility is inversely proportional to temperature and diverges at the phase transition temperature, exhibiting symmetry about the transition temperature. Due to the limitations of our simulation, there may be some issues, and the effects might not be very prominent. However, we can interpret the “peak” in the graph as the diverging point, indicating a change in the overall properties of the system after the phase transition. Near the phase transition point, the system becomes highly sensitive to external fields.
2D simulation
Here, we present the variations in the cluster behavior of the Ising model under different parameter changes in the two-dimensional case.
Manipulation
Manipulation
Experiments
Experiments
Varying T
Varying T
When we vary the temperature, particularly when transitioning from high temperatures to low temperatures, we can observe a convergence of spins, leading to the formation of cluster phenomena in the case of ferromagnetic materials (J>0).
As shown in the graph, where we selected temperatures ranging from 1.1 to 10 with intervals of 2.5, we can observe that the regions with neighboring spins pointing in the same direction become smaller and more scattered. This indicates that the system tends towards disorder at higher temperatures, but it also becomes more symmetric.
As shown in the graph, where we selected temperatures ranging from 1.1 to 10 with intervals of 2.5, we can observe that the regions with neighboring spins pointing in the same direction become smaller and more scattered. This indicates that the system tends towards disorder at higher temperatures, but it also becomes more symmetric.
Varying d
Varying d
As we adjust the lattice size, we find that larger lattices require more time for spins to reach complete uniformity within a given simulation time step. This is reasonable because as we increase the lattice size, the Monte Carlo algorithm demands more iterations to traverse every spin in the system.
In the graph provided, we examined the system by varying the lattice size from 5 to 30 with intervals of 7. It is evident that as we increase the lattice size, the system becomes less susceptible to becoming ordered. It needs lower temperatures to induce disorder in the system.
In the graph provided, we examined the system by varying the lattice size from 5 to 30 with intervals of 7. It is evident that as we increase the lattice size, the system becomes less susceptible to becoming ordered. It needs lower temperatures to induce disorder in the system.
Varying j
Varying j
When J > 0, indicating a ferromagnetic material, the larger the value of J, the less time it takes for the spins to reach complete magnetization. On the other hand, when J < 0, representing an antiferromagnetic material, at very low temperatures, the spins in the system tend to align in the opposite direction to their surrounding spins.
In our study, we explored the range of parameter j from -2.5 to 3 with intervals of 1.5. It is observed that when j is negative, neighboring spins tend to point in opposite directions, and this phenomenon becomes more pronounced as the absolute value of j increases. On the other hand, when j is positive, the system becomes more prone to orderliness, and this tendency becomes stronger as the magnitude of j increases.
In our study, we explored the range of parameter j from -2.5 to 3 with intervals of 1.5. It is observed that when j is negative, neighboring spins tend to point in opposite directions, and this phenomenon becomes more pronounced as the absolute value of j increases. On the other hand, when j is positive, the system becomes more prone to orderliness, and this tendency becomes stronger as the magnitude of j increases.
Concluding remarks
Concluding remarks
After running the Monte Carlo simulation, we obtained the configuration of the system. First, we calculated the heat capacity of the system and found a “jump,” indicating the occurrence of a phase transition. Considering the characteristics of the system, we then calculated both the magnetization vector (“magnetization”) and the magnetic susceptibility to observe the effects of the second-order phase transition on the system from both a phenomenological and a property perspective, focusing on the breaking of symmetry.
Moreover, we can observe that when the interaction coefficients are positive, indicating a ferromagnetic material, they are positively correlated with the occurrence of cluster motion in the system. However, in our case, the effect is reversed for the lattice size. In the antiferromagnetic material, each spin eventually aligns itself in the opposite direction to its neighboring spins when the temperature decreases .
In our next step, we will thoroughly investigate the intricacies of the fully connected Ising model, potentially integrating mean-field theory to enhance our understanding. Subsequently, we will explore practical applications of the Ising model in real-world scenarios, such as crystal structures or spin glasses. This exploration will be accompanied by an analysis of the energy landscape, enabling us to delve into the fascinating realm of symmetry-related phenomena.
Moreover, we can observe that when the interaction coefficients are positive, indicating a ferromagnetic material, they are positively correlated with the occurrence of cluster motion in the system. However, in our case, the effect is reversed for the lattice size. In the antiferromagnetic material, each spin eventually aligns itself in the opposite direction to its neighboring spins when the temperature decreases .
In our next step, we will thoroughly investigate the intricacies of the fully connected Ising model, potentially integrating mean-field theory to enhance our understanding. Subsequently, we will explore practical applications of the Ising model in real-world scenarios, such as crystal structures or spin glasses. This exploration will be accompanied by an analysis of the energy landscape, enabling us to delve into the fascinating realm of symmetry-related phenomena.
Keywords
Keywords
◼
Ising model
◼
Collective pattern
◼
Phase transition
Acknowledgment
Acknowledgment
I would like to express my sincere gratitude to my mentors Jon, Hatem, Shivam, and Xerxes for their invaluable assistance in this project. I am also grateful to Johannes and Bowen for their guidance and support in conceptualizing and completing the coding aspects. I would also like to thank Wolfram Research and Stephen Wolfram himself for their help with this project. It has been a wonderful experience to meet everyone in a foreign land, and it is truly an unforgettable journey.
References
References
◼
Yaroslav Mikhaylovich, SECOND ORDER PHASE TRANSITIONS, https://ocw.tudelft.nl/wp-content/uploads/lecture6_reading.pdf
◼
Richard Fitzpatrick, The Ising model, 2006-03-29, https://farside.ph.utexas.edu/teaching/329/lectures/node110.html
◼
Wikipedia contributors. (2023, July 11). Ising model. In Wikipedia, The Free Encyclopedia. Retrieved 15:03, July 12, 2023, from https://en.wikipedia.org/w/index.php?title=Ising_model&oldid=1164815476
◼
Wikipedia contributors. (2023, June 21). Partition function (statistical mechanics). In Wikipedia, The Free Encyclopedia. Retrieved 15:04, July 12, 2023, from https://en.wikipedia.org/w/index.php?title=Partition_function_(statistical_mechanics)&oldid=1161281168
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Exploring spin dynamics: an Ising model study
by Hengyuan Xu
Wolfram Community, STAFF PICKS, July 12, 2023
https://community.wolfram.com/groups/-/m/t/2959065
by Hengyuan Xu
Wolfram Community, STAFF PICKS, July 12, 2023
https://community.wolfram.com/groups/-/m/t/2959065