Anomalous Dynamics and Equilibration in the Classical Heisenberg Chain

The search for departures from standard hydrodynamics in many-body systems has yielded a number of promising leads, especially in low dimension. Here we study one of the simplest classical interacting lattice models, the nearest-neighbour Heisenberg chain, with temperature as tuning parameter. Our numerics expose strikingly different spin dynamics between the antiferromagnet, where it is largely diffusive, and the ferromagnet, where we observe strong evidence either of spin super-diffusion or an extremely slow crossover to diffusion. This difference also governs the equilibration after a quench, and, remarkably, is apparent even at very high temperatures.

The search for departures from standard hydrodynamics in many-body systems has yielded a number of promising leads, especially in low dimension. Here we study one of the simplest classical interacting lattice models, the nearest-neighbour Heisenberg chain, with temperature as tuning parameter. Our numerics expose strikingly different spin dynamics between the antiferromagnet, where it is largely diffusive, and the ferromagnet, where we observe strong evidence either of spin superdiffusion or an extremely slow crossover to diffusion. This difference also governs the equilibration after a quench, and, remarkably, is apparent even at very high temperatures.
The focus of this work is the classical Heisenberg spin chain, for which the nature of hydrodynamics has provoked extensive debate. Based on the lack of integrability, it has been argued that ordinary diffusion holds for both spin and energy [55][56][57][58][59][60][61]. However, there have also been proposals of anomalous behaviour [62][63][64][65], including an argument for logarithmically enhanced diffusion [65]. Ref. [61], in contrast, has argued from a theory of non-abelian hydrodynamics that each component of the spin follows a separate, ordinary diffusion equation.
In this paper, we present a systematic numerical study of the dynamical correlations and equilibration dynamics over a wide range of temperatures, T < |J| to T = ∞. We find ordinary diffusion of both spin and energy at T = ∞ and ordinary diffusion of energy at all (nonzero) temperatures in both the ferromagnetic (FM) and antiferromagnetic (AFM) chains [66]. Most strikingly, we find a qualitative difference between ferromagnetic and antiferromagnetic models at finite temperatures. This manifests as a temperature-dependent finite-time dynamical exponent in the spin correlations of the ferromagnetic chain, which departs from the diffusive exponent α = 1/2; whereas the antiferromagnetic chain displays behaviour compatible with spin diffusion at all tem-peratures studied. This deviation is apparent even at high temperatures, where the correlation length is still of the order of a single lattice spacing -far from the low-temperature regime where the distinction between quadratic ferro-and linear antiferromagnetic spin-wave spectra may play a role. We have thus identified a, hitherto perhaps unappreciated, fundamental difference between the dynamics of the FM and AFM models.
The observed behaviour of the ferromagnet could be interpreted as anomalous diffusion with a temperaturedependent exponent, or alternatively as a crossover at remarkably large timescales, rendering diffusion in practice unobservable experimentally for a wide range of temperatures. Intriguingly, at low temperatures where we obtain the best fit to a single power-law, we observe the KPZ exponent almost perfectly across three decades in time. In addition, the spacetime profiles of correlation functions closely follow the KPZ scaling form. This establishes intermediate time KPZ scaling at low temperatures in the FM Heisenberg model, even if ultimately followed by a crossover to normal diffusion at very long times.
As a related phenomenon, we study equilibration dynamics after quenches from an XY to a Heisenberg chain. Equilibration is shown to proceed via a power-law approach to the equilibrium value, with an exponent determined by that observed in the corresponding unequaltime equilibrium correlation function, again displaying anomalous finite-time exponents in the case of the FM.
Model.-We consider the periodic-boundary classical Heisenberg chain, with Hamiltonian for unit length classical spins S i ∈ S 2 . Here J = 1 for the FM chain, and J = −1 for the AFM chain. The dynamics is given by the classical Landau-Lifshitz equation of motion,Ṡ which we solve numerically [66]. In equilibrium, we probe the spin-spin correlations and the energy correlations where E j = −JS j · S j+1 is the bond energy, and E = E is the internal energy density. We use internal energy and temperature interchangeably, via E(T ) = T − coth(1/T ) [66,67]. Also, the (equal-time) spin correlation length is which, as a function of E, is the same for the Heisenberg (1) and XY chains (11) [66]. Both of these correlation functions are symmetric under parity and time-reversal. To evaluate these correlations for a given E, we first construct an ensemble of 20,000 initial states drawn from the canonical ensemble of H at the temperature T (E) [66,68]. Each state is evolved in time, cf. (2), with snapshots stored at intervals of ∆t = 10J −1 . The correlation function at a fixed time difference t is calculated by averaging over 1000 consecutive snapshots for every initial state.
Hydrodynamics and Scaling Functions.-The hydrodynamic theory posits an asymptotic scaling form for the correlations of the conserved densities, with a scaling exponent α and universal function F. The exponent is, in principle, independent of the precise form of F, and may be extracted by fitting the autocorrelation function, A(t) = C(0, t), to a power-law, Ordinary diffusion corresponds to an exponent of α = 1/2, and a Gaussian scaling function, where χ = dx C(x, t) = j C(j, t), and D is the diffusion constant. This scaling function may be obtained directly by solving the ordinary diffusion equation.
The most well-known anomalous scaling is the KPZ universality class, with an exponent α = 2/3. There is no analytic form for the scaling function, but it is tabulated in [69].
Even if the the asymptotic behaviour is diffusive, one might have finite-time corrections. The lowest order correction to a diffusive autocorrelation function from nonabelian hydrodynamics [61] is of the form The AFM displays ordinary spin diffusion, with the diffusive power-law observable after a comparatively short time t ≈ 10 3 . The FM does not exhibit diffusion, at any finite temperature, over the timescales of our simulations.
The autocorrelations of the FM are, for these timescales, well-approximated by a power-law (7), with superdiffusive exponents ( Fig. 1(b)). One may also fit a crossover of the form (9). Adopting this point of view, we may extract a crossover time t × (E) after which we would predict the system to show diffusion, via the effective exponent with the crossover defined, arbitrarily, by α eff (t × ) = 0.505 (i.e., the time after which α eff is sufficiently close to 1/2). The estimated crossover times obtained for the FM are orders of magnitude larger than the AFM, reaching t × ≈ 10 7 J −1 at low-to-intermediate temperatures ( Fig. 1(d)).
In Fig. 2 we show the scaling collapses at these temperatures. The AFM is clearly consistent with a diffusive collapse (8); the FM is not. Since the autocorrelation function is well-fit by an anomalous power-law (7), we use these exponents to perform the scaling collapse in the FM. This collapses the correlations rather well, though there is some noise in the tails.
Moreover, despite the noise, one may observe that the tails of the correlations decay faster than a Gaussian. This suggests that an enhancement of the diffusion constant alone (whether of the form (7), or a crossover (9)) is not the correct picture. KPZ Scaling.-We thus examine the possibility that we are observing a crossover from KPZ scaling. Indeed, the numerical evidence at low temperature is remarkably strong, shown in Fig. 3. The correlations up to t = 10 4 collapse onto the KPZ function for E = −0.8 and E = −0.9. Beyond this time the noise in the tails is too great to reliably distinguish the form of the spatial decay, but the scaling exponent measured by the autocorrelation function is consistent with α = 2/3 up to the final time t = 10 5 . Moreover, at E = −0.8, there are apparently no finite-time corrections to the power-law decay A(t) ∼ t −2/3 , for three decades in time.
Equilibration Dynamics.-In addition to our equilibrium simulations examining unequal-time correlations, we perform equilibration simulations probing the relaxation to thermal equilibrium after a quench. This allows us to test whether the anomalous behaviour, and the distinction between FM and AFM, are also observable in out-of-equilibrium dynamics.
We initially prepare the system in a thermal state of the XY chain, (11) for unit length classical rotors S i ∈ S 1 . At time t = 0, we quench the system, and evolve under the dynamics (2) of the Heisenberg chain. We examine the relaxation of the following observables: which measures the energy attributed to the µ th spin components; and which measures the total magnitude of the µ th spin components. These are natural measures of the anisotropy, which characterises the relaxation from the initial state, satisfying S z i = 0 ∀i, to a (quasi-)thermal state of the isotropic Heisenberg chain. The equilibration of the energy fluctuations is measured using the heat capacity, where we take the spatial variance before the ensemble average to obtain a time-dependent quantity. As in the equilibrium simulations, we average over an ensemble of 20,000 states, initially drawn from the canonical ensemble of H XY . We expect that the equilibration dynamics will be similarly hydrodynamic, since establishing the new global equilibrium requires the transport of conserved densities over long distances [1]. The relaxation of an observable O is therefore expected to follow a power-law where O eq is the thermal value of the observable in the Heisenberg chain. These simulations exhibit complementary aspects of the same broad phenomenology observed in equilibrium. Fig. 4 shows the equilibration dynamics at E = −0.5. The extracted (anomalous) equilibration exponents have qualitatively similar dependence on energy as those extracted from equilibrium correlation functions [66]. The energy fluctuations, as measured by the heat capacity, always equilibrate diffusively. In the AFM, E µ and Q µ also show diffusive equilibration. In the FM, however, the equilibration of E µ and Q µ is anomalous. It should be noted that, although E µ has dimensions of energy, it is, like Q µ , a measure of the magnetic anisotropy, and therefore equilibrates anomalously in the FM, rather than tracking the diffusive behaviour of the energy fluctuations.
Thus, as in equilibrium, we observe a striking difference between the FM and AFM, with only the former displaying anomalous exponents.
While our simulations do not allow us to rule out a potential crossover to diffusive equilibration at even longer times, the observables can reasonably be described to have (fully) equilibrated with these anomalous exponents, in particular when considering a realistic experimental situation in which resolution and time scales might be limited.
Discussion & Conclusions.-We have conducted a detailed numerical study of the equilibrium and out-ofequilibrium dynamics of the classical Heisenberg chain, with the largest system sizes, simulation times, and range of temperatures we are aware of so far for this model. We find that, although ordinary diffusion is expected at infinite time, the FM exhibits a long-lived regime that is well-described by an effective superdiffusive, temperature-dependent scaling exponent, with remarkably clean KPZ-like behaviour at low temperature. The AFM, by contrast, swiftly evinces ordinary diffusion at all temperatures. The existence of such large intermediate scales is obviously relevant to experiments probing anomalous diffusion -anomalous behaviour might be the only accessible experimental regime even when the longer-term behaviour is diffusive.
A possible explanation of the intermediate-time regime and the stark difference between FM and AFM cases could be the presence of integrable ferromagnetic models, such as the continuum Landau-Lifshitz model [70][71][72][73] and the lattice FM model with log(1 + S i · S i+1 ) interactions [54,[73][74][75][76][77]. The Heisenberg FM studied in the present work could arguably be considered to be increasingly similar to either of these integrable models at lower temperatures. This would be consistent with our finding that the low-temperature FM behaviour is closer to KPZ. Recent work [12][13][14][15] suggests the perturbative stability of KPZ scaling in systems which are close to integrability and preserve the rotational symmetry.
Nevertheless, it is remarkable that we observe an anomalous regime even at near-infinite temperatures, where the correlation length (5) is short, e.g., already less than a single lattice spacing at E = −0.3. Intuitively, this regime does not seem to be close to either the continuum model or the integrable log-interaction model. This points to the need for a better understanding of how proximity to integrable points might play a role in the physics of the Heisenberg FM, especially at high temperatures.

CONTENTS
In this supplementary we provide further details of our simulations, and further evidence for our conclusions. In S-I we provide the exact thermodynamics of the XY and Heisenberg chains. S-II contains an account of our numerical methods for the construction of thermal states and time evolution. In S-III we show evidence for spin diffusion at infinite temperature, and in S-IV we show the diffusion of energy at all temperatures. Finally, in S-V we provide more information about our equilibration simulations, and provide the extracted anomalous exponents.

S-I. EXACT THERMODYNAMICS
For reference, we provide here the exact thermodynamics of the classical Heisenberg [67] and XY chains. The derivation is given with open boundary conditions, since this is simpler -but the boundary effects vanish in the thermodynamic limit.
Consider first the XY chain. The partition function is where I n denotes a modified Bessel function of the first kind, and in the second line we have transformed the coordinates to ϕ i = φ i − φ i+1 . Similarly, the partition function of the Heisenberg chain is where, in the second line, we rotate the coordinate system of S i+1 such that the z-axis aligns with S i . From these expressions, all of the thermodynamic quantities may be calculated. In particular, taking the thermodynamic limit and setting J = 1, the internal energy density and specific heat are: for the XY chain, and: for the Heisenberg chain. The equal-time two-point correlation function is given by and so the correlation length is which, as a function of internal energy, is the same for both the XY and Heisenberg chains. In the antiferromagnet, the above is replaced with the staggered correlator.

Construction of Thermal States
Our initial thermal states are constructed using a heatbath Monte Carlo method [68], where we use the fact that we can precisely invert the thermal probability distribution for a single spin in a magnetic field.
The general method is known as inverse transform sampling. Let X ∈ [a, b] ⊆ R a random variable on some real interval, with probability density function p(X). The cumulative distribution function (CDF) is i.e., the probability that a randomly sampled X is less than or equal to x. Then the random variable Y = F −1 X (u), where u is uniformly random over [0, 1], has the same probability distribution as the original variable X since, by construction, The specific problem is to invert the CDF. For a Heisenberg spin, this can be done analytically. Letting h = hẑ, the CDF for S z is (S9) Then, using F z (F −1 z (u)) = u, we find that we can sample S z as The S x and S y components have random direction in the plane, and magnitude set by the condition |S| = 1. For a magnetic field of arbitrary direction, one need only appropriately rotate the sampled spin. For XY spins, the inversion cannot be performed analytically. In this case, we let h = hx, and sample the angle φ, where (S x , S y ) = (cos(φ), sin(φ)), by numerically solving the integral equation for randomly generated u. Again, the generalisation to arbitrary magnetic field direction is via a rotation.
To sample a state from the canonical ensemble, we randomly generate the first spin S 1 . We then sweep through the chain, generating S i+1 from the thermal distribution of the effective field h = JS i . For open boundary conditions, cf. the coordinate transformations used to calculate the partition functions, this exactly samples the canonical ensemble. For periodic boundary conditions we perform an additional 1000 sweeps through the chain.

Time Evolution
For the equilibration simulations, we integrate the equations of motion with the standard fourth-order Runge-Kutta (RK4) method, using a fixed timestep of ∆t = 0.002J −1 . This method conserves the magnetisation to machine precision, and for a system size L = 16384 and a final time t f = 4096J −1 the error in the energy density is limited to ∼ 10 −10 .
For the equilibrium simulations, we use a system size of L = 8192, but a much longer final time t f = 1.1×10 5 J −1 . To achieve such times, we use the discrete-time oddeven (DTOE) method, with a larger timestep of ∆t = 0.05J −1 . The method consists of updating the odd spins for a timestep ∆t by exactly solving the equations of motion with the even spins held fixed, and then vice versa. This method conserves the energy to machine precision, but the error in the magnetisation is suppressed only as O(∆t 2 ). However, the error does not grow with time -for the timestep chosen the magnetisation error is ∼ 10 −5 .

S-III. SPIN DIFFUSION AT INFINITE TEMPERATURE
There has been some recent controversy over the nature of the hydrodynamics at T = ∞ in the Heisenberg chain, with [65] claiming logarithmically enhanced diffusion and [61] arguing for ordinary spin diffusion. Here we provide our own contribution to this debate: we see no evidence for logarithmically enhanced diffusion at long times, which would predict t 1/2 A(t) → 0 as t → ∞, see Fig. S1(a).
In Fig. S1(b) we show the diffusive scaling collapse of the spin correlations from t = 2000 to t = 10 5 . Note that the ferromagnet and antiferromagnet are indistinguishable at infinite temperature.

S-IV. ENERGY CORRELATIONS
We have reported in the main text that the energy correlations are found to be diffusive for both the FM and the AFM. We show the evidence for this in Fig. S2, where we plot the Gaussian width of the energy correlations C E (x, t) as a function of time. We find that, except at the lowest temperatures, they are well-fit by the diffusive power-law.
At low temperatures, ballistically propagating spinwave modes persist to intermediate times. This makes the observation of diffusion in our simulations rather difficult -even at longer times where the ballistic modes have decayed -because, over their lifetime, the ballistic modes increase the width of the correlations much faster than diffusion. By the time this effect is negligible, the width is comparable to the system size, and finite size effects take over. We show the ballistic collapse in Fig. S3.  S4. Anomalous equilibration exponents αE and αQ for the observables E µ and Q µ in the FM. The measured equilibrium exponent is also shown for comparison. The discrepancy at low-temperature is probably due to timescales -the equilibrium exponent is extracted over the range t = 10,000 to t = 100,000, whereas the equilibration exponents are extracted up to t = 4096.

S-V. EQUILIBRATION DYNAMICS
We provide here some further details of the equilibration simulations -in particular, how we determine the thermal values, and the temperature dependence of the (finite-time) anomalous exponents.
We begin from a thermal state of the XY chain, with every spin confined to the plane S z = 0, and evolve towards a quasi-thermal state of the Heisenberg chain. Recall from the main text that we measure the degree of anisotropy with the observables E µ (t) = −J S µ i (t)S µ i+1 (t) (S12) and Q µ (t) = S µ i (t) 2 . (S13) For a state with energy density E, these observables are constrained by µ Q µ = 1 and µ E µ = E. Their Heisenberg equilibrium values are thus determined by isotropy, to wit, Q µ eq = 1/3 and E µ eq = E/3. There is a caveat: at finite size there is a small, but conserved, total magnetisation, which prevents Q µ and E µ from attaining their precise equilibrium values. However, this correction may be calculated exactly. Given, at system size L, the q = 0 component of the static structure factor, the asymptotic values are: Q z → 1/3 − ∆/3, Q → 1/3 + ∆/6, E z → E/3 + J∆/3, E → E/3 − J∆/6, (S15) where we have defined the averages of the in-plane observables, Q = (Q x + Q y )/2 and E = (E x + E y )/2. We take the average of the two in-plane components to account for the finite-size magnetisation spontaneously breaking the rotational symmetry of the initial XY state -which means we can read off the initial values exactly from the sum rules.
To examine the equilibration of energy fluctuations, we consider the heat capacity. Recall that the heat capacity can be estimated from a thermal ensemble as where E is the energy of the state and the angle brackets denote the ensemble average. Since the energy of each state is conserved by the Hamiltonian dynamics, this is time-independent. However, we define the heat capacity of a single state as where the variance is taken over the spatial distribution of the energy. In equilibrium, this is equal to the ensemble-based definition (S16), but it is not conserved by the dynamics. The initial value, of course, is the heat capacity (S3) of the XY chain, except that, since the correspondence between internal energy and temperature is different in the two chains, we must multiply the above by (T XY /T H ) 2 , with T XY and T H the temperatures that correspond to E. The equilibrium value C eq is then given by the heat capacity of the Heisenberg chain (S4). As mentioned in the main text, the equilibration simulations probe a different aspect of the underlying phenomenology: Q µ and E µ equilibrate diffusively in the AFM, but anomalously in the FM. The heat capacity always equilibrates diffusively.
The anomalous exponents obtained from the equilibration in the FM are shown in Fig. S4.