Classical many-body chaos with and without quasiparticles

We study correlations, transport and chaos in a Heisenberg magnet as a classical model many-body system. By varying temperature and dimensionality, we can tune between settings with and without symmetry breaking and accompanying collective modes or quasiparticles. We analyse both conventional and out-of-time-ordered spin correlators (`decorrelators') to track the spreading of a spatiotemporally localised perturbation -- the wingbeat of the butterfly -- as well as transport coefficients and Lyapunov exponents. We identify a number of qualitatively different regimes. Trivially, at $T=0$, there is no dynamics at all. In the limit of low temperature, $T=0^+$, integrability emerges, with infinitely long-lived magnons; here the wavepacket created by the perturbation propagates ballistically, yielding a lightcone at the spin wave velocity which thus subsumes the butterfly velocity; inside the lightcone, a pattern characteristic of the free spin wave spectrum is visible at short times. On top of this, residual interactionslead to spin wave lifetimes which, while divergent in this limit, remain finite at any nonzero $T$. At the longest times, this leads to a `standard' chaotic regime; for this regime, we show that the Lyapunov exponent is simply proportional to the inverse spin-wave lifetime. Visibly strikingly, between this and the `short-time' integrable regimes, a scarred regime emerges: here, the decorrelator is spatiotemporally highly non-uniform, being dominated by rare and random scattering events seeding secondary lightcones. As the spin correlation length decreases with increasing $T$, the distinction between these regimes disappears and at high temperature the previously studied chaotic paramagnetic regime emerges. For this, we elucidate how, somewhat counterintuitively, the ballistic butterfly velocity arises from a diffusive spin dynamics.

In parallel the classical versions of OTOCs (alternatively dubbed as decorrelators) have been applied to study the spatiotemporal chaos in classical many-body systems [37][38][39][40][41], in part with the goal to elucidate their significance in the more conventional realm of classically chaotic models, but also to understand the semi-classical limit of many-body quantum chaos. In case of chaotic classical many-body systems with short range correlations, the study of the decorrelators clearly reveals two complementary aspects of the butterfly effect-(1) the exponential temporal growth of a localised (in real space) infinitesimal difference in the initial conditions characterised by the Lyapunov exponent, λ, and (2) its ballistic spread characterised by the butterfly speed, v B [37,38].
A natural question then pertains to the dependence and relation between the above chaos time-scales, λ −1 , and lengthscales, v B λ −1 , on the thermodynamic and dynamic properties of the system. This assumes particular importance with regard to two central results about chaos in quantum-many body systems where it has been shown that: (1) in maximally chaotic system, the Lyapunov exponent is universally bounded by the absolute temperature as λ max = 2πk B T / , [8] and (2) under selected circumstances, the diffusion constant, D, of a conserved charge and the chaos time and length scales are related via D ∼ v 2 B /λ. The above two questions remain equally important in the context of classical many-body systems, where it has recently been shown that the maximal bound is generically 'violated' in the disordered phase of spin-rotation symmetric classical spin systems and thermalised fluids where λ ∝ √ T is observed [37][38][39]41]: the limits of temperature T → 0 and spin S → ∞ do not commute. In particular, in a classical spin liquid that remains disordered down to T = 0, the above behaviour was found down to the lowest temperatures as well as arXiv:2011.04700v1 [cond-mat.stat-mech] 9 Nov 2020 D S ∼ v 2 B /λ for the spin diffusion [38]. The interplay of symmetries and finite temperature however is much richer in many-body systems including the possibility of spontaneous symmetry breaking at low temperatures through a thermal phase transition that is accessed by lowering the temperature.
In this paper, we concern ourselves with the study of the spatiotemporal chaos in a paradigmatic many-body system of classical spins that is tuned through a spontaneous symmetry breaking magnetic phase transition at finite temperatures, or that at least exhibits a divergent correlation length at low temperatures. In particular, we consider short range interacting spin rotation invariant Heisenberg Hamiltonians for classical spins sitting on a d-dimensional hypercubic lattice (with d = 1, 2, 3): where S i are unit vectors encoding classical three component spins at each lattice site, i, J is the nearest-neighbour exchange constants for spins joined by bond ij , and we study both ferromagnetic, J < 0, and antiferromagnetic, J > 0, cases. The dynamics we consider is precessional, i.e. given by the Landau-Lifshits equations of motion, Eq. 12 [42,43], which are generated either by canonical Poisson brackets in the classical case, or commutators for quantum systems. Indeed, an advantage of this model system lies in such a direct connection of the classical model to the quantum system in the form of a semi-classical limit, where observables tracking the chaotic properties can easily be constructed to apply to both cases [37] providing a potentially more direct connection of classical and quantum chaos than in other situations.
We find that in practice, even though the Heisenberg magnet in d = 1 has no true long-range order at any nonzero temperature [44], that it is most convenient to focus our numerical efforts on this case: thanks to its divergent correlation length in the limit of zero temperature, it does exhibit a regime where for practical purposes, long-range order and a finite spin-wave velocity exist [45], while large linear system sizes, up to L = 4 × 10 4 spins, are numerically tractable. However, we also present numerical data on systems in d = 2, 3 with a similar number of spins but correspondingly smaller linear sizes. Indeed, most of our results are generally consistent across dimensions and in absence (d = 1, 2) or presence (d = 3) of a true thermal phase transition. The situations where the presence or absence of the phase transition is of essential importance -as in the case of the temperature dependence of the butterfly velocity (Fig. 12) -are discussed at their respective places.
We provide a systematic understanding of the following questions. Are there distinct regimes for chaos in accordance with these phases, and if so, what are their properties? How can we characterise them through observables? What role do quasiparticles, present in systems with spontaneously broken symmetry, play? A first step towards answering these questions was the study of the temperature dependence of chaos in absence of phase transitions in a classical spin system [38]. There are of course numerous precursors to this work which have considered spatiotemporal chaos in classical many-body systems, often focussing on high/infinite temperatures, or without a notion of temperature at all [46][47][48][49][50][51][52]. More recently, many-body chaos was studied close to thermal phase transitions in a scalar field theory [40], and in a classical XXZ spin system in two dimensions [41], where qualitatively different behaviour of the chaos quantities, v B and λ, was found at or below/above the phase transition.
Using a combination of direct numerical simulations and mode-coupling calculations for the spin-waves, we reveal the features of the short time emergent integrability at low temperatures ultimately giving away to chaos at intermediate times paving way for the long-time thermalisation.
The chaotic behavour is quantitatively characterised via the decorrelator, the specific measure at the centre of this investigation, which we define in Sec. III A, Eq. 4. This is part of a detailed introduction to the system and the observables to be studied, to which Sec. III is devoted.
Our narrative then proceeds from low to high temperatures, as outlined in the abstract. The low-T ordered regime is the subject of Sec. IV. There, we study in detail the spatiotemporal behaviour of the decorrelator which we use to characterise the many-body chaos. We find that for short times, it is well described by non-interacting spin wave theory. In particular, we find ballistic propagation of an initially localised perturbation, and an initial powerlaw decay of the decorrelator consistent with this ballistic spreading of spin-waves. Remarkably, the butterfly velocity, i.e. the speed with which the light-cone advances, continues to be given by the quasiparticle velocity from linear spin wave theory, even when the decorrelator overall is dominated by the expontential growth characteristic of spatiotemporally chaos.
Furthermore, at low temperatures we find a regime of a "scarred" decorrelator. This is sandwiched between the short-time integrable and the late time chaotic regimes. It arises through repeated scatterings of the long-lived, welldefined propagating quasiparticles, which seed 'secondary' lightcones. The superposition of many such secondary lightcones then generates the 'diffusive' core of late-time chaos. This diffusive core grows parametrically more slowly than the the (primary) ballistic lightcone, inside which it is located.
In Sec. V, we then present the behaviour of the butterfly velocity, characterising the spatial propagation of chaos, and the Lyapunov exponent, characterising the exponential temporal growth, as a function of temperature in d = 1, 2, 3. We find that both characteristically change behaviour at the (finite-size) thermal phase 'transition'. Whereas the Lyapunov vanishes as a powerlaw in temperature below the transition, the velocity saturates to a finite value both in the paramagnetic high-temperature and the magnetically ordered low temperature regime.
In Sec. VI we then contrast the behaviour of the ordered magnet to the previously studied case of the spin-liquid, and discuss the high-temperature paramagnetic phase, where diffusive and ballistic behaviour co-exists. Most notably, we find a direct relation, in the form of a proportionality, between the Lyapunov exponent and a combination of a characteristic velocity of the quasiparticles, and their associated scat- tering rates. From kinetic theory, this leads to the relation D ∼ v 2 qp τ qp ∼ v 2 B /λ: the quasiparticle velocity plays the role of the butterfly velocity. Finally, for the high-temperature short-range correlated paramagnetic regime, we show how the diffusive nature of the spin transport is consistent with the ballistic spreading of the decorrelator front. This section also makes contact with our previous study of chaos in a cooperative paramagnet on the kagome lattice [38], which exhibits scaling forms for chaotic observables down to T = 0, and which is relevant to the high-temperature paramagnetic phase of the hypercubic Heisenberg models studied here. We conclude with a discussion Sec. VII.
For the benefit of the reader interested in a qualitative impression of our central results without having to wade through the copious details, Sec. II provides a summary of these in the form of what we hope is a visually compelling survey of the behaviour of decorrelators in the various regimes.

II. OVERVIEW OF THE REGIMES
Here, we present the different regimes of many-body chaos that we have identified. To unclutter the presentation, we concentrate on the case of the d = 1 Heisenberg ferromagnet, and present corresponding figures for antiferromagnets, as well as higher d, in later sections. All figures display the spatiotemporal behaviour of the decorrelator, Eq. 4, on different timescales, and at different temperatures.
Integrable regime: The left panels of Fig. 1 show the shorttime spread of the perturbation evaluated in two different ways. The top one displays the results of linear spin-wave theory, a completely non-interacting theory without scattering or chaos. The lower panel displays the result of our numerical simulations. The two plots agree in considerable detail, exhibiting the following features. First, there is a clear lightcone, advancing with velocity v B ∼ |J|, beyond which the signal vanishes rapidly. The amplitude decays as a powerlaw t −1 within this light-cone. Second, there is considerable structure in the decorrelator, which reflects the properties of the spin waves across the entire spectrum, as the localised initial perturbation excites all of them.
Scarred regime: At longer times, rare scattering event become visible, middle panels of Fig. 1. In these, secondary lightcones appear at random points in spacetime, which have ballistically propagating edges in both directions. As an increasing number of these overlap as time progresses, a distinctive scarred appearance emerges.
Diffusive regime: Eventually, the number of secondary lightcones becomes so large that a given point in space is under the influence of many of them, leading to a statistical average. This is determined by the 'random walk' of the lightcones, and therefore takes on a standard diffusive shape, with equipotentials of the decorrelator following trajectories x 2 ∼ t. Note that the crossover between the scarred and the diffusive regime can take place at rather late times: for T = 0.04 at t = 1000, scars are still plainly visible even after averaging over 10 3 initial states.
This diffusive regime grows at the expense of the other two as T is raised. This connects to the analysis of the spin liquid regime presented in Ref [38], where this is the only regime present across all temperatures.

III. MODEL AND OBSERVABLES
In this section, we first introduce the central object of study for much of this work, namely the decorrelator. We then describe the thermal properties of our model system, followed by setting the stage for the analysis of its dynamics.

A. The decorrelator
We start by introducing the decorreletor: It measures the divergence of two copies in response to a perturbation applied to one of them, which can be chosen to be spatiotemporally localised, as in the wingbeat of the butterfly. This is an instance of an out of time ordered correlator, and was first introduced in Ref. [37] to study the classical many-body chaos in a spin chain at infinite temperature.
Specifically, we label the two copies of the same spin system I and II, which both evolve under the same Hamiltonian through the equation of motion given by Eq. 12. We consider an infinitesimal difference in the initial condition (i.e., at time t = 0) of the two copies only at site i = 0: where 0 < ε 1, and The decorrelator D(i, t) [18,[37][38][39]53], defined as the local difference of the spin configuration of the two copies after a time t (where · · · T denotes averaging over a set of initial conditions chosen from a canonical ensemble of spin configurations at temperature T ), measures the temporal growth and spatial spread of the difference in initial condition. In this sense δS i,t is the difference field that evolves spatially outward from i = 0 with time. It will be useful to introduce the Fourier modes for the various low energy fields. In particular, where δS k is the Fourier mode with with momentum k such that the initial condition (Eq. 2) becomes The decorrelator in Eq. 4 is now given by which contains the entire spatiotemporal information. We shall also find it useful to introduce the space averaged decorrelator defined by which, due to the averaging over the spatial information, such as the beating patterns discussed above, conveniently allows to isolate the temporal evolution. This space-averaged decorrelator is similar to the distance function studied in Ref. [54] in the context of cellular automata. We note that the question of the dynamics of a single misaligned spin in a ferromagnetically ordered background was investigated several years ago in Refs. [55] and [56]. While these studies did find ballistic spread in the magnetically ordered state, the high-temperature phase was assumed to be purely diffusive. However, as shown in Ref. [37] the decorrelator displays ballstic spreading even at infinite temperature.

B. The classical Heisenberg magnet
The model we study is the classical Heisenberg model, the Hamiltonian of which is given by Eq. 1. In this work we shall exclusively consider nearest neighbour interactions on a set of bipartite hypercubic lattices with varying dimensionality: (1) one dimensional chain, (2) two dimensional square lattice, and (3) three dimensional cubic lattice. We consider both ferromagnetic and antiferromagnetic interactions.
As the system is spin rotation invariant the three components of the total magnetisation are conserved along with the total energy E = J ij S i · S j . The fate of these conservation laws influence the thermodynamic and transport properties of the system across the entire temperature regime. We also note that due to the conservation of both energy and magnetisation during the dynamics, generically any observable computed from dynamical trajectories, such as the decorrelator described above, is a function of the total magnetisation S T and energy E (and any other conserved quantity). A given ensemble, e.g. a thermal ensemble of initial states at temperature T , then corresponds to averages over the corresponding (thermal) distributions of the conserved quantities. However a micro-canonical ensemble at fixed S T and E contains equally valid information about the dynamical properties of the model. In particular, this implies that a finite-size system with a finite magnetisation (even if the thermodynamic limit would show a vanishing magnetisation) can shed light on the dynamics in a symmetry-broken state.

The thermodynamics: phases and phase transitions
Our central tuning parameter is temperature. This allows us to access regimes with short-range correlations, and (at least effectively) long-range order.
At high temperatures, a spin rotation invariant thermal paramagnet is always realised, where the spin correlations are short ranged. At low temperature in d = 3, there is a magnetic ordering transition at a critical temperature, T c [57]. By contrast, for d ≤ 2 a thermal phase transition at a nonzero temperature is absent [44]. Instead, the correlation length diverges in the limit T → 0 in d = 1, 2. As a result, at sufficiently low temperatures, the correlation length will always exceed the size of any finite lattice. In such a regime, the behaviour of the system in many respects resembles that of a symmetry-broken phase-i.e., a uniform magnetised phase for the ferromagnetic interaction and a Neel order phase for the antiferromagnetic interactions.
We illustrate this behaviour in Fig. 2 for both the ferromagnet and the antiferromagnet model in d = 1, 2, 3 by considering the behaviour of the spin correlation length ξ versus temperature T . Note that, whereas both two and three dimensional systems show a rather well-defined feature, in one dimension there is a much smoother onset of the correlations. We will return to this point in Section V.
In order to use the decorrelator in a meaningful way at a given, common, temperature for both copies, I and II, of the system, we need to make sure that the energy difference of the two copies due to the local perturbation is consistent with arising from a typical thermal fluctuation. To ensure this, the energy difference should be small compared to the temperature: This condition is fulfilled with the following order of limits This is done in our calculations for both the linearised and non-linear equations [37].

Precessional dynamics and effective hydrodynamics
We are interested in the real-time many-body dynamics of a spin system initialised in a state representative of a particular temperature. The dynamics of the classical spin system is generated by the spin-Poisson bracket where f and g are functions of the spins. This leads to the equation of motion which is just the precession of the spins in the local exchange field due to its neighbours. The dynamics of the classical Heisenberg and related models have been extensively studied in the literature both numerically [58][59][60][61][62][63] as well as using hydrodynamic approaches [64,65]. At high temperatures spin-spin dynamical correlators are generally found to be diffusive [59,60,64] whereas below the transition the dynamical spin structure factor reveals characteristic spin wave features both in 3D with a true thermodynamic phase transition [61], and as a crossover in 2D [62,63]. For 3D, in particular, the hydrodynamic theory of spin waves [64] predicts that for the Heisenberg ferromagnet at long wavelengths, the spin-wave dispersion is given by where ρ s is the spin stiffness constant, m ord is the equilibrium magnetisation and w is the strength of spin-wave scattering. For the antiferromagnet, the low energy spin-waves in hydrodynamics obey where χ s and D denotes the spin susceptibility and diffusion constant respectively. The ballistic propagation of the linear spin-waves is correctly captured within the spin-wave theory as shown in the Appendix A. At this point we note that at least in one dimension the long-time behaviour of the dynamic correlations may be dictated by non-linear effects and need to be studied within the theory of non-linear fluctuating hydrodynamics [66]. Above T c , the conservation of the total magnetisation S T (Eq. 9) in the symmetric phase and energy E implies that for this short range interacting system both these currents are conserved leading to corresponding diffusion equations, with diffusion constants D s and D E respectively. The diffusion of these conserved quantities then captures the long-time hydrodynamics of the high temperature paramagnetic phase [67][68][69].
In the following, we will in particular discuss how this 'standard time-ordered' diffusive dynamics is related to the dynamics of an OTOC.

Dynamics of difference field
The decorrelator directly depends on the dynamical evolution of the difference field δS i = S II i − S I i defined above. Using the equation of motion (Eq. 12) and writing S II i = S I i +δS i (Eq. 2), we obtain the equation of motion for the difference, δS i,t , as where we have dropped the superscript for clarity, S I i (≡ S i ). This shows that the δS i evolves in the background of the dynamic spin field S i .
In the numerical simulations we will mainly consider the limit ε → 0, such that the second line in the equations of motion for the difference field drops out, and all quantities in the simulations are defined without factors of ε. This we refer to as the 'linearised' decorrelator, which is not to be confused with the linear spin-wave analysis: the former preserves the chaotic nature of the dynamics, while the latter does not.

IV. CHAOS AT LOW TEMPERATURE
We start our analysis of spatiotemporal chaos at low temperatures where spin rotation symmetry is spontaneously broken in an ordered state. This goes along with the emergence of long-lived Goldstone mode. The items we discuss are existence and nature of a'short-time' integrable regime of noninteracting wavepacket propagation and a long-time 'hydrodynamic' chaotic regime. A central finding is the fact that the Lyapunov exponent of this chaos is directly related to the spin-wave scattering rate. In between those two regimes, the scarred regime mentioned above appears.
For a ferromagnet (antiferromagnet on a bipartite lattice), the ordering in question is uniform (Néel) order. The follow-ing discussion is centred on the ferromagnetic case. The analysis for the antiferromagnet, which is largely analogous albeit somewhat more complicated in terms of calculations, is relegated to the Appendix C.
As remarked above, despite the lack of true long-range order in thermodynamic limit for d = 1, 2, for finite systems the rapidly increasing correlation length at low temperature leads to the effective appearance of local order capable of supporting long-lived elementary excitations. Our discussion therefore proceeds in terms of spin waves about an ordered background, appropriate at least for time scales short compared to the spin-wave life-and scattering time and length-scales smaller than the correlation length.

A. Equation of motion at low T
Our approach is to cast the equations of motion, Eq. 15 into the form of a linear term corresponding to the propagation of an 'integrable' wavepacket, subject to non-linear 'scattering' processes. The analysis of the latter will underpin our central result linking the spin-wave lifetime and the Lyapunov exponent, Eq. 30.
Deep inside the symmetry broken phase, we can expand in small deviations from the collinear ordering pattern, where n i represents the direction of collinear order whereas L i is the spin-wave amplitude (with n i · L i = 0) and the latter are the gapless Nambu-Goldstone modes describing the lowest energy long wavelength excitations about the ground state. For the ferromagnet (antiferromagnet) they have quadratic, ω ∼ k 2 (linear, ω ∼ k) dispersion (Appendix A). This suggests that the low energy dynamics of the difference field, δS i , is best understood in terms of the interaction of spin waves. To this end we use Eq. 16 in the evolution equation (Eq. 15), to obtain For the nearest neighbour ferromagnet, J < 0, we can choose without loss of generality n i =ẑ (and hence L ⊥ẑ). Therefore from Eq. 17, we get, after Fourier transforming (see Appendix B for details), For the first linear term in d = 1, 2, 3 dimensions and the matrix Z is given by Eq. A5. This term corresponds to the linearised equations of motion and accounts for the free ballistic propagation of δS k . The non-linear second term represents the scattering of δS k with the dynamic spin-waves with the scattering determined by the matrix While the detailed forms are given by Eq. B2 and A6, we note that these scatterings imply that the different Fourier modes of difference field scatter from the dynamic spin-wave modes and this results in the mode-coupling route to chaos at low temperatures as we shall see below. From the explicit forms of the scattering kernels, it is clear that the scattering vanishes as k → 0 and hence the long-wavelength modes are more long-lived, as expected for Goldstone modes.

B. The linearised de-correlator and emergent integrability
The full solution can be expressed in integral form as This can then clearly be used as the starting point of the modecoupling [70] expansions for the difference field at low temperatures.
The free solution is obtained by setting A k,q = 0. This is given by the first term of Eq. 21 as where the free propagator G 0 k (t) is given by Eq. A8. In this case each momentum mode is independent. However they interact with different spin-wave modes eventually leading to the coupling of different modes of decorrelations. We now first develop the form of the decorrelator from the free solution followed by inclusion of scattering further down.
The explicit form of the free solution is obtained from Eq. 22 and is thus where η = ( x ) 2 + ( y ) 2 = ε 1 − |L 0 (0)| 2 and hence proportional to magnetisation and φ = tan −1 (L y (0)/L x 0 (0)). Thus in the thermodynamic limit, we have (see Appendix B for details) where m T is the magnetisation and F(i, t) is a dimension dependent function given by with i ≡ (i 1 , i 2 , · · · i d ) denotes the position on a d = 1, 2, 3 dimensional cubic lattice and J ν (t) denotes Bessel function of first kind.
24 for a one-dimensional cut along the cartesian coordinate axes in dimensions d = 1, 2, 3 in the left panels of Fig. 3, and compare to the corresponding results of the numerical simulations (see E) of the full spin dynamics in the right panels. The striking similarity of the two provides crucial evidence of the emergent integrability at low temperatures in the symmetry broken phase. We note that the time and length scales shown here are chosen to be below the inverse Lyapunov time, λ −1 and below the correlation length, ξ −1 , and, thus the approximation of non-interacting spin waves in an ordered background is expected to be good as potential scattering between spin waves will be negligible. At this linear order the solution in Eq. 23 has the same structure as the linear spin-wave solution (Appendix A), but with a localised (in real space) initial condition. Therefore the short-time behaviour of the de-correlator is nothing but akin to the spread of an initially localised spin-wave packet. This, in the linear-spinwave regime is a ballistic phenomena without any exponential amplification in conformity with the full numerical calculation. We further note the dimensionality dependence of the striking beating structures visible in Fig. 3, and also reflected in the full numerical dynamics, expected for interference of non-interacting spin waves.
We provide a more detailed comparison of the above shorttime physics in the one-dimensional case via constant time slices plotted in Fig. 4. Besides the fine-structure visible within the emergent light-cone of the decorrelator, we also observe a propagating peak -the primary packet of decorrelationtravelling out at twice unit speed (in units of |J| which has been chosen to be unity in Fig. 3 and the rest of this work). Notice that this is also the case at T = ∞, however the mechanism to generate the chaos and the butterfly velocity is different as there are no spin waves at infinite temperature as shown by the analysis of the two-point spin correlators which are purely diffusive.
Finally by using the asymptotic form of the Bessel function given by we obtain powerlaw decay ∼ t −d of Eq. 24 in agreement with the full numerical simulations as an intermediate asymtotic behaviour until scattering become important. To understand the ballistic nature and the emergence of a light-cone, one needs to consider the asymptotic scaling along rays of fixed velocity v = x/t. Using stationary phase arguments, see App. D, this can be seen to result in power-law decay for v i ≤ 2J, and exponential decay for v i > 2J. This sharp dis- This free solution is expected to properly capture the short time behaviour (t λ −1 ), where no significant exponential growth of the decorrelator has taken place and hence λ ≈ 0. In Sec. V and Fig. 11, we demonstrate that the time-window for the validity of the linear solution expands and tends to diverge as the temperature is decreased to zero. We provide a more detailed picture of this short time integrability at short times and low temperatures via constant space slices in app. B, Fig. 14.
C. Long-time chaos through scattering with spin-waves We now turn to the effects of the spin waves, L, scattering with the difference field, δS i . At times much longer than the Lyapunov time, λ −1 , chaos sets in which is expected to arise due to the scattering of spin waves. Indeed, we establish that it is the scattering time of the spin waves which directly determines the Lyapunov exponent.
At the outset we note that the light-cone velocity remains unchanged even at late times when scattering cannot be neglected. This can be understood as follows. As the packet of δS i propagates outward it will eventually scatter off spinwaves in the background field, either being reflected back inside the light-cone or splitting into two spin-waves propagating forwards along the light-cone and backwards inside the light-cone respectively. Thus, the leading edge, determining the light-cone velocity, will always be dominated by the remaining weight of the initial peak that was never reflected propagating outwards at the initial velocity given by Eq. 24. In other words infinitesimally near the light-cone edge, the linear form of the de-correlator in Eq. 24 always remain valid and hence correctly predicts the light-cone velocity.
To demonstrate this, we show the late time light-cone for the one-dimensional ferromagnet in Fig. 5 for a system with L = 40000 at T = 0.04 for times up to t = 20000. We note that here the inverse Lyapunov time is t λ = 1/λ ≈ 170, thus, we are probing the dynamics deep into the chaotic temporal regime (see below), and for spatial scales considerably larger than the spin-spin correlation length ξ ≈ 20. Remarkably, we still observe a perfectly linear light-cone with a velocity determined by the spin wave velocity as for the short time dynamics discussed above.
However well within the light-cone repeated scattering results in the rapid amplification of the signal at times greater than λ −1 due to the non-linear contributions denoted by the mode-coupling term in Eq. 18. This temporal growth is most directly seen in the space averaged de-correlator, i.e., I(t) as given by Eq. 8. In fact, it is possible to show from 24 that within linear theory I(t) = ε 2 , i.e., a constant. Thus, any deviation, in particular, an exponential growth of this, is a direct indicator of the chaotic regime where scattering is important. This is shown in the main panel of Fig. 6 which for long times clearly exhibits exponential growth, I(t) ∼ e 2λt with the Lyapunov exponent, λ. Indeed, we generally find that the Lyapunov exponent based on the summed decorrelator converges to the exponential form on times of order 1/λ.
In contrast, as mentioned above, the upper left inset shows the decorrelator at the initial site of perturbation, i.e., D(i = 0, t) where the initial t −1 powelaw decay is clearly visible till t λ ∼ 1/λ. After t * the exponential growth takes over with the same Lyapunov exponent, λ, with which I(t) grows.
Finally, the lower right inset shows the short time behaviour of the summed decorrelator on a linear scale, showing only slow initial growth, which a-posteriori justifies our treatment of spin-waves as (almost) non-interacting in this regime.

Mode-coupling theory and the emergence of chaos
The mode-coupling theory for δS k introduced above (Eq. 21) can provide crucial insight into the temporal aspects of the late-time chaos in terms of the properties of the spinwaves. This section is devoted to the derivation of the result linking the Lyapunov exponent, to the spin-wave scattering rates via Eqs. 30 and 31.
Our analysis starts by expanding the integral Eq. 21 iteratively to obtain (using Eq. 22 and Eq. 6) a mode coupling expansion of δS k (t). The details are given in Appendix B. From this, we readily obtain the leading order contribution to the space-averaged de-correlator (Eq. 8) given by where · · · refer to higher order terms. The first term is indeed the constant contribution of the free decorrelator as discussed above.
The first order correction to this is due the scattering as denoted by the second term on the right hand side of the above expression. This scattering term contains two separate contributions corresponding to the two terms in Eq. 20. For the first contribution which is proportional to explicit calculation shows that this term is proportional to α β for α = β. This stems from the general antisymmetric structure of L q (Eq. B3). Hence on taking the average over the thermalised initial condition this term vanishes, since α β T ∝ δ αβ . However, the second term gives a non-zero contribution of the form which then is of the same form as the second term in Eq. 24 albeit summed over the lattice points. We note that in deriving the above expression we have not assumed free-spin waves. At second order, there are three classes of terms which can be schematically written as OO, OM and MM. Again due to the averaging the cross terms vanish. While the third term is nothing but the higher order version of Eq. 29, the first term is the sub-leading contribution with respect to Eq. 29. This provides the basis for neglecting all terms proportional to O. The dynamics of the difference field (Eq. 18) now becomes exactly equivalent to that of interacting spin-waves for a ferromagnet in Eq. A3.
This concretely shows that the same coupling of the modes that gives rise to the spin-wave lifetime also leads to chaos. The respective forms are obtained by re-summing the series in Eq. 27 within the above approximation (and similar series for the spin-waves). However, note that, in the regime where spin wave scattering leads to chaotic behaviour, the individual modes of the exponentially growing difference field, δS i , cease to be well-defined in the long-time limit. It is therefore not possible to follow the exponential growth of a particular k mode in isolation from that of the others. Then, it is the combined effect of all interacting modes which is measured in the summed decorrelator, and its form is suggestively written as where τ k is the lifetime of the kth mode. At long times the actual Lyapunov exponent defined by I(t) ∼ e 2λt will therefore be set by the lifetime of the shortlived spin waves: Recall that the limit of → 0 has been taken in the numerical calculations to open a sufficiently long time-window (in fact infinitely long for the linearised calculations as described in App. E) before the decorrelator saturates due to the finite phase-space volume of the unit sphere. The above mode coupling theory therefore connects the chaos time-scale and the lifetime of the quasiparticles in a classical system. Similar results have been proposed recently in context of quantum many-body systems, particularly Fermi-liquids [71].
We emphasize that λ is then dominated by the short scattering time-scales which are necessarily away from the ordering momentum, k = 0, where τ k diverges for the Goldstone modes and hence the present finite λ is not inconsistent with the very long-lived Goldstone modes at the longest wavelengths.

D. Cross-over between integrability and chaos: the scarred regime
We next address the nature of the cross-over between shorttime integrable and the long-time chaotic regimes. Note that already in Fig. 5, some fine-structure is visible in the decorrelator at long times. To make this more visible we study the scaled decorrelator via D scal (i, t) = D(i, t)/I(t) which removes the exponential temporal growth and thus enhances the spatial structure at given time. In Fig. 7, the scaled decorrelator exhibits a scarred appearance, due to the coexistence of a multitude of secondary light-cones at different spatial and temporal scales. We interpret these as originating at points where a given spin-wave with a well defined velocity and appreciable weight scatters of other spin-waves and get reflected, or "splits", resulting in a backward-propagating, or two spin waves propagating in opposing directions. Note that the spin wave can 'survive' a scattering event, in the sense that it splits into a part which continues to propagate in its original direction, while a second part of the signal is seeded which propagates in the opposite direction.
In turn, these products of the scattering event can individually scatter further, thus seeding further light cones; eventually, these overlap and merge, giving rise to the diffusive regime at the longest times. This has a natural interpretation in terms of the scars: with the increase in the number of lightcones, a given point in space receives contributions propagating outwards or inwards, whose statistics are those of a random walk, whence the diffusive nature.
An alternate way to visualise the crossover is to track the maximum weight of the difference field in space-time plane; for the one dimensional ferromagnet, this is shown in Fig. 8. For short times, the maximum of the decorrelator tracks the lightcone. However at later times, the scattering events nucleate growth of the decorrelator and hence increase the number of subsequent scattering events, so that the maximum then moves closer to the origin, where more lightcones overlap. In Fig. 8, this is clearly visible for the higher temperature T = 0.1 where scattering is stronger. In contrast at T = 0.04 we observe free propagation for a significantly longer time, and subsequently repeated scattering events during which the propagation direction of the main peak in the signal is fully reversed. This in turn supports that the life-time of the spinwaves is indeed larger than the scattering time, such that spinwaves undergo a random walk at fixed speed interrupted by reversal of the propagation direction when a scattering event occurs. Thus, in this regime we expect the emergence of a diffusive regime.
We demonstrate this directly in the lower panel of Fig. 7, where we show contours to the scaled decorrelator and diffusive fits to the contour lines of the form x 2 ∼ t. The results show large fluctuations, strongly enhanced by the chaotic growth, even for very large times due to the increasing life time of the spin waves at decreasing temperatures.
We note in passing that this cross-over is not the only manifestation of the existence of quasiparticles. There exists in addition a feature in the decorrelator at the initial site of the perturbation, which shows powerlaw decay in time for t 1/λ, before the exponential growth takes over. This is absent in the case of the kagome spin liquid discussed further down [38].
At this point, we note that in practice we cannot probe the scarred behaviour in dimensions larger than 1 due to the very long time and system sizes required to make this scattering and the emergent diffusive core visible.

E. Bipartite Antiferromagnet
The situation for the nearest-neighbour bipartite antiferromagnet turns out to be quite close to that of the ferromagnet. Nonetheless, there are some striking differences between the two, and we summarise the general features emphasising these differences in the following paragraphs. The details of the underlying calculations are relegated to Appendix C for the case of a Neel state, where A and B represent the two sublattices.
The free decorrelator (expressions in Eqs. C18 and C19) is compared to the numerical solutions of the full spin dynamics in Fig. 9 for short times. As in the ferromagnetic case, the two agree quantatively, both in the detailed spatiotemporal interference patterns and the ballistically propagating peaks. We would like to note that for the antiferromagnet the free decorrelator at the origin-the initial localtion of the difference field-decays much more slowly compared to the ferromagnet, a behaviour which is expected from the long time scaling in the stationary phase calculation, see App. D. However, the slow decay is soon swamped by the exponential amplification of the decorrelator as a whole. Indeed, these non-linear effects yield late time chaos like in the ferromagnet, albeit with different temperature dependence as we show below.
F. The integrable regime, and cross-over to chaos, in d > 1 The visually most striking difference between FM and AFM is the shape of the free decorrelator, plotted in Fig. 10 for d = 2 for LSW and for the full numerical solutions. The FM (left panels) displays a square shape with bright peaks at the corner. This structure is readily inferred from the factorised form of Eq. 25. By contrast, the AFM (right panels) free decorrelator comes in the shape of an isotropic circle.
This difference is explained in terms of the origin of contributions to the difference field in k-space. For the ferromagnet, the explicit expression of the Bessel function as well as the stationary phase solution of the integral expression (Eq. B9) shows that the integral is dominated by the weight of lattice dependent modes away from the ordering vector. This is related to the fact that the ferromagnetic spin waves have a quadratic dispersion. However for the antiferromagnet, the stationary phase solution (Eq. C18 and C19) is dominated, due to the linearly dispersing spin waves, by the contribution from modes near the magnetic ordering vector. While we do not have straightforward closed asymptotic expression for the integral form, the isotropy nonetheless follows from the emergent rotational symmetry near these soft modes. Indeed a stationary phase approximation shows that the isotropic butterfly speed, v B = 2J √ d, which matches very well with the numerical calculation as well as the full integral solution of the free decorrelator, App. D.
This naturally raises the question about the persistence of the square wavefront at long times for the nearest neighbour ferromagnet. We note that upon increasing temperature, the square gets replaced by an isotropic circular shape as the order underpinning the existence of the quasiparticles ceases. We note in passing that for the classical spin liquid at low temperatures [38], the decorrelator is also circular.
Indeed, we cannot exclude that that nonlinear effects at long times could restore isotropy. If this were the case the chaotic long-time decorrelator at low temperatures would also become isotropic. A direct numerical verification of this is beyond our present numerical capacities even for d = 2 on account of the rapidly growing length/timescales. At any rate, note that the decorrelator at these scales will be dominated by the exponentially growing core. This however only grows diffusively, in contrast to the ballistic spreading of the free decorrelator. A parametrically separated coexistence between a square form of the latter, and a circular one of the former, is hence another possibility.

V. TEMPERATURE DEPENDENCE OF CHAOS SCALES IN d = 1, 2, 3
We now collate our numerical results for the temperature dependence of the central chaos scales, the Lyapunov exponent and the butterfly velocity, for the hypercubic lattices in d = 1, 2, 3, considering both the ferro-and the antiferromagnet.

A. Lyapunov exponent
The Lyapunov exponent, Fig. 11, was obtained from a fit to the time-dependence of the decorrelator. The data shown results from considering its behaviour at the initial site D(x = 0, t) ∼ e λt . We account for the initial powerlaw decay in the ordered regime (see Fig. 6) by only fitting after the cross-over time t λ ∼ 1/λ where the behaviour is clearly exponential. We have corroborated these values by comparison to the spatially integrated decorrelator I(t) ∼ e λt , which does not show the initial power law decay, yielding consistent exponents.
Generally, the Lyapunov exponent exhibits two distinct regimes. In the high-temperature, short-range correlated regime, there is little temperature dependence. This is not at all surprising: the state of a system with a bounded local energy spectrum changes only little when the temperature is raised well above this bandwidth. In this regime, λ ∼ |J|, only weakly dependent on the dimension, is thus The second, T -dependent regime, is entered as correlations start to develop on a scale set by J, with a crossover to a welldeveloped low-temperature power-law regime, λ(T ) ∼ T α with α > 0. One may try to understand this decrease of the Lyapunov exponent at lower temperaturs from the reduction of the available phase space volume. As the energetic constraints become active around the cross-over temperature the entropy of the system is reduced, fluctuations decrease overall, and the dynamics slows down. Indeed, even for the cooperative paramagnet in d = 2 for which the correlation length saturates to a small value in the low-T limit, there is a powerlaw λ(T ) ∼ T 1/2 [38]. Noticeably, in Fig. 11, while the high to low temperature crossover is smooth for λ in d = 1, 2, the transition for d = 3 leads to the appearance of a kink for both the FM and the AFM.
The actual value of the exponent of the low-temperature scaling, α, itself depends both on the dimension and the sign of the exchange constant J. Thus, unlike the thermodynamic static properties, including the position of the transition, which are equivalent for FM and AFM classically, their dynamic chaotic behaviour and temperature scaling is different. This is of course not unexpected, in view of their differing respective (quadratic and linear) low energy dispersions.
Given the correlations are longer-ranged, and the long- wavelength quasiparticles increasingly well-developed (and strictly long-lived for d = 3 below the transition), a slowing down of the dynamics compared to the cooperative paramagnet is again unsurprising: the power α in Fig. 11 is greater than 1/2 in all cases. However Eq. 31 implies that chaos is dominated by the short spin-wave lifetimes, i.e. not from the vicinity of the Goldstone modes (Eq. 13 and 14). This, we believe, is also the reason why an evaluation of τ k for longwavelengths fails to account for the observed value of α in Fig. 11.

B. Butterfly velocity
The temperature dependence of the butterfly velocity is shown in Fig. 12. This is obtained from fits to the position along lattice axes where the decorrelator exceeds a threshold D 0 , e.g. x thr = vt with D(x thr , t) > D 0 for D 0 = 10 −10 . The resulting velocity is independent of the threshold for sufficiently small thresholds in the range we can probe, and consistent with the propagation velocity of the main peak (see Fig. 4). We however note that for the FM in the ordered regime it will be direction dependent, i.e. the speed will be minimal along lattice-axes, and maximal along the body-diagonal due to the cubic form of the wave-fronts.
In the high temperature regime, note that a light-cone arises despite the absence of spin waves and their 'ballistic' propagation. In fact the dynamic spin-spin correlator in this regime is diffusive [37,38]. As for the Lyapunov exponent, in this regime the butterfly speed is essentially constant and determined by the strength of the exchange coupling.
In contrast, in the low-temperature regime, the spin-waves as shown above are well defined quasi-particles even deep into the chaotic regime. In particular the velocity calculated in the previous section from the free de-correlator agrees well with the butterfly velocity for the lowest measured temperatures, particularly for the antiferromagnet where the wavefront is isotropic and along the body-diagonal for the FM where we d, e.g. 2, 2.82 and 3.46 in d = 1, 2, 3. For d = 3 there is a gradual hardening of the butterfly velocity across the transition temperature. This is more prominent for the antiferromagnet, where it seems to follow the pronounced increase of the spin stiffness, which determines the spin-wave velocity via Eq. 14. Ref. [41] also studied the Lyapunov exponent and butterfly velocity across a BKT and Ising transition in two-dimensional XXZ spin models, with similar results as ours detailed above. Intriguingly, Ref. [41] did not detect any sharp signature of phase transitions in the these measures of chaos.

VI. COMPARISON BETWEEN ORDERED MAGNET AND SPIN LIQUID
The role of ordering with the concomitant appearance of quasiparticles can be crisply juxtaposed to the situation in the case of the classical kagome spin-liquid [38], where no ordering occurs down to T = 0 on account of the geometric frustration of the spin interactions. Both the Lyapunov exponent and the light-cone velocity show smooth crossovers from the hightemperature paramagentic regime (common to both cases) to the low-temperature cooperative paramagnetic, or spin-liquid, regime. In this regime, both quantities vanish algebraically with temperature, obeying the following relation with the spin diffusion constant [38]: The coexistence of the diffusion constant, D, with the ballistic butterfly speed v B may at first sight seem somewhat surprising. Below, we first provide a simple picture for how these two different types of behaviour can coexist naturally. This yields a remarkable connection between the physics of chaos, and that of hydrodynamics, linking spatio-temporal chaos time-and length-scales to the transport coefficient of a conserved charge.
We then contrast this expression with the situation in a system with quasiparticles, where the diffusion constant within an kinetic theory set-up appropriate for this situation, is given in terms of a characteristic velocity associated with the quasiparticles, v qp , and their scattering rate, λ qp , Collecting our previous results implies that the butterfly velocity is simply replaced by the characteristic speed of ballistic propagation of the quasiparticles, with λ qp related to the Lyapunov exponent via Eq. 31. Then, from our above discussion we conclude that while for antiferromagnet v qp indeed represents the velocity of the long wavelength spin waves as calculated from the hydrodynamic theory, for the ferromagnet modes away from the ordering modes set this scale.

A. Coexistence of diffusive and ballistic correlators
For an insight into the form of Eq. 33, it is convenient to define the average, and difference, combinations of the two copies, S I T , S II T introduced to define the decorrelator in Eq. 4. Both the total difference (Eq. 2) as well as the total average magnetisation are conserved. Therefore the corresponding densities of δS i and S + i are expected to diffuse (it is straightforward to see from Eq. 12 that the dynamics of the two densities, though coupled with each other, are also local) with diffusion constants (say) D − and D + respectively. These diffusions can be studied through the respective two two-point correlators and and we expect Indeed for the kagome spin liquid, these two diffusion constants are plotted as a function of temperature in Fig. 13 (left  panel).
For our protocol, the diffusion of δS it suggest that the difference field evolves as where is n i,t is a vector that encodes magnitude and direction of the fast local fluctuations in the difference field. This describes the motion of the difference field throughout the system starting from i = 0 at t = 0 such that for all times. For a chaotic system, the growth of this difference field at a given location, i, is expected to be exponential, exp(λt), for x t e iωt e −iqx C±(x, t), via S±(q, ω) ∼ 1 ω 2 +κ(q) 2 , with κ(q) = Dq 2 depending quadratically on q around q = 0. Right Panel: Ratio v 2 b /λ of the square of the butterfly velocity to the Lyapunov exponent versus temperature T . The Lyapunov is defined from arbitrarily long times in the limit → 0. While the random directions of n i,t ensure that Eq. 41 is obeyed, i.e. the decorrelator, is expected to have a form where λ is the Lyapunov exponent and by definition the butterfly speed, v B = D − λ.
The constants D ± refer to the diffusion of two different conserved quantities, and they need not be equal in general. This is most prominent in case that the diffusion is due to separate set of objects (particles/excitations) that carry one of the two charges of symmetric and staggered spin-rotation symmetry, O ± (3). However, in case the same objects carry both the charges, then a Wiedemann-Franz type law can emerge, but now relating the two diffusion constants, would be expected, where η is a constant which may depend on the energy scales, as well as details, of the system under consideration. In this case, the spin diffusion constant is as advertised above, Eq. 33. The numerical verification of these ideas is shown in the right panel of Fig. 13 for the classical kagome spin liquid which exhibits no ordering at any T , thereby allowing fits over a large range of T .

B. Diffusion and chaos with and without quasiparticles
We now return to the observation that a similar relationship between the chaos and diffusion holds both for the symmetry broken case, Eq. 34, as well as the spin liquid, Eq. 33.
In case of the ordered phase, we obtained a relation between the Lyapunov exponent and a time scale λ ∼ τ −1 , Eq. 31, where τ is an appropriately defined single-particle lifetime of the low energy quasi-particles, the spin waves. These spin waves propagate with a velocity v qp ∼ J, which subsumes the ballistic propagation of the perturbation wave front, v B ∼ v qp . In this situation, the kinetic theory of dilute gases yields a diffusion constant, D + ∼ v 2 qp τ ∼ v 2 B /λ. Note that, despite the visual similarity between Eq. 34 and Eq. 33, their underlying physics differs fundamentally: there is no straightforward quasi-particle description for the spin liquid, whose low-energy sector with its huge ground state degeneracy is completely unlike that of the ordered magnet with its emergent integrability and long-lived quasiparticles. In this sense, the emergence of the butterfly velocity described in the previous subsection is an entirely separate, and remarkable, feature of many-body chaos at low temperature.
In passing, we note that for the symmetric combination, i.e., the average magnetisation density, the fast local fluctuations, however cannot grow exponentially due to the already large background present in the form of a local spin-length. Therefore while the constraint (that can be derived from the equation of motion for the spins) indicates that the symmetric part contains the same information about the chaotic properties, e.g. the butterfly velocity and the Lyapunov, as the decorrelator, we note that the signal might in practice be hidden under the magnetic fluctuations and impossible to extract from the symmetric part. However, as noted in Ref. 53, in classical spin systems without spin conservation (say a XYZ model) the initial difference of total magnetisation between the two copies grows exponentially.

VII. SUMMARY AND OUTLOOK
In summary, we have presented a study of chaos and its temperature dependence in a family of model Hamiltonian many-body system which allows considerable variation in terms of choice of lattice, dimension and interaction. This provides access to different thermodynamic phases, such as a disordered [paramagnetic], ordered [(anti)ferromagnetic], as well as the critical regime separating these. This has in turn permitted us to identify a number of distinct regimes characterised by different natural degrees of freedom, transport mechanisms, and concomitant velocity and time scales.
Our combined numerical and analytical investigations concretely connect the signatures of many-body chaos in classical spin systems, the Lyapunov exponent λ and the butterfly velocity v B , with the velocity v qp and scattering time τ qp associated with the spin waves, i.e. the quasiparticles in the phase which spontaneously breaks spin rotation symmetry. These relations are directly manifested in the de-correlator, which quantifies how two copies with weakly, and locally (in real space), perturbed initial conditions diverge in time and space.
At low T and short times (t < λ −1 ) integrable behaviour emerges whose butterfly velocity and power-law temporal decay are concretely captured within linear spin-wave theory. Interestingly, the shape of the wave-front may depend on the interactions, being 'hypercubic' in shape for the nearest neighbour ferromagnet in d > 1, while the antiferromagnet exhibits isotropic 'hyperspherical' behaviour. This is quite a tangential point to the present study, but besides some early work [56], we are not aware of a systematic study of the properties of freely propagating disturbances across lattices as a function of Hamiltonian parameters. Indeed it would be an interesting question how the symmetries of the considered lattices and further-neighbour couplings affect the observed wave-fronts. We also leave open the fate of the shape of these wave-fronts at longer times when interactions become relevant, which might be sufficient to restore isotropy even on these lattices.
The integrability, however, is only approximate and thus a transient phenomenon, even though the corresponding timescale, the spin-wave lifetime, τ k (where k is the spinwave momentum), grows as the temperature is decreased. The short time integrable behaviour thus gives way to fully developed chaos as witnessed by the exponential temporal growth of the decorrelator in addition to its ballistic spread through an intermediate scarred region. We interpret this a consequence of well defined quasi-particles, namely the spin-waves, transporting weight of the decorrelator ballistically while undergoing repeated scattering events. This results in random-walk like behaviour with diffusive core of the decorrelator on top of the exponential chaotic temporal growth. The novel scarred regime still holds plenty of interest for further study, for instance regarding the detailed mechanism leading to the generation of the secondary lightcones, as well as their statistical distribution in what appears to be a regime dominated by rare events.
In the ordered low temperature regime, the Lyapunov exponent also vanishes as a power law in temperature, albeit with a different exponent which in turn depends on both dimension and sign of the coupling; at the same time, the light-cone velocity gets subsumed by the velocity of the ballistically propagating quasi-particles, and also saturates to a finite value at low temperatures. We also find that both the Lyapunov and the light-cone velocity show characteristic features at the phasetransition, where the Lyapunov changes from esentially constant in the paramagnetic regime to a powerlaw decay, and the velocity shows a minimum for the ferromagnet and a characteristic stiffening for the d = 3 antiferromagnet reflecting the emergent spin stiffness.
We note that the same spin-wave scattering is responsible for the thermalisation of the weakly interacting gas of spinwaves. This notion of a dilute thermalised gas of spin-waves then forms the basis of the kinetic theory of transport at low temperature. Indeed the above picture is generic and forms one of the central pillars of low temperature transport theory in symmetry broken systems, and our work ties this in with the nature of the concomitant many-body chaos.
Such a low energy transport theory can be contrasted to the low temperature cooperative paramagnet, where quasiparticles are absent but diffusive transport persists. There, the transport coefficients are directly determined by the chaos timescales and lengthscales, with both Lyapunov exponent and light-cone velocity exhibiting a power-law temperature dependence.
The present observation of the importance of chaos time/length scales for hydrodynamics in presence and absence of quasiparticles quantitatively indicates an intriguing and important role of many-body chaos in the transport of strongly correlated systems which we think are of much broader significance in both classical and quantum many-body systems. Indeed, the connection between classical and quantum chaos presents one of the most fascinating aspects of the field of many-body dynamics at present. While for the quantum setting a number of results -under the headings of bounds for chaos and relatedly, Planckian transport -have been established, their fate in the semiclassical limit is a subject of current investigation [72][73][74].
We would like to end with two open interesting questions pertaining to two well-known frameworks of transport in the context of symmetry breaking and thermal phase transitions. The first pertains to the connection between the present observations and the elegant theory of non-linear fluctuating hydrodynamics [66] for the dynamic correlation function which characterises the long distance and time scaling of the low energy modes in the broken symmetry phase. Application of such approaches to classical spin systems in one dimensions are rather recent [65] and their connection to chaos remains to be understood. The second is related to the similar connection between our present approach and the theory of dynamic critical phenomena [64]. While the characterisation of chaos near a critical point has been recently addressed in model systems [40,41], concrete connection to the different forms of coarse grained hydrodynamics and their possible relationship with many-body chaos remains open. In either case the systematic understanding of the mode-coupling there for the difference field, as developed here, and its similarity with the spin-wave dynamics can serve as a starting point of such an understanding.
of the spin-waves as

Nearest Neighbour Ferromagnet
For the ferromagnet, J < 0, the ground state is given by n i =ẑ (andn ⊥ L). The transverse fluctuations are described by and where γ k is given by Eq. 19. The linear solution is obtained by considering the bare Green's function which, written explicitly (in real space), is of the form Therefore M k,q represents three spin-wave modes scattering, i.e., the leading scattering term. Within a k-dependent relaxation time approximation, we can re-write Eq. A3 as where τ k is the lifetime of a spin-wave. This is similar to the Landau-Lifshits-Gilbert equation.

Nearest neighbour bipartite antiferromagnet (Neel order)
For an Néel state (Eq. 32), as in the ferromagnetic case, Eq. A1 gives and Now, defining we get, after Fourier transforming, From this, again similarly to the FM case, for the free spinwaves, we get This gives where Therefore we have The above equations can be solved to get A k and B k . From this we find that it is useful to define which diagonalizes the propagator in spin-space. Thus we get and We re-write Eq. A32 as To understand the effect of scattering, it is useful to re-write Eq. A18 as describes the free evolution and is the scattering matrix with elements given bȳ This describes the scattering of the spin waves and leads to their finite lifetime.
Appendix B: Decorrelator for the Ferromagnet From Eq. 17, for the ferromagnet with n i =ẑ, the equation of motion for the difference field reads The Fourier transformation (Eq. 5) of the above equation leads to (with L k = (L x kx + L y kŷ )) Eq. 18, with in Eq. 20, and M k,q is given by Eq. A6. The decorrelator can then be calculated from Eq. 7.

Free solution
Neglecting scattering, the explicit form of the linear equation (first term on the right and side of Eq. 18) of motion of δS i is given by whose solution is summarised in Eq. 22. The initial conditions are obtained as follows. From Eqs. 2 and 3, we have Eq. 6 such that Therefore, using Eq. 23, the decorrelator reads which can be re-written as Eq. 24 by using and associated properties of the Bessel function of the first kind, J ν (x). In Eq. B9 we have also used the relation in thermodynamic equilibrium, |L 0 (0)| 2 = 1 − m 2 T .

Mode coupling and scattering of difference field with spin waves
The mode-coupling expansion for δS k (t) is obtained by iterating Eq. 21 and is given by From Eq. B11, the expansion of the decorrelator is obtained via an expansion of δS k (t)·δS k (t) whose perturbation series is given by 27. Note that in deriving Eq. 27 we have used Appendix C: Decorrelator for the bipartite antiferromagnet Next, we discuss the spread of decorrelations for the bipartite antiferromagnet. From Eq. 17, using Eq. 32 for a bipartite antiferromagnet, we get for sublattice A and for sublattice B. Introducing the symmetric and antisymmetric modes for the difference field similar to the antiferromagnetic spin-waves in the section above, Fourier transforming yields The various scattering terms are given by with L ± q being given by equations similar to Eq. B3 with L ± q appropriately used for the two spin wave modes (see above); and M ± , N and R are given by Eqs. A20 and A22.

The free de-correlator
Neglecting the scattering terms in Eq. C4, the equation of motion for the free decorrelator is Again due to the structure of the matrix Z (see Eq. A5), the longitudinal component (along the ordering direction) does not evolve. Explicitly solving Eq. C11 for the transverse components we get δS 0 where Γ k and ρ k are defined in Eqs. A26 and A29 respectively. Now, similarly to Eq. A31, we define ∆ k as withγ k and ρ k being defined by Eqs. A24 and A29. The free solution can now be written as where G 0 k is the free propagator defined in Eq. A35 and We note that in case of the antiferromagnet, we have used the initial conditions as given by Eq. 6 and B6. This breaks sublattice symmetry by having the perturbation initially at sublattice A.
For sublattice A and B the free correlators are then given by and We note that the initial perturbation was put on sublattice A, as is evident from the above expressions and n T is the Neel order parameter.

Scattering and chaos
Incorporating the scattering, the equation of motion in terms of the ∆ k introduced in Eq. C14 can be written as where χ k , given by Eq. A37, controls the free evolution and Ξ k,q is the scattering kernel (that couples the different modes, similar to Eq. 18) given by with Ξ 1 k,q given by Eq. A38 and where now we havē M ± k,q ,N k,q andR k,q are given by Eqs. A40 and A42. The calculation for ∆ k now proceeds very similarly as in the ferromagnet, leading to the solution, cf. Eq. 21 : where the free solution ∆ 0 k (t) is given by Eq. C15. This can be expanded in terms of the free solution as A diagrammatic representation similar to Eq. B11 is obtained for both δS as well as the decorrelator (not shown).
We begin with the I 0 term. This is already not fully trivial. Firstly, the dispersion is not a smooth function. However, evaluating the contribution due to the non-differentiable part at the centre of the Brillouin zone d d k e i|k|t ∼ 1/t d (D5) we find this contribution to be subleading. In addition, the points of stationary phase are not isolated, but rather form a d − 1-dimensional manifold defined via i cos(k i ) = 0. In 1D this does not change things and we have the same scaling as for the FM, e.g. t −1 for the decorrelator. In higher dimensions this complicates things. Furthermore, while in 2D generically the critical points have a nonvanishing Hessian, for special points, e.g. (k x , k y ) = (π, 0) the Hessian vanishes.
We consider these contributions in turn. A (d-1) manifold of normal saddles gives a contribution of the form In 2D we also have higher order saddles of the form dk x dk y e i(kx−π) 2 k 2 y t ∼ log(t)/ √ t (D7) Thus, in 2D and 3D we obtain additional contributions dominant compared to the normally expected t −d/2 scaling, and there is no clear powerlaw behaviour for the AFM in 2D and 3D.
Next we consider the I 1 and I 2 terms. The additional presence of ρ(k) (1/ρ(k)) does not change the asymptotic scaling for those critical points where ρ(k) (1/ρ(k)) attains a finite value. Thus, we do not need to revisit the discussion of the d − 1 manifold where cos(k i ) = 0, and only need to consider the effect on isolated critical points.

Comparison
Thus, for the FM we generically obtain powerlaw decay in time, whereas for AFM we obtain a constant contribution to the decorrelator in 1D which will dominate at long times, and complicated crossovers between different scalings in 2D and 3D.
In addition, the stationary phase calculation suggests a cubic wavefront for the FM and an isotropic one for the AFM as borne out in the explicit calculations and full numerical calculations.
needs to be kept in mind when comparing to the expressions explicitly containing these factors. We then average 10 3 such trajectories to compute the thermal average.
We note that since the equations of motion preserve the magnetisation, the difference between the two copies is limited at low temperature by the ordered moment, as only the transverse component of the spins can decorrelate which has length of order √ 1 − m 2 . Similarly, in the linearised equations of motion the growth will be in the transverse components as well. Thus, in principle, one may distinguish the longitudinal (parallel to the ordered moment) and the transverse (perpendicular to the ordered moment) components of the decorrelator, which thus acquires a tensorial structure. For the results presented in this work, we have not addressed this additional structure, considering the sum over all components as defined in Eq. 4. This is indeed dominated by the transverse components.
This ties into another advantage of the approach using the linearised equations of motion. While the exponential growth is still limited to the transverse components, it is not limited in size or time, such that we obtain a clean exponential growth up to arbitrary long times, rather than only over a finite window as in the full non-linear equations of motion between two copies.

Role of finite magnetisation
We discuss two additional complication due to the finite magnetisation. Firstly, when comparing to a perturbed copy, one may be concerned that our perturbation, Eq. 2, does not preserve the order parameter. Of course, for small ε this change is also small. Still, we performed simulations with a modified prescription, where two neighbouring spins are rotated around the order parameter by the angle ε instead, with no difference in the results.
Secondly, the equations of motion result in a precession of all spins around the magnetisation vector. For the FM the order parameter is the total magnetisation and thus constant in time. However, for the AFM, the order parameter is the staggered moment which is not exactly preserved under the dynamical evolution. While the dynamics of the staggered moment slows down as system size increases, to vanish in the thermodynamic limit, it is present on finite systems. We have therefore computed all observables for the AFM in a static frame, and a dynamic frame defined by a rotation which is chosen such that the staggered moment remains constant. This is primarily important when distinguishing the longitudinal and transverse components which depend on the orientation of the order parameter, rather than the total magnetisation.