Discrete Laplacian Thermostat for Spin Systems with Conserved Dynamics

A well-established numerical technique to study the dynamics of spin systems in which symmetries and conservation laws play an important role is to microcanonically integrate their reversible equations of motion, obtaining thermalization through initial conditions drawn with the canonical distribution. In order to achieve a more realistic relaxation of the magnetic energy, numerically expensive methods that explicitly couple the spins to the underlying lattice are normally employed. Here we introduce a stochastic conservative thermostat that relaxes the magnetic energy while preserving the constant of motions, thus turning microcanonical spin dynamics into a conservative canonical dynamics, without actually simulating the lattice. We test the thermostat on the Heisenberg antiferromagnet in d=3 and show that the method reproduces the exact values of the static and dynamic critical exponents, while in the low-temperature phase it yields the correct spin wave phenomenology. Finally, we demonstrate that the relaxation coefficient of the new thermostat is quantitatively connected to the microscopic parameters of the spin-lattice coupling.


I. INTRODUCTION
The impact of symmetries and conservation laws on the dynamics of physical systems cannot be overstated, and spin systems are no exception. Not only conservation laws can change the low-temperature dispersion relations, but they can also radically change the dynamical critical exponents [1]. The most effective method to numerically study spin systems with symmetries and conservation laws is to microcanonically integrate the reversible equations of motion [2,3], a technique called Spin Dynamics (SD) by Landau and coworkers, who advanced it very significantly [4][5][6][7][8][9][10]. As the energy is conserved, in order to thermalize the system SD draws the initial conditions from a canonical ensemble at temperature T using Montecarlo. Although SD provides excellent results, one can raise an issue, which is both conceptual and practical. Consider a microcanonical simulation of a particles system, as in standard Molecular Dynamics (MD); despite the inevitable simplifications, one can argue that MD is conceptually the same as the actual physical dynamics of an isolated system. On the other hand, microcanonical SD has not quite the same conceptual standing as microcanonical MD: in an actual spin system, where the spins belong to the atoms of an underlying lattice, thermal relaxation occurs mostly through the exchange of energy (but not of magnetization) between the spins and the nuclei; by excluding the lattice from the simulation, and including the temperature only through the initial conditions, SD takes a (clever) shortcut, which has however no actual physical counterpart, as in most physical systems it is quite hard to isolate the spins from the lattice. While within MD one can consider a subsystem A as the heat bath acting upon an adjacent subsystem B, in most spin systems the heat bath is provided by the lattice, which is instead absent in SD. This issue has also practical implications; for example, if we want to change the temperature during a spin simulation or if we want to study the effects of a quench, SD has a problem, as T is fixed once and for all at the beginning of the simulation by the initial conditions. It was precisely to deal with this problem that Tauber and Nandi devised an interesting hybrid method, in which microcanonical SD is alternated with canonical Kawasaki Montecarlo (KM) at temperature T , giving rise to a SD-KM-SD-KM-. . . dynamical sequence [11,12]. Notice that although KM is conservative, it is not a true dynamics, as conservation is achieved by swap moves, rather than being dynamically generated by the symmetry through Noether's theorem; which is precisely why KM needs to take turns with SD in the method of [11,12].
A more fundamental approach is to explicitly take into consideration the interaction between spins and lattice by running in parallel a spin dynamics and a molecular dynamics simulation, an approach that we will call SD+MD [13,14]. Even though in this microcanonical dynamics the total energy is conserved, there is energy exchange between spin and lattice, so that the magnetic energy relaxes. SD+MD is the most complete and realistic numerical method to simulate magnetic systems, but it is also significantly more expensive than SD from a computational point of view, as it needs to update the positions and the momenta of the nuclei, in addition to the spins. It would be useful to have a method as simple and economic as SD, but which includes the relaxational effects of the spin-lattice coupling. We note that such methods exist for the case where interactions do not conserve the total magnetization. The Landau-Lifshitz-Gilbert equation [15] and other Langevin-type equations [16,17] relax the spin energy through multiple dissipative terms, violating the conservation law. However, analogous mechanisms for the conservative case, where the total magnetization is constant, are still lacking. Here, we fill that gap and present a new stochastic thermostat that turns microcanonical spin dynamics into a conservative canonical dynamics: the novel numerical method updates only the spins, hence having a low computational cost -similar to SD -and yet it thermalizes the magnetic energy, as if the spins were coupled to an underlying lattice. arXiv:2212.09647v2 [cond-mat.stat-mech] 26 Jun 2023 The work is organized as follows. In Section II, we introduce the new conservative thermostat. In Section III, we test the new method on the Heisenberg antiferromagnet and obtain the correct critical exponents, spin wave dispersion law and spin wave softening. Finally, we show in Section IV that the relaxation coefficient of the stochastic thermostat can be qualitatively and quantitatively connected to the microscopic parameters of the spin-lattice coupling.

II. DISCRETE LAPLACIAN THERMOSTAT
In this Section, partly taking inspiration from the dynamical mesoscopic equations of conserved fields, we will devise a way to write a conservative thermostat for a discrete spin dynamics. As we shall see, the key objects to achieve this result will be the discrete Laplacian and the incidence matrix; for a discussion of both these quantities in the context of graph theory see [18].

A. Microcanonical Spin Dynamics
We consider a system of N spins, σ µ i , with i = 1, . . . , N and µ = 1, . . . , d, obeying the Poisson brackets, where ε µνρ is the Levi-Civita antisymmetric tensor. Microcanonical SD amounts to integrate the reversible equations of motion, which naturally conserve the energy H. We consider the case in which also the total magnetization, is a constant of motion, which, of course, amounts to require, As we discussed in the Introduction, the only way to thermalize a microcanonical SD simulation is to start the numerical experiment from an initial condition previously thermalized at temperature T with a canonical non-conservative dynamics, typically Montecarlo [9]. We want to change that; our aim is to add to the microcanonical SD dynamics (2) some new irreversible relaxational terms, i.e. a stochastic thermostat, which relaxes the energy H, while at the same time conserving the magnetization M .

B. Inspiration from dynamical field theory
Within dynamical field theory there is a standard method to achieve the equilibration of a field φ(x, t), subject to the constraint that its space integral is a constant of motion; this method consists in adding to the reversible parts of the dynamics, an irreversible relaxational force and a stochastic noise linked to each other by a kinetic coefficient proportional to the Laplacian [1]. More precisely, we can write, where the noise correlator satisfies the equilibrium condition, and where (crucially) the kinetic coefficient Γ is given by, The positive parameter λ is usually called transport coefficient; notice that also the kinetic coefficient Γ is positive, as the continuous Laplacian is a negative-definite operator. To see why this method works, it is sufficient to go to Fourier space, where −λ∇ 2 → λk 2 , so that the space integral of the field -namely the mode φ(k, t) at k = 0is automatically conserved, both by the relaxational force and by the noise. Using the standard terminology of Hohenberg and Halperin [1], this method is used in Model B (spinodal decomposition), in Models E and F (superfluid helium), in Model G (quantum antiferromagnet), and in Model J (isotropic quantum ferromagnet). Moreover, a conserved noise with the form of (7)-(8) is used in the context of non-equilibrium theories, in particular in the case of the conserved KPZ equation [19]. We want to take inspiration from this mesoscopic continuous case to devise a conservative relaxational dynamics that works also in the microscopic discrete case. We stress however that we are not discretizing equations (6)-(8) in any concrete way; we will simply try and export the physical mechanism of conservation used in the continuous case to the discrete setup.

C. The role of the discrete Laplacian
Transposing to the discrete case the idea behind the mesoscopic conservative dynamics (6)- (8) does not seem too difficult, given that there exists a very well-known discrete version of the Laplacian operator; this matrix, that we shall call Λ ij , is defined in the following way [18], where n ij is the adjacency matrix defining the lattice's topological structure: n ij = 1 if two sites interact with each other, n ij = 0 otherwise; we will assume a symmetric interaction network, so that also the Laplacian is symmetric. Notice that -at variance with its continuous counterpart -the discrete Laplacian is a positive-definite matrix, i.e. Λ ∼ −∇ 2 (to make this correspondence dimensionally consistent we should include the square of the lattice spacing; however -as we have already remarkedwe are not pursuing an actual discretization of the continuous case, but just using it as a conceptual guideline). One of the crucial features of the Laplacian is the fact that it has a zero mode, namely, We can exploit this property and directly mimic equations (6)- (8), so to achieve a relaxational dynamics of the discrete spins, which at the same time enforces the conservation law. We propose to do this by writing the following canonical stochastic equations, where -in order to achieve thermal equilibrium -the dimensionless relaxation coefficient λ and the noise ξ i must satisfy the Fluctuation-Dissipation (FD) relation, Let us show that dynamics (11)-(12) conserves the total magnetization; we recall that {H, M } = 0 and that but from (12) we see that the random variable i ξ i has zero mean and zero variance, hence it must be identically zero for each realization of ξ i , finally proving that, On the other hand, the new irreversible relaxational term proportional to the 'force', ∂ σ H, pushes the spins to relax towards the configuration that minimizes the Hamiltonian, so that the energy converges to its equilibrium value at temperature T . Indeed, a standard Fokker-Planck argument [20] shows that the equilibrium probability distribution generated by (11)-(12) is the Gibbs-Boltzmann canonical ensemble, FIG. 1. Schematic view of the adjacency matrix n, incidence matrix D, and discrete Laplacian Λ in a very simple lattice.
We therefore have a new canonical stochastic dynamics, in which the microscopic spins are thermalized at temperature T and the total magnetization is conserved. And yet, we still have some more work to do.

D. Incidence matrix noise
Up to now the method simply mimics the continuous case, but of course the real problem is how to produce a noise ξ i whose correlator is proportional to the discrete Laplacian, Λ ij , as required by equation (12). Let us see how we can solve this problem.

The complicated way
If we insist working exclusively in the space of sites, the first possibility that comes to mind is to try and find the matrix D ij whose square (over the sites) is the Laplacian, and then define the conservative noise as, ξ i = k D ik ϵ k , with ⟨ϵ k ϵ l ⟩ ∼ δ kl , hence giving ⟨ξ i ξ j ⟩ ∼ Λ ij , as desired.
Although at first sight this seems feasible, it is in fact a dead end. The problem with this method is that solving equation (16) is far from straightforward: the matrix D ij heavily depends on the specific nature of the lattice, because one needs to go through the explicit form of the eigenvalues and eigenvectors of the discrete Laplacian to find it; in fact, one finds, where w q i are the eigenvectors of the discrete Laplacian matrix Λ, with λ q the corresponding eigenvalues. This form of D ij is problematic, because the spectrum of the Laplacian is known analytically only for a limited number of regular lattices, and even in these cases only with periodic boundary conditions, while we would like to have a method that works irrespective of the specific topology of the lattice, and of its boundary conditions. Moreover, even in those cases where the Laplacian spectrum is exactly known and the matrix D ij can be calculated, its mathematical expression is extremely cumbersome, even for the simplest lattices. An even more serious problem is that of locality: in general, the matrix D ij in (17) is non-zero even for non-interacting sites, namely for pairs of sites for which the adjacency matrix is zero, n ij = 0; this means that the conserved noise, ξ i = k D ik ϵ k , connects sites that were not directly interacting in the original Hamiltonian, which seems unnatural, to say the least.
We want to develop a method that is local and that employs noting more complicated than the plain adjacency matrix itself, n ij , and that possibly entails no calculations what-so-ever, neither hard, nor easy. Once again, field theory comes to our rescue.

The simple way
In the continuous case, the fact that the noise variance is proportional to the Laplacian suggests that, in some way, the noise must be proportional to a gradient, ξ ∼ ∇. Fortunately, a simple discrete equivalent of the gradient does exist in graph theory: interestingly, it is a matrix defined in the space of sites and links, rather than of sites only. Let us label the sites of the lattice with {i, j, . . . } and the links with {a, b, . . . }. The incidence matrix, D ia , is constructed as follows [18]: after arbitrarily assigning a direction to each link a, we set D ia = +1 if i is at the end of a, D ia = −1 if i is at the origin of a, and D ia = 0 if site i does not belong to a (Fig.1); note that, by construction, we have, A little reflection immediately shows in what sense the incidence matrix is the discrete equivalent of the gradient: the 'derivative' of a discrete set of variables {σ i } over link a can now be written as, in view of which, relation (18) simply expresses the obvious, i.e. that the derivative of a constant is zero. Notice also that the arbitrariness in the definition of D ia due to the arbitrary choice of the directions of the links, reflects the inevitable arbitrariness in the definition of the derivative on a general discrete lattice.
We can now state the crucial property of the incidence matrix, namely that its square over the links is equal to the discrete Laplacian [18], What we have to do now seems clear: we need to switch from a noise defined on the sites, to a noise defined on the links. More precisely, on each link a we define a standard δ-correlated Gaussian noise, ϵ a , with variance, so that the site noise acting on each spin i can finally be constructed as, which gives discrete flesh to the idea that the conserved noise is proportional to a gradient, ξ ∼ ∇. Let us compute the variance of this new noise, so that we recover exactly the desired expression, equation (12). Moreover, we see from equation (22) that, within this construction, the site noise, ξ i , is the sum of all the link noises, ϵ a , incident on that site; because by construction i D ia = 0, from (22) we have that, which makes it even more apparent the fact that the new noise conserves the total magnetization in (11).

E. Summary of the Discrete Laplacian Thermostat
We have finally derived a closed set of equations specifying the canonical stochastic dynamics of a spin system with conserved magnetization. Because of the rather lengthy derivation, we summarize here the new canonical equations, where Λ ij is the discrete Laplacian and D ia is the incidence matrix. We call this method Discrete Laplacian Thermostat (DLT). It is important to stress that for λ = 0 the DLT canonical dynamics (25)- (27) becomes identical to microcanonical SD (2), where the energy does not relax; hence, we expect the relaxation coefficient λ to be related to the inverse of the time-scale of energy thermalization. Apart from this, we will show that the relaxation coefficient does not impact on the qualitative features of the system: both the static and dynamic universality classes are unchanged, and the classic spin-wave phenomenology is correctly reproduced by DLT.
We conclude this Section with a notational clarification. After equations (6), (7) and (8), one could have expected us to call λ the 'transport coefficient', as in the field theory context [1]. However, the terminology 'transport coefficient', as well as 'kinetic coefficient', belongs to the very specific context of hydrodynamics, which is a coarse-grained theory. Here, on the other hand, we are dealing with microscopic dynamic equations. Moreover, under coarse-graining, the microscopic spins may in general give rise to both conserved and non-conserved hydrodynamic fields, depending on the specific model [1,21], so that the microscopic parameter λ could contribute at the coarse-grained level both to the transport coefficient of a conserved field and to the kinetic coefficient of a non-conserved field. Therefore, we prefer to adopt a more neutral name for the microscopic parameter λ, and call it 'relaxation' coefficient.

III. TESTING DLT IN THE HEISENBERG ANTIFERROMAGNET
We numerically test DLT on the classical Heisenberg antiferromagnet. The Hamiltonian is, where J > 0 and n kl corresponds to a d = 3 cubic lattice of side L with PBC. The Hamiltonian is rotationally invariant in the absence of an external field, so that {H, M } = 0 and the global magnetization is conserved, M = 0. The order parameter is the non-conserved staggered magnetization, Ψ = i π i σ i , where π i = ±1 is the parity of site i. By plugging Hamiltonian (28) into the DLT equation (25), and after using Poisson's relations (1), we obtain, where ξ i is given by (26) and (27). In order to fix the norm of the spins one could use a Lagrange multiplier, which would however slow down the simulation; instead, we use a single-particle potential suppressing norm fluctuations (see Appendix A1 for more details). We have performed simulations of system sizes from L 3 = 216 up to L 3 = 8000 (the lattice spacing is ℓ = 1). We work in units such that k B = 1 and ℏ = 1; we also choose units such that J = 1, hence we are effectively measuring time in units of J −1 .

A. Numerical integrator
We need to make an important technical remark here. Standard, microcanonical SD employs a non-stochastic fourth-order Runge-Kutta integrator with time-step ∆t = 5·10 −4 ; of course, in any reversible microcanonical dynamics, in which energy must be conserved, it is important that the integrator is highly accurate, lest the conservation of energy is violated, which does not bode well for a microcanonical dynamics. But if we add a stochastic thermostat to the dynamics, the energy is no longer conserved, so that the tiny inaccuracies in the integration of the reversible part, which would cause a violation of energy conservation in the microcanonical case, become now irrelevant compared to the relatively huge fluctuations of the energy caused by the irreversible stochastic term; hence, what would normally happen when we switch from a non-stochastic microcanonical dynamics to a stochastic canonical one, is that we should also switch from a nonstochastic highly accurate integrator to a stochastic one.
But in our case, we want to be able to precisely recover the microcanonical SD dynamics in the limit λ → 0; while the case at exactly λ = 0 could be studied by switching back to a non-stochastic integrator, this is not possible for small values of λ, when inaccuracies in the integration of the reversible term are not irrelevant compared to the energy fluctuations caused by the irreversible stochastic term. Therefore, the deterministic term must still be integrated with high accuracy, to correctly reproduce also the case of small relaxation coefficient. This is the reason why, even though we are dealing with a stochastic differential equation, we employ the same non-stochastic fourth-order Runge-Kutta integrator as in standard SD.
On the other hand, by construction, both the reversible term and the irreversible thermostat satisfy exactly the constraintṀ = 0 (see Sections II C and II D), independently from the accuracy of the integrator, simply thanks to the antisymmetric form of the equations; hence, the conservation law of the magnetization is immune from all this.

B. Conservation, transition, susceptibility
In Fig.2a we show that DLT conserves the magnetization with very high precision, while relaxing the energy, both when starting from random initial conditions, Ψ ∼ 0, M ∼ 0, and when starting from ordered initial conditions, |Ψ| = N , M = 0; to avoid slowing down due to coarsening we used ordered initial conditions in the rest of our study. Fig.2b shows that energy thermalization is quicker the larger is the relaxation coefficient, λ; the energy is not a critical variable, hence its thermalization time, τ E , is a harmless microscopic scale of the system; we find τ E ∼ λ −1 , hence confirming the expectation that the relaxation coefficient is the inverse of the energy time scale (had we not set J = 1, we would have τ E ∼ (Jλ) −1 ).
The antiferromagnet has a continuous phase transition at T c = 1.446 [22,23], which DLT captures well (Fig.2c). Notice that the low-T linear behaviour of the modulus of the staggered magnetization, 1 − Ψ/N ∼ T , predicted by the theory [24], is also correctly reproduced by DLT.
We test the static critical behaviour by studying the susceptibility, which satisfies the finite-size scaling relation, χ = ξ γ/ν g(L/ξ), where g(x) is a scaling function; we probe the scale-free regime by selecting at each size L the temperature T c (L) at which χ is maximal (see Appendix A2 for a detailed description of the procedure), hence we obtain ξ ∼ L, and therefore χ ∼ L γ/ν ; the fit to the theoretical exponent [25] is quite satisfactory (Fig.2d).

C. Spin waves and their dispersion relation
The low temperature regime of antiferromagnets is dominated by spin waves [26], which are the quintessential consequence of the symmetry and conservation law; hence it is important to check how DLT performs in relation to them. The transverse scattering function S T (k, ω) of the staggered magnetization -which is computed following [7] -is reported in Fig.2e: as expected, above the critical temperature there is a simple diffusive peak, while two spin-wave peaks at ±ω c emerge for T < T c .
In the spin wave phase, the characteristic frequency ω c depends on the wavevector k according to the exact dispersion relation (see Appendix B for the full derivation), which is fairly well reproduced by DLT, considering that no fitting parameters what-so-ever are used (see Fig.2f).

D. Critical dynamics and the exponent z
As fundamental as the existence of spin waves is the emergence of critical slowing down at the phase transition: in the bulk, the relaxation time of the order parameter at k = 0 diverges at the critical point as a power law of the correlation length, τ ∼ ξ z [1]. An exact renormalization group calculation of the dynamic critical exponent yields z = 1.5 for the Heisenberg antiferromagnet in d = 3 [27], a result that has been confirmed both by numerical simulations [10] and by experiments on perovskites [28,29].
Notice that z = 1.5 differs very significantly from the value z ≈ 2 of the universality class of standard nonconserved ferromagnets (as the Ising model), also called Model A class in the classification of Hohenberg-Halperin [1]. The profound reason for this difference is the connection between symmetry and conservation law: if the critical dynamics of the antiferromagnet is studied through a Montecarlo method, one obtains z ≈ 2; interestingly, it is not only standard non-conserved Montecarlo that fails to yield the correct dynamical critical exponent [30], but also Kawasaki Montecarlo (KM) fails: although KM conserves M , this conservation has no effects on the dynamics of Ψ, hence giving z ≈ 2; as we have already noted, it is not the conservation per se, but the deep dynamical connection between symmetry and conservation law that yields the correct universality class.
To calculate the critical exponent z we use dynamical scaling [31], according to which the relaxation time has the form, where g(x) is a scaling function and all the dependence on the temperature T goes into ξ(T ) (the precise definition of relaxation time is reported in Appendix A3). We now define the largest relaxation time, τ , as that corresponding to the lowest mode, k = 2π L ; from (31) we obtain, τ = L z g(ξ/L) .
If for each size L we work in the scale-free regime, namely at the pseudo-critical temperature, T c (L), where χ is maximal (see Section III B and Appendix A2), we have which is the relation we test, using sizes ranging from L 3 = 1000 to L 3 = 8000. The result of DLT is quite satisfactory (Fig.3): we find z = 1.47 ± 0.07 for λ = 0.1, z = 1.56±0.08 for λ = 0.2, and z = 1.50±0.10 for λ = 1.0, all values consistent with the exact result, z = 1.5. The tests of thermalization that we used to check that all simulations are actually at equilibrium are reported in Appendix A4. Hence, the critical exponent z does not depend on the relaxation coefficient λ, which therefore does not impact on the dynamic universality class of the model; the role of the relaxation coefficient is simply to shift the prefactor of τ , which is reasonable, given the identification of λ −1 as a non-critical time scale of the system.

E. Spin wave softening and damping
A feature over which the relaxation coefficient does have a nontrivial impact is the blunting of spin waves (Fig.4): the larger is λ, the weaker is the spin wave peak (SW damping [32]) and the lower its characteristic frequency (SW softening [33]).
Damping and softening are typically caused by the coupling of the magnetic degrees of freedom (magnons) with the lattice vibrational degrees of freedom (phonons) [34], an interpretation confirmed by numerical simulations [35,36]. Hence, the fact that λ is connected to softening suggests to try and link the relaxation coefficient of DLT to the microscopic details of the magnon-phonon system. This is what we do in the next Section.

IV. CONNECTION BETWEEN DLT AND SPIN-LATTICE DYNAMICS
The canonical DLT method has just one extra parameter compared to microcanonical SD, namely the relaxation coefficient λ, which -as we have seen -is essentially the inverse of the energy relaxation time; this makes sense, as in SD the energy does not relax, of course, and the limit λ → 0 reproduces exactly the SD case. We also have found from the numerical simulations that λ has an impact on the softening of spin waves, while no spin-wave softening can be observed in microcanonical SD. In actual physical systems, that is in systems where the spins are coupled to an underlying atomic lattice, both these phenomena (energy relaxation and spin-wave softening) are obviously caused by the very interaction between spins and lattice. Therefore, it seems natural to seek a connection between the relaxation coefficient λ of DLT and the physical parameters of the actual spin-lattice coupling.
To do that, we need a model of the interaction between spins and atomic lattice.

A. The simplest SD+MD dynamics
In the very simplest case the coupling between the spins and the underlying lattice can be described by the Hamiltonian [37,38], where r kl = |r k − r l |, while, is the deviation of nucleus k from its equilibrium position. The first term in H TOT is the spin exchange interaction, while the second and third terms describe the vibrational excitations of the lattice. At variance with microcanonical SD, the spins' positions are not fixed and the spin exchange interaction, J(r), depends on the distance. Within a microcanonical simulation of (34) (the SD+MD method described in the Introduction [13,14]), the total magnetization M is conserved, as well as the total energy H TOT of the system; however, there is an energy exchange between phonons and magnons, so that the spin magnetic energy is not conserved and it relaxes to its thermal value. Hence, we can interpret phonons as a conservative thermostat coupled to magnons, which is exactly the role of DLT.
To explore this analogy further, following [32,38,39] we expand the spin exchange coupling to first order, where, becomes the spin-spin coupling constant at the leading order, while, plays the role of the spin-lattice coupling constant; notice that α is positive, because the exchange interaction, J(r), decreases with distance. Finally, in (36) we have defined, We can now write the equations of motion for all the degrees of freedom, These equations conserve the total magnetisation thanks to the fact that the unit vectorê ik is antisymmetric. Compared to purely microcanonical SD, however, the spin equation, (40), has an extra spin-phonon term proportional to α, which is responsible for relaxing the magnetic energy.

B. Marginalizing the lattice degrees of freedom
In order to find an effective canonical dynamics for the sole spin variables, we follow the rather standard path [20] of marginalizing the spin-phonon term over the dynamics of the phonon variables, {u i , p i }, using thermalized initial conditions. To do this, one must first rewrite the dynamical equations using the eigenvectors of the discrete Laplacian, which is somewhat lengthy and cumbersome; for this reason, the full details of the calculation are reported in Appendix C.
Once this marginalization is performed, a relaxational term and a noise term emerge, linked to each other by the FD relation and both conserving the total magnetization, thus acting as an effective thermostat very similar to DLT. We can write the formal solution of equation (42) treating the spin term as an external driving force; then, the solution of (42) can be plugged into (40), and after some manipulations and reasonable approximations (described in Appendix C), we end up with, The first term at the r.h.s. is the usual reversible force coming from Poisson's brackets structure, where, is the standard Heisenberg Hamiltonian, as in (28). The second term in (43) is a relaxational force, which takes the form (see Appendix C), where R µν ij (t−t ′ ) is a dimensionless memory kernel, whose explicit form is specified in Appendix C (see equation (C36)) and that has the following crucial property, which guarantees the conservation of the total magnetization in (43). The third term, ξ i , is a noise, with variance, Using the same line of thought of Section II C, it is clear that this noise too conserves the total magnetization, because, thanks to (46) and (47), we have that i ξ i = 0 for each realization of the noise. All brackets, ⟨·⟩, in the relations above indicate an average over the initial conditions for phonons, as they are random variables drawn from the Gibbs-Boltzmann distribution. We also observe that the FD relation is satisfied, since in both the relaxational and the noise term the same kernel R µν ij (t−t ′ ) appears.

C. A simple dimensional argument
The derivation of equations (43)-(47) -derivation that we provide in great detail in Appendix C -is lengthy and cumbersome. Hence, let us give to the reader a very simple dimensional argument to calculate the effect of the spin-lattice coupling on the dynamics of the spins, and show how DLT quite naturally emerges in this context.
If we go back to equation (40), we see that the first term at the r.h.s. is the same as in pure SD dynamics, therefore it will not be affected by the marginalization over the phonon variables, while the second term, proportional to the spin-lattice coupling constant α, is the one that -containing the phonon degrees of freedom -will give rise both to the relaxational force, Ξ, and to the noise, ξ; let us concentrate on the latter, even though, of course, dimensionally the two are the same. By dimensional analysis of the spin-phonon term in (40) we get, where -again -brackets, ⟨·⟩, indicate an average over thermal phonons. By using equipartition of the phonons potential energy, we can rewrite (48) as, We immediately see that, dimensionally, this relation is the same as (47); moreover, the only reasonable way to generalize (49) to different sites, i, j, different coordinates, µ, ν, and different times, t, t ′ , without changing the dimensions, is indeed to introduce a dimensionless kernel R µν ij (t − t ′ ), which gives exactly equation (47); moreover, even if we did not do the explicit calculation of this kernel, we would still know that the conservation of magnetization that holds in the original equations (40) must survive marginalization over the phonons, hence we must have that i R µν ij = 0. Now that we have dimensionally justified the form of the effective noise variance (47), we can compare it with the DLT variance, that we report here for clarity, Such comparison may suggest that we are almost done: we have obtained a non-Markovian non-isotropic version of the DLT noise and one may be tempted to identify the DLT relaxation coefficient λ with (α 2 /ℏK); but that would be a dimensional mistake, because the kernel R µν ij (t − t ′ ) is dimensionless, while the combination Λ ij δ µν δ(t − t ′ ) has the dimensions of the inverse of a time.
To make progress, we introduce the memory time scale of the non-Markovian kernel, where we can forget about its possible ij and µν dependence as long as we are only proceeding through dimensional analysis. The advantage of having singled out τ m is that we can rewrite (47) as, where now the operator (R µν ij /τ m ) has the right dimensions to play the role of a 'fatter' δ-function, so that a direct comparison between (52) and the DLT original noise (50) finally gives, which is a quantitative relation between the relaxation coefficient of DLT and the microscopic parameters of the spin-lattice coupling. These results can actually be obtained analitically, and we refer the reader to Appendix C for the details. There, we show that -under some reasonable approximationsthe kernel is proportional to the discrete Laplacian, so that the effective noise after the phonons marginalization becomes indeed a non-Markovian version of the DLT noise, Then, the kosher way to proceed after this is to discretize the spin dynamics using a time interval ∆t ≫ τ m , so that the dynamics becomes effectively Markovian over time scales comparable to the discretization scale; in this way we finally obtain a Markovian noise that can be directly compared with the time-discrete version of the DLT noise (50), thus giving again equation (53) (see -as usual -Appendix C for the details). The relation we found between the DLT relaxation coefficient λ and the parameters of the spin-lattice systemequation (53) -looks very sound at the qualitative level: when K goes to infinity (i.e. for infinitely stiff -namely fixed -lattice) or when α goes to zero (no spin-lattice coupling), we obtain λ → 0, correctly recovering microcanonical SD. Hence, the relaxation coefficient of DLT wraps into a single quantity all the parameters of the (possibly very complicated) interaction between spins and lattice: λ is larger the stronger is the phonon-magnon coupling, and because this coupling is responsible for blunting spin waves [36], we finally understand why softening and damping within DLT are the stronger the larger is λ. But can we trust relation (53) at the quantitative level?
First, let us discuss what happens when λ is either too small or too large. As we have said, if there is no magnon-phonon interaction we obviously obtain λ = 0; but we know that 1/λ is the energy relaxation time, which does not diverge in real antiferromagnets. On the other hand, as we have seen from numerical simulations, too high a value of λ would completely suppress spin waves, which are a quintessential feature of real antiferromagnets. These considerations suggest that it may be possible to actually obtain a reasonable estimate of λ by plugging into (53) the microscopic parameters of actual magnetic materials found in the literature [13,16]. We can estimate the strength of the phonon-magnon coupling, α, from equation (38), in particular by using the characteristic amplitude and width of standard forms of the function J(r) (see -for example - Table A1 of [13]), getting α ∼ 1.0 − 2.0 × 10 −11 J/m. The lattice stiffness K can be derived from sound velocities, atomic masses, and lattice spacings for magnetic systems, which gives K ∼ 40 − 50 N/m. Estimating the value of the memory kernel time scale, τ m , is more challenging, as it depends on the phonon dynamics -itself coupled to the spins; it is known that the characteristic time scale of the phonon-phonon interaction is of the order of 1ps, while the spin-phonon interaction has a characteristic scale of 100ps (see [13]). It is, therefore, reasonable to expect that τ m should be within this range. Hence, we obtain that the value of the (dimensionless) relaxation coefficient in real materials should lie in the range 0.02 < λ < 10, which is consistent with the values of the relaxation coefficient used in our DLT simulation, where λ was chosen between 0.1 and 1 (see Figs. 2, 3 and 4).
In conclusion, the connection between the effective canonical dynamics provided by DLT and the actual microcanonical dynamics of systems with spin-lattice coupling, seems to be deeper than a generic numerical shortcut: given a certain material, with a certain set of parameters describing its microcanonical SD+MD dynamics, it should be possible to calculate the DLT relaxation coefficient λ through (53) and run an effective canonical dynamics appropriate for that specific material.

V. CONCLUSIONS
We have introduced a new thermostat that is at the same time conservative and stochastic, thus preserving the true dynamics of the spins, yet making it canonical. The potential applications of DLT seem promising. First of all, remaining at the equilibrium level, the method should be tested in other dynamical universality classes for which symmetries and conservation laws are important, as Model J (Heisenberg ferromagnet) and Models E/F (superfluid Helium) [1]. But DLT seems particularly promising for the study of out-of-equilibrium systems. As we anticipated in the Introduction, the method could be employed in the context of aging studies [11], as quenches and dynamic changes of temperature become straightforward within the canonical dynamics. Moreover, DLTor, more precisely, the conservative noise that we have introduce within DLT -can be used for the study of inherently out-of-equilibrium systems, such as those that violate detailed balance; in these cases, Metropolis-like Montecarlo simulations -whose dynamics is built upon the assumption of detailed balance -cannot be employed. One outstanding example in this class of systems is the conserved KPZ equation [19], which has non-thermal fluctuations. DLT noise could be an interesting new way to simulate such equations in real space.
Of course, DLT is not free of limitations. As with any other thermostat, and despite the connection between the relaxation coefficient and the spin-lattice parameters, a substantial part of the physical information regarding the microscopic degrees of freedom is lost, resulting in an effective simplified canonical dynamics. Moreover, we have seen how DLT can deal with respecting one conservation law, that of the total order parameter; it remains to be seen whether, and how, DLT could be adapted or generalized to the case of several conservation laws simultaneously holding in the dynamics.
An interesting issue for future studies is the possible emergence of dynamical crossovers related to the weak violation of the conservation law. In our -much -simplified microscopic model for spins and lattice we assumed that the only relevant interaction was the exchange interaction, which is isotropic, namely invariant under rotations in the internal space of spins, thus leading -through Noether's theorem -to the exact conservation law of the global magnetization. Within DLT, this exact conservation law is empowered by the fact that the relaxational force and the noise variance are proportional to the discrete Laplacian, Λ ij . However, anisotropic interactions, such as spin-orbit or dipole-dipole interaction, can in some cases be relevant, thus violating the symmetry and the conservation law. The interesting point is that we can easily generalize DLT to this cases, by adding a non-conservative component to the relaxation coefficient, namely, leaving all the other relations untouched. Now the effective friction, η, dissipates the total magnetization, hence it can be used to describe those physical cases in which the original conservation law is hindered by all sorts of non-symmetric interactions. When η is very large compared to the conservative coefficient λ, the symmetry and conservation laws are completely washed away, and one recovers the physics of overdamped non-conservative systems [15]. But when anisotropic interactions are weak, the value of η will be non-zero but small, so that the relative magnitude of the conservative relaxation coefficient, λ, vs the non-conservative one, η, may give rise to interesting finite-size dynamical crossovers. Such crossovers have been already studied in the field-theoretical context using the renormalization group [40], so that the possibility to analyze at the microscopic level the interplay between conserved and non-conserved relaxation seems a promising direction of investigation in the general field of magnetism. Finally, we notice that there might be a broader spectrum of applications of DLT beyond the realm of physical magnets. In biophysics there is growing experimental [41,42], numerical [43] and theoretical [44,45] evidence that conservation laws are important in determining the collective dynamical properties of some natural systems. More precisely, it has been noted that, although the mechanisms of biological imitation between neighbours in a biological group can be described by borrowing the standard concepts of ferromagnetism [46], an overdamped dynamics as in the standard non-conserved case (model A class [1]) does not reproduce some key experimental traits, as information propagation across the group and dynamic scaling [41,42]. On the other hand, it has emerged that the dynamics of underdamped mode-coupling systems, that is systems where a conservation law is important (as the Model G class [1] of the antiferromagnet studied here), is more suited to the study of collective behaviour and it has a better agreement with biological experiments [45]; essentially, the behavioural inertia of the biological individuals forces the interaction rules to comply to a second order dynamics that can only emerge as a result of a conservation law [41]. In the light of this, a conservative thermostat tailored on spin dynamics could become a relevant numerical tool also in physical biology.

The norm constraint
The DLT equation of motion reads, The first term (cross product) automatically conserves the norm of each spin. The second and third terms (relaxational and noise terms, respectively), however, would produce some deviations from the initially fixed norm if no countermeasure is taken. In order to fix the norm of the spins one could use a Lagrange multiplier; however, the equation to work out the Lagrange multiplier would need to be solved numerically at each time step, hence slowing down significantly the simulation. We therefore use a different method: we introduce a soft constraint on the norm by adding to the Hamiltonian the potential, so that the spin force becomes, It is easy to check that the reversible term is not affected by the soft constraint and that the total magnetization M remains constant. Hence, the soft constraint does not modify the reversible dynamics and it conserves the total magnetization, independently of the value chosen for A. For A → ∞, one would recover a strictly fixed norm, but of course for very large A we would need a much higher precision to integrate the equation. To set A, therefore, we just need a large value that preserves in practice the norm, but does not demand an extreme precision in the simulations. Fortunately, because nothing essential in the dynamics actually depends on A, finding this balance is not too difficult. We have fixed A = 100 and we have carefully checked that this value does not play a relevant role in the results presented throughout the main text and the SI. In Fig. 5 we can see that at this value of A, norm, energy, and order parameter are not far from their asymptotic, A = ∞, values.

Identification of the scaling regime
The critical temperature of the Heisenberg antiferromagnet in d = 3 is T c ∼ 1.446J [22,23]. However, the actual phase transition only occurs in the thermodynamic limit, where ξ ∼ (T − T c ) −ν → ∞ and the static susceptibility χ ∼ (T − T c ) −γ → ∞. In a finite system, the correlation length ξ cannot diverge, as it is limited by the the system's size, L. Finite-size scaling [48] states that near-criticality the susceptibility χ (reported in Fig.6a) obeys the relation, where F (y) is a dimensionless scaling function. This scaling form implies that if we plot χ/L γ/ν vs L(T − T c ) ν , we should have a collapse of the data, provided that we use the correct critical exponents of the Heisenberg antiferromagnet in d = 3 [25], that is γ/ν = 1.97 and ν = 0.70. This collapse is shown in Fig. 6b, which is already quite convincing. But to have a harder check, we need to work out the critical exponents in a more direct way than the collapse. The most effective way to explore the critical point at finite size is to keep fixed the scaling variable, y = L(T − T c ) ν , thus defining a size-dependent pseudo-critical temperature, T c (L); the easiest way to do that is to locate the point y c where F (y) has its maximum, which corresponds to locating the temperature T c (L) where χ(T, L) has its maximum, at each given L (see Fig.(6)a). In this scaling regime we have ξ(T c (L)) ∼ L, and because the scaling variable y c is kept constant by following the maximum of χ at each L. Equation (A5) is tested in Fig.2d in the main text; the best fit to a log-log representation of the data gives γ/ν = 1.92 ± 0.07, quite close to the correct value.

Relaxation time
Given the dynamic correlation function of the order parameter in the momentum-frequency domain, C(k, ω), we can obtain the static correlation function, C 0 (k) as, which gives the following normalization condition, A convenient and very robust definition of the characteristic frequency ω k is [31], The characteristic frequency, ω k , is naturally the inverse of the relaxation time, τ k , which -working out (A8) in the time domain -is defined by, which is the relation we use in our simulations to calculate the relaxation time (the transverse scattering function S T (k, ω) in the main text is the transverse part of C(k, ω)). The advantage of calculating τ k through equation (A9) is that no a priori fitting form for C(k, t) is needed; moreover, being an integral equation that uses the entire temporal range of C(k, t), this estimate is considerably more robust than crossing the correlation function with an arbitrary constant.

Tests of thermalization
We test that the analysis to obtain the exponent z has been done for a thermalized system. To do so, we check that the relaxation time τ does not depend on the duration of the simulation. We compute τ for trajectories of increasing duration t max (see Fig. (7)). The value of τ grows with t max until it reaches a stationary plateau. This procedure allows us to say that, when the duration is t max ≥ 10τ , we are sure that the system has thermalized; we repeat this test for all temperatures T and sizes L. In the main text, we report the dispersion relation for spin waves at low temperatures. In this section, we show how to derive this result for the case of the classical Heisenberg antiferromagnet in d = 3. Let us consider the microcanonical equation of motion, with ∂H ∂σ i = J j n ij σ j . From (B1), one obtains, Multiplying both sides by the parity of the site π i = ±1 (neighboring spins have opposite parity), we can rewrite Eq. (B2) in terms of the local staggered magnetizations, ψ i = π i σ i , In the first term of the r.h.s, the indexes j and k have the same parity, while the opposite holds for the second term. Eq. (B3) can then be rewritten as In the low temperature regime, i.e. T ≪ T c , the system is deeply polarized. We can thus expand the local staggered magnetization around the global polarization direction, that is, withn · π i = 0, |π i | ≪ 1, and where we have i π i = 0. The spin wave expansion is an expansion in terms of the {π i }. Because of the constraint, |ψ i | = 1, we have that ψ i ∥ = 1 − π 2 i . To first order, then, ψ i ∥ ≃ 1, and the equation of motion becomes, d 2 π i dt 2 =J 2 jk n ij [n ik (n × (n × π i ) +n × (π k ×n)) − n jk ((n × π j ) ×n + (π k ×n) ×n)] =J 2 jk [n ij n ik (−π i + π k ) − n ik n jk (−π j + π k )] .

(B7)
If we now define (for a cubic lattice) n c = k n ik = 2d, we get, where, we remind, Λ ij = −n ij + n c δ ij is the discrete Laplacian. In Fourier space, this equation can be easily expressed in terms of the eigenvalues of the Laplacian operator The resulting dispersion relation reads, and taking q = (q, 0, 0), one obtains, This expression is precisely the one reported in Eq. (30) and plotted in Fig.2f of the main text.
For a simple cubic lattice we have, with ℓ the lattice spacing. Due to the PBC, q can only assume discrete values, q = 2π L (n x , n y , n z ), with L the lattice size, n x , n y , n z integers from 0 to n max = (L/ℓ − 1). The change of basis using the eigenvectors of the discrete Laplacian gives, Hence we have, To rewrite the interaction energy we define the auxiliary quantities, and so we get, hence, Finally, combining all terms, we end up with, We will now write the equations of motion for both the spins and the lattice vibrational degrees of freedom. In the new basis, the variables {Q q } are not directly coupled one with the other; we can therefore solve explicitly their equations and plug back the solutions in the equations for the spins. The result is an equation of motion for the spin variables only, where the effect of the phonons is the appearance of a thermal bath composed by a "noise" and a "relaxational" term, both of which preserve the conservation of the global magnetization. This thermal bath always satisfies the Fluctuation Dissipation (FD) relation, but the specific memory kernel will depend on the details of the dynamics. In the following sections, we are going to derive it for both an isolated system, and in the case where the nuclei, but not the spins, are in contact with a (standard) thermal bath.

Isolated system a. Equations of motion
The outline for deriving a thermal bath out of many oscillatory degrees of freedom of an isolated system can be found in [20]; in the case of ferromagnetic systems, a similar computation has been done in [49] within a mesoscopic field theory framework. We start by writing the equations of motion for phonons and spins, The plan is now to solve Eq.(C20) and plug the result back into Eq.(C21). The solution of the homogeneous equation is given by, where Q q (0) and P q (0) are the initial conditions at time t = t 0 = 0, and, The complete solution is therefore given by, Substituting into the equations of motion for the spins we get, where the effective Hamiltonian is given by, the relaxational term is given by (square brackets indicate functional dependence), and the noise term is given by, Equation (C21) is the same equation of motion (written in Fourier's space) as Eq. (40) of the main text, used to find the dependence of λ on the microscopic parameters. We notice, however, that the spin-lattice interaction term in the r.h.s. of Eq. (C21), which depends on Q q , is generating not only the noise (as in the simplified argument of the main text), but also the relaxational term. This is what we should expect, since they always emerge together in actual physical systems.

b. Noise
In the context of the microcanonical dynamics discussed in this section, if we want to describe the properties of the system at a given temperature T , the natural choice is to consider the initial conditions as drawn with a canonical distribution e −βH /Z. We note that the total Hamiltonian of Eq. (C19) can be rewritten as, We therefore have for the canonical averages, We can exploit these expressions to compute the properties of the effective noise ξ i . From Eq.(C28)) we get, and, from Eqs. (C31) and (C33), Hence, the {ξ i } are random variables with a Gaussian distribution, zero mean and non-trivial correlations both in space, as noises in different sites i and j are not independent, and in time. We can indeed identify in Eq. (C35) the following non-Markovian, dimensionless memory kernel, The noise correlator can then be expressed in a more compact form as, Equation (C37) is the same eq. as (47) of the main text, and already looks quite similar to (55), written in the main text using dimensional analysis. Some more work, however, is needed to recover the same thermostat as DLT.
Using the above equation we find, We then find that the relaxational term can be written as, Hence, the same memory kernel appears in the relaxational term and in the noise correlator, Eq.(C37). Together with the factor k B T , this fully recovers the FD relation. We recall here that the FD relation must be satisfied in order to have a thermal bath, which preserves the equilibrium probability distribution of the system.

d. Conservation of the total Magnetization
Let us prove here that the dimensionless kernel appearing in both the noise and relaxational terms preserves the total magnetization along the dynamics. Let us consider Eq. (C25), the sum over i yields the derivative of the total magnetization of the system, We now show that each term of the r.h.s. is zero. Recalling the definition of f β σ in Eq (C16), we note that, From this expression we get, where the last passage derives from summing an anti-symmetric object. Thanks to this relation, and recalling (C36), we conclude that, which implies that the second and third term in (C41) are zero. For what concerns the first term in (C41), we notice that it can be rewritten as, where the first term in the r.h.s. gives 0, as we knew from the very beginning, while the second vanishes due to Eq. (C43). The total magnetization is therefore strictly conserved during the dynamical evolution.

Approximations for recovering DLT
Up to this point, no approximation has been done, and both the conservation of the total magnetization and the FD relation are satisfied exactly. The equation of motion for the spin that we have obtained is Eq. (C46) looks similar in structure to the equation of motion introduced in the main text in the DLT scheme: there is a Hamiltonian precession term, a conservative noise term, and a conservative relaxational term linked to the derivative of the Hamiltonian. However, there are several differences. First of all, the derivative of the Hamiltonian in the relaxational term involves the total Hamiltonian, which includes the phononic degrees of freedom making the equation itself not closed in the spin variables. The precession term involves the effective Hamiltonian (C26), instead of a simple spin exchange interaction. The relaxation term is proportional to something which, while having a zero mode, is not the discrete Laplacian. Finally, the memory kernel is not proportional to a Dirac delta distribution and the equation is therefore not Markovian. We address here the first three issues, while leaving the fourth to a separate section at the end of the SI.

a. Approximation for H TOT in the relaxational term
Eq.(C46) involves the derivative of the total Hamiltonian with respect to the spin variables σ j . In the simpler case of a colloidal particle coupled to a bath of harmonic oscillators (see [20]), this term does not depend on the degrees of freedom of the thermal bath, leaving a closed equation for the colloidal particle. In the case of spin systems, though, this derivative includes the phononic degrees of freedom, so we have to make an approximation to close the equation. From the expression of the total Hamiltonian (C19) and the definition of H eff we get The phononic variables Q β −p (t) on the r.h.s. are random variables, since they depend on the initial condition Q β −p (0), which we have drawn from the canonical distribution. If we assume that their dynamics is faster than the one of the spin variables we can approximate their time dependent value with the expected value with respect to their marginal equilibrium distribution at fixed σ. Hence, the second term of Eq. (C47) can be approximated with its average. Since the initial canonical distribution is invariant under the microcanonical dynamics, this average value is zero (see Eq. (C30)), thus giving

The effective Hamiltonian approximation
The effective Hamiltonian H eff appears now in both the precession and the relaxational term, and the equation of motion becomes, Eq. (C49) involves a Hamiltonian, which is a function of the spins only, and it is therefore self-consistent. In the main text, though, we used an equation for the spin dynamics where only H 0 appears, in place of H eff . In general, corrections to H 0 might be potentially relevant for the equilibrium probability distribution of the slower variables. However, this is not likely what happens in our case. Indeed, if we consider the difference between H 0 and H eff we have that, and using, Eq. (C16), The above expression is quite complicated, but we can see that both the ground state of the ferromagnet (all the spins aligned) and the ground state of the antiferromagnet (all spins anti-aligned with their neighbors) yield ∆H = 0 . Hence, we argue that the correction does not change any fundamental feature of the system and it is safe to neglect this contribution. This is also consistent with the fact that in the starting Hamiltonian we neglected all the 4-spins (or more) contributions keeping only 2-spins interactions. With this approximation, we hence have,

. Introducing the Laplacian
We now focus on the site dependence of the memory kernel. As discussed in the previous sections, the kernel exactly has a zero mode, which ensures the conservation law for the global magnetization. Adopting some additional approximations, we can show that it is proportional to the discrete Laplacian, with the appropriate q-dependent coefficient. Let us call, the site-dependent part of the memory kernel. Using Eq. (C43) we find, It is also convenient to introduce the following auxiliary quantity, First, we assume isotropy with respect to the spin vectorial components, so that this matrix becomes proportional to the identity in that space, By applying the vectorial identity we can rewrite this quantity in a more useful form, We then approximate the above expression with its expected value with respect to the stationary spin distribution, and we consider spin correlations to be negligible unless the scalar product is taken between two spins on the same site (i.e. for i = j and h = k in the first term and h = j and i = k in the second term), A four-point correlation function of the spins appears in Eq. (C59), namely, where no site dependence is present in the l.h.s because i and h are always first neighbours (see (C54)) and the resulting expected value is therefore the same for all the first neighbours pairs. The time dependences appear only as a time difference, because expected values are taken with respect to the stationary distribution. The actual shape of this function is not important: we only assume that the time dependence is much slower than the phononic dynamics, so that we can assume it to be constant in time and approximate it with its value in 0, namely C (4) (t − t ′ ) ≃ C (4) (0) = 1. Now we can use the above approximations in formula (C54) and obtain, We remind thatê hi = −ê ih , thus finding, The indexes i and j identify a well defined direction α = α(i, j), so the last equation reads, Lastly, we make the assumption of isotropy in the q-space, hence we end up with, With this expression, we finally recover the Laplacian matrix appearing in the DLT scheme. We now can reduce the memory kernel introduced before to, where, And we can rewrite, after all the approximations discussed so far, Eq.(C49) as, with the noise correlation, This dynamics is now nearly identical to DLT. The only remaining difference is that the dimensionless memory kernel R(t − t ′ ) is not Markovian, as we hinted in our shortened derivation in the main text using dimensional analysis. The time dependence of the kernel is in general non-trivial and it depends on the details of the phononic spectrum, and on the phonon dynamics. So far, we have considered the case where the system is isolated and the dynamics is fully microcanonic. In the next section, we compute the time dependence of the memory kernel when the phononic degrees of freedom (but not the spins) are coupled with an external thermal bath.

Phonons coupled to an external thermal bath a. Equations of motion
In this Section we perform the same computation as in Sec. C 2 a, but considering this time the vibrational degrees of freedom coupled to an external thermal bath. The spin variables, on the other hand, still follow Hamiltonian equations of motion and feel the presence of the thermal bath only via their interaction with the lattice. We use for the phonons a standard Gaussian white thermal bath, with friction coefficient η and temperature T . Instead of Eq. (C20) we therefore have, where the bath noise satisfies, There are three relevant (inverse) timescales in this equation, i.e.
The solution of the homogeneous equation is not relevant, since it only gives an exponentially decaying transient, and for large enough times that system forgets the initial conditions. Setting the initial condition at t 0 = −∞ the solution is therefore given by, The above equations can then be inserted back in the equations of motion for the σ variables, Proceeding as in Sec. C 2 a we find, where H eff is the same defined in Eq. (C26), but the relaxational and the noise terms now have a different expression, (C77) and, Let us start by considering the noise term ξ i defined in (C78). In the microcanonical treatment of Sec. C 2 a, the noise was a fluctuating variable due to the presence of the initial conditions, which were drawn randomly with a canonical distribution. Here, on the other hand, the initial conditions are not relevant and the fluctuating nature of this term is due to the presence of the stochastic variable ζ β q representing the external thermal bath coupled to phonons. Taking the average of Eq.(C78) over the thermal bath at fixed σ we get, and, for the noise correlator, approximations, to the following effective equations of motion for the spins only, 1 where ξ i is a Gaussian noise, with zero mean and correlation given by, R(t − t ′ ) is a dimensionless non-Markovian memory kernel dependent on the details of the phonons. For an isolated system, we have, while for a system where phonons (but not spins) are in contact with a thermal bath, we have,

Discrete-time Markovian approximation
Because of the memory kernel in Eq.(C87), what we have obtained so far is a non-Markovian thermostat, while the DLT presented in the main text is Markovian. In this section we show that, under appropriate conditions, we can recover a Markovian dynamics by discretizing the continuous stochastic process over a sufficiently long time scale. Every stochastic differential equation must ultimately be brought back to its discrete time version for it to make sense, so we rewrite (C87) as, where dw i = ξ i dt is the noise on scale dt. The microscopic timescale dt has to be interpreted as the minimal physical time increment over which the system's dynamics takes place. The derivation of the previous sections shows that on this minimal timescale the noise correlator has memory, i.e. ⟨dw i (t) · dw i (t ′ )⟩ ∼ R(t − t ′ )dt 2 . We can however ask whether, when considering the dynamics on scales larger than this minimal increment, a Markovian process can be recovered. Let us assume that the spins vary in time slowly compared with the time span of R(t − t ′ ). In this case, we immediately see that a simplification occurs in Eq. (C91), because the quantity ∂H 0 /∂σ j only depends on the spin variables and can be brought out of the integral. If we then define, Eq. (C91) can be rewritten as, Consistently with the above assumption, it is reasonable to look at the spin dynamics over a lower time resolution, because changes in the spins are negligible on the minimal scale dt. We can then consider a discrete time increment ∆t ≫ τ m , which is much larger than the decay time of R, but also much smaller than the time scale over which the spins vary significantly. The dynamical evolution of the spins on this coarse-grained timescale can be easily obtained by integrating Eq. (C93) between t and t + ∆t, and considering the spins to be constant in this integration interval. We get, where, Eq. (C94) is now the lower time resolution version of Eq. (C91). Memory is not present anymore in the relaxational term and ∆w i is the coarse-grained noise over scale ∆t. From (C95) we find, where t 1 and t 2 are multiples of the new time resolution ∆t. Now, since ∆t is much larger than the decay time of R, if t 1 ̸ = t 2 Eq. (C96) gives 0; on the other hand, if t 1 = t 2 , all the times where the function is significantly different from 0 are inside the integration range, thus yielding the aforementioned τ m . Therefore, we have We then recover white uncorretaled noise (i.e. ∆w µ i (t) is a Wiener process on scale ∆t), together with the correct FD relation between its correlation and the relaxational coefficient.
Equation (C94) is therefore Markovian and it has the very same structure on scale ∆t as the DLT equation used in the main text. This can be seen explicitly by discretizing Eq. (11) of the main text on scale ∆t, which becomes with ⟨∆w µ i (t 1 )∆w ν j (t 2 )⟩ = 2k B T ℏ −1 λΛ ij δ µν ∆t δ t1,t2 .
By comparing these two equations with Eqs. (C94) (C97) we can finally identify the relaxational coefficient λ as, which is (finally) equation (53) of the main text.