# Universal long-time relaxation on lattices of classical spins: Markovian behavior on non-Markovian timescales

## 1 Introduction

We perform a numerical study of the long-time behavior of infinite temperature correlation functions defined on an infinite lattice of classical spins as:

$$G(t)=S_k^x(t)_n\text{cos(}𝒒\mathbf{}𝒓\text{kn}\text{)}S_n^x(0),$$
(1)   
where $`S_k^\mu `$ is the $`\mu `$th ($`x`$,$`y`$ or $`z`$) spin component on the $`k`$th lattice site; $`𝒓`$<sub>kn</sub> is the translation vector between the $`k`$th and the $`n`$th sites; and $`𝒒`$ is a wave vector commensurate with the lattice periodicity. We consider three types of lattices: a simple one-dimensional chain, a two-dimensional square lattice, and a three-dimensional cubic lattice. In each case, the dynamical evolution of the system is driven by the nearest-neighbor interaction represented by the Hamiltonian

$$=\underset{k,n}{}[J_xS_k^xS_n^x+J_yS_k^yS_n^y+J_zS_k^zS_n^z],$$
(2)   
where $`J_\mu `$ are coupling constants. With such a Hamiltonian, the timescale of the individual spin motion referred to below as the “mean free time” can be given by the time $`\tau =[\frac{1}{3}NS^2(J_{x}^{}{}_{}{}^{2}+J_{y}^{}{}_{}{}^{2}+J_{z}^{}{}_{}{}^{2})]^{1/2}`$, where $`N`$ is the number of the nearest neighbors (twice the number of the lattice dimensions).

In the context of inelastic neutron scattering, correlation functions (1) are called “intermediate structure factors”. If $`q=0`$, Eq.(1) can also represent the free induction decay in nuclear magnetic resonance (NMR).

In this work, we provide extensive numerical evidence that the generic long-time behavior of $`G(t)`$ has one of the following two functional forms: either

$$G(t)e^{\xi t},$$
(3)   
or

$$G(t)e^{\xi t}cos(\eta t+\varphi ),$$
(4)   
where the constants $`\xi `$ and $`\eta `$ are of the order of $`1/\tau `$.

It is important to realize that, if the above functional form of the long-time behavior is, indeed, generic (i.e. independent of the specific details of interaction), then this property is very likely related to the randomness generated by the spin dynamics. At the same time, the problem cannot be reduced to the Markovian paradigm of “a slow variable interacting with a fast equilibrating background” — the characteristic timescale $`\tau `$ in Eqs.(3,4) is not ”slow”. It is, in fact, the fastest natural timescale of the problem. Therefore, whatever is the ultimate explanation of that behavior, it will certainly be a step beyond the standard theory of Brownian-type motion.

Our interest in the long-time behavior of the correlation functions (1) was originally motivated by two isolated pieces of evidence supporting the oscillatory behavior (4) in quantum (spin 1/2) systems: (i) experiments on NMR free induction decay in CaF<sub>2</sub>  and (ii) the results of numerical diagonalization of spin 1/2 chains . In the both cases, quantities analogous to the one defined by Eq.(1) have been measured or computed, and the results look very similar to the plots shown in the left column of Fig. 1.

We came to recognize the importance of a detailed study of the classical limit, when, in an attempt to explain the long-time relaxation in quantum spin systems, we developed a theory that turned out to be simultaneously applicable to classical spins. That theory is presented in a different paper which has been written simultaneously with the present one. The present paper is mainly numerical: it is not intended to be a brief exposition of Ref.. Below, we only provide the summary of the results from Ref..

The theory developed in Ref. describes long-time relaxation as a correlated diffusion in finite volumes. In the classical case, those finite volumes correspond to the spherical surfaces on which the tips of classical spin vectors move, while in the quantum case, the finite volumes originate from a more sophisticated construction in Hilbert space. The overall structure of such a treatment has noticeable parallels with the theory of Pollicott-Ruelle resonances in classical chaotic systems. A definite prediction from the correlated diffusion description is that the functional form of the long-time relaxation should be given by Eqs.(3,4).

The important part of the above theory is not the diffusion description itself but the reason why it is applicable, given the “non-Markovian” relaxation timescale. The theory is based on the fairly strong conjecture that, for a broad class of many-body systems, a formal extension of the Brownian-like description applies to the the long-time behavior of the ensemble average quantities, even when the problem exhibits no separation of the timescales between the slow and the fast motions.

Before proceeding with the description of the simulations, it should be mentioned that, for the classical spin systems at infinite temperature, the long-time behavior of the $`q`$-dependent correlation functions (1) decaying on the timescale of $`\tau `$ has never been addressed. The closest to this subject was the work of de Alcantara Bonfim and Reiter, who considered the Heisenberg spin chain and focused on the long-time behavior of correlation functions (1) with small $`q`$. Those correlation functions, however, are not typical for our purposes, because, as a consequence of the total spin conservation, they decay on the timescale, which is much longer than the characteristic timescale of one-spin motion. In that situation, the hypothesis of spin diffusion would lead to the prediction of nearly exponential decay with the decay constant proportional to $`q^2`$. The results of de Alcantara Bonfim and Reiter did not cover the range of values, which would be sufficient to confirm or rule out the exponential character of the long-time decay. However, those results (in line with others,) indicated that, if, the spin diffusion regime exists for classical spin chains, the approach to that regime is anomalously slow.

The present work includes one example of the Heisenberg interaction, just to show that this case does not appear to be special with respect to the long-time property (3,4).

## 2 Simulations

Our computational strategy was similar to that of Müller. Namely, we did not deal with very large systems but, instead, performed an ensemble averaging over a large number of finite, but not too small, lattices having periodic boundary conditions. The finite size effects were then controlled by varying the size of the lattice.

For a given lattice size, many computational runs have thus been performed. Each of them started from completely random initial conditions (corresponding to the infinite temperature) and generated the evolution of the system over a time interval two orders of magnitude longer than the mean free time $`\tau `$ (see Table 1 for specific numbers). The correlation functions were then obtained by averaging the data within each run and over different runs.

The following algorithm has been used in order to simulate the evolution of the system.

At each time step, the spins were advanced sequentially in such a way that, if the spin number $`k`$ interacted with the spin number $`n`$, and the $`k`$th spin was advanced first, then the new coordinates of the $`n`$th spin were computed based on the local field created by the already advanced $`k`$th spin.

The procedure for advancing a given ($`k`$th) spin to the next point along the discrete time grid consisted of two steps.

Step 1: The coordinates of the $`k`$th spin were changed by $`\delta `$$`𝑺`$<sub>k</sub> according to the straightforward discretization of the equations of motion, i.e.

$$\delta 𝑺_k=[𝑺_k\times 𝒉_k]\delta t,$$
(5)   
where $`𝒉`$<sub>k</sub> was the local field equal to $`_n^{\text{n.n.}}[J_xS_n^x\widehat{𝒆}_x+J_yS_n^y\widehat{𝒆}_y+J_zS_n^z\widehat{𝒆}_z]`$ (sum over the nearest neighbors of the $`k`$th spin); and $`\delta t`$ was the discretization time step.

Step 2: The higher order errors that changed the length of the spin vector were eliminated. This was done by contracting the spin component perpendicular to the local field, so that it took the absolute value it had before Step 1.

The whole manipulation could not change the spin projection parallel to the local field and, therefore, the energy of interaction of that spin with its neighbors. Since the next spin was advanced in the newly updated local field, the energy of the whole system was conserved exactly during the entire integration. Apart from insuring meaningful behavior of the computed trajectories, the exact conservation of energy substantially improved the convergence of the algorithm with respect to the limit $`\delta t0`$.

In our simulations, the discretization time steps (given in Table 1) were admitted as sufficient when their further reduction appeared to have no effect on the computed correlation functions. We also checked that the averaging over a larger number of much shorter runs led to results consistent with the longer runs we used.

The time lengths of the runs indicated in Table 1 were chosen to optimize the resulting efficiency of averaging: too short runs (of the order of the time length of the computed correlation function) did not make many independent contributions to the correlation functions, while too long runs did not improve the quality of the averaging proportionally to their length.

## 3 Results

The results of our simulations for different dimensions, interaction constants, and wave numbers are presented in Fig. 1, together with the long-time theoretical fits based on either Eq.(3) or (4). Since we aimed at demonstrating the exponential character of the long-time behavior, it was natural to use a logarithmic scale for $`G(t)`$. However, because in the half of the cases, $`G(t)`$ was also oscillating, we chose to show the logarithmic plots for the absolute value of $`G(t)`$, which explains the cusps in the left column of plots in Fig. 1. Those cusps correspond to the points where $`G(t)`$ crosses zero. The reason that the cusp minima do not reach $`\mathrm{}`$ is that a discrete grid was used.

The selection of parameters for the simulations was subject to certain practical constraints: The computed correlation functions could be considered reliable only within quite a limited range along both the $`t`$\- and Log($`G(t)`$)- axes. Therefore, we avoided the correlation functions that reached too-small values too fast, or, on the contrary, decayed too slowly.

From our experience, the finite size effects were least pronounced for correlation functions with $`q=0`$. For this reason, $`q`$ was chosen to be zero in six of eight examples presented in Fig. 1.

Since the long-time behavior of the correlation functions was only marginally accessible with our computational resources, we present the results in substantial detail, thus making clear the uncertainties associated with insufficient ensemble averaging and finite size effects. Each of the correlation functions presented in Fig. 1 was computed four times: two statistically independent averaging results for each of two different lattice sizes. “Two statistically independent averaging results” means that, in each case, the same number of sample runs was performed but two different sets of random numbers were used for setting the initial orientations of spins. Each frame in that figure thus contains a superposition of four plots. The time interval where these plots do not deviate from each other can be considered as representing the limit of infinite lattice size with sufficient ensemble averaging.

The finite size effects are not evident in any of the examples shown in Fig. 1 — in every case, the plots representing different simulation outcomes do not deviate from each other before the statistical fluctuations for each of the two lattice sizes become apparent.

Our experience indicates that further improvement in the accuracy of the computed correlation functions would simultaneously require a finer discretization, much more extensive ensemble averaging and, probably, larger system sizes, i.e. much greater computational effort.

Summarizing the evidence, we observe that, with the marginal exception of Fig. 1(d), in every other case presented, there is an interval, covering at least one decade of the values of $`G(t)`$, where the simulation results agree with the long-time fits (3) or (4). We would also like to point out that, while in all cases the exponential behavior becomes pronounced quite early, the long-time exponents in Figs. 1(e,g) describe almost the entire correlation functions.

Thus our numerical results lend strong support to the idea expressed in Ref., that a discrete spectrum of well-separated exponents describes the long-time behavior of the correlation functions considered, with the slowest of those exponents responsible for the asymptotic functional form given by Eq.(3) or (4). It is also very likely, though slightly less reliable, that in the examples presented, our simulations revealed the slowest exponents. The reservation here is for the possibility that even slower exponents could enter the long-time expansion of $`G(t)`$ with anomalously small coefficients.

## 4 Conclusions

In conclusion, we have presented a numerical evidence that the long-time behavior of the correlation functions considered is exponential with or without the oscillatory component. Leaving the details to Ref., here we just mention that our best hope for the theoretical explanation of that behavior is associated with the strong chaotic properties of the spin dynamics. Those properties are likely to be quite generic, i.e. present in other systems.

## 5 Acknowledgements

The author is grateful to A. J. Leggett for numerous discussions, and to C. M. Elliott and R. Ramazashvili for helpful comments on the manuscript. This work has been supported in part by the John D. and Catherine T. MacArthur Foundation at the University of Illinois, by the Drickamer Endowment Fund at the University of Illinois and by the Foundation of Fundamental Research on Matter (FOM), which is sponsored by the Netherlands Organization for the Advancement of Pure Research (NWO).

Note added to proof: After submitting this paper, the author became aware of a closely related work of T. Prosen, in which the decay of the correlation functions defined on the lattices of quantum spins was interpreted in terms of Policott-Ruelle resonances.

TABLE CAPTIONS:

Table 1:

Simulation parameters and the long-time fits corresponding to the plots presented in Fig. 1. The numbers in the right four columns characterize each of the four data sets superimposed in the corresponding frame.

FIGURE CAPTIONS:

Figure 1:

Correlation functions $`G(t)`$ of the form (1) for 1D chain, 2D square lattice, and 3D cubic lattice. The interaction coefficients and the wave numbers (in the units of inverse lattice spacing) are indicated above the plots. The main frame of each figure shows the logarithmic scale of the absolute value of $`G(t)/G(0)`$, while the inset frame shows the direct plots of $`G(t)/G(0)`$ (with neither the logarithm nor the absolute value being taken). Within each frame, the simulation results are presented by almost indistinguishable superposition of data for two lattice sizes, and each size is represented by two statistically independent averaging results; therefore, two solid lines for the larger size (Size 1) and two dash-dotted lines for the smaller size (Size 2). The spread of the four lines indicates the computational uncertainty — it becomes visible only in the lower right corner of each frame. The dashed lines in each figure are the long-time theoretical fits of form (3) or (4). The numbers relevant to each data set are given in Table 1.