A novel approach to the study of critical systems

We introduce a novel approach to study the critical behavior of equilibrium and non-equilibrium systems which is based on the concept of an instantaneous correlation length. We analyze in detail two classical statistical mechanical systems: the XY model and the Ising model, and one of the prototype models of Self-Organized Criticality: the forest fire model (FFM). The proposed method can both capture the critical behavior of the XY model and the Ising model and discriminate between the nature of the phase transition in the two scenarios. When applied to the FFM, it gives surprising results, suggesting that the model could be critical despite displaying broken scaling in the distribution of cluster sizes.


INTRODUCTION
The concept of criticality is widely used in many disciplines, spanning finance [1][2][3][4], meteorology [5][6][7], neuroscience [8][9][10][11] and physics [12][13][14][15]. A system in a critical state is usually characterized by scale invariance and self-similarity, and by the development of strong instabilities which are caused by the emergence of longrange temporal or spatial interactions. In statistical physics, the term criticality indicates the behavior of a system near a critical point, which is typically associated with a phase transition between two different states. A classic example of phase transition is the behavior of magnets near a critical temperature T c , which separates an ordered state at low temperatures (T < T c ) from a disordered one at high temperatures (T > T c ). When a system is in a critical state, it is highly susceptible to external perturbations, and it is characterized by the emergence of long-range correlations between its constituent components. This high susceptibility is a direct consequence of the self-similarity of the correlation function, which emerges from microscopic interactions and leads to the presence of strong correlations on all scales of the system. Formally, correlations are described by the covariance between two microscopic physical quantities. In statistical physics, the correlation function is usually defined as the difference between the canonical ensemble average ... of the scalar product between two random variables s 1 and s 2 (usually spins or particles) at positions r 0 and r 0 + r and their uncorrelated average product: C(r) = s 1 (r 0 )s 2 (r 0 + r) − s 1 (r 0 ) s 2 (r 0 + r) Introducing the external control parameter X, at a critical point X c and in the thermodynamic limit one expects to find scale-invariance in the correlations, which corresponds to a power-law behavior of the correlation function: Eq. 2 implies that correlations behave in the same way for any arbitrary rescaling of the distance by a factor µ, i.e. if r → µr then one still has C(r|X c ) ∼ (µr) −η ∼ r −η .
The fact that correlations are present at all scales translates in long range correlations and the resulting critical behavior of the whole system. Away from the critical point, correlations typically decay as an exponential function, and the characteristic length of the exponential is referred to as the correlation lengthξ. The typical functional form that is assumed for the correlation function near a critical point is whereξ indicates the typical length over which two agents are correlated and depends on the control parameter of the system X and the system size L. This length is limited by L and diverges in the thermodynamic limit in correspondence of the critical value of the control parameter X c , i.e.ξ(X c , L → ∞) → ∞, giving Eq.2. It is clear then that the correlation length act as a parameter that describes the typical extension of correlations inside a system and therefore represents the most reasonable quantity to look at when one investigates the critical behavior of a physical system. However, it is essential to observe that measuring a diverging correlation length is not enough to determine if a system is in a critical state, because it does not convey any information about the scaling behavior of the system. In other words, one could observe a divergent correlation length even in a system that is not scale-invariant and therefore not critical, as will be discussed in the next sections.
In this paper, we introduce a new method to investigate the critical behavior of a system. This method is still based on the study of the correlation function, but introduces a new correlation length that is no longer a parameter of the system, but a stochastic variable which distribution is able to catch at the same time the scaleinvariance of the system, the asymptotic behavior of the correlation length and the universal properties of the model.

THE INSTANTANEOUS CORRELATION LENGTH FORMALISM
The method proposed is based on the instantaneous correlation length introduced in [16].
For simplicity, we consider models defined on a 2D lattice from which we sample N independent lattice configurations S 1 , S 2 , . . . S N during the time evolution. The classic estimate of the correlation lengthξ goes as follows: for each configuration S t one computes the two-point correlation function C t (r) between two spins s 1 and s 2 at positions r 0 and r o + r. Assuming translational invariance, one has: where the average . . . is taken summing over all the possible pairs of spins and values of r 0 at a time t. Iterating this procedure for different configurations, one obtains an ensemble of correlation functions {C 1 , C 2 , . . . , C N }, which can be used to compute the time-averaged correlation function C(r) for which the following functional form is usually assumed near a critical point: whereξ(X) is the correlation length which depends on a control parameter X. In correspondence of the critical value of the control parameter X = X c , the correlation lengthξ(X c ) diverges in the thermodynamic limit, and the correlation function decays algebraically. Now we introduce the instantaneous correlation length formalism. Assuming that the system size is sufficiently large to give reasonable statistics for the instantaneous correlation function C t (r), one can fit the instantaneous correlation length ξ t using the same functional form that is used in Eq. 5. Doing this, one obtains an ensemble of instantaneous correlation lengths {ξ 1 , ξ 2 , . . . , ξ N }. Each ξ t is a measure of how critical a single configuration S t is. If the system is far from a critical point, one expects ξ t to be always small because the correlation function will decrease exponentially fast. On the other hand, as the system approaches X c there will be an increasing fraction of configurations with a big correlation length, which corresponds to a power-law behavior of the correlation function. Although one expects the ensemble averaged correlation lengthξ(X) and the average instantaneous correlation length ξ to scale in the same way, it is essential to stress the fact that they are two distinct mathematical objects: the first being a parameter of the ensemble averaged correlation function and the second being a stochastic variable. Indeed, the strength of this new approach lies in the fact that we can now use the ensemble of ξ t to compute not only the average correlation length ξ , but also the distribution of correlation lengths P (ξ). P (ξ) is an entirely new physical object and, as we will see, contains plenty of information about the critical behavior of the system under analysis.

THE DISTRIBUTION OF THE INSTANTANEOUS CORRELATION LENGTHS
Using P (ξ), it is possible to determine whether a system is at a critical point or not. This can be done by looking at the conditional probability P (ξ|X), which should become scale invariant in correspondence of X c . Assuming simple scaling, one expects: for ξ and ξ c bigger than a constant lower cut-off ξ 0 . In Eq.6, ξ c represents an upper cut-off that diverges in the thermodynamic limit, G ξ ξc is a universal scaling function and τ is a critical scaling exponent. In general, the upper cut-off scales as ξ c ∼ aL β , where a is a nonuniversal metric factor and β is related to the universal spatial dimension of the observable [17,18]. From Eq. 6 one can compute the n th moment as Imposing normalization (n = 0) one gets τ ≥ 1. If one absorbs the non-universal constant a in the definition of G(u) and assumes that integral in Eq. 7 converges in zero, then in the thermodynamic limit the average correlation length is given by In the following sections, we will study the behavior of P (ξ) in two traditional statistical mechanical systems, the Ising Model and the XY Model, and to one of the prototype models of Self-Organized Criticality, the forest Fire Model. We conclude this section observing that if Eq. 6 holds, then it automatically allows the introduction of the new critical exponent τ .

ISING MODEL
The Ising model is a mathematical model of ferromagnetism that was invented by Wilhelm Lenz in 1920 and solved for the first time in one dimension by Ernst Ising in 1925 [19,20]. The model consists of N interacting twostate spin variables σ i ∈ [−1, 1] which represent adjacent magnetic dipoles. The energy that is associated with a given macro-configuration σ is given by where the first sum is over all pairs of adjacent spins < i, j >, J ij is the interaction strength, and h i is the external magnetic field. The two-dimensional square lattice Ising model was solved in 1944 by Onsager [21] in the case of no external field (h i = 0) and assuming periodic boundary conditions and constant interaction strength along the x-axis (J x ) and the y-axis (J y ). The 2D Ising Model is central in statistical physics because it is one of the simplest statistical models to exhibit a phase transition between an ordered phase (low temperatures) and a disordered phase (high temperatures). In the case of isotropic interactions J x = J y , the critical value of the temperature T = T c that marks the phase transition is given by If one looks at the lattice at different temperatures, it can be noted that the high-temperature phase is characterized by disorder, because the entropy introduced in the system by the temperature destroys long-range correlations, which results in random configurations with roughly half of the spins up and half of the spins down and no emergent complex structures. On the other hand, at low temperatures, most of the spins will be able to align in order to minimize the energy, giving rise to ordered configurations. At the critical point, the correlation function decays as a power-law, and we are in the presence of long-range interactions and the creation of fractal structures, as can be seen in Fig.1. This behavior is reflected by the fact that the correlation length becomes proportional to the system size L and diverges in the limit L → ∞. In terms of the instantaneous correlation length formalism, if we assume that the upper cut-off scales like L, i.e. β = 1, then in order to maintain the linear relationship between the system size L and ξ, we should have τ = 1 in Eq. 8, meaning that the constant of proportionality would be given by the integral of the universal function G(u). Under these assumptions and in the large L limit, Eq. 8 reduces to As τ = 1 implies that the limit lim u→0 G(u) = 0 [18]. In order to verify if the theory is correct, we need to evaluate whether there is a value of the control parameter X = X c such that P (ξ|X c ) becomes scale invariant. For the 2D Ising model, this corresponds to the critical temperature given by Onsager's solution in Eq. 10. In our simulations, we used the Wolff Algorithm [22] in order to reduce the critical slowing down and for each system size L we sampled 10 6 independent configurations to estimate P (ξ) at T c . In the following, we will use ξ 0 = 1 and ξ c = L √ 2 as the reference upper cut-off, as this is the maximum physical distance between two points on a square lattice with periodic boundary conditions. Plotting P (ξ)ξ as a function of ξ ξc , we can perform a data collapse in correspondence of T c (Fig.2). The resulting curve corresponds to the universal scaling function G(u) which according to Eq. 8, can be used to compute the proportionality constant between ξ and L. In Fig. 3 it is shown how the integral of the universal function G(x) converges to a value that is consistent with the estimated gradient of the line ξ = mL. In summary, when applied to the 2D Ising model, our method was able to identify the critical temperature as the T for which P (ξ|T ) becomes scale invariant and to capture the scaling behavior of ξ , which is consistent with the classical theory. In addition to these two wellknown results, we were able to introduce a new critical exponent for the Ising model, i.e. τ = 1, and to relate the rate of growth of the correlation length to the universal function G(u).

XY MODEL
The two-dimensional XY-model is a paricular case of the Heisenberg model, which was introduced in 1928 [23] as a model for ferromagnetism. Similarly to the Ising Model, it consists of a system of spins in a lattice with the difference that the individual spins can rotate in any direction and are not constrained to take only two values. The energy of the model is given by (12) where the first sum is over pairs of adjacent spins < i, j >, J ij is the interaction strength, and θ i is the angle that a spin σ i makes with respect to some arbitrary direction in the lattice plane. As for the Ising model, in our simulation we keep the interaction strengths constant J ij = J and apply periodic boundary conditions. A typical realization of the model is represented in Fig. 4. The two-dimensional version of this model is of particular interest because at high temperatures correlations decay exponentially fast, while at low temperatures they decay with a power-law, even though in both cases the overall magnetization is zero. This peculiar transition is named after Kosterlitz and Thouless who first discovered it in 1973 [24]. The XY-model is a relevant case to discuss in this context because of the behaviour of the correlation length, which diverges even for finite systems at temperatures below the Kosterlitz-Thouless temperature T KT 0.892J [25][26][27]. In the XY model, the two-point correlation function is defined as [28] C(r) = cos(θ(r 0 ) − θ(r 0 + r)) In our simulations we used the Wolff algorithm [22] sampling 10 5 independent configurations to estimate P (ξ) at T KT and used ξ 0 = 1 and ξ c = L √ 2 . As for the Ising model, it is possible to perform a data collapse for P (ξ) in correspondence of T KT and for τ = 1 and ν = 1. Although the Ising model and the XY model share the same exponents, we can observe in Fig. 5 that in the XY model, ξ is able to exceed the system size L. This is in line with the theory, which predicts a pure power-law in two dimensions in correspondence of T KT [29]. The presence of the Kosterlitz-Thouless phase transition and the behavior of the correlations is summarized in Fig. 7, where we plot the conditional probability P (ξ > L) at different temperatures. As one lowers the temperature, the fraction of correlation lengths that exceed the system size goes from 0 to 90%, which corresponds to the pure power-law decay of correlations at T < T KT .

FOREST FIRE MODEL
The last model we consider is one of the prototype models of Self Organized Criticality: the Drossel−Schwabl Forest Fire Model (FFM) [30]. This model is different form the Ising Model and the XY model because it entails a dissipative dynamics and does not have an external control parameter, like temperature, that can be fine-tuned in order to reach a critical state. The dynamic involves the occupation of empty sites on a 2D grid with new trees (planting steps) and the removing of entire clusters of trees (burning steps). The creation of new trees and the removal of clusters results in the typical patchy appearance of the lattice, which is characterized by the presence of patches of different densities (Fig. 8). The way we implement the FFM follows [31][32][33][34] and is concisely summarized by the following pseudo-code:  To estimate P (ξ), we collected 10 6 independent configurations after a transient of 5 · 10 6 burning steps. From Alg. 1 it is clear that two parameters must be considered: the number of trees that one tries to plant θ, and the system size L. In order to reach a critical state one would like to have both L and θ infinitely large, although there is not a clear rule about how to tune θ for a finite system, and in the literature different authors have used quite a large span of θ values for the same systems size L [34,35]. Despite the model being introduced as critical, it was subsequently realized that the observed power-law in the distribution of clusters sizes displayed deviations from perfect scaling for large system sizes [34,35], implying that the model is not critical in the sense of being scale-free [34], and that all proposed scaling laws seem to be just transient [31]. The correlation length was first studied in [36] for systems sizes L and θ up to L = 512 and θ = 2048, finding that ξ ∼ θ ν , with ν = 0.56. The authors also studied the connected correlation function finding ν c = 0.58, and attributed this discrepancy between the two exponents to numerical error. Another estimate for larger system sizes was given in [37], where the authors used up to L = 17408 and θ = 10 4 finding ν = 0.541 and ν c = 0.576 to be statistically inconsistent, and therefore concluding that the model presents two different diverging correlation lengths. This finding points in the same direction as the lack of scaling observed in the distribution of cluster sizes. However, as it was noted in [31], there seem to be small deviations from a power-law in Fig.1 of [37], meaning that the estimate of ν would be unreliable and therefore not suitable to confirm the presence of multiple diverging correlation lengths. Now we want to apply the instantaneous correlation length formalism to investigate whether P (ξ) displays broken scaling as one should expect from a non-critical model. A similar approach was adopted in [16], where the critical exponent was obtained by fitting the tail of P (ξ). However, the tail includes contributions from the universal function G(ξ, L) and therefore that estimate of the critical exponent is spurious.

Critical Behavior in the Forest fire Model
As we discussed in the previous section, it is not clear how to tune the system size L and θ. In previous studies on the correlation length, the standard procedure consisted in keeping the system size L fixed and looking at the behavior of the correlation length as a function of θ [36,37]. Following this approach, it turns out that it is impossible to perform a data collapse for P (ξ), which agrees with the general lack of scaling observed in the literature so far. The same broken scaling can be observed keeping θ fixed and changing the value of L. If we consider the correlation length as a surface in the space of parameters ξ(θ, L), to keep one of the two dimensions fixed corresponds to two different ways of crossing this surface. In particular, increasing the systems size L without a suitable re-scaling of the parameter θ could lead to a different statistical behavior of the system, although most observables like the average density of trees or the average cluster size seem to be quite robust for a wide range of θ at a fixed L. Even though there are infinite ways of coupling θ and L, it is sensible to choose θ L 2 = k for a constant k (k = 10 − 3 in our simulations). In this way, for different system sizes, one tries to plant the same fraction of trees, which seems to be reasonable if one wants to assure statistical consistency at different values of L. This particular path choice is shown in Fig.  9. Surprisingly, coupling the value of θ and L in this way makes a data collapse for P (ξ) possible (Fig. 10), making P (ξ) the first scale-invariant distribution observed in the Forest Fire Model so far. As for the Ising Model and the XY Model, we used ξ 0 = 1 and ξ c = L β √ 2 , but this time we found τ = 1 and β = 1.123 ± 0.038 (Fig.  11), which corresponds to ν = 0.561 ± 0.019 with a 95% confidence bound. This measurement is consistent with the exponents that have been computed for the two-point and the connected correlation lengths in previous studies [36,37]. We conclude this section observing how the broken scaling in the distribution of cluster sizes P (S) found in [34,35] is not affected by the choice of keeping fixed the ratio θ L 2 . This means that the distribution of cluster sizes is not scale-invariant, although the distribution of correlation lengths is scale-free. Therefore, even though the clusters grow in a non-critical and non-scale-free way, there seems to be some global order in terms of the correlations, which is highlighted by the scale invariance of P (ξ). This is a highly non-trivial result and an aspect that surely requires further investigations. Finally, we observe that the correlations in the FFM seem to grow at a higher rate than in the Ising model and in the XY model (β = 1.123). This is likely due to the burning mechanism, which introduces long range correlations in the system as a consequence of the simultaneous removal of sites that belong to the burning cluster.

CONCLUSIONS
The instantaneous correlation length formalism that we have introduced was able to reproduce the well-known results about the critical behavior of the Ising Model and the XY Model, proving that P (ξ) can be used to identify the presence of a phase transition and to estimate the asymptotic behavior of the correlation length. Furthermore, the introduction of P (ξ) allowed us to define a new critical exponent τ , which happens to be equal to 1 for all the three models discussed in the paper. When applied to the Forest Fire Model, this method allowed to identify a coupling of the two parameters L and θ for which P (ξ) is scale invariant. The scale invariance of P (ξ) was unexpected as it is the first scale-free distribution observed in the model so far, and this opens once again the debate about the criticality of the Forest Fire Model. In particular, we observe that the FFM shares the same critical exponent τ = 1 of the Ising model and the XY model but displays an algebraic growth of the correlation length β = 1.12 which could be the reason behind the broken scaling observed in the distribution of cluster sizes P (S). From a theoretical perspective, all systems that present a critical exponent τ = 1 share a very elegant property, namely the fact that constant of proportionality and the system size dependence are described by the integral of the universal function G(u). In the case of the Ising model and the XY model, we found τ = 1 and β = 1. This means that all the details of the two models are contained in the integral of the universal function G(u), which is characteristic of the model under analysis and becomes the only relevant quantity to distinguish between the critical behavior of correlations for the Ising model and the XY model. In the appendix, it is discussed in more detail the relationship between the ensemble correlation length and the instantaneous correlation length, and how it is possible to obtain the classic critical exponent for the correlation length starting from the instantaneous correlation length formalism. Finally, we observe how the presented method could be easily applied to the study of real-world phenomena, such as brain activity or rain precipitation, as the estimate of P (ξ) only requires to collect different images of the system during its time evolution. The study of P (ξ) in real-systems could be a useful tool to assess the scaleinvariance of the systems under examination and to contribute to a more accurate characterization of their critical behavior.
ACKNOWLEDGMENT LP gratefully acknowledges an EPSRC-Roth scholarship from the Department of Mathematics at Imperial College London, the High-Performance Computing facilities provided by the Research Computing Service, and Gunnar Pruessner for very helpful conversations.

AUTHOR CONTRIBUTIONS
Both authors discussed the results of the numerical simulations and contributed to the final version of the manuscript. L. P. performed the numerical simulations and wrote the paper.
Appendix: Critical exponent of the correlation function As it is well known from classical statistical mechanics, the correlation function of the 2D Ising model is characterized by a critical exponentη = 0.25 [29]. It is therefore natural to investigate whether it is possible to recover this critical exponent employing the formalism we have hereby introduced. It is worth to stress the fact that although we assume the same functional form for the instantaneous correlation function and the classic one, the instantaneous values of ξ and η represent two different mathematical quantities with respect to their traditional counterpart. The crucial point is that we expect the standard correlation lengthξ and ξ to scale in the same way, even though the two quantities are defined differently. In particular, ξ is a variable that is related to how correlated a single configuration is, and it is not bounded by the system size L. Regarding the critical exponent η, since it is a constant, it is not expected to scale with the system size, and we expect it to converge to a value that could be different fromη = 0.25 because the two quantities are averaged differently. This is confirmed by our simulations, which show that the distribution of η is not scale-invariant and that the mean value of η tends to η = 0.34 as L increases (Fig.13). However, it is still possible to estimate the ensemble critical exponentη = 0.25 and, at the same time, check the accuracy of our method. In order to do so, one can use the parameters estimated via fit for each configuration i and reconstruct the correspondent correlation function C i (r). If the error that we do in fitting C i (r) is negligible, we should be able to compute the classical correlation function averaging over all configurations, and hence recoverη = 0.25. Indeed, plotting C(r)r 0.25 vs 2r L for different system sizes we can perform a data collapse (Fig.12), meaning that the fitting error is negligible and that we can safely recover the ensemble critical exponentη = 0.25. * l.palmieri16@imperial.ac.uk