Computation of the dynamic critical exponent of the three-dimensional Heisenberg model

TL;DR

This paper accurately computes the dynamic critical exponent of the three-dimensional Heisenberg model using advanced equilibrium and out-of-equilibrium techniques, providing results consistent with previous studies and experiments.

Contribution

It presents the first high-precision computation of the dynamic critical exponent for large lattices in the 3D Heisenberg model using both equilibrium and out-of-equilibrium methods.

Findings

Dynamic critical exponent z ≈ 2.03–2.04

Consistent static critical exponents η and ν

Agreement with previous small-lattice and theoretical estimates

Abstract

Working in and out of equilibrium and using state-of-the-art techniques we have computed the dynamic critical exponent of the three dimensional Heisenberg model. By computing the integrated autocorrelation time at equilibrium, for lattice sizes $L\le 64$, we have obtained $z=2.033(5)$. In the out of equilibrium regime we have run very large lattices ($L\le 250$) obtaining $z=2.04(2)$ from the growth of the correlation length. We compare our values with that previously computed at equilibrium with relatively small lattices ($L\le 24$), with that provided by means a three-loops calculation using perturbation theory and with experiments. Finally we have checked previous estimates of the static critical exponents, $\eta$ and $\nu$, in the out of equilibrium regime.

Tables

Table 1. Table 1: Integrated correlation time of 𝑴 2 superscript 𝑴 2 \mathitbf{M}^{2} , τ int , 𝑴 2 subscript 𝜏 int superscript 𝑴 2 \tau_{\mathrm{int},\mathitbf{M}^{2}} for c = 10 𝑐 10 c=10 , as a function of the lattice size, L 𝐿 L . We also report the length of the run at equilibrium, n sweep subscript 𝑛 sweep n_{\mathrm{sweep}} , in units of τ int , 𝑴 2 subscript 𝜏 int superscript 𝑴 2 \tau_{\mathrm{int},\mathitbf{M}^{2}} . For each lattice size we have performed 50 intial conditions. Notice that all the reported runs satisfy n sweep > 10000 ​ τ int , 𝑴 2 subscript 𝑛 sweep 10000 subscript 𝜏 int superscript 𝑴 2 n_{\mathrm{sweep}}>10000~{}\tau_{\mathrm{int},\mathitbf{M}^{2}} . Sokal

$L$	${\tau }_{\mathrm{int},{𝑴}^2}$	$n_{\mathrm{sweep}}/{\tau }_{\mathrm{int},{𝑴}^2}$
8	24.84(1)	4122383
12	53.11(4)	1928074
16	93.57(8)	1094368
24	211.2(3)	484848
32	378.3(5)	270542
48	860(3)	58333
64	1545(9)	13782

Table 2. Table 2: Hardware features of the CPUs and GPUs.

CPU/GPU	CPU Intel	GPU Geforce	GPU Tesla
model	Core i7	GTX 1080 G1	K80
Cores	20	2560	4992
Core clock	2.26 GHz	1.86 GHz	0.88 GHz
Total memory	24 GB	8 GB	24 GB
Memory	-	-	480 GB/s
bandwidth

Equations

$$\mathcal{H}=-\sum _{<\text{mathitbf}r,\text{mathitbf}{r}^{\prime }>}\text{mathitbf}{S}_{\text{mathitbf}r}\cdot \text{mathitbf}{S}_{\text{mathitbf}{r}^{\prime }}.$$

$$C_{\mathcal{O}}(t)\equiv ⟨\mathcal{O}(s)\mathcal{O}(t+s)⟩-⟨\mathcal{O}(t){⟩}^2,$$

$${\rho }_{\mathcal{O}}(t)\equiv C_{\mathcal{O}}(t)/C_{\mathcal{O}}(0).$$

$${\tau }_{\mathrm{int},\mathcal{O}}=\frac{1}{2}+\sum _{t=0}^{\infty }{\rho }_{\mathcal{O}}(t).$$

$${\tau }_{\mathrm{int},\mathcal{O}}=\frac{1}{2}+\sum _{t=0}^{c{\tau }_{\mathrm{int},\mathcal{O}}}{\rho }_{\mathcal{O}}(t),$$

$${\tau }_{\mathrm{int},\mathcal{O}}\sim L^z(1+O(L^{-\omega })).$$

$${\tau }_{\exp ,\mathcal{O}}\equiv \underset{t\to \infty }{\mathrm{limsup}}\frac{-t}{\log {\rho }_{\mathcal{O}}(t)},$$

$$\rho_{\cal O}(t)=t^{-p_{\cal{O}}}f_{\cal{O}}\Big{(}\frac{t}{\tau_{\mathrm{exp},{\cal O}}},\frac{\xi(L)}{L}\Big{)}\,,$$

$$\text{mathitbf}M=\sum _{\text{mathitbf}x}{\text{mathitbf}S}_{\text{mathitbf}x}.$$

$$e(t)=\frac{⟨\mathcal{H}{⟩}_t}{V}.$$

$$e(t)=e_{\infty }+C{t}^{(d-1/\nu )/z}\left(1+A{t}^{-\omega /z}\right),$$

$$C(r,t)=\frac{1}{V}\sum _{\text{mathitbf}x}⟨{\text{mathitbf}S}_{\text{mathitbf}x}{\text{mathitbf}S}_{\text{mathitbf}r+x}{⟩}_t,$$

$$C(r,t)=\frac{1}{r^a}f\left(\frac{r}{\xi (t)}\right)$$

$$C(r,t)\sim \frac{1}{r^{d-2+\eta }}=\frac{1}{r^{1+\eta }},$$

$$I_k(t)={\int }_0^{L/2}dr{r}^{k}C(r,t),$$

$${\xi }_{k,k+1}(t)\equiv \frac{I_{k+1}(t)}{I_k(t)}.$$

$${\xi }_{k,k+1}(t)\sim t^{1/z}\left(1+A_k{t}^{-\omega /z}\right).$$

$$\chi (t)=\frac{1}{V}⟨{\text{mathitbf}M}^2{⟩}_t,$$

$$\chi (t)=\int d^{3}xC(|\text{mathitbf}x|,t).$$

$$\chi (t)=4\pi I_2(t).$$

$$\chi (t)\sim t^{(2-\eta )/z}\left(1+A{t}^{-\omega /z}\right),$$

$$\chi (t)\sim {\xi }_{k,k+1}(t{)}^{2-\eta }\left(1+C_k\xi (t{)}^{-\omega }\right).$$

$$C(r,L)=\frac{A}{r^{1+\eta }}+\frac{A}{(L-r{)}^{1+\eta }}.$$

$$C(r,t)=\frac{a_1}{r^{a_2}}\exp (-a_3{r}^{a_4}).$$

Full text

I Introduction

The study of the dynamics in and out of equilibrium in a critical phase is of paramount importance since it permits to extract the critical exponents of the system, hence, to characterize its universality class. In the last decades a great amount of work, analytical, numerical and experimental, has been devoted to study these issues.

One of the studied systems has been the three dimensional (isotropic) classical Heisenberg model. The dynamic critical exponent, $z$, has been computed using field theory by studying its Model A dynamics (pure relaxational dynamics in which the order parameter is not conserved).Hohen:77 ; Folk:06 ; Tauber:17 A three-loop computation reported in Ref. Antonov:84 provided $z=2.02$. 111 $z=2+c\eta\,\,,\,\,c=0.726-0.137\epsilon+O(\epsilon^{2})\,,$ where $\eta$ is the anomalous dimension of the field (from static) and $\epsilon=4-d$, $d$ being the dimensionality of the model. The equilibrium dynamics of this model was studied by means of numerical simulations in Ref. Peczak:93 and $z=1.96(6)$ was reported. The authors were aware that this exponent was slightly below the analytical computation of Ref. Antonov:84 and discuss in the paper different systematic bias. For example, a relatively narrow range of the lattice sizes and despite the accuracy of their values for the correlation times, a more precise determination of these times were needed to study the corrections to the scaling presented in the model.

From the experimental side, the situation is complicated due to the crossover from the Heisenberg universality class to the dipolar one which induces a change from $z\sim 2.5$ (Heisenberg with conserved magnetization and reversible forces, model J Hohen:77 ; Folk:06 ; Tauber:17 ) to $z\sim 2$ (dipolar). Chow:80 ; Hohenemser:83 In particular using PAC 222Perturbed Angular Correlations of $\gamma$ ray spectroscopy. Hohenemser et al. found $z=2.06(4)$ Hohenemser:83 for Ni and $z\simeq 2$ for Fe; Dunlap et al. Dunlap:80 using ESR 333Electron Spin Resonance. found $z=2.04(7)$ for EuO; and $z=2.09(6)$ was found by Bohn et al for EuS Bohn:84 using inelastic neutron scattering. It seems that the interplay of spin dipoles with orbital angular momentum or dipolar interactions breaks the conservation of the magnetization on these materials, producing a crossover between Heisenberg model J ($z\sim 2.5$) and Heisenberg model A ($z\sim 2$). Chow:80 ; Hohenemser:83

Recently, Pelissetto and Vicari Peli:16 have used the value provided by field theory in the scaling analysis of their numerical data to study the off-equilibrium behavior of three-dimensional $O(N)$ models driven by time-dependent external fields and assigned it an error of 0.01, so $z=2.02(1)$, to take into account the uncertainty on the extrapolation to $\epsilon=1$ of the three-loop-expansion result.

Consequently, it is of paramount importance to obtain an accurate value for this dynamic critical exponent, in order to be used in future numerical analysis and experiments, and also to check the accuracy of the three-loops analytical computation. The main goal of our study is to improve the value of $z$ using numerical simulations by studying the dependence of the integrated correlation time in the equilibrium regime and the behavior of the correlation length, susceptibility and energy with the simulation time in the out of equilibrium region. In this way, we can compare the performance of equilibrium and out-of-equilibrium methods in the computation of the dynamic critical exponent.

Nowadays, a great amount of work has been devoted to study numerically the dynamics of disordered systems, see for example Refs. Parisi:99 ; Hasen:07 ; Janus1 ; Janus2 ; Lulli:16 ; Janus3 ; Janus4 . In general a sudden quenched is performed to work in the off-equilibrium, yet, in other studies the models have been simulated at equilibrium. For example in the three dimensional diluted Ising model both approaches gave the same dynamic critical exponent. Hasen:07

In the equilibrium part of the paper we compute the integrated correlation time, avoiding some of the problems which appear in the computation of the exponential one (e.g. assume that the autocorrelation function is a single exponential function Sokal ; Salas ).

We also study the correlation length in the out-of-equilibrium regime. In the last two decades, this observable has played an important role both in numerical simulations Janus1 ; Janus2 and experiments out of equilibrium Orbach:14 in spin glasses. Due to this, powerful numerical techniques has been developed to compute this observable with high accuracy which has allowed a precise determination of the dynamic critical exponent just at the critical point as well as inside the critical spin glass phase. Janus3 We apply these techniques to the three dimensional (non disordered) Heisenberg model. In addition to the computation of the dynamic critical exponent, we have checked the consistency of previous and very accurate determinations of the static critical exponents ($\nu$ and $\eta$) in the out-of-equilibrium regime. Our starting point will be the (very precise) critical temperature computed in Ref. Balles:96 and the static critical exponents reported in Refs. Campos:02 ; Hasen:11 .

We have also measured the dynamic critical exponent from the decay of the energy at criticality. This decay has also been studied in the past in finite dimensional spin glass Janus1 and recently has played a central role together with the behavior of the correlation length in the analysis of the Mpemba effect, a striking memory effect. Janus4

The structure of the paper is the following. In the next section we introduce the model and the observables. In section III we describe our numerical results: in Sec. III.1 we report the equilibrium determination of $z$ via the integrated autocorrelation time; in Sec. III.2 we study the dependence of the correlation length with time and the computation of $z$ out of equilibrium; in Sec. III.3 and Sec. III.4 the correlation function and the energy have been studied (respectively). Sec. IV is devoted to the conclusions. Two appendices close the paper, one to describe our implementation of GPU and the last one to describe how we have computed the statistical error of the exponents with highly correlated data.

II The model and Observables

The Hamiltonian of the three dimensional Heisenberg model is

[TABLE]

$\mathitbf{S}_{\mathitbf r}$ is a classical three component spin on the site $\mathitbf{r}$ of a three dimensional cubic lattice with volume $V=L^{3}$ and periodic boundary conditions. Without loss of generality we will assume that the spins are unit vectors. The sum runs over all pairs of nearest neighbors spins. We have simulated this model with the standard Metropolis algorithm 444In the equilibrium simulations, Sec. III.1, we have used the Metropolis algorithm proposing a random spin in the unit sphere. In the out of equilibrium simulations we have used the Metropolis algorithm in the standard way AmitMartinMayor : we modify the original spin by adding a random vector and normalizing the final vector to the unit sphere. The magnitude of the random vector is selected in order to maintain an acceptance between 40% and 60%. Both versions of the Metropolis algorithm belong to the same dynamic universality class.Sokal and we have run in CPU (smaller time simulations) and GPU (for larger time simulations). Details on the simulations can be found in the appendix A.

II.1 Equilibrium

We address the problem of the computation of the dynamic critical exponent in the equilibrium regime by means the computation and further analysis of the integrated aucorrelation time as a function of the lattice size.

We compute for a given observable ${\cal O}(t)$, the autocorrelation function (we follow Refs. Madras ; Sokal ; AmitMartinMayor ):

[TABLE]

and its normalized version

[TABLE]

The integrated autocorrelation time is given by

[TABLE]

In a run with $N$ measurements, the number of independent measurements of the observable ${\cal O}$ is just $N/(2\tau_{\mathrm{int},{\cal O}})$. Madras ; Sokal ; AmitMartinMayor If the number of measurements is finite, for large times $t$ the noise will dominate the signal in $\rho_{\cal O}(t)$ and to bypass this problem we use the following self-consistent method to compute the integrated time

[TABLE]

where $c$ is usually taken to be 6 or bigger. Madras ; Sokal ; AmitMartinMayor

At the critical point the integrated autocorrelation time of a long distance observable diverges with the size of the system Hasen:07

[TABLE]

where $z$ is the dynamic critical exponent and $\omega$ is the leading correction-to-scaling exponent (the leading irrelevant eigenvalue of the theory).

Another time, the exponential correlation time, is defined as

[TABLE]

which also depends on the observable used to define $\rho_{\cal O}(t)$. The exponential autocorrelation time controls the approach to the equilibrium.Sokal

Once we have defined the exponential correlation function we can write the general scaling form for the correlation function Sokal

[TABLE]

where $\xi(L)$ is the equilibrium correlation length computed on a system of size $L$. Integrating Eq. (8) in time, we obtain the integrated correlation time and that $\tau_{\mathrm{int}}\sim\tau_{\mathrm{exp}}^{1-p_{\cal{O}}}$. Both times are proportional if and only if $p_{\cal{O}}=0$ and in this situation $z$ is the same for both times. Otherwise, $p_{\cal{O}}\neq 0$, and $\tau_{\mathrm{exp}}$ and $\tau_{\mathrm{int}}$ will provide with different dynamic critical exponents. See also the discussion of Ref. Hasen:07 .

In this paper we will use the slowest mode provided by the non local operator ${\cal O}=\mathitbf{M}^{2}$, where the magnetization ${\mathitbf M}$ is defined as

[TABLE]

II.2 Out of equilibrium

We have focused on only one local observable, the energy, defined as

[TABLE]

We denote the average over different initial conditions at the Monte Carlo time $t$ by $\langle(\cdots)\rangle_{t}$. The renormalization group predicts AmitMartinMayor ; Tauber:17 , at the critical point, the following behavior for this observable:

[TABLE]

where $d$ is the dimensionality of the space (three in this study) and $\nu$ is the critical exponent which controls the divergence of the equilibrium correlation length. The other two exponents $z$ and $\omega$ has been defined in the previous subsection.

One of the main observables on this paper is the correlation function defined as:

[TABLE]

satisfying, at criticality, the following scaling law Tauber:17

[TABLE]

which defines the dynamic correlation length, $\xi(t)$. As we approach the equilibrium regime, $\xi(t)$ reaches its equilibrium value.

At the $d=3$ critical point and in equilibrium, one should expect

[TABLE]

$\eta$ being the anomalous dimension of the field.

The correlation length $\xi(t)$ can be estimated by computing Janus1 ; Janus2

[TABLE]

by means of

[TABLE]

We focus in this work on $\xi_{2,3}$. On spin glasses was measured $\xi_{1,2}$ with a correlation function decaying like $1/r^{0.5}$. Janus1 ; Janus2 In our case, to decrease the weight of the smallest distances we have resorted to compute higher values of $I_{k}$. In the appendix B we describe the detailed procedure we have used to compute the integrals and how we have estimate the statistical error associated with $\xi_{k,k+1}(t)$. The dependence of the dynamic correlation length with time is

[TABLE]

The magnetic susceptibility is given by

[TABLE]

or equivalently by

[TABLE]

In the regime of large $\xi(t)$ we recover rotational invariance and we obtain

[TABLE]

The temporal dependence of $\chi(t)$ is

[TABLE]

which can be rewritten as

[TABLE]

III Numerical results

In this section we report the computation of the integrated correlation time at equilibrium. After this analysis, we describe our results in the out of equilibrium regime. In particular, we consider the short and long time behavior of correlation length and the long time behavior of the correlation function and that of the energy. The data are obtained after a sudden quench from $T=\infty$ to $T=1/\beta_{c}$. All the numerical simulations were performed at $\beta_{c}=0.693001$. Balles:96

III.1 Equilibrium

To obtain the dynamic critical exponent in the equilibrium regime, we compute the integrated correlation time of $\mathitbf{M}^{2}$ when the numerical simulation has reached the equilibrium. We follow the methodology described in Sec. II.1, using the self consistent windowm algorithm with a window size given by $c\tau_{\mathrm{int}}$. We have analyzed the correlation functions with $c=6$, 8, 10 and 12 and we have checked that the $c=10$ data are fully compatible with that of $c=8$ and 12. We report in the following $c=10$ integrated autocorrelation times.

In Table 1 we report the values of $\tau_{\mathrm{int},\mathitbf{M}^{2}}$ and other parameters of the performed runs. In order to improve the statistics on $\tau_{\mathrm{int},\mathitbf{M}^{2}}$ we have performed 50 independent runs (initial conditions). We have computed the statistical error on the integrated autocorrelation times by using the jackknife method over the independent runs. jackknife ; Young

We have fitted $\tau_{\mathrm{int},\mathitbf{M}^{2}}$ to Eq. (6) using $8\leq L\leq 64$ obtaining $z=2.033(5)$ and $\omega=2.7(3)$ with a $\chi^{2}/\mathrm{d.o.f.}=0.36/3$.555In this case, the data for different lattice sizes are not correlated and we can safely use the diagonal covariance matrix. We report this fit and the numerical data in Fig. 1. Fitting the data using only a power law (i.e. neglecting the correction-to-scaling term) we obtain a good fit only for $L\geq 24$ obtaining $z=2.026(4)$ with $\chi^{2}/\mathrm{d.o.f.}=0.28/2$. Both reported values are fully compatible.

Finally we present a scaling analysis of the $\rho(t)$ function at equilibrium and at the critical point to show that $\tau_{\mathrm{exp}}$ and $\tau_{\mathrm{int}}$ are proportional and therefore both times diverge with the same dynamic critical exponent $z$. Fig. 2 shows the scaling law of the correlation function $\rho(t)$ as a function of $t/\tau_{\mathrm{int}}(L)$ instead of $t/\tau_{\mathrm{exp}}(L)$, as stated in Eq. (8). Scaling in the new variable holds if and only if $\tau_{\mathrm{exp}}\propto\tau_{\mathrm{int}}$ and this is the case apart from small scaling corrections on the $L\leq 16$ data induced by the term $\xi(L)/L$ in the scaling function $f_{\cal{O}}$ of Eq. (8), see also the inset of this figure for a detailed and more quantitative view of this effect.

In Ref. Peczak:93 the (biggest) exponential correlation time was computed for the magnetization $\sqrt{\mathitbf{M}^{2}}$. However, the computation of this exponential time is very involved in the case the autocorrelation function $\rho(t)$ does not show a single exponential decay  Salas ; Peli:16 : in Fig. 2 we have plotted a single exponential decay and the correlation function clearly departs from this behavior.

III.2 Correlation length: Shorter times

We report in Fig. 3 the behavior of $\xi_{23}$ as a function of time for different lattice sizes that we have been able to thermalize. The $\xi_{23}$-plateaus obtained for the largest times are a clear evidence that the numerical simulations have reached the equilibrium.

In order to extract the dynamic critical exponent by using Eq. (17), we need to work in the out of equilibrium regime, avoiding the transient regime and the equilibrium domain. Therefore, we need to check the following points:

• •

We need to avoid the transient regime between the power law behavior (pure out of equilibrium regime) and the plateau one (equilibrium one).

• •

Eq. (17) holds after the initial transient of the dynamics. Therefore, we need to fix a minimum time $t_{\mathrm{min}}$.

• •

Finally, in order to avoid finite size effects, we need to compare the correlation lengths for different lattice sizes.

In this section we have performed numerical simulations with $t<10240$ (in order to avoid the transient and equilibrium regimes). We have simulated 4000 random initial conditions for $L=128$ and $L=200$, and 5325 initial conditions for $L=250$.

In Fig. 4 we have plotted, to check finite size effects, the differences among the correlation lengths of $L=128$ and $200$ and that of $250$. In this figure we can see that the data of the $L=200$ and $L=250$ lattices are compatible in the statistical error for $t<8100$.

From the previous discussion we must fit the data for $\xi_{23}$ in the time interval given by $t\in[t_{\mathrm{min}},8100]$ using the $L=250$ data: $t_{\mathrm{min}}$ being the smallest value of the Monte Carlo time that provides a good $\chi^{2}/\mathrm{d.o.f.}$ (e.g. $\sim 1$) by fitting the data to Eq. (17).

In Fig. 5 we show the behavior of $\xi_{23}$ for the largest lattice size simulated in this time regime, $L=250$. By fitting $L=250$ data in the interval $t\in[100,8100)$ we have obtained $z=2.04(2)$ and $\omega=2.2(4)$ with $\chi^{2}/\mathrm{d.o.f.}=287/796$. Furthermore, we can report that a fit neglecting the contribution of the correction-to-scaling term provides $z=2.012(13)$ with $t\in[400,8100)$ and $\chi^{2}/\mathrm{d.o.f.}=361/768$. All two reported values are statistically compatible.

We have computed the statistical error on the $z$-exponent by means of the jackknife method. jackknife ; Young As described in the appendix B, we compute the $\chi^{2}$ using a diagonal covariance matrix (neglecting the correlations of the data), but we use a jackknife procedure to take into account the (important) different correlations among the data. Hence, in the following all $\chi^{2}$ are computed assuming a diagonal covariance matrix; we refer the reader to the appendix B for a discussion of the interpretation of this diagonal $\chi^{2}$ and for more details on the procedure we have followed to take into account the correlation among the data (in time or in distance, see below) and the way we have computed the statistical errors on the values of the critical exponents. 666The same fit performed with the help of Gnuplot Gnuplot (with a diagonal covariance matrix) provides a $z=2.012$ with an asymptotic error of $0.00075$. In order to obtain the right statistical error, we need to divide this asymptotic error by $\sqrt{\chi^{2}/\mathrm{d.o.f.}}=0.687$ Young , obtaining the final value of $z=2.012(1)$. Notice that the computed error discarding correlations among the different times is a factor 13 times smaller.

Having computed $\xi_{23}$ and $I_{2}\propto\xi_{23}^{2-\eta}$ we can, as a check, estimate the $\eta$ exponent. Fig. 6 shows $I_{2}$ as a function of $\xi_{23}$. We can compute $\eta$ using the time interval $t\in[100,8100)$ obtaining $\eta=0.029(20)$ and $\omega=0.8(4)$ with $\chi^{2}/\mathrm{d.o.f.}=342/796$. In this case the $\omega$-exponent is similar to the equilibrium value $\omega\simeq 0.78$. Guida:98 ; Hasen:01 We can improve the value of $\eta$ by fixing the correction-to-scaling exponent $\omega$ to the equilibrium value $\omega=0.78$, obtaining $\eta=0.044(7)$ with $\chi^{2}/\mathrm{d.o.f.}=666/757$ ($t\in[50,8100)$). Our value compares very well (but with 20 times more error) with that computed at equilibrium: $\eta=0.0378(3)$. Campos:02 ; Hasen:11 Finally, neglecting the correction-to-scaling term, we have found $\eta=0.043(6)$ with $\chi^{2}/\mathrm{d.o.f.}=676/768$ ($t\in[400,8100)$). All three reported values are statistically compatible.

III.3 Correlation function for larger times

In Fig. 7 we plot $C(r,t)$ for different times using $L=128$ data (200 initial conditions) and very long times. One can see the crossover of the dynamic correlation function between the off-equilibrium regime and the equilibrium one. In appendix B we provide more details about the functional form of $C(r,t)$ in the out of equilibrium regime. We can also check that we have reached the equilibrium regime by plotting the behavior of $\xi_{23}(t)$ (see $L=128$ curve of Fig. 3). This non-local observable has clearly reached its equilibrium (plateau) value. We can safely assume that for $t>4\times 10^{5}$ we have thermalized the $L=128$ lattice and we can try to extract the value of the the anomalous dimension by averaging the correlation function above this time.

The analytical behavior at the critical point in this regime (large $L$) is given by Eq. (14). Having in mind that we are using periodic boundary conditions, we can write the following improved equation to fit our numerical data

[TABLE]

By fitting the data of Fig. 8 to this functional form, we obtain $\eta=0.026(4)$ (by using only $t>4\times 10^{5}$, $r\geq 16$ and $\chi^{2}/\mathrm{d.o.f}=44/48$) in a good statistical agreement with the value drawn from equilibrium studies $\eta=0.0378(3)$. We have followed the method described in appendix B in order to obtain the error in the $\eta$ exponent. 777The same fit, assuming no correlation among the different values of the correlation function, provides an error of 0.0013, three times smaller than that obtained in our procedure.

III.4 Energy for larger times

We have analyzed the behavior of the energy at criticality in order to compute the ratio of critical exponents $(d-1/\nu)/z$. To analyze this behavior, we have run $L=128$ (153 initial conditions, i.c. in the following), $L=160$ (600 i.c.), $L=200$ (684 i.c.) and $L=250$ (684 i.c.) for longer times $t<102400$.

Firstly, we study in Fig. 9 the effect of a finite size lattice on the values of energy as a function of time. From this figure one can see that it is safe to take fits only in the range $t<48000$ in order to avoid finite size effects (at least in the precision of our simulation).

In Fig. 10 we show the results for the largest lattice $L=250$. We have fitted the $L=250$ data to a power law, in the time interval $t\in[1000,48000)$ obtaining $z=2.034(22)$ and $e_{\infty}=-0.989505(17)$, with a diagonal $\chi^{2}/\mathrm{d.o.f.}=985/939$. We have fixed in the fit the value $\nu=0.7117(5)$. Campos:02 ; Hasen:11 The really small error bar of the $\nu$ exponent has not a measurable effect in the final error bar of $z$. To finish the analysis of the energy, we have also checked corrections to scaling for this observable and we have found that the exponent $\omega_{\mathrm{eff}}=2\times 0.78$ describes very well the numerical data obtaining $z=2.13(7)$ and $e_{\infty}=-0.989525(22)$ with $\chi^{2}/\mathrm{d.o.f}=980/947$.

IV Conclusions

By performing in and out equilibrium numerical simulations we have computed the dynamic critical exponent $z$.

The most accurate value has been computed in the equilibrium regime by studying the integrated correlation time as a function of the lattice size: $z=2.033(5)$. We have found a correction-to-scaling exponent $\omega=2.7(3)$. In addition we have provided strong numerical evidences about the proportionality of the integrated and exponential correlation times.

We have also computed the $z$ exponent in the out-of-equilibrium regime obtaining $z=2.04(2)$ and $\omega=2.2(4)$ which is similar to the $\omega$-exponent computed at equilibrium. Moreover, we have checked the consistency of the computed critical exponents at equilibrium with the out of equilibrium ones with and without considering corrections to scaling. The (equilibrium) value of $\nu$ provides us, by monitoring the energy, with a another dynamic exponent estimate ($z=2.034(22)$) fully compatible with the previous ones.

Furthermore, our value of $z$ has improved the statistical precision of that computed in numerical simulations performed at equilibrium in relatively small lattices ($L\leq 24$). Peczak:93 Our computed values match very well with that obtained in experiments and with the exponent computed using field theoretical techniques Antonov:84 (although in this framework it is very difficult to assign an uncertainty to this estimate).

Acknowledgements.

We thank L. A. Fernandez, M. Lulli, A. Pelissetto, V. Martin-Mayor, J. Salas, J. A. del Toro and D. Yllanes for discussions. This work was partially supported by Ministerio de Economía y Competitividad (Spain) through Grant No. FIS2016-76359-P, by Junta de Extremadura (Spain) through Grant No. GRU10158 and IB16013 (partially funded by FEDER). We have run the simulations in the computing facilities of the Instituto de Computación Científica Avanzada (ICCAEx) and in the CETA-Ciemat thanking Dr. A. Paz for his support.

Appendix A Details of the numerical simulations and GPU parallel implementation

We have simulated the Heisenberg model using the Metropolis Algorithm on CPUs and GPUs (see Table 2). We have simulated $L=128$, 160, 200 and 250 for more than 10000 random initial conditions. The GPU code has been programmed in CUDA C. nvidia11a The original C code which simulates the Heisenberg model has been parallelized in three parts:

Computation of the nearest neighbors of each spin: the C code has a loop which goes sequentially through all the spins one by one. However, in the GPU code each spin has associated an execution thread and all the nearest neighbors of every spin are computed at once. 2. 2.

Metropolis Algorithm: in the sequential C code we can find several loops in the Metropolis part. So, the parallel GPU code reduces meaningfully the execution time especially in large systems ($L\sim 200$). Moreover, the lattice has been divided using a checkboard scheme (Fig. 11). Lulli:15a In this way, the Metropolis algorithm has been executed first of all in the “white” spins and after that in the “black” ones. 3. 3.

Random numbers: to have high quality random numbers is mandatory in Computational Physics. Initially, we have used the CURAND random numbers which are part of the CUDA C distribution.  nvidia11a The problems with the CURAND random numbers have appeared when we have performed long simulations using a huge quantity of random numbers. To avoid these problems we have used Congruential Random Numbers. Knuth:98a

Making use of the GPU Tesla K80 we have achieved a speedup of 22 which represents an important reduction of the execution time.

Appendix B Details of the analysis of the computation of the correlation length

We describe the different steps we have followed in order to compute $\xi(t)$ and its associated exponent $z$. Janus1 ; Janus2 ; ThesisDavid ; Michael:94 The important point of this approach is to avoid the use of the full covariance matrix since this matrix is frequently singular (see for example ThesisDavid ; Seibert:94 ). Thus, our procedure is the following:

We compute using the jackknife method over the set of the initial conditions, the statistical error of $C(r,t)$, denoted as $\sigma[C(\Lambda,t)]$. 2. 2.

To compute $I_{k}$ we introduce a cutoff to have a good control of the signal to noise ratio of $C(r,t)$ for large values of $r$ (see also Fig. 7).

• •

We compute the cutoff $\Lambda$ using the condition $\sigma[C(\Lambda,t)]=4C(\Lambda,t)$.

• •

For a fixed $t$ and $r_{\mathrm{min}}<r<\Lambda$ we fit the correlation function to the functional form given by

[TABLE]

with $r_{\mathrm{min}}$ is the minimum value of $r$ which provided, for $C(r,t)$, a good fit (e.g. $\chi^{2}/\mathrm{d.o.f.}\sim 1$) to Eq.(24). In Fig. 12 we report the dependence of the exponents $a_{2}$ and $a_{4}$ with the Monte Carlo time. Notice that $a_{2}$ converges to the equilibrium value (see Eq. (14)) given by $1+\eta=1.0378$ and $a_{4}\simeq 1.8$.

• •

We compute the integral in Eq. (15) using the numerical values of the correlation $C(r,t)$ for $r<\Lambda$ and using the values provided by the fit (Eq. (24)) for $\Lambda<r<L/2$.

• •

Using the previous procedure, we compute the statistical error of $\xi(t)$ using again the jackknife method over the set of the initial conditions. The time interval for the fit is decided by imposing a diagonal $\chi^{2}/\mathrm{d.o.f}\sim 1$.

• •

The jackknifed $\xi$’s are used to compute the jackknifed values of $z$ and this allows us to compute the statistical error of the dynamic critical exponent using the standard deviation in the jackknife method. Notice that for extracting $z$ on each jackknife block, we use the diagonal covariance matrix. However, the jackknife procedure reproduces with high accuracy the effect of the correlations among the different times.

Notice that the diagonal $\chi^{2}/\mathrm{d.o.f.}$ has not a rigorous interpretation as that of the full (non diagonal) one. One can show (see the detailed analysis of this procedure carried out in section B.3.3.1 of Ref.  ThesisDavid ) that the diagonal $\chi^{2}/\mathrm{d.o.f.}$ behaves as if there were a small number of degrees of freedom, hence, one can not compute confident limits as usual.

Finally, in Ref. Lulli:16 was shown that the error bars are essentially equal (using this jackknife procedure, neglecting the correlations among the data) to those obtained taking into account all the statistical correlations among the data.