Quantum Brownian Motion for Magnets

Spin precession in magnetic materials is commonly modelled with the classical phenomenological Landau-Lifshitz-Gilbert (LLG) equation. Based on a quantized spin+environment Hamiltonian, we here derive a general spin operator equation of motion that describes three-dimensional precession and damping and consistently accounts for effects arising from memory, coloured noise and quantum statistics. The LLG equation is recovered as its classical, Ohmic approximation. We further introduce resonant Lorentzian system--reservoir couplings that allow a systematic comparison of dynamics between Ohmic and non--Ohmic regimes. Finally, we simulate the full non-Markovian dynamics of a spin in the semi--classical limit. At low temperatures, our numerical results demonstrate a characteristic reduction and flattening of the steady state spin alignment with an external field, caused by the quantum statistics of the environment. The results provide a powerful framework to explore general three-dimensional dissipation in quantum thermodynamics.

The continued miniaturisation of critical components in consumer electronics and neighbouring technologies will require a deeper understanding of thermal noise and general thermodynamic principles beyond the classical macroscopic world. Quantum thermodynamics [1][2][3] has emerged as a field addressing the conceptual challenges related to the exchange of energy and information at the nanoscale. Recent advances include studies of heat transport in quantum systems [4][5][6][7][8][9][10][11], the characterisation of memory effects in their dynamics [12][13][14][15][16][17][18][19], and clarification of the impact of quantum coherence and correlation on thermodynamic processes [20][21][22][23][24]. The establishment of a generalised thermodynamic framework, valid for nanoscale systems that strongly couple to environmental modes, is well under way [25][26][27][28][29][30][31][32][33], and for magnetic molecules an environment-induced renormalisation of the anisotropy has been predicted [34]. Two open quantum systems models have served as the workhorse for many of these conceptual studies; the Caldeira-Leggett model for quantum Brownian motion [8,35,36] and the spinboson model of a spin (or many spins) coupled to a one-dimensional harmonic bath [5,6,33,[37][38][39]. These describe a very wide range of physical situations extending to studies of quantum effects in bio-chemical reactions [40], where they are used to model exciton-phonon interactions [41].
Until now few nanoscale technologies have required the use of advanced open quantum systems techniques. But advances in engineering magnetic materials for magnetic hard drives at unprecedented length and timescales [42] are likely to require a more detailed picture of spin dynamics including memory and quantum signatures. Here we introduce a three-dimensional open quantum system model to characterise the quantum Brownian motion of spins in magnetic materials.
Magnetic behaviour has been studied extensively based on the classical phenomenological Landau-Lifshitz-Gilbert (LLG) equation [43][44][45][46][47] ∂M ∂t = γM × B eff − η G ∂M ∂t , (1) * janet@qipc.org which is routinely solved with micromagnetic and atomistic simulations. Here M is the magnetic moment, γ is the gyromagnetic ratio and B eff is an effective magnetic field which includes the external field [106], exchange and anisotropy effects, as well as stochastic magnetic noise b ∝ √ T stemming from an environment at temperature T that was added by Brown [107] in 1963 [48]. The final term on the right of (1) is the so-called "Gilbert damping" term and the positive constant η G is the damping parameter [108], which is often rewritten as η G = η/|M ||γ| with a unit-free η.
Gilbert damping is not derived from microscopic principles, but chosen as the simplest term that could serve to align the magnetic moment with the applied field [43]. As we will see, it contains no memory which is increasingly seen as a limitation [49,50]. Advances in engineering magnetic materials at the nanoscale and manipulating them on ultrafast timescales indicate that a theory beyond the classical LLG equation is required [47]. Early attempts have pursued a path integral derivation of a quantum spin dynamics equation [51], as well as other conceptually related classical and quantum derivations [52,53]. These derivations were not directly applied to the calculation of magnetization dynamics or steady states, nor have they been connected to recent generalizations of Gilbert damping that include inertial terms [49,50,[54][55][56][57][58][59] or provided an assessment of quantum effects.
Here we go further and develop a comprehensive and quantum-thermodynamically consistent theory suitable to describe the quantum dynamics of spins in magnetic materials including non-Markovian damping, coloured noise and quantum zero-point fluctuations. Unlike the conceptually pioneering Caldeira-Leggett model that has few experimental realisations, the developed threedimensional quantum spin model is directly applicable for atomistic spin dynamics simulations [47,60], ultrafast magnetism experiments [61], and systems exhibiting anisotropic damping [62].
The paper is organised as follows: In section I the general quantum spin dynamics equation for spin operator precession in three dimensions is derived. For the simplest, Ohmic, coupling this equation is found to reduce to the memory-free LLG equation. In section II we introduce Lorentzian couplings as a systematic method for exploring non-Markovian dynamical regimes in gen-eral open quantum systems, including spins. Finally, in section III we detail a numerical method to simulate non-Markovian dynamics, and present results for a single classical spin that illustrate the differences between spin dynamics and steady states arising with non-trivial memory, coloured noise, and quantum bath statistics in comparison to those obtained with the memory-free LLG equation. Conclusions and open questions are discussed in section IV.

I. Quantum spin dynamics equation A. System+environment Hamiltonian
We begin by introducing the quantized Hamiltonian describing the different contributions to the total energy of the system, consisting of spins as well as environmental degrees of freedom (e.g. electrons and phonons), given byĤ whereĤ S is the bare spin Hamiltonian operator which captures the spin energy in external fields and interactions between spins,Ĥ R is the environmental or reservoir Hamiltonian, andV int is the interaction between the spins and the reservoir. We chooseĤ S as the sum of the interaction with a homogeneous external field B ext [109] and the exchange interaction between three-dimensional spin vector operatorsŜ 3 ) at sites n of a lattice [110], Here J (nm) is the exchange tensor for spin pairs (n, m) [111], which can include the Dzyaloshinskii-Moriya interaction [47]. It is straightforward to include additional energetic terms in the bare spin Hamiltonian, such as energies associated with magnetic anisotropy. Instead of the magnetic moment M used in Eq. (1), we will here work with the spin angular momentum S proportional to M , M = γS, where γ is the gyromagnetic ratio. In the following we will assume the gyromagnetic ratio γ = −g e µ B / = −1.76 · 10 11 s −1 T −1 for an electron. The reservoir Hamiltonian is commonly modelled as a set of harmonic oscillators [35,64], and we here follow the continuous reservoir approach by Huttner and Barnett [64], taking the reservoir Hamiltonian aŝ It describes a continuous frequency reservoir at each lattice site n, whereΠ ω are (three-dimensional) momentum and position operators of the reservoir oscillator with frequency ω. The position operatorsX (n) ω physically represent variations in the environment to which the spin at site n responds, see illustration Fig. 1, as for example, in magnon-phonon mediated loss [65]. Unlike most system+environment Hamiltonians which assume one-dimensional coupling, we here take the spinreservoir interaction to be of the three-dimensional form This coupling allows angular momentum transfer as well as energy transfer between the spins and the environment. Here C (n) ω is a three-dimensional coupling tensor and a function of frequency ω. At each ω, the coupling tensor determines the strength of the coupling of each spin to its reservoir oscillators at frequency ω, thus acting as a frequency filter. As we shall see, the choice of the coupling C (n) ω will determine the damping of the spin dynamics as well as the stochastic noise experienced by the spins.
For readers concerned about time-reversal symmetry ofV int in (5), we note thatX (n) ω should be interpreted as an effective magnetic field seen by the spins due to their interaction with the environment, which has the same time symmetry asŜ (n) [112]. Indeed, the theory of magnetic materials developed here on the basis of the system+environment HamiltonianĤ is analogous to macroscopic QED, an effective medium theory which successfully describes quantum electromagnetism in dielectric materials [64,66,67]. Instead of trying to give a fully microscopic description that accounts for every light-matter interaction in the material, macroscopic QED characterizes electromagnetic materials in terms of two measured frequency dependent susceptibilities. The quantum Hamiltonian is then written in terms of these susceptibilities, and can be used in applications from predicting the Lamb shift to the Casimir effect [66,68].

B. Equations of motion
Having set up the full Hamiltonian (2) of the spins and environment degrees of freedom allows the study of the spins' reduced state dynamicsρ S (t) = tr R [ρ SR (t)] of the total statê where the reservoir stateρ R = e −βĤ R /tr[e −βĤ R ] is thermal at some inverse temperature β = 1/k B T . In what follows it will be more convenient to work instead in the Heisenberg picture where the state is stationary,ρ SR (0), while the time dependence of an opera-torÔ(t) is governed by the commutator dÔ(t)/dt = (i/ )[Ĥ,Ô(t)]. Expectation values at time t can then be obtained as FIG. 1. Illustration of Hamiltonian model. In addition to precessing in an external field Bext, and coupling to its spin neighbours m with strength J (nm) , each spinŜ (n) couples to its environmental mode (phonons, electrons) at frequency ω with a coupling function C (n) ω . All environmental modes are assumed to be thermal at the same temperature T . ω ,k ] = i δ nm δ jk δ(ω − ω ), we obtain the following equations of motion for the spin operatorsŜ is a magnetic field operator generated by the reservoir oscillator positionsX (n) ω at site n. In turn, the equations of motion for these operators are i.e. the reservoir oscillators are driven by the motion of the spins, with the (transposed) coupling tensors C (n) T ω governing the degree of driving for each of the continuum of oscillators. We assume retarded boundary conditions so that the reservoir responds only to the past behaviour of the spins. The retarded Green , and Eq. (9) can then be solved exactly bŷ Here,â (n) ω andâ (n) † ω are (vectors of) bosonic ladder operators with their components obeying [â . Classically these correspond to the two integration constants for the differential equation (9) which set the initial amplitude and velocity of the oscillator.
Substituting the reservoir solutions (10) into the equations of motion for the spins (8), we obtain the first result: The Heisenberg-Langevin equation that governs three-dimensional quantum spin dynamics under the influence of memory and coloured quantum noise is The termb (n) (t) is a Hermitian magnetic noise operator for site n, which plays the role of the stochastic noise first described by Brown [48]. Here it arises from the spin's interaction with its reservoir. As we will see below, the bath noise can be coloured and contain quantum zeropoint fluctuations. In addition to the coloured noisê b (n) , a kernel tensor K (n) (t − t ) appears in Eq. (11), which captures the damping of the spins. It arises from the coupling tensor C ω and is given by where τ = t − t . Here Θ is the Heaviside function which makes the spin's dynamics at time t, see (11), a function of the spin's state at previous times t 0 ≤ t < t [113].
The three-dimensional spin dynamics equation (11) describes the evolution of spin operators and explicitly includes memory of the past dynamics (non-Markovianity). This contrasts with previous derivations of quantum spin dynamics in the form of a master equation [69] which can be solved numerically. But to obtain the master equation a range of simplifying assumptions where made, including as weak spin-environment coupling, as well as the Markov and secular approximations. These approximations are quite strong and may not always be justified for a given spin system. We further remark that a different method of including (classical) coloured noise is based on the Miyazaki-Seki approach [70]. Similar to (8), this model assumes that the equation of motion of the spins is coupled to a stochastic equation of motion for the lattice enabling transfer of energy and angular momentum between the lattice and the spin systems. Different coupling potentials, such as harmonic and Morse potentials, have been considered and the spin-lattice coupling impact on the magnetisation has been characterised [71].
The above spin+environment Hamiltonian and its resulting equations of motion share many similarities with the well-known Caldeira-Leggett model for harmonic quantum Brownian motion [26,35] and the spinboson model [33,37,39]. A key difference is the threedimensional nature of the reservoir interaction in (5) which leads to the cross product in (11). The spin-boson model is recovered as a special case for spin 1/2 operators and rank-1 coupling tensors, see Appendix A1.

C. Fluctuation-Dissipation Relation andŜ 2
The presence of the kernel in (11) gives rise to spin damping that can significantly differ from Gilbert damping, as explored further in section II. Previous generalisations to Gilbert damping have considered the spins' interactions with the lattice and electron motion [54,[71][72][73], and included inertial terms [49] and further memory terms within a damping kernel [55][56][57][58]. First measurements have recently confirmed the presence of inertial corrections [50,59]. But these studies have not considered how the noise may have to change from the standard white magnetic noise that is included in the effective magnetic field in the LLG equation (1).
A separate strand of reasoning has argued that to extend the LLG equation to the quantum setting requires that the stochastic noise included in the effective field B eff should have a quantum distribution rather than the classical Boltzmann one [60,[74][75][76][77]. These studies partially reproduce experimentally measured magnetization over temperature curves. But the noise considered vanishes at low temperatures, whereas the full Bose-Einstein distribution famously has non-zero quantum fluctuations even at T = 0 K. Other pioneering investigations have discussed the effect of coloured noise and non-Markovian effects on the ultrafast demagnetization rate [70], while not requiring the fluctuation-dissipation theorem (FDT) to hold and not including quantum effects.
The quantum spin dynamics equation (11) brings these pieces, the damping kernel, quantum statistics and coloured noise, together in a quantum thermodynamically consistent framework. To see this we now prove that the reservoir's two signatures -the stochastic field b (n) and memory kernel K (n) -always fulfil the (quantum) fluctuation-dissipation theorem, for any choice of coupling tensor C where {. , .} is the anti-commutator. Inserting (12), one finds that the only non-trivial contributions to this expression come from the bosonic ladder operator expectation values {â , is related to the imaginary part of the Fourier transform of the damping kernel K (n) (τ ), i.e.
Hence one obtains for all choices of the interaction tensor C (n) ω in Eq. (5) the quantum fluctuation-dissipation theorem (FDT) where the last line is the well-known classical hightemperature (or low frequency) approximation, and the first line is the low-temperature limit where the Bose-Einstein distribution of the reservoir oscillator modes becomes important [116].
Finally, the spin dynamics equation (11) leaves the square of the spin operators a constant of motion, since see Appendix A3. HereB is an effective magnetic field at site n see Eq. (8), which commutes with all components ofŜ (n) . This confirms that the spin length |Ŝ is a constant in time, e.g. for a spin-1/2Ŝ = 2σ withσ j the Pauli matrices, one has a constantŜ

II. Coupling functions
The general quantum spin dynamics equation (11) is specified completely by the coupling tensor C ω in (5) that sets the interaction between spins and reservoirs. Here we will show that the standard LLG equation (1) arises as a special case of (11), for a particular choice of coupling C ω . We will further introduce a class of Lorentzian coupling functions which allow a systematic exploration of spin dynamics behaviours beyond the LLG equation.
For simplicity from here on we will drop the lattice site superscripts and consider isotropic coupling tensors C ω = C ω 1 3 with C ω a scalar function, and similarly for the corresponding kernels K(τ ) = K ω (τ ) 1 3 and power spectraP(ω) =P (ω)1 3 . These choices are appropriate for 3D materials in which all spins couple to the environment in the same manner in all spatial directions. For other materials, such as 2D layers, non-isotropic coupling tensors can be considered.

A. Ohmic coupling and recovery of LLG equation
In the open quantum systems literature, coupling functions that are linear in frequency are referred to as "Ohmic", while those proportional to higher and lower powers of ω are called super-and sub-Ohmic, respectively. For the magnetic system considered here, Ohmic coupling means where η G is a positive constant with units of kg/A 2 m 2 s, as in Eq. (1), and + is an infinitesimal positive constant which we take to zero at the end of the calculation. The corresponding kernel (13) is close to instantaneous, from which we can see that a positive value of + ensures the response kernel K Ohm is causal. For Ohmic coupling (18), the damping term in (11) reduces to i.e. Ohmic coupling results in a damping term proportional to a first order time derivative of the system variable. In the magnetic case this isŜ, while for quantum Brownian motion it would be the oscillator position [35]. Inserting (20) into the general quantum spin dynamics equation (11), we see that it recovers (the quantum version of) the LLG equation (1) with η G being identified as the damping parameter. This becomes the well-known classical equation in the limit of large spins. Consequently we will also refer to the Ohmic coupling function (18) as "LLG coupling" [117]. Given a coupling function it is straightforward to find the corresponding power spectrum (16a) governing the time correlations of the quantum noise. For C Ohm ω the quantum power spectrum is where we have taken the limit + → 0. The Ohmic quantum power spectrum (21) is proportional to the damping parameter η G , and tends to a linear (diverging) function of frequency in the low temperature limit ω k B T , see Fig. 2g+h. However, power spectra diverging with frequency are unphysical. Coupling any system to environmental modes must go to zero at large enough frequencies, such as the interaction of spins with lattice phonons or the effect of conduction electrons scattering off the spins. A common way of integrating this physical information into Ohmic coupling is to introduce a cut-off, i.e. an upper frequency to which the coupling grows linearly, and after which it is set to 0 or decays to 0 in either algebraic or exponential form [38]. A different approach was taken in [60,74] for the modelling of spin dynamics, where a semi-quantum Ohmic power spectrum similar to (21) was chosen, but with the zero-point noise subtracted ensuring that at T = 0K it vanishes for all frequencies.
In the high temperature limit k B T ω, Eq. (21) reduces to the classical power spectrum which is frequency independent (white) noise. Thus Ohmic coupling plus the high temperature approximation to the power spectrum recovers the LLG equation (1) commonly used to simulate magnetic materials over a wide range of temperatures.
As a first step to unravel how the dynamics predicted by the quantum spin equation (11) can deviate from that predicted by the LLG equation (1), one can expand the spin vector (operator) at time t in Eq. (11) around the end point of integration to arbitrary order, The 0-th order coefficient, κ 0 does not appear in the expansion of the damping operator (23), asŜ(t) ×Ŝ(t) = 0. For the quantum operatorsŜ(t) this is ensured by the angular momentum commutation relations, which make this cross product anti-Hermitian. As outlined in Appendix A2, when using any truncated form of expansion (23) one must employ an explicitly Hermitian form of Eq. (11), and hence any anti-Hermitian contributions drop out. The higher expansion coefficients κ m>0 are proportional to the m th one-sided moment of the memory kernel where we have assumed that the dynamics has been running for some time longer than the kernel decay time, for which one can replace the initial time t 0 by −∞.
For the Ohmic kernel only the first moment is nonzero and corresponds to the (negative) damping parameter, while κ Ohm m>1 = 0. This complete lack of higher moments, which would maintain a certain degree of memory in the dynamics, shows that Ohmic coupling dynamics can only be an approximation to any real dynamics. For example within magnetism, the memory-free form of damping (20) is known as Gilbert damping [43], and is almost universally used to describe magnetization dynamics through the LLG equation (1). However, "inertial" corrections to such dynamics have been proposed [49] and their presence was recently confirmed experimentally [50].

B. Lorentzian coupling
Here we provide a tool to systematically study dynamics beyond the Ohmic case, allowing one to include memory and coloured noise effects in a manner consistent with the quantum fluctuation dissipation theorem (16a). We consider the class of Lorentzian coupling functions where A is a coupling amplitude, with the following properties: i) for small ω, C Lor ω grows linearly with ω and can be approximated by an Ohmic coupling function, ii) at large ω, C Lor ω smoothly decays to zero, and iii) at some intermediate frequency ω 0 , C Lor ω has a resonant peak with some width Γ. This peak characterises the confined range and relative strength of systemenvironment interaction with two parameters. Alternative "peaks" such as Gaussians or top hat functions could be considered, but here we chose the Lorentzian shape due to the fact that many expressions can be solved analytically and, as we will demonstrate in section III, Lorentzian couplings allow us to efficiently simulate non-Markovian dynamics.
We call the above functions "Lorentzian coupling" since the corresponding damping kernel in the frequency domain is the widely studied Lorentzian response where the imaginary part is obtained from (15) and the real part is determined using the usual Kramers-Kronig relations. In the time-domain the Lorentzian memory kernel is where ω 1 = ω 2 0 − Γ 2 4 and Γ/2 can now be interpreted as the kernel decay rate. For the coupling function (26), the collective response of the environment is thus equivalent to a single harmonic oscillator of resonant frequency ω 0 and damping rate Γ/2 [78]. From the quantum FDT (16a) it follows that the corresponding power spectrum is which takes its largest values at frequencies close to ω 0 and tends to zero as ω −3 at large ω. This power spectrum differs from classical Ohmic noise (22) in two important respects. Firstly, the quantum mechanical treatment means that the low temperature noise is not proportional to temperature, and does not vanish at zero temperature. Secondly, even in the high temperature limit, the noise spectrum is frequency dependent (coloured), unlike the white noise of Eq. (22). The presented theory thus captures both of these aspects in a consistent quantum thermodynamic framework.
Unlike Ohmic coupling (18) which depends on a single parameter η G , the Lorentzian coupling function (26), and hence its kernel and spectrum, depends on three parameters. These allow a systematic study of different regimes of the environment response. Specifically, the memory time of the environment can be continuously varied by changing ω 0 and Γ, which can lead to very different magnetic behaviour. Beyond the spin dynamics explored here, the proposed Lorentzian coupling may also be a useful tool for the characterisation of quantum Brownian motion of a variety of systems, including oscillators and free particles.

C. Two coupling regimes
To better understand the relation between the Ohmic and Lorentzian coupling functions, one may consider their kernel moments κ m in expansion (23). In contrast to the Ohmic case, for the Lorentzian kernel (28) all κ m are non-zero and given by The first relevant two moments are and when comparing to the Ohmic case, one finds that the first Lorentzian moment can be identified with minus the damping parameter η G , i.e. κ Lor For a material with a given damping parameter η G this fixes one of the Lorentzian parameters, i.e.
which now only depends on the two parameters ω 0 and Γ. For a specific material these may be approximately determined through information contained in the density of states of the environment to which the spins couple [79].
Inserting expansion (23) with Lorentzian moments (30) in the quantum spin dynamics equation (11) one can distinguish two different dynamical situations. Ohmic regime: When the resonant frequency ω 0 and the damping rate Γ of the reservoir coupling is much larger than the spin operators' typical frequency of motion, ω 0 ω, each successive term in expansion (23) is smaller by an extra factor of (ω/ω 0 ). In the limit of infinite ω 0 but finite damping parameter η G , the Lorentzian damping term in (11) thus tends to the Ohmic one (20), Non-Ohmic regime: When the environmental frequency ω 0 is comparable to typical spin motion frequencies, ω 0 ≈ ω the Ohmic approximation to the Lorentzian kernel begins to fail. The first deviation in (11) is a new term proportional to κ Lor 2 . I.e. in addition to the Gilbert damping term, one adds the term γ 2 κ Lor 2Ŝ (t) × ∂ 2 tŜ (t) containing a second time derivative of the spin operator S. By analogy with the classical equation of motion for a massive body, this is known as an "inertial" modification to the spin dynamics [49]. As the ratio κ 2 /κ 1 has the dimensions of time, one may introduce an "inertial timescale" τ in [49], which for the Lorentzian is A large inertial time τ in indicates the presence of non-Markovian dynamics, i.e. dynamics that has a certain degree of memory. For a high quality factor resonance ω 0 Γ the inertial time is (half) the kernel's decay time, τ d = 2/Γ. For increasing resonance width Γ, the inertial time τ in decreases and memory effects become less important. In magnetic systems this timescale determines the time over which nutation oscillations are observed in the precession of the spin. Such inertial corrections to standard magnetism have very recently been observed for the first time in ultrafast experiments on thin films [50]. Curiously, the inertial timescale becomes negative when Γ > ω 0 , a fact that may be the subject of future investigation.
Similarly to the kernel expansion (23), one can also expand the Lorentzian power spectrum in frequency, see Appendix A4. Generally, only the odd moments κ 2m+1 appear in the power spectrum. This implies that if a kernel only has non-trivial first (κ 1 ) and second (κ 2 ) moments, while higher moments vanish (κ m>2 = 0) then the power spectrum will still be given by the (quantum) Ohmic one (21), with η G = −κ 1 . When third or higher moments are non-zero, then the power spectrum of the noise will deviate from the Ohmic case at all temperatures.
To summarise, here we have demonstrated that Lorentzian coupling functions, kernels and power spectra provide a systematic framework to explore system dynamics that arises from inertial terms and other memory effects, while recovering the standard Ohmic limit whenever the Lorentzian resonance frequency ω 0 is much larger than the typical system frequencies.

D. Unit-free variables and Lorentzian parameters
In section III we perform semi-classical simulations of the dynamics of Eq. (11), for the Lorentzian kernel of Sec. II B. For this purpose we re-write expressions in In (e-h) the Lorentzian power spectrum is compared to the quantum LLG approximation (cyan squares), and its high temperature white noise limit (magenta dash-dotted).
the operator equation (11) in terms of a unit-free set of quantities. The time coordinate t is replaced with the unit-free coordinate ω L t, where ω L = |γ||B ext | is the Larmor frequency. In addition the spin operatorŜ with largest eigenvalue S 0 is re-written in terms of a unit-free operatorŝ with largest eigenvalue 1,Ŝ = sign(γ) S 0ŝ . The sign of the gyromagnetic ratio is included in the definition ofŝ so thatŝ aligns with the magnetic field, whatever the sign of γ. From Eq. (11) we can see that the damping kernel has dimensions of magnetic field squared divided by angular momentum, which leads us to identify a unit-free damping kernel Similarly, the unit-free coupling function c ω is defined through C ω = |B ext | S −1/2 0 c ω and the unit-free spectral functionsp throughP = |B ext | 2 S −1 0 ω −1 Lp . Looking at the Lorentzian kernel K(t − t ) in (27), the pulling of dimensions can be achieved by setting the kernel amplitude to A = |B ext | 2 S −1 0 α where α now is a frequency. For the simulations we choose the Lorentzian parameters ω 0 , Γ and α to be independent of spin length S 0 . Through (5) this implies an interaction energy scaling ofV int ∝ S 0 √ A ∝ √ S 0 which sets the scaling of the interaction versus self-energy tô V int /Ĥ S ∝ 1/ √ S 0 . Apart from it being implied by dimensional analysis, such scaling is heuristically plausible in many physical contexts. E.g. it is similar to the increasing ratio of the surface (where reservoir interaction occurs) to volume (self-energy) for decreasing system sizes. For microscopic systems for whichV int is no longer small in comparison toĤ S a thermodynamic treatment beyond the weak coupling limit [26,27,81], a limit tacitly assumed in standard thermodynamics, may be required.
Similarly for Ohmic coupling leading to the LLG equa-tion (1), the above scaling choice amounts to choosing η G ∝ A ∝ 1/S 0 implying that η = η G γ 2 S 0 = αΓω 2 L /ω 4 0 is assumed to be independent of S 0 . In physical situations where this assumption is not justified, one may instead choose η and α to depend on the spin length S 0 .
For a spin in an external field B ext the typical frequency of the dynamics is set by the Larmor frequency, ω L . In the simulations discussed in section III we will use the following two sets of Lorentzian parameters, all expressed in terms of ω L , In the second rows we have also listed the equivalent unit-free Gilbert damping η, the inertial timescale τ in , and the memory kernel decay time τ d for the Lorentzian kernel. Figs. 2-5, show plots obtained with Lorentzian coupling functions (26) with parameter Sets 1 and 2, which are shown in blue and red, respectively. The Ohmic LLG approximation is shown in magenta when the classical reservoir (22) is considered, and in cyan when the quantum reservoir (21) is considered. The external field is set to B ext = 10 T e z with e z the unitvector in z-direction throughout. Parameter Set 1 has been chosen to have a resonant frequency ω 0 much larger than the characteristic spin precession frequency ω L (Ohmic regime). Consequently we can truncate the series given in Eq. (23) to leading order and recover the Ohmic form (20) typically considered in magnetism theory. The validity of this approximation is demonstrated in the top row of Fig. 2. Fig. 2a) shows that the Lorentzian coupling function c Lor ω is well approximated by LLG (Ohmic) coupling for the relevant frequency range, while Fig. 2c) shows that the kernel is approximately instantaneous on the timescale ω −1 L , in line with Ohmic damping for which κ Ohm m>1 = 0. Fig. 2e) shows that at high temperature (T = 200 K) and for relevant frequencies ω ≈ ω L , the power spectrump Lor is well approximated by the quantum Ohmic (LLG) power spectrum (21), and its classical limit (22). Fig. 2g) shows that at lower temperatures (T = 1 K) quantum noise becomes important. Here the Ohmic approximation (21) remains valid, while its classical limit (22) is invalid.
Parameter Set 2 is chosen such that it has a resonant frequency ω 0 comparable to the precession frequency ω L (non-Ohmic regime). Here it is inaccurate to truncate the series (23) and the damping will be fundamentally non-Ohmic. To directly compare with Ohmic dynamics generated by Lorentzian coupling with Set 1, both parameter sets have been chosen to correspond to the same unit-free Gilbert damping parameter, η = 0.02. The failure of the Ohmic approximation is demonstrated in the bottom row of Fig. 2. Fig. 2b) shows that the linear approximation to the coupling function fails in the relevant frequency range. For this set of parameters the damping kernel now exhibits significant memory and Fig. 2d) shows that the response persists over a timescale of several ω −1 L . This memory kernel implies, through the FDT (16), a coloured quantum noise power spectrump Lor qu . As shown in Fig. 2h), this coloured Lorentzianp Lor qu (red) differs from the LLG counterpart, p Ohm qu (cyan). Furthermore, Fig. 2f) shows that in the high temperaturesp Lor qu (red) also differs very significantly from the LLG power spectrump Ohm cl (magenta). The presence of both memory and coloured quantum noise for the Lorentzian with parameter Set 2 are both signatures of a thermostat that substantially deviates from the classical Ohmic assumptions and, as we will see in the next section, leads to markedly different short time dynamics and steady state of s z .

III. Semi-classical spin dynamics simulations
The general spin dynamics equation (11) is an operator equation for quantum spins in a lattice, each interacting with neighbouring spins and with a bosonic reservoir. Solving the quantum dynamics using, for example, Lorentzian coupling, kernel and spectrum, is rather difficult without approximations, even numerically, and such exploration is left for future work.
To make progress here, we will numerically solve the full non-Markovian for a semi-classical version of Eq. (11), while including coloured quantum noise and memory effects arising from the coupling to the environment. It replaces the quantum spin operatorŜ with a classical spin vector S, and the quantum stochastic noise field vectorb with a stochastic classical field vector b with statistics that obey the quantum fluctuationdissipation theorem (16a). This semi-classical approach is currently used in the theory of molecular and ionic dynamics [82,83], and was perhaps first applied by Koch [84] to include the effects of quantum fluctuations in Josephson junctions. It has been justified through an expansion of a forward-backward path integral [85,86] (note the remark of Caldeira and Leggett on pg. 589 of [35]), and is valid when the potential energy can be expanded in the path integral to first order in the deviations from the average path. The validity of applying this approach to the decay of metastable states was investigated in detail in [87].
Here we simulate a single spin allowing us to illustrate the impact of memory effects and the reservoir's quantum statistics on the spin dynamics and steady state. The simulation details presented below can readily be extended to multiple interacting spins and could be integrated in sophisticated atomistic spin dynamics simulations such as those used in [47,60].
A. How to simulate coloured noise and memory kernel Here we detail how to efficiently simulate non-Markovian dynamics that arises as a result of Lorentzian coupling (26) for the example of spins vectors. Numerical implementation of (11) requires both -the integration of the kernel with the spin state of previous time steps and the inclusion of coloured noise as follows.
Dropping the spin index and for simplicity assuming any isotropic kernel K(τ ) = 1 3 K(τ ), the three vector components b j (t) for j = 1, 2, 3 of the magnetic noise (12) are implemented as [88] b where ξ j (t ) is standard white Gaussian noise for the jth component, which is delta correlated ξ j (t) ξ k (t ) = δ jk δ(t−t ). The "coloured noise" comes from choosing F (t − t ) as the Fourier transform of the square root of the power spectrum associated with the kernel through (16), i.e.
which can be implemented using a fast Fourier transform. To simulate the effect of a Lorentzian damping kernel (27) we numerically integrate [89] the following set of first order coupled differential equations for the spin vector S and two dummy vectors V and W : The integrated values of the dummy vectors and the spin are separated by the time step dt. Solving these equations is equivalent to solving the integro-differential equation (11) for a Lorentzian kernel, see Appendix A5, but is numerically more straightforward to implement.
B. Single trajectories for different couplings and noises. We wish to illustrate on a single trajectory level, the differences between the dynamics predicted by (11)   first, because the dynamics is intrinsically stochastic, trajectories will naturally differ and cannot readily be compared. However, looking at the noise generation in Eqs. (35) and (36), one can see that the same white noise ξ j (t) for j = 1, 2, 3 may be used as a seed to create comparable "stochastic" noise for different power spectraP (ω). Fig. 3 shows the stochastic short time dynamics of a single classical spin for two pairs of spin length and temperature, S 0 = 1 /2 at T = 1 K (left panel) for a single electron, and S 0 = 200 /2 at T = 200 K (right panel) for a mesoscopic cluster of spins with a combined larger effective spin.
The dynamics is obtained according to Eq. (11) for Lorentzian coupling (26) with S 0 -scaling A = |B ext | 2 S −1 0 α, for parameter sets Set 1 (top panel, blue) and Set 2 (bottom panel, red), and with the quantum coloured noise given by (29). For comparison we also show the short time dynamics according to the LLG equation (1) with S 0 -scaling η G = AΓ/ω 4 0 = |B ext | 2 S −1 0 ω −2 L η with the Gilbert damping parameter η = 0.02 common to both Lorentzian parameter sets, see (34). That implies that the top and bottom LLG plots are identical. For the standard LLG equation two types of noise are considered -high-temperature classical noise (magenta) see Eq. (22), and quantum noise (cyan) see Eq. (21). Since the same white noise time series is used as a seed for producing the stochastic noise for all traces, we can compare them directly. We will here focus on s z , the component of s = sign(γ)S/S 0 aligned with the external field B ext .
Three features stand out in Fig. 3: i) as expected from section II, the dynamics generated with Eq. (11) with Lorentzian Set 1 (top, blue) closely matches the dynamics obtained with the LLG equation (1) with quantum noise (cyan) for both spin-temperature pairs, ii) the quantum statistics of the reservoir (cyan) at low temperatures (left) introduces differences to the LLG dynamics compared to the LLG dynamics obtained with classical noise (magenta), and iii) memory effects that are present for Lorentzian Set 2 (bottom, red) result in significantly different dynamics from that arising with the memory-free Lorentzian Set 1 (top, blue).
We remark that due to the spin/temperature ratio being the same for the two spin-temperature pairs, the LLG equation with classical noise (magenta) integrates to exactly the same dynamics in left and right panel, see Appendix A8. This scaling relation ceases to be true for the LLG equation with quantum noise (cyan). Another difference to note is that in Fig. 3a-d) the dynamics for Lorentzian parameter Set 1 (blue) varies more rapidly in time than for Lorentzian parameter Set 2 (red). This is due to the high frequency content of Set 1's power spectrum, see Fig. 2e+g).
Finally, the spin component s x (green) and the spinvector length |s| (black) are shown for the Lorentzian coupling with Set 1 and Set 2 in the top and bottom panels of Fig. 3, respectively. The plots of |s| show that the numerical integration of Eq. (11) indeed leads to a constant spin-vector length |s| = 1, i.e. no renormalisation is required.
C. Ensemble-averaged s z trajectories. Fig. 4 shows the ensemble averaged s z over time, averaged over 500 stochastic trajectories. We now highlight two important features in Fig 4. Firstly, at low temperatures (left) the quantum statistics of the reservoir (blue, red, cyan) results in a much depleted value of s z , roughly at around 0.28, in comparison to that obtained with the LLG equation with classical noise (magenta), ca 0.85. This indicates that for this particular choice of spin length and temperature the quantum character of the reservoir strongly affects the value of s z , as further discussed below, and the classical hightemperature limit taken in (16b) would not be appropriate. For the high temperature T = 200 K + larger spin pair S 0 = 200 /2 (right), the difference between classical and quantum statistics of the reservoir can be neglected and s z settles at 0.85 independent of whether the spin dynamics integration was done for Eq. (11) with either Set 1 (blue) or Set 2 (red), or for Eq. (1) with either classical (magenta) or quantum noise (cyan).
Secondly, for the large spin-temperature pair (right), there is clear evidence of a much quicker relaxation to steady state (by a factor of a third) for Lorentzian Set 2 (red) compared to the other plots (blue, cyan, magenta). This is a non-Markovian effect that arises because the memory kernel for Set 2 has an appreciable memory over time, see Fig. 2d), while the other memory kernels are (close to) instantaneous. This quicker equilibration occurs because the non-Markovian kernel leads to a smoother dynamics which in turn is more quickly sampled by the dynamical system. D. Steady state s z as a function of temperature.
. Secondly, for simulations that include the full quantum noise (cyan, blue, red) in the dynamics of the small spin at low temperatures (left), we observe reduced s z values in the range 0.2-0.4 at T = 0 K, i.e., well below the classical value of 1. This arises because the power spectrumP Lor qu , given through the quantum FDT (16a), includes quantum fluctuations which remain even for T → 0 K. The steady state curves with quantum noise also show a characteristic "flattening" compared to the steep decay of the steady state with temperature for classical noise (magenta). Qualitatively speaking, this quantum zero point noise, when compared to classical noise, is as if thermal noise is "on" even at 0 K. I.e. taking the Larmor frequency as the relevant frequency, and settingP Ohm qu (0 K) =P Ohm cl (T cl ) for the Ohmic coupling, for example, defines a classical temperature of T cl = 6.7 K "equivalent" to the quantum zerotemperature case. For S 0 = 1 /2 the classical statistical physics steady state value at T cl is s z stat phys ≈ 0.3. This indeed is of comparable size to the s z values obtained with quantum noise at T = 0 K. The corresponding steady state value for the large spin (S 0 = 200 /2) is s z stat phys ≈ 0.995 ≈ 1, see Fig. 5c).
Generally, for the classical temperature T cl = ω L 2k B which is "equivalent" to the quantum zero-point noise, one obtains s z stat phys = coth 2S0 − 2S0 , which only depends on the spin length S 0 while being independent of the field strength |B ext |. With increasing S 0 this function rises very sharply from ≈ 0.3 to 1. For example for spin length S 0 = 5 /2, the quantum zero temperature value is s z stat phys ≈ 0.8 and its decay with increasing temperate is shown in Fig. 6b) in Appendix A7. The middle panel, Fig. 5b), gives an alternative illustration of the steady state value for the small spin as a function of environment temperature T . It shows the same plot as panel a), but with the y-axis rescaled as m(T ) = s z (T ) / s z (0) . All plots now start at 1 at T = 0 K, independent of whether the dynamics was integrated with quantum or classical noise. Interestingly, the overarching behaviour of the resulting curves bears some resemblance with heuristically rescaled magnetization curves that match experimental data [90,91]. Running high-end atomistic simulations of Eq. (11), instead of (1), for multiple interacting classical spins would answer if the quantum power spectrum's impact on their low temperature magnetisation behaviour as well as their Curie temperature is the reason for the apparent rescaling. We note that larger numerical uncertainties arise for the quantum noise because an additional scale is present in comparison to classical noise, see Appendix A8. Error bars obtained from an ensemble of simulations are indicated at one temperature value in a) and c) for all four s z curves in Fig. 5.

IV. Conclusions and open questions
We have derived a general quantum spin dynamics equation, Eq. (11), capable of describing threedimensional precession and damping. The terms arising from the reservoir interaction are treated in a quantum thermodynamically consistent manner, by tracing the origin of both the memory kernel, K(τ ), and the stochastic noise,b(t), to a single coupling function, C ω . Secondly, Lorentzian coupling functions were proposed and shown to provide a systematic means to investigate different dynamical regimes -from Ohmic to non-Ohmic dynamics which is subject to memory and coloured noise. We showed that only in the Ohmic regime, the standard LLG equation with Gilbert damping, widely used in magnetism, is recovered. Finally, we provided details of how to include Lorentzian memory and coloured noise in numerical simulations of open sys-tem dynamics. For the example of a single spin vector, we illustrated that a non-Ohmic Lorentzian kernel leads to a faster equilibration time of s z in comparison to the Ohmic (LLG) regime. We also discussed the steady state differences that arise when the full quantum thermostat with quantum zero-point noise is employed, in contrast to classical white noise.
The above three ingredients provide a complete framework for the simulation of damped threedimensional precession including memory and coloured noise. It can readily be adapted in atomistic spin dynamics simulations [47,60] that solve the dynamics of millions of interacting spins.
The theory presented here will be a useful tool for investigating non-Markovian behaviour, opening up a number of avenues for future research at the intersection of quantum thermodynamics, magnetic materials and beyond. For example, it is an open question to clarify the connection between the three-dimensional precession described by the spin equation (11) and rotational Brownian motion. The orientation of a nonsymmetric rotating body behaves analogously to the three-dimensional spin vector, and the motion and viscosity of a gas surrounding a rotating body simultaneously act on its motion while obeying the FDT as discussed in recent work [92,93].
Within magnetism, for particular materials of interest, detailed models of the coupling functions C ω can be developed that are based on an understanding of the interactions between spins, phonons, and electrons in the material [79]. Coupling to optical modes may further be included to describe, for example, whispering gallery photon-magnon coupling which leads to an effective Gilbert damping term that can take either sign [94]. A direct experimental characterisation of a material's damping kernel K that determines memory and noise in (11) may be attempted, for example with high field experiments such as those recently reported in [50]. Of particular interest are dynamical features beyond the inertial kernel approximation, which will also modify the noise spectrum at larger temperatures.
While we here discussed scalar couplings to the environment in depth, Eq. (11) does hold for any real 3x3 matrix C ω describing the spin-environment interaction in three dimensions. Anisotropic coupling tensors can be implemented, suitable for describing magnetization dynamics within thin films [62], where one direction is coupled differently to environmental modes than the other two. One simplification of our three-dimensional model is to choose a coupling tensor such that spins interact with only one-dimensional environmental modes. This reduces the theory to the spin-boson model, see Appendix A1, whose quantum thermodynamic properties have been discussed very extensively, recently for example in [33].
Microscopic heat transport in spin systems can also be analysed by allowing non-equilibrium situations where individual reservoir modes at frequencies ω and for spins n are thermal -but at different temperatures. This will result in spin dynamics that shuffles energy from one reservoir mode to another, and could result in two-and more-temperature models. For example, the possibility of different phonon modes, each with their own temperature, to couple with different strengths to electrons has recently been analysed in [95] for a magnetic system excited by an ultra-short laser pulse. Furthermore, in deriving the FDT we have assumed bosonic environmental modes but it would be insightful to identify changes to the properties of equation (11) that arise when the spins couple directly to electrons, or fermionic modes in general [14,15,96].
Beyond the quantum character of the reservoir, it will be important to numerically solve the full quantum dynamics according to Eq. (11), including spin operators interacting with neighbouring spin operators. Advanced quantum numerical methods such as Hierarchical Equations Of Motion (HEOM) [97], and the recently proposed time-evolving matrix product operators (TEMPO) method [98] will be required to efficiently describe the time evolution of even just a sin-gle quantum spin coupled to a non-Markovian environment. For multiple interacting spins at low temperatures one can expect entanglement between the spins being present during the short-time dynamics, and even in steady state [99][100][101]. Unfortunately, evaluating such properties will very quickly become a numerically hard problem, requiring advanced numerical techniques such as density-matrix renormalisation group (DMRG) [102] to find realistic approximate solutions. Vice versa, solving (11) within the classical spin vector approximation while including a full quantum power spectrum for the environmental modes, may prove insightful and numerically tractable in the context of finding suitable models for noise in quantum computing hardware, such as superconducting qubits that are held in the mK range [103]. The results may also inform implementations of Young's double slit experiment with a levitated single magnetic domain nanoparticle using the Einsteinde Haas effect [104,105]. but one could also use a continuum description, as in micromagnetics [44].
[111] Note that tensors and vectors are set in calligraphic and bold font, respectively, and that scalar products between vectors are indicated with ·, while a tensor followed by a tensor or a vector is to be understood as matrix multiplication.
[112] Note that in contrast to what is typically done in Caldeira-Leggett type models [35] no counter term has been included here. In any case, coupling to a spin would result in a term proportional to S 2 ∝ 1, which will incur an offset in the overall Hamiltonian that does not affect the dynamics.
[113] The Fourier transformK (n) (ω) of the kernel K (n) (τ ) automatically satisfies the Kramers-Kronig relations [63], connecting the dissipative and reactive parts of the response kernel, as is required for any causal response.
[114] The autocorrelation function of two reservoir operators A and B † in the thermal reservoir stateρR is defined as the expectation value of the Hermitian operator, In the classical case A and B † commute at all times, removing the need for this distinction.
[115] The Fourier transform is here defined as [116] The power spectrum given in (16) is the correct general version for any kernel K (n) (t), fulfilling the full quantum FDT [63]. For a Gilbert damping kernel a power spectrum proportional to ω/(exp ( ω/kBT ) − 1) was given in [60,[74][75][76]. This is missing the quantum ground state contribution of ω/2, which acts as stochastic noise on the spin system even at zero temperature.
[117] Note that although (20) is not an explicitly Hermitian operator, Appendix A2 shows that Eq. (11) is nevertheless equivalent to a Hermitian equation of motion.
R is a decoupled two-dimensional reservoir that can be dropped from the dynamics.
A2. Hermiticity of the quantum spin dynamics equation The quantum spin dynamics equation (11) is not written in an explicitly Hermitian form. The integral term containing the damping kernel K (n) includes an operator product that does not equal its conjugate transposê Nevertheless equation (11) is Hermitian, as it is the time integral that commutes withŜ This can be verified from an observation that Eq. (11) is simply a re-written form of the explicitly Hermitian equation (8). Any confusion can be avoided through re-writing Eq. (11) in an equivalent but explicitly Hermitian form   does not affect the evolution of the spin operator, even though the operator cross product is non-zero. A consequence of this result is that the zeroth order term in the expansion of the damping operator (23) does not contribute to the evolution of the spin operator. We note that Eq. (41) implies that only the sum of all the terms in the damping kernel expansion (23) commutes with the spin operator. When using a truncated form of the expansion (23) we must therefore use the explicitly Hermitian equation of motion (42).
A3.Ŝ 2 is a constant of motion of Eq. (11) To evaluate the derivative of (Ŝ (n) (t)) 2 we first express Eq. (11) in explicitely Hermitian form (42). Dropping site index and time for simplicity, we find where we have applied the angular momentum commutation relations, interchanged indices, and used the anti-symmetric property of jkm . The final line follows from the fact that the spin and the effective magnetic field commute.

A4. Lorentzian power spectrum expansion
Similar to the damping kernel term expansion (23), in moments (30) and time-derivatives, the Lorentzian power spectrum (29) can be expanded in powers of frequency ω, as where we have kept the quantum coth unexpanded. The κ Lor m are the same coefficients as those given in (30). For small frequencies ω the first term in the series (46) dominates and the power spectrum takes the (quantum) Ohmic form where comparison with (21) again shows that −κ Lor 1 is the effective Gilbert damping constant.
Beyond the Ohmic regime, one can see in (46) that only the odd moments κ Lor 2m+1 contribute. Therefore the inertial term κ Lor 2 , which is the first deviation of the damping kernel from Ohmic behaviour, does not change the quantum fluctuations in (16). Only when the third order time derivative of the spin operator contributes significantly to equation (11), will memory effects begin to colour the spectrum away from the (quantum) Ohmic form (21).

A5. Set of equations for kernel simulation
Here we show that the simulation of the kernel in Eq. (11) can be achieved by numerically integrating a set of first order coupled differential equations. We assume a single spin and rewrite Eq. (11) as where we have defined V (t) = γ t t0 dt K(t − t ) S(t ).
Furthermore defining W (t) = dV (t) dt , now leads to a differential equation for W (t): where we have assumed K(0) = 0 andK(0) = 0. Expressing K(t − t ) through its Fourier transformK(ω), choosing a Lorentzian kernel (27) and considering the expression Z(t) := dW (t) dt +ΓW (t)+ω 2 0 V (t), we obtain Rearranging gives as stated in the main text. (Note that the assumption K(0) = 0 andK(0) = 0 is fulfilled for the Lorentzian kernel, (28), since the Heaviside function Θ(τ ) = 1 for τ > 0, and zero elsewhere.) However, for the quantum Ohmic power spectrum the components of the stochastic noise can be written as Clearly, in the quantum case the temperature T now appears separately from spin length S 0 , thus introducing an additional scale in comparison to the classical case. Moreover the fact that the frequency integration for the stochastic field does not simplify as in (57) means that relaxation to the steady state at low temperatures (where the coth x cannot be approximated as 1/x) will be much more noisy than in the high temperature case. Thus in our simulations, this additional scale leads to larger uncertainties in the steady state results, as seen in Fig. 5a).
Finally, we note that for the integration of the quantum Ohmic power spectrum in (58) we have introduced a frequency cut-off ω c by hand, which is necessary at low temperatures to avoid the integral diverging. At low temperatures, this cut-off will set an additional, somewhat artificial, scale of the problem. Importantly, such cut-off is not required for the Lorentzian coupling since the power spectrum (29) decays at high frequencies, even at low T .