Ageing and non-equilibrium critical phenomena in Monte Carlo simulations

It is considered the non-equilibrium critical evolution of statistical systems which displays some features, such as ageing and violation of the fluctuation-dissipation theorem. We review here some theoretical results of computations that have been obtained in recent years for universal quantities, such as the exponents determining the scaling behaviour of dynamic response and correlation functions and the fluctuation-dissipation ratio, associated with the non-equilibrium critical dynamics, with particular focus on the 3D pure and diluted Ising models with Glauber dynamics. We analyse an influence of critical fluctuations, different non-equilibrium initial states and presence of nonmagnetic impurities in spin systems on two-time dependence of correlation and response functions on characteristic time variables as waiting time tw and time of observation t – tw with t > tw. We discuss the obtained values of non-equilibrium exponents for autocorrelation and response functions and values of the universal long-time limit of the fluctuation-dissipation ratio X∞. Analysis of simulation results show that the insertion of disorder leads to new values of X∞ with X∞diluted > X∞pure.


Introduction
The collective behaviour of statistical systems close to critical points is characterized by an extremely slow dynamics which, in the thermodynamic limit, doesn't give them to relax to an equilibrium state after a change in some thermodynamic parameters. The non-equilibrium evolution in this case displays some of the features typically observed in glassy materials, such as ageing and violation of the fluctuation-dissipation theorem (FDT) [1][2][3][4][5]. It can be detected through determination of dynamic susceptibilities and correlation functions of the order parameter, the scaling behaviour of which is characterized by universal exponents, scaling functions, and amplitude ratios. This universality allows one to calculate these quantities in suitable simplified models [6,7] and Monte Carlo methods are a natural way for this analysis.
We consider here some of the theoretical results of computations that have been obtained in recent years for universal quantities, such as the fluctuation-dissipation ratio, associated with the non-equilibrium critical dynamics, with particular focus on original results of authors for the 3D pure and diluted Ising models with Glauber dynamics.
The investigation of critical behaviour of disordered systems remains one of the main problems in condensed-matter physics and excites a great interest, because all real solids contain structural defects. The structural disorder breaks the translational symmetry of the crystal and thus greatly complicates the theoretical description of the material. The influence of disorder is particularly important near the critical point where behaviour of a system is characterized by anomalous large response on any even weak perturbation. The description of such systems requires the development of special analytical and numerical methods. The effects produced by weak quenched disorder on critical phenomena have been studied for many years [8][9][10][11][12][13][14]. According to the Harris criterion [8] the disorder affects the critical behaviour only if α, the specific heat exponent of the pure system, is positive. In this case a new universal critical behaviour, with new critical exponents, is established. In contrast, when α < 0, the disorder appears to be irrelevant for the critical behaviour. Only systems whose effective Hamiltonian near the critical point is isomorphic to the Ising model satisfy this criterion (α > 0) for drastic change of universality class of the critical behaviour with introduction to system of whatever concentration of short-range correlated quenched defects [15].
A large number of publications is devoted to the study of the critical behaviour of diluted Ising-like magnets by the renormalization-group (RG) methods, the numerical Monte Carlo methods, and experimentally (for a review see [16]). The ideas about replica symmetry breaking in the systems with quenched disorder were presented in Refs. [17,18]. A refined RG analysis of the problem has shown the stability of the critical behaviour of weakly disordered threedimensional systems with respect to the replica symmetry breaking effects [19,20]. All obtained results confirm the existence of a new universal class of the critical behaviour, which is formed by diluted Ising-like systems. However, it remains unclear, whether the asymptotic values of critical exponents are independent of the rate of dilution of the system, how the crossover effects change these values, and whether two or more regimes of the critical behaviour exist for weakly and strongly disordered systems. These questions are the subjects of heated discussions [21,22].
In the present paper, we numerically investigate the non-equilibrium critical dynamics with a non-conserved order parameter (model A) [23] in the 3D pure and site-diluted Ising systems with spin concentrations p = 0.8 and 0.6.

Non-equilibrium critical dynamics, main peculiarities and characteristics
Statistical systems with slow dynamics have recently attracted considerable theoretical and experimental interest, in view of the rich scenario of phenomena they display: dramatic slowing down of relaxation processes, hysteresis, memory effects, etc. After a perturbation a system with slow dynamics does not achieve equilibrium even after long times and its dynamics is not invariant either under time translations or under time reversal, as it should be in thermal equilibrium. During this never-ending relaxation ageing occurs: two-time quantities such as response and correlation functions depend on characteristic time variables as waiting time t w and time of observation t − t w with t > t w not via t − t w only. Important, that decays for these quantities as functions of time of observation t − t w are slower for larger waiting times t w .
At variance with one-time quantities like the order parameter converging to asymptotic values in the long-time limit, two-time quantities clearly demonstrate the ageing. Ageing was known to occur in such disordered and complex systems as glassy materials [1][2][3] and only in the last ten years has attention been focused on simpler systems such as critical ones, whose universal features can be rather easily investigated by using different methods and which might provide insight into more general cases [24][25][26].
Consider a system with a critical point at temperature T c , order parameter S(x, t), and prepare it in some initial configuration which might corresponds to an equilibrium state at a given temperature T 0 . At time t = 0 bring the system in contact with a thermal bath with temperature T s not equal T 0 . The relaxation process is expected to be characterized by some equilibration time t eq (T s ) such that for t ≫ t eq (T s ) equilibrium is attained and the dynamics is stationary and invariant under time reversal, whereas for 0 < t ≪ t eq (T s ) the evolution depends on the specific initial condition. Upon approaching the critical point T s = T c the equilibration time diverges t eq ∼ |T − T c | −zν , where z is the dynamic critical exponent, ν is the exponent for correlation length, and therefore equilibrium is never achieved. To monitor the time evolution we where < ... > stands for the mean over the stochastic dynamics, and the linear response (susceptibility) R x (t, t w ) to a small external field, applied at time t w , which is defined by relation Causality implies R(t, t w > t) = 0. According to the general picture of the relaxation process, one expects that for t > t w ≫ t eq (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 the corresponding equilibrium quantities, related by the fluctuation-dissipation theorem (FDT) The FDT suggests the definition of the so-called fluctuation-dissipation ratio (FDR) with t > t w . For t > t w ≫ t eq (T s ) the FDT yields X(t, t w ) = 1. This is not generically true in the ageing regime. The asymptotic value of the FDR is a very useful quantity in the description of systems with slow dynamics, since X ∞ = 1 whenever the ageing evolution is interrupted and the system crosses over to equilibrium dynamics, i.e. t eq (T s ) < ∞. Conversely, X ∞ ̸ = 1 is a signal of an asymptotic non-equilibrium dynamics. Moreover, X ∞ can be used to define an effective non-equilibrium temperature T eff = T /X ∞ , which might have some features of the temperature of an equilibrium system, e.g. controlling the direction of heat flows and acting as a criterion for thermalization. What is known in general about X ∞ ? As a consequence of the fluctuation-dissipation theorem X ∞ (T > T c ) = 1. On the other hand, on the basis of general scaling arguments for the phase ordering regime [27], it has been shown that X ∞ (T < T c ) = 0. These results are expected to be actually independent of the specific system and of its microscopic details. For T = T c there are no general arguments constraining the value of X ∞ and therefore it has to be determined for each specific model. In table 1 we report some of the values that has been found either by exact solutions or by means of Monte Carlo (MC) simulations (a more complete table can be found in [24]). Clearly, X ∞ (T = T c ) depends on the model and on the space dimensionality d. Nevertheless it has been argued on the basis of scaling arguments [27,29] that X ∞ (T = T c ) should be a universal quantity associated with the critical dynamics.
At present time, it is well known that two-time dependence for autocorrelation and response functions for systems starting from high-temperature initial state with m 0 = 0 (or m 0 ≪ 1) satisfies the following scaling forms:  [34]. A R and A C are non-universal amplitudes which are fixed by the condition f R,C (0) = 1. With this normalization f R,C are universal. From these scaling forms the universality of X ∞ follows as an amplitude ratio (an example of M (t) behaviour is given in Fig.7). The initial rise of magnetization is changed to the well known decay M (t) ∼ t −β/zν for t ≫ t cr . The critical exponents θ and θ ′ depend on the dynamic universality class [23] and have been calculated by the RG method for a number of dynamic models such as the model with a non-conserved order parameter [22,[34][35][36] (model A), the model with an order parameter coupled to a conserved density [37] (model C), and the models with reversible mode coupling [38] (models E, F, G, and J).
It can be distinguished three stages of the non-equilibrium relaxation process. The first quasiequilibrium stage is observed for small time separation t − t w ≪ t w with t w ≫ 1 where ageing does not exist and the dynamic evolution of the correlation and response functions exhibits a stationary part and does not depend on waiting time with C = C(t−t w ) and R = R(t−t w ). The second ageing regime is realized for times t − t w ∼ t w ≫ 1 where the correlation and response functions are derived by relations and therefore at different waiting times do not superpose and are characterized by different slopes for each t w (in (7) it was used the relation 2β/(νz) = d/z − a − 1). At long time separations with t − t w ≫ t w ≫ 1 the scaling functions in (7) decay as power laŵ where exponent c a = d/z −θ ′ is the same which describes time dependence of the autocorrelation function in the short-time regime of non-equilibrium behaviour [22,39,40]. At this short-time dynamics stage ageing effects are not developed too. Scaling analysis of response function R(t, t w ) behaviour in short-time dynamics regime predicts that c r = c a . Renormalization-group investigations of non-equilibrium critical behaviour in d-dimensional pure systems with n-component order parameter and weakly dilute Random Ising Model (RIM)  [41] and [42]. The asymptotic values of the FDR X ∞ were calculated with the use of ε-expansion method (ε = 4−d) in two-loop approximation for pure systems [41] ( and in one-loop approximation for diluted Ising model [42] It was obtained X ∞ 3DIs = 0.429(6) for the three-dimensional Ising (ε = 1, n = 1) model, (8) for XY (ε = 1, n = 2) model, and X ∞ 2DIs = 0.30(5) for the two-dimensional Ising (ε = 2, n = 1) model which are in good agreement with Monte Carlo results for these models given in Table 1. For 3D diluted Ising model (that is the only physical relevant case) it was obtained value X ∞ 3DRIM ≃ 0.416. To this order it is not clear whether disorder really changes X ∞ in a sensible way or not. In any case this could not be safely stated from low-order computations since the √ ε-expansion is known to be not well-behaved at d = 3 [16,20,43]. Results of Monte Carlo investigations of non-equilibrium critical dynamics for the 3D diluted Ising model will be presented in this article below.
If initial state of system is characterized by magnetization m 0 ̸ = 0 the renormalizationgroup analysis of non-equilibrium dynamics for model A predicts that the correlation C(t, t w ) and response R(t, t w ) functions display the following scaling behaviours after a quench to T = T c [34,44] where modification of these scaling relations in comparison with (6) is connected with new time scale t m set by the initial value of the magnetization m 0 and which displays a universal dependence on it where the universal scaling exponent κ > 0 is given, in terms of static and dynamic equilibrium and non-equilibrium exponents, by κ = θ + a + β/(νz). The two-time quantities C(t, t w ) and R(t, t w ) are homogeneous functions of the three time scales t − t w , t w , and t m . In particular, when t w < t ≪ t m , which is always for case with m 0 = 0, the scaling form of C and R becomes as in equation (6) with F C,R (x, 0) = f C,R (x). In the opposite case with large times compared to t m , i.e. t m ≪ t w < t the scaling form of C(t, t w ) and R(t, t w ) becomes [44]: where new exponentθ = −βδ/(νz) = −(1 + a + β/(νz)) andF C,R are universal scaling functions related to the large-y behaviour of F C,R (x, y). In the ageing regime which is realized for times t − t w ∼ t w ≫ t m the correlation and response functions are derived by relations with the scaling functionsF C,R (t/t w ) which decay at long time separations limit with t − t w ≫ t w ≫ t m as power lawF where exponent ϕ = d/z − a + βδ/(νz).
In paper [44] it was investigated the non-equilibrium behaviour of the d-dimensional Ising model with purely dissipative dynamics during its critical relaxation from a magnetized initial state. The universal scaling forms of the two-time response and correlation functions of the magnetization were derived within the field-theoretical approach and the associated scaling functions were computed in first order of the ε-expansion. It was shown that ageing behaviour is clearly displayed and the asymptotical universal fluctuation-dissipation ratio is characterized by relation which gives X ∞ 3DIs ≃ 0.78 for the three-dimensional Ising (ε = 1, n = 1) model and X ∞ 2DIs ≃ 0.75 for the two-dimensional Ising (ε = 2, n = 1) model. The obtained results were confirmed by Monte Carlo simulations of the two-dimensional Ising model with Glauber dynamics, from which it was found X ∞ M C = 0.73(1) [44].

Monte Carlo simulations of the 3D pure Ising models with relaxation from high-temperature and low-temperature initial states
The one from simplest non-trivial models in which ageing occurs is the Ising model in ddimensions evolving with a purely dissipative dynamics after a quench to the critical point.

Its Hamiltonian on an hypercubic lattice is given by
where J > 0 is the short-range exchange interaction between spins S i fixed at the lattice sites, and assuming values of S i = ±1. We performed our Monte Carlo simulations using the heat bath updating rule [45], simulating a large Ising spin system on cubic lattice with linear size L = 128 with periodic boundary conditions. It were computed the magnetization M (t) and the two-time autocorrelation function C(t, t w ) where the angle brackets stand for an average over initial configurations and MC realizations.
The averaging of C(t, t w ) is carried out on 3000 MC runs for every t w .
For case when we have simulated the dynamics of the systems starting from the hightemperature initial state with magnetization value m 0 = 0.02, the response function and fluctuation-dissipation ratio were calculated using the relations [32] where S W i = tanh(J ∑ m̸ =i S m /T ) and The averaging of R(t, t w ) and X(t, t w ) is carried out on 90000 MC runs for every t w .
For case when we have simulated the dynamics of the systems starting from the lowtemperature initial state with magnetization value m 0 = 1 we calculated the integrated response function [46]: with response function determined by (2). In the large time limit χ(C) = ∫ 1 C X(q)dq and the fluctuation-dissipation ratio can be defined as The MC obtained time dependencies of the autocorrelation function C(t, t w ) and the response function R(t, t w ) from observation time t − t w for different initial non-equilibrium states and t w values are presented in Fig.1 and 2. These functions demonstrate the ageing effects that is the slowing down of time correlations and decreasing of response with increasing of system age t w .
In the ageing regime the time dependence of correlation functions is characterized by scaling relations (7) and (14). These scaling relations are rather well displayed by our data, as shown in Fig. 3. For time intervals with (t − t w )/t w ≫ 1 the values of exponents c a = 1.333 (40) and c r = 1.357 (18) were determined which demonstrate a good agreement with each other and with value c r = 1.362 (19) [39] obtained by short-time dynamics method.
The obtained data for fluctuation-dissipation ratio for case with m ≪ 1 are presented in Fig. 5. The asymptotical value of FDR X ∞ = 0.380 (13) was obtained by linear extrapolation of data for X(t w /(t − t w )) to limit t w /(t − t w ) → 0.
For completely ordered initial state with m 0 = 1 the asymptotical value of the FDR X ∞ = 0.77(6) was obtained as realization of limit dependence T χ(C) (Fig. 6) in accordance with relation (22). This value of X ∞ is in excellent agreement with theoretical field value X ∞ ≃ 0.78 [44].

Monte Carlo simulations of the 3D diluted Ising models with relaxation from high-temperature initial state
The Hamiltonian of diluted by nonmagnetic impurity atoms ferromagnetic Ising model is given by where random occupation numbers p i take the values 0 or 1, and p i equal to 1 if site of lattice contains spin and 0 otherwise. We considered the cubic lattice with periodic boundary conditions. Let us denote, N s = pL 3 is the number of spins in cubic lattice with linear size L and p is the spin concentration. We performed our Monte Carlo simulations also with the use of the heat bath updating rule for Ising spin system with the spin concentrations p = 0.8 and 0.6 on cubic lattice with L = 128.  Figure 5.
It were computed the magnetization and the two-time autocorrelation function where the square brackets stand for additional average over disorder configurations. The averaging of C(t, t w ) was carried out on over 1000 samples with different disorder configurations with 15 MC runs for each sample.
In this part of investigations we have simulated the non-equilibrium critical dynamics of the systems starting from the high-temperature initial state only with magnetization value m 0 = 0.01 for p = 0.8 and m 0 = 0.005 for p = 0.6. The response function and fluctuation-dissipation ratio were calculated using the relations where S W i = tanh(J ∑ m̸ =i S m /T ). The averaging of R(t, t w ) and X(t, t w ) was carried out on over 5000 samples with different disorder configurations with 15 runs for each sample.
We present in Fig. 7 the time evolution of the magnetization M (t) from different initial states with magnetization m 0 = 0.02 at T c = 4.5114(1) for pure Ising model, m 0 = 0.01 at T c = 3.4995(2) for samples with spin concentration p = 0.8, and m 0 = 0.005 at T c = 2.4241(1) for samples with p = 0.6 [21,22] which are characterized by power-law dependence in form M (t) ∼ t θ ′ in initial stage of evolution. The time scale of a critical initial increase of the magnetization is t cr ∼ m . The initial rise of magnetization is changed to the wellknown decay m(t) ∼ t −β/zν for t ≫ t cr . However, in the limit of m 0 → 0, the time scale t cr goes to infinity. The times with t < t cr can be considered as time intervals of non-equilibrium behaviour of considered systems. The values of exponent θ ′ calculated with the use of corrections to scaling are following: θ ′ = 0.106(4) for p = 1, θ ′ = 0.127 (16) for p = 0.8, and θ ′ = 0.186 (39) for p = 0.6 [22,49,50].
The time dependence of autocorrelation C(t, t w ) and response R(t, t w ) functions is plotted in Fig. 8 for different values of waiting time t w and spin concentrations p. These curves demonstrate that the ageing effects are increased with growth of defect concentrations.
To check the scaling dependence of C(t, t w ) and R(t, t w ) given by relations (7) we plotted t 2β/(νz) C(t, t w ) and t 1+2β/(νz) R(t, t w ) versus t/t w with the use of exponent values: z = 2.191(21), 2β/ν = 1.016(32) for p = 0.8 [22] and z = 2.663 (30), 2β/ν = 0.924(80) for p = 0.6 [49,50]. The obtained data are presented in Fig. 9. These functions demonstrate the collapse of curves for different t w with fixed spin concentration p into single curve with universal scaling dependence. Note, that systems with different p are characterized by different scaling functions F C,R (t, t w , p).  [22,49,50] by short-time dynamics method, but are in a poor agreement with value c a = 1.05(3) from paper [51] obtained as independent on p characteristics of non-equilibrium behaviour of the autocorrelation function. Reasons for this discrepancy have been discussed in detail in our paper [22]. At later investigations we have computed the fluctuation-dissipation ratio in compliance with relation (25). The obtained data are plotted in Fig. 10 for different spin concentrations p. For analysis we used time interval from t − t w ∼ t w to t − t w ≫ t w . In contrast to pure Ising model the data of X(t, t w ) for the diluted Ising model are characterized by an explicit dependence from t w . At the beginning, we calculated the asymptotical values of X(t w , p) for different t w from plot in Fig. 10 in the limit t w /(t − t w ) → 0, and then, using the obtained values of X(t w , p), we made extrapolation 1/t w → 0 to gain the asymptotical fluctuation-dissipation ratios X ∞ (p). Dependence FDR X(t w , p) as function of 1/t w for different p. Asymptotical value of X ∞ (p) can be obtained in limit 1/t w → 0 The results are presented in Fig. 11. The final values of the asymptotical fluctuation-dissipation ratio are X ∞ (p = 1) = 0.380 (13), X ∞ (p = 0.8) = 0.413 (11), and X ∞ (p = 0.6) = 0.446 (8). We are giving in table 2 the values for critical universal characteristics like exponents c a , c r and the asymptotical fluctuation-dissipation ratio X ∞ , obtained at present paper for different spin concentrations p = 1.0, 0.8, and 0.6, and values of the critical exponents 2β/ν, z, and θ ′ for these spin concentrations which were used in procedure of c a and c r determination and were calculated by Monte Carlo simulations of non-equilibrium critical dynamics in Refs. [22,39,49,50].

Conclusions
In this paper we have reviewed some theoretical results of computations in the range of nonequilibrium phenomena that have been obtained in recent years for universal quantities, such as the exponents determining the scaling behaviour of dynamic response and correlation functions and the fluctuation-dissipation ratio, associated with the non-equilibrium critical dynamics. In original part, we have considered the results of our MC simulations for the 3D pure and diluted Ising models with Glauber dynamics. The influence of critical fluctuations, different nonequilibrium initial states and site-quenched disorder on two-time dependence of the correlation and response functions at non-equilibrium critical regime were investigated for these models.
Analysis of the time dependencies for the autocorrelation and response functions at ageing   (8) regime of evolution showed that the ageing effects are increased with growth of defect concentrations. The obtained values of the autocorrelation critical exponents c a and the fluctuation-dissipation ratio X ∞ demonstrate that these systems belong to different classes of non-equilibrium critical behaviour. Analysis of simulation results show that the insertion of disorder leads to new universal FDRs with X ∞ strong diluted > X ∞ weak diluted > X ∞ pure .