# Multifractality can be a universal signature of phase transitions

###### Abstract

Macroscopic systems often display phase transitions where certain physical quantities are singular or self-similar at different (spatial) scales. Such properties of systems are currently characterized by some order parameters and a few critical exponents. Nevertheless, recent studies show that the multifractality, where a large number of exponents are needed to quantify systems, appears in many complex systems displaying self-similarity. Here we propose a general approach and show that the appearance of the multifractality of an order parameter related quantity is the signature of a physical system transiting from one phase to another. The distribution of this quantity obtained within suitable (time) scales satisfies a -Gaussian distribution plus a possible Cauchy distributed background. At the critical point the -Gaussian shifts between Gaussian type with narrow tails and Lvy type with fat tails. Our results suggest that the Tsallis -statistics, besides the conventional Boltzmann-Gibbs statistics, may play an important role during phase transitions.

Phase transitions are ubiquitous in physical systems. A first general theory of phase transitions is the mean field theory, from which in the vicinity of the critical point the free energy can be expanded in a power series of a universal quantity called “order parameter” Plischke2006 (). This quantity fluctuates around zero in a “disordered” phase above the critical point, and is non-zero in an “ordered” phase below the critical point representing the broken symmetry. Unfortunately the predictions from this theory are often not correct when comparing with experiments. To remedy this, much progress has been made in the theory of phase transitions during last fifty years with the landmark of the discovery of some universal scaling laws and the introduction of the renormalization group method Huang1987 (); Binney92 (); MEFisherRMP1998 (). With them the seemingly complex phase transitions in a variety of systems can be quantified by a few critical exponents which represent the self-similar characteristics of systems during transitions. Recently, such ideas of statistical physics have been applied to many complex systems in various fields and have extensively expanded the understanding to these systems KwapienPhysRep12 (). Yet, these studies indicate that in many cases a few scaling exponents representing the self-similarity fail to quantify certain property of systems. Instead, one needs a large number (spectrum) of scaling exponents to clarify the characteristic of systems which represents a higher level of complexity. Examples include the processes of diffusion-limited aggregation (DLA) AmitranoPRB91 (), turbulence and chaos JensenPRL85 (); MuzyPRL91 (), Human heartbeat dynamics Plamennature1999 (), climate change AshkenazyGRL03 () and many financial quantities such as stock prices KwapienPhysRep12 (). This characteristic has been called “multifractality”. A natural question is thus whether such behaviour can be seen in phase transitions. During last two decades people have verified the multifractality of certain quantities on some specific transitions such as random field transitions BenePRA88 () and Anderson transitions at the critical point CastellaniJPhysA86 (); MirlinPRL06 (); RodriguezPRL09 (); RodriguezPRL10 (). However, a general interpretation is still absent on how the multifractality appears with phase transitions. Especially, is the multifractality a universal property of phase transitions as those scaling laws? To address these questions, we argue that one should investigate the property of systems which appears universally in phase transitions. The only ideal quantity to our knowledge is the order parameter.

Another important progress recently in statistical physics is non-extensive Tsallis -statistics TsallisJSP88 (); TsallisPRL95 (), which is a generalization of the conventional Boltzmann-Gibbs statistics. Specifically, when the effective number of degrees of freedom is small, one obtains the Tsallis -statistics following the same arguments to obtain the Boltzmann-Gibbs statistics BarangerPhysicaA02 (); HanelEPL11 (). In this frame some variables such as the entropy and the energy, may be non-extensive when systems have long-range correlations TsallisBJP99 (). It is well-known that the correlation length is infinite at the critical point of a continuous transition. Thus, systems near the critical point are ideal candidates to observe such non-extensive behaviour. Just as the limit distribution of the sum of independent and identically distributed variables is Gaussian or Lvy-stable distribution, in -statistics where the long-range correlation is present we observe a -Gaussian distribution UmarovMJM08 (). The distributions related to Tsallis -statistics have been observed in many experiments, such as cold atoms in dissipative optical lattices DouglasPRL06 (), superdiffusion LiuPRL08 (), solar plasma dynamics Burlaga09 (), spin glass relaxation PickupPRL09 (), tissue radiation response Sotolongo-GrauPRL10 (), and financial signals BorlandPRL02 (); KwapienPhysRep12 ().

Here we show that, at temperatures near the critical point the
distribution of an order parameter related quantity satisfies a
-Gaussian distribution plus a possible Cauchy background. At the
critical point, the distribution shifts between the
Lvy regime and the Gaussian
regime NakaoPhysLettA00 (); DrozdzEPL09 (), triggering the
multifractal behaviour and signalling the phase transition.

Distribution and the measuring quantity

The -Gaussian distribution of interest is a symmetric distribution and has the following form:

(1) |

where is a scale parameter related to the variance of the distribution, indicates the peak position, and is a normalization parameter. Letting one recovers the Gaussian distribution. When , the distribution is in the Gaussian regime. The signal is monofractal and can be described by a single exponent. The distribution falls into the Lvy regime with infinite variance when , and the critical value is TsallisPRL95 (). In this regime decays asymptotically as that of a Lvy-stable distribution. Specifically the Cauchy distribution corresponds to where is a scale parameter. It has been shown that the signal with the -Gaussian distribution in the Lvy regime is multifractal NakaoPhysLettA00 (); KantelhardtPhysicaA02 (); DrozdzEPL09 (). The multifractal behaviour is manifested by two attractors at and , where is the singularity strength and is the singularity spectrum DrozdzEPL09 (). Further, the multifractality may survive even when the time correlation is destroyed by shuffling the data DrozdzEPL09 ().

We hypothesize a universal measuring quantity which possibly satisfies a -Gaussian distribution as follows. It contains only the order parameter and does not depend on details of a specific system, thus should be a dimensionless quantity. We further utilize the “universal” self-similar property of certain physical quantities of the system near the phase transition. This implies that such property does not change when the measuring scale alters. Combining these ideas we construct the desired quantity as the ratio of two order parameters measured at different spatial scales. Its distribution where and are order parameters measured at scale and scale , respectively.

When far away from the critical point the distribution of an order parameter should be Gaussian at different spatial scales. Its mean is zero in a disordered phase and non-zero in an ordered phase. Thus in the disordered phase satisfies a Cauchy distribution. In the ordered phase, we mark the mean and the fluctuations of the order parameter at the scale “” (scale “”) as and [ and ], respectively. For sufficiently large systems and . Thus the ratio is:

(2) | |||||

From equation (2) we find that in the ordered phase
to the first order of is a Gaussian. Therefore, in two
extreme cases, satisfies a -Gaussian distribution. We next
investigate how behaves when the system is near the critical
point from some tunable examples.

Example tunable models

We now focus on some tunable examples of finite size spin systems to verify our hypothesis. For all models we obtain the data through Monte Carlo simulations and apply periodic boundary conditions to systems of different sizes. When the size one obtains the property of a macroscopic system. The system at the scale “1” is the original spin system, and the block spin system at the scale “” is constructed using the standard coarse-graining operation in statistical physics Huang1987 (); Binney92 (); MEFisherRMP1998 () (see the Appendix .1).

We investigate two well-known models: (1) The 2-dimensional nearest-neighbour Ising model with the Hamiltonian , where represents nearest neighbours, ferromagnetic parameter and the spin . This model has a second order phase transition at . The temperature is in units of where is the Boltzmann constant. Its order parameter is the magnetization (per site) and we obtain the time series using the Wolff algorithm WolffPRL89 (). The ratio where and are order parameters measured at the same Monte Carlo step (MCS) “”. (2) The 3-dimensional nearest-neighbour Ising glass model which has a second order phase transition with very slow dynamics near the critical point BhattPRB88 (); ChenPRL10 (). Its and the Hamiltonian where we set satisfying a Gaussian distribution with zero mean and unit variance, and the spin . Each configuration of is one sample. At each temperature after the equilibrium we simulate 150 samples with at least two million MCS per sample. The order parameter is the spin overlap where is a replica of , is the number of spins. The order parameter time series are obtained using the Metropolis algorithm Landau2005 ().

We also provide another model as a representing scenario for general cases of phase transitions. The example is the 2-dimensional nearest-neighbour -state Potts model where is a positive integer and each spin takes values . This model has a second order transition for and a first order transition for WuRevModPhys82 (). Specifically, when it returns to the Ising model. Its , the Hamiltonian and we set . The order parameter is where is the index of a spin, is the value that the spin can take and indicates a thermal average WuJStatPhys88 (). The order parameter series are obtained using the Metropolis algorithm Landau2005 (). Due to the spin symmetry, we obtain utilizing all where is the index of the block spin, is the index of any original spin within the block spin . We have fixed the value of in the calculation.

In the simulations we set the rescaling factor for the Ising model
and the Potts model, and for the glass model.

near the critical point

Results of for the Ising model near the critical point are shown in Fig. 1a and in the Appendix .2. When approaching the critical point from the ordered phase (), a Gaussian fit to gradually does not work well, and it is worse when . In these situations the testing Gaussian fit underestimates at the peak, overestimates on two shoulders, and diminishes faster than the exponential decay compared to the decay of which contains power-law tails. Instead, a -Gaussian fit coincides very well with the central part of at all temperatures. We observe an increasing value of the parameter with increasing temperature. At the -Gaussian enters Lvy regime with infinite variance, and we take this point as our critical temperature for the finite system with the size . Thus from this specific example, we find that near the critical point the order parameter ratio is our desired quantity.

When the ratio is very far away from the peak , we find that separates from such -Gaussian fit. At all temperatures we observe that fat tails of decay as , suggesting possible “universal” origin. Since in the disordered phase follows a Cauchy distribution, we fit these fat tails with this distribution and take them as background. The Cauchy fits work well even when (see the Appendix .2). As shown in Fig. 1a, the summation of a -Gaussian distribution and a suitable Cauchy background at each temperature of interest provides good fit to in all ranges of , i.e., where is due to the Cauchy tails and is the probability of the Cauchy contribution. This implies that the central part and tails part of may be independent. By reversing our procedure we can deduce the origin of the Cauchy part (see the Appendix .3). We find that the distribution of or which contributes to the Cauchy part outside the cross points with the central part of achieves a local maximum and is symmetric about zero, indicating a disorder-like behaviour. It decays with approximately exponential tails for large values of and decays faster for the larger system.

Similar considerations can be done on both the Ising glass model and
the Potts model. We find that the ratio distributions of them both
share the similar behaviour with that of the Ising model. The Gaussian
distribution could not fit the central part of at temperatures
near the critical temperature . This is true even for the
5-state Potts model which has a first order phase transition. In
Fig. 1b we show the fits at for both
models. When approaching we observe Cauchy distributed fat
tails where the ratio is far away from the peak . Further,
combining a -Gaussian distribution and a suitable Cauchy background
provides a good fit to at all ranges of . Nevertheless, for
the Potts model, one has to choose a suitable thermal averaging
interval (in units of MCS) for the order parameter. For
the 5-state Potts model with the size we take
at all temperatures.

General scenario

As shown in Fig. 2a, taken with randomly chosen , e.g., for the 5-state Potts system with the size , may be asymmetric thus would not follow a -Gaussian distribution. With the Potts model as an example, here we propose an approach which can be applied to general physical systems.

First, from the construction of our measuring quantity we hypothesize that the possible “universal” -Gaussian distribution is relevant to the self-similarity of certain physical quantities near the transition. To verify this, we investigate the self-similarity of the distribution of the order parameter at different spatial scales and find that it is broken for small values of (see the Appendix .1). Correspondingly we obtain asymmetric . The self-similarity of becomes better for larger values of . When for the 5-state Potts model with the size the distributions at two scales are almost identical after the rescaling. Correspondingly becomes symmetric.

We thus next quantify how responds when the value of alters. To do this, we calculate the skewness of the central part of . As examples shown in Fig. 2a, the considered region is within two dashed lines. We fix the temperature and study how the skewness changes with different values of (see the Appendix .4 for relevant details). For different sizes of systems we find that the skewness is closer to zero for larger value of . For sufficiently large the skewness stays at around zero, while is symmetric and its central part follows a -Gaussian distribution. The value of seems larger for larger size of systems and it is comparable with the characteristic time scale of the order parameter (see the Appendix .5). Specifically, for the 2-state Potts model is around the number of spins (see the Appendix .4), thus the results of the 2-state Potts model is consistent with that of the Ising model. Starting from , some specific short-range properties of the Potts model are lost and the value from a -Gaussian fit to keeps constant in a broad range of (see Fig. 2b). Thus we consider the value of obtained with chosen in such range as the suitable value of . When further increases, however the value of will increase. In this situation the order parameter itself is gradually not a good quantity to quantify the system. The above are remarks how we should choose the order parameters.

Finally we evaluate the effect of the Cauchy noise which appears for all models. It is pronounced where is far away from and its contribution is . We thus consider as a reliable quantity to measure the strength of the Cauchy noise. In the suitable range of , the of the Cauchy noise is invariant (see Fig. 2b and the Appendix .4). We then study how it depends on the system size for fixed and suitable . As examples shown in Fig. 3, for both the 2-state Potts model from our general approach and the Ising model we find that it decays in power-law with the system size. This implies that for a macroscopic system the Cauchy contribution may go to zero and is -Gaussian in all ranges of . We also study how the of the Cauchy noise depends on the parameter for fixed and . By comparing the results of two equivalent models mentioned above, we find that from the former contains weaker Cauchy noise. Further, we also obtain better -Gaussian fit to of the Potts model which is constructed with fewer points (see Fig. 1).

Critical behaviour and multifractality

We now subtract the Cauchy noise and investigate the rest of , whose behaviour is related to macroscopic systems. In Fig. 4 we show for all three models how the value of from a -Gaussian fit depends on the reduced temperature , where is obtained from the conventional methods. For each finite system with the size , we search for its corresponding critical temperature . We employ this on different sizes of systems and all models. From the finite size scaling we have . We take the values of as those from conventional methods and find good power-law behaviours between and (see inset of Fig. 4). The critical exponents for the Ising model, the Ising glass model and the 5-state Potts model are 1.16, 1.22 and 1.61, respectively. Thus the critical temperatures from our method is consistent with those from conventional considerations.

We next quantify the multifractality in signals using the multifractal detrended fluctuation analysis (MFDFA) PengPRE94 (); ChenPRE02 (); KantelhardtPhysicaA02 (); ChenPRE05 () (see the Appendix .6), which is considered to be a reliable method to quantify the multifractal characteristic in a nonstationary series KwapienPhysRep12 (). With this method we can calculate the -th order fluctuation function , where is the size of the observing window. Varying indicates the multifractality. We also calculate widely-used the singularity strength and the singularity spectrum . When spreads from one point to a variety of points in the X-Y plane, the multifractal behaviour appears.

In Fig. 5 with the Ising model as an example we show the multifractal behaviour of the shuffled series of the order parameter ratio where the time correlation is absent. Similar behaviour is also seen on other models (see the Appendix .7). We find that the order parameter ratio is monofractal when the corresponding is far below . When we approach from the ordered phase, the slope of gradually shifts from small scales of for certain values of . Yet at large scales all fluctuation functions still maintain the identical slope. At the multifractality appears at all scales we observe. Nevertheless, the range of where the multifractality presents shrinks when continuing increasing the value of . We further show the singularity spectrum v.s. . We find that the spreading range of the spectrum is small for small , indicating a Gaussian-like monofractal behaviour. However when approaches , rapidly increases and then slowly expands with increasing .

The multifractal behaviour shown above should not be system specific
since we have considered the shuffled series. To further verify this,
we generate some -Gaussian distributed artificial signals and set
the same parameters for the artificial signals and the testing model
(see the Appendix .8). We find that the results are
almost identical for the artificial signals and the model. Thus such
multifractal behaviour is determined by the two parameters and
of the -Gaussian distribution. To observe the
multifractality, it requires . When , the range where
the multifractality presents is controlled by . For fixed
value of , the larger size of system, the smaller value of
(see the Appendix .8). Correspondingly the larger system
presents multifractal behaviour in broader range of .

Conclusion

With tunable examples we have shown that physical systems
near the phase transition present much more complex self-similar
behaviour represented by the multifractality of the ratio of two order
parameters measured at different spatial scales. At all temperatures
around the critical point the distribution of the ratio follows a
non-extensive -Gaussian distribution plus a possible Cauchy
background which decays in power-law with the system size and may
disappear for a macroscopic system. The -Gaussian distribution
enters the Lvy regime at the critical point, and
triggers the multifractality at all scales we observe. We have
proposed a general approach which relates only to the broken symmetry
yielding zero and non-zero order parameters in different phases as
well as the self-similar characteristics of systems near the
transition. Thus it should be applicable to other physical
systems. The multifractality appears for the order parameter
obtained within suitable range of (time) scales where certain
short-range properties of the specific system are lost and shows
good spatial self-similarity. In this situation the long-range
correlation in the order parameter prevails and the ratio of order
parameters follows a non-extensive q-Gaussian distribution. Our
results suggest that the Tsallis q-statistics may play an important
role in phase transitions.

This work is supported by the National Natural Science Foundation of China (Grant no. 11275184), Anhui Provincial Natural Science Foundation (Grant no. 1208085MA03), and the Fundamental Research Funds for the Central Universities of China (Grants no. WK2030040012 and 2340000034).

## References

- (1) Plischke, M. & Bergersen, B. Equilibrium Statistical Physics, 3rd Edition (World Scientific, New Jersey, 2006).
- (2) Huang, K. Statistical Mechanics, 2nd Edition (John Wiley & Sons, New York, 1987).
- (3) Binney, J.J., Dowrick, N.J., Fisher, A.J. & Newman, M.E.J. The Theory of Critical Phenomena: An Introduction to the Renormalization Group (Oxford Univ. Press, New York, 1992).
- (4) Fisher, M.E. Renormalization group theory: Its basis and formulation in statistical physics. Rev. Mod. Phys. 70, 653-681 (1998).
- (5) Kwapien, J. & Drozdz, S. Physical approach to complex systems. Phys. Rep. 515, 115-226 (2012).
- (6) Amitrano, C., Coniglio, A., Meakin, P. & Zanetti, M. Multiscaling in diffusion-limited aggregation. Phys. Rev. B 44, 4974-4977 (1991).
- (7) Jensen, M.H., Kadanoff, L.P., Libchaber, A., Procaccia, I. & Stavans, J. Global universality at the onset of chaos: Results of a forced Rayleigh-BÃ©nard experiment. Phys. Rev. Lett. 55, 2798-2801 (1985).
- (8) Muzy, J.F., Bacry, E. & Arneodo, A. Wavelets and multifractal formalism for singular signals: Application to turbulence data. Phys. Rev. Lett. 67, 3515-3518 (1991).
- (9) Ivanov, P.Ch., et al. Multifractality in human heartbeat dynamics. Nature 399, 461-465 (1999).
- (10) Ashkenazy, Y., Baker, D.R., Gildor, H. & Havlin, S. Nonlinearity and multifractality of climate change in the past 420,000 years. Geophys. Res. Lett. 30, 2146-2149 (2003).
- (11) Bene, J. & Szepfalusy, P. Multifractal properties in the one-dimensional random-field Ising model. Phys. Rev. A 37, 1703-1707 (1988).
- (12) Castellani, C. & Peliti, L. Multifractal wavefunction at the localisation threshold. J. Phys. A: Math. Gen. 19, L429-L432 (1986).
- (13) Mirlin, A.D., Fyodorov, Y.V., Mildenberger, A. & Evers, F. Exact relations between multifractal exponents at the Anderson transition. Phys. Rev. Lett. 97, 046803 (2006).
- (14) Rodriguez, A., Vasquez, L.J. & Romer, R.A. Multifractal analysis with the probability density function at the three-dimensional Anderson transition. Phys. Rev. Lett. 102, 106406 (2009).
- (15) Rodriguez, A., Vasquez, L.J., Slevin, K. & Romer, R.A. Critical parameters from a generalized multifractal analysis at the Anderson transition. Phys. Rev. Lett. 105, 046403 (2010).
- (16) Tsallis, C. Possible generalization of BoltzmannâGibbs statistics. J. Stat. Phys. 52, 479-487 (1988).
- (17) Tsallis, C., Levy, S.V.F., Souza, A.M.C. & Maynard, R. Statistical-mechanical foundation of the ubiquity of Lvy distributions in nature. Phys. Rev. Lett. 75, 3589-3593 (1995).
- (18) Baranger, M. Why Tsallis statistics? Physica A 305, 27-31 (2002).
- (19) Hanel, R. & Thurner, S. When do generalized entropies apply? How phase space volume determines entropy. Europhys. Lett. 96, 50003 (2011).
- (20) Tsallis, C. Nonextensive statistics: Theoretical, experimental and computational evidences and connections. Brazilian Journal of Physics 29, 1-35 (1999).
- (21) Umarov, S., Tsallis, C. & Steinberg, S. On a -central limit theorem consistent with nonextensive statistical mechanics. Milan J. Math. 76, 307-328 (2008).
- (22) Douglas, P., Bergamini, S. & Renzoni, F. Tunable Tsallis distributions in dissipative optical lattices. Phys. Rev. Lett. 96, 110601 (2006).
- (23) Liu, B. & Goree, J. Superdiffusion and non-Gaussian statistics in a driven-dissipative 2D dusty plasma. Phys. Rev. Lett. 100, 055003 (2008).
- (24) Burlaga, L.F. & Ness, N.F. Compressible ”turbulence” observed in the heliosheath by Voyager 2. ApJ 703, 311-324 (2009).
- (25) Pickup, R.M., Cywinski, R., Pappas, C., Farago, B. & Fouquet, P. Generalized spin-glass relaxation. Phys. Rev. Lett. 102, 097202 (2009).
- (26) Sotolongo-Grau, O., Rodriguez-Perez, D., Antoranz, J.C. & Sotolongo-Costa, O. Tissue radiation response with maximum Tsallis entropy. Phys. Rev. Lett. 105, 158105 (2010).
- (27) Borland, L. Option pricing formulas based on a non-Gaussian stock price model. Phys. Rev. Lett. 89, 098701 (2002).
- (28) Nakao, H. Multi-scaling properties of truncated LÃ©vy flights. Phys. Lett. A 266, 282-289 (2000).
- (29) Drozdz, S., Kwapien, J., Oswiecimka, P. & Rak, R. Quantitative features of multifractal subtleties in time series. Europhys. Lett.88, 60003 (2009).
- (30) Kantelhardt, J.W., et al. Multifractal detrended fluctuation analysis of nonstationary time series. Physica A 316, 87-114 (2002).
- (31) Wolff, U. Collective Monte Carlo updating for spin systems. Phys. Rev. Lett. 62, 361-364 (1989).
- (32) Bhatt, R.N. & Young, A.P. Numerical studies of Ising spin glasses in two, three, and four dimensions. Phys. Rev. B 37, 5606-5614 (1988).
- (33) Chen, Z. & Yu, C.C. Comparison of Ising spin glass noise to flux and inductance noise in SQUIDs. Phys. Rev. Lett. 104, 247204 (2010).
- (34) Landau, D.P. & Binder, K. A Guide to Monte Carlo Simulations in Statistical Physics, 2nd Edition (Cambridge University Press, Cambridge, UK, 2005).
- (35) Wu, F.Y. The Potts model. Rev. Mod. Phys. 54, 235-268 (1982).
- (36) Wu, F.Y. Potts model and graph theory. J. Stat. Phys. 52, 99-112 (1988).
- (37) Peng, C.-K., et al. Mosaic organization of DNA nucleotides. Phys. Rev. E 49, 1685-1689 (1994).
- (38) Chen, Z., Ivanov, P.Ch., Hu, K. & Stanley, H.E. Effect of nonstationarities on detrended fluctuation analysis. Phys. Rev. E 65, 041107 (2002).
- (39) Chen, Z., et al. Effect of nonlinear filters on detrended fluctuation analysis. Phys. Rev. E 71, 011104 (2005).

Appendix

## .1 Self-similarity in distributions of the order parameter

When a system is in the vicinity of the critical point, many physical quantities of it display self-similar (scaling) properties. The distribution of the order parameter also shares this characteristic. In Fig. S1 we show an example for the 2-dimensional 5-state Potts model with the size . We construct the associated block spin system with size utilizing the standard coarse-graining procedure. For a -dimensional spin system with lattices, we can transform it into a block spin system with block spins. Each of block spins contains spins. A coarse-graining operation is reached when we substitute each of block spins by a single spin with the value determined by the majority rule. For example, for the Ising model the block spin is +1 if there are more spins up than down, and vice versa. In particular, when the amount of spins up and spins down is equal, we assign the block spin +1 or -1 randomly. In this way we have constructed the Ising system at the scale “”. The block spin system for the Potts model can be constructed similarly.

Nevertheless, for the Potts model we do not see good scaling when the thermal averaging interval for the order parameter is too small (see Fig. S1). In this situation contains certain short-range information specifically related to the Potts model. It is well-known that the critical properties of a physical system do not depend on the short-range details, but on the characteristics of long-range fluctuations. Such short-range information is not self-similar and may diminish with larger where the time correlation in the order parameter becomes weaker. (For example, see the Appendix .5.) For sufficiently large value of , the distribution of the order parameter for the original system and that for the block spin system become similar, i.e., . Interestingly, when such scaling relation is effective, the distribution of the order parameter ratio becomes symmetric and follows a -Gaussian distribution, as we show in the paper.

## .2 Cauchy noise in distributions of the order parameter ratio

In Fig. S2 we present the Cauchy distributed fat tails in distributions of the order parameter ratio . For the Ising model we observe very good fit to the tails of down to , while the length of the order parameter ratio is 8-16 million. We also note that at high temperatures with close to 2, the Cauchy background may be still pronounced. In an example shown in Fig. S2 for the system at and with the size , we find that a single -Gaussian fit works well for the central part, however it still underestimates at positions of far away from the peak .

The above behaviour has also been found in the Ising glass model and the Potts model. For the former at each temperature we have 150 samples each with 2 million points after the equilibrium. For the latter at each temperature the length of the order parameter ratio is 5.8 million. As examples in Fig. S2 we show such behaviour at the critical point which we have defined in the paper. For both models good fits to the tails of down to are observed.

## .3 Origin of the Cauchy noise

By reversing our procedure we can deduce the origin of the Cauchy distributed fat tails. To see this, we define a quantity , where is our Cauchy fit and is the distribution of the model from the Monte Carlo simulation. At each Monte Carlo step with the order parameter ratio , we take as the probability of this data point belonging to the Cauchy noise. We did this at each and thus could obtain the corresponding and which belong to the Cauchy noise. Such considerations can be done in all ranges of or in partial range of , e.g., we can choose the range of which is outside the cross points of the central part and the Cauchy part of . As shown in the orange and brown dashed lines of Fig. S3 we find that the distribution of or which is associated with the Cauchy part of outside the cross points with the central part of achieves a local maximum and is symmetric about zero, then decays with approximately exponential tails for large values of , indicating a disorder-like behaviour. Further, the tails decay faster for the system with larger size. This implies that the fat tails in may go to zero when the system size . Such behaviour is also true when the distribution of the order parameter is asymmetric.

## .4 Asymmetric and symmetric with different

For the -state Potts model with we find that the distribution of the order parameter ratio is not symmetric when the thermal averaging interval is small, as shown in Fig. 2a of the paper. Further, such behaviour diminishes with larger values of . To characterize this, we can calculate the skewness of around the peak . (In practice, the width of the region we choose to calculate is around , where is the scale parameter of the testing -Gaussian fit to .) As shown in Fig. S4a, such asymmetry seems larger at lower temperature with smaller value of the parameter . Further, it is also stronger for the Potts system with more possible spin states. This is manifested in Fig. S4b, where we fix the -Gaussian fit parameter . We find that — the minimal value of the thermal averaging interval to obtain symmetric — increases with increasing value of the state number . Specifically, we note that is consistent with the results of the Ising model which is related to the 2-state case of the Potts model. Since for the system has a second order phase transition, such behaviour is not related to the order of the transition.

For sufficiently large thermal averaging interval , becomes symmetric, and it would keep symmetric when continuing increasing the value of , as shown in Fig. 2b of the paper. Further, the value of the parameter in a -Gaussian fit to also keeps almost invariant in a broad range of . Here we show some examples in Fig. S5. In this suitable range the strength of the Cauchy noise, marked by the of the noise, also keeps constant.

## .5 Time correlation in the order parameter time series

The order parameters which we obtain through Monte Carlo simulations may be correlated in time. Such correlation may be stronger when the system is in the vicinity of the critical point. As an example, here we investigate the 5-state Potts model near the phase transition. We fix the temperature and vary the size of systems. To obtain the time correlation function , for different sizes of systems we fix the thermal averaging interval . The order parameter of the Potts model depends on the spin value and the spin index . We also fix and for each original spin we calculate the time correlation . We show the average of all original spins in Fig. S6. When the time decay is small decays exponentially, i.e., where is the characteristic time scale. Comparing with the results we show in Fig 2b, we find that the distribution of the order parameter ratio becomes symmetric when the thermal averaging interval is comparable or larger than the characteristic time scale of the order parameter .

## .6 Multifractal detrended fluctuation analysis (MFDFA)

The multifractal detrended fluctuation analysis (MFDFA) is an efficient and reliable method available to quantify the multifractality in a non-stationary series. For a series with the length , its procedure is the following:

(1) Calculate the profile

where is the mean of .

(2) Divide into non-overlapping parts of length . The maximal value of should be smaller than to avoid statistically unreliable results. We perform this step from both the beginning and the end of the signal. Thus we totally obtain parts.

(3) For each part with index we calculate the local trend in this part by a least-square fit of the series with a polynomial function, where is the order of the polynomial function. We then subtract this local trend and determine the variance:

for and

for

(4) Average over all parts we obtain the -th order fluctuation function:

For we choose

(5) Determine the scaling of the fluctuation function: . If is a constant, the signal is monofractal; otherwise it is multifractal.

The singularity spectrum can be further calculated from: where the singularity strength . If the signal is monofractal, we find that and .

## .7 Multifractal behaviour of the 5-state Potts model

As another example here we investigate the multifractal behaviour of the 5-state Potts model in the vicinity of the critical point. We set the thermal averaging interval at all temperatures. After obtaining the series of the order parameter ratio, we first shuffle the data and then apply the MFDFA. The results of the fluctuation function versus the observing window are shown in Fig. S7. We find that, below the critical point (in the ordered phase) the system is monofractal since the slope of is identical for different values of . When approaching the critical point from the ordered phase, the slope of starts to shift at small scales, indicating a multifractal behaviour within these scales. Such region is broader when the system is closer to the critical point. At the critical point we observe multifractal behaviour at all scales we measure. When continuing increasing the temperature, the system is in the disordered phase. The value of the parameter from the -Gaussian fit is increasing while the region where the multifractal behaviour presents shrinks. Such behaviour is similar to that of the Ising model. Nevertheless, for the Potts model we find that the scale parameter of the -Gaussian fit is increasing with increasing temperature and fixed , while it is not sensitive to the temperature for the Ising model if the measured is far below 2.

## .8 Characteristics of the multifractal behaviour

Here we investigate the characteristics of the multifractal behaviour we show in the paper. To do this, we generate some -Gaussian distributed artificial signals and we set the same parameters for the artificial signals and for the testing model — the Ising model. For the model we have first shuffled the data thus eliminated the dynamics in signals. As shown in Fig. S8 we find that the results are almost identical for the artificial signals and the model when the parameters and are the same. This verifies that the multifractal behaviour we present in the paper is not a system specific behaviour but should be applicable to other physical systems.

We then investigate how the multifractal behaviour depends on the size of a specific system. Such dependence should be manifested from the values of two parameters and . To see this, we investigate the Ising model with different sizes and . We first shuffle the data and eliminate the dynamics in signals. We then apply the MFDFA and find that with the same parameter , the difference in the observed multifractal behaviour is owing to different values of the scale parameter (see Fig. S9). When is fixed, the larger system with smaller value of shows multifractal behaviour in broader region than the smaller system shows.