Behavior of Early Warnings near the Critical Temperature in the Two-Dimensional Ising Model

Among the properties that are common to complex systems, the presence of critical thresholds in the dynamics of the system is one of the most important. Recently, there has been interest in the universalities that occur in the behavior of systems near critical points. These universal properties make it possible to estimate how far a system is from a critical threshold. Several early-warning signals have been reported in time series representing systems near catastrophic shifts. The proper understanding of these early-warnings may allow the prediction and perhaps control of these dramatic shifts in a wide variety of systems. In this paper we analyze this universal behavior for a system that is a paradigm of phase transitions, the Ising model. We study the behavior of the early-warning signals and the way the temporal correlations of the system increase when the system is near the critical point.


Introduction
Complex systems are those in which the agents or elements that compose the system interact non-linearly and in such a convoluted way that it is impossible to describe the behavior of the system in terms of the simpler behavior of its components.Intuitively, it should be easy to define precisely what a complex system is and under which conditions complexity emerges.However, so far there is no general agreement of what complex systems are.In fact, the very definition of complexity is a main research topic [1][2][3][4].Complex systems share certain general properties which in some fashion describe them: emergence, selforganization, homeostasis, entangled properties on multiple scales, ability to efficiently transmit and process information, etc. [1,2,4].In this work we focus on a property that is apparently common to many complex systems: the existence of critical thresholds [5].A large number of complex systems display behaviors related to criticality and phase transitions [6][7][8].When a physical system is in a critical point it acquires unique properties, one of which is scale invariance.Theoretically, at that point the correlation length diverges, which in practice means that the correlation length becomes very large when compared to the scales of interaction of the system.These properties may be the key to the robustness and adaptability of complex systems.One of the main properties that most of complex systems share is the fact that they seem to lie on the boundary between order and chaos [9].It is in this regime where complex systems acquire robustness and where they can adapt to environmental changes.If the system has excessive order, it tends to be too rigid to react promptly to the pressures of the environment, and is thus incapable of responding to its requirements.It is not capable of evolving.On the other hand, if the system behaves too randomly, then it is too fragile and any environmental perturbation will impair fundamental properties of the system.Living on the edge between order and chaos, complex systems are capable of being adaptive to changes in a robust way.In the physical sciences, this regime between order and chaos is frequently associated with phase transitions and criticality.Although the relation between robustness, adaptability and criticality in complex systems is quite evident, a second way in which the complex systems are related to the phase transition phenomena is through the appearance of catastrophic shifts.It is known that many complex systems have critical thresholds at which the system changes dramatically from one stable state to another.Usually these changes are abrupt and dramatic relative to the long-range time scales of the system.These shifts are driven by small perturbations (again in terms relative to the typical scales of the system), follow paths with hysteresis, and are related to the critical slowing down phenomenon [5].This kind of transition is present in a large and diverse set of complex systems, which include ecosystems [10][11][12], species population [13], biodiversity, the global economy [14], physiological systems [15][16][17], climate systems [18], social systems, etc.Although the particular features of these systems are quite diverse, all these transitions display characteristic signals which in principle are independent of the particularities of the system.In other words, the dynamics of systems near critical points exhibit universal properties.There has recently been an increased interest in understanding how a complex system behaves in the vicinity of catastrophic shifts, in part to predict and possibly control the timing and evolution of such transitions.Usually early warnings and critical transitions are studied either in model systems, for which it is possible to describe the behavior of the system in terms of a mathematical model, or in actual complex systems in which the description is based in a rich collection of data, but where a full mathematical description is not available [5].In this work we focus on a very well known physical system, the Ising model [19].We choose this model because the phase transition associated to it is very well known and the presence of critical slowing down is a well studied phenomenon [20].For the study of biological systems, the Ising model has been applied to the study of human brain [21][22][23], cancer [24], protein folding [25], ion channels [26], statistical genetics [27] and cardiac activation [28].The Ising model has also been used successfully to study collective phenomena in social systems [29], specially in order to simulate crowd dynamics and opinion formation.It is a system on which we can delicately control the parameters and the location of the critical point, but that is not trivial in terms of its dynamics, i.e., it is not possible to fully understand the system in terms of dynamical equations in a trivial way.The main goal of the present work is to analyze the dynamics of the Ising model in terms of the classic early warning properties described in literature.

The Ising model
The Ising model is a statistical physics model for ferromagnetism.It is paradigmatic both for systems in which cooperative phenomena play an important role and for the study of physical phase transitions.The definition of the system is very simple.Consider a lattice of N sites with a spin state σ defined on each site.We let each of the spins take one of two possible orientation values, denoted by σ = ±1.There are thus 2 N possible configurations of the system.Each of these spin sites interact with its nearest neighbors with an interaction energy given by where the first summation runs only through neighboring spins and J ij represents the coupling strength between spins i and j.If this coupling is positive then the neighboring spins will tend to align parallel to each other, since this minimizes the energy.B i represents the external magnetic field acting on site i and µ is the magnetic moment.In this work we focus only in the case where J ij = constant and B i = 0 (no external magnetic field).The probability that the system is in a given configuration depends on the energy of the configuration and the value of the parameter T , which is identified as the temperature of the system.This probability is given by the Boltzmann distribution where β = (kT ) −1 , k is the Boltzmann constant and the normalization Z β is the partition function.It is possible to measure the order present in the system through the total magnetization, defined as A well-known fact is that if it is defined on a 1-dimensional lattice, in which each spin has only two nearest neighbors, the system will have no phase transition.However, for lattices in 2 or more dimensions the system goes through a phase transition when T is equal to a critical value T c .Below the critical value, the system undergoes spontaneous magnetization and all the spins tend to align towards either the +1 state or the −1 state.For temperatures higher than T c , the system becomes paramagnetic, where the total magnetization of the system is zero on average.The presence and size of clusters of equally aligned spins is also markedly different in these two regimes: when T is lower than T c , large resilient clusters form, while above T c only small clusters can survive momentarily.If the temperature is high enough, all the clusters are completely destroyed.In the critical point (T = T c ), however, clusters are continually formed and destroyed in a wide range of scales, with the distribution of cluster sizes following a power law. Figure 1 shows a typical spatial configuration for each of the three regimes of temperature for a 2-dimensional system.Black squares represent spins with σ = +1 and white ones represent those with σ = −1.
While we are interested in the dynamics of the Ising model, the model described thus far does not incorporate dynamics since there is no kinetic term in the Hamiltonian.A kind of artificial dynamics can be imposed on the system through Monte Carlo simulations using the Metropolis algorithm [30].This algorithm iteratively generates successive spin configurations.While these configurations do not represent the time evolution of a system of spins, it is possible to associate them to such a system in contact with a heat reservoir through Glauber dynamics [31].We will analyze the successive configurations obtained with the simulation assuming that they represent the evolution of a correlated system.We are not interested in the Ising system as a model for ferromagnetism, but rather as a nearest neighbor interaction system in which two competing effects are acting: on the one hand, the short-range interactions between neighbors, while on the other the stochastic fluctuations caused by temperature.We are thus not interested in recovering the usual thermodynamic properties through the simulation.While the critical slowing down phenomenon near the critical transition has been previously seen with the Metropolis algorithm [20], the focus has always been on trying to avoid it since the main goal of Monte Carlo simulations is to obtain independent configurations of the system.The presence of critical slowing down is an indication that the successive configurations obtained through the simulation are correlated, i.e., not independent.Usually, the critical slowing down is estimated through relaxation times and it is considered a deficiency rather than a feature of the system.It is thus typically used to determine how many simulation steps need to be skipped in order to obtain independent configurations.In this work we take an opposite view and consider critical slowing down as a feature of the system which can be used to detect criticality.Based on this we simulate the Ising model through the Metropolis algorithm and consider correlation between successive configurations as an indication of early warnings of an oncoming critical threshold.

The Metropolis algorithm
The goal of the Metropolis [30] Monte Carlo simulation is to generate a large number of different and independent configurations of the system in order to have enough statistical sampling to estimate average values of thermodynamic properties.In this algorithm, new configurations are generated from a previous state using a transition probability which typically depends on the energy difference between the initial and final states.The exact form of the probability comes from considering detailed balance in the master equation for the transition probabilities.Following the Boltzmann distribution, the probability of the system being in a state n is: where E n is the energy of the state, k is the Boltzmann constant, T is the temperature and Z is the partition function.The transition probability from state n to state m is then given by where Given a previous state of the system, defined by a specific set of values for all spin sites σ i , the algorithm proceeds as follows: 1. Randomly choose a site i; An iteration or a simulation time unit has passed when every spin in the system has had a chance to flip.We randomly choose the order in which the spin sites are selected for the process described above.It is common that when the system temperature is near the critical point T c , relaxation times are computed in order to skip several iterations and avoid correlations between successive configurations.In this analysis we keep all the iterations because we are precisely interested in the correlations in the system and how these correlations behave as T approaches T c .Particularly, we analyzed an Ising system defined over a 2-dimensional square lattice so that each spin has 4 nearest neighbors.We choose the units of the system such that T c ≈ 2.27 which means that k = 1 and J ij = 1.The lattice has a finite size of 100 × 100 sites with periodic boundary conditions.Although the finiteness of the lattice can affect the value of T c , we are not interested in simulating the system precisely at the critical point but rather in wider regions near criticality.We simulate the system for three temperature regimes: T < T c , T ≈ T c and T > T c .For temperatures lower than T c we choose as initial condition a configuration where all the spins are aligned with the +1 state.We do this in order to avoid as much as possible the metastable configurations typical of this regime (metastable states are described later).For both regions where T ≈ T c and T > T c , we choose an initial condition where all the initial spin states are randomized, so that the average total magnetization is M = 0. We sampled the system from T = 1.42 to T = 3.12, changing the temperature in intervals of ∆T = 0.05.For each of these temperatures we run an ensemble of 1000 simulations of 5000 iterations each (where the first 1000 iterations are discarded in order to remove the initial transient).We then computed the total magnetization of the system as a function of time and analyzed the properties of the magnetization temporal fluctuations considering them as a time series.It is important to note that we are not dynamically changing the temperature; its value is fixed for every run of the simulation.Thus our system is not being driven towards the critical threshold.This can also be thought as the rate of change of temperature being very small compared to the dynamic scale of the system.Under this assumption, we can consider that we are sampling separate instants of this very slow transition through the critical threshold.Figure 2 shows a typical magnetization time series for each of the three temperature regimes.For each temperature value we computed the ensemble average of several parameters that can function as early warnings in order to better understand how the system behaves when it is near the critical region.

Early Warnings
The variety of early warning signals that have been proposed in the literature is overwhelming.There recently has been an effort to summarize, compare and contextualize the different techniques that have thus far been explored [5,32].It is possible to broadly classify early warning signals in two main categories: the metric-based indicators which essentially detect subtle changes in the statistical properties of the time series, and the model-based indicators which detect changes in the time series dynamic fitted by a reasonable model.In this work we will focus only on the former approach.Among the metric-based estimators we can distinguish two main lines of thought.
The first one is based on the non-stationarity of the time series.If we think of the time series as a stochastic process, it is possible to estimate the probability distribution of events.Theoretically, if the system approaches a critical point then due to critical slowing down several moments of the probability distribution will change.The second moment of the distribution, the variance, will diverge because near a critical threshold a system recovers very slowly from perturbations, which in principle allows the system to drift across the boundaries of different states [33].Depending on the particularities of the system, it is possible that the fluctuations become asymmetric if the system approaches configurations with an unstable equilibrium.This will produce changes in the third moment of the distribution, the skewness [34].Changes in the fourth moment of the distribution are also possible, because near criticality the system will visit extreme states more often [32].
The second approach in the family of metric-based estimators is to analyze changes in the memory of the time series via its temporal correlations.Changes in the temporal correlation of events are strongly linked with the situation in which the system approaches a tipping point.If the system recovers slowly from perturbations, as happens near the critical threshold, then it is expected that the temporal correlations of the system for short time scales will increase.This effect can be calculated by means of the autocorrelation function for small time lags τ , specially through the autocorrelation at lag 1, C(τ = 1) [35].The correlations for long time scales are also modified when the system approaches criticality.These correlations are associated with slow oscillations in the fluctuations and with long range memory effects in the system.This kind of effects are visible when an analysis for all the temporal scales is performed.Some of the effects most commonly used early warnings are the Power Spectral Density (PSD) analysis [18] and the Detrended Fluctuation Analysis (DFA) [36].The PSD is connected to the amount of correlations present in the system.Particularly, it is well known that the PSD of scale invariant series obey a power law [7,8].It is important to mention that in this work we have focused on temporal early warnings.However, given the rich spatial behavior of the Ising model, it would also be interesting to study early warnings in the spatial domain [37,38].It is our intention to analyze this kind of early warnings and their relation with the magnetization cluster patterns formed when the Ising system approaches the critical region in a future work.

Results and Discussion
We carried out simulations of the Ising model dynamics for a range of temperatures.We followed the total magnetization of the system through time as the simulation evolved, and in this way we constructed a magnetization time series for each simulation.We repeated the experiment 1000 times for each temperature in order to estimate the ensemble behavior of the system.For each temperature we computed the ensemble average of the early-warning signals described previously, the first four moments of the probability distribution, the autocorrelation at lag 1 and the power spectral density.
The first moment of the distribution is the mean and corresponds to the average magnetization of the system.Figure 3 shows the ensemble average for the magnetization as a function of temperature (it is important to mention that for low temperatures the global magnetization can also converge towards −1; we have thus in general used only the positive values).It is clear that our simulation follows the usual behavior of an Ising model.The three regimes can be clearly identified in the simulations.a) For T > T c , the stochastic contribution to the fluctuations dominates the neighbor spin interactions.In this regime clusters of magnetization are very small, and practically all spin sites change orientation independently of the state of the neighboring spins.The system is governed by stochastic fluctuations and any attempt to transmit information through the system will fail, since correlations are quickly destroyed by the high stochasticity of the system.b) On the other hand, for T < T c , the main contribution to the dynamics comes from the short range interaction between neighbors.In this regime large clusters of magnetization are formed, which is shown by the fact that magnetization does not quite reach a value of 1. Again, any attempt of transmit information fails, essentially because any spin flip is promptly outweighed by the dominating neighbor interaction.c) Finally, at the critical point, T = T c , the two opposing effects are in balance.As we mentioned before, in this regime magnetization clusters of various sizes are formed and the corresponding size distribution of the clusters follows a power law, which is a very strong indication that the system is spatially scale invariant.The effects of perturbations in this regime are stable and robust.[39] succinctly summarized this phenomenon as follows: "For all other temperatures, one can disturb the system locally and the effect of the perturbation will influence only the local neighborhood.However at the transition temperature, the local distortion will propagate throughout the entire system.The effect decays only algebraically rather than exponentially.Although only 'nearest neighbor' members of the system interact directly, the interaction effectively reaches across the entire system.The system becomes critical in the sense that all members of the system influence each other."The fundamental emergent property of the system near criticality is its capacity to transmit information over scales comparable to the entire scale of the whole system.While the average magnetization of the system is zero, it is located precisely at the frontier between the zero magnetization regime and the magnetized one.
Changes in the higher moments of the distribution constitute important early warnings.Figure 4 shows the ensamble behavior of the temporal variance of the magnetization and the way it changes as a function of temperature.There are two ways to approach the critical value: either reaching it from temperatures below T c , or coming from higher temperatures.It is evident that the response of the system is different in these two cases.At the critical point the variance is maximal, but its increase is quite different depending on which side of the critical point the system is coming from.In the high temperature regime the variance grows smoothly as the temperature diminishes.In this regime the variance will constitute an excellent early-warning signal, for its changing behavior will indicate that the system is approaching a critical region even if we are not aware beforehand of the existence of a criticial point.On the other hand, for the low temperature regime the variance increases abruptly, so it is difficult to use the variance as an early-warning signal when reaching the critical temperature from below.It is also interesting that the ensamble standard error of the mean (defined as the standard deviation divided by the square root of the sample size and shown with error bars in the figures) becomes larger when the system is near the critical point.
The third moment of the distribution, the skewness, is shown in Figure 5.This moment is related to the asymmetry of events in the time series, and we expect that this asymmetry to appear only in the low temperature regime.The reason for this is that in this regime the system has the possibility of becoming trapped in a meta-stable state where the average magnetization is lower than the expected value.These meta-stable states are produced when instead of one large cluster of magnetization, two stable clusters with opposite alignments are formed.A large perturbation is required for the system to spring out of this state.If the temperature is too low, the stochastic perturbations are too weak compared to the neighbor interactions and the system can spend a long time in meta-stabilty.However, near the critical point the stochastic perturbations are large enough to carry the system away from such meta-stability.This effect increases the skewness of the probability distribution.We thus conclude that skewness is an effective early-warning signal if the system is going towards the critical point from below.When approaching T c from higher temperatures it is not possible to produce asymmetries in the distribution and thus the skewness is not sensitive to the approaching criticality.
Figure 6 shows the kurtosis of the distribution.It is clear that this moment is a good early warning for this system regardless of the direction in which the system approaches criticality.Kurtosis behaves quite symmetrically and increases smoothly when coming from either direction.The only difference is whether the distribution is more strongly peaked or not than the reference normal distribution (which has a kurtosis of 3).When approaching criticality from lower temperatures we have a strongly peaked or leptokurtic distribution, while when doing so from higher temperatures we have a flattened or platykurtic distribution.
In terms of correlations, it is well known that simulations of the Ising model with the Metropolis algorithm display critical slowing down.It is expected that early warnings based on correlation estimations are useful for this system.Figure 7 shows the ensemble behavior of the autocorrelation at lag τ = 1 as a function of temperature.As expected the autocorrelation of the system for very short temporal scales increases when the temperature approaches the critical value T c .The increase is not symmetrical; it is faster when the system approaches the critical value from lower temperatures.However, in both cases it is possible to use this information as an early warning.The autocorrelation is almost 1 for temperatures near the critical value, which means that the system configuration of a given iteration is highly dependent on the previous configuration from which it was obtained.The increase of variance and short term correlations is related to the form of the potential driving the dynamics of the system.They are different if we approach the critical threshold from lower temperatures than from higher temperatures because the form of the potential, basically the way in which it flattens out, is different in these two regions.For this short time scale the system exhibits a clear memory effect.While the short term autocorrelation function is a well known early-warning signal, long term correlations are a better representation of the interesting properties of a critical system.The autocorrelation function is related with the PSD through the Wiener-Khinchin theorem, provided that the time series is a stationary random process [40].In order to explore the behavior of the long range correlations present we analyzed the PSD of the system.Figure 8 shows the memory effect for the whole range of scales.In panel (a) we can observe the evolution of the PSD for temperatures smaller than the critical value.The temperature increases from bottom to top, with the critical value corresponding to the topmost curve.Panel (b) shows the corresponding PSD evolution for the high temperature regime.Again, the topmost curve corresponds to the critical temperature, and this time temperature increases from top to bottom.In both panels the PSD curves for the different temperatures have vertically shifted for clarity.It is clear that a power law appears in the PSD at the critical point.Power laws have been previously connected to criticality, specially with temporal scale invariance [7,8].As we have mentioned, at the critical point the Ising model exhibits spatial scale invariance as well as fractal structure in the sizes of the magnetization clusters formed by the system.It is remarkable that the system also displays temporal scale invariance.Temporal scale invariance means that the time series is statistically the same at all temporal scales.This property is related to long range correlations and long range memory in the system [8].As soon as the system's temperature departs from the critical value, either to lower or higher temperatures, the low frequency part of the PSD flattens out.A flat PSD is characteristic of an uncorrelated system, where the fluctuations are white noise.We observe that the flat region of the PSD becomes wider as temperature gets further away from the critical value.For temperatures that are very far from the critical one, we can expect that the PSD will flatten out for all frequencies.For every temperature, there is a crossover frequency at which the power law changes from a power law to the flat profile, an uncorrelated system.This crossover is connected to the temporal scale at which the long range correlations cease to be relevant.States of the system that are separated by a time interval equal or greater than this scale should thus be independent.This behavior is known in the early warning literature as spectral reddening [18].It is a property of the PSD which appears when the system approaches a critical threshold.However, in order to be used as an early warning signal it must first be quantified.The way in which this is usually done is through the spectral exponent or through the spectral ratio between low and high frequencies.However, there are systems for which the power law behavior appears only near the critical threshold, the Ising model being a prime example.We thus consider that the application of the spectral exponent is limited.The crossover frequency where the power spectrum behavior changes is directly related to the range of the temporal correlations, in turn related with the temporal scale at which the system becomes independent from its past states.When a system is near a critical transition, the correlation length increases, and scale invariance appears.In this state the crossover frequency appears at a scale comparable to the scale of the system.Thus, by comparing the scale at which the crossover occurs to the scale of the system, it is possible to estimate how far the system is from a critical transition.It is important to notice that this split in the PSD profile has been reported in other systems in which criticality is a desired property and the appearance of the the two-part spectrum is considered as an indication that the system is losing the properties associated with criticality, i.e., robustness and adaptability.For instance, in physiological time series, and in particular for cardiac interbeat intervals, the PSD exhibits a power law when the heart of the subject is deemed healthy [41,42].The underlying hypothesis is that the heart is a complex system that has evolved to be both robust and adaptable, and it has achieved this by approaching spatial and temporal scale invariance.One can think that these properties allow the heart to be resilient against environmental changes while at the same time being able to adapt to the typical efforts to which a heart is subjected.This resilience and criticality are lost over time and due to illness, and it has been reported that when this happens the PSD stops being a power law and a two part spectrum appears [43,44].
The fact that criticality and critical transitions are not the same spawns an interesting question: what is the connection of the critical point in the Ising model with bifurcations?It is certainly possible to analyze the Ising model in terms of a bifurcation.By considering a mean field approximation it is possible to neglect fluctuations in the Ising model.Under this assumption the average magnetization satisfies the following rate equation: which can be expanded as: This expression is precisely the normal form of a supercritical pitchfork bifurcation, which is typical in physical problems that have symmetry.Under this view it is clear that Figure 3 can be thought as the upper half of a pitchfork bifurcation diagram (if we imagine the symmetric negative magnetization branch), and that the stable fixed point at M = 0 for T > T c becomes an unstable fixed point when T goes below T c.The negative derivative of the free energy can be related to an effective force driving the magnetization dynamics.Expanding the free energy as a power series in the magnetization gives: which is known as the Landau expansion.Below the critical temperature the free energy has two minima and one maximum, while for temperatures greater than the critical value there is a stable minimum at M = 0.According to this analysis, the system is bistable for temperatures lower than the critical value.
When the system acquires a nonzero magnetization the symmetry of the system has been broken.

Conclusions
In this work we have analyzed the behavior of the Ising system in search of the most common earlywarning signals.We have selected the Ising model because it is a far from trivial example of a system in which a very well known critical point is present.When the system is near the critical point it exhibits the critical slowing down phenomenon and our numerical experiments demonstrate the deep connection of this behavior with several early-warning signals in the evolution of the magnetization time series.We concentrate on early-warning signals related to the statistical properties of the time series probability distribution.We have shown that the change in these properties is asymmetrical, i.e., it is not the same when the critical point is approached from different directions.This effect restricts the applicability of some of the statistics-based early-warning signals, depending on the system at hand.On the other hand, we also analyzed early-warning signals that are based on the estimation of the correlation properties of the system.These early-warning signals seem more robust in terms of asymmetry, but the change in the correlations is slower than the change in the statistical properties.It is important to notice that the in the regions far from criticality this system tends toward states in which either the correlations are nil (high temperatures) or trivial (low temperatures), but certainly this should not necessarily be the case for real complex systems which are not under our control.However, we believe that it is enlightening to understand how these early warnings behave under precisely controlled circumstances, specially because the main goal of early warning theory is to detect critical thresholds in systems where there is no a priori knowledge of whether a critical transition is near or even exists at all.Analyzing the correlations behavior by means of the Power Spectral Density in all the temporal scales available to the system lets us interpret the spectrum reddening as a characteristic way in which the scale invariance properties of critical systems disappear when the system departs from criticality.For systems poised at criticality, a power law in the PSD appears, and this power law is an indication that the system possesses scale invariance in the time domain.When the system departs from criticality the PSD splits in two distinct regions.The low frequency regime is characterized by a flat spectrum, while the high frequency domain exhibits a power law.The crossover frequency that divides these two regions depends on how far the system is from the critical point.When the temperature is far away from the critical value, the spectrum flattens out completely, reflecting a total lack of correlations.At criticality, a single, full power law appears, i.e., the system acquires long term memory.We propose that the frequency at which this crossover occurs can be used as an early-warning signal that can indicate to what extent the temporal correlations have been lost, as well as provide a measure of the loss of temporal scale invariance in the system.

2 . 4 .
Calculate the change in energy ∆E if spin site i were to be flipped; 3.If ∆E is negative, then flip the spin of site i.If on the other hand ∆E is positive, generate a uniformly distributed random number between 0 and 1 and flip the spin only if this random number is less than exp(−∆E/kT ); Choose another site and go back to step 1.

Figure 1 .
Figure 1.Spatial configurations in the Ising model.Typical spatial configurations for a 2-dimensional Ising model.Three regimes are shown: a) T < T c , b) T ≈ T c and c) T > T c .Black squares represent spins with σ = +1 and white one correspond to σ = −1.

Figure 2 .
Figure 2. Total magnetization as a function of time in the Ising model.Typical behavior of the total magnetization time series in a 2-dimensional Ising model.Three regimes are shown: a) T < T c , b) T ≈ T c and c) T > T c .It is important to notice the change of scale between plots.

Figure 3 .
Figure 3. Temporal mean as a function of temperature.Ensemble behavior of the mean as a function of temperature.The mean corresponds to the total magnetization of the system.Three regimes are shown, T < T c , T ≈ T c and T > T c .Note that the mean can also approach −1 at low temperatures; we only show here the positive values.

Figure 4 .
Figure 4. Temporal variance as a function of temperature.Ensemble behavior of the variance as a function of temperature.Three regimes are shown, T < T c , T ≈ T c and T > T c .

Figure 5 .
Figure 5. Absolute values of temporal skewness as a function of temperature.Ensemble behavior of the skewness as a function of temperature.Three regimes are shown, T < T c , T ≈ T c and T > T c .

Figure 6 .
Figure 6.Temporal kurtosis as a function of temperature.Ensemble behavior of the kurtosis as a function of temperature.Three regimes are shown, T < T c , T ≈ T c and T > T c .

Figure 7 .
Figure 7. Temporal auto correlation at lag 1 as a function of temperature.Ensemble behavior of the autocorrelation function for lag τ = 1 as a function of temperature.Three regimes are shown, T < T c , T ≈ T c and T > T c .

Figure 8 .
Figure 8. Power Spectral Density as a function of temperature.Ensemble behavior of the Power Spectral Density as a function of temperature.Panel (a) shows the behavior of the PSD for temperatures T ≤ T c .Temperature increases from bottom to top, with T c corresponding to the topmost curve.Panel (b) shows the behavior of the PSD for temperatures T ≥ T c .Temperature increases from top to bottom, with T c corresponding to the topmost curve.