Hysteresis-Loop Criticality in Disordered Ferromagnets–A Comprehensive Review of Computational Techniques

Nomenclature

1  Introduction

Disordered ferromagnetic materials are the subject of intense theoretical and experimental research due to their physical characteristics associated with the domain structure [1–3]. Recent research surveys [4–7] elucidate that these materials, especially low dimensional samples and ferromagnetic-antiferromagnetic heterostructures, are considered promising materials for modern technology applications. Disordered ferromagnets possess intricate domain structure leading to the hysteresis behavior [8,9] in the external magnetic field; see also recent review [10] and references therein. Hysteresis behavior is vital for various applications of magnetic materials [11]. Therefore more attention is devoted to predicting and optimizing hysteresis properties using different numerical tools [12]. During the reversal processes by slow ramping of the external field over the hysteresis loop, the moving domain walls interact with the structural and magnetic disorder [13], leading to stochastic changes of the magnetization. The magnetization fluctuations are experimentally measured as Barkhausen noise (BHN), for example, in bulk materials [14–17] and thin films [18–21], as well as in samples with varied thickness [22,23]. Considering the Ising model, one of central pillars of statistical physics [24], theoretical investigations focusing to the dynamic phenomena on hysteresis use the field-driven spin models with random bond [25] or random field disorder [26–28]. These investigations revealed that the collective nature of the magnetization changes with avalanche-like behavior on the hysteresis are associated with a disorder-induced critical point [29–32].

This out-of-equilibrium critical dynamics is characterized by long-range temporal correlations and multifractal features of the BHN [33]. Among different types of models of disordered systems, the Ising model with random fields [34–37] appeared as the most attractive for theoretical investigations of the in- and out-of-equilibrium criticality. It was also recognized as an appropriate model for weakly disordered antiferromagnetic materials in an external field, where structural disorder induces local random fields [38–40]. For more details regarding recent studies on antiferromagnetic systems, see [41–44]. It should be noted that the 100 years of the Ising model, celebrated this year, have shown its relevance to many different phenomena in complex systems [45]. However, certain limitations of the model in describing the magnetism of solid materials are evident, for example, taking into account more complex domain walls in domain structures that are associated with hysteresis phenomena, distinguishing hysteresis properties along different easy axes and the influence of other types of disorder that do not break the local rotational symmetry. In addition to the theoretical interpretation of the behavior of disordered ferromagnetic systems, using the random-field Ising model in interpreting experimental Barkhausen noise proved crucial for some systems [17,46]. In a more general context, statistical physics and random field theory were recently applied for spatial data modelling, as described in the book [47].

The model treats a system of Ising spins located at lattice sites and quenched impurities interpreted by randomly distributed magnetic fields. Among the nearest neighboring spins the ferromagnetic interaction occurs, in addition to the interaction with the external magnetic and a quenched random field. The zero-temperature variant of the model proposed in [48] captures the essential features of avalanching dynamics while neglecting thermal effects. Triggered by the external magnetic field changes, the system relaxes in spin-flipping avalanches, reflected by the magnetization jumps. The stochastic process of magnetization reversal occurs along the hysteresis loop, characterizing the collective response of the system. The occurrence of the out-of-equilibrium critical dynamics on the hysteresis loop together with its dependence on pertinent physical parameters, system’s spatial dimensionality and shape, represent essential challenges for controlling and predicting the reversal processes in these materials.

Recently, it has been recognized that besides hysteresis-loop criticality in disordered ferromagnets, numerous other complex systems exhibit a similar avalanche-like response to the external driving forces. Further examples of such systems span from earthquakes [49–53] and propagation of interfaces in various kinds of random media [54–58] to brain networks and neuronal activity [59–62], dislocations in crystal structures [63–66], invading imbibition fronts in porous media [67,68], response of the mechanically pressured wooden materials [69] to financial booms [70] and epidemics [71], all having in common that the underlying systems evolve through the metastable states due to an avalanche type of relaxation. Despite being different at a microscopic level, these systems have some universal features of their critical states in common. This universality implies the possibility of extending the findings and methodologies developed on the systems that are tractable for experimental realizations, such as disordered ferromagnets, to characterize such phenomena in systems that are elusive to controlled experiments.

Here, we provide a comprehensive review of numerical approaches to model the criticality of the hysteresis loop in disordered ferromagnets based on extended nonequilibrium (NEQ) Random Field Ising Model (RFIM) with different physical parameters and driving regimes inspired by experimentally achievable conditions. The presented overview of the results shows how these modeling approaches, supported by advanced computational techniques using the inherent scale invariance of the dynamics, adequately describe the experimentally observed hysteresis behavior and provide new physical insights into the dynamic critical phenomena of disordered ferromagnetic systems. The developed computational techniques can be adjusted to analyze nonequilibrium critical dynamics of ferromagnetic/antiferromagnetic bilayers [72,73] and various quantum [74–76] and classical complex systems, in particular, nanonetworks [77] and higher-order self-assembled nanostructures [78–82] that are increasingly interesting to modern technology.

Challenges in Numerical Modeling of Disordered Ferromagnets Hysteresis Behaviour

Occurrence of extreme events such as infinite or system-spanning avalanches underpins the system criticality, characterized by the scale invariance of the created events whose features are power-law distributed and diverging correlation length [34]. Over time, several models were developed to describe this complex phenomenon and answer whether the model under study exhibits a critical behavior [83]. Unraveling this is not an easy task due to the multitude of mutually intertwined factors to which the underlying systems are sensitive, such as the degree of disorder or the type and rate of driving, presence of the demagnetizing fields and thermal effects, sample geometry, which is of particular relevance for applications, and other. In the following, we briefly describe these factors that represent challenges for computational modeling of the hysteresis-loop critical phenomena.

As the extensive studies [28] have shown, the critical behavior of the zero-temperature (ZT) RFIM depends on the dimension d of the studied system [31,84]. In the range of dimensions 2≤d<6, ZT NEQ RFIM displays a nontrivial critical behavior [31], while for d≥6, the behavior of the system is described by mean-field approximation [35,84,85]. In the nonequilibrium field-driven version, the system evolves through metastable states that represent local, and not global energy minima, in contrast to the equilibrium critical behaviour studied by the RFIM [86–88].

Criticality of the equilibrium and nonequilibrium version of the model for systems with dimension d≥3, displayed in the matching of many avalanche properties, lead to the presumption that they share the same universality class [37,87,88]. However, that is not the case for the 2D model since the ferromagnetic phase, whose existence is demonstrated in the case of nonequilibrium model [89,90], is not possible in the thermodynamic limit of the equilibrium version [91]. In an attempt to describe this nontrivial critical behavior, both perturbative [92–94] and nonperturbative renormalization group approaches [95,96] were used. Recently some advances have been made on understanding the principles of universality [36,97], dimensional reduction [98] and supersymmetry [99] in the equilibrium version of the model. On the other hand, dimension and geometry influence of underlying lattice are shown to have different impact on the nonequilibrium critical behavior [100–103].

Particular attention in the recent studies of nonequilibrium systems is devoted to the extreme, catastrophic events in which most of the system’s constituents change their state [104,105]. Understanding the particular conditions under which these events occur would reveal essential details regarding the mechanism of some natural phenomena such as earthquakes, snow avalanches, cracks in materials, etc. [27,106–108], that are of immense importance due to the possible significant consequences they might cause.

In recent years, due to the increased interest in the practical applications of thin ferromagnetic materials, studies of nonequilateral systems with different geometry aspects [109] and thin systems [110–114] emerged, with some aspects in resemblance to the criticality of spin systems situated on a complex network topology [115–118]. Disordered ferromagnetics are mainly used as memory materials in the form of thin films [119–122] and nanowires [123–126], putting forward the experimental investigations of critical dynamics of BHN [127–129]. The ZT NEQ version of the RFIM was very suitable for interpreting actual experiments since the thermal fluctuations are negligible in most of them, and the underlying dynamics is closer to the one in externally driven ferromagnets. Driving these systems by the external magnetic field at finite rates [130] provides a more realistic framework for the numerical interpretation of experimental measurements conducted on ferromagnetic strips and ribbons [21,131,132] and has been of increasing interest in recent experimental studies [17,22,114,133].

Besides dimensionality, one may ask to which extent the model describes the influence of a variation of the geometry of the underlying lattice [109], including the reconsideration of the universality classes within the RFIM [100,134]. Recently, a competitive conjecture was put forward that the universality classes in the NEQ model are determined by the topology of the lattice emphasizing the role of its coordination number in addition to its spatial dimensionality [102,134–136]. This conjecture was to some extent confirmed by contrasting the criticality of the model on triangular [100], quadratic [89] and hexagonal (i.e., honeycomb) [103] 2D lattices. The absence of critical behavior for the latter lattice introduced doubt in the presumed existance of a single universality class for all periodic 2D lattices in the case of ZT NEQ RFIM.

Another significant aspect, particularly for experimental studies, is the presence of unwanted external noise in the recorded data requiring the inevitable procedure of using the discrimination threshold. In the experimentally obtained signals, avalanches are identified as signal parts outside the imposed threshold region (i.e., the region between lower and upper threshold levels). Recent experimental [137] and theoretical [138–141] studies delivered some relevant results regarding the consequences of the implementation of the finite detection threshold when analyzing the original signal, which leads to correlated bursts of activity by separating the avalanche events into subavalanches. Effects of thresholding have also been considered in some other types of signals, e.g., in the context of fracture [142–144], and argued to be of importance in seismicity [145,146].

Driving mode represents one of the essential factors influencing the avalanche-like response. Theoretically, the three types of driving are recognized: adiabatic, in which the external magnetic field is increased for the exact amount needed to trigger the least stable spin and kept unchanged until the avalanche stops; quasistatic, in which the incrementing of the external magnetic field is performed in fixed steps until the conditions for the initiation of the avalanche are met, from which moment is kept constant as long as the avalanche propagates; and a finite rate driving during which, throughout the whole simulation time, the external magnetic field is increased at a constant rate. According to the value of the driving rate, the regimes of slow, intermediate and fast driving are identified [130]. The underlying dynamics in the finite driving regime is profoundly influenced by the interplay of the system’s disorder and the driving rate at which the magnetic field is incremented. In this type of driving, multiple avalanches may simultaneously propagate, making the analysis and interpretation of the obtained data of the avalanching dynamics far more complex.

The different ways of driving have a significant impact on the system’s behavior, with one of the most noticeable effects being the time/space profile of avalanche evolution. Adiabatic driving, during which only one avalanche is active at a time, is conducted under precisely defined conditions realizable in the numerical simulations but not in the actual experiments. A little bit closer to realistic is the quasistatic driving type, in which the regimes with small (adiabatic-like) and large field increments can be identified, promoting the temporal and occasionally spatial merging of avalanches [147]. This amalgamation of avalanches is expressed the most at the fast regime of finite rate driving protocol, comprising an overall system-spreading activity [148] without the possibility of distinguishing the contribution of individual avalanches [130,147–149]. This type of driving is even more realistic and comparable to the experimental situations. In one of the previous studies [150], it has been shown that due to the merging of avalanches (swelling) the overall duration of avalanches can be extended; at the same time, the merged avalanche will appear in the distribution while the merging avalanches vanish (merging absorption). Similarly, for the avalanches overlapped in time (and spatially separated), the overlapped avalanche appears, whereas the overlapping avalanches disappear from the distribution (temporal absorption). Until now, many experiments have been conducted on disordered ferromagnetic materials using the finite driving rate protocol for example, in [14–16,133]. Works in [26,150–153] attempt to describe theoretically the observed critical dynamics phenomena with the finite driving rate and its impact on the hysteresis loop and coercive field was considered in [18,19,154–156].

The study of the effects that the type of driving inflicts on the system is further augmented by employing the newly proposed stochastic driving, realized by stochastic increments of the external magnetic field in each time step of the system’s evolution, regardless of the current state of the ongoing system’s activity. This type of driving opens the possibility to extend the findings and perform the more realistic interpretation of seismic activity [50], owing to its importance due to the unwanted ramifications. In some recently conducted experiments on the propagation of a crack line in a random environment [57,142,157,158], a seismic-like behavior has been demonstrated, also found in the studies of compressional fracture [159]. Stochastic driving represents so far the most realistic scenario, allowing much closer benchmarking with the results of the experimental studies and the possibility to extend to the numerous models developed in an attempt to explain and interpret this complex avalanche dynamics [143,160–163].

Another aspect concerns the thermal effects and the presence of the demagnetization fields in the system. These physical factors bring challenges, especially in understanding the intricate interplay of various parameters, interpreting and implementing the solution numerically, and comparing the numerical results with experiments. As shown in [164], the extended form of RFIM can be an excellent alternative for the micromagnetic modeling, which is particularly relevant for applications. The stochastic nature of the thermally triggered spin flips tends to impede accurate avalanche detection by contaminating the signal, producing an excess of minor activity events, and limiting the spread of large avalanches. Nonetheless, the demagnetizing field limits the amount of underlying spin activity events, and the maximum length of magnetization jumps by counteracting the external field and introducing non-local interactions. In a recent study [164], an appropriate algorithm to study these combined effects was developed. Notably, it demonstrates that the non-zero temperature modifies the remanent magnetization and coercive field, leading to the shrinkage of the hysteresis loop and the amplification of minor activity events. By introducing extended linear segments, the higher demagnetizing coefficient modifies the loop shape and changes the multifractal nature of the magnetization fluctuations. Similar to the analogous ZT dynamics, the statistics of well-identified intermediate-range activity events are controlled by the same scaling exponents.

In the next section, we describe the most general version of the NEQ RFIM to study the critical behavior on the hysteresis loop of disordered ferromagnetic systems, assuming different driving conditions and the aforementioned physical factors, followed by the section presenting the appropriate computational techniques. The rest of the review is structured so that each section presents the results obtained in numerical studies of particular hysteresis phenomena, emphasizing the influence of specific physical factors.

2  Random Field Ising Model (RFIM): Definition, Parameters and Simulation Approaches

The Random Field Ising Model describes systems of Ns mutually interacting classical Ising spins si=±1 exposed to external magnetic field H and a quenched local random magnetic field hi. The spins are situated at the sites i of some underlying finite lattice with space dimensionality d, shape and size specified in each model realization. Mostly employed are the two-dimensional (2D), and three-dimensional (3D) lattices, however lattices of other dimensions are also used (for d>3, see, e.g., [31]). Frequently utilized are finite lattices cut out of a corresponding infinite lattice with translation symmetry like the quadratic lattice of size (Lx,Ly) in 2D, cubic lattice of size (Lx,Ly,Lz) in 3D, and analogously hypercubic lattice for d>3. In all three of the foregoing lattice examples the number znn of nearest neighbors for each site equals 2d, however, this number can be different for other types of lattice elementary cell, like znn=6 for the triangular [100], and znn=3 for the honeycomb [103] 2D lattices. Finally, the model can be situated on a lattice without translation symmetry, such as Bethe lattice, and/or translation symmetry may be absent due to the presence of irregularly distributed vacancies (lattice sites not ‘occupied’ by a spin), and/or irregular interfaces consisting of spins with fixed value.

In general, the interaction of each spin si→ with the external magnetic field is accomplished through the coupling −μ0mi→⋅Hi→(t) between the magnetic dipole moment mi→=γi(s)ℏsi→ associated with si→, and the external magnetic field Hi→(t) at the spin’s site i at the current moment (of time) t; here, γi(s) is the gyromagnetic ratio for the spin si→, usually taken to be the same for all spins in the model, γi(s)=γ(s). Because the RFIM is not a vector spin model, one can assume here that some (arbitrary) direction is chosen in space such that for all spins the projection s→i⋅u→0 of si→ onto the unit vector u→0 of this direction equals si enabling in the RFIM the incorporation of the coupling between spins and the external magnetic field as −Hi(t)si; here, Hi(t)=H→i(t)⋅u→0 is the projection of H→ onto u→0, and the units are appropriately rescaled for simplicity (which is allowed in models).

Analogous form −hisi is taken for the coupling between the spin si and the random field hi at the spin’s site i. This local random field {hi}i=1Ns is postulated as quenched (meaning that its values are fixed in time), and could be considered to originate from various types of real systems’ imperfections not specified in the RFIM. Besides, this field is random, meaning that its values hi at different lattice sites are chosen out of some type of zero-mean random distribution ρ(h) with finite standard deviation R measuring disorder in the RFIM systems, which could be site-dependent. Here, the common approach is to use Gaussian distribution, ρ(h)=1R2πexp⁡(−h2/2R2), with the same disorder R for all sites. The RFIM studies, performed with another choice of the random filed distribution ρ(h), exist but are relatively rare-see, e.g., the case of parabolic and uniform distributions in [87], Laplace and double Gaussian distribution in [97], and the case of Gaussian distribution with R modulated by the presence of crystal grains [165].

Independent on the choice of the distribution ρ(h) is the question whether there is some correlation between the values of random field at different lattice sites. When the values are uncorrelated, then ⟨hihj⟩=Ri2δij, where δij is the Kronecker delta function and Ri is the value of disorder at the site i. This is the case in majority of past RFIM studies although the model itself doesn’t exclude possible correlations.

Besides the interaction of spins with magnetic fields (external and random), the spins also mutually interact. The strongest type of interaction between the spins is the (electrostatic in nature) exchange interaction between pairs of spins, giving the contribution −Jijsisj from each pair si and sj, to the overall system Hamiltonian through the coupling constant Jij. When Jij>0 the coupling is ferromagnetic, for Jij<0 the coupling is antiferromagnetic, while for Jij=0 the spins si and sj are exchange decoupled. The exchange interaction is typically considered as short-ranged, meaning that it is the most prominent for the nearest neighbor spins and significantly weaker for further neighbors (e.g., exponentially decaying with the inter-site distance). On this ground, frequently used is the nearest-neighbor approximation in which all pairs of spins except nearest neighbors are considered as exchange decoupled.

Another type of spin-spin interaction is the dipole-dipole interaction between the magnetic dipole moments m→i and m→j associated with pair of spins s→i and s→j. In general, the dipole-dipole interaction is of the form

−μ0(γ(s))2ℏ24πrij5[3(si→⋅rij→)(sj→⋅rij→)−si→⋅sj→rij2],

where rij→=ri→−rj→ is the position vector of site i relative to site j. In the RFIM this expression is simplified to

−Jdipole3cos2θij−1rij3sisj,

where Jdipole=μ0(γ(s))2ℏ2/4π, and θij is the angle between rij→ and u→0. Being inversely proportional to the cube of inter-site distance rij, the dipole-dipole interaction is not (like the exchange interaction) limited to near spins, but instead, it affects all pairs of spins, i.e., it is long-ranged. Yet, it is typically weaker than the exchange interaction, and therefore often omitted in the model analysis.

Besides, the spins located at the boundaries of finite lattices generate inside the sample a long-range (effective) demagnetizing magnetic field which acts against the external field. In the RFIM, this field is taken as a homogeneous field −JDM and its coupling with individual spins as JDMsi, where M=(∑isi)/N is the actual magnetization of the system and JD≥0 is the demagnetization coefficient (factor). The adopted form of demagnetizing field is by all means approximate, so in simulations and analyses it is most appropriate to assign such values to JD that conform to the values theoretically found for the homogeneously magnetized samples (e.g., JD=1/3 for cubic and spheroidal samples).

Therefore, the most general form of Hamiltonian for the RFIM spin systems in the external magnetic field with a time profile Hi(t) reads

ℋ=−∑{ij}(Jij+Jdipole3cos2θij−1rij3)sisj−∑i=1Ns(Hi(t)+hi−JDM)si,(1)

where the summation in the first term is performed over all pairs of distinct spins, so that each spin si is under the influence of the effective magnetic field

hieff(t)=∑j≠i(Jij+Jdipole3cos2θij−1rij3)sj+Hi(t)+hi−JDM,(2)

giving contribution −si(t)hieff(t) to Hamiltonian (1). If si(t)hieff(t)≤0, the spin si(t) is field-unstable at the moment t, and the value of Hamiltonian (1), henceforth simply the system energy, will be reduced by flipping of si (i.e., change of si to the opposite value −si). Otherwise, if si(t)hieff(t)>0, the spin si(t) is field-stable, so its flipping would increase the system energy.

Although absent in Hamiltonian (1), the lattice-spin interactions are indirectly present in the RFIM through the thermal fluctuations of spins manifested when the system temperature Θ>0 in which case the model is called thermal. Otherwise, when Θ=0, the thermal fluctuations of spins are absent and the model is athermal or zero-temperature. This Θ=0 version of RFIM is of twofold importance because in many real-world systems, to which the RFIM can be applied, the thermal fluctuations are negligible in the relevant range of temperatures, and also because the athermal RFIM version is much simpler than the thermal, both for simulations and analyses.

In the RFIM, the state {si}i=1Ns of the spin system changes through the flipping of individual spins. In the nonequilibrium (NEQ) model version, which is the subject of this review article, the individual spins flip due to appropriate change of their effective magnetic field and/or thermally if Θ>0. Thus, in the athermal NEQ RFIM version, only the first cause applies, and this is usually done by parallel updating of orientation of spins, meaning that all spins which are field-stable at the moment t will remain unchanged and the field-unstable spins will be flipped at the next moment t+Δt (in simulations, both t and its increment Δt are discrete, usually with Δt=1). Due to such local dynamical rule and parallel updating, the system evolves deterministically, traversing (almost exclusively) through the nonequilibrium states in a tendency to reduce its energy. Other (and less used) way, employed in the athermal model version, is to perform sequential field-stability testing and flipping of spins one by one, in which case the system’s evolution depends on the random order in which the spins are tested.

Flipping of spins, together with a possible change of the external magnetic field, modifies the effective magnetic field hieff of the system’s spins. Thus, when JD>0, the magnetization M is (typically) altered due to spin flipping which causes a change in hieff for all spins through the term −JDM in (2). Also, the dipole-dipole interaction term in (2) becomes different for all (especially near) spins, provided that this type of interaction is included in the model. However, the contribution coming from the change in the exchange-coupling term in (2) is the most prominent, or the only one when JD=0, Jdipole=0, and constant H at the next moment t+Δt.

In thermal (i.e., Θ>0) model, together with the preceding field-stability testing, the spins are also checked for thermal flipping. Thus, in [166], the changes ΔH of the external magnetic field are sandwiched between a preselected number Nth of Monte-Carlo sweeps sequentially applied, so that in each sweep a spin si is randomly chosen and flipped with the probability min(1,e−ΔEi/Θ), where Θ>0 is the temperature, and ΔEi is the change in energy proposed by flipping of si. By this rule, the selected spin si is certainly flipped if it is field-unstable (ΔEi<0), however, if it is field-stable (ΔEi>0), it may be also flipped with probability e−ΔEi/Θ. As a consequence, the central role is played by the time-scale of thermal flipping, set by the number Nth regardless of the system size, while the time-scale for the field changes is introduced indirectly, relative to the thermal flipping time-scale.

To check the stability at the current moment for all (and not only selected) spins, a different approach is proposed in [164]. There, the thermal flipping is tested on a set of spins, randomly selected at each moment t and containing a preselected fraction c=ΔNs/Ns of spins. At the next moment t+Δt, each spin outside this set will be flipped if it is field-unstable, while each spin si from the set will be flipped with the (Boltzmann-type) probability of thermal flipping probability

pi(th)=exp⁡(−sihieff/Tr)exp⁡(−sihieff/Tr)+exp⁡(sihieff/Tr),(3)

where the parameter Tr=Θ/Tc stands for the temperature relative to some temperature Tc that is characteristic in the model for the underlying system (e.g., Tc=J for the system of purely ferromagnetic spins with Jij=J). Hence, for sihieff<0 (i.e., if the selected spin si is field-unstable) it is likely (but not certain) that it will be flipped, whereas for sihieff>0, the selected spin si has a small nonzero chance for flipping.

As stated, the external magnetic field varies in time in the NEQ RFIM. This is realized either in some deterministic variation pattern (e.g., adiabatic, quasistatic, finite driving rate), or by stochastic external field increments (see, e.g., Reference [149]). In the adiabatic and quasistatic driving protocols, employed in the athermal model versions, the external field is changed only if all spins in the system are field-stable. In the adiabatic case the field is changed in the exact amount that destabilizes the least stable spin (and flips it at the next moment), while in the quasistatic driving the field increment is constant, ΔH=const, possibly causing flipping of one or more spins, especially for the large systems. Flipping of spin(s) may destabilize some spins (most often the nearest neighbors of the flipped spin); their flipping may further destabilize other spins, and so on, leading to a time-series of flipping in the form of an avalanche. In the case of adiabatic models with JD=0, Jdipole=0, and the exchange coupling limited only to the nearest neighbors, each avalanche is nucleated by flipping of a single spin and afterward spread over a cluster of (connected) spins. In other adiabatic cases, or when the driving is quasistatic, multiple avalanche nucleations at spatially distant locations may happen. This leads to the onset of avalanches that simultaneously propagate over several spin clusters initially space-separated and later possibly merging in space. Whether due to such, or to individual avalanches, the system response is expressed in terms of the numbers n+(t) and n−(t) of spins that flip up and down at the moment t, giving the response signal of the system, V(t)=n+(t)−n−(t), and the corresponding change in magnetization ΔM(t)=2V(t) during one time-step at the moment t. Each subsequence of the V(t) sequence, consisting of consecutive terms such that the underlying n+(t) and n−(t) satisfy n+(t)+n+(t)>0, represents a system activity event, which is a generalization of an individual avalanche.

In addition to adiabatic and quasistatic driving (which tends to adiabatic when ΔH→0 for any finite system), there is the finite rate driving protocol in which the external magnetic field, regardless of the system activity, increases/decreases along the rising/falling part of the magnetization curve at some constant driving rate Ω=ΔH/Δt (given by the field increment/decrement ΔH when Δt=1 and also tending to adiabatic when Ω→0). This type of driving, enhancing spin flipping and thus facilitating activity propagation due to perpetual modifications of spins effective magnetic field, is of considerable importance because it is encountered in a large number of experiments performed on magnetic systems.

Compared to athermal model, the system response is essentially different in the thermal model regarding events duration. Thus, while all events in athermal model are of finite duration this is almost never the case in the thermal model due to the thermal flipping of spins. For this reason, unless some threshold is introduced, the decomposition of system response into activity events separated in time is impossible, making the adiabatic and quasistatic driving unrealistic. So, in [164], the external magnetic field is incremented in each step, whereas in [166] the field is incremented after Nth (elemental) time-steps, while meanwhile the thermal flipping of spins is performed. Note that in pure ferromagnetic case, due to thermal flipping and other factors influencing the effective magnetic field, some spins may be back-flipped, i.e., changed from +1 to −1 on the rising part of magnetization curve (and opposite on the falling part). Such back-flipping is possible when JD>0 even in athermal model and impossible when JD=0 and Jdipole=0.

The analysis of RFIM systems evolution is commonly performed on their response signal decomposed into events; the exception is the analysis of power spectra P(f) [167] giving the variation with frequency f of the released power frequency density which is performed without any decomposition on the entire selected part or on the whole response signal. In the athermal version of ferromagnetic systems each activity event gives a part of consecutive non-zero values in the response signal, V(t)>0 along the rising part of the magnetization curve (and V(t)<0 along the falling part). These parts are taken as events which are the subject of further statistical analysis. In other cases, such simple decomposition is impossible, and the decomposition of the response signal into events is accomplished with the aid of some (suitably chosen) discrimination threshold Vth. Provided that Vth≥0 is selected, the response signal is decomposed into longest subsequences {V(ti),V(ti+1),…,V(te)} consisting of signal values |V(t)|>Vth, all positive (or all negative), and registered in consecutive moments between the initial moment ti and ending moment te of the subsequence in question. Each such subsequence is considered as an event, and for any moment t that is outside all events |V(t)|≤Vth. Each event can be characterized by several parameters like its duration T=te−ti+1, size S=∑te≤t≤tiV(t)Δt, energy E=∑te≤t≤tiV2(t)Δt, and amplitude A=max{|V(ti)|,|V(ti+1)|,…,|V(te)|}, and for each such parameter X the corresponding distribution DX(X) is collected. These distributions can be either integrated (i.e., collected from the entire signal), or windowed (i.e., collected in a selected window (Hmin,Hmax) of the external magnetic field). Besides, the distributions can be classified according to the event’s type (e.g., nonspanning and different types of spanning events). Among possible distribution shapes, of particular importance are those of the type

DX(X)=X−a𝒟X(X),(4)

representing a power-law with cutoff specified by the power-law exponent a and the so-called cutoff function 𝒟X(X) which (beside X) depends on some parameters for simplicity not shown in notation. The analytic form (4) describes a pure power-law, DX(X)∝X−a, only in some limits of the cutoff function parameters leading to 𝒟X(X)→const. Outside such limits the cutoff function rapidly tends to zero for large X (and possibly for very small X), while in between, a region of moderate X-values may exist such that 𝒟X(X)≈const. If so, this region is called the scaling region, and it could be used for determination of the value of the power-law exponent as a equals the gradient (i.e., slope) of the scaling region data presented in a log-log plot; other (and preferred) method is by the collapsing of data which can be applied when the cutoff function has universal scaling properties as will be demonstrated in the Results sections. For the later convenience, here we quote the common RFIM exponents: τ, α, ε, and μ for the windowed distributions of event’s size, duration, energy, and amplitude, respectively; additional NEQ RFIM exponents will be introduced in the Results sections.

Besides, other quantifiers of response signal can be defined, like the magnetization as a function of the external magnetic field, M=M(H), and analogously the susceptibility dM/dH, the correlation functions, average event shape, and the distributions of various types of waiting time, like the external waiting time Text measuring the separation in time of consecutive events.

Integral part of model is the specification of boundary conditions which are either closed or open on each system boundary. Closed conditions mean that the spin values are periodic along the corresponding direction, e.g., s(ix,iy,iz)=s(ix+Lx,iy,iz) in the case of (Lx,Ly,Lz) 3D cubic lattice and closed boundary condition along x− direction. Closed boundary conditions are important because they enable faster convergence towards the infinite systems. On the other hand, they are less realistic than the open boundary conditions for which there are no spins outside considered lattice.

When studying system evolution, initial conditions have to be specified that include the initial state for all spins and the external magnetic field. If suitable, terminal conditions are specified as well.

Finally, let us explain the role of averaging. Averaging is performed in order to collect more reliable statistics of system response. However, for fully deterministic variants of the model repeating the simulation under identical conditions is useless because it gives exactly the same results. In such cases, the so-called quenched averaging is performed in which the simulation is repeated using different configurations of the random field corresponding to the same choice of disorder and the same values of the remaining model parameters. Otherwise, the simulations may be repeated using the same configuration of the random field but with randomized other thermal ingredients, e.g., the different sets of spins tested for thermal flipping.

The theoretical analyses of the NEQ RFIM was focused so far on the possible critical behavior of this model. Thus, for the adiabatically driven ZT model situated on the (hyper)cubic lattices, the renormalization group (RG) analyses [28,37,84] have revealed a mean-field behavior for dimensions d≥6 and a nontrivial critical behavior for 3≤d<6, but without the final RG conclusions about the 2D case.

Having explained the main features of various versions of the NEQ RFIM, let us briefly state the main characteristics of the equilibrium model. In this RFIM version at each external (commonly homogeneous) magnetic field of interest is only the ground spin state (or states in the case of frustration when more than one state have minimal energy) at Θ=0, whereas for Θ>0, the corresponding thermal distribution of spin states is analyzed. Hence, the temporal evolution of the underlying spin system remains out of the scope together with all associated features (events, power-spectrum, waiting times, etc.). For this reason, the equilibrium version is less informative than the nonequilibrium version, especially for the analyses of response of the real-world systems. Nevertheless, the equilibrium version is very important in theoretical studies revealing many conclusions that are not attainable within the nonequilibrium version.

As mentioned in the Introduction, the focus of this review is on extensive computational techniques and simulations of the field-driven spin reversal dynamics, avalanches statistics, finite-size scaling, correlations and multifractal analysis of Barkhausen noise signals, that are developed to study the nonequilibrium criticality of the hysteresis loop. They are detailed throughout the results presented in different sections. More precisely, in the following section, we give the fundamental aspects of the simulations, which are then adapted and elaborated in each section to analyze different driving modes, dimensionality and sample shapes, athermal and thermal fluctuations and demagnetizing effects. See also the program flow in Section 6.1, which refers to the model with demagnetizing fields. In addition, a systematic comparison of the analysis of simulated data and data collected in the experiment with the nanocrystalline sample in [168] is presented; see Section 7.1. It demonstrates the advantages and disadvantages of presented numerical experiments compared to laboratory experimental data.

3  Simulational Methods in the RFIM

Even in the simplest model version, the RFIM simulations are computationally very challenging. Realization of the finite-size scaling analysis, and avoiding the finite-size effects, demand simulations of very big systems (e.g., with 109 spins) which require large amounts of computer memory and long running time even for a single run, let alone for repeated simulations with different realizations of the random magnetic field necessary for quenched averaging. Efficient algorithms are therefore a must regarding both memory (RAM and storage) space and execution time.

An essential breakthrough in simulations is achieved by the sorted-list algorithm and bit-per-spin algorithm detailed in [169] for the adiabatically driven athermal NEQ RFIM (with Jij=J>0 for nearest neighbors and zero otherwise, Jdipole=JD=0, and homogeneous H). In this case, the external magnetic field H is to be increased only after all spins become stable and that by the amount triggering the least stable spin. To find it, the easiest tactic is to check the values of the effective magnetic field hieff at all lattice sites i, find the site i0 with the least negative value of hieff, increase H by ΔH=−hi0eff, and flip the spin si0. Then, iteratively until the system becomes stable, for all spins flipped at the previous moment, flip at the next moment their nearest neighbors that became unstable. Despite being simple, the preceding algorithm, named the brute-force method, is extremely slow because its total running time scales as O(Ns2); therefore it is not used in simulations of larger systems (with, e.g., more than 105 spins) even by modern (sequential) computers.

Traversing the whole lattice in search for the least stable spin is avoided in the sorted-list algorithm. To this end, the array {hi}i=1Ns of the random field values, generated before the start of the simulation, are sorted in descending order, {hπ(1)>hπ(2)>…>hπ(Ns)}, where π:{1,…,Ns}→{1,…,Ns} is the permutation such that π(k) points to the index of the element in the original array occupying k-th position in the sorted array h∘π. Also is introduced a pointer (integer) array {Pn↑}n↑=0znn and all its elements initially set to 1. After the system becomes stable at the current value of H, each element Pn↑ of this (nondecreasing) array is updated (increased) so to point to the smallest index k in the sorted array, i.e., to the largest element in the sorted array h∘π, such that sπ(k) would flip if it had n↑ upward oriented nearest neighbors, and the external field for that n↑ set to Hnew(n↑)>H equal to Hnew(n↑)=−[(2n↑−znn)J+hπ(k)]. Thereupon, the external magnetic field is updated to the new value Hnew=min{Hnew(n↑):n↑=0,1,…,znn}, and the new avalanche nucleated by flipping the least stable spin caused by H→Hnew.

Having the total running time that scales as O(Nslog⁡Ns), the sorted-list algorithm is, so far as is known, the fastest algorithm used for simulations of the adiabatically driven athermal NEQ RFIM. However, the algorithm is memory-hungry because, besides the arrays used for storing the spins and the random magnetic field, it requires an additional integer array {πi}i=1Ns keeping track about the permutation used in sorting the original array {hi}i=1Ns, requiring for storage at least 64 bits of memory per each element.

The bit-per-spin algorithm greatly reduces memory demands because it eliminates the need for the array {hi}i=1Ns. Besides this, memory is additionally saved by storing each spin in only one bit of memory, also possible in any other algorithm regarding the storage of Ising spins. Although the algorithm’s historical name refers to bit-per-spin memory storage, the main aspect of its optimization is generating the random field values on the fly, resulting in at least 96 times reduced memory demands, but at the cost of approximately two-fold increase in the execution time. Thus, the array {hi}i=1Ns is not generated before the start of the simulation, nor the memory for its storage is reserved. Instead, provided that the system is stable at the external field Hold, a new value of the external field, Hnew is found. To this end, using the complementary error function erfc=(2/π)∫x∞exp⁡(−t2)dt, one firstly calculates the probabilities

P↓(n↑,H)=12erfc(H+(2n↑−znn)JR2),(5)

that at the value H of the external field a spin with n↑ up neighbors is pointing down, and consequently the probability

Pnone(Hold;Hnew)=∏n↑=0znn[P↓(n↑,Hnew)P↓(n↑,Hold)]Nn↑,(6)

that no spin has flipped between Hold and Hnew. After that, the equation

rn=Pnone(Hold;Hnew),(7)

is numerically solved for a chosen random number rn (uniformly distributed between zero and one), using in this solving the recipe from [169], for example. In the next step, the rates for spin flipping for each n↑ at the Hnew are calculated according to

Γ(n↑,Hnew)=Nn↑ρ(Hnr(n↑,Hnew))P↓(n↑,Hnew),(8)

where Nn↑ is the number of spins with n↑ up neighbors (taken before H is set to Hnew), and ρ(Hnr(n↑,Hnew)) stands for the Gaussian distribution ρ(h)=exp⁡(−h2/2R2)/R2 of the random field at h=Hnr(n↑,Hnew), where Hnr(n↑,Hnew)=Hnew+(2n↑−znn)J. For these rates Γ(n↑,Hnew), the total rate

Γ(Hnew)=∑n↑=0znnΓ(n↑,Hnew)(9)

is calculated and n↑ selected using another random number uniformly distributed between zero and Γ(Hnew). Finally, randomly searching the lattice, a down spin with n↑ up neighbors is found and flipped. This way a new avalanche is nucleated and its propagation at Hnew is followed like in the brute force method until the ongoing new avalanche dies. The whole preceding procedure (for finding Hnew and propagating the emerged avalanche) is iteratively repeated, until all spins become stable in the up state.

The bit-per-spin algorithm enabled large-scale simulations in the last decade of the XX century standing as a breakthrough in the numerical studies of the model [31,170]. Together with the sorted-list algorithm, it was later extended to be applicable for the quasistatic [147], finite-rate [130,148,171,172], and stochastic driving [149] along the entire hysteresis loop. Nevertheless, their extension was impossible for the RFIM versions modelling the systems at finite temperatures and/or systems having a more complex inter-spin interaction (e.g., nonhomogenous exchange interaction, combined ferromagnetic-antiferromagnetic layers, vacancies, interfaces, etc.). In these cases, the single run execution time is possible to be reduced by converting the brute-force method from the sequential to parallel code, which is particularly convenient for the shared-memory systems using Message Passing Interface (MPI)–a standard designed to function on parallel computing architectures. Note, however, that the sequential code provides the shortest run-time per thread, and that the run-time acceleration usually rapidly saturates with the number of involved parallel threads. This, together with the repetition of simulations with different random field configurations requested for the quenched averaging, means that the number of involved parallel threads has to be compromised to achieve the least overall execution time of the whole set of simulations. An example of the flowchart of a parallel algorithm is given in Section 6.1.

Besides sophisticated simulation algorithms, the RFIM studies also have computational challenges regarding the extraction of the relevant statistics and their subsequent analysis. As an illustrative example, let us describe the algorithm proposed in [89] for the collapsing of the suitably scaled distributions. These distributions are of the power-law type and are represented by their histograms collected in the logarithmic bins which reduces the random scattering at the side of large events. Thus, for a given set of discrete histograms, one firstly finds the union set X of all their ordinates, next for each x in X the (weighted) mean of the interpolated curves, and then the width around this mean curve, resulting in a merit function (i.e., width w). Various types of interpolation and weighted widths can be used, which is important for noisy data having different uncertainties. The simplest choice is the linear interpolation, no weights, and width w=D2/Np, where D2 is the sum of squared distances from the mean curve, and Np is the number of degrees of freedom (crudely equal to the number of points in X).

4  Adiabatic Collective RFIM Dynamics with Spin-Flipping Avalanches

4.1 Zero-Temperature (ZT) Nonequilibrium (NEQ) RFIM on Three-Dimensional (3D) Simple Cubic Lattice

Previous theoretical and numerical studies [28,35,84,89] have shown that the adiabatically driven ZT NEQ RFIM systems, situated at the equilateral (hyper)cubic lattices with dimension d≥2, exhibit in the thermodynamic limit (L→∞) a dynamical critical behavior at the (d-dependent) value of critical disorder Rc(d), discriminating two domains of disorder. This value, denoted as Rc in what follows, separates the ferromagnetic phase that exists for low disorders R<Rc, at which the infinite avalanche appears causing a jump in magnetization, from the paramagnetic phase occurring for high disorders R>Rc, in which no infinite avalanche is created, so the magnetization curve is smooth.

In finite systems, the role of infinite avalanches is played by the spanning avalanches [27] spreading along at least one of the system’s spatial dimensions and, therefore, classified according to the number of dimensions they span as 1d, 2d, and 3d spanning avalanches (further classified in [27] into critical and subcritical 3d spanning avalanches [173]). As shown in the comprehensive work investigating the ZT NEQ RFIM on finite equilateral 3D cubic lattices [108], the distributions of the number Ns(R;L) of all spanning avalanches per single run in the systems of size L, illustrated in the main part of Fig. 1a, scale as

Figure 1: The total number of spanning avalanches per single run, Ns(R;L), and their scaling collapse (10) are displayed by symbols in the panel (a) and its inset. The model curves (11) that best fit the Ns(R;L) data are shown by the full lines in the main part of this panel. (b) The histograms of spanning field Hsp for 1d, 2d and 3d spanning avalanches obtained for system with size L=128 and disorder R=2.22. Presented data are averaged over 40,000 different realizations of the random field. Figure is replotted from Reference [108], combining parts of Figs. 1 and 2

Ns(R;L)=LθN~s(rL1/ν),(10)

see in the inset of Fig. 1a, where N~s(rL1/ν) is the universal scaling function employing the reduced disorder r=(R−Rc)/R, and θ is the exponent whose value is, after considering the significantly wider range of the employed values of L, with great accuracy refined in [108] to θ=0. Additionally, it was shown in [108] that the distributions Ns(R;L) are well-described by the model function

Nsmod(R;L)=Hexp⁡{−[(R−Rceff(L)]2/2W2(L)}+0.5erfc{[(R−Rceff(L)]/W(L)},(11)

depending on parameters H, Rceff(L)>Rc, and W(L) which specify height, center, and width of the Gaussian in the first addend of the Eq. (11). Besides, Rceff(L) and W(L) in the second addend specify the inflection point and the width of the complementary error function taking values 0.75 and 0.25 for R=Rceff(L)−W(L) and R=Rceff(L)+W(L), and 0.5 at the inflection point R=Rceff(L). For large L, the width becomes very small tending to zero as W(L)∼L−1/ν, while Rceff(L)→Rc according to [Rceff(L)−Rc]/Rceff(L)∼L−1/ν. As long as R is significantly below Rceff(L)−W(L) only one (3d) spanning avalanche appears per run, so that Ns(R;L)≈1. Therefrom, the number Ns(R;L) increases when R grows and reaches its maximum at R=Rceff(L), followed by a rapid decrease towards zero as W(L) is small when L is large. Owing to such behavior, Rceff(L) is singled out as a characteristic value of disorder that can be considered as the effective critical disorder for the systems of size L such that (in simplified terms) the spanning avalanches are absent for R>Rceff(L) or (more precisely) unlikely for R surpassing Rceff(L)+W(L). As a consequence, three domains of disorder with distinct scaling exist for finite systems with negligible W(L), namely the domain of low (or subcritical) disorders for R<Rc, the domain of transitional disorders for Rc<R<Rceff(L), and the domain of high (or supercritical) disorders for R>Rceff(L), taking into account that for smaller systems the borders of the domains are not very sharp due to finite W(L).

Results published in [108] are obtained in by-all-means extensive numerical simulations, averaged across up to 200,000 distinct realizations of the random magnetic field for up to 36 various values of disorder on the systems with lattice size up to L=2048 containing almost 1010 spins. These results showcase that for disorders R>Rceff(L)+W(L), surpassing the effective critical disorder, the systems’ properties are almost independent on the lattice size L, practically obeying the scaling predictions that hold in the thermodynamic limit. However, significant size-dependence is documented in the other two disorder domains, subcritical and transitional, evident in all characteristic features of systems’ behavior.

In the subcritical and transitional domains of disorder, the single-run magnetization curves have jump(s) that appear at the spanning field Hsp (i.e., the value of the external field at which a spanning avalanche is triggered) whose values are determined by the chosen RFC. As the distribution of Hsp values (due to varying employed RFCs) is not sharp but of finite width, see Fig. 1b, the positions of magnetization jumps are smeared, cf. bottom-right inset in Fig. 2b. Therefore, the corresponding part of the averaged magnetization curve is slanted, see top-right inset in Fig. 2b, gradually turning into vertical when L→∞, cf. right inset in Fig. 2a. On the other hand, for finite systems in the supercritical domain of disorder, both single-run magnetization curves and their averaged counterparts, exemplified in the right inset of Fig. 2c, are smooth functions of the external field H.

Figure 2: Main panels show the scaling collapses (12) of the averaged magnetization curves in the domains of disorder: (a) below critical, (b) transitional, between the critical and the effective critical, and (c) above the effective critical disorder, for a constant value of L|r|ν and the values of critical exponents and parameters presented in Table 1. The corresponding left insets show the collapses (13) of the averaged magnetic susceptibility curves, while the right insets show examples of the averaged magnetization curves. Bottom-right inset of panel (b) shows magnetization jumps in the two single-run magnetization curves, illustrating that more than one spanning avalanche (together with the associated Hsp-values) may appear in a single run. Presented data are replotted from Reference [108], combining parts of Figs. 6, 8 and 9

Averaged magnetizations and magnetic susceptibility curves for finite systems of lattice size L, according to [27,108], follow the scaling forms

MR,L(H)=|r|βℳ~(h^Reff/|r|βδ,1/L|r|ν),(12)

χR,L(H)=|r|β−βδ𝒳~(h^Reff/|r|βδ,1/L|r|ν),(13)

where h^Reff=H−Hceff(R) is the effective reduced magnetic field, measuring the shift of the data along the H-axis by the value of the effective critical field Hceff(R) pinpointing the pertinent susceptibility maxima; Mceff(R)=M(Hceff(R)) is the effective critical magnetization, ℳ~ and 𝒳~ stand for the universal scaling function, different in each of the disorder domains, and β, βδ and ν are the standard RFIM exponents [28]. Scaling collapses of averaged magnetizations and susceptibilities in all three domains of disorder, each having parameters L, R chosen so that the value of L|r|ν is constant and allowing for the estimation of critical parameters and exponents as the best collapsing values (shown in Table 1), are presented in Fig. 2.

The distributions of avalanche parameters are cutoff-ending power-laws of the distributed parameter which depend on the underlying values of disorder R and lattice size L. Considered as the generalized homogeneous functions of their arguments, these distributions should follow finite-size scaling predictions in three equivalent forms, two of which, for the integrated size distributions and according to [27,108], read

DS(int)(S;R,L)=L−τ′/σνD~S±(S/L1/σν,1/L|r|ν),(14)

and

DS(int)(S;R,L)=S−τ′𝒟S±(Sσ|r|,1/L1/ν|r|),(15)

where τ′ is the pertinent critical exponent of the integrated distribution of avalanche size, τ′=τ+σβδ, and σ is the cutoff exponent describing the scaling Smax∝|r|−1/σ [28]. These forms should hold in all domains of disorder, with the domain-specific universal scaling functions D~S± and 𝒟S±, predicting data collapsing for any set of size distributions satisfying the condition L1/νr=const. The results of numerical simulations [108] supported the preceding scaling conjectures in the supercritical domain of disorder, as is illustrated in Fig. 3a by the data collapsing (14) in the left, and (15) in the right inset. Additionally, these simulations also showed that the DS(int)(S;R,L) distributions become almost L-independent for R>Rceff(L)+W(L), see Fig. 3b, meaning that they obey the scaling

Figure 3: (a) Presented integrated size distributions are all from the supercritical domain of disorder and correspond to the (L,R) pairs from the legend satisfying the condition Lrν=const. While the main part of this panel shows the original distributions, the left inset shows their scaling collapse according to (14), and right inset according to (15). (b) Integrated distributions DS(int)(S;R,L) of avalanche size S in the main panel are obtained for several values of L quoted in legend and the same disorder R=2.3 surpassing Reff(L)+W(L) for L≥128, except for L=64. Thus, being in the supercritical domain, all L≥128 distributions mutually overlap, while the case L=64 manifests a bump before the large-S cutoff due to the presence of spanning avalanches. Inset illustrates the same effect for the susceptibilitity curves dM/dH. (c) Attempts to collapse the integrated size distributions of nonspanning avalanches using (16) for the system with size L=360 in a broad disorder range encompassing all three domains. Collapsing is successful only in the supercritical domain, while in the remaining two domains it is possible only for the (L,R) parameters satisfying L|r|ν=const. The pertinent inset illustrates that the R-dependent position on the Sσ|r|-scale of the maximum of the scaled distribution Sτ′Dns(int)(S;R) attains its minimum value at R=Rc and saturates to 1 in the supercritical domain. (d) The main panel shows the collapsing (14) of the integrated size distributions DSall(int) of all avalanches for disorders from the subcritical domain, left inset shows the collapsing (19) of DSns(int) for the nonspanning avalanches, and right insets the collapsing (20) of various types of integrated size distributions DS1d(int), DS2d(int), and DS3d(int) for the spanning avalanches of type S1d, S2d, and S3d, respectively; the distributions correspond to the (L,R) pairs from the legend satisfying Lrν=const, and all collapses are obtained for τ′=τ+σβδ=2.03 and fractal dimension Df=2.98. Presented data are replotted from Reference [108], combining parts of Figs. 3, 10, 12 and 14

DS(int)(S;R)=S−τ′𝒟∞(Sσ|r|),(16)

for infinite systems, following from (15) in the L→∞ limit, and similarly for magnetization and susceptibility

MR(H)=|r|βℳ~∞(h^Reff/|r|βδ),(17)

χR(H)=|r|β−βδ𝒳~∞(h^Reff/|r|βδ),(18)

following from (12) and (13).

On the other hand, according to the same Reference [108], in the subcritical and transitional domains of disorder, the L-independent scaling (16) is violated, while the L-dependent scaling conjectures (14) and (15) remain in effect, see Fig. 3c. In the main part of Fig. 3d we illustrate the collapsing (14) in the subcritical disorder domain of the integrated size distribution of all avalanches DSall(int)=DSns(int)+∑αsp=13DSαsp(int), where DSns(int)(S;R,L) is the integrated size distribution of nonspanning avalanches, and DSαsp(int)(S;R,L) are the integrated size distributions of 1d, 2d, and 3d types of spanning avalanches (i.e., αsp=1,2,3). These distributions follow the scaling

DSns(int)(S;R,L)=L−Dfτ′𝒟~Sns(int)(SL−Df,|r|L1/ν),(19)

and

DSαsp(int)(S;R,L)=L−Dfτ′𝒟~Sαsp(int)(SL−Df,|r|L1/ν),(20)

predicting their collapsing when L1/ν|r|=const; here, Df=2.98 is the common value of fractal dimension providing the best distributions’ collapsing for all spanning and nonspanning avalanches [108], predicted in [48] to be Df=1/σν. Qualitatively similar behavior holds for the distributions of the remaining avalanche parameters (duration, energy and amplitude).

The integrated correlation functions, measuring the probability per spin that flipping of the spin initiating the avalanche will cause the flipping of spins at the distance x away within the same avalanche, scale as [28,31].

GR(int)(x)∼1xd+β/ν𝒢¯±(x|r|ν).(21)

Here d is the system dimension, and ν stands for the correlation length exponent describing the divergence of the correlation length ξ∼|r|−ν when r→0, presented in panel (d) of Fig. 4. In the supercritical disorder domain, correlation functions are monotonically decreasing with the inter-spin distance x in the cutoff region, while in the remaining domains, the onset of spanning avalanches causes the occurrence of characteristic bumps, as is displayed in the main panels (a)–(c) of Fig. 4. Scaling collapses of integrated correlation functions, satisfying the condition |r|L1/ν=const, are shown in pertinent insets.

Figure 4: Integrated correlation functions GR(int)(x) vs. inter-spin distance x shown in the main panels, and the pertinent collapses (21) shown in the insets in all three disorder domains, subcritical in (a), transitional in (b), and supercritical in (c). Divergence of the correlation length ξ when r→0, that is when R→(Rc)+, is displayed in panel (d) on the log-log scale in the main panel part, and on the lin-lin scale in the inset. Presented data are replotted from Reference [108], combining part of Fig. 15 and Fig. 17

4.2 Adiabatically Driven ZT NEQ RFIM on Two-Dimensional (2D) Lattices

The question of whether the adiabatically driven ZT NEQ RFIM on 2D quadratic lattices exhibits nontrivial critical behavior remained open for almost 20 years. This puzzle was initiated by the Mermin-Wagner theorem [174] which proved the absence of coexistence of two ferromagnetically ordered phases in isotropic Heisenberg models [175] in d≤2 dimensions. After that, the existence of two phases at low temperatures in 3D RFIM has been established in [92] using an exact renormalization-group (RG) flow down to the zero-temperature zero-field fixed point. Next, in 1989 the existence of an ordered phase at finite temperatures and weak fields for EQ RFIM in 3D was demonstrated, while in 2D it was rigorously proved that the ferromagnetic ordering is absent in the thermodynamic limit [91,176]. Additionally, many similarities between the EQ [86,177–179] and NEQ ZT RFIM [27,106] were observed regarding their criticality (matching of exponents, scaling functions, and spatial structures of avalanches), including the studies of four-dimensional systems [180] and at and beyond the upper critical dimension [181,182], leading to the conclusion that both models are rather likely to be in the same universality class for d≥3 dimensions [87,88]. Although not obtained for the NEQ model, these findings suggested that infinite 2D systems in the NEQ model may also lack ferromagnetic ordering because the two phases, evidenced in finite 2D systems, might vanish in the thermodynamic limit, i.e., when the system becomes big enough. The question of how large the NEQ system should be to lose/retain two ferromagnetically ordered phases was addressed in [183] where it was shown that the two phases ‘survive’ the thermodynamic limit provided they exist in 2D systems with size greater than the ‘breakup length’ Lb≈exp⁡(A/R2) with A=2.1±0.2 for the Gaussian distribution of the random magnetic field having disorder R.

Complementary to the preceding theoretical issues stood the question of whether the real 2D disordered ferromagnetic samples exhibit critical behavior until being experimentally confirmed in studies [184,185], opening a wide venue for further fundamental experimental and theoretical investigations of the 2D disordered ferromagnets and their applications.

4.2.1 The Case of Quadratic Lattices

Numerical simulations reported in [89,90], performed for the system sizes up to L=131072, significantly exceeding the ‘breakup length’ Lb and having up to ≈1.7×1010 spins, provided a numerical evidence for the critical behavior of the adiabatically driven ZT NEQ RFIM on the 2D quadratic lattices. The values of the non-universal critical parameters and of the universal critical exponents, obtained in [89,90], are quoted in Table 2 including the value of critical disorder Rc=0.54 corroborated in Fig. 5 implying Lb≈104.

Figure 5: Effective critical disorder Rceff(L) vs. system size L (symbols) and the power-law prediction [Rceff(L)−Rc]/Rceff(L)∼L−1/ν (continuous line) with Rc=0.54 and ν=5.14. Top inset shows Bray-Moore scaling ξ∼exp⁡[a~/[Rceff(L)−Rc]2], Rc∼0; see [186]. Bottom inset presents the windowed size distributions for R=0.38 to 0.55, L=65,536, and more than 600 RFC for each R. Their shapes in the bottom inset indicate the onset of spanning avalanches for R≤0.54. This is Fig. 6 from [89]

In Fig. 6a, we presented (the rising part of) the magnetization curves M(H) obtained in [89,90] for the system size L=131,072 and several values of disorder, ranging from R=0.52 (which is below the critical disorder Rc=0.54) up to R=0.76 surpassing the effective critical disorder Rceff(L) (= 0.605 for L=131,072). Due large value of L, one can take that the M(H) curves for disorder above the effective critical disorder Rceff(L) obey the scaling

MR,L(H)=|r|βℳ~(h^Reff/|r|βδ),(22)

which follows from (12) in the L→∞ limit, and analogously for the susceptibility curves

χR,L(H)=|r|β−βδ𝒳~(h^Reff/|r|βδ),(23)

and Eq. (13). The collapsing of M(H) curves according to (22) is presented in the main panel of Fig. 6b, and collapsing for the susceptibility curves and Eq. (23) in the inset. Besides, the distributions of avalanche parameters also indicate criticality by manifesting the power-law shapes. Thus, in Fig. 7 we show the distributions of size and duration in panels (a) and (b), respectively. In the bottom inset of panel (a) presented are the integrated size distributions for the five values of L and the corresponding effective critical disorder Rceff(L). The collapsing for R>Rceff(L) size distributions according to

Figure 6: (a) Rising part of the magnetization curves MR(H) for disorders R=0.52−0.76 and system size L=131072 where the effective critical disorder is Rceff=0.605. For R<Rceff, the value Hsp of the external magnetic field at which the spanning avalanche occurs varies with the random field configuration (RFC). Main panel shows magnetization curves for single RFC sorted in increasing disorder R. Inset displays magnetization curves averaged over 30 RFC; while this number is small, visible are the steps appearing due to spanning avalanches and stochastic nature of Hsp. (b) Collapsing (22) of four magnetization curves and (23) of susceptibility curves χ(H) are shown in the main panel and inset, respectively, for disorders R>Rceff from the legend. In this figure, combined are Figs. 1 and 2 from [89]

Figure 7: (a) The main panel shows the collapsing of the integrated distributions of avalanche size collected at R>Rceff for L=65,536, top insets show the collapsing of the corresponding windowed distributions, whereas the bottom inset shows the (non-scaled) size distributions collected for the values of L quoted in the legend at the corresponding effective critical disorder Reff(L). (b) Main panel and inset like for (a), but for the duration distribution. The number of RFC used in quenched averaging of the data ranged between 30 for the largest and 8000 for the smallest L. In this figure, combined are Fig. 3 from [89] and Figs. 6 and 7 from [90]

DS(int)(S;R)=S−(τ+σβδ)𝒟(Sσ|r|),(24)

following from Eq. (15) in the L→∞ limit, is presented in the main part of panel (a). Although the DS(int)(S;R) curves collapse well for large avalanches, they show a noticeable branching for small sizes, originating from the avalanches triggered near the end of the magnetization curve when almost all (except trapped) spins are already flipped, so that small avalanche sizes are more probable than Eq. (24) predicts. As the top inset shows, the branching disappears for the windowed distributions

DS(w)(S;R)=∫Hceff+rβδh0′Hceff−rβδh0′DS(S;R,H),(25)

collected in complementary windows (−rβδh0′,rβδh0′) covering the central parts of magnetization curves where the scaling

DS(S;R,H)=S−τ𝒟(Sσ|r|,h′|r|−βδ),(26)

of the size distribution DS(S;R,H) collected at the external field H applies, and h′=H−Hc−br; for details, see [90]. As is illustrated in panel (b), analogous behavior is found for the duration distribution, as well as for the (here not shown) distributions of avalanche energy and amplitude; see [90] for the latter two and also for the joint distributions for pairs of avalanche parameters. Values of the corresponding exponents are quoted in Table 3.

The behavior of the correlation function GR(int)(x) and the correlation length ξ(r,h′), depending on the reduced disorder r and reduced magnetic field h′, is presented in Fig. 8 for h′=0. The data from the main part and the bottom inset of panel (a) show the divergence ξ∼|r|ν of the correlation length ξ with the reduced disorder r, while the top inset shows the collapsing of the integrated correlation function GR(int)(x) according to the scaling prediction

GR(int)(x)∼x−(d+β/ν)𝒢±(int)(x|r|ν),(27)

where x is the distance between the spins flipped in the same avalanche. Complementary, the panel (b) shows the collapsing of the correlation function GR,H(x), giving correlations at the external field H, which scales as

GR,H(x)=1d−2+η𝒢±(xξ(r,h′)),(28)

where η is the exponent named anomalous dimension [28,31] (here, η=1).

Figure 8: (a) The main part shows power-law divergence ξ∼|r|ν of the correlation length ξ with reduced disorder r for reduced magnetic field h′=0. In the bottom inset, the same data are shown against the disorder R on a linear scale. Top inset: scaling collapse (27) of the correlation function GR(int)(x) for disorders R=0.64−0.90 and system size L=131,072. The curves are averages of 30 RFCs for each R. (b) The scaling collapse (28) of the correlation function GR,H(x) for the avalanches at reduced field h′=0. The collapse is obtained for η=1 and correlation lengths from the left panel. Inset: the same collapse on the lin−log scale illustrates the applicability of the approximation GR,H(x)∼exp⁡(x/ξ)/xd−2+η which is used in determination of the correlation length ξ. In this figure are combined Figs. 4 and 5 from [89]

4.2.2 Spanning Avalanches on Quadratic Lattices

Below the effective critical disorder Rceff(L), system response is dominated by the spanning avalanches studied in [107] in the case of adiabatically driven ZT NEQ RFIM on 2D equilateral quadratic lattices. The number per single run, N2(R,L), of 2d spanning avalanches (i.e., the avalanches that span the system along both dimensions), shown in Fig. 9a against R for lattice sizes L=1024−16,384, scale as

N2(R,L)=N~2(rL1/ν),(29)

enabling their collapsing presented in the panel (b), and indicating that limL→∞N2(R,L)=U(Rc−R), where U(x) is the unit step function. On the other hand, the number per single run N1(R,L) of 1d spanning avalanches (spanning the system along only one dimension and shown in the inset of left panel) scales as

N1(R,L)=L−θ1N~1(rL1/ν),(30)

with θ1=0.08±0.02, meaning that they become irrelevant in the L→∞ limit. The number per single run Nns(R,L) of the remaining (i.e., nonspanning) avalanches scales as

Nns(R,L)=L2N~ns(rL1/ν),(31)

and is presented in the inset of Fig. 9b. Due to scale invariance, the clusters of spins flipped during a spanning avalanche are fractals; their fractal dimension Df is smaller for 1d than for the 2d spanning avalanches as is illustrated in Fig. 10.

Figure 9: (a) Main part and inset show the number of spanning avalanches per single run, N2(R,L) for 2d and N1(R,L) for 1d spanning avalanches, respectively, against disorder R for the system sizes L quoted in legend. (b) Collapsing (29) of 2d spanning and (31) of nonspanning avalanches shown in the main part and inset, respectively. In this figure, combined are Figs. 2, 4, 5 and 13 from [107]

Figure 10: Plots of a 2d spanning avalanche (left) and a 1d spanning avalanche (right) for L=4096 and R=0.68 [107]; the time scale of spin flipping is shown by the color legend; not affected spins are white. Fractal dimensions are: Df=1.9828 for the 2d spanning avalanche, and Df=1.9125 for the 1d spanning avalanche. Presented figure is replotted from Figs. 7 and 8 from Reference [107]

Spanning avalanches cause bumps in distributions and jumps in magnetization ΔM realized at some value of the external field Hsp, called the spanning field, whose distribution is determined by the RFC. This is shown in [107] and illustrated in Fig. 11.

Figure 11: (a) (Scaled) size distributions of all avalanches (circles), nonspanning avalanches (full line), and spanning avalanches in the inset. (b) Magnetization jumps ΔM shown in inset, and their collapsing in the main panel. (c)–(e) Distribution of spanning field Hsp for three different values of disorder R; for details, see [107]. In this figure, combined are Figs. 14, 18 and 19 from [107]

4.2.3 The Case of Triangular and Hexagonal Lattices

Besides quadratic, the ZT NEQ RFIM can be studied on other 2D lattices with translation symmetry, but with different topology of elementary cell, e.g., with different number of nearest neighbors given by the coordination number znn. Two examples are the triangular and hexagonal lattices which, compared to quadratic lattice with four nearest neighbors, have six and three nearest neighbors for each spin, respectively; see Fig. 12.

Figure 12: Triangular (left) and hexagonal (right) lattices. The right panel is Fig. 1 from [103]

The studies [100] on triangular and [103] on hexagonal lattices gave partial answers to the conjecture introduced in [135] that the coordination number znn, and not solely the lattice dimension (as believed beforehand), plays the key role in determining the universality class of the ZT NEQ RFIM critical behavior. Thus, the study [100] showed that the ZT NEQ RFIM exhibit critical behavior on the triangular lattices (where znn=6) in accordance with the general expectation that this should be the case for the lattices with znn≥4. However, the values of some critical exponents and nonuniversal critical parameters reported for this lattice, see in Table 4, differ outside the error bars from the corresponding ones reported in [89,90] for the quadratic lattices, suggesting that the universality classes for the triangular and quadratic lattices might be different. As an illustration of the critical behavior of the model on triangular lattice, in Fig. 13, we show the graph of Rceff(L) vs. L, implying that Rc=0.86 for the triangular lattice (left panel), and the collapsing of the windowed size distributions (right panel) obtained for the exponents from Table 4. Similar study [103] revealed the absence of critical behavior on the hexagonal lattice having znn=3.

Figure 13: (a) Effective critical disorder Rceff(L) for the triangular lattice vs. system size L (symbols), and the power-law prediction [1−Rc/Rceff(L)]∼L1/ν (straight line) for Rc=0.86 and ν=5.26. (b) The scaling collapse of the integrated distributions of avalanche size for triangular lattice with L=65,536. In this figure are combined Figs. 2 and 8 from [100]

4.3 Crossover from 3D to 2D ZT NEQ RFIM Systems

The dimensional crossover from three to two spatial dimensions has been studied so far both experimentally [22,188] and in equilibrium models [189–192], where it has been established that the systems with constant thickness l, and two diverging spatial dimensions L→∞, behave in the asymptotic limit essentially as 2D systems, showing a critical temperature Tc(l) that shifts from the critical temperature of the planar 2D system, Tc(l=1)=Tc2D, to the critical temperature of the bulk 3D system, Tc(l→∞)=Tc3D, with a first approximation for this crossover function Tc(l)−Tc3D∼l−1/ν3D for large l, where ν3D is the correlation length exponent in the 3D case.

The 3D to 2D dimensional crossover in ZT NEQ RFIM was studied in [110–112] on nonequilateral 3D L×L×l cubic lattices with quadratic L×L base and thickness l in adiabatic regime with closed boundary conditions on base, and open conditions on thickness. For each l, one can define the critical disorder Rc(l) for infinite (i.e., L→∞) systems as the value of disorder such that for R<Rc(l) there is a finite jump ΔM in the magnetization curve M(H) tending to zero when R tends to Rc(l). Additionally, for R>Rc(l) the curve M(H) is smooth and has finite susceptibility χ(H)=dM/dH for any H, while at R=Rc(l) the M(H) curve is still smooth but with infinite susceptibility at some value Hc(l) of the external magnetic field, called the critical field for the infinite lattices of thickness l. In finite L×L×l systems the biggest change in magnetization, and therefore the maximum value of susceptibility, is caused by the spanning avalanches which (roughly speaking) appear at disorders R≤Rceff(l,L), where Rceff(l,L) is the effective critical disorder (for given thickness l and base size L) which tends to Rc(l) in the L→∞ limit. At R=Rceff(l,L) the susceptibility attains its maximum value at the value Hceff(l,L) of the external field called the effective critical field which tends to Hc(l) in the L→∞ limit. The analyses in [110] resulted in the analytical prediction

Rcth(l,L)=Rc3D[1−Δl1/ν3D−(A−Δ)l1/ν2DL1/ν2Dl1/ν3D]−1,(32)

for the effective critical disorder Rceff(l,L); here A=0.63±0.18 is an adjustable parameter, Δ=1−Rc3D/Rc2D, while ν3D and Rc3D are the correlation length exponent and critical disorder in 3D model (and analogously for ν2D and Rc2D and 2D model), see Eq. (7) in [110]. Extension of this analysis, presented in [112], led to

Hcth(l,L)=Hc3D+Θcrossl1/ν3D+B−Θcrossl1/ν3D(lL)1/ν2D(33)

predicting the effective magnetic field Hceff(l,L), where B=0.20±0.07, Θcross=0.68±0.07, and Hc3D is the critical field in 3D model. From (32) and (33) follow the expressions

Rcth(l)=Rc3D[1−Δl1/ν3D],(34)

Hcth(l)=Hc3D+Θcrossl1/ν3D,(35)