# Universality in the off-equilibrium critical dynamics of the 3⁢𝑑 diluted Ising model

## Abstract

We study the off-equilibrium critical dynamics of the three dimensional diluted Ising model. We compute the dynamical critical exponent $`z`$ and we show that it is independent of the dilution only when we take into account the scaling-corrections to the dynamics. Finally we will compare our results with the experimental data.

PACS numbers: 05.50.+q, 75.10.Nr, 75.40.Mg

The issue of Universality in disordered systems is a controversial and interesting subject.

Very often in the past it has been argued that critical exponents change with the strength of the disorder . While, on a deeper analysis, it has turned out that those exponents were “effective” ones, i.e. they are affected by strong scaling corrections. So, when one studies the critical behavior of a disordered system it is mandatory to control the leading correction-to-scaling in order to avoid these effects that could modify the dilution-independent values of the critical exponents. For instance, in Ref.  the equilibrium critical behavior of the three dimensional diluted Ising model was studied. The authors showed that taking into account the corrections-to-scaling it was possible to show that the static critical exponents (e.g. $`\nu `$ and $`\eta `$) and cumulants were dilution-independent. These numerical facts supports the (static) perturbative renormalization group picture: all the points of the critical line (with $`p<1`$) belong to the same Universality class (with critical exponents given by the random fixed point) . Their final values of the exponents  were in very good agreement with the experimental figures (see below).

We will show that an analogous effect also happens in the off-equilibrium dynamics of the diluted ferromagnetic model and we will take it into account in our data analysis, in order to get the best estimate of the critical dynamical exponent.

The critical dynamics of the diluted Ising model has been studied experimentally in Ref.  using neutron spin-echo inelastic scattering on samples of $`\mathrm{Fe}_{0.46}\mathrm{Zn}_{0.54}\mathrm{F}_2`$ (antiferromagnetic diluted model) and has been compared with the results obtained in pure samples ($`\mathrm{FeF}_2`$. For the pure model a dynamical critical exponent $`z=2.1(1)`$ was found (in good agreement with the theoretical predictions based on one-loop perturbative renormalization group (PRG) ) whereas in the diluted case the exponent $`z=1.7(2)`$ was computed (three standard deviations away of the analytical prediction based on (one-loop) PRG that provides $`z2.34`$ ). Furthermore, the dynamical exponent was computed in the framework of the PRG up to two loops and it was obtained $`z=2.237`$  and $`z=2.180`$  (the experimental value is at 2.5 standard deviation of the two loops analytical result).

In the experiment  were measured critical amplitudes 100 times smaller than those computed in the pure case. It is clear that a more precise experiment on this issue will be welcome. We should point out that the critical dynamics of a diluted antiferromagnet is the same as of a diluted ferromagnet.

A numerical study of the on-equilibrium dynamics in diluted systems was performed in 1993 by Heuer . He measured the equilibrium autocorrelation functions for different concentrations and lattice sizes. The autocorrelation time ($`\tau `$) depends on the lattice size ($`L`$) via the formula $`\tau L^z`$ (neglecting scaling corrections). He found that all the data, for concentrations not too close to 1, were compatible, for large $`L`$, with the assumption of a single dynamical exponent, different from the one of the pure fixed point and similar to the analytical estimate of Ref.  ($`z2.3`$). The final value reported by Heuer was $`z=2.4(1)`$.

The main goal of this work is to check Universality in the critical dynamics of diluted models (i.e. whether the dynamical critical exponent is dilution independent) in the off-equilibrium regime . To do this we monitor scaling corrections in the same way it was done in the static simulations . Therefore we will also obtain the value of the corrections-to-scaling exponent for the dynamics. Our motivation to study the off-equilibrium dynamics instead of the equilibrium one is based on two reasons. The more important reason is that the experimental data was obtained in the off-equilibrium regime and the second one is that (in general) it is easier to simulate systems in the off-equilibrium regime. Moreover it will be possible to confront our $`z`$ computed in the off-equilibrium regime with that obtained at equilibrium .

The relevance of the corrections-to-scaling is twofold. The first one is that the scaling-corrections are very important in the right determination of the static (equilibrium simulation) critical exponents . In some models the corrections-to-scaling change the anomalous dimension of the order of 10 % (see for example Ref. ). The second one is that the correction-to-scaling exponent can be (and it has been) computed in a real experiment .

We have studied the three dimensional diluted Ising model defined on a cubic lattice of size $`L`$ and with Hamiltonian

$$=\underset{<ij>}{}ϵ_iϵ_jS_iS_j,$$
(1)   
where $`S_i`$ are Ising spin variables, $`<ij>`$ denotes sum over all the nearest-neighbor pairs and $`ϵ_i`$ are uncorrelated quenched variables which are 1 with probability $`p`$ and zero otherwise.

We have measured, at the infinite volume critical point and for several concentrations $`p`$, the non-connected susceptibility, defined by

$$\chi =\frac{1}{L^3}\underset{ij}{}\overline{S_iS_j},$$
(2)   
where the brackets stand for the average over different thermal histories or initial configurations and the horizontal bar for an average over the disorder realizations. The indices $`i`$ and $`j`$ run over all the points of the cubic lattice. In practice we use a large number of disorder realizations ($`N_S=512`$) each with a single thermal history, what amounts to neglect the angular brackets in Eq.(2). This procedure is safe and does not introduce any bias.

With the notation of the book of Ma  we can write, for instance, the following equation for the response function, under a transformation of the dynamical renormalization group (RG) with step $`s`$,

$`G(𝒌,\omega ,𝝁)=`$ (3)   
$`s^{2\eta }G(s𝒌,s^z\omega ,𝝁^{}\pm \left({\displaystyle \frac{s}{\xi }}\right)^{y_1}𝒆_1+O(s^{y_2})),`$ (4)   
where $`\omega `$ is the frequency, $`𝒌`$ is the wavelength vector, $`z`$ is the dynamical critical exponent, by $`𝝁`$ we denote all the parameters of the Hamiltonian, $`𝝁^{}`$ is the fixed point of the renormalization group transformation, $`\xi `$ is the static correlation length and finally $`y_1`$ is the relevant eigenvalue (equals to $`1/\nu `$: $`y_1`$ is the scaling exponent associated with the reduced temperature), $`𝒆_1`$ is its associate eigenvector and $`y_2`$ is the greatest irrelevant eigenvalue ($`y_2<0`$) of the renormalization group transformation (we have assumed that the system posses only one relevant operator).

Using Eq.(4) and considering the leading scaling-corrections for a very large system  at the critical temperature, we can write the dependence of the susceptibility on the Monte Carlo time as

$$\chi (t,T_c(p))=A(p)t^{\frac{\gamma }{\nu z}}+B(p)t^{\frac{\gamma }{\nu z}\frac{w}{z}},$$
(5)   
where $`t`$ is the Monte Carlo time, $`T_c(p)`$ is the critical temperature, $`A(p)`$ and $`B(p)`$ are functions that depend only on the spin concentration, $`\gamma `$ is the exponent of the static susceptibility, $`\nu `$ is the exponent of the static correlation length, $`z`$ is the dynamical critical exponent and finally $`wy_2`$ is the correction-to-scaling exponent. Hereafter we denote $`w_dw/z`$. We recall that $`w`$ corresponds with the biggest irrelevant eigenvalue of the RG in the dynamics, in principle $`w`$ will be different from the leading correction in the static (that we will denote by $`w_s`$. In addition, an analytical correction-to-scaling comes from the non singular part of the free energy and gives us a background to add to Eq.(5). In our numerical simulations we can neglect this background term (i.e. we will show that $`\gamma /(\nu z)\omega /z0.50`$).

Moreover, Eq.(5) is valid for times larger than a given “microscopic” time and for times (in a finite lattice) less than the equilibration time (that is finite in a finite lattice).

To study numerically the present issue we have simulated $`L=100`$ systems for different spin concentrations $`p=1.0,0.9,0.8,0.65,0.6,0.5`$ and $`0.4`$ at the critical temperatures reported in Ref. . The Metropolis dynamics  provides our local dynamics. We have checked in all the simulations that we were in an off-equilibrium situation: for the volumes and times we have used, the non-connected susceptibility is far from reaching its equilibrium plateau (in a finite system). For completeness we also report the numerical estimate of the critical exponents for the random fixed point, where all the systems with $`p<1`$ should converge for large length scales : $`\gamma =1.34(1)`$, $`\nu =0.6837(53)`$, $`\eta =0.0374(45)`$ and $`w_s=0.37(6)`$ (Ref. , using PRG, provides $`\omega =0.372(5)`$ in the massive scheme and $`\omega =0.39(4)`$ in the minimal subtraction one). It is worth noting that experimentally the best estimate of the susceptibility exponent is $`\gamma =1.33(2)`$ .

At this point we can recall the one-loop prediction of the PRG for the $`\nu `$ and $`\eta `$ exponents: $`\nu =\frac{1}{2}+\frac{1}{4}\sqrt{\frac{6ϵ}{53}}`$ and $`\eta =ϵ/106`$ , where $`ϵ=4d`$, $`d`$ being the dimensionality of the space. If we substitute $`ϵ=1`$ we obtain the following (“naive”) estimates: $`\eta =0.0094`$ and $`\nu =0.5841`$. Obviously the previous naive estimates are far from the numerical and experimental values of the critical exponents. This would also imply that even the one-loop PRG estimate of the dynamical critical exponents will stay far from the true value.

Notice also that the anomalous dimension exponent ($`\eta `$) takes near the same value either at the pure and at the random fixed point. One can argue that this holds using the arguments provide in Ref.  using a $`ϵ^{}`$-expansion (where $`d2+ϵ^{}`$. This fact, assuming the naive dynamical theory (Van Hove theory or conventional theory) , implies that the dynamical critical exponent $`z=\gamma /\nu =2\eta `$ is the same for both, diluted and pure system, to first order in $`ϵ^{}`$. We will show that this is not the case for our diluted model. The Van Hove theory was used in  to interpret the experimental data.

An analytical estimate of the value of the dynamical critical exponent has been taken from Ref.  where a dynamical $`\sqrt{ϵ}`$-expansion ($`ϵ4d`$) was done: $`z=2+\sqrt{6ϵ/53}+O(ϵ)`$, that in three dimension becomes $`z2.34`$, where we have neglected the terms $`O(ϵ)`$. We can recall the two loops computation: $`z=2.237`$  and $`z=2.180`$ . One of the results of this work should be about the reliability of the previous estimates of $`z`$ (the first and second term of an $`\sqrt{ϵ}`$-expansion).

With all these ingredients we can analyze our numerical data for the dynamical non-connected susceptibility and check whether or not the Universality, based on renormalization group arguments, holds.

In the first plot (Fig. 1) we show the numerical data in a double logarithm scale. The slope gives, neglecting the corrections-to-scaling, the ratio $`\gamma /(\nu z)`$. It seems that all the lines behave in a power law but with different slopes \[i.e. different exponents $`\gamma /(\nu z)`$\]. This fact could call for non universality in this model (i.e. critical exponents vary along the critical line). In addition, if we take into account the main result from the static  which states that the static critical exponents (e.g. $`\nu `$ and $`\eta `$) do not depend on the dilution degree, we obtain a dynamical exponent that depends on the dilution, violating the prediction of the dynamical perturbative renormalization group . In fact, following the RG flow (for $`p<1`$) we should always end at the same random fixed point and so, for large scales (in time and space) $`z`$ is not expected to depend on the dilution degree.

In the previous analysis we have not taken into account the scaling corrections. However, we are able to monitor the leading scaling-corrections given by the exponent $`w`$. We succeeded in fitting (using the MINUIT routine ) all our numerical data to Eq.(5) for $`0.5p0.8`$. We have 10 parameters to fit : $`A(p)`$ and $`B(p)`$ for four dilutions ($`p=0.8`$, $`0.65`$, $`0.6`$, $`0.5`$), $`\gamma /(\nu z)`$ and $`\omega /z`$, these last exponents assumed dilution independent.

In this way we have computed the functions $`A(p)`$ and $`B(p)`$ in Eq.(5) and $`\gamma /(\nu z)`$ and $`\omega /z`$. By fitting the data using $`t4`$ we have obtained a very good fit (with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=33.8/34`$, where d.o.f stands for degrees of freedom) and the following values for the dynamical critical exponent and the leading dynamical scaling-corrections

$$z=2.62(7),\omega =0.50(13),$$
(6)   
where we have used the value of the static critical exponents $`\gamma =1.34(1)`$ and $`\nu =0.6837(53)`$ .

In order to check the stability of the previous fit we have tried a new fit using only times $`t8`$. The fit again is very good (with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=29.7/30`$) and

$$z=2.58(7),\omega =0.72(16),$$
(7)   
Clearly the fit is very stable since both exponents are compatible inside the error bars (one half standard deviation in $`z`$ and one standard deviation in $`\omega `$). Therefore we take, as our final values, $`z=2.62(7)`$ and $`\omega =0.50(13)`$.

In Fig. 2 we show our results for the amplitudes $`A(p)`$ and $`B(p)`$ (using the results of the fit with $`t4`$; $`t=4`$ plays the role of the microscopic time for this model and algorithm, see the previous discussion). The main result of these fits is that the numerical data can be well described using dilution independent exponent (both dynamical and static), while the value of the dilution only enters in the non-universal amplitudes, $`A(p)`$ and $`B(p)`$. This fact clearly supports Universality in this model.

From Fig. 2 we can compute the value of the dilution in which there is not (leading) scaling-corrections (one kind of “perfect Hamiltonian” for this dynamical problem). For $`p0.63`$ we obtain $`B(p)0`$ and so with this dilution it is possible to measure dynamical critical exponents \[e.g. $`\gamma /(z\nu )`$ from the growth of the susceptibility, $`(d1/\nu )/(z\nu )`$ from the relaxation of the energy, etc.\] neglecting the underlying (leading) scaling-corrections. This dilution could be a good starting point in order to monitor the sub-leading scaling-corrections.

Systems with spin concentrations $`p=0.9`$ have also been simulated, but the data from these runs are not been included in the previous analysis, because they can not be well fitted with the formula of Eq.(5). We can explain this fact assuming that for this dilution the system is in the cross-over region, for the lattice and times we used. Also in the static studies a similar effect was found and only for $`p0.8`$ was possible to obtain final values (for exponents and cumulants) dilution independent .

In order to convince the reader of the goodness of our fits we plot in Fig. 3 the non-connected susceptibility divided by just the correction-to-scaling factor $`[A(p)+B(p)t^{w_d}]`$. If universality holds (i.e. all the critical exponents, dynamical and static, are dilution independent) all the data points (corresponding to four dilution degrees) should collapse on a straight line in a double logarithm scale. It is clear from this figure that it is what happens. The equation of the curve is $`t^{\frac{\gamma }{\nu z}}`$ with $`\frac{\gamma }{\nu z}=0.748`$.

We have shown that it is possible to describe the off-equilibrium numerical data assuming critical exponents (dynamical as well as static) independent on the dilution for a wide range of dilutions. This supports the predictions of the (perturbative) renormalization group for the statics as well as for the dynamics. So, the (perturbative) RG scenario that predicts that all the points on the critical line (for $`p<1`$) belong to the same Universality class is very well supported by numerical simulations.

We have found that our estimate of the dynamical critical exponent $`z=2.62(7)`$ is incompatible with the experimental value $`z=1.7(2)`$. Further numerical and experimental studies should be done in order to clarify this discrepancy.

We can compare the value of the dynamical critical exponent computed off- and on-equilibrium. The Heuer’s estimate was $`z=2.4(1)`$ and the difference with our estimate $`z=2.62(7)`$ is $`z_{\mathrm{off}\mathrm{eq}}z_{\mathrm{eq}}=0.22(12)`$, i.e. $`1.8`$ standard deviations. The conclusion is that both estimations are compatibles in the error bars. In any case, it will be interesting to compute $`z`$ on-equilibrium by controlling the scaling corrections.

Moreover, our estimate is not compatible with that of PRG to order $`\sqrt{ϵ}`$ in the $`\sqrt{ϵ}`$-expansion ($`z=2.34`$). The comparison with the two loops estimates of $`z`$  is still worse. One possible explanation for this disagreement could be the lack of Borel summability the diluted model shows . We remark again that the one loop PRG estimates of the static critical exponents was very bad (see below).

Another interesting issue is to compare the dynamical scaling-corrections and the static ones. Unfortunately our statistical precision is unable to solve this issue. For instance, taking the values of the $`t4`$ we obtain $`\omega \omega _d=0.13(6)`$ that is compatible with zero assuming two standard deviations. If we take the values of the $`t8`$ fit we obtain $`\omega \omega _d=0.35(17)`$. We will devote further work (analytical and numerical) in order to discern if the leading dynamical scaling-correction corresponds to the leading static scaling-correction.

We wish to thank H. G. Ballesteros, D. Belanger, L. A. Fernández, Yu. Holovatch, V. Martín Mayor and A. Muñoz Sudupe for interesting discussions.