Influence of initial states and structure defects on non-equilibrium critical behavior of 3D and 2D Ising models

Results of Monte Carlo description of features of non-equilibrium critical behavior in three- and two-dimensional Ising models conditioned by presence of structural quenched disorder are presented. It was revealed that the increase in the concentration of defects is accompanied by the strengthening of the aging effects manifested as the slowing down of correlation and relaxation processes in structurally disordered systems as compared to the pure system. Non-equilibrium initial states begin to increasingly more strongly influence peculiar features and characteristics of the system’s evolution. It was shown that limiting fluctuation-dissipation ratio values X∞ characterizing the degree of system departure from equilibrium satisfy inequality X∞ < 1 and depend on the universality class of non-equilibrium critical behavior to which they belong: one of these classes corresponds to the high-temperature, and the other to the low-temperature initial state of the system. Values of these non-equilibrium characteristics in 3D Ising model demonstrate belonging to universality classes of critical behavior of weakly and strongly diluted systems, but non-equilibrium characteristics in 2D Ising model depend continuously on concentration of defects owing to the crossover effects in the percolation behavior.


Basic concepts and model representations of the non-equilibrium behavior
One of the peculiarities arisen in describing the critical behavior of systems is the critical slowing down effect. It is associated with an anomalous increase in the system relaxation time t rel when approaching the second-order phase transition temperature T c . As a result, the system being at the critical point appears in no condition to reach equilibrium throughout the entire relaxation process. Therefore, at times t t rel , extraordinary non-equilibrium phenomena characteristic of systems with slow dynamics, such as aging, fluctuation-dissipation theorem (FDT) violation, and the effect of various initial non-equilibrium system states arise in the behavior of systems [1].
The aging effects appearing in the non-equilibrium relaxation stage of the system with slow dynamics are characterized by the existence of two-time dependencies for such functions as the autocorrelation function C(t, t w ) and the response functions R(t, t w ) on the waiting time t w and observation time t − t w determined for the spin system with spin density S(x) by the relations International where ... characterizes the statistical averaging over different realizations of initial spin configurations, h is a weak magnetic field applied to a system at the time point t w . The waiting time t w is defined by time between a sample preparation and the beginning of measurement of its characteristics. During t, t w t rel , the effect of initial states of the system and aging effects manifests itself in the temporal system behavior, which is characterized by both translation symmetry breaking in time and slowing down relaxation and correlation processes with increasing the system "age" t w .
It is supposed that the non-equilibrium behavior of a system is realized via its transition at the starting instant t 0 from the initial state at temperature T 0 to the state with temperature T s differing from T 0 . The accompanying equilibration process is characterized by relaxation time t rel (T s ), and equilibrium corresponding to temperature T s is reached in times t t rel (T s ), while the system dynamics prove stationary and invariant with respect to time reversal. However, in times t t rel (T s ), the evolution of the system depends on its initial state. In this connection, the non-equilibrium behavior of the system depends on whether it evolves from a high-temperature T 0 > T s or a low-temperature T 0 < T s initial state.
According to general concepts about non-equilibrium processes, it is expected that for times t > tw t rel (T s ) C(t, t w ) = C eq (t − t w ) and R(t, t w ) = R eq (t − t w ), where C eq and R eq are corresponding equilibrium values related by the fluctuation-dissipation theorem (FDT) The consequences of the FDT provide theoretical grounds of various experimental methods for the measurement of radiation scattering and absorption by matter.
In the case of non-equilibrium behavior of systems for t, t w t rel (T s ), the generalized FDT assumes the form with introduction of the quantity X(t, t w ) which is called as the fluctuation-dissipation ratio (FDR): For times with t > t w t rel , the FDT establishes that X(t, t w ) = 1. However, in the general case, X(t, t w ) = 1 for times with t, t w t rel . The asymptotic FDR value introduced as is an important universal characteristic of non-equilibrium processes in the various systems. Let us consider the general properties of X ∞ and its dependence on the system's quenching temperature T s . For systems demonstrating dissipative slow dynamics from the hightemperature initial state with initial magnetization m 0 = 0 (m 0 1), the free propagator R 0 (q, t, t w ) and correlator C 0 (q, t, t w ) as the bare response function and correlation function are characterized by following expressions [1,2]: where C is the correlator characterizing the influence of the starting conditions. In (6)(7)(8)(9), parameters q and τ ∼ (T − T c )/T c are wave vectors and the reduced temperature, correspondingly. The FDR in momentum space is given by the expression and we obtain in the Gaussian approximation The asymptotic FDR in the Gaussian approximation is characterized by the following values This means that only the zero order-parameter mode with q = 0 is characterized by aging effects at the critical point; in other words, it does not relax to equilibrium state and the FDT for this mode is violated.
Non-equilibrium dynamics of systems evolving from a low-temperature initial state with initial magnetization m 0 = 0 is characterized by following expressions for the bare response function and correlation function [1,3]: where t m = B m m −k 0 is a new time-scale with k = 2 and B m = 2/3 in the Gaussian approximation.
At the critical point τ = 0 the FDR can be presented in the Gaussian approximation as with x = t w /t m = 2t w m 2 0 /3. The FDR in (15) is unrelated to t, i.e., X 0 (q, t, t w ) ≡ X 0 (t w /t m , q 2 t w ), which is a distinctive feature of the FDR in the Gaussian approximation. X 0 (x = 0, y) corresponds to case with evolution from the high-temperature initial state.
For the zero mode with q = 0 (y = 0), the FDR as a function of variable x is characterized by the relation X 0 (x, y = 0) = 4 5 1 + 3 5 Its value grows monotonically from X 0 (x = 0, y = 0) = 1/2 at x = 0 to X 0 (x → ∞, y = 0) = 4/5 as x → ∞. Thus, the FDR for the global order parameter (complete magnetization) in the longtime limit (x → ∞) reaches the value of X ∞ 0 = 4/5 if initial magnetization m 0 = 0 in the case of evolution from a magnetized low-temperature state, and X ∞ 0 = 1/2 if initial magnetization m 0 = 0 in the case of evolution from a high-temperature state. The value of X ∞ 0 = 4/5 invariably remains independent of a concrete m 0 = 0 value.
Strong fluctuation effects accompanying second-order phase transitions result in fluctuation corrections to values of these characteristics of non-equilibrium critical behavior considered in the Gaussian approximation.
Consequently, we can generalize that for states of the system with temperature T s > T c , X ∞ (T s > T c ) = 1. For a low-temperature ordered phase with T s < T c , X ∞ (T s < T c ) = 0. These results are unrelated to specific properties of individual systems. In the case of T s = T c , there are no general arguments defining the value of X ∞ (T s = T c ), which necessitates its derivation for each individual statistical model, and X ∞ (T s = T c ) is a universal quantity associated with the universality class of the model's critical dynamics.

Influence of defects on critical behavior
The renormalization-group [4,5], numerical [6][7][8][9], and experimental [10] investigations of the critical dynamics of structurally disordered systems have clearly established that the presence of uncorrelated structural defects and defects with effects of long-range correlations in the studied systems leads to new types of critical behavior and to a significant enhancement of critical slowing down effects in comparison with pure systems. In this respect, specific features of a non-equilibrium behavior, such as aging effects, should undoubtedly find a more striking manifestation in structurally disordered systems with new universal values of the fluctuationdissipation ratio.
Well-known that phenomenological approach to description of an universal critical behavior of systems with disorder is determined by the Hamiltonian where φ(x) is the n-component order parameter, V (x) characterizes a potential of random field of defects. Under the assumption that spatial distribution of the potential V (x) is specified by a Gaussian distribution, statistical properties of V (x) are determined by only first moment with V (x) = 0 and second moment V (x)V (y) ∼ u(x − y). In the case of uncorrelated point-like defects [11], where v ∼ V 2 c imp , c imp is a concentration of defects. The model with ε D -extended defects parallel oriented relative to one another and random distributed in sample perpendicularly to defect direction [12] is characterized by The model of long-range correlated defects with random orientations was proposed in the work of Weinrib and Halperin [13] is determine by The free energy of disordered system is given by relation where . . . indicates an averaging over disorder configurations, . . . 0 means statistical averaging with the Hamiltonian of pure system. The specific heat of pure system can be introduced as where α 0 is critical exponent of the specific heat for pure system. With the use of scaling analysis, we can determine the free energy F 0 ∼ τ 2−α 0 and the correlator . Then asymptotic behavior of the free energy of disordered system on temperature can be presented in the following form where ϕ is crossover critical exponent, which characterizes influence of defects on critical behavior.
As it is seen from (23), influence of defects on critical behavior for ϕ > 0 is relevant and leads to a new type of critical behavior with new values of critical exponents, while for ϕ < 0 influence of defects on critical behavior is irrelevant and critical behavior of disordered system is characterized by values of critical exponents for pure system.
For uncorrelated point-like defects, d d x u(x)G(x/ξ) ∼ const that leads to ϕ = 2−Dν 0 = α 0 . Therefore, the presence of point-like defects in sample is relevant for critical behavior if critical exponent of the specific heat for pure system α 0 is positive. This statement forms a content of well-known Harris criterion [14].
This paper devotes to study of influence of quenched point-like defects on properties of nonequilibrium critical behavior and we must note that only three-dimensional Ising-like systems are characterized by positive value of the specific heat critical exponent α 0 : 0.1091(24) (RG calculations), 0.1109(15) (MC calculations), 0.1099(7) (experiment). The 2D Ising model is a marginal case with the specific heat C v ∼ ln |τ | and suitable exponent α 0 = 0.
In case of systems with ε D -extended or isotropic long-range correlated defects, the crossover exponent ϕ = α 0 + ε D ν 0 and ϕ = α 0 + (D − a)ν 0 = 2 − aν 0 , respectively, and defects affect on critical behavior of more wide class of systems than Ising-like systems. We recommend for familiarity with question of influence of ε D -extended and isotropic long-range correlated defects on critical behavior the articles [15,16].
Investigations into the critical behavior of diluted Ising type magnets with the use of RG methods, MC numerical simulations, and experimental techniques confirm the existence of a new universality class of critical behavior, but the dependence of asymptotic values of critical exponents on the degree of their dilution and the influence of crossover effects in weakly and strongly disordered systems await elucidation and remain subjects of ardent discussion. It is worthwhile to note that analytical RG methods applied to study the critical behavior of impurity systems are suitable only for weakly diluted magnetic materials with concentration of defects c imp = 1 − p 1. As a system becomes increasingly more diluted with nonmagnetic impurity atoms at spin concentrations p (s) c and form a fractal-like structure with effective long-range spatial correlation in impurity distribution [16].
A change in effects of order parameter fluctuation scattering from impurity atoms must give rise to new fixed points for the vertices of interaction between order parameter fluctuations.
is characterized by a new type of critical behavior of three-dimensional Ising models, corresponding to the region of the strong structural disorder.
Presence of point-like defects in the 2D Ising model does not changes static critical exponents but logarithmic corrections may become important in the vicinity of the critical point. It has been found analytically within bond disorder model for the specific heat C v ∼ ln | ln |τ | and correlation length ξ ∼ |τ | −1 (1 + g ln(1/|τ |)) 1/2 with g 4.8(1 − p)/p [17] as well as for the magnetic susceptibility χ ∼ |τ | −7/4 ln 7/8 |τ | [18]. This prediction for weak bond disorder was confirmed numerically in a number of papers [19][20][21]. In this paper we study the 2D site-diluted Ising model, for which the phase diagram contains two crucial points, the pure Ising fixed point and the percolation fixed point at p (s) c 0.5927. This model was investigated numerically in a number of papers [22][23][24][25]. Some of authors claim in [22,23] that the static critical exponents β, γ vary with defect concentration, but the ratios β/ν, γ/ν do not change. Others [24,25] concluded that defects lead only to logarithmic corrections as for weak bond dilution. At the same time, the studies of the critical dynamics of the model reported in [26,27] demonstrated that, near the spin percolation threshold (for the systems with the spin concentrations p ≤ 0.85), the dynamic critical exponent z, determining the temperature dependence of the relaxation time, exhibits the dependence on the concentration of defects with the violation of the standard form of dynamic scaling and takes the form This suggests a rather nontrivial effect of defects on the characteristics of the non-equilibrium behavior of the 2D Ising model.

Non-equilibrium relaxation of the magnetization in 3D and 2D Ising models
In this part, we present the results of numerical Monte Carlo study of the effects of initial states and structural defects on relaxation properties of the magnetization M (t) in the 3D and 2D Ising models. The Hamiltonian of a site-diluted Ising model with external magnetic field h is given by the expression where J > 0 is the short-range exchange interaction between spins S i = ±1 fixed at the sites of a lattice, p i are the random occupation numbers indicating the presence of quenched structurally uncorrelated disorder in the system: p i = 1 for the i site with a spin, and p i = 0 for the site with a nonmagnetic impurity atom. Structural defects are canonically distributed in the system according to the distribution function P ( where p = p i defines the spin density in the system. The position of structural defects is fixed for a particular impurity configuration.
The simulation was performed on the simple cubic lattice for the 3D Ising model and on the planar square lattice for the 2D Ising model with imposed periodic boundary conditions. N s = pL D characterizes the number of spins in the lattice with linear size L.
To calculate characteristics of the non-equilibrium critical behavior, the dynamic process of one-spin flips was numerically simulated within the Monte Carlo statistical method. The dynamic process of one-spin flips was implemented using the Metropolis algorithm and the heat bath algorithm [28]. For an Ising model with two possible spin states S j = ±1, the probability of a spin flip with the use of the heat bath algorithm can be represented in the form Non-equilibrium behaviour of the magnetization relaxed from an initial state with magnetization m 0 is characterized at the critical point by the following scaling dependence The curves M (t) clearly show the significant both qualitative and quantitative differences in magnetization relaxation from a high-temperature initial state with m 0 1, a low-temperature completely ordered state with m 0 = 1, and intermediate states with 0.1 ≤ m 0 ≤ 0.6. We can see that the relaxation curves for systems started from initial states m 0 = 1 asymptote to the relaxation curve from a low-temperature initial state with m 0 = 1. In this case, for systems with m 0 1, a characteristic increase in the magnetization is observed at the non-equilibrium evolution stage, which is described by the power law M (t) ∼ t θ with θ = 0.111(4) for p = 1.0, θ = 0.127(16) for p = 0.8, and θ = 0.167 (18) for p = 0.6. At times t > t cr ≡ t m ∼ m −k 0 , this evolution stage is followed by the mode characterized by the power-law time dependence of the magnetization M (t) ∼ t −β/zν . In the case of system evolution from the initial ordered state with m 0 = 1, the time dependence of the magnetization at the critical point is immediately defined by     defects on relaxation of the 3D Ising model. We can see that the presence of structural defects slows down the system relaxation, leading to an increase in the relaxation time with the defect density. In Figs.1-2 we also demonstrate the results of numerical tests of predictions of the time scaling dependence (27) for the magnetization M (t, t m ) as a function of initial values of the magnetization m 0 for the 3D pure Ising model and systems with p = 0.8 and p = 0.6. For the curves of the scaling function F M (x) = t β/zν M (t, t m ) of the variable x = tm k 0 , a "collapse" of the data obtained for various m 0 is observed in the unified universal curve whose shape depends on the system spin density p. The universal curves for the scaling function F M (x) in the loglog plot are characterized by a linear initial part corresponding to the power-law dependence F M (x) x 1/k . By the slope of these linear parts, the following exponents k for the 3D Ising model were calculated: k = 2.77(2) for p = 1.0, k = 2.83(3) for p = 0.8, and k = 3.08(7) for p = 0.6.
The similar behavior are demonstrated in Fig.4 by the curves of the critical magnetization relaxation M (t) and its scaling function F M (t/t m ) = M (t)t β/zν for the 2D site-diluted Ising model. Differences against the 3D Ising model are determined by smaller values of exponent β/zν which is equal, for example, 0.046(1) for p = 0.8 while β/zν = 0.224 (10) for the 3D Ising model with the same dilution p = 0.8. It leads to increase of the critical slowing-down effects in the 2D Ising model and gives possibility better to study non-equilibrium behavior of this model numerically.  Figure 6. Dependence of the exponent z on spin density p at logarithmic scale | ln(p − p s c )|.
Particularly, we show in Fig.5 the procedure of taking into consideration of the logarithmic corrections for description of the time dependence of the magnetization in the case of evolution from the low-temperature initial state with m 0 = 1 using these corrections in the form: where g = 4 [17,30]. In result, we calculated influence of the logarithmic corrections on values of the exponents β/(νz) and z which demonstrate logarithmic dependence of the exponent z on the spin concentrations (Fig.6). The relation (24) is confirmed in the case of the spin concentrations p ≤ 0.85. However, the values of the exponents β/(νz) and z for systems with p = 0.95 and p = 0.9 coincide within the statistical errors with the values for the pure system β/(νz) = 0.05784 (30) and z = 2.161 (11). This result confirms that the values of the dynamic critical exponent z for weakly disordered systems with p ≥ 0.9 remain constant and the critical dynamics of the 2D site-diluted Ising model belongs to the universality class of the pure model [26,27]. The value of z for the pure system is in a good agreement with the value z = 2.1667(5) from the work [32]. We must note that the logarithmic corrections provide a marked contribution to the values of z.
Analysis of the scaling functions F M (tm k 0 ) for the 2D Ising model gives to calculate the exponent 1/k which has the values 1/k 0.245(11) for spin concentrations p = 1.0 ÷ 0.9 and 1/k 0.228(4)÷0.160(5) for p = 0.85÷0.7 with dependence on p by a logarithmic law | ln(p−p s c )|.

Aging in two-time behavior of the autocorrelation and the response functions.
Calculation of the FDR Aging effects during the non-equilibrium process demonstrate such two-time dependent quantities as the autocorrelation function C(t, t w ), and the linear response function R(t, t w ) to a weak field applied at time t w . The functions C(t, t w ) and R(t, t w ) can be defined as where the angle brackets stand for averaging over the different realizations of initial state and the square brackets are for averaging over the different impurity configurations. In simulation of system's dynamics with the use of the thermal bath algorithm, the response function in (30) is calculated using the following relations [35] where S W i (s) = tanh(J j =i p j S j /T ). The linear response function R(t, t w ) corresponding to definition (30), cannot be directly measured experimentally. The generalized susceptibility in the form of the integral response function (thermostatic susceptibility) can be measured experimentally and calculated numerically. The quantity ∆S (32) is calculated in the simulation process from initial state at t = 0 until the time t w , where S W i (s) is the same as in (31). Currently, it is known that, in the case of evolution from an initial state with m 0 = 0, the non-equilibrium critical behavior of the autocorrelation function, the response function, and dynamic susceptibility is characterized by the following scaling dependencies [3,33,34]: where the exponents θ and a are related to the critical exponents of the system under consideration: a = (2 − η − z)/z, θ = θ − (2 − z − η)/z. The scaling functions F C (t w /t, t/t m ), F R (t w /t, t/t m ), and F χ (t w /t, t/t m ) are finite at t w → 0 and t/t m → 0, A C , A R , and A χ are the non-universal amplitudes whose values are fixed by the conditions F C (0, 0) = 1, F R (0, 0) = 1, F χ (0, 0) = 1.
The quantities in (33) are generally homogeneous functions of three time scales tt w , t w , and t m . When the following their correspondence t w < t t m is satisfied, which is realized in the case of evolutions from a high-temperature initial state with m 0 = 0, dependencies (33) transform to the relations corresponding to this case [2,36].
In the opposite case with times t m t w < t, which is realized in the case of evolution from a low-temperature completely ordered initial state with m 0 = 1, the scaling dependencies (33) take the form where the new exponentθ = −βδ/(νz) = −(1 + a + β/(νz)) is introduced, andF C,R,χ are universal scaling functions obtained from F C,R,χ (t w /t, t/t m ) in (33) at t/t m → ∞.
In the aging regime realized at observation times t − t w ∼ t w t m , in which the twotime dependence of functions C(t, t w ), R(t, t w ), and χ(t, t w ) is most pronounced, the scaling  dependence of these functions on waiting times t w and t is given by relations with scaling functionsF C,R,χ (t/t w ), which decrease at the long-term stage of their variation with t − t w t w t m according to the power law where the exponent φ = d/z − a + βδ/zν [3]. For the high-temperature initial state with m 0 = 0, with Results of calculations of two-time dependencies of the autocorrelation function C(t, t w ) and the dynamic susceptibility χ(t, t w ) on observation times t − t w for various spin concentrations  p and waiting times t w with evolution from the high-temperature and low-temperature initial states are presented in Figs.7-9. The aging effects clearly manifest themselves at the observation times t − t w ∼ t w and are characterized by the slowing down of these functions with increase of the waiting time t w . The aging effects increase in the structurally disordered systems in comparison to those in the pure system.
Comparison of the autocorrelation function plots in Figs.7-8 with different initial states shows that the C(t, t w ) are decreased more slowly for the case of evolution from high-temperature initial state with m 0 1 than from the low-temperature initial state with m 0 = 1. The role of defects is the most clearly pronounced in the strong slowing down of the correlation effects in structurally disordered systems as compared to the pure system for the case of evolution from the low-temperature initial state (Fig. 8). We attribute these pronounced changes in the behavior of the autocorrelation function to the pinning of domain walls by structural defects. The pinning occurs in the course of non-equilibrium changes in the domain structure of the system at its transition from the single-domain state at T 0 = 0 to a multidomain fluctuation-induced structure arising at the critical temperature T c .
To confirm the scaling dependence for autocorrelation function C(t, t w ) and susceptibility χ(t, t w ) at the aging regime with t − t w ∼ t w (33-35), we draw the plots of t   2β/zν w χ(t, t w ) versus (t − t w )/t w in Figs.10-12 for cases with different initial states with the use of exponent β/zν which values were determined from study of relaxation properties of the magnetization. These plots demonstrate the collapse of the data obtained for different values of t w on the universal curves corresponding to different spin densities p and characterized by the scaling functions F C (t/t w ) and F χ (t/t w ) for case with evolution from the high-temperature initial state with m 0 1. However, the strong slowing down of the correlation effects in structurally disordered systems as compared to the pure system which was revealed for the case of evolution from the low-temperature initial state leads to violation of usual scaling behavior of function F C (t/t w ) on t/t w with non-availability of collapse of the data obtained for different values of t w on the universal curves. Putting into operation the scaling function F C (t/t µ w ) on t/t µ w with exponent µ > 1 restores the collapse of the data obtained for different values of t w on the universal curves as it is demonstrated in Fig.13. This phenomenon is classified as superaging and characterized by the exponent µ = 2.30(6) for weakly disordered 3D Ising systems and by µ = 2.80(7) for the 3D Ising systems with strong disorder. We revealed that the exponent µ for the 2D Ising model demonstrates the dependence on spin density p: µ(0.95) = 6.15(5), µ(0.9) = 6.25 (5), µ(0.85) = 6.50(5), µ(0.8) = 6.75, µ(0.75) = 7.00(5), and µ(0.7) = 7.50 (5). For scaling function F χ (t/t µ w ), the canonical aging with µ = 1 was revealed for both pure and diluted 3D and 2D Ising systems.
At the next stage of our studies, we calculated the fluctuation-dissipation ratio according to the relation [37,38]:  Table 1. The asymptotic values of the FDR X ∞ HT for case of system evolution from the hightemperature initial state and X ∞ LT the low-temperature initial state for 3D and 2D Ising models Dependencies T χ on C are presented in Fig.14 for the evolution of systems from different initial states. We calculated the FDR for each value of the waiting time t w , and then the asymptotic value of the FDR X ∞ was determined by extrapolation X(t w → ∞). The asymptotic values of the FDR are presented in Table 1 for case of system evolution from the high-temperature X ∞ HT and the low-temperature X ∞ LT initial states for the 3D and 2D Ising models. For the 3D Ising model similarly to exponents, calculated values of X ∞ HT (p) could consider as universal characteristics of three different classes of non-equilibrium critical behavior for pure,