Magnetization in the 1D Ising model and the Metropolis algorithm: a glimpse into statistical physics and Monte Carlo simulations

Magnetization in the 1D Ising Model and the Metropolis Algorithm A Glimpse into Statistical Physics and Monte Carlo Simulations Arpan Dey M1 Physics (IDIL) Université de Montpellier Acknowledgments I would like to express my sincere thanks to Dr. Jean-Charles Walter, Prof. Andrea Parmeggiani, and Linda Delimi for their guidance and support. This document is a pedagogical report on the one-dimensional Ising model, and was prepared during my bibliographic studies in advance of my internship at the Laboratoire Charles Coulomb (L2C), Université de Montpellier. Table of Contents Abstract 1 The 1D Ising Model: Analytical Treatment • The Model • The Partition Function • Correlation • Mean Field Approximation • Exact Solution of Magnetization 2 Monte Carlo Simulations • Statistical Physics • Ergodicity and the Postulate of Equal A Priori Probabilities • Principles of Monte Carlo Simulations • The Metropolis Algorithm • Space-Time Evolution of 1D Ising System • A Brief Look at the 2D Ising Model Discussion and Conclusion References Abstract In this report, we present a comprehensive exploration of the one-dimensional (1D) Ising model using a combination of exact analytical derivations and Monte Carlo simulations. The 1D Ising model is the simplest nontrivial lattice model that is exactly solvable, yet interesting enough (shows collaborative effects at low temperatures) to be relevant in biological physics, among other things. It is an ideal starting point for developing a proper understanding of numerical sampling techniques. Analytically, the 1D model exhibits no true thermodynamic phase transition at any finite temperature; we get a relatively high, nonzero magnetization at sufficiently low temperatures, which is a consequence of the finite size effect. The report bridges the conceptual foundation of statistical mechanics with practical computational approaches, demonstrating how simple stochastic rules can reproduce ensemble behavior consistent with thermodynamic predictions. Through comparison of exact transfermatrix results, mean-field approximations and Metropolis Monte Carlo data, we highlight how fluctuations and correlations – which are ignored in mean field theory – play a crucial role even in the simplest spin system. To extend beyond equilibrium averages, we also visualize the space-time evolution of the spin configurations in the 1D Ising system at low, intermediate and high temperatures. These discrete 2D plots reveal the transition from long-lived ordered domains to rapidly fluctuating disordered states as temperature increases, illustrating the competition between coupling 𝐽 and thermal agitation. Finally, we also briefly discuss the 2D Ising model. This report is brief and unoriginal, but unique in its pedagogical approach: instead of merely reproducing known results, it tries to use the 1D Ising model as a minimal laboratory for understanding ergodicity, equilibration and emergent order in stochastic systems. While higherdimensional models exhibit true phase transitions, the 1D case serves as an ideal testing ground for numerical methods, offering insight into how equilibrium arises from microscopic randomness. The 1D Ising Model: Analytical Treatment The Model We consider the 1D Ising model – essentially a long one-dimensional lattice with 𝑁 lattice sites, where each lattice can have a spin of either +1 or −1 (𝑠𝑖=±1). The system has 𝑍2 symmetry, which means if we flip all the spins (𝑠𝑖→−𝑠𝑖), the physics describing the system does not change. We assume a spin can interact only with adjacent spins (nearest neighbors), and the strength of interaction is characterized by the coupling constant 𝐽 (𝐽>0). We also assume periodic boundary conditions to avoid edge effects (we want all the spins to have two adjacent neighbors – we do not want the first and last spins to have one neighbor; so we are essentially talking about a “ring” instead of a “chain”). Apart from simplifying the mathematics, this ensures the system is uniform and can be treated to be a part of an infinite system; if we have a very large (infinite) system, the edge effects would not matter very much. Figure 1: A visualization of 1D Ising spins (Alice Heiman) Figure 2: A visualization of periodic boundary conditions on a 1D Ising chain of length N (https://doi.org/10.1088/1742-6596/2543/1/012006) Naturally, spins tend to align with their neighbors; if you have spin +1 and your neighbors on both sides have spin −1, you would want to flip your spin to −1 too. The Hamiltonian for the 1D Ising model is given by: 𝐻=−𝐽∑𝑠𝑖𝑠𝑖+1 𝑁 𝑖=1 Since the model is one-dimensional, we consider the spins of two adjacent lattice sites 𝑖 and 𝑖+1. The periodic boundary conditions can be enforced by setting 𝑠𝑁+1≡𝑠1. Each interaction between two adjacent spins contributes an energy: 𝐸𝑖,𝑖+1=−𝐽𝑠𝑖𝑠𝑖+1 If the spins are aligned (either both −1 or both +1), we have an energy contribution of −𝐽 from that bond, and if the spins are opposite, we get an energy contribution of +𝐽. Clearly, the energetically-favorable move would be to flip the spins in a way that they align with their neighbors; this ensures the system reaches the most stable state with the minimum energy. At very low temperatures, thermal fluctuations are insufficient to override this tendency and randomly flip spins – hence the system is ordered (all spins are aligned) and the magnetization 𝑚 (the overall average orientation) is nonzero and large. At high temperatures, thermal fluctuations break the order and introduce randomness in the system. As temperature 𝑇→∞, the magnetization drops to zero. It is now important to talk about the finite size effect. If we had an infinite chain with all spins aligned, there are infinitely many possible lattice sites where the spin may flip due to thermal fluctuations – even if the temperature is very slightly above zero. Thus, long range order is impossible for an infinite 1D Ising chain except when the temperature is exactly zero. So there is no true phase transition for any 𝑇>0. Now, if we had a finite-sized 1D Ising chain (size 𝑁) there are 𝑁 possible lattice sites, and the energy cost of flipping a spin that is already aligned with its neighbor is +𝐽−(−𝐽)=2𝐽. At sufficiently low (but still nonzero) temperatures, the probability of such a flip – which is proportional to exp(−2𝐽/𝑘𝐵𝑇)=exp(−2𝛽𝐽),𝛽=(1/𝑘𝐵𝑇) – is negligibly small. Hence, for a finite-sized 1D Ising chain, we have a nonzero magnetization at sufficiently low temperatures, and the magnetization smoothly drops to zero with increasing temperature, characteristic of a secondorder phase transition. The Partition Function The partition function of a system gives us information about the possible states the system can be in, and the corresponding probabilities. We will now intuitively derive the partition function for the 1D Ising model (in absence of any external magnetic field) using the idea of domain walls. A domain wall is present only between two opposite spins. Clearly, the energy of the chain if all 𝑁 spins were aligned (all spins either −1 or +1) is −𝐽𝑁, and if there are domain walls, each domain wall would contribute an extra energy of 2𝐽. Hence, if we assume there are 𝑘 number of domain walls in our 1D Ising chain, the energy of that configuration would be: 𝐸(𝑘)=−𝐽𝑁+2𝐽𝑘 The energy of any configuration with 𝑘 domain walls would be of the above form. The corresponding Boltzmann weight of such configurations would be: 𝑒−𝛽𝐸(𝑘)=𝑒−𝛽(−𝐽𝑁+2𝐽𝑘) The corresponding partition function would be: 𝑍= ∑ 2(𝑁 𝑘)𝑒−𝛽(−𝐽𝑁+2𝐽𝑘) 𝑁 𝑘=0 𝑘 𝑒𝑣𝑒𝑛 =2𝑒𝛽𝐽𝑁 ∑ (𝑁 𝑘)𝑥𝑘 𝑘 𝑒𝑣𝑒𝑛 for 𝑥=𝑒−2𝛽𝐽. We now recall the binomial theorem: (𝑥+𝑦)𝑁=∑(𝑁 𝑘)𝑥𝑘𝑦𝑛−𝑘 𝑁 𝑘=0 For 𝑥=𝑥,𝑦=1 and 𝑥=−𝑥,𝑦=−1, we respectively have: (1+𝑥)𝑁=∑(𝑁 𝑘)𝑥𝑘 𝑁 𝑘=0 and (1−𝑥)𝑁=∑(𝑁 𝑘)(−𝑥)𝑘= 𝑁 𝑘=0 ∑(𝑁 𝑘)(−1)𝑘𝑥𝑘 𝑁 𝑘=0 Now, we add (1+𝑥)𝑁 and (1−𝑥)𝑁: (1+𝑥)𝑁+(1−𝑥)𝑁=∑(𝑁 𝑘)𝑥𝑘 𝑁 𝑘=0 +∑(𝑁 𝑘)(−1)𝑘𝑥𝑘 𝑁 𝑘=0 =∑(𝑁 𝑘){1+(−1)𝑘}𝑥𝑘 𝑁 𝑘=0 =2 ∑ (𝑁 𝑘)𝑥𝑘 𝑁 𝑘=0 𝑘 𝑒𝑣𝑒𝑛 for even values of 𝑘. This implies: ∑ (𝑁 𝑘)𝑥𝑘 𝑁 𝑘=0 𝑘 𝑒𝑣𝑒𝑛 =(1+𝑥)𝑁+(1−𝑥)𝑁 2 Using this result in our partition function, and using 𝑥=𝑒−2𝛽𝐽, we get: 𝑍=2𝑒𝛽𝐽𝑁 ∑ (𝑁 𝑘)𝑥𝑘 𝑘 𝑒𝑣𝑒𝑛 =𝑒𝛽𝐽𝑁{(1+𝑒−2𝛽𝐽)𝑁+(1−𝑒−2𝛽𝐽)𝑁} Some trivial algebra gives: 𝑍=𝑒𝛽𝐽𝑁{𝑒−𝛽𝐽𝑁(𝑒𝛽𝐽+𝑒−𝛽𝐽)𝑁+𝑒−𝛽𝐽𝑁(𝑒𝛽𝐽−𝑒−𝛽𝐽)𝑁} =(𝑒𝛽𝐽+𝑒−𝛽𝐽)𝑁+(𝑒𝛽𝐽−𝑒−𝛽𝐽)𝑁=(2cosh𝛽𝐽)𝑁+(2sinh𝛽𝐽)𝑁 Now, we will show that as 𝑁→∞, the second term in the above expression becomes negligible. We know that: cosh𝑥=𝑒𝑥+𝑒−𝑥 2, sinh𝑥=𝑒𝑥−𝑒−𝑥 2,tanh𝑥=sinh𝑥 cosh𝑥 Since 𝑒𝑥+𝑒−𝑥>𝑒𝑥−𝑒−𝑥 for any real 𝑥>0, we see that cosh𝑥>sinh𝑥>0, which implies sinh𝑥 cosh𝑥=tanh𝑥<1. We may write: 2sinh𝑥 2cosh𝑥=tanh𝑥 Now, we rearrange the above expression and raise it to the 𝑁-th power: (2sinh𝑥)𝑁=(2cosh𝑥)𝑁(tanh𝑥)𝑁 Since tanh𝑥<1, (tanh𝑥)𝑁 decays exponentially as 𝑁→∞, which means under the limit 𝑁→∞: (2sinh𝑥)𝑁≪(2cosh𝑥)𝑁 Thus, for sufficiently large 𝑁, we get: 𝑍≈(2cosh𝛽𝐽)𝑁 The free energy per lattice site is given by: 𝑓=𝐹 𝑁=−𝑇 𝑁ln(𝑍)=−𝑇ln(2cosh𝛽𝐽) The above method works only in absence of any external magnetic field. If there is some constant nonzero external field ℎ, each spin on each lattice site is already biased. We can use the transfer matrix approach in this case. The new Hamiltonian with the external field ℎ is: 𝐻=−𝐽∑𝑠𝑖𝑠𝑖+1−ℎ∑𝑠𝑖 𝑁 𝑖=1 𝑁 𝑖=1 The second term encodes the effect of the external field ℎ, which wants to align all the spins and hence lower the energy (which explains the minus sign). Like before, we have periodic boundary conditions 𝑠𝑁+1≡𝑠1. Now, we define a symmetric matrix 𝑇𝑠,𝑠′ (𝑠,𝑠′≡𝑠𝑖,𝑠𝑖+1={±1} are the spins): 𝑇𝑠,𝑠′=exp[𝛽𝐽𝑠𝑠′+𝛽ℎ 2(𝑠+𝑠′)] We essentially split the effect of the field ℎ symmetrically between adjacent spins, and the products of the 𝑇 matrices, summed over all possible spin configurations {±1} gives the partition function: 𝑍= ∑ ∑ … ∑ 𝑇𝑠1,𝑠2𝑇𝑠2,𝑠3…𝑇𝑠𝑁,𝑠1 𝑠𝑁=±1𝑠2=±1𝑠1=±1 Using the result (𝑇𝑁)𝑎,𝑏=∑𝑇𝑎,𝑖1𝑇𝑖1,𝑖2…𝑇𝑖𝑁−1,𝑏𝑖1,𝑖2,…,𝑖𝑁−1 , and putting 𝑎=𝑏=𝑠1, we get: 𝑍= ∑(𝑇𝑁)𝑠1,𝑠1 𝑠1=±1 =𝑇𝑟(𝑇𝑁) where 𝑇𝑟(𝑇𝑁) refers to the trace of the matrix 𝑇𝑁, in other words the sum of its diagonal elements. There are four possible combinations of 𝑠 and 𝑠′, which gives the four elements of the matrix 𝑇𝑠,𝑠′: • When 𝑠=𝑠′=+1, 𝑇++=𝑒𝛽𝐽+𝛽ℎ • When 𝑠=+1,𝑠′=−1, 𝑇+−=𝑒−𝛽𝐽 • When 𝑠=−1,𝑠′=+1, 𝑇−+=𝑒−𝛽𝐽 • When 𝑠=−1,𝑠′=−1, 𝑇−−=𝑒𝛽𝐽−𝛽ℎ We now write down the transfer matrix 𝑇 and calculate its eigenvalues: 𝑇=(𝑒𝛽𝐽+𝛽ℎ 𝑒−𝛽𝐽 𝑒−𝛽𝐽 𝑒𝛽𝐽−𝛽ℎ) Setting the determinant of 𝑇, det(𝑇)=0, we obtain the characteristic equation of the eigenvalue problem: 𝜆2−𝜆(𝑒𝛽𝐽+𝛽ℎ+𝑒𝛽𝐽−𝛽ℎ)+(𝑒2𝛽𝐽−𝑒−2𝛽𝐽)=0 𝜆2−𝑇𝑟(𝑇)𝜆+det(𝑇)=0 The corresponding eigenvalues are: 𝜆±=𝑇𝑟(𝑇)±√(𝑇𝑟(𝑇))2− 4det (𝑇) 2 Now, the trace of the matrix 𝑇 is the sum of its diagonal elements: 𝑇𝑟(𝑇)=𝑒𝛽𝐽+𝛽ℎ+𝑒𝛽𝐽−𝛽ℎ=𝑒𝛽𝐽(2cosh𝛽ℎ) Now, we calculate (𝑇𝑟(𝑇))2− 4det (𝑇): (𝑇𝑟(𝑇))2−4det(𝑇)=4𝑒2𝛽𝐽cosh2𝛽ℎ−4(𝑒2𝛽𝐽−𝑒−2𝛽𝐽) =4[𝑒2𝛽𝐽{cosh2𝛽ℎ−1}+𝑒−2𝛽𝐽] Using cosh2𝑥−sinh2𝑥=1, we get: √(𝑇𝑟(𝑇))2−4det(𝑇)=2√𝑒2𝛽𝐽sinh2𝛽ℎ+𝑒−2𝛽𝐽 The eigenvalues are: 𝜆±=𝑒𝛽𝐽cosh(𝛽ℎ)±√𝑒2𝛽𝐽sinh2𝛽ℎ+𝑒−2𝛽𝐽 The partition function is: 𝑍=𝑇𝑟(𝑇𝑁)=𝜆+ 𝑁+𝜆− 𝑁 Now, we set the external field ℎ=0, and using cosh(0)=1 and sinh(0)=0, we get: 𝜆+=𝑒𝛽𝐽+𝑒−𝛽𝐽=2cosh𝛽𝐽 and 𝜆−=𝑒𝛽𝐽−𝑒−𝛽𝐽=2sinh𝛽𝐽 Thus, the partition function for ℎ=0 is: 𝑍=𝜆+ 𝑁+𝜆− 𝑁=(2cosh𝛽𝐽)𝑁+(2sinh𝛽𝐽)𝑁 For 𝑁→∞ and substituting 𝛽ℎ→𝑥, we have: 𝑚=1 𝑁𝜕ln𝜆+ 𝑁 𝜕𝑥 =𝜕ln𝜆+ 𝜕𝑥 𝑚=1 𝜆+𝜕 𝜕𝑥{𝑒𝛽𝐽cosh𝑥+√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽} 𝑚=1 𝜆+{𝑒𝛽𝐽sinh𝑥+ 𝑒2𝛽𝐽sinh𝑥cosh𝑥 √𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽} 𝑚= 𝑒𝛽𝐽sinh𝑥 𝑒𝛽𝐽cosh𝑥+√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽 +𝑒2𝛽𝐽sinh𝑥cosh𝑥 √𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽{𝑒𝛽𝐽cosh𝑥+√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽} 𝑚= 𝑒𝛽𝐽sinh𝑥√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽+𝑒2𝛽𝐽sinh𝑥cosh𝑥 √𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽{𝑒𝛽𝐽cosh𝑥+√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽} 𝑚= sinh𝑥{𝑒𝛽𝐽√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽+𝑒2𝛽𝐽cosh𝑥} 𝑒𝛽𝐽cosh𝑥√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽+𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽 𝑚= sinh𝑥{𝑒𝛽𝐽√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽+𝑒2𝛽𝐽cosh𝑥} √𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽{𝑒𝛽𝐽cosh𝑥+√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽} 𝑚= 𝑒𝛽𝐽sinh𝑥{𝑒𝛽𝐽cosh𝑥+√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽} √𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽{𝑒𝛽𝐽cosh𝑥+√𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽} 𝑚= 𝑒𝛽𝐽sinh𝑥 √𝑒2𝛽𝐽sinh2𝑥+𝑒−2𝛽𝐽 Putting back 𝑥→𝛽ℎ, we finally get the exact analytical solution for the magnetization: 𝑚= sinh𝛽ℎ √sinh2𝛽ℎ+𝑒−4𝛽𝐽 We now compare the magnetization 𝑚 of the 1D Ising model obtained from the mean-field approximation and the exact transfer-matrix solution as a function of temperature for different external fields ℎ. The mean-field values of 𝑚 have been obtained by fixed-point iteration of 𝑚=tanh(2𝐽𝑚+ℎ 𝑇), while the exact results were computed directly from the analytic expression for 𝑚 derived above. Figure 9: A comparison of exact and mean-field plots of |𝑚| vs 𝑇 for different values of ℎ (1D Ising model) We see that the mean-field curves (dashed lines) decay more slowly with temperature because the mean-field approach neglects thermal fluctuations and spin-spin correlations, effectively assuming that each spin interacts with an averaged background field. Thus, it is no surprise that the overall order (magnetization) predicted by the mean-field approximation is an overestimation. In contrast, the exact 1D solution (solid lines) correctly shows that magnetization diminishes more rapidly since fluctuations destroy long-range order at any finite temperature 𝑇>0. For ℎ=0, the exact magnetization remains strictly zero. Algebraically, this is because sinh𝛽ℎ=0 at ℎ=0; physically this reflects the unbroken Z2 symmetry of the system – configurations with spins −1 and +1 are equally probable, and their net contribution is zero. The system as a whole can be thought of as flipping between both the spin states infinitely often, and does not prefer a direction overall. The mean field approach ignores these fluctuations, and the divergence between the exact curve and the mean-field curve for ℎ=0 in the above plot can be attributed to the fact that we started with a nonzero value of the initial magnetization. In the plot below, we compare the magnetization 𝑚 as a function of external field ℎ at different temperatures, for both the exact 1D Ising solution (solid lines) and the mean-field solutions (dashed lines). At low temperatures, the mean-field curves show an abrupt jump in magnetization near ℎ=0, reflecting how mean field theory predicts the system to be strongly ordered and capable of switching collectively between the two states. The exact 1D solution, however, changes smoothly with ℎ, since thermal fluctuations prevent any truly collective flip in one dimension. The disparity between the two theories decreases at higher temperatures – when thermal agitation dominates over interactions and the spins are effectively independent. The slight horizontal shift or asymmetry of the mean-field curves at low temperatures can be attributed to the effect of the initial magnetization bias on the fixed-point iterations. Notice that the exact solution remains symmetric about ℎ=0, reflecting the unbroken Z2 spin-flip symmetry of the system. Figure 10: A comparison of exact and mean-field plots of 𝑚 vs ℎ for different values of 𝑇 (1D Ising model) Monte Carlo Simulations Statistical Physics Statistical mechanics tries to understand large (macroscopic) systems in terms of the dynamics of the microscopic constituents, and with the help of probabilistic assumptions. Although we can use quantum mechanics to study the behavior of a single microscopic particle (within the limits set by the uncertainty principle), it becomes practically impossible to apply quantum mechanics to explain the behavior of large, classical systems. Statistical mechanics makes it possible to effectively understand the macroscopic properties observed in such large systems in terms of the dynamics of the microscopic constituents. A major portion of statistical mechanics is concerned with systems in equilibrium, simply because it is much easier to study systems that are in equilibrium with their surroundings. However, most of the systems around us are not in perfect equilibrium. In fact, if everything around us was in perfect equilibrium with the surroundings, they could not exchange energy and information with the surroundings, and nothing would occur. It would be a very boring world! In classical mechanics, we study the trajectory of a particle – the path it follows through space over time. Let us consider a particle in a one-dimensional space (represented by the 𝑥 axis). At every point in time 𝑡, the particle will be at some position 𝑥. The state of the particle can be represented by a point on position-vs-time space. If we plot the points corresponding to all the states of the particle as it moves through space over time, we will trace out the trajectory of the particle. In statistical mechanics, however, it is more useful to study the dynamics of a system in phase space. Consider the same scenario: a particle moving in a one-dimensional space. For this particle, the phase space will be a two-dimensional momentum-vs-position space. Every point on the phase space represents a microstate of the particle. To completely determine the state of a particle, we need to know its position and momentum (or velocity, since momentum can be calculated by multiplying the velocity with the mass of the particle) at each point in time. For a single particle, we can obtain this information just by looking at the trajectory of the particle in real space. The trajectory is a position-vs-time curve, the slope of which gives the velocity of the particle. However, when we have a macroscopic system comprising a large number of microscopic particles in three dimensions (the kind of systems statistical physics deals with), it is extremely difficult to calculate the states of the microscopic particles (microstates) of the system as it evolves, if we work in real space. Phase space is defined in a way that each point on the phase space corresponding to a particle automatically contains information about the position and momentum of that particle, and effectively represents the complete state of the particle. Of course, the phase space is an abstract mathematical space, and is not necessarily three dimensional. However, it is an extremely powerful formulation when dealing with systems made of a large number of particles. Now, we briefly discuss the idea of microstate and macrostate in statistical mechanics. Let us consider a box that is divided into two partitions; let the left partition be A and the right partition be B. We also have four distinguishable particles colored red, green, blue and yellow. Say we want a configuration in which there are three balls in the left partition A and one ball in B. We can achieve this configuration in four different ways, as shown below. Figure 11: Four microstates corresponding to the same macrostate (Source: Chemistry LibreTexts) Let us assume the color is a microscopic property of the balls, which cannot be directly measured by us. In other words, macroscopically the four cases shown in the above figure are the same. This defines the macrostate of the system: three balls in A, one in B. If we talk about the microstate, however, the four cases above are distinct. Thus, there can be more than one microstate corresponding to the same macrostate. For a real system, like a box containing a large number of gas particles, examples of macroscopic properties are temperature and volume. A given value of temperature, at a given instant of time, is a macrostate of this system. Macroscopic properties can be measured by us, and they are defined on the entire system. It is meaningless to ask: What is the temperature of that particular gas molecule? We can only ask: What is the temperature of this box of gas right now? Using statistical mechanics, we can show how these macroscopic properties arise from the behavior of the microscopic particles in the system. Temperature, for instance, is a collective effect of the vibration of the microscopic particles; the more vigorously the particles vibrate, the more would be the temperature of the system. This depends on the microscopic properties of the particles, like their positions and momenta (in other words, the microstate). We now briefly discuss the concept of entropy, which plays a central role in statistical physics. It is commonly said that the more the disorder, the more the entropy. However this simple definition of entropy can be deceptive without a deeper understanding. Entropy can be defined as the number of ways a particular state can be achieved. This does not mean just the number of ways the different parts of the system can be arranged (spatial configurations). We also have to take into account the velocities and energies of all the parts of the system at a given point in time. According to the second law of thermodynamics, the entropy of a closed system (a system that does not exchange energy and information with its surroundings) must always either remain the same or increase. In our universe, the overall entropy is increasing with time. It might seem that the entropy increases in our universe because the universe is expanding and more and more spatial configurations become possible within the volume of the universe due to the expansion. However, entropy does not depend solely on the number of different spatial configurations possible. Even if the universe were contracting, the overall entropy could have increased due to the fact that every object in the universe could emit radiation, which contains energy, thus increasing the number of possible energy configurations of the universe, since the radiation (or photons) would spread throughout the universe and disperse energy everywhere. Equilibrium refers to the state with the maximum entropy. It is a natural tendency of systems that are not in equilibrium to try to reach an equilibrium state (highest entropy state) over time. If we forcibly try to prevent a system from reaching equilibrium (the highest entropy state), we must be aware of the states of the particles of the system to achieve this feat (the information of the states in our brains can be thought of as a factor that contributes to the increase of the overall entropy), and we must supply some external energy (do some work) to the system to prevent it from attaining an equilibrium state. This system can no longer be regarded as a closed system, since there is some energy input from the surroundings. The second law of thermodynamics states that the entropy of closed systems must not decrease; however, for open systems, a local decrease in entropy is perfectly allowed, as long as the total entropy of the system plus surroundings increases. A common scenario often employed to illustrate the idea of entropy is that if we have an untidy room, it is in a high entropy state, whereas a neat and clean room will have a low entropy. It is perfectly possible for us to clean the room and decrease the entropy of the room, but we have to spend a lot of energy to do the cleaning, which contributes to the increase of the overall entropy. Now let us properly define entropy. We may define entropy as a measure of disorder in a nontechnical fashion, but that is not exactly so. It is true that an untidy room is in a higher entropy state as compared to a neat and clean room, but this is not a very accurate description of the concept of entropy. Let us consider a better example. If we have five differently colored balls and two boxes, the lowest entropy state would be achieved if we keep all five balls in one box and thus, no balls in the other. There are only two ways in which we can achieve this configuration: either put all the balls in the first box keeping the second box empty, or vice versa. This indicates low entropy. Also, this state is clearly not an equilibrium configuration, because the “concentration” of the balls is higher on one side. If we decide to keep two balls in one box and three in the other, by taking into account the five different colors, the number of ways we can achieve this state goes up drastically, indicating a higher entropy. And in this case, a near-equilibrium state is attained. This is a beautiful illustration of how increasing entropy corresponds to the system moving toward equilibrium. Figure 12: The ordered arrangement on the left becomes disordered over time due to increase in entropy (Jack Westin) Since everything tends to reach an equilibrium state, we may argue (within the limits of statistical mechanics) that the probability of finding a system in equilibrium is the maximum. And since entropy increases over time as the system evolves, and the equilibrium state has the maximum entropy, it seems there must be a positive correlation between entropy of a state and the probability of finding a system in that state (and if this state is the equilibrium state, the corresponding entropy is maximum). Based on this assumption, Boltzmann arrived at a mathematical expression for the entropy of a state that we find the system in: 𝑆=𝑘ln(𝑊)+𝐶 Here 𝑆 is the entropy of the state, 𝑊 is the number of possible ways we can arrive at this state (more precisely, the number of microstates corresponding to the given macrostate), and 𝑘 and 𝐶 are constants. Based on further arguments, it was showed that 𝐶=0 and the value of the constant 𝑘 was also determined. It is now known as Boltzmann’s constant 𝑘𝐵. Thus, we can write: 𝑆=𝑘𝐵ln(𝑊) Although 𝑊 is sometimes referred to as “thermodynamic probability,” it is clearly not a fraction between 0 and 1, and is nothing like ordinary probability. The more the number of possible ways allowed for the state we are concerned with, the higher the value of 𝑊 (it is more probable to find the system in that state), and the higher the value of the corresponding entropy 𝑆. It should be noted that the above equation is applicable only for systems in equilibrium or near-equilibrium with their surroundings. From 𝑆=𝑘𝐵ln(𝑊), we see that the thermodynamic probability (number of allowed microstates) associated with a given macrostate is: 𝑊=𝑒𝑆/𝑘𝐵 In classical thermodynamics, entropy (more precisely, the change in entropy 𝑑𝑆) is defined as: 𝑑𝑆≡𝛿𝑄𝑟𝑒𝑣 𝑇 where 𝛿𝑄𝑟𝑒𝑣 is the heat exchanged in a reversible process. Although 𝛿𝑄 is not easily measurable since it depends on the path taken by the process between states, the ratio 𝛿Q/𝑇 (for reversible process) is path independent, and depends only on the initial and final states. This property allows entropy 𝑆 to be defined as a thermodynamic state variable. The connection to energy follows from the first law of thermodynamics (which is basically energy conservation): 𝑑𝑈=𝛿𝑄−𝛿𝑊 Here, 𝑑𝑈 is the change in internal energy, 𝛿𝑄 is the heat exchanged and 𝛿𝑊 is the work done during the process. For a simple compressible system, the mechanical work is 𝛿𝑊=𝑃𝑑𝑉 (𝑃 is pressure, 𝑑𝑉 is the change in volume of the system during the process), so the heat exchanged can be written as: 𝛿𝑄=𝑑𝑈+𝑃𝑑𝑉 In many statistical-mechanical models, including the Ising model, the volume is fixed (𝑑𝑉= 0), which reduces the first law of thermodynamics to simply 𝛿𝑄=𝑑𝑈, and the entropy change reduces to: 𝑑𝑆=𝑑𝑈 𝑇 This relation shows that changes in entropy are directly linked to changes in internal energy when the system is held at constant volume. Using 𝑊=𝑒𝑆/𝑘𝐵, an increase in energy 𝑑𝑈 leads to a change in the number of accessible microstates that scales as: 𝑑𝑊∝e𝑑𝑈/𝑘𝐵𝑇 Now, we consider the probability of a specific microstate with energy 𝐸 of a system, when the system is in thermal contact (can exchange heat) with a heat reservoir at temperature 𝑇. In this case, the total entropy is the sum of the entropy of the system and that of the reservoir. When the system gains an energy 𝐸, the reservoir loses the same amount of energy, and its entropy decreases by: Δ𝑆𝑟𝑒𝑠=−𝐸 𝑇 We have seen that the expression for the thermodynamic probability for a given state is 𝑊= 𝑒𝑆/𝑘𝐵. Here, the thermodynamic probability 𝑃(𝐸) of the system being in a given microstate with energy 𝐸 would be proportional to the number of accessible microstates of the reservoir with entropy change Δ𝑆𝑟𝑒𝑠=−𝐸/𝑇 (since the system microstate is fixed). Thus, we get: 𝑃(𝐸)∝𝑒−𝐸/𝑘𝐵𝑇=𝑒−𝛽𝐸 where 𝛽=1/𝑘𝐵𝑇. This exponential suppression of higher-energy configurations is the Boltzmann weight, which underlies the canonical ensemble and governs the equilibrium statistics of the Ising model and many other systems. Ergodicity and the Postulate of Equal A Priori Probabilities A fundamental assumption in statistical mechanics is that every possible microstate is equally likely for a system in equilibrium. Fundamentally, we cannot assume a particular microstate to be more probable than any other microstate without precisely explaining this preference. For a system in equilibrium with its surroundings, we do not expect to see any particular possibility being preferred over any other possibility at a fundamental level. Essentially, we feel every possibility must occur in reality somewhere at some point in time – something the essence of which is beautifully captured in the popular saying: “Everything not forbidden is compulsory.” In statistical mechanics, we have two extremely fundamental ideas that we now briefly discuss: ergodicity and the postulate of equal a priori probabilities. Let us begin with ergodicity. According to the ergodic hypothesis in statistical physics, the time average of a macroscopic parameter of a system in equilibrium must be equal to the ensemble average of the same parameter of the same system in equilibrium. Suppose we want to measure a macroscopic parameter A of a box of gas in equilibrium with its surroundings. Let us construct a hypothetical scenario in which we measure A for this system a very large number of times, and at infinitesimal intervals of time. Then we calculate the average of all these values. This is the time average of the parameter A. Next, let us consider another hypothetical picture in which we create a very large number of identical copies (an ensemble) of this system in space. When we say “identical copies,” we mean all these copies must be in the same macrostate. As we have already discussed, there can be many possible microstates corresponding to the same macrostate; thus all these identical copies of our system are allowed to be in different microstates, as long as they correspond to the same macrostate. Now imagine we measure the parameter A of all these copies of the system at once, and calculate the average of all the obtained values; this gives the ensemble average of A for this system. According to the ergodic hypothesis, both these average values (the time average and the ensemble average) should be equal. One might wonder what is the need of even considering the ensemble average, since we have said that all the “identical copies” have the same macrostate, and we are measuring a macroscopic parameter A. Should all these “identical copies” not have the same value of A? It is important to realize that certain macroscopic parameters of a system are functions of the underlying microstate, and can vary across microstates corresponding to a macrostate that is defined using some other macroscopic parameter. Ensemble-averaging becomes important when we are measuring such a parameter. For example, say we have a box of gas where the energy is fixed. There are many possible microstates corresponding to this same energy macrostate. We, however, are interested in measuring the pressure. The pressure of this system will be different for these different microstates, even though all of them correspond to the same energy. In some microstates, the particles will be colliding with the walls of the box more frequently; in some other microstates, the particles will be colliding among themselves in the middle of the box more frequently, and so on. The pressure would be different for each microstate, and the average pressure can be determined by averaging the pressure measurement outcomes obtained for each microstate (ensemble average). Let us now think about the implications of the apparently obvious-sounding claim that the time and ensemble averages must be equal. When we consider the identical copies (ensemble) of the system, we take into account all the possible, unique microstates of the system. And this number is very large (essentially infinite). Imagine all the possible combinations of position and momenta of all the microscopic gas molecules in a box of gas made of a large number of molecules; this is a huge number (again, essentially infinite). This means in the ensemble, we have one copy of each microstate corresponding to the given macrostate. The ensemble average is the average of all the values obtained for each microstate, counted only once. And when we are making a large (infinite) number of measurements at very small intervals of time to calculate the time average, we get the same result (according to ergodicity). This means as the system evolves over a sufficiently long period of time, it accesses all the possible, unique microstates (that correspond to the given macrostate). And importantly, the system will never revisit a microstate it has already visited without first visiting all the available microstates at least once. Considering the number of microstates to be infinite, this implies that it is equally likely to find a system in equilibrium in any one of the possible microstates that correspond to the given macrostate. This, in fact, is the postulate of equal a priori probability, which implicitly assumes ergodicity. Principles of Monte Carlo Simulations Monte Carlo methods are extremely useful to effectively simulate statistical systems by sampling high-dimensional configuration spaces and generating a representative ensemble of only those configurations that contribute significantly to the average value of an observable. This way, it becomes possible to study the average behavior of the system without having to solve it analytically (which is often impossible for large and complex systems). Ergodicity is one of the principles on which Monte Carlo simulations work. When we run a Monte Carlo simulation, we do not sum over every possible configuration explicitly. Instead, we walk through the configuration space using random updates that accept or reject moves according to some specified rule. The goal is to generate configurations 𝐶 with energies 𝐸(𝐶) distributed according to the Boltzmann weights: 𝑃(𝐶)∝𝑒−𝛽𝐸(𝐶) Then we measure observables as averages over the sampled configurations (we assume the sampled configurations are uncorrelated); the averages represent the full ensemble average. For this to work, it must be possible to reach all relevant configurations on the configuration space given enough time. In context of Monte Carlo simulations, “time” usually refers to number of simulation steps, and not real time. The idea is that given enough (but finite) number of Monte Carlo steps, the random walk on the configuration space must be able to visit any configuration (microstate) with a nonzero Boltzmann weight. This is ergodicity. We now illustrate the basic idea of Monte Carlo simulations. We iteratively propose a small, random change in a configuration 𝐶𝑖, resulting in a trial configuration 𝐶𝑖+1 𝑡 that is either accepted or rejected. This is a Monte Carlo move. If it is accepted, the trial configuration becomes the new configuration: 𝐶𝑖+1=𝐶𝑖+1 𝑡. If it is rejected, the initial configuration does not change: 𝐶𝑖+1=𝐶𝑖. Notice that in this scheme each step depends only on the previous step (and not the steps prior to that), and the resulting set of configurations for 𝑖=1,2,…,𝑛 is referred to as a Markov chain in the configuration space. The Markov process is described by the master equation: 𝑑𝑃𝐴(𝑡) 𝑑𝑡 =∑[𝑃𝐵(𝑡)𝑊(𝐵→𝐴)−𝑃𝐴(𝑡)𝑊(𝐴→𝐵)] 𝐴≠𝐵 Here, 𝑃𝐴(𝑡) is the probability of the system being in configuration 𝐴 at time 𝑡 (in general this is a function of time) and 𝑊(𝐴→𝐵) is the transition rate from configuration 𝐵 to configuration 𝐴. In the master equation, we look at the time evolution of 𝑃𝐴(𝑡), and this is the difference between 𝑃𝐵(𝑡)𝑊(𝐵→𝐴) and 𝑃𝐴(𝑡)𝑊(𝐴→𝐵) summed over all possible configurations 𝐵≠ 𝐴. Notice that 𝑃𝐵(𝑡)𝑊(𝐵→𝐴) is the transition rate of the system moving to configuration 𝐴 from a configuration 𝐵≠𝐴 at time 𝑡, multiplied by the probability of the system actually being in such a configuration 𝐵 at time 𝑡. This favors the configuration 𝐴, and thus has a positive effect on 𝑑𝑃𝐴(𝑡)/𝑑𝑡. The term 𝑃𝐴(𝑡)𝑊(𝐴→𝐵) captures the transitions away from 𝐴, and impact negatively on 𝑑𝑃𝐴(𝑡)/𝑑𝑡. The master equation is simply saying that the time evolution of the probability of the system to be in configuration 𝐴 at any given time is the difference between the influx to and the outflux from the configuration 𝐴 on the configuration space. We have the conditions 𝑊(𝐴→𝐵)≥0 and ∑𝑊(𝐴→𝐵)=1 . The first condition ensures that the Markov chain on the configuration space is not stuck in configuration 𝐴, and always has some nonzero probability of moving to a different configuration that is physically allowed (for instance, that does not violate any conservation law). Ergodicity is defined only on accessible configurations; for inaccessible configurations, 𝑊(𝐴→𝐵) can technically be exactly zero. The second condition ensures that the sum of all the possible transitions from 𝐴 to 𝐵≠𝐴 is 1, which makes sense because 𝑊 represents a transition probability (more precisely, transition probability per unit time, or the transition rate). 𝑊(𝐴→𝐵) can be further decomposed into a trial probability 𝑇(𝐴→𝐵) corresponding to a proposed change of the configuration 𝐴 to 𝐵, and an acceptance probability 𝐴(𝐴→𝐵) of this proposal. Thus, we may write: 𝑊(𝐴→𝐵)=𝑇(𝐴→𝐵)𝐴(𝐴→𝐵) statistically possible for longer chains, and independent runs can end up in different configurations. This increases the standard deviation in the measured magnetization. We also observe that for longer chain lengths 𝑁, on average the magnetization magnitude |𝑚| becomes smaller for the entire temperature range. This is because in one dimension, domain walls (boundaries separating regions of opposite spin alignment) cost only a finite energy, and so the probability of spins flips increases with the system size 𝑁. As 𝑁 grows, the likelihood of multiple domains appearing within the same chain rises, reducing the net magnetization even though local correlations remain strong. Hence, larger systems (for example the green curve, 𝑁=128) exhibit weaker overall order and the corresponding |𝑚| values lie below those of shorter chains (like the blue curve, 𝑁=32). The system sizes were chosen to be powers of two (𝑁=32,64,128) for numerical convenience – it is not a mathematical requirement. Now, we plot the magnetization 𝑚 as a function of the external field ℎ (for a fixed value of temperature 𝑇=1.5) for different lengths of the chain (𝑁). We observe that the magnetization is nearly identical for different chain lengths because the external field explicitly breaks spinflip symmetry and suppresses spontaneous reversals (hence the number of lattice sites where a spin flip is possible – and hence the length of the chain – becomes irrelevant). In contrast, the |𝑚|-vs-𝑇 curves at zero field showed finite-size dependence since longer chains can host more domain walls. Figure 15: 𝑚-vs-ℎ plot (at 𝑇=1.5) for different lengths of 1D Ising chain using Metropolis algorithm We now choose one specific case and look at plots of 𝑚 as a function of 𝑇 (for a fixed ℎ) and then as a function of ℎ (for fixed 𝑇) for the 1D Ising model using three methods: • Exact analytical solution (transfer matrix) • Mean field approximation (fixed point iteration) • Monte Carlo simulations (Metropolis algorithm) Figure 16: |𝑚|-vs-𝑇 plots (for ℎ=0.2 and 𝑁=128) using three methods (1D Ising model) Figure 17: 𝑚-vs-ℎ plots (for 𝑇=1.5 and 𝑁=128) using three methods (1D Ising model) The plots of magnetization 𝑚 versus temperature 𝑇 and field ℎ show the expected behavior for the 1D Ising model. Across both plots, the mean-field approximation deviates noticeably from the exact transfer-matrix results. In contrast, the Metropolis Monte Carlo data align almost perfectly with the exact analytical curves, confirming the accuracy and reliability of the approach. Space-Time Evolution of 1D Ising System We now briefly study the evolution of the 1D Ising system in space-time. The plots below show space-time snapshots of the 1D Ising model, simulated using the Metropolis algorithm for 𝑁= 100 spins and 100 Monte Carlo sweeps. Each snapshot is a two-dimensional map: the horizontal axis represents the position of the spins along the one-dimensional lattice, while the vertical axis (increasing downward) shows the progression of Monte Carlo sweeps – effectively, the evolution of the spin configuration over time. Each pixel therefore corresponds to the state of a particular spin at a particular iteration. The color coding is discrete: yellow denotes a spin 𝑠𝑖=+1, while green represents 𝑠𝑖=−1. At low temperatures (here we put 𝑇=0.1) the 𝐽/𝑇 ratio is high (since 𝐽=1 and 𝑇=0.1, 𝐽/𝑇=10), which means the coupling constant 𝐽 dominates thermal agitation, favoring alignment of neighboring spins. This causes the system to evolve toward an ordered state, where large, uniform domains of the same spin orientation form and persist over many sweeps (see the plot below). Thermal fluctuations are insufficient to frequently flip spins, so magnetization remains roughly constant in time, reflecting the collective alignment of spins into one of the two possible orientations. Figure 18: Space-time evolution of 1D Ising system at a low temperature under Metropolis algorithm At intermediate temperatures (for the plot below we used 𝑇=1) the 𝐽/𝑇 ratio is intermediate (since 𝐽=1 and 𝑇=1, 𝐽/𝑇=1), and the coupling constant 𝐽 and thermal agitation are of comparable strength (for an appropriate choice of 𝐽). In this regime, neither the tendency for neighboring spins to align nor the randomizing influence of temperature dominates. As a result, the system fluctuates dynamically between local order and disorder, producing a roughly equal distribution of yellow (𝑠𝑖=+1) and green (𝑠𝑖=−1) regions in the snapshot. In contrast, at high temperatures (we used 𝑇=100), the 𝐽/𝑇 ratio is very small (here 𝐽/𝑇= 0.01), and thermal energy overwhelms the spin-spin coupling. Individual spins flip frequently due to strong thermal agitation, destroying any long-range correlations that might form. The result is a disordered configuration where the fractions of −1 and +1 spins become roughly equal on average. This balance corresponds to a magnetization that fluctuates around zero – “uniformly distributed” when viewed over time. The rapid random flipping ensures this distribution remains uniform over time. Figure 19: Space-time evolution of 1D Ising system at an intermediate temperature under Metropolis algorithm Figure 20: Space-time evolution of 1D Ising system at a high temperature under Metropolis algorithm We set the external magnetic field ℎ=0 to isolate the intrinsic collective behavior of the spins from biases introduced by an external field. With no bias, neither spin orientation is favored, allowing the system to spontaneously evolve toward order or disorder purely as a function of the temperature 𝑇 and coupling strength 𝐽. Including an external field would introduce a global bias that drives the spins preferentially in one direction, obscuring the spontaneous ordering behavior. In summary, these space–time snapshots vividly demonstrate how temperature governs the balance between order and randomness in the Ising model (and in many other systems!) – at low 𝑇 interactions dominate and yield large coherent domains, whereas at high 𝑇 entropy wins and randomness takes over. It is extremely important to understand that the space-time snapshots above only illustrate the qualitative evolution of order and disorder in the 1D Ising model under Metropolis dynamics; there is no true phase transition in the 1D Ising model. As we have discussed, in one dimension with short-range interactions, no stable long-range order can exist at any nonzero, finite temperature (no matter however small), and the critical temperature is 𝑇𝑐=0. The apparent “domains” observed at low but finite temperatures are therefore transient and arising from short-range correlations rather than genuine spontaneous symmetry breaking. Only at exactly 𝑇=0 does the system exhibit true long-range order, with all spins aligned in one of the two ground states. A Brief Look at the 2D Ising Model We now briefly discuss the two-dimensional Ising model. In two dimensions, there is a genuine competition between energy (which favors spin alignment) and entropy (which favors disorder), leading to a true phase transition at a finite, nonzero critical temperature. Below this temperature, the system exhibits spontaneous symmetry breaking and stable long-range order, while above it, thermal fluctuations dominate and the system becomes disordered. The 2D Ising model therefore provides the simplest setting in which collective ordering and critical behavior emerge in a nontrivial manner. In the 2D Ising model, we consider spins 𝑠𝑖=±1 located on the sites of a square lattice. Each lattice site can take two possible values: +1 (up spin) or −1 (down spin). Even for relatively small lattice sizes, the number of possible microscopic configurations becomes astronomically large. For example, a lattice of size 𝑁=10 contains 10X10=100 spins, which results in 2100 possible configurations! This highlights the necessity of using statistical and numerical methods, such as Monte Carlo simulations, to efficiently study the macroscopic behavior of the system. Figure 21: Random Ising spins on a 2D square lattice (Jay Shen) Unlike the 1D case, the 2D Ising model exhibits spontaneous magnetization at sufficiently low temperatures. Here, a genuine competition exists between energy (promoting alignment of neighboring spins) and entropy (favoring disordered configurations). Energetic considerations dominate over entropic effects at sufficiently low temperatures, and then the system can spontaneously choose one of the two symmetric ordered states (either all spins up or all spins down), resulting in a nonzero magnetization even in the absence of an external magnetic field. This phenomenon corresponds to a true phase transition at a finite (nonzero) critical temperature. The macroscopic behavior of the Ising model is governed by the minimization of the Helmholtz free energy 𝐹=𝑈−𝑇𝑆, which encodes the competition between internal energy 𝑈 and entropy 𝑆. At low temperatures, the entropic contribution 𝑇𝑆 is small, and the system primarily minimizes its internal energy. For the 2D Ising model, the internal energy 𝑈 arises from spinspin interactions between nearest neighbors and is minimized when large regions of aligned spins are formed (domains). As a result, the system settles into an ordered phase with a nonzero magnetization. At high temperatures however, the entropic term dominates the free energy, and in order to minimize the overall free energy, maximizing the entropy becomes more important than minimizing internal energy. The entropy 𝑆 measures the number of microscopic spin configurations consistent with a given macroscopic state, and is maximized when spins fluctuate independently and randomly. At large 𝑇, the system explores the large number of disordered spin configurations even though they have higher energy, and large regions where spins are aligned are never formed (thermal fluctuations destroy long-range order and the magnetization vanishes). At an intermediate temperature 𝑇=𝑇𝑐 (the critical temperature), neither energy nor entropy fully dominates; their competition leads to a qualitative change in the equilibrium state of the system. This manifests as a phase transition between the ordered and disordered phases, characterized by the spontaneous breaking of spin-flip symmetry and the emergence of critical fluctuations spanning multiple length scales. The Hamiltonian of the 2D Ising model is defined analogously to the 1D case, with nearestneighbor interactions on the square lattice (and an optional external field). Since now we are on a 2D square lattice, each spin interacts with four nearest neighbors (in contrast to the 1D case where each spin has only two neighbors). Specifically, a spin 𝑠𝑖,𝑗 located at lattice site (𝑖,𝑗) interacts with the spins at (𝑖−1,𝑗), (𝑖+1,𝑗), (𝑖,𝑗−1) and (𝑖,𝑗+1). This increased coordination number strengthens the tendency toward collective ordering and plays an important role in allowing a finite temperature phase transition in two dimensions. For a single spin 𝑆𝑖,𝑗, the local contribution to the Hamiltonian (no external field) can be written as: 𝐻(𝑠𝑖,𝑗)=−𝐽(𝑠𝑖,𝑗𝑠𝑖−1,𝑗+𝑠𝑖,𝑗𝑠𝑖+1,𝑗+𝑠𝑖,𝑗𝑠𝑖,𝑗−1+𝑠𝑖,𝑗𝑠𝑖,𝑗+1) This is the energy associated with the initial configuration that depends on 𝑠𝑖,𝑗. Now, if we flip the spin 𝑠𝑖,𝑗→−𝑠𝑖,𝑗: 𝐻(−𝑠𝑖,𝑗)=+𝐽(𝑠𝑖,𝑗𝑠𝑖−1,𝑗+𝑠𝑖,𝑗𝑠𝑖+1,𝑗+𝑠𝑖,𝑗𝑠𝑖,𝑗−1+𝑠𝑖,𝑗𝑠𝑖,𝑗+1) This is the energy associated with the new configuration that depends on the new spin at site (𝑖,𝑗) (now −𝑠𝑖,𝑗). Thus, the energy difference for this move would be: Δ𝐸(𝑠𝑖,𝑗→−𝑠𝑖,𝑗)=2𝐽𝑠𝑖,𝑗(𝑠𝑖−1,𝑗+𝑠𝑖+1,𝑗+𝑠𝑖,𝑗−1+𝑠𝑖,𝑗+1) where 𝐽>0. The total Hamiltonian of the system is obtained by summing this interaction over all lattice sites. To preserve the physical symmetries of the model, we must choose appropriate boundary conditions. In particular, the global Z2 spin-flip symmetry (changing all up spins to down and all down spins to up does not change the physics of the system) and translational invariance (no matter which lattice site we are on, all spins follow the exact same physics) are naturally maintained by imposing periodic boundary conditions. In 1D, this simply meant turning the 1D line into a ring. In 2D, it means turning the 2D square lattice into a toroidal lattice (see the figure below). Figure 22: A torus can be constructed from a flat two-dimensional surface by joining its edges as shown (Sophia Potoczak Bragdon) Without periodic boundaries, spins at the edges of the lattice would have fewer neighbors than spins in the bulk of the lattice, explicitly breaking translational symmetry and introducing boundary effects that can distort bulk behavior. In numerical implementations, periodic boundary conditions are enforced by using modulo arithmetic. For a 2D 𝑁X𝑁 lattice, neighbors are accessed using indices such as (𝑖±1) 𝑚𝑜𝑑 𝑁 and (𝑗±1) 𝑚𝑜𝑑 𝑁, ensuring that spins at one edge of the lattice correctly interact with spins on the opposite edges, thus allowing all lattice sites to be treated on equal footing and ensuring that the simulated system effectively represents an infinite system. We now employ the Metropolis algorithm in the same manner as before to plot the average magnetization as a function of inverse temperature 𝛽 (𝛽=1/𝑘𝐵𝑇, in the simulations we set 𝑘𝐵=1, which corresponds to a choice of units for temperature. Since 𝑘𝐵 is a constant, this rescaling does not affect the qualitative behavior of the curves). We start the simulation from low values of 𝛽 and move on to higher values, which corresponds to a high value of temperature, hence avoiding the problem of long equilibration times. This means we are essentially cooling the system and sampling configurations according to their Boltzmann weights as it cools. The plot below shows the magnetization versus 𝛽 curves for four different lattice sizes (10X10, 20X20, 40X40 and 80X80). In this simulation, we allowed the system to equilibrate for 6400 steps, and collected measurements averaged over 20 independent runs for 10000 steps. Figure 23: |𝑚|-vs-𝛽 plots for ℎ=0 and four different lattice sizes using Metropolis algorithm (2D Ising model) The monotonic rise of magnetization with 𝛽 reflects the system becoming progressively more ordered as thermal agitation weakens. At a critical temperature 𝑇=𝑇𝑐, the system transitions from disordered state (zero or near-zero magnetization) to an ordered state (high magnetization). This is spontaneous magnetization (we have no external field). In the above plot, visually 𝑇𝐶 is between 𝛽=0.4 and 𝛽=0.5. Notice that even at very low 𝛽 (very high temperatures), the measured magnetization is not exactly zero. This is due to finite-size fluctuations (hence the effect is more pronounced for smaller lattices, the blue and orange curves in the above plot). A perfectly zero magnetization would require an unrealistically exact cancellation of spins (for example, a perfectly alternating pattern, where all the four nearest neighbor sites of every +1 spin is occupied by −1 spins, and vice versa). In a finite system however, random thermal fluctuations inevitably produce small local clusters and imbalances, and these residual fluctuations lead to a small but nonzero value of magnetization even deep in the disordered phase. For high temperatures (𝑇>𝑇𝐶, small 𝛽), entropy dominates and spins fluctuate almost independently, producing a disordered phase with magnetization close to zero and no persistent domains. For low temperatures (𝑇<𝑇𝐶, large 𝛽), the interaction energy dominates and neighboring spins tend to align, forming large, coherent domains and a nonzero magnetization; one of the two spin orientations (either +1 or −1) effectively dominates the lattice (up to rare thermal excitations). The sharp crossover region around the transition becomes increasingly steep as the lattice size increases. This is because increasing the lattice size reduces finite-size effects, allowing the system to more closely approach the thermodynamic limit. It is important to keep in mind that in this plot we are looking at the average magnetization of the system. In small systems, thermal fluctuations can affect a large fraction of spins, flipping the system back and forth easily and smoothing out the transition. As the system size increases, relative fluctuations become less important and the system increasingly favors either the ordered or disordered phase (since as system size increases, it becomes increasingly difficult to flip a large group of spins when the system, as a whole, starts favoring order or disorder). As a result, the crossover sharpens and approaches an abrupt transition in the thermodynamic limit. It is also observed that the error bars are larger for the larger lattice curves (red and green in the above plot), and in general are noticeable in the transition region (they are vanishingly small in the ordered and disordered regimes for all lattice sizes). This is because of critical fluctuations and critical slowing down. At 𝑇=𝑇𝐶, the correlation length diverges, which means no matter how large the system is, the correlation length becomes comparable to the system size near the critical temperature, and so effectively the full system becomes strongly correlated as it approaches criticality. As a result, configurations change very slowly, reducing the effective number of independent measurements (and hence making the error bars larger). Figure 24: Space-time evolution of 2D Ising system, left snapshot corresponds to low temperature, middle snapshot corresponds to critical temperature, right snapshot corresponds to high temperature (https://doi.org/10.1371/journal.pcbi.1000314) The correlation length becoming comparable to the system size (for systems of all sizes) near the critical temperature 𝑇𝑐 simply means that the system exhibits fluctuations on all length scales, leading to spin configurations with no characteristic domain size. Near 𝑇𝑐, the system is at the point where energy and entropy balance each other almost perfectly. In other words, creating a small ordered cluster or a large ordered cluster costs nearly the same free energy, because the energetic gain from alignment is almost canceled by the entropic cost. As a result, fluctuations of all sizes become equally likely, and near criticality, small, medium and very large clusters all coexist. In this critical regime, space-time snapshots (see the figure above) of the spin configuration would show complex, interpenetrating clusters of up and down spins with no single characteristic domain size. These are fractal-like patterns, because they are visually self-similar over a large range of scales (since the correlation length of the system diverges at criticality, the system becomes scale-invariant and is neither ordered nor completely disordered). Discussion and Conclusion The 1D Ising model, despite its simplicity, is a powerful pedagogical system for connecting analytical reasoning, computational methods and physical intuition. By combining exact theory with Monte Carlo simulations, we have captured the complete behavior of the model for the chosen observable (which is magnetization here). The excellent agreement between the Monte Carlo and the exact analytical results for the 1D Ising model illustrates the strength of the Metropolis algorithm. It is important to note that the Metropolis algorithm is a local Monte Carlo update scheme, and many more advanced algorithms exist for simulating the Ising model and other systems (for example, algorithms that flip entire spin domains at once and significantly reduce critical slowing down). While the agreement between the Metropolis Monte Carlo and the exact analytical methods is remarkable, the mean-field approximation, while qualitatively correct, overestimates order by neglecting spatial correlations. Finite-size effects also emerge naturally; in 1D, domain walls cost finite energy and therefore long-range order cannot persist in the thermodynamic limit. These subtle effects encapsulate the essence of statistical mechanics – how local rules give rise to macroscopic regularity and its breakdown. We have also briefly touched upon the 2D Ising model, and analyzed the behavior of the magnetization as a function of inverse temperature 𝛽 using the Metropolis algorithm – the emergence of collective ordering and a sharp crossover associated with a true phase transition. The main focus of this report was the 1D Ising model, and although the 1D Ising model is relatively simple compared to more complex, higher-dimensional systems, it establishes a robust foundation for exploring those topics. The methods and reasoning developed in this report – especially the understanding of concepts like detailed balance and ergodicity – can be directly generalized to more complex systems, making this work both a conclusive summary and a conceptual starting point for future studies in statistical physics and computational modeling.