Many-body chaos and anomalous diffusion across thermal phase transitions in two dimensions

Chaos is an important characterization of classical dynamical systems. How is chaos linked to the long-time dynamics of collective modes across phases and phase transitions? We address this by studying chaos across Ising and Kosterlitz-Thouless transitions in classical XXZ model. We show that spatio-temporal chaotic properties have crossovers across the transitions and distinct temperature dependence in the high and low-temperature phases which show normal and anomalous diffusions, respectively. Our results also provide new insights into the dynamics of interacting quantum systems in the semiclassical limit.


Introduction
Chaotic systems are described by a growth rate, the maximum Lyapunov exponent λ L (> 0), of perturbation to the initial condition. In recent years, a quantum Lyapunov exponent, and a butterfly velocity v B for ballistic spread of local perturbation, computed from the out-oftime-order commutator (OTOC) have emerged as important measures for chaos in quantum many-body systems having some well-defined semiclassical limit [1][2][3][4]. In these quantum systems where the Lyapunov exponent can be extracted, it has been perceived as a rate for early-to-intermediate-time thermalization. Here, thermalization refers to the emergence of statistical mechanical description during dynamical evolution of the system. However, it is still controversial [5] whether chaos in isolated interacting classical systems is essential for thermalization.
The recent interest in chaos in quantum many-body systems stems from the proof of a remarkable upper bound, 2πk B T / , for λ L [6] and a relation, D ∼ v 2 B /λ L , between diffusion coefficient D, a quantity related to transport, and λ L and v B in certain strongly interacting systems [7]. Such bounds have been phenomenologically conjectured for transport scattering rate [8,9], or more conventional relaxation rates extracted from usual time-ordered correlation functions. However, these conjectured bounds have not been concretely established. As a result, even if we set aside the issue of actual role of chaos in thermalization, the temperature (T ) dependence of λ L , and other quantities related to chaos, can serve as a fundamental characterization of phases and phase transition as various chaotic fixed points, like certain non-Fermi liquids and Fermi liquids [10][11][12][13]. The former are highly chaotic with λ L ∼ T at low temperature, whereas the weakly interacting Fermi liquids show λ L ∼ T 2 , essentially dictated by quasiparticle decay rate. Motivated by these, here we ask whether many-body chaos can be used to classify phases and finite-T phase transitions in classical systems with intrinsic dynamics. An example of such systems is interacting classical Heisenberg spins on a lattice with precession dynamics.
However, typically, chaos is a probe of the short-and intermediate-time behaviour. On the other hand, the long-time dynamical properties of interacting many-body systems, in symmetry broken and unbroken phases, and across phase transitions, are mostly characterized by the properties of the collective low-energy excitations, hydrodynamic and critical modes. How are the short-time chaotic properties of many-body systems related to their long-time dynamics? To address these questions, we look into the connection of chaos with transport, characterized in terms of usual dynamical spin-spin correlations in the spin system. Classical many-body chaos under Poisson-bracket spin dynamics has been studied in some recent works for various spin models at high [14] and low temperatures [15], as well as across a spin-glass transition in a zero-dimensional spin model [16]. Also, there have been similar studies in other dynamical models like Burgers hydrodynamics [17], and in the classical limit of a relativistic field theory [18]. There are many earlier studies of maximum Lyapunov exponent as well as full Lyapunov spectrum in spin models, molecular systems, and across classical phase transitions [19][20][21][22][23][24][25][26][27][28][29][30]. The conventional Lyapunov spectrum analysis, in principle, can reveal the spatiotemporal structure of chaos [31][32][33]. However, Lyapunov spectrum analysis typically makes the direct information about the spatial structure somewhat obscure. In this regard, a spatiotemporal correlation function like OTOC [14,34] provides much more transparent way by explicitly revealing spatial spread of a local perturbation. Thus, in contrast to earlier works [19][20][21][22][23][24][25][26][27][28][29][30], using a classical version of OTOC [14], we establish the detailed temperature dependence of both temporal (λ L ) and spatial (v B ) characteristics of chaos across two classic thermal phase transitions in two dimensions. The models are described by a microscopic spin dynamics that is directly connected with the quantum dynamics in the semiclassical limit. We show that the chaotic properties, in general, are rather impervious to the nature of transport, namely whether the system exhibits diffusion or anomalous diffusion.

Model, dynamics and OTOC
We study the classical XXZ model on a square lattice described by the Hamiltonian where S r = (S x r , S y r , S z r ) are unit length spin vectors on lattice site r with total N sites and J is the coupling between spins on the nearest neighbor bonds along δ = ±x, ±ŷ directions. The anisotropy ∆ ≥ 0 can be varied to change the nature of the finite-T phase transition. For example, the system undergoes a transition at a non-zero temperature T KT in the Kosterlitz-Thouless (KT) [35] universality class for ∆ < 1 (easy plane), and an Ising transition for ∆ > 1 (easy axis), while the isotropic (∆ = 1) point does not have any finite-T transition. The chaotic properties of the isotropic point have been studied by Bilitewski et al. [36] in an independent work.
We study chaotic properties of the model in Eq.(1) using the classical OTOC [14], along with more conventional dynamical spin correlation function S r (t) · S r (0) , for the Poisson bracket dynamics The Poisson bracket of two functions f ({S r }) and g({S r }) is defined as {f, g} = rijk is the effective field on the spin at r. For ∆ = 1, apart from total energy, the dynamics conserves S z total = r S z r . This hydrodynamic mode is expected to lead to diffusive behaviour for dynamical spin correlation function at long times.
To characterize the chaotic properties of the model of Eq.(1), we use a classical version of OTOC, or the so-called cross correlator or decorrelator, introduced in ref. [14], where a and b denote two copies of the initial configuration (t = 0), with b slightly perturbed from a at r = 0 such that S br (0) = S ar (0) + δS 0 δ r,0 . The small perturbation, |δS 0 | ≈ ε, is chosen to be orthogonal to S a0 (0), i.e. δS 0 · S a0 (0) = 0. More specifically, following Ref. [14], we generate the initial configuration {S br (0)} for replica b by rotating S a0 slightly about a unit vectorn = (ẑ×S a0 )/|(ẑ×S a0 )| such that the perturbation at t = 0 becomes δS 0 = ε(n×S a0 ). The perturbation only preserves the normalization of the spin S b0 at r = 0 up to O(ε), i.e. S 2 b0 1 + O(ε 2 ). For the particular choice of perturbation, the connection of D(r, t) with the quantum out-of-time-ordered commutator [S r (t),n · S 0 (0)] 2 in the semi-classical limit, where commutator is replaced by Poisson bracket, has been discussed in Ref. [14]. We note that D(r = 0, 0) identically zero, and, since δS 0 ⊥ S a0 , D(0, 0) = 0 too as S ar (0) · S br (0) = 1 at any r. Thus, D(r, t) starts from zero for any r due to the special choice of the initial orthogonal perturbation. The averaging . . . is over initial equilibrated spin configurations {S ar (0)} drawn from a thermal distribution ∝ e −H({Sar(0)}) /T (Boltzmann constant k B = 1). Starting from the slightly different initial conditions as discussed above, the two copies are time evolved independently via spin-precession dynamics of Eq.(2). The classical OTOC D(r, t) measures the amount of de-correlation at (r, t) between the configurations in the two replicas or the trajectories, which are almost completely correlated at t = 0. The classical OTOC differs at O(ε 2 ) for r = 0, t = 0 from the more conventional measure of spatiotemporal divergence of two trajectories [37], (δS r (t)) 2 = (S ar (t) − S br (t)) 2 , which we denote as trajectory divergence for brevity. Starting at (δS 0 (0)) 2 ∼ O(ε 2 ) at t = 0, the trajectory divergence is expected to grow exponentially at r = 0 as ε 2 e 2λ L t over a Lyapunov time window t ∼ λ −1 L ln ε −2 in a chaotic system. We show that (δS r (t)) 2 and the classical OTOC of Eq.(3), both have an early-time regime 0 ≤ t ≤ t 0 , where (δS 0 (t)) 2 and D(0, t) initially change non-exponentially to O(ε 2 ), and then grows exponentially for t > t 0 . In fact, both (δS 0 (t)) 2 and D(0, t) initially decrease before start increasing with a power-law time dependence till t 0 . Nonetheless, t 0 is found to be closely connected with the chaos time scale 1/λ L , as we discuss later.
The main motivations for studying the spatio-temporal OTOC [Eq.(3)] in the model of Eq.(1) are twofold. One, as stated in the introduction, is to dynamically characterize the thermodynamic phase diagram of a spin model with well-known two-dimensional (2d) phase transitions. The XXZ model allows to tune the relative contribution of various hydrodynamic, low-energy and critical modes in the dynamics by changing temperature and anisotropy, and thus to probe the potential role of these collective modes on chaos. The second motivation comes from the fact that the spin precession dynamics [Eq.(2)] can be obtained as a classical large-S (spin) limit of the Heisenberg equation of motion for the quantum XXZ model. Hence, chaotic properties of such classical model can give useful insights even about the quantum model. The results from the classical dynamics could be particularly relevant near finite-temperature continuous phase transitions, where quantum effects for the dynamics are generally believed to be irrelevant [38] due to divergent length and time scales.
In quantum systems, truly chaotic behaviour, namely the exponential growth of OTOC [34], can only be observed certain large-N models, e.g. Sachdev-Ye-Kitaev (SYK) and related models dual to black holes [1,3,4,7,10,11,39], other large-N theories [40][41][42], and weakly interacting systems with semiclassical quasiparticle dynamics [43,44]. In these models the exponential growth can be observed over a parametrically long time window between t ∼ λ −1 L and λ −1 L ln N or λ −1 L ln(1/ ) for large N or the semiclassical ( → 0) limits, respectively. The large-N models are either infinite range or have a large local Hilbert space. In contrast, short-range quantum models with finite local Hilbert space, and without any semiclassical limit, typically do not show any exponential growth regime in OTOC [34,45,46]. This lack of exponential growth is either simply due to the absence of chaotic growth or else due to very short, and thus unresolvable, temporal window of the growth. It is an outstanding unresolved question whether such quantum systems can show chaos. However, as shown in Ref. [47], the semiclassical limit, though sufficient, may not be a necessary condition to observe the exponential growth. In particular, even for a short-range quantum model with finite local Hilbert space and without any obvious semiclassical limit, the exponential growth may be ascertained through a suitably defined spatially integrated OTOC if v B /λ L 1, where is a microscopic length scale. Based on our calculations in the classical limit, we identify a possible temperature regime in the XXZ model where such a condition could be satisfied, and thus the exponential growth may be observed even in the quantum limit. Moreover, as remarked earlier, quantum effects typically become unimportant near finite-temperature continuous phase transitions. Thus, one can naively conjecture that even short-range quantum models with finite local Hilbert space may show chaos due to effective coarse graining of degrees of freedom near the transitions. However, chaos is only short and intermediate-time property and maybe unaffected by such critical coarse-graining at long time scale. Nevertheless, the exploration of this possibility will require the simulation of real-time dynamics of the quantum XXZ model across the 2d transitions and is beyond the scope of this paper. Here we only study the chaotic properties across the phase transition in the large-S limit of the XXZ model.
In the context of the interrelation between chaos and the dynamics of collective mode, we particularly focus on the dependence of chaos on the nature of transport in the presences of conserved quantities. Interesting interplay between operator spreading characterized via OTOC and diffusion due to conserved modes have been explored in quantum systems [48,49], albeit in the toy models of random unitary circuits [50,51]. However, these toy models are non chaotic from the perspective of exponential growth of OTOC, though they can be classified as quantum chaotic based on other diagnostics, like entanglement growth [52]. For the chaotic quantum systems of strongly interacting diffusive metal [7,39] built from solvable large-N SYK model, the OTOC exhibits exponential growth with a ballistic light cone, i.e. D(r, t) ∼ exp [λ L (t − r/v B )], with v 2 B /λ L exactly equal to the energy diffusion constant. Moreover, in weakly-interacting diffusive metal, D(r, t) ∼ exp [λ L t(1 − (r/v B t) 2 )] with charge diffusion constant D = v 2 B /4λ L . Similar relation between spin diffusion constant and v 2 B /λ L has been deduced numerically in the interacting classical spin-liquid regime of a frustrated spin system [15]. We find the functional form, , with ν varying between 1 to 2 from low to high temperature, to be a good description for the OTOC close to the ballistic chaos front [34] for both easy-plane and easy-axis anisotropies. The relation D = v 2 B /4λ L is violated either qualitatively or quantitatively even at high temperatures. Moreover, we find the evidence of anomalous diffusion at low and intermediate temperatures.
As mentioned in the introduction, there are many earlier studies [19,20,24,25,29] on Lyapunov exponent across phase transitions in classical lattice spin models. In the context of OTOC and spatio-temporal evolution of chaos across phase transitions, more recent results on OTOC in O(N ) models [18,41,42] are directly relevant for our work. As discussed later in detail, the temperature dependence and finite-size scaling of the butterfly speed v B close to the transitions could be related with dynamical critical exponent z [53]. Ref. [18] performed numerical simulation of high-temperature classical dynamics of 2+1d relativistic quantum field theory with O(1) order parameter. Unlike the large-S classical limit in the quantum XXZ model, taking the classical limit of the 2+1d O(1) field theory is somewhat more involved [18,54]. The high-temperature dynamics in 2+1d O(1) model is relevant across finite-temperature 2d Ising phase transition in the model. The dynamics is relevant for 2d transverse field Ising model [55], but, it does not conserve the order parameter. This is unlike the XXZ dynamics in the easy-axis case considered here. The non-conserved order parameter dynamics in the O(1) model falls in the Model C category among the dynamical universality classes [38] and has z = 2 [18,56,57] for 2d Ising transition. In contrast, the dynamics in Eq.(2) conserves the Ising order parameter, i.e. the z component of spin for ∆ > 1, and expected to be in the Model B or D dynamical universality class with z = 4 − η with anomalous exponent η = 0.25 for 2d Ising transition [35]. Refs. [41,42] obtained temperature dependence of λ L and v B in the ordered and disordered phases, and in the quantum critical regime of 2+1d O(N ) model in the large N approximation. The dynamics of the model for N = 2 is more appropriate for 2d quantum rotors and planar antiferromagnets [55] and the finite-temperature transition is expected to have a dynamical exponent z ≈ 2 [35,56]. The large-N approximation, unlike our direct numerical Monte Carlo and spin dynamics simulation in the XXZ model, cannot appropriately describe KT transition in the 2d O(2) model. In contrast, the dynamics [Eq. (2)] in the ferromagnetic XXZ model for the easy-plane case ∆ < 1 is described by the Model E dynamics [35,58,59], where the Poisson bracket terms between planar spin components, i.e. the order parameter, and conserved z component of spin are important. In this case, one expects a dynamical exponent z = 1 in 2d [58,59]. We discuss these points further in the context of our results for v B (T ). Our main results are summarized schematically in a phase diagram in Fig.1. Before describing the results in detail, we give an overview of our main results below. 1. We show that λ L (T ) has a crossover across both KT and Ising transitions, clearly distinguishing low-and high-temperature phases. In particular, we find λ L ∼ T 0.5 and λ L ∼ T 2.5−3 above and below the transitions. 2. The spatio-temporal evolution of the OTOC exhibits ballistic spreading of perturbation in the form of a linear light-cone throughout the temperature range for both easy-plane and easy-axis anisotropies, as shown in Figs.2(a),(b), above and below T KT , for ∆ < 1. Unlike typical quantum systems [34,42,50], we do not find any signature of broadening of the ballistic propagation front of OTOC, even close to the phase transitions. However, we find that there is a delay t 0 in the onset of the light-cone. The time scale t 0 increases with decreasing temperature and seems to diverge for T → 0, like 1/λ L . 3. We find the butterfly speed v B has a non-monotonic temperature dependence, showing a minimum at the transitions. This is the only sharp signature of the phase transition detectable via many-body chaos. 4. Contrary to λ L (T ), sharp signatures of the phase transitions is seen in τ (T ), the time scale extracted from the temporal decay, , above the transition for the auto-correlation function of the planar components of spins. This implies that the chaos time-scales 1/λ L , t 0 are unrelated to the relaxation time τ . We find a power-law decay C xy (t) ∼ 1/t α (α < 1) for the easy-plane case below T KT . 5. We show clear evidence of anomalous diffusion below and close to the transitions for the easy-axis case. We find sub-diffusive to super-diffusive (α > 1) crossover across Ising transition T c for correlation function of the conserved out-of-plane component, The correlation function C zz (t) shows oscillatory behaviour, expected from spin waves, below T KT for the easy-plane case. On the contrary, for both ∆ > 1 and ∆ < 1, C zz (t) always exhibits diffusive behaviour at high temperatures with α ≈ d/2 = 1, as expected for two dimensions (d = 2). We also corroborate the high-temperature diffusive behaviour by computing the dynamical correlation function S zz (q, is Fourier transform of z-component of spins at time t. However, we find that the actual diffusion coefficient D extracted from S zz (q, t) is, in general, either quantitatively or qualitatively different fromD = v 2 B /4λ L . We find spin diffusion constant D D only at infinite temperature for the easy-plane case in the XXZ model.
The above results indicate that there is no qualitative difference between KT and Ising transitions in terms of many-body chaos, at least for the range of anisotropies and temperature we access within our simulations. However, the dynamical spin-spin correlations show qualitatively very different behaviors in the KT and Ising ordered phases, within the time scale over which the perturbation spreads throughout the entire system for the system sizes studied. These imply that, relation between chaos and transport is much more intricate for phases with anomalous diffusion, unlike that in the high-temperature phase well above the transitions, where the diffusive behavior of spin correlation can be linked with the ballistic spread of chaos [14,15].

Results
We study the model Eq.(1) with J = 1 and periodic boundary condition for two values of anisotropy, ∆ = 0.3 (easy plane) and 1.2 (easy axis), for square lattices with N = L 2 sites, with L = 32, 64, 128. We generate 10 4 initial equilibrated configurations at each T via Metropolis Monte Carlo (MC) simulations, and time evolve the configurations via Eq.(2) using fourth-order Runge-Kutta method with time step ∆t = 0.005. As already mentioned, we look into two types of correlation functions -(1) The dynamical spin correlation functions, C xy (t), C zz (t), S zz (q, t), and (2) The classical OTOC of Eq.(3).

Thermodynamics
We first characterize the thermodynamic phases from MC simulations. In particular, we estimate the KT and Ising transition temperatures T KT 0.74 for ∆ = 0.3 and T c 0.96 for ∆ = 1.2, respectively [see Appendix A]. We mainly focus close to the phase transitions and carry out the calculations below and above the transitions for a range of temperatures 0.5 T 2.0, in the KT and Ising-ordered phases as well in the paramagnetic phase. In the easy-plane case (∆ = 0.3), one expects the dynamics in the low-temperature phase to be controlled by gapless spin waves [58,60] which lead to algebraic spatial correlation S r (0)·S 0 (0) ∼ r −η , where the exponent η = T /(2πρ s ) and ρ s the spin stiffness (see Appendix A). The KT transition occurs due to vortex-antivortex unbinding, resulting into a vortex plasma phase for T T KT [35,61,62], where the dynamics is expected to be dictated by the motion of free vortices. We obtain T KT from the universal Nelson-Kosterlitz jump criterion [63] (Fig.9, Appendix A). The statics and dynamics are qualitatively different for easy-axis anisotropy ∆ = 1.2. We obtain the two-dimensional (2d) Ising transition temperature T c from divergence of specific heat and vanishing of the order parameter m z = (1/N ) r S z r (Appendix A). The spin waves of the Ising-ordered phase have a gap ∆ 0 = 4(∆−1) (Appendix A), and the spatial correlation decays exponentially with distance for all temperatures except at T c . Below we investigate how these well-known static and dynamic properties of the model influence transport and chaos in the two cases.

Many-body chaos
We demonstrate the growth and the spread of initial perturbation at r = 0 via D(r = xx, t) for a 1d cut along x direction at T = 2.0 > T KT and T = 0.5 < T KT in Fig.2 (a) and (b), respectively, for ∆ = 0.3. It is evident that, both below and above T KT , the chaos has a ballistic spread like a light cone. As evident from Fig.2(a) for T = 2, and as we have observed even for T close to T KT (not shown), the chaos front across the light cone remains sharply defined and we do not see any evidence for broadening of the front, unlike the diffusive broadening seen for quantum systems with short-range interactions and finite local Hilbert space [34,42,50]. We observe the same phenomena for easy-axis anisotropy ∆ = 1.2 as shown in Fig.10, Appendix B.
To get a better look at the spatio-temporal evolution of the perturbation, we plot D(x, t) as a function of t for a few x in Fig.2(c). Due to the choice of the orthogonal perturbation, D(0, t) starts from zero and initially becomes negative [ Fig.2(d)] over an early-time regime, followed by a power-law growth (linearly with t, not demonstrated) till t 0 , before it starts growing exponentially from a value D(0, t 0 ) ε 2 . As evident from Fig.2(c), the exponential growth ensues at a later time for x = 0. Lyapunov exponent: To quantify spatio-temporal profile of chaos, we define a generalized Lyapunov exponent Using the above, we obtain a light cone from the locus of λ L (x, t) = 0, i.e. where the generalized Lyapunov exponent crosses zero or D(x, t) = ε 2 , as plotted in Figs.2(a),(b). At low temperature T = 0.5 [ Fig.2(b)], the tip of the light cone at x = 0 gets rounded, and, more importantly, shifts to a later time t 0 , compared to that at T = 2.0 [ Fig.2(a)] (also see Fig.3(c)). This clearly suggests a temperature-dependent delay t 0 in the onset of the light cone. We also find similar time scale from (δS x (t)) 2 as shown in Fig.12, Appendix B. As mentioned in Sec.2, the quantity (δS x (t)) 2 starts from ε 2 at x = 0, t = 0. But, it initially decreases with time at x = 0, just like D(0, t) in Figs.2(c),(d).  We extract the Lyapunov exponent λ L (T ) as a function of temperature by fitting D(0, t) ∼ ε 2 e 2λ L t in the exponential growth regime, e.g. in Fig.2(c). The results are shown in Figs. 3(a),(b) across the KT and Ising transitions, respectively, for different system sizes. A smooth crossover around the transitions can be clearly seen indicating a change of temperature dependence of λ L . We find λ L ∼ T 2.86 for T ≤ T KT in the KT phase for ∆ = 0.3 and λ L ∼ T 2.45 in the Ising ordered case for ∆ = 1.2. For the latter, the spin-wave spectrum has a gap ∆ 0 0.8 (Appendix A), and we expect [41,42] an activated T dependence, λ L ∼ e −∆ 0 /T , possibly with a power-law pre-factor, at low temperatures T ∆ 0 . However, for relatively high temperature T 0.5, studied here, we presumably capture only the power-law prefactor ∼ T 2.45 in Fig.3(b). In both the easy-axis and easy-plane cases, λ L (T ) is consistent with a √ T dependence above the transitions, albeit over a limited range of temperature. At very high temperature T 2, λ L eventually saturates to the infinite temperature value ∼ 1 [ Fig.3(a),(b)], which we calculate separately by doing spin-dynamics simulations starting with completely random initial configurations {S ar (0)}. The √ T dependence for the Lyapunov exponent has also been seen for the classical spin liquid phase of a frustrated spin system [15]. The results in Figs. 3(a),(b) indicate no significant effect of L and critical slowing down on λ L , unlike that observed across liquid-gas critical point [30]. However, the results imply that the individual phases can still be distinguished in terms of λ L (T ) [10][11][12][13]. Thus many-body chaos indeed could be an additional tool to characterize dynamics in phases in classical systems and may give new insights not contained within traditional static and dynamical properties. The crossover in chaos across KT transition has been studied earlier [19,20,22], either for smaller system sizes or with a different dynamics. Similar crossover in λ L (T ) has been reported for the Ising transition with various different types of dynamics [18,23,24,29]. Butterfly speed: We next show the temperature dependence of the the butterfly speed v B (T ) in Fig.3(d). The light cones, i.e. the locus of λ L (x, t) = 0, at a few temperatures, e.g., as shown in Fig.3(c), are fitted using t = x/v B + t 0 with v B and t 0 as fitting parameters (see Figs.11(a),(b), Appendix B for more details). The speed v B exhibits a non-monotonic temperature dependence, having a broad minimum around the KT and Ising transition temperatures. A non-monotonic behavior in v B (T ) has been observed in Ref. [18] for the finite-temperature 2d Ising transition in the classical limit of O(1) model. However, in contrast to our results, there v B (T ) shows a maximum at the transition for O(1) model. This implies that chaotic properties are dependent on the details of the dynamics even close to critical points. As discussed in the Appendix.D, one can obtain dynamical scaling laws for OTOC across finite temperature transitions with diverging length and time scales, as in the case of quantum phase transition [53]. Based on these scaling laws, or even just simple dimensional argument [41], v B ∼ ξ/ξ z = ξ 1−z , with dynamical exponent z ≥ 1. Similarly, for ξ L, i.e. close to the transitions, v B ∼ L 1−z , giving the finite-size scaling of v B . As mentioned earlier, z = 1 for the easy-plane case [58,59], thus the weak system size dependence of v B (T ) (∼ L 0 ) in the KT phase (T ≤ T KT ) for ∆ = 0.3 in Fig.3(d) is consistent with the scaling law. However, the same features in v B (T ) are seen around T c for the easy-axis case ∆ = 1.2 [ Fig.3(d)] where one expects z = 4 − η [35,38] and much stronger dependence of v B on |T − T c | and L. This discrepancy could be due to the fact that the easy-axis anisotropy ∆ = 1.2 studied here is not large enough to access true critical regime expected over a narrow range of T around T c . Also, dynamics for such a large z ≈ 4 becomes extremely slow and thus it may be difficult to capture the asymptotic critical dynamics within our simulations times. We note that Ref. [18] also finds very weak temperature and system-size dependence for v B (T ) close to T c in 2d O(1) model, where z = 2 and stronger variations with |T − T c | and L are expected for v B (T ). We keep more detailed analysis of the scaling laws for OTOC and v B for future studies.
At high temperature we find that the quantityD = v 2 B /4λ L to be temperature independent, as discussed later [ Fig.8], suggesting v B ∼ T 0.25 , similar to that found in a classical spin liquid phase [15]. The v B (T ) in Fig. 3(d) suggests faster spread of chaos at low temperatures. This could be due to well-defined spin-wave excitations in the low-temperature KT and Ising-ordered phases. In this regime v B increases at lower temperature whereas λ L → 0 as T → 0, implying a large v B /λ L ( 1, the lattice spacing). This feature may persist even for quantum XXZ model with small S and thus one maybe able to observe [47] the exponential growth in the quantum limit for such a regime dominated by weakly interacting spin waves.
The delay time t 0 extracted from the light cones [ Fig.3(c)] is shown as a function of temperatures in Figs.4(a) and (b). The existence of the delay time and the linear form of the light cone for t > t 0 are further corroborated by plotting λ L (x, t) as a function of x/(t − t 0 ) in Figs.5(a),(b). The λ L (x, t) for different t collapses on a single curve near the light cone λ L (x, t) = 0. We also find that, sans the region deep inside the light cone, D(x, t) can be fitted with a ballistic form ε 2 exp[λ L t{1 − (x/v B (t − t 0 )) ν }] for t > t 0 for both easy-plane and easyaxis cases, as discussed in the Appendix.C. The exponent ν changes from 2 to 1 going from high to low temperatures across the transitions, as shown in Figs.14(b),(c) in Appendix.C. Surprisingly, as shown in Fig.4(a) and (b), t 0 , which characterizes early-time regime prior to exponential growth, roughly follows the temperature dependence of λ −1 L , especially at high temperature. Naively, one would expect t 0 ∼ 1/J ∼ 1, a microscopic time scale. But, t 0 tends to diverge for T → 0 [Figs.4(a),(b)]. Thus interaction effects, that lead to chaotic growth, presumably influence the pre-chaotic non-exponential growth regime of D(x, t) too.
We have also separately computed (δS i x (t)) 2 = (S i ax (t) − S i bx (t)) 2 /2 for i = x, y, z. All the components give the same λ L (t) and v B (T ) and exhibit qualitatively same behaviour, unlike the planar and out-of-plane auto-correlation functions C xy (t) and C zz (t) that we discuss below.

Dynamical spin-spin correlations, diffusion and anomalous diffusion
To understand the possible connection of growth and spread of chaos to transport and dynamical correlations, we first look into C xy (t) and C zz (t) at various temperatures, as shown in Fig.6(a) for ∆ = 0.3. For T > T KT , C zz (t) ∼ 1/t α , i.e. C zz (t) exhibits a power-law decay at long times (t 10). The exponent α ≈ 1 [ Fig.6(b)] is consistent with the expected diffusive behaviour for C zz (t). Above the KT transition, C xy (t) decays exponentially with a time scale τ [see Fig.15 (a), Appendix E]. The evidence of critical slowing down can be observed in τ (T ), as shown in Fig.4(a). Since C xy (t) approaches a finite value C xy (t → ∞) = C ∞ xy (see Appendix E) in the long time limit below the transitions, we plot C xy (t) = C xy (t) − C ∞ xy in Fig.6(a) for T < T KT . C ∞ xy = 0 for T < T KT due to strong finite-size effect [64].C xy (t) shows a power law decay at long times with α < 1, as shown in Fig.6(b). However, we could extract α only close to T KT due to large error in estimating C ∞ xy for T 0.5. Qualitatively, the power-law decay of C xy (t) is expected from non-interacting gapless spin waves in the KT phase giving rise to a temporal spin correlation C xy (t) ∼ t −η [60], implying α η = T /(2πρ s ). However, our results for α close to T KT does not quantitatively match the non-interacting spin-wave results for η(T ) plotted in Fig.6(b). This could be due to coupling between planar and out-of-plane components via spin-wave interactions or finite time ( 100) accessed in our simulations. The asymptotic power law may set in at very long times, as well known for Heisenberg chain [65]. To verify whether a relatively steady power-law regime is reached for C xy (t), we also look into a time-dependent exponent or local logarithmic slope α(t) = d ln( C xy (t))/d ln(t) [ Fig.16(a), Appendix E]. α(t) has a clear drift towards a larger value. This implies that, below T KT , long-time asymptote for C xy (t) is not reached for the time scales over which the chaos spreads ballistically through our finite-sized systems. We note that the transient power-law regime [20 t 80, Fig.16(a), Appendix E], though partially overlaps with the non-exponential growth regime t < t 0 ∼ 5 − 35 [ Fig.4(a)], but extends far beyond the latter and continues deep inside the light cone. Thus the observed power-law (α < 1) regime persists over the time window of the ballistic spread of chaos over the system sizes considered here. We could not extract any power-law exponent for C zz (t) below T KT since it becomes small and oscillatory at low temperatures (not shown).
We did similar calculations of C xy (t) and C zz (t) for the easy-axis case (∆ = 1.2). As in the easy-plane case, C zz (t) [ Fig.15(b), Appendix E] shows a diffusive power law with α 1 at high temperatures (T T c ), as shown in Fig.6(b). However, α decreases approaching the transition indicating a sub-diffusive behaviour. In this regime, the local slope α(t) only shows slight drift with t as shown in Fig.16(b), Appendix E. Below T c , again we obtaiñ C zz (t) = C zz (t) − C ∞ zz [ Fig.15(b), Appendix E] by subtracting the t → ∞ value of C zz (t). In the Ising case, C ∞ zz = 0 below T c because of the symmetry breaking. We find a surprising crossover from sub to super-diffusive scaling of C zz (t), with α(T ) increasing rather sharply from α 1 to a value greater than one below T c , as shown in Fig.6(c). The analysis of α(t) The temperature dependence of α across the transitions. For the easy-plane anisotropy, α ≈ 1 is extracted from C zz (t) above T KT . Below T KT , α < 1 has been extracted fromC xy (t) and α < 1 . For T < T KT , α is compared with the exponent η = T /2πρ s , expected from non-interacting gapless spin waves. For the easy-axis case (c), α has diffusive → sub-diffusive → super-diffusive (α > 1) crossovers from high to low temperature across T c . In this case, α is extracted fromC zz (t). The vertical dotted lines mark the transitions.
[ Fig.16(b), Appendix E] indicates a steady exponent corroborating the super-diffusive power law. The latter may again be an intermediate-time behavior, but it happens over the same times scale over which the chaos spreads in the system. We do not have any good understanding of the anomalous sub and super-diffusive behaviour across the transition and in the Ising ordered phase. For the latter, the phenomena may be arising from some additional conserved mode emerging in the ordered phase and nonlinear coupling between the hydrodynamic modes as happens in one dimensional XXZ model [66]. The planar correlation C xy (t) decays exponentially with the decay time τ (T ) [ Fig.4(b)] for T > T c . However, in contrast to the easy-plane case [ Fig.4(a)], τ (T ) sharply decreases approaching the Ising transition. We note that the planar components are not the critical modes for ∆ > 1 and hence the associated relaxation time does not necessarily need to show critical slowing down. C xy (t) has strongly oscillatory behaviour for T < T c (not shown). Finally, to probe further the high-temperature diffusive phase and the relation betweeñ D = v 2 B /4λ L , which superficially looks like a diffusion constant from dimensional ground, and the actual spin diffusion coefficient D, we compute the dynamical structure factor (or its Fourier transform) For computing the above from spin dynamics simulation, we rewrite above expression as S zz (q, t) = S z q (t)S z −q (0) . Here S z q (t) is Fourier transform of z-component of spins at time t. We obtain S z q (t) from a configuration {S z r (t)} at time t, and average over configurations to obtain S zz (q, t) in Eq.(5). If the system exhibits normal diffusion, for large wavelength q → 0 we expect S zz (q, t) to decay exponentially in time, i.e., e −κ(q)t where the decay rate κ(q) = Dq 2 . We choose momenta q = q x along x direction, close to q = 0 and fit S zz (q, t) with the exponential form to get κ(q) for a given q. Fig.7(a)  S zz (q, t) ∼ e −Dq 2 t is evident, implying diffusive behaviour. The exponential decay of S zz (q, t) as a function of t for small momenta is shown in Fig.7(b) for the same parameter values. We extract the diffusion constant D for a range of temperatures where κ(q) could be fitted via Dq a with a ≈ 2. As we show in Fig.8(b) by plotting a(T ), κ(q) follows the quadratic dependence in the high temperature regime, for both easy-plane and easy-axis cases, where the auto-correlation exponent α 1 [ Figs.6(b),(c)]. An example of near-quadratic dependence of κ(q) is shown for ∆ = 0.3 and T = 2.0 in Fig.7(c). At lower temperatures, we find κ(q) to deviate from the diffusive form and S zz (q, t) to exhibit oscillatory behaviour as a function of both q and t (not shown). The oscillatory behaviour is expected [58] below the transition and even slightly above it, due to spin-waves.
The diffusion constant D calculated for both easy-plane and easy-axis anisotropies in the high-temperature diffusive regime (α 1) is plotted in Fig.8, and compared withD, extracted from Fig.3. In the easy-axis case, both D andD are independent of temperature at high temperatures (1 T /T c 2) and closely approach their infinite-temperature values, however, D D . In the easy-plane case, D varies substantially with temperature even at high temperature and slowly approaches its infinite-temperature value D ∞ . This behaviour is unlike that ofD, which varies little with temperature and coincides with its infinite-temperature valueD ∞ . Moreover, we findD ∞ D ∞ for the easy-plane case. Nevertheless, our results suggest thatD is quite distinct from the actual diffusion constant D in general, unlike that in a correlated classical spin-liquid state [15] or in quantum systems like strongly or weakly interacting diffusive metals [7,39,43,44]. For the strongly interacting metals,D is related to the energy diffusion constant implying that the chaos or scarambling directly controls the thermal diffusion [7,39]. To this end, our results raise interesting questions about actual physical process governingD in the semiclassical limit of 2d XXZ model. For example, it would be interesting to compute the energy diffusion constant for the spin model in future and see whether it corresponds more closely toD, rather than the spin diffusion constant that we calculate here.  (Fig.6). Note that for simplicity of notation, we refer to D as diffusion constant even when a deviates substantially from the diffusive value 2.

Conclusion
We have studied here the OTOC, and the dynamical spin correlations in the semiclassical limit of the 2d XXZ model. In particular, we have tuned the anisotropy in the model to study the dynamical properties across KT and 2d Ising transitions via simulation of spin precession dynamics as a function of temperature. Thus we obtain a dynamical phase diagram in terms of chaos and spatio-temporal spin correlations for the classical 2d XXZ model. We have computed temperature dependence of the Lyapunov exponent and butterfly speed, which show crossover across the transitions and no effect of critical slowing down. Only relatively sharp signature of the transitions is exhibited by a non-monotonic temperature dependence of butterfly speed having a minimum at the transition.
Overall, we find chaotic growth and spread, and the dynamical correlations above the transitions at high temperature in the easy-plane and easy-axes cases very similar. However, the dynamical spin-spin correlations are very different at intermediate times in the low-temperature KT and Ising ordered phases, and close to the transitions, although chaos still spreads ballistically in these regimes. This leads us back to the question on the connection of chaos and transport. A simple, albeit heuristic, way [67] to obtain ballistic light-cone for chaos from diffusive transport is to take the separable ansatz D(x, t) ≈ ε 2 e λ L t C zz (x, t) for the OTOC and plug in the diffusive form C zz (x, t) ∼ e −x 2 /4Dt /t d/2 . This leads to a velocity-dependent generalized Lyapunov exponent [14,34,36] λ L (x, t) = λ L 1 − (x/v B t) 2 , and naturally gives rise to the relation D ∼D = v 2 B /4λ L . However, we findD to be very different from D, except for the easy-plane case at infinite temperature. More importantly, the simple ansatz clearly fails in the phases exhibiting anomalous diffusion, like Ising ordered phases in the XXZ model. Our observation of the anomalous diffusion over a large range of temperature for easy-axis anisotropy in the semiclassical limit of 2d XXZ model is intriguing and it would be good to get a proper understanding of these phenomena and their possible connections to chaos.
We reveal an early-time pre-exponential regime in the form of a temperature-dependent overall delay t 0 (T ) in the onset of the light cone. The time scale t 0 is presumably connected with chaos time scale and originates from the same many-body interactions that give rise to chaotic growth. This result suggests the possibility of extracting new dynamical regime with a suitable choice of the correlation function. Such a regime may even exist for non-chaotic systems, e.g. integrable or fully quantum ones, where there is no exponential growth inside the light cone [34].
In the absence of any good analytical understanding of many-body chaos in classical systems, our results call for the development of a theoretical framework [68] to compute the classical OTOC along the line of that done for quantum systems, in a suitable large-N limit [16], or in some perturbative regime [36,68], like at low temperatures with weakly interacting spin waves. Such a theory may give rise to new insights into the dynamics of interacting classical systems as well as quantum systems in the semiclassical limit. It would be desirable to obtain hydrodynamic description of the OTOC in such classical spin systems with Hamiltonian dynamics or some related tractable toy models, e.g. with random classical Liouvillian dynamics, along the line of those developed for quantum systems [43,[48][49][50][51].

A.1 Spin stiffness and KT transition
To obtain the Kosterlitz-Thouless (KT) transition temperature T KT for ∆ = 0.3, we calculate the spin stiffness ρ s , which measures the rigidity of the spin configuration to small twist or rotation of the spins along some direction. For the KT transition the relevant spin stiffness is obtained by twisting the planar components (S x r , S y r ) = (S r cos φ r , S r sin φ r ) of the spins. The spin stiffness for the model of Eq.1 in the main text can be obtained as We calculate ρ s (T ) via MC simulations and apply the Nelson-Kosterlitz universal jump criterion [63] ρ s (T KT )/T KT = 2/π to obtain T KT , as shown in Fig. 9(a).

A.2 Two dimensional (2d) Ising transition
To estimate the 2d Ising transition temperature T c for ∆ = 1.2, we calculate the magnetization and the specific heat per site T c is obtained from the divergence of c V shown in Fig. 9(b), as well as, from m z (T ) which continuously goes to zero at T c [ Fig. 9(b)(inset)].

Easy-plane anisotropy
In this case, the spin-wave dispersion is obtained by expanding the dynamical equations [Eq.2, main text] around the T = 0 ground state, which corresponds to all spins aligned along a direction (sayx) in the xy-plane. In this case, S x r 1 S y r , S z r , and Eq.2 (main text) gets reduced to Using Fourier transformation S i r (t) → S i q (ω), it is straightforward to obtain the spin-wave dispersion for q → 0 for the square lattice Hence the spin waves are gapless in the easy-plane case.

Easy-axis anisotropy
In this case, the all the spins can be taken to be aligned alongẑ direction, and S z r 1 S x r , S y r . Following the same method as the easy-plane case, we obtain the spin-wave dispersion with a spin-wave gap ∆ 0 = 4J(∆ − 1).

B Classical OTOC
We show the OTOC D(x, t) for a cut along x direction for ∆ = 1.2 at two temperatures, T = 1.6 > T c and T = 0.74 < T c in Fig.10 ballistically, as indicated by the light cones λ L (x, t) = 0, and the light cones start at a finite time t 0 , which increases with decreasing temperature. We extract the light cones, e.g. in Fig.3(c) (main text), by finding the locus t D 0 (x) of λ L (x, t) = 0 or D(x, t) = D 0 = ε 2 , as shown in Fig.11

B.1 Trajectory divergence and decorrelation
In all our calculations and results shown in the main text, we consider the decorrelation function [Eq.(3)] or the classical OTOC, which essentially measures how uncorrelated two spin configurations a and b are. As discussed earlier, the copy b differs from a initially only at a single site (r = 0), namely S br (0) = S ar (0) + ε(n × S a0 )δ r,0 . One can also consider the trajectory divergence, which is defined (δS r (t)) 2 = (S br (t) − S ar (t)) 2 . The decorrelation function and trajectory divergence differ at O(ε 2 ) initially at r = 0, i.e.
This initial-time difference eventually relaxes in the intermediate time regime where chaos sets in and we get the same exponential growth and ballistic spread characterized by λ L and v B , respectively, from both of these quantities.
To compare the decorrelation function and the trajectory divergence, we plot (δS x (t)) 2 in Fig.12 as a function of t for a few x, for T = 2.00, 0.50 and ∆ = 0.3. (δS x (t)) 2 exhibits behaviour very similar to D(x, t) [ Fig.2(c),(d)]. By definition, (δS 0 (t)) 2 starts from ε 2 and decreases over an early-time regime, followed by a power-law growth till t 0 , before it starts growing exponentially from a value D(0, t 0 ) ε 2 . The exponential growth occurs at a later time for x = 0. As shown in Fig. 13  C Velocity-dependent Lyapunov exponent As shown in Fig.5(a),(b), the generalized Lyapunov exponent λ L (x, t) [Eq.(4)] at different t can be collapsed into a single curve as a function of a velocity v = x/(t − t 0 ) for t > t 0 over a relatively large range around v = v B for both outside (v > v B ) and inside (v < v B ) the light cone. However, this range shrinks progressively with decreasing temperature. Over this range, the velocity-dependent Lyapunov exponent λ L (v) can be fitted well with a ballistic form, as shown, for example, in Fig.14(a) for easy-plane anisotropy ∆ = 0.3 at T = 0.96. The deviation from the scaling for v 1.5v B in the inset of Fig.14(a) and in Fig.5 is due to the numerical precision. The extracted values of the exponent ν are plotted as a function of T for ∆ = 0.3 in Fig.14(b). Here to obtain ν we have used a linear approximation λ L (v) λ L ν(1 − v/v B ) to the non-linear fitting form in Eq. (11) for v ∼ v B and fitted with ν as the fitting parameter for fixed values of λ L and v B obtained from Fig.3. The resulting linear fit and the non-linear function are compared with the data for ∆ = 0.3, T = 0.96 in Fig.14(a). We find ν ≈ 1.9 at high temperature, but ν decreases towards ∼ 1 with decreasing temperature T 0.6. We could not get a reliable goodness of fit at lower temperature since the fitting range shrinks substantially for both v > v B and v < v B for T 0.6.
To assess the goodness of the fits we obtain the errorbars [ Fig.14 generated with different initial conditions at temperature T into multiple groups. Based on these errorbars on D(x, t) and χ 2 fitting of the data for λ L (v) with the linear approximation to Eq.(11) mentioned above, we obtain the error in ν, and estimate the goodness of fit using [69] where N data is number of data points over the fitting range. A healthy fit is defined as 0.01 P P max , where P max is slightly less that 1 [69]. The errors in the estimate of ν and P values are indicated in Fig.14(b) for the easy-plane case ∆ = 0.3. We have not carried out detail error analysis for easy-axis anisotropy ∆ = 1.2, but ν in this case is shown as a function T in Fig.14(c).
To verify the possibility of the broadening of the chaos front around v ≈ v B , we have also tried to fit [34,42] with p as the fitting parameter. We find p 0, consistent with the absence of broadening [34] and the fact that the λ L (v) is more or less linear for v > v B [Figs.5,14(a)], over the range of v where the scaling collapse works.

D Dynamical scaling law for OTOC
One can obtain dynamical scaling laws for OTOC and the butterfly speed v B across finitetemperature phase transitions with diverging length scale, as in the case of quantum critical point (QCP) [53]. To this end we consider F(r, t) = 1 − D(r, t) = S ar (t) · S br (t) .
We can write a scaling form for F(r, t) by applying scaling transformation F(r, where z is the dynamical exponent, Φ a universal scaling function, and ξ is the correlation length that diverges for T → T c for Ising transition in the easy-axis case and for T ≤ T KT at the KT transition and KT phase for easy-plane anisotropy. Since, F(r, 0) = 1 for any L, r and ξ, ∆ F = 0. Choosing b = ξ, we obtain F(r, t) = Φ(L/ξ, r/ξ, ξ −z t).

E Spin auto-correlation functions
We calculate the spin auto-correlation functions C ii (t) = (1/N ) r S i r (t)S i r (0) from the spin dynamics simulations starting with thermal initial conditions at different temperatures, as discussed in the main text. We mainly look into the planar correlation C xy ≡ C xx + C yy and out-of-plane correlation C zz . More specifically, we computẽ  where C ∞ ii = C ii (t → ∞) = (1/N ) r S i r (0) 2 , as we have S i r (t → ∞)S i r (0) = S i r (t → ∞) S i r (0) = S i r (0) 2 for thermal initial conditions. As shown in Fig.9(a)(inset), in the easyplane case, C ∞ xy = 0 below T KT since S x/y i (0) = 0 for strong finite-size effect [64]. In the easy-axis case, C ∞ zz = 0 below T c due to spontaneous symmetry breaking [ Fig.9(b)(inset)]. In Fig.15(a), we show that C xy (t) decays exponentially with a relaxation time τ (T ) for T > T KT and ∆ = 0.3. Similar exponential decay is observed for the easy-axis case ∆ = 1.2. However, the relaxation time τ (T ) increases approaching the transition for ∆ = 0.3, whereas it decreases for ∆ = 1.2, as shown in Figs.4(a),(b) in the main text. We show [ Fig.15(b)] thatC zz (t) exhibits a power-law decay for t 5 across the Ising transition. The exponent α changes from a diffusive value 1 at high temperature to sub-diffusive values (< 1) close to T c , and finally to super-diffusive values (> 1) at low temperature [ Fig.6(c)]. To verify whether C xy (t) and C zz (t) have attained a steady power-law behaviour within our finite simulation time ( 100), we obtain a time-dependent exponent α(t) = d ln(C(t))/d ln t (C = C xy , C zz ) for T < T KT in the easy-plane case [ Fig.16(a)] and across T c for the easy-axis case [ Fig.16(b)]. These suggest that a steady power-law exponent is achieve for ∆ = 1.2 except at the lowest temperature studied, whereas the exponent shows perceptible drift towards a larger value for ∆ = 0.3.