Spin dynamics and relaxation in the classical-spin Kondo-impurity model beyond the Landau-Lifschitz-Gilbert equation

The real-time dynamics of a classical spin in an external magnetic field and locally exchange coupled to an extended one-dimensional system of non-interacting conduction electrons is studied numerically. Retardation effects in the coupled electron-spin dynamics are shown to be the source for the relaxation of the spin in the magnetic field. Total energy and spin is conserved in the non-adiabatic process. Approaching the new local ground state is therefore accompanied by the emission of dispersive wave packets of excitations carrying energy and spin and propagating through the lattice with Fermi velocity. While the spin dynamics in the regime of strong exchange coupling J is rather complex and governed by an emergent new time scale, the motion of the spin for weak J is regular and qualitatively well described by the Landau-Lifschitz-Gilbert (LLG) equation. Quantitatively, however, the full quantum-classical hybrid dynamics differs from the LLG approach. This is understood as a breakdown of weak-coupling perturbation theory in J in the course of time. Furthermore, it is shown that the concept of the Gilbert damping parameter is ill-defined for the case of a one-dimensional system.


I. INTRODUCTION
The Landau-Lifshitz-Gilbert (LLG) equation 1-3 has originally been considered to describe the dynamics of the magnetization of a macroscopic sample. Nowadays it is frequently used to simulate the dynamics of many magnetic units coupled by exchange or magnetostatic interactions, i.e., in numerical micromagnetics. 4 The same LLG equation can be used on an atomistic level as well. [5][6][7][8][9] For a suitable choice of units and for several spins S m (t) at lattice sites m, it has the following structure: It consists of precession terms coupling the spin at site m to an external magnetic field B and, via exchange couplings J mn , to the spins at sites n. Those precession terms typically have a clear atomistic origin, such as the Ruderman-Kittel-Kasuya-Yoshida (RKKY) interaction [10][11][12] which is mediated by the magnetic polarization of conduction electrons. The non-local RKKY couplings J mn = J 2 χ mn are given in terms of the elements χ mn of the static conduction-electron spin susceptibility and the local exchange J between the spins and the local magnetic moments of the conduction electrons. Other possibilities comprise direct (Heisenberg) exchange interactions, intra-atomic (Hund's) couplings as well as the spin-orbit and other anisotropic interactions. The relaxation term, on the other hand, is often assumed as local, α mn = δ mn α, and represented by purely phenomenological Gilbert damping constant α only. It describes the angular-momentum transfer between the spins and a usually unspecified heat bath.
On the atomistic level, the Gilbert damping must be seen as originating from microscopic couplings of the spins to the conduction-electron system (as well as to lattice degrees of freedom which, however, will not be considered here). There are numerous studies where the damping constant, or tensor, α has been computed numerically from a more fundamental model including electron degrees of freedom explicitly [13][14][15] or even from first principles. [16][17][18][19][20][21] All these studies rely on two, partially related, assumptions: (i) The spin-electron coupling J is weak and can be treated perturbatively to lowest order, i.e., the Kubo formula or linear-response theory is employed. (ii) The classical spin dynamics is slow as compared to the electron dynamics. These assumptions appear as well justified but they are also necessary to achieve a simple effective spin-only theory by eliminating the fast electron degrees of freedom.
The purpose of the present paper is to explore the physics beyond the two assumptions (i) and (ii). Using a computationally efficient formulation in terms of the electronic one-particle reduced density matrix, we have set up a scheme by which the dynamics of classical spins coupled to a system of conduction electrons can be treated numerically exactly. The theory applies to arbitrary coupling strengths and does not assume a separation of electron and spin time scales. Our approach is a quantum-classical hybrid theory 22 which may be characterized as Ehrenfest dynamics, similar to exact numerical treatments of the dynamics of nuclei, treated as classical objects, coupled to a quantum system of electrons (see, e.g., Ref. 23 for an overview). Some other instructive examples of quantum-classical hybrid dynamics have been discussed recently. 24,25 The obvious numerical advantage of an effective spinonly theory, as given by LLG equations of the form (1), is that in solving the equations of motion there is only the time scale of the spins that must be taken care of. As compared to our hybrid theory, much larger time steps and much longer propagation times can be achieved. Opposed to ab-initio approaches 16,17,26 we therefore consider a simple one-dimensional non-interacting tightbinding model for the conduction-electron degrees of freedom, i.e., electrons are hopping between the nearestneighboring sites of a lattice. Within this model approach, systems consisting of about 1000 sites can be treated easily, and we can access sufficiently long time scales to study the spin relaxation. An equilibrium state with a half-filled conduction band is assumed as the initial state. The subsequent dynamics is initiated by a sudden switch of a magnetic field coupled to the classical spin. The present study is performed for a single spin, i.e., we consider a classical-spin Kondo-impurity model with antiferromagnetic local exchange coupling J, while the theory itself is general and can be applied to more than a single or even to a large number of spins as well.
As compared to the conventional (quantum-spin) Kondo model, 27,28 the model considered here does not account for the Kondo effect and therefore applies to situations where this is absent or less important, such as for systems with large spin quantum numbers S, strongly anisotropic systems or, as considered here, systems in a strong magnetic field. To estimate the quality of the classical-spin approximation a priori is difficult. [29][30][31] For one-dimensional systems, however, a quantitative study is possible by comparing with full quantum calculations and will be discussed elsewhere. 32 There are different questions to be addressed: For dimensional reasons, one should expect that linearresponse theory, even for weak J, must break down at long times. It will therefore be interesting to compare the exact spin dynamics with the predictions of the LLG equation for different J. Furthermore, the spin dynamics in the long-time limit can be expected to be sensitively dependent on the low-energy electronic structure. We will show that this has important consequences for the computation of the damping constant α and that α is even ill-defined in some cases. An advantage of a full theory of spin and electron dynamics is that a precise microscopic picture of the electron dynamics is available and can be used to discuss the precession and relaxation dynamics of the spin from another, namely from the electronic perspective. This information is in principle experimentally accessible to spin-resolved scanning-tunnelling microscope techniques [33][34][35][36] and important for an atomistic understanding of nano-spintronics devices. 37,38 We are particularly interested in the physics of the system in the strong-J regime or for a strong field B where the time scales of the spin and the electron dynamics become comparable. This has not yet been explored but could become relevant to understand real-time dynamics in realizations of strong-J Kondo-lattice models by means of ultracold fermionic Yb quantum gases trapped in optical lattices. 39,40 The paper is organized as follows: We first introduce the model and the equations of motion for the exact quantum-classical hybrid dynamics in Sec. II and discuss some computational details in Sec. III. Sec. IV provides a comprehensive discussion of the relaxation of the classical spin after a sudden switch of a magnetic field. The reversal time as a function of the interaction and the field strength is analyzed in detail. We then set the focus on the conduction-electron system which induces the relaxation of the classical spin by dissipation of energy. In Sec. V, the linear-response approach to integrate out the electron degrees of freedom is carefully examined, including a discussion of the additional approximations that are necessary to re-derive the LLG equation and the damping term in particular. Sec. VI summarizes the results and the main conclusions.

II. MODEL AND THEORY
We consider a classical spin S with |S| = 1/2, which is coupled via a local exchange interaction of strength J to the local quantum spin s i0 at the site i 0 of a system of N itinerant and non-interacting conduction electrons. The conduction electrons hop with amplitude −T (T > 0) between non-degenerate orbitals on nearest-neighboring sites of a D-dimensional lattice, see Fig. 1. L is the number of lattice sites, and n = N/L is the average conduction-electron density.
The dynamics of this quantum-classical hybrid system 22 is determined by the Hamiltonian Here, c iσ annihilates an electron at site i = 1, ..., L with spin projection σ =↑, ↓, and s i = 1 2 σσ c † iσ σ σσ c iσ is the local conduction-electron spin at i, where σ denotes the vector of Pauli matrices. The sum runs over the different ordered pairs ij of nearest neighbors. B is an external magnetic field which couples to the classical spin.
To be definite, an antiferromagnetic exchange coupling J > 0 is assumed. If S was a quantum spin with S = 1/2, Eq. (2) would represent the single-impurity Kondo model. 27,28 However, in the case of a classical spin considered here, there is no Kondo effect. The semiclassical single-impurity Kondo model thus applies to systems where a local spin is coupled to electronic degrees of freedom but where the Kondo effect absent or suppressed. This comprises the case of large spin quantum numbers S, or the case of temperatures well above the Kondo scale, or systems with a ferromagnetic Kondo coupling J < 0 where, for a classical spin, we expect a qualitatively similar dynamics as for J > 0.
We assume that initially, at time t = 0, the classical spin S(t = 0) has a certain direction and that the conduction-electron system is in the corresponding ground state, i.e., the conduction electrons occupy the lowest N one-particle eigenstates of the non-interacting Hamiltonian Eq. (2) for the given S = S(t = 0) up to the chemical potential µ. A non-trivial time evolution is initiated if the initial direction of the classical spin and the direction of the field B are non-collinear.
To determine the real-time dynamics of the electronic subsystem, it is convenient to introduce the reduced oneparticle density matrix. Its elements are defined as expectation values, in the system's state at time t. At t = 0 we have ρ(0) = Θ(µ − T (0)). The elements of ρ(0) are given by where Θ is the step function and where U is the unitary matrix diagonalizing the hopping matrix T (0), i.e., U † T (0)U = ε with the diagonal matrix ε given by the eigenvalues of T (0). The hopping matrix at time t is can be read off from Eq. (2). It comprises the physical hopping and the contribution resulting from the coupling term. Its elements are given by Here δ ii = 1 if i, i are nearest neighbors and zero else. There is a closed system of equations of motion for the classical spin vector S(t) and for the one-particle density matrix ρ(t). The time evolution of the classical spin is determined via (d/dt)S(t) = {S, H class. } by the classical Hamilton function H class. = H . This equation of motion is the only known way to consistently describe the dynamics of quantum-classical hybrids (see Refs. 22,41,42 and references therein for a general discussion). The Poisson bracket between arbitrary functions A and B of the spin components is given by, 43,44 where the sums run over x, y, z and where ε αβγ is the fully antisymmetric ε-tensor. With this we find This is the Landau-Lifschitz equation where the expectation value of the conduction-electron spin at i 0 is given by and where J s i0 t acts as an effective time-dependent internal field in addition to the external field B. The equation of motion for s i t reads as d dt

S(t)
where the sum runs over the nearest neighbors of i. The second term on the right-hand side describes the coupling of the local conduction-electron spin to its environment and the dissipation of spin and energy into the bulk of the system (see below). Apparently, the system of equations of motion can only be closed by considering the complete one-particle density matrix Eq. (3). It obeys a von Neumann equation of motion, as is easily derived, e.g., from the Heisenberg equation of motion for the annihilators and creators. As is obvious from the equations of motion, the realtime dynamics of the quantum-classical Kondo-impurity model on a lattice with a finite but large number of sites L can be treated numerically exactly (see also below). Nevertheless, the model comprises highly non-trivial physics as the electron dynamics becomes effectively correlated due to the interaction with the classical spin. In addition, the effective electron-electron interaction mediated by the classical spin is retarded: electrons scattered from the spin at time t will experience the effects of the spin torque exerted by electrons that have been scattered from the spin at earlier times t < t.

III. COMPUTATIONAL DETAILS
Eqs. (5), (7), (8) and (10) represent a coupled nonlinear system of first-order ordinary differential equations which can be solved numerically. By blocking up the von Neumann equation (10), the differential equations are written in a standard formẏ = f (y(t), t), where y(t) is a high-dimensional vector, such that an explicit Runge-Kutta method can be applied. A high-order propagation technique is used 45 which provides the numerically exact solution up to 6-th order in the time step ∆t. For a typical system consisting of about L = 10 3 sites this implies that ∼ 10 6 coupled equations are solved.
We consider a one-dimensional system with open boundaries consisting of L = 1001 sites and a local perturbation at the central site i 0 of the system, see Fig. 1. For a half-filled tight-binding conduction band the Fermi velocity v F = 2T roughly determines the maximum speed of the excitations and defines a "light cone". 46,47 This means that finite-size effects due to scattering at the system boundaries become relevant after a propagation time t max ∼ 500 (in units of 1/T ). A time step ∆t = 0.1 is usually sufficient for reliable numerical results up to t max , i.e., about 5000 time steps are performed. The computational cost is moderate, and calculations can be performed in a few hours on a standard desktop computer.
Assuming, for example, that B = (0, 0, B), the Hamiltonian is invariant under rotations around the z axis. It is then easily verified that not only the length of the spin |S| = 1/2 is conserved but also the total number of conduction electrons, the z-component of the total spin, as well as the total energy, Despite the fact that the model does not include a direct (e.g., Coulomb) interaction among the conduction electrons, the average occupation numbers of the basis of one-particle states in which the hopping matrix T (t) is diagonal at time t are not conserved. This is due to the effective retarded interaction mediated by the classical spin. Hence, the system is not integrable, unlike a free fermion gas. The conservation of the above-mentioned global observables serves as a sensitive check for the accuracy of the numerical procedure.

IV. COUPLED SPIN AND ELECTRON DYNAMICS
A. Spin relaxation where a non-zero but small polar angle ϑ = π/50 is necessary to slightly break the symmetry of the initial state and to start the dynamics.
For the same setup, the Landau-Lifschitz-Gilbert equation would essentially predict two effects: first, a precession of the classical spin around the field direction with Larmor frequency ω L = B z , and second, a relaxation of the spin to the equilibrium state with S parallel to B = Be z for t → ∞. Both effects are also found in the full dynamics of the quantum-classical hybrid model. The frequency of the oscillation of S x (t) that is seen in Fig. 2 is ω L , and S z is reversed after a few hundred time units. The precessional motion is easily explained by the torque on the spin exerted by the field according to Eq. (7). The explanation of the damping effect is more involved: Even for the high field strength considered here, the spin dynamics is slow as compared to the characteristic electronic time scale such that it could be reasonable to assume the electronic system being in its instantaneous  ground state at any instant of time and corresponding to the configuration of the classical spin. This, however, would imply that the expectation value s i0 t of the local conduction-electron moment at i 0 is always strictly parallel to S(t) and, hence, there would not be any torque mediated by the exchange coupling J on S(t).
In fact, the direction of s i0 t is always somewhat behind the "adiabatic direction", i.e., behind −S(t): This is shown in Fig. 3 for a field of strength B = 0.1 where the spin dynamics is by a factor 10 slower, compared to Fig. 2, and for different stronger exchange couplings J. Even in this case the process is by no means adiabatic, and the angle γ(t) between S(t) and s i0 t is close to but clearly smaller than γ = π at any instant of time. This non-adiabaticity results from the fact that the motion of the classical spin affects the conduction electrons in a retarded way, i.e., it takes a finite time until the local conduction-electron spin s i0 t at i 0 reacts to the motion of the classical spin.
This retardation effect results in a torque J s i0 t × S(t) = 0 exerted on the classical S(t) in the +z direction which adds to the torque due to B and which drives the spin to its new equilibrium direction. Hence, retardation is the physical origin of the Gilbert spin damping.
With increasing time, the deviation of γ(t) from the instantaneous equilibrium value γ = π increases in magnitude, i.e., the direction of s i0 t is more and more behind the adiabatic direction, and the torque increases. Its magnitude is at a maximum at the same time when the oscillating x, y components of S(t) are at a maximum (see Fig. 2). The z-component of the torque does not vanish before the spin has reached its new equilibrium position S(t) ∝ e z . Fig. 3 shows results for different J. Generally, nonadiabatic effects show up if the typical time scale of the dynamics is faster than the relaxation time, i.e., the time necessary to transport the excitation away from the location i 0 where it is created initially. Roughly this time scale is set by the inverse hopping 1/T . One therefore expects that, for fixed T , a stronger J implies a stronger retardation of the conduction-electron dynamics. The results for different J shown in Fig. 3 in fact show that the maximum deviation of γ(t) from the adiabatic direction γ = π increases with increasing J (for very strong J the dynamics becomes much more complicated, see below). This results in a stronger torque on S(t) in z direction and thus in a stronger damping. The picture is also qualitatively consistent with the LLG equation as the Gilbert damping constant α increases with J.
The spin (almost) reverses its direction after a finite reversal time τ which is shown in Fig. 4 as a function of J. Calculations have been performed for an initial direction of the classical spin with S z (0) = −(1/2) cos ϑ with two different polar angles ϑ 1,2 . If ϑ is sufficiently small, the results for different ϑ are expected to differ by a J and B independent constant factor only. As Fig. 4 demonstrates, the ratio τ 1 /τ 2 is in fact nearly constant. For weak J and up to coupling strengths of about J 30, we find that the reversal time decreases with increasing J. With increasing J, the retardation effect increases, as discussed above, and the stronger damping results in a shorter reversal time. The prediction of the LLG equation for the reversal time of a single spin τ can be derived analytically 48 and is given by However, down to the smallest J for which τ can be calculated reliably, our results for the full spin dynamics do not scale as τ ∝ 1/J 2 as one would expect for weak J assuming that α ∝ J 2 (see discussion in Sec. V B). The B dependence of the reversal time is shown in Fig.  5. With increasing field strength, the classical spin S(t) precesses with a higher Larmor frequency ω L ≈ B around the z axis. Hence, non-adiabatic effects increase. The stronger the field, the more delayed is the precessional motion of the local conduction-electrons spin s i0 t . This results in a stronger torque in +z direction exerted on the classical spin. Therefore, the relaxation is faster and the reversal time τ smaller. For weak and intermediate field strengths, τ is roughly proportional to 1/B. This is consistent with the prediction of the LLG equation, see Eq. (14).
In the limit of very strong fields one would expect an increase of the reversal time with increasing B since the field term will eventually dominate the dynamics, i.e., only the precessional motion survives which implies a diverging reversal time. In fact, for a field strength exceeding a critical strength B c , which depends on J, there is no full relaxation any longer, and τ = ∞. This strong-B regime cannot be captured by the LLG equation and deserves further studies.
The strong-J regime is interesting as well. For coupling strengths exceeding J ≈ 30 the reversal time increases with J (see Fig. 4). Eventually, the reversal time must even diverge. This is obvious as the dynamics is described by a simple two-spin model in the limit J = ∞ which cannot show spin relaxation. The corresponding equations of motion are obtained by from Eqs. (7) and (9) by setting t = 0 and B = 0: Note that we have | s i0 t | = 1/2 for J → ∞. The equations are easily solved by exploiting the conservation of the total spin S tot = S(t) + s i0 t . Both spins precess with constant frequency ω 0 = JS tot around S tot . Their components parallel to S tot are equal, and their components perpendicular to S tot are anti-parallel and of equal length. However, the two-spin dynamics of the J = ∞ limit is not stable against small perturbations. Fig. 6 shows the classical spin dynamics of the full model (with T = 1 and B = 0.1) for a very strong but finite coupling J = 100. Here, the motion of the classical spin gets very complicated as compared with the highly regular  Fig. 2). In particular, the z-component of S(t) is oscillating on nearly the same scale as the x and y components. This characteristic time scale ∆t ≈ 4 of the oscillation corresponds to a frequency ω ≈ 1.5 which differs by more than an order of magnitude from both, the Larmor frequency B = 0.1 and from the exchange-coupling strength J = 100. Note that the oscillation of S z (t) is actually the reason for the ambiguity in the determination of the precise reversal time and gives rise to the error bars in Fig. 4 for strong J.
We attribute the complexity of the dynamics to the fact that the torque due to the field term and the torque due to the exchange coupling are of comparable magnitudes, see second and third panel in Fig. 6. It is interesting that even for strong J, where one would expect S(t) and s i0 t to form a tightly bound local spin-zero state, there is actually a small deviation from perfect antiparallel alignment, i.e., γ(t) = π, as can be seen in the fourth panel of Fig. 6. This results in a finite z-component of the torque on S(t) which leads to a very fast reversal with S z (t) ≈ +0.5 at time t ≈ 2.1. Contrary to the weakcoupling limit, however, the z-component of the torque changes sign at this point and drives the spin back to the −z-direction. At t ≈ 3.2, however, the z-component of S(t) once more reverses its direction. Here, the torque due to the exchange coupling vanishes completely as S(t) and s i0 t are perfectly antiparallel (see the first zero of γ(t) in the fourth panel). The motion continues due to the non-zero field-induced torque. This pattern repeats several times. S(t) mainly oscillates within a plane including and slowly rotating around the z-axis.
Eventually, there is a perfect relaxation of the classical spin for large t but in a very different way as compared to the weak-coupling limit. While the deviation from γ = π is small at any instant of time as for weak J, the most apparent difference is perhaps that the new ground state is approached with an oscillating behavior of γ(t) around γ = π, i.e., s i0 t may run behind or ahead of S(t) as well.
Let us summarize at this point the main differences between the quantum-classical hybrid and the effective LLG dynamics of the classical spin: For weak J and B, the qualitative behavior, precessional motion and relaxation, is the same in both approaches. Quantitatively, however, the LLG equation is inconsistent with the observed J dependence of the reversal time when assuming α ∝ J 2 . The B dependence of 1/τ is linear as expected from the LLG approach. For strong J, the spin dynamics qualitatively differs from LLG dynamics and gets more complicated with a new time scale emerging. Absence of complete relaxation, as observed in the strong-B limit, is also not accessible to an effective spin-only theory.

B. Energy dissipation
To complete the picture of the relaxation dynamics of the classical spin, the discussion should also comprise the dynamics of the electronic degrees of freedom. The spin relaxation must be accompanied by a dissipation of energy and spin into the bulk of the electronic system since the total energy and the total spin are conserved quantities, see Eqs. (12) and (13), while conservation of the total particle number, Eq. (11), is trivially ensured by the particle-hole symmetric setup considered here where the average conduction-electron number at every site is time-independent: σ n iσ t = 1.
The total energy is given by E tot = H , see Eq. (13), and is a sum over different contributions, E tot = E B (t) + E hop (t) + E int (t), namely the interaction energy with the field the kinetic (hopping) energy of the conduction-electron system and the exchange-interaction energy The time dependence of those contributions is shown in Fig. 7 for J = 5 and B = 1.
The top panel of Fig. 7 shows that the system releases the interaction energy |2BS| of the classical spin in the external field by aligning the spin to the field direction. In the long-time limit, this energy is stored in the conduction-electron system: The average kinetic energy per site (L = 1001) increases by the same amount as shown in the second panel (note the different scales). The exchange-interaction energy changes with time but is the same for t = 0 and t → ∞ (see third panel). The total energy is constant (second panel).
The relaxation of the classical spin in the external field B implies an energy flow away from the site i 0 into the bulk of the conduction-electron system such that locally, in the vicinity of i 0 the system is in its new ground state. To discuss this energy flow, it is convenient to consider the total energy as a lattice sum E tot = i e i (t) over the total energy "density" defined as where the sum over j runs over the nearest neighbors of site i. The time dependence of the energy density in the vicinity of i 0 and at distances 50 and 100 is shown in the fourth panel of Fig. 7.
At any site in the conduction-electron system, the energy density increases from its ground-state value, reaches a maximum and eventually relaxes to the energy density of the new ground state. Since the latter is just the ground state with the reversed classical spin, the new ground-state energy density is the same as in the initial state at t = 0. As can be seen in the fourth panel of Fig. 7, there is also a slight spatial oscillation of the ground-state energy density which just reflects the Friedel oscillations around the impurity at i 0 .
Complete relaxation means that the excitation energy is completely removed from the vicinity of i 0 and transported into the bulk of the system. That this is in fact the case can be seen by comparing the energy density at different distances from i 0 . It is also demonstrated by Fig. 8 which visualizes the energy-current density which symmetrically points away from i 0 : The total energy of the excitation flowing through each pair of sites i 0 ±∆i is constant, i.e., the time-integrated energy flux, dt e i (t) is the same for all i.
As is seen in Fig. 7 (fourth panel) and Fig. 8, there is a considerable dispersion of the excitation wave packet carrying the energy. For example, at i 0 + 100 it takes more than four times longer, as compared to i 0 + 1, until most of the excitation has passed through (note the logarithmic time scale).
The broadening of the wave packet, due to dispersion, is asymmetric and bound by an upper limit for the speed of the excitation which is roughly set by the Fermi velocity v F = 2T . This Lieb-Robinson bound 46,47 determines the "light cone" seen in Fig. 8.

C. Spin dissipation
The same upper speed limit, given by the Fermi velocity of the conduction-electron system, is also seen in the spatiotemporal evolution of the conduction-electron spin density s i t . This is shown in Fig. 9 for a different magnetic field strength B = 0.1 where the classical spin dynamics is slower. Apparently, the wave packet of excitations emitted from the impurity not only carries energy but also spin. It symmetrically propagates away from i 0 and, at t ≈ 300, reaches the system boundary where it is reflected perfectly. Up to t = 500 there is hardly any effect visible in the local observables close to i 0 that is affected by the finite system size.
Snapshots of the conduction-electron spin dynamics are shown in Fig. 10 for the initial state at t = 0 and for states at four later times t > 0 which are also indicated by the arrows in Fig. 9. At t = 0 the conduction-electron system is in its ground state for the given initial direction of the classical spin. The latter basically points into the −z direction, apart from a small positive x-component (ϑ = π/50) which is necessary to break the symmetry of the problem and to initiate the dynamics. This tiny effect will be disregarded in the following.
From the perspective of the conduction-electron system, the interaction term JSs i0 acts as a local external magnetic field JS which locally polarizes the conduction electrons at i 0 . Since J is antiferromagnetic, the local moment s i0 points into the +z direction. At half-filling, the conduction-electron system exhibits pronounced antiferromagnetic spin-spin correlations which give rise to an antiferromagnetic spin-density wave structure aligned to the z axis at t = 0, see first panel of Fig. 10.
The total spin S tot = 0 at t = 0, i.e., the classical spin S is exactly compensated by the total conductionelectron spin s tot = i s i = −S in the ground state. This can be traced back to the fact that for a D = 1dimensional tight-binding system with an odd number of sites L, with N = L and with a single static magnetic impurity, there is exactly one localized state per spin projection σ, irrespective of the strength of the impurity potential (here given by JS = 0.5Je z ). The number of ↑ one-particle eigenstates therefore exceeds the number of ↓ states by exactly one.
Since the energy of the excitation induced by the external field B is completely dissipated into the bulk, the state of the conduction-electron system at large t (but shorter than t ≈ 500 where finite-size effects appear) must locally, close to i 0 , resemble the conduction-electron ground state for the reversed spin S = +0.5e z . This implies that locally all magnetic moments s i t must reverse their direction. In fact, the last panel in Fig. 10 (left) shows that the new spin configuration is reached for t = 250 at sites with distance |i − i 0 | 100, see dashed line, for example. For later times the spin configuration stays constant (until the wave packet reflected from the system boundaries reaches the vicinity of i 0 ). The reversal is almost perfect, e.g., s i0 t=0 = 0.2649 → s i0 t≥250 = −0.2645. Deviations of the same order of magnitude are also found at larger distances, e.g., i = i 0 −100. We attribute those tiny effects to a weak dependence of the local ground state on the non-equilibrium state far from the impurity at t = 250, see right part of the last panel in Fig. 10.
The other panels in Fig. 10 demonstrate the mechanism of the spin reversal. At short times (see t = 60, second panel) the perturbation of the initial equilibrium configuration of the conduction-electron moments is still weak. For t = 80 and t = 100 one clearly notices the emission of the wave packet starting. Locally, the antiferromagnetic structure is preserved (see left part) but superimposed on this, there is an additional spatial structure of much longer size developing. This finally forms the wave packet which is emitted from the central region. Its spatial extension is about ∆ ≈ 300 as can be estimated for t = 250 (last panel on the right) where it covers the region 200 i 500. The same can be read off from the upper part of Fig. 9. Assuming that the reversal of each of the conduction-electron moments takes about the same time as the reversal of the classical spin, ∆ is roughly given by the reversal time times the Fermi velocity and therefore strongly depends on J and B. For the present case, we have τ 1 ≈ 150/T which implies ∆ ≈ 150 × 2 = 300 in rough agreement with the data.
In the course of time, the long-wave length structure superimposed on the short-range antiferromagnetic texture develops a node. This can be seen for t = 100 and i ≈ 40 (fourth panel, see dashed line). The node marks the spatial border between the new (right of the node, closer to i 0 ) and the original antiferromagnetic structure of the moments and moves away from i 0 with increasing time.
At a fixed position i, the reversal of the conductionelectron moment s i t takes place in a similar way as the reversal of the classical spin (see both panels in Fig.  9 for a fixed i). During the reversal time, its x and y components undergo a precessional motion while the z component changes sign. Note, however, that during the reversal | s i | gets much larger than its value in the initial and in the final equilibrium state.

A. Perturbation theory
Eqs. (7) and (9) do not form a closed set of equations of motion but must be supplemented by the full equation of motion (10) for the one-particle conduction-electron density matrix. This implies that the fast electron dynamics must be taken into account explicitly even if the spin dynamics is much slower. Hence, there is a strong motivation to integrate out the conduction-electron degrees of freedom altogether and to take advantage from a much larger time step within a corresponding spin-only timepropagation method. Unfortunately, a simple effective spin-only action can be obtained in the weak-coupling (small-J) limit only. 13,14 This weak-coupling approximation is also implicit to all effective spin-only approaches that consider the effect of conduction electrons on the spin dynamics. 49 In the weak-J limit the electron degrees of freedom can be eliminated in a straightforward way by using standard linear-response theory: 50 We assume that the initial state at t = 0 is given by the conduction-electron system in its ground state or in thermal equilibrium and an arbitrary state of the classical spin. This may be realized formally by suddenly switching on the interaction J(t) at time t = 0, i.e., J(t) = JΘ(t) and by switching the local field from some initial value B ini at t = 0 to a final value B for t > 0. The response of the conduction-electron spin at i 0 and time t > 0 ( s i0 t = 0 for t = 0) due to the time-dependent perturbation J(t)S(t) is up to linear order in J. Here, the free (J = 0) local retarded spin susceptibility of the conduction electrons Π (ret) (t, t ) is a tensor with elements where α, β = x, y, z. Using this in Eq. (7), we get an equation of motion for the classical spins only, which is correct up to order J 2 . This represents an equation of motion for the classical spin only. It has a temporally non-local structure and includes an effective interaction of the classical spin at time S(t) with the same classical spin at earlier times t < t. In the full quantum-classical theory where the electronic degrees of freedom are taken into account exactly, this retarded interaction is mediated by a nonequilibrium electron dynamics starting at site i 0 and time t and returning back to the same site i 0 at time t > t . Here, for weak J, this is replaced by the equilibrium and homogeneous-in-time conduction-electron spin susceptibility Π (ret) (t − t ). Compared with the results of the full quantum-classical theory, we expect that the perturbative spin-only theory breaks down after a propagation time t ∼ 1/J at the latest.
Using Wick's theorem, 50 the spin susceptibility is easily expressed in terms of the greater and the lesser equilibrium one-particle Green's functions, G > ii,σσ (t, t ) = −i c iσ (t)c † iσ (t ) and G < ii,σσ (t, t ) = i c † iσ (t )c iσ (t) , respectively: Assuming that the conduction-electron system is characterized by a real, symmetric and spin-independent hopping matrix T ij (as given by the first term of Eq. (2)), G > and G < are unit matrices with respect to the spin indices. They are easily expressed as explicit functions of T (see Ref. 51, for example). With tr (σ α σ α ) = 2δ αα , we find for a conduction-electron system at inverse temperature β and chemical potential µ.
Using Eq. (24) we have computed the spin susceptibility for the ground state (β = ∞) of the conductionelectron system. This fixes the kernel in the integrodifferential equation Eq. (22) which is solved numerically by standard techniques. 52 We again consider the For weak coupling, up to J = 1 (Fig. 11), there is an almost perfect agreement between the results of the exact quantum-classical dynamics (full lines) and the linearresponse theory (dashed lines) up to the maximum propagation time t max = 500. We note that, compared to the full theory, there is a tiny deviation of the linearresponse result for the z component of S(t) visible in Fig. 11 for times t 10. Hence, on this level of accuracy, t 1 ≈ 10 sets the time scale up to which the linear-response theory is valid. This may appear surprising as this implies t 1 J = 10 for the "small" dimensionless parameter of the perturbation theory. One has to keep in mind, however, that even if the perturbation is "strong", its effects can be rather moderate since only non-adiabatic terms ∼ S(t) × S(t ) contribute in Eq. (22).
For J = 3, see Fig. 12, damping of the classical spin sets in much earlier. Visible deviations of the linearresponse theory from the full dynamics already appear on a time scale that is almost two orders of magnitude smaller as compared to the case J = 1. A simple reasoning based on the argument that the dimensionless expansion parameter is t 1 J fails as this disregards the strong enhancement of retardation effects with increas- ing J, which have been discussed in Sec. IV A. These effects make the perturbation much more effective, i.e., lead to a torque, which is exerted by the conduction electrons on the classical spin, growing stronger than linear in J.
We conclude that linear-response theory is highly attractive formally as it provides a tractable spin-only effective theory. On the other hand, substantial discrepancies compared to the full (non-perturbative) theory show up as soon as damping effects become stronger. Note that, with increasing time, these deviations must diminish and disappear eventually since both, the full and the effective theory, predict a fully relaxed spin state for t → ∞ -see lower panel of Fig. 12, for example. At least for simple systems with a single classical spin, as considered here, this implies that the effective theory provides qualitatively reasonable results.

B. Landau-Lifschitz-Gilbert equation
To derive the LLG equation, the linear-response theory must be further simplified: 14,53 As is obvious from Eq. (22), the spin-susceptibility Π (ret) (t−t ) can be interpreted as an effective retarded self-interaction of the spin. We assume that the electron dynamics is much faster than the spin dynamics. On the time scale of the spin dynamics, the self-interaction then takes place almost instantaneously, i.e., the memory kernel Π  (22) is peaked at t ≈ t. We can therefore approximate S(t ) ≈ S(t) + (t − t)Ṡ(t) under the integral in Eq. (22). This immediately yields: where Note that this is a necessary step to arrive at a constant damping parameter which can again be justified by noting that Π (ret) (t) is peaked at t = 0. Before proceeding, let us stress that Eq. (27) is, or is equivalent to, the standard expression used for computing the Gilbert damping constant in various studies: After Fourier transformation, one ends up with which has also been derived, e.g., in Refs. 14,54 in different contexts. The frequency-dependent spin correlation Π (ret) (ω) in Eq. (29) can be obtained explicitly as the Fourier transform of Π (ret) (t) given by Eq. (24). A straightforward calculation yields: where f (ω) = 1/(exp(βω) + 1) is the Fermi function, and for L → ∞ is the local one-particle spectral function.
Here we have assumed periodic boundary conditions, i.e., spatial homogeneity, such that the hopping matrix T is diagonalized by Fourier transformation: with U ik = L −1/2 exp(ikR i ). The eigenvalues of T are given by the tight-binding dispersion ε(k) of the conduction-electron Bloch band. Within our simple tight-binding model for the conduction-electron system, Eqs. (29) and (30) are equivalent with Kambersky's breathing Fermi-surface theory, related torque-correlation models and scattering theory and have frequently been used for ab initio as well as model computations of the Gilbert damping constant. 15,18,20,21,[55][56][57][58] Let us remark that Eq. (30) demonstrates that α < 0. This results from the convention for the coupling of the magnetic field to the spin, namely H = H(B = 0) − BS [see Eq. 2)], which has been adopted here. As a consequence, the precession of S(t) around B is described by a left-hand helix,Ṡ = S(t) × B, and thus α must be negative to describe damping.

C. Ill-defined Gilbert damping
The above discussion shows that Eq. (27) represents the fundamental definition of the Gilbert damping constant α and that the limit t → ∞ is crucial to recover the LLG equation in its standard form. The existence of the long-time limit, however, decisively depends on the longtime behavior of the retarded spin-correlation function Π (ret) (t). Starting from Eq. (24), this is easily computed as Π (ret) (t) = Θ(t)Im 1 L 2 k,p e −iε(k)t 1 + e −β(ε(k)−µ) e iε(p)t e β(ε(p)−µ) + 1 .
(33) Fig. 13 gives an example for the time-dependence of Π (ret) (t). The calculations have been done at half-filling, β = ∞ for L = 10 4 sites. We note that the susceptibility is in fact peaked at t = 0. For long times, it oscillates with frequency ω Π = 4 and decays as 1/t. This is an important observation as it implies that the limit in Eq. (27) does not exist and that, therefore, the damping constant α is ill-defined.
To analyze the physical origin of the divergent integral, we rewrite the spin susceptibility in Eq. (33) as  (ω) ≡ f (ω)A loc (ω), part of the spectral density, see Eq. (31).
For functions with smooth ω dependence, the Fourier transform generically drops to zero exponentially fast if fact that the LLG approach, by construction, recovers the correct final state where the spin is parallel to the field. In fact, quantitative deviations are found during the relaxation process. The LLG approach is based on first-order perturbation theory in J and on the additional assumption that the classical spin is slow. To pinpoint the source of the deviations, we have numerically solved the integro-differential equation that is obtained in firstorder-in-J perturbation theory and compared with the full hybrid dynamics. The deviations of the perturbative approach from the exact dynamics are found to gradually increase with the propagation time (until the proximity to the final state enforces the correct long-time asymptotics). This is the expected result as the dimensionless small parameter is Jt. However, with increasing J the time scale on which perturbation theory is reliable decreases much stronger than 1/J due to a strong enhancement of retardation effects which make the perturbation more effective and produce a stronger torque.
Generally, the perturbation can be rather ineffective in the sense that it produces a torque ∝ S(t) × S(t ) which is very weak if the process is nearly adiabatic. This explains that first-order perturbation theory and the LLG equation is applicable at all for couplings of the order of hopping J ∼ T . For the present study this can also be seen as a fortunate circumstance since the regime of very weak couplings J T is not accessible numerically. In this case the spin-reversal time scale gets so large that the propagation of excitations in the conduction-electron subsystem would by affected by backscattering from the edges of the system which necessarily must be assumed as finite for the numerical treatment.
For the one-dimensional lattice studied here, a direct comparison between LLG equation and the exact quantum-classical theory is not meaningful as the damping constant α is ill-defined in this case. We could argue that the problem results from the strength of the van Hove singularities in the conduction-electron density of states which dictates the long-time behavior of the memory kernel of the integro-differential equation which is given by the equilibrium spin susceptibility. As the type of the van Hove singularity is characteristic for all systems of a given dimension, we can generally conclude that the LLG approach reduces to a purely phenomenological scheme in the one-dimensional case. However, it is an open question, which will be interesting to tackle in the future, if this conclusion is still valid for systems where the Coulomb interaction among the conduction electrons is taken into account additionally.
There are more interesting lines of research which are based on the present work and could be pursued in the future. Those include systems with more than a single spin where, e.g., the effects of a time-dependent and retarded RKKY interaction can be studied additionally. We are also working on a tractable extension of the theory to account for longitudinal fluctuations of the spins to include time-dependent Kondo screening, and the competition with RKKY coupling, on a time-dependent mean-field level. Finally, lattice rather than impurity variants of the quantum-classical hybrid model are highly interesting to address the time-dependent phase transitions.