Threshold-induced correlations in the Random Field Ising Model

We present a numerical study of the correlations in the occurrence times of consecutive crackling noise events in the nonequilibrium zero-temperature Random Field Ising model in three dimensions. The critical behavior of the system is portrayed by the intermittent bursts of activity known as avalanches with scale-invariant properties which are power-law distributed. Our findings, based on the scaling analysis and collapse of data collected in extensive simulations show that the observed correlations emerge upon applying a finite threshold to the pertaining signals when defining events of interest. Such events are called subavalanches and are obtained by separation of original avalanches in the thresholding process. The correlations are evidenced by power law distributed waiting times and are present in the system even when the original avalanche triggerings are described by a random uncorrelated process.

When a threshold is imposed on a signal, the events of interest are connected bursts of activity above the threshold. Each such burst is a subavalanche of some underlying avalanche, comprising entire activity of the system at the current stage of its evolution. A subavalanche, selected by thresholding, begins at the moment of time when the signal exceeds the threshold, and ends when the signal falls below it. The difference between these two moments is taken as duration T, and the area between the imposed threshold and the portion of the signal above it as the size of subavalanche, T s 0 th , see Fig. 1a, where t is the time measured from the start t s of the subavalanche.
Once a threshold is applied, a portion of the signal will remain below it. This implies an introduction of the concept of waiting time T w , describing the time interval between two consecutive subavalanches, selected by the imposed threshold. Provided the start and end of each avalanche are known, like in simulations or in experiments after some minimal threshold is imposed, one can differentiate two kinds of waiting times. The internal waiting time is the waiting time T int between two consecutive subavalanches thresholded from the same avalanche (avalanche j in Fig. 1a), while the external waiting time T ext is the time between two consecutive subavalanches thresholded from two different avalanches (avalanches i and j in Fig. 1a).
A further distinction can be made among different types of contribution to the external waiting time T ext (i, j; V th ); thus, see Fig. 1a, where T end (i; V th ) is the time taken by the avalanche i to end (i.e. fall from V th to zero), T mid (i, j; V th ) is the time spent by a whole sequence of consecutive avalanches that lie between the avalanches i and j and remain below V th (note that this sequence may be empty), and T ini (j; V th ) is the time taken by the avalanche j to rise from zero to V th . In the main panel of Fig. 1b  Thresholding of RFIM signal. The signal is obtained in simulations of 3D system with size L = 1024 and 40 random field configurations for each disorder R. (a) For a (blue) part of a train of avalanches (shown in bottom, and zoomed in top panel) and the imposed threshold V th (red line), we illustrate: the determination of size S and duration T of a subavalanche (starting at the moment t s and ending at t e = t s + T) taken out of the avalanche i, the internal waiting time T int between two subavalanches of avalanche j, and the contributions T end , T mid , and T ini to the external waiting time T ext between avalanches i and j, see Eq. (1). (b) Distributions D(T) of duration (main panel), and distributions D(S) of size (inset) of subavalanches selected by thresholds from a wide range shown in legend. (c) 〈S〉 T shown against T for the thresholds in legend, where 〈S〉 T is the average size of subavalanches with duration T; variation of exponent γ S/T with V th is shown in inset. (d) γ S/T vs V th data, obtained for various disorders R (see legend), collapse onto a same curve when presented against V th r, where the reduced disorder r = (R − R c )/R measures a distance to the critical disorder R c of the model. Inset shows how γ S T / (0) (i.e. the exponent γ S/T taken for V th = 0) depends on the reduced disorder r. c eff for the given system size (see Methods for more details). The distributions are collected for the subavalanches extracted above threshold V th from the avalanches triggered in a zero-centered window of external magnetic field in which the response signal can be considered as stochastically stationary. We found that in a wide range of thresholds, both types of distribution follow power-laws , terminated by the cutoff scaling functions g T (x) and g S (x) for duration and size, respectively. The cutoff time T 0 , and the cutoff size S 0 , decrease when the threshold V th increases. The values of exponents, pertaining to these distributions are: τ T = 1.64 ± 0.02 for duration, and τ s = 1.38 ± 0.03 for size of subavalanches. These values are obtained using 0 0 S , where σ T and σ S are the cutoff exponents whose values are close to 1 for all the distributions being analysed.
The average size 〈S〉 T of subavalanches with duration T is shown against T in the main panel of Fig. 1c for a family of curves, corresponding to different values of threshold V th . The graph demonstrates that 〈 〉 ∼ γ S T T S T / in a broad range of thresholds, and that the power-law exponent γ S/T varies with threshold. As can be seen in inset, when the threshold increases, the exponent γ S/T decreases (from the value 1.77 for very low thresholds, to the value of 1.44 for very high threshold levels), forming some sort of plateau, like in the crack-line propagation model 30 .
In order to gain a more complete insight into the variation of γ S/T with V th , we show in the main panel of Fig. 1d its values against V th r, i.e. threshold multiplied by r, where r is the reduced disorder r ≡ (R − R c )/R, measuring the distance to the critical disorder R c of the model. The data obtained for different disorders collapse onto a same curve, suggesting that a joint plateau is formed at the value γ = . ± .
1 49 0 02 . The existence of plateau may be considered as an important feature of the model, because the plateau value γ S T / (pl) remains stable under variation of both threshold and disorder unlike, for instance, the value γ S T / (0) of exponent γ S/T for zero threshold, which (seeming linearly) changes with disorder -see the inset of for the exponent of the waiting time, and terminated by the cutoff scaling function g w (x), taken for T w /T w,0 , where T w,0 is the cutoff waiting time. In contrast to the cutoff time T 0 which decreases with threshold, the cutoff waiting time T w,0 increases with V th . This is shown in the bottom inset of Fig. 2a, where one can see that the relation t h is satisfied with δ ≈ 1.30 ± 0.02. Regarding the shape of distributions D(T w ), one can see that the power law part appears with the increase of threshold, indicating the onset of correlations due to avalanches that are partially hidden below the detection threshold. Opposite to that, one can notice that the power law gradually vanishes for the very low threshold levels due to very small value of the cutoff waiting time T w,0 . This is illustrated in the top inset of Fig. 2a, where the power law reduces to the (approximately) exponential cutoff scaling function g w (T w /T w,0 ). The distribution of waiting times D(T w ) for the zero threshold, should be given by a delta function limit of the exponential distribution with vanishing cutoff waiting time T w,0 , due to specific pattern of driving explained in Methods, giving no time separation between consequtive avalanches.
Finally, in Fig. 2 panel b, we present the cutoff time T 0 and the cutoff waiting time T w,0 against threshold V th for a family of curves obtained for different disorders R. For small thresholds, the cutoff waiting time grows as a power law ( , and more rapidly than that for larger V th . This leads to the conclusion that if the threshold is so high that only the avalanches from the cutoff of the "true" distribution are "observable", then a very large separation between the time scales of the waiting times and avalanche durations may ensue 32 . Pertaining scaling collapse of the cutoff waiting times, corresponding to different disorders, is obtained by multiplying the threshold axis by r, as is shown in the inset.

Scaling theory of avalanche correlations. Scaling properties of the threshold induced correlations can
be predicted for signals emitted by a wide range of systems that respond to stationary external driving in stochastically stationary trains of avalanches. To this end, let V(t) be a response signal that continuously varies with (continuous) time t, and let V(t) > 0 at any moment of time when the system is active, while V(t) = 0 otherwise. Such signal is a sequence of avalanches, separated by intervals of time when the system is quiet, each avalanche being a continuous burst of values V(t) > 0 taken at any t s < t < t e between the moments t s and t e when the avalanche starts/ends, and therefore V(t s ) = V(t e ) = 0. Essentially the same can be said for discrete signals sampled at discrete moments of time, provided that they are continuously (e.g. linearly) extrapolated to all moments of time throughout the signal duration, like it is done here with the RFIM signals.
The avalanches can be classified into types so that all avalanches of the same type have the same profile f(t′) ≡ V(t s + t′) with respect to time t′ measured from their start. Using this equivalence relation, one can obtain the set of all avalanche types I, and introduce a one-parameter family of scaling transformations where x is an exponent specified by the type of the involved system 5 . For this family of transformations one can prove that where T i is duration, and , and also for the following types of waiting times: all for the subavalanches extracted above the threshold V th out of avalanche of type i, and for the corresponding x end t h describing the scaled avalanche ˆi S b , extracted above the scaled threshold b V x th . Furthermore, if one scales the whole portion of response signal, starting with an avalanche (of type i) that surpasses threshold V th , and ending with successive avalanche (of type j) above the same threshold, then analogous expression hold for the waiting time T mid (i, j; V th ), spent by the avalanches that lie between these avalanches and remain below V th , namely: which, combined with Eqs (1) and (5), gives the same type of scaling  Next, let dp(i; λ) be an elementary probability of obtaining an avalanche of type i in the response signal under observation conditions specified by some appropriate multiparameter λ = (λ 1 , λ 2 , …λ n ), like λ = (h′, r, 1/L) for the RFIM signals (L is system size, r is the reduced disorder, and h′ is the reduced magnetic field, see Methods). Having at our disposal this probability, one can express the distribution D T (T; V th , λ) of subavalanches that are obtained under conditions λ and have the duration T above the threshold V th . Thus, where w is a probability exponent, and ζ = (ζ 1 , are the scaled conditions under which the type ˆi S b is observed. Starting from this expression, which is in fact a generalized scaling hypothesis, one can obtain the scaling laws The foregoing general predictions can be tested in the case of any response signal for which can be expected that the assumptions, used in their derivation, are reasonably satisfied. The first step towards that in the case of RFIM signal is to specify the observation conditions λ, and next to express the generic exponents x, y, w, and ζ in the terms of standard RFIM exponents. Here, as we already mentioned, the observation conditions are λ = (h′, r, 1/L), while for the exponents x and w, and for the multiexponent ζ = (ζ h , ζ r , ζ L ), we found that h r L where σ, ν, z, α, β and δ are the standard RFIM exponents 3,9,10 . In Fig. 3 we present the collapsing for distributions of duration and of various types of waiting times, all for the subavalanches above thresholds. The subavalanches are taken from a family of response signals observed under conditions which are aligned according to the collapsing requirements together with the corresponding collapsing predictions. Thus, in panel a, the data are scaled in agreement with  power-law parts of waiting time distributions gradually vanish as the value of the threshold decreases, down to the lowest threshold when the distribution of waiting times turns approximately exponential. This implies that even though the avalanches are triggered by a random process, applying a finite threshold implicitly introduces underlying temporal correlations in a given signal. On the other hand, distributions of the subavalanche durations and sizes follow power law with cutoffs decreasing when the threshold level increases. Thus the present paper verifies the fact that the scaling predictions obtained for the previously studied crack-line propagation model 30,32 also hold in the case of RFIM signals, suggesting their general validity in accordance with the scaling hypothesis and derived scaling forms proposed in the previous section. In our analysis we have also identified different contributions of waiting times (T int , T ext , T ini , T end , T mid ), which all follow the same scaling form Eq. (20). This form predicts that rescaled distributions of each type of waiting time, obtained for system parameters adjusted according to the collapsing requirements, all collapse onto a single curve. As the proposed general scaling theory predicts, we anticipate that the derived scaling forms of the threshold induced correlations, provided that the requirements of the theory are fulfilled, should hold for any response signal originating from system that, as a response to slowly changing external conditions, relaxes in an avalanche-like intermittent way. For such response signals, one can say that, in general, the external waiting times are affected by the two mechanisms: (i) implicit or explicit application of a finite threshold resulting in low levels of avalanche activity not being detected, and (ii) effects due to details of the implementation of the external driving. While we focus here on the first mechanism, a more detailed study of the interplay between the two, considering the joint effects of finite thresholds and driving rates, would be an interesting avenue of future work.
Additionally, we have found that the exponent γ S/T , obtained from the scaling of the average avalanche size with duration, is affected by the level of the applied theshold, with theoretically expected value of 1.77 recovered only in the limit of very low thresholds. The effective value of γ S/T deviates from this value in a way that it initially decreases and then reaches the plateau where it remains stable for a wide range of threshold and disorder values. This result is of special importance for analysis of experimental results where setting a finite detection threshold is an inevitable procedure in order to be able to perform the analysis of the recorded signal.
Qualitatively speaking, our results meet very well to the ones obtained for the formerly studied crack-line propagation model 30,32 . Given that our numerically generated data are noise-free, we anticipate the possibility to observe in real experiments the additional effects that the presence of noise could impose on the thresholded signal 28 . Thus, our work calls for further and more detailed investigation of experimental data in order to unravel the true nature of underlaying mechanisms causing the onset of temporal correlations in these systems.

Methods
Simulations of the Random Field Ising Model. In order to investigate the threshold induced temporal correlations in the Random Field Ising Model (RFIM) we performed numerical simulations of its athermal (T = 0) variant in nonequilibrium adiabatic regime. The athermal RFIM describes a system of N ferromagnetically coupled classical Ising spins S i = ±1, located at sites i of some underlying lattice. The spins are influenced by a homogeneous external magnetic field H, and by some local magnetic field h, whose values h i vary randomly from site to site due to random distribution of quenched impurities that generate that field. Therefore, it is considered that each instance of RFIM system is specified by the configuration of values 1 that the quenched random field h takes at the sites i of that system, and that these values remain frozen throughout any evolution of the system.
In the basic RFIM case, presented here, the ferromagnetic coupling between spins extends only to the nearest-neighbors, so the effective magnetic field acting on spin S i is i j j i eff where S j are the nearest neighbors of S i , and J > 0 is a ferromagnetic coupling constant. Hence, the system Hamiltonian reads and in this expression the first sum refers only to the nearest-neighbor spins S i and S j , the second term describes the coupling of spins with the external field H, while the last term gives the coupling with quenched random field h. At any site i, its value h i is taken randomly from some zero-centered distribution, same for all sites, and for any two different sites i ≠ j the values h i and h j are chosen independently, so that the expected value of their product is 〈h i h j 〉 = 0. For generating the values of the quenched random field, here we use a Gaussian distribution and take its standard deviation = 〈 〉 R h 2 as a measure of disorder in the system. In the nonequilibrium athermal RFIM, the system evolves according to the following flipping rule: each spin S i remains stable while its sign equals the sign of the effective field h i eff at its site; otherwise, S i becomes unstable, and flips at the next moment t + 1 of discrete time. Thus, the flipping of each spin influences the effective field for all of its nearest neighbors. All those neighbors that become unstable will flip at the next moment of time, which in turn may cause the flipping of their neighbors, making an avalanche which lasts until all spins become stable.
Once all spins become stable, the only way to trigger a new avalanche is to change the external field H, and in this way drive the system by a sequence of H-increments, forming a driving pattern that is set in advance. Typically, the changes between two consecutive moments of time are small, resembling the usual real-world situation with two well separated time scales: fast one for spin flipping, and slow scale for the external field. In the limiting regime of infinitely slow (i.e. adiabatic) driving, the external field is kept constant during (any) avalanche. After the avalanche dies, H continues to change (i.e. increase or decrease following the current direction of driving pattern) until it reaches exactly the value that triggers only the least stable spin. Note, however, that because all spins during the foregoing change remain stable, and therefore unaltered, the overall change of H is allowed to be done in a single jump, which is utilized in computer simulations for better efficiency. The consequence of such driving pattern is that the next avalanche is triggered immediately after the previous one has ended.
Together with the driving pattern, one also needs to specify some initial and final conditions. Here, we take that initially H = −∞ and all spins are −1, and then we gradually increase H until all spins become +1. At each moment of time t, we register the number of spins V(t) flipped at that moment, and in this way collect system's response along the whole rising part of the saturation hysteresis loop. Note that if one repeats the run with the same sample (i.e. same configuration of random field), in the same driving regime, and with the same initial conditions, the system response will be the same because the flipping rule is deterministic. Therefore, reliable SCIEntIfIC RepoRTS | (2018) 8:2571 | DOI:10.1038/s41598-018-20759-6 avalanche statistics are collected by repeating the whole procedure many times using different random field configurations (quenched or sample averaging).
Our RFIM simulations are done in the nonequilibrium adiabatic regime with J = 1, and with closed boundary conditions on 3D lattices L × L × L of linear size L = 1024. For the analysis of the effects of the of imposed threshold V th , we have used the parts of signal where the signal can be considered to be approximately stationary. This is fulfilled in the narrow window of external field H, taken around the coercive value (i.e. the value at which the system magnetization = ∑ = M S i N i 1 is zero). One fragment of such signal is shown in the bottom panel of Fig. 1a. The number of spins flipped at a given moment of time t is taken as the signal value V(t) at that moment, and the time is measured from the window start. The statistics, collected by quenched averaging for given disorder R, are described using the scaling variables: reduced disorder r = (R − R c )/R and reduced magnetic field h = H − H c − b r r, where R c is the critical disorder, H c is the critical value of magnetic field H, and b r is the rotational parameter accounting how the effective critical value H r ( ) c eff of the external field (i.e. the value of H at which the maximum of susceptibility occurs), shifts with reduced disorder r 3,10 . In this paper we have confined our study to disorders that are above the effective critical disorder R c eff pertaining to the underlying lattices 33,34 . For this value of disorder, precisely defined in the quoted references, it is likely for the spanning avalanches to occur. The spanning avalanche is the avalanche that spans the finite system along at least one of its dimensions and therefore plays the role of infinite avalanche, causing the jump of magnetization in infinite systems below the critical disorder. In three-dimensional RFIM systems, these avalanches have different distributions than remaining (i.e. non-spanning) avalanches 35 , and violate the scaling assumptions given in section Scaling theory of avalanche correlations. Data availability. The datasets generated during and analysed during the current study are available from the corresponding author on reasonable request.