Study of Quasi-Static Magnetization with the Random-Field Ising Model

: The topic of this paper is modeling based on Hamiltonian spin interactions. Preliminary studies on the identiﬁcation of quasi-static magnetizing ﬁeld in a magnetic system were presented. The random-ﬁeld Ising model was then used to simulate the simpliﬁed ferromagnetic structure. The validation of algorithms and simulation tests were carried out for the 2D and the 3D model spaces containing at least 10 6 unit cells. The research showed that the response of a slowly driven magnetic system did not depend on the external ﬁeld sweep rate. Changes in the spatial magnetization of the lattice were very similar below a certain rate of the external ﬁeld change known as the quasi-static boundary. The observed di ﬀ erences in obtained magnetization curves under quasi-static conditions stemmed from the random nature of the molecular ﬁeld and the avalanche-like magnetization process


Introduction
The concept of a cellular automaton was introduced by Ulam and Von Neumann [1,2]. At present, cellular automata are used as a scientific tool in studies of complex dynamic systems, microstructures, and many other applications [3][4][5][6][7]. Cellular automaton (CA) is an example of how simple rules and local interactions can lead to very diverse and complicated behaviors [8,9]. Usually, cellular automata are simple and regular grids. Each node of the grid is in a defined discrete state. The states of the nodes are updated synchronously at discrete moments and the state of each node at the next moment is a function of the state of its neighbors in the current moment. The Ising model, being a magnetic system based on statistical mechanics, is a magnificent example of a cellular automaton that was developed much earlier than the aforementioned cellular automaton concept [10,11]. The model is widely used in computer physics for the analysis of magnetic hysteresis [12,13] phase transitions [14,15], dynamics of magnetization [16,17], Barkhausen noise [13,15,17,18], and energy loss in magnetic materials [19], as well as for defect and inhomogeneity studies in magnetic structures [20][21][22][23][24][25].
Magnetic modeling using the interactions of so-called "spins" arranged in a lattice of regular unit cells is the most widely used approach. A single spin can be interpreted as a magnetic moment of an atom in a crystal lattice, crystal grain, magnetic domain, or quantized amorphous space [26]. Modeling of spontaneous magnetization processes including hysteresis, displacements of domain walls, and the Barkhausen effect are some of the common applications of the Ising model [26][27][28]. The vast majority of works use the random-field Ising model (RFIM) based on the Hamiltonian as the energy function with an additional term representing the random molecular magnetic field assigned to each spin [13,14,29,30].
However, the RFIM does not explicitly take into account the time scale. Therefore, the simulation results are qualitatively and quantitatively consistent only for time-independent processes such as the primary magnetization curve or the static magnetic hysteresis loop [28].
Studies of hysteresis in soft magnetic materials have always been a part of scientific work. Various magnetic, physical, structural, thermodynamic, and many other parameters are identified on the basis of the shape and area of the hysteresis loop. Basically, loops can be classified as static or dynamic [12,31], but the magnetizing field must be variable over time to obtain both hysteresis loops. The static loop is interesting because it illustrates the magnetizing process without time-dependent effects. In practice, the change in material magnetization is quite complex and there is no strict interpretation when the process can be considered static, even when an extremely low frequency sweep rate is used [6,32]. Therefore, this approach shall only be considered as a quasi-static approximation of a boundary static process [32]. Many works on testing magnetic materials arbitrarily assume, without any explanation, that the frequency sweep rate is low enough that, e.g., eddy currents may be neglected. The real limitation of the static condition results very often from the technical limitations of the measuring instruments. In this case, one cannot assess how large the discrepancy of such an approximation is. In this context, the paper provides a contribution in the field of slowly driven magnetic systems through the identification of the quasi-static limiting of magnetization processes in magnetic materials. Coliaori et al. [19] showed the dependency of field sweep rate in the RFIM on the coercive field of magnetic thin films and pointed out a direct link between the theory of loss separation for bulk materials and a simple magnetic system with disorder. Thus, the aim of this work is to check if RFIM, as a mathematical model of ferromagnetism, shows a tendency to reach the limit of quasi-static magnetization under general conditions.

Random-Field Ising Model
The model space is defined as a set of unit cells associated with the S zi spin operators. Spins denoted as S zi are surrounded by j = 2 D adjacent spins S zj . The S zi spin is coupled with S zj through exchange interactions. Spin operators take values from the set of S zi , S zj = { S − z = −1, S + z = +1} and change the values to the opposite by spin flip.
A set of periodically arranged unit cells creates either a spatial or a flat model lattice with spins placed in all sites of the lattice. The number of spins in the model space depends on the form of the unit cell, the number of space dimensions D and the lattice length L for each dimension. The unit cell is here referred to as a site in the lattice. The model space is therefore defined as the D-dimensional lattice of regular unit cells limited in each direction by its pertinent length L D , and comprising N sites with localized spins (Equation (1)).
Other magnetic structures require consideration of both other forms of unit cells and other numbers of localized spins coupled by exchange interactions [20,27]. The use of regular cells instead of more complicated forms is a typical simplification of the model space. In other cases, when the structure is not a crystal lattice or has a very heterogeneous form, one can quantize the space with regular cells and estimate equivalent parameters. The fundamental Hamiltonian form H of the Ising model takes into account both the exchange interactions between the nearest-neighbored spins in the lattice and the applied external field H. The Hamiltonian of the RFIM (Equation (2)) imposes quenched disorder in the form of a random molecular field h i assigned to all sites i of the lattice [26].
The value of the h i field in Equation (2) has an impact only on triggering the spin S zi flip induced either by adjacent spins S zj or by the external field H [28]. Parallel polarization of all spins in the lattice occurs under the influence of exchange interactions in the Weiss molecular field [12,13]. In the RFIM, the subsequent values of the h i field are random variables and depend on the parameters of the random distribution described by Equation (3).
The set of random values h i ∈ {h 0 , h 1 , . . . h N-1 } is generated on the basis of the Gaussian distribution function described by the relationship in Equation (3) and is additionally subject to the limitations defined by Equation (4). The mean value of the set h i in the space defined by N spins tends to zero (Equation (4)). Variance of random fields h i is a function of the disorder R at the temperature T~0 K.
The range of changes in the parameter R is limited by the value of critical disorder R c [13,28]. The sequence of random h i ∈ {h 0 , h 1 , . . . h N-1 } for R > R c may contain h i elements that do not guarantee the convergence of the simulation process (non-ergodic sequence of random variables) [29]. A more comprehensive analysis of the disorder parameter in the RFIMs was carried out in works [15,17,18,26,28]. If the effective field h eff (Equation (5)) acting in a given site i changes the sign to the opposite, then flipping of the S zi spin occurs and its value is also changed to the opposite.
The flip condition (Equation (6)) may occur if there is either a change in the configuration of adjacent spins S zj for a given instantaneous value of field H, or a change in the instant value of the magnetizing field H with a fixed configuration of spins S zj .
The initial equilibrium of the model is achieved at a temperature close to absolute zero (with no thermal fluctuations, accidental spin excitations, changes in spontaneous magnetization saturation): M s (T = 0, H = 0). Evolution from the initial equilibrium state to the subsequent equilibrium or metastable states takes place in the magnetization process, forced by a homogeneous, external H field. The estimation of average magnetization during metastable magnetization states (H, m = const) is based on the relationship in Equation (7), Where S zi + denotes spins that have been flipped and S zi spins that remain in their initial configuration. To estimate the area of the model space where unstable changes in magnetization are induced, the correlation function (Equation (8)) of x number of spin flips denoted as x and random field disorder R was used [13,26,28].
Algorithms 2020, 13, 134 4 of 12 The estimation of the G value for a given R indicates the dynamics of magnetization and shows changes in the number of correlated spin flips (avalanches) triggered by a single spin flip.

Numerical Implementation
Numerical implementation of the RFIM involves the algorithms for generation of pseudo-random numbers and collective spin flipping.
To determine the sequence of random fields h i , the three-stage random number generator was used. In the first step, a start sequence of five pseudo-random values initiated the relevant generator. The initial set was generated by a linear congruential generator for a given seed value [33]. Next, a sequence of pseudo-random numbers from the range (0,1) was generated on the basis of the Tausworth algorithm [34]. In the third stage, one random variable with a Gaussian distribution was calculated from the two previously generated independent variables of the even distribution by means of the Box-Muller method [33,35]. The output value was scaled according to the determined standard deviation specified by the disorder R. All S zi spins were set to −1 and random values h i are assigned to them. The initial value of the external field H min < min(h i ) was determined.
Avalanche spin flipping was triggered by increasing gradually the external field with a step ∆H. The condition (Equation (5)) of spin flip was checked by the brute force linear search method [18]. The algorithm is not very effective numerically, but allows one to determine explicitly the instantaneous changes in spatial magnetization and map other parameters of the model. The space was searched as long as all spins for the given external field H could be flipped. In the next step, the external field H i+1 = H i + ∆H was again increased and the subsequent step of the search algorithm was executed. The procedure was repeated until all spins in the entire model space were flipped (Equation (6)).

Results and Discussion
Numerical simulations using RFIM were carried out to verify the Gaussian pseudo-random number generator and spin flip procedure. Subsequently, the dynamics of the magnetization process of 2D and 3D structures was examined.
Model verification included both the discrepancies in the average value of the generated pseudo-random set h i according to dependencies 3 and 4 and the correlation of the spin flip process G (x, R) depending on the set of pseudo-random parameter values using the correlation functions (Equations (8) and (9)).
Analyzing the changes in the average value of the h i parameters against the size of the model space N for various values of the disorder R (Figure 1), one can assume that these changes are negligibly small in comparison with the h i values for the model space comprising above N = 10 4 spins. However, if the value N falls below 10 4 , one should either consider the non-equilibrium magnetization state resulting from the non-zero mean h i or reduce the discussed problem by increasing the size of the model space so that the number of spins would exceed 10 4 . Thus, the reduced impact of the non-equilibrium initial state occurs in simulations of two-dimensional spaces (D = 2) for the values of the parameter L > 100, and L > 22 in the case of the three-dimensional model (D = 3).
Another tested parameter of the model was the correlation of spin flips in areas of the model where avalanche processes occurred. The effect of avalanches, which introduce discontinuous changes in the magnetization of the model space, was determined on the basis of Equation (9). It can therefore be assumed that the discrete states of magnetization, recorded in the course of simulation, also strongly depend on the correlation of spin flips. The value of spin flip correlation can be estimated by determining the distance x (avalanche size) defined as the number of sites of inverted spins relative to the spin starting the avalanche process. Based on the normalized number of classes, range width, and class limits, the histogram of correlated inversions was determined in accordance with the relationships found in Equations (8)  Analyzing the changes in the average value of the hi parameters against the size of the model space N for various values of the disorder R (Figure 1), one can assume that these changes are negligibly small in comparison with the hi values for the model space comprising above N = 10 4 spins. However, if the value N falls below 10 4 , one should either consider the non-equilibrium magnetization state resulting from the non-zero mean hi or reduce the discussed problem by increasing the size of the model space so that the number of spins would exceed 10 4 . Thus, the reduced impact of the non-equilibrium initial state occurs in simulations of two-dimensional spaces (D = 2) for the values of the parameter L > 100, and L > 22 in the case of the three-dimensional model (D = 3).  Another tested parameter of the model was the correlation of spin flips in areas of the model where avalanche processes occurred. The effect of avalanches, which introduce discontinuous changes in the magnetization of the model space, was determined on the basis of Equation (9). It can therefore be assumed that the discrete states of magnetization, recorded in the course of simulation, also strongly depend on the correlation of spin flips. The value of spin flip correlation can be estimated by determining the distance x (avalanche size) defined as the number of sites of inverted spins relative to the spin starting the avalanche process. Based on the normalized number of classes, range width, and class limits, the histogram of correlated inversions was determined in accordance with the relationships found in Equations (8) and 9). A graph of the spin flip correlation function was created by accumulating number of flipped spins in class intervals for a given distance x of correlated spin interactions ( Figure 2) The high variability of the spin flip correlation function indicated that avalanches were nucleated frequently [13,30], but their propagation was often quenched by fluctuations of random fields hi. The size of avalanches depended on the value of disorder introduced to the model space. Test results of the developed pseudo-random number generator (PRNG) showed that the numerical performance of the generator did not introduce significant limitations associated with simulation time. The implementation efficiency of the PRNG algorithm at 10 7 numbers per second was sufficient for the RFIM model, because the simulation did not require the generation of pseudorandom numbers ad hoc. Considering the limitations related to the mean value of the hi set, it should be stated that the model space with 10 6 spins guarantees the correct convergence of computer simulations [26,27].
The developed algorithm of a three-stage pseudo-random generator (Section 3) led to the formation of clusters. These were areas where the directly coupled spins had random values of hi and the same sign. Clusters can introduce a weak effect of collective interactions. It is very likely that all spins in a cluster will be flipped in the same simulation step. Clusters are able to boost occurring avalanches or trigger the new ones. Figure 3 shows the distribution of random fields with the same sign over the 2D model space. Additionally, the exemplary clusters are highlighted. Test results of the developed pseudo-random number generator (PRNG) showed that the numerical performance of the generator did not introduce significant limitations associated with simulation time. The implementation efficiency of the PRNG algorithm at 10 7 numbers per second was sufficient for the RFIM model, because the simulation did not require the generation of pseudo-random numbers ad hoc. Considering the limitations related to the mean value of the h i set, it should be stated that the model space with 10 6 spins guarantees the correct convergence of computer simulations [26,27].
The developed algorithm of a three-stage pseudo-random generator (Section 3) led to the formation of clusters. These were areas where the directly coupled spins had random values of h i and the same sign. Clusters can introduce a weak effect of collective interactions. It is very likely that all spins in a cluster will be flipped in the same simulation step. Clusters are able to boost occurring avalanches or trigger the new ones. Figure 3 shows the distribution of random fields with the same sign over the 2D model space. Additionally, the exemplary clusters are highlighted. The impact of both the disorder and the range of variation of the disorder parameter R on magnetization were examined. The magnetizing process depending on the disorder parameter R was depicted in graphs of meta-stable magnetization changes dM/dH (Figures 4 and 5) and hysteresis loops ( Figure 6).  The impact of both the disorder and the range of variation of the disorder parameter R on magnetization were examined. The magnetizing process depending on the disorder parameter R was depicted in graphs of meta-stable magnetization changes dM/dH (Figures 4 and 5) and hysteresis loops ( Figure 6). The impact of both the disorder and the range of variation of the disorder parameter R on magnetization were examined. The magnetizing process depending on the disorder parameter R was depicted in graphs of meta-stable magnetization changes dM/dH (Figures 4 and 5) and hysteresis loops ( Figure 6).  . A greater magnetic order of the magnetic system and thus a smaller standard deviation of the random field component in each site of the model space affected the exchange energy of the localized spins in the considered site. The ratio of the strength of exchange interactions to the random field hi was greater, and therefore the energy needed to overcome the magnetic order in each site was also greater. Therefore, the coercive field was the highest when the influence of the random field component in the model space was the lowest. A magnetic system with higher disorder tends to lower magnetic susceptibility and reduces the coercive field as well [19]. Hence, smooth hysteresis . A greater magnetic order of the magnetic system and thus a smaller standard deviation of the random field component in each site of the model space affected the exchange energy of the localized spins in the considered site. The ratio of the strength of exchange interactions to the random field hi was greater, and therefore the energy needed to overcome the magnetic order in each site was also greater. Therefore, the coercive field was the highest when the influence of the random field component in the model space was the lowest. A magnetic system with higher disorder tends to lower magnetic susceptibility and reduces the coercive field as well [19]. Hence, smooth hysteresis A greater magnetic order of the magnetic system and thus a smaller standard deviation of the random field component in each site of the model space affected the exchange energy of the localized spins in the considered site. The ratio of the strength of exchange interactions to the random field h i was greater, and therefore the energy needed to overcome the magnetic order in each site was also greater. Therefore, the coercive field was the highest when the influence of the random field component in the model space was the lowest. A magnetic system with higher disorder tends to lower magnetic susceptibility and reduces the coercive field as well [19]. Hence, smooth hysteresis loops significantly depend on the symmetrical shape of the h i distribution rather than the nearest neighbor interactions.
Figures 7a-d and 8a-d respectively present mapping of the 2D and the 3D metastable states of magnetization due to spin flipping. Single-colored areas (the colors are attributed randomly) refer to regions in which the magnetization processes were avalanche-like. Flip of spins in these areas occurred under conditions of the temporarily determined value of the external field H. Each single spin flip for a given field value was marked with the same color of pixels/cubes on the 2D and 3D maps.
Algorithms 2020, 13, x FOR PEER REVIEW 8 of 12 loops significantly depend on the symmetrical shape of the hi distribution rather than the nearest neighbor interactions. Figures 7a-d and 8a-d respectively present mapping of the 2D and the 3D metastable states of magnetization due to spin flipping. Single-colored areas (the colors are attributed randomly) refer to regions in which the magnetization processes were avalanche-like. Flip of spins in these areas occurred under conditions of the temporarily determined value of the external field H. Each single spin flip for a given field value was marked with the same color of pixels/cubes on the 2D and 3D maps.  Simulations of the quasi-static process were performed in the same way as earlier tests of numerical implementation efficiency of the RFIM model. A single-spin flip method and a sorting algorithm were used to identify spins meeting the conditions described by the relationships 5 and 6.
The quasi-static magnetization process, almost independent of the external field step rate, is depicted in Figure 9. Hence, it is possible to achieve a similar course of metastable changes in magnetization triggered by an external field with a different rate of change. Similar shapes of avalanches presented in Figure 7c,d also may prove the quasi-static magnetization process. The Simulations of the quasi-static process were performed in the same way as earlier tests of numerical implementation efficiency of the RFIM model. A single-spin flip method and a sorting algorithm were used to identify spins meeting the conditions described by the relationships 5 and 6.
The quasi-static magnetization process, almost independent of the external field step rate, is depicted in Figure 9. Hence, it is possible to achieve a similar course of metastable changes in magnetization triggered by an external field with a different rate of change. Similar shapes of avalanches presented in Figure 7c,d also may prove the quasi-static magnetization process. The observed damping of the metastable changes in magnetization was generally consistent with experimental data [12] (Vol.II, Ch.3). The small increment of the external field H had a negligible impact on the occurring avalanche-like changes in magnetization. The avalanche of spin flips caused discrete changes in magnetization called the Barkhausen effect [15,27]. The greatest impact of the spin-flip avalanche on the change of magnetization was observed in simulations with moderate disorder R∈<1.0, 2.0> of random field h i .
Algorithms 2020, 13, x FOR PEER REVIEW 10 of 12 observed damping of the metastable changes in magnetization was generally consistent with experimental data [12: Vol.II,Ch.3]. The small increment of the external field H had a negligible impact on the occurring avalanche-like changes in magnetization. The avalanche of spin flips caused discrete changes in magnetization called the Barkhausen effect [15,27]. The greatest impact of the spin-flip avalanche on the change of magnetization was observed in simulations with moderate disorder R∈<1.0, 2.0> of random field hi.

Conclusions
The RFIM model was applied for the numerical studies on slowly driven magnetic systems. Preliminary research on the identification of a quasi-static magnetizing processes in 2D magnetic systems was presented. The paper proves that the a magnetic system driven with a field sweep rate below the quasi-static limit reaches convergent response. Hence, if the response of the system does not depend on the magnetizing field sweep rate, the shape of the magnetizing waveform should not matter. Further development of the model should be focused on the implementation of an explicit time scale and implementation of the long-range magnetic interactions. The quasi-static boundary may be quantified in the expanded model with time-dependent spin interactions and a time-varying external magnetizing field.

Conclusions
The RFIM model was applied for the numerical studies on slowly driven magnetic systems. Preliminary research on the identification of a quasi-static magnetizing processes in 2D magnetic systems was presented. The paper proves that the a magnetic system driven with a field sweep rate below the quasi-static limit reaches convergent response. Hence, if the response of the system does not depend on the magnetizing field sweep rate, the shape of the magnetizing waveform should not matter. Further development of the model should be focused on the implementation of an explicit time scale and implementation of the long-range magnetic interactions. The quasi-static boundary may be quantified in the expanded model with time-dependent spin interactions and a time-varying external magnetizing field.