Prerelaxation in quantum, classical, and quantum-classical two-impurity models

We numerically study the relaxation dynamics of impurity-host systems, focusing on the presence of long-lived metastable states in the non-equilibrium dynamics after an initial excitation of the impurities. In generic systems, an excited impurity coupled to a large bath at zero temperature is expected to relax and approach its ground state over time. However, certain exceptional cases exhibit metastability, where the system remains in an excited state on timescales largely exceeding the typical relaxation time. We study this phenomenon for three prototypical impurity models: a tight-binding quantum model of independent spinless fermions on a lattice with two stub impurities, a classical-spin Heisenberg model with two weakly coupled classical impurity spins, and a tight-binding quantum model of independent electrons with two classical impurity spins. Through numerical integration of the fundamental equations of motion, we find that all three models exhibit similar qualitative behavior: complete relaxation for nearest-neighbor impurities and incomplete or strongly delayed relaxation for next-nearest-neighbor impurities. The underlying mechanisms leading to this behavior differ between models and include impurity-induced bound states, emergent approximately conserved local observables, and exact cancellation of local and nonlocal dissipation effects.


I. INTRODUCTION
A small system ("impurity") in an excited state and coupled to a large bath ("host") at zero temperature is usually expected to relax over time and to approach its ground state.After the initial excitation of the impurity, the excess energy is dissipated via the impurity-host coupling and through the coupling of the host degrees of freedom into the bulk of the host system (see Fig. 1).For a system with a macroscopically large number of degrees of freedom subject to the principles of thermodynamics, this dissipation process is irreversible.This picture of generic relaxation dynamics explains the interest in exceptional cases, where the system is trapped in a metastable state that does not decay on timescales exceeding by far the typical intrinsic timescales governing the microscopic degrees of freedom.
Incomplete or delayed relaxation and metastability in impurity-host models [1,2] is closely related to incomplete or delayed thermalization of extended lattice models.In both cases, much of the interest in 1.Typical structure of an impurity model.The impurities are coupled to a much larger host system.After an initial local excitation the full system, impurities plus host, is expected to relax locally in its ground state in the vicinity of the impurities after the excitation energy is dissipated into the bulk.metastable states is due to their promise for controlling non-equilibrium dynamics and for related functionalities [3].Compared to notoriously difficult lattice models, models with single or few, initially excited impurities embedded in a large host represent an interesting class of comparatively simple systems that may hold a key to the understanding of metastability.Here, we report on metastable states in the real-time dynamics of three different prototypical system-bath models, an uncorrelated quantum, a classical and a quantum-classical hybrid model.In all three cases, the exact dynamics is numerically accessible on long time scales.
In the recent decades, much progress has been made in understanding the thermalization of generic macroscopically large quantum systems, the main paradigm of quantum-statistical physics and the foundation of thermodynamics [4,5].Here, an important concept is the eigenstate thermalization hypothesis [6][7][8][9][10] as reviewed, e.g., in Ref. [11].One route to non-thermal states of quantum-lattice models is via integrability.For (onedimensional) systems with a large number of conserved local observables, the long-time dynamics may result in a state described by a generalized Gibbs ensemble [12][13][14][15].Another route is provided via disorder, either on the single-particle level or via many-body localization [16,17].This is similar to classical Hamiltonian dynamics [18]: It is known that ergodicity and the equivalence between long-time and ensemble averages of observables can be broken in the case of a large number of integrals of motion.Violations of ergodicity are found for integrable systems [18] but also for systems parametrically close to integrability [19][20][21][22] or in systems with glassy dynamics [23,24].
For impurity-host systems (or open quantum systems) the situation is not very different: After an initial local excitation of the impurity, one generally expects a relaxation of the reduced density matrix of the impurity to its (canonical) thermal state if the impurity-host coupling is weak [2,[30][31][32][33].On the other hand, in the case of band gaps or finite band widths, incomplete relaxation and residual dissipationless dynamics may occur [34,35].Relaxation to non-thermal states may be found in the gapless case for a sufficiently strong impurity-bath coupling [35][36][37].
Recently, a metastable state and incomplete spin relaxation have been observed in a system consisting of a classical impurity spin that is exchange coupled to a spinful Su-Schrieffer-Heeger model at an edge site [38].Here, the topological state of the host and the associated presence or absence of a protected edge mode are found to control the relaxation of the classical spin.This is a prime example for an impurity system with pre-relaxation dynamics, analogous to pre-thermalization in quantum-lattice models.However, the dynamical decoupling of the impurity and the stabilization of the excited state on a long time scale is due to a gapped spectrum for two-particle excitations, which blocks further energy dissipation.A similar effect has been observed for a classical spin locally coupled to a one-dimensional half-filled Hubbard model [39].In this case the Hubbard-U and the narrow spectrum of (quantum) spin excitations control pre-relaxation and metastability.Fast but incomplete relaxation to a metastable intermediate excited state, followed by extremely slow complete relaxation is also known from the decay of a local doublon excitation in the Hubbard model at large U [40,41], or for a magnetic doublon [42] in the strong-J limit of the Kondo lattice.
Here, we study the exact real-time dynamics by numerical integration of the fundamental equations of motion for three different prototypical impurity models, a quantum, a classical, and a quantum-classical hybrid model.All share equivalent geometries, namely a onedimensional lattice model serving as the bath and two additional impurities which are locally coupled to nearestneighbor (n.n.) or to next-nearest-neighbor (n.n.n.) sites of the lattice.Specifically, we study (i) a tight-binding quantum model of independent spinless fermions on a lattice with two stub impurities, (ii) a classical-spin Heisenberg model with two weakly coupled classical impurity spins, and (iii) a tight-binding quantum model of independent electrons with two classical impurity spins.The long-time relaxation dynamics initiated by a local excitation of the impurities can by studied numerically for large lattices in all three cases, and in all cases we find qualitatively very similar results: There is complete relaxation to a time-independent final state in the case of n.n.impurities while there is incomplete relaxation or pre-relaxation for n.n.n.impurities.In all three cases the effect can be understood after a thorough theoreti-cal analysis.However, it turns out that the uncovered mechanisms are very different.The different systems are discussed separately in Secs.II, III, and IV.Our conclusions are summarized in Sec.V.

II. STUB IMPURITY MODEL
We start the discussion by considering a tight-binding model with two stub impurities.The host system is given by non-interacting spinless fermions on a one-dimensional chain of L sites with open boundaries.The nearestneighbor hopping T = 1 sets the energy scale.The two impurities (a and b) are given by two additional sites or orbitals coupling via a hybridization of strength V to the host sites i a and i b .We will consider n.n. or n.n.n.sites i a , i b located at the center of the chain.A sketch of the system is shown in Fig. 2. The Hamiltonian consists of three terms: the Hamiltonian of the host, the impurity sites, with on-site energy ϵ f = 0, and the host-impurity hybridization Here, c i annihilates a fermion at site i, and f a annihilates a fermion at impurity site a.We study the real-time dynamics of the system after a quantum quench of the hybridization from zero to a finite value V .At time t = 0, the state of the host is assumed to be prepared in its non-degenerate ground state with N = L/2 fermions, i.e., half-filling: Sketch of the stub impurity model of spinless fermions.Two fermionic impurity sites or orbitals are coupled via a hybridization V to a one-dimensional lattice with nearest-neighbor hopping T .T = 1 sets the energy scale.The host system is at half-filling.
where k runs over the one-particle eigenstates of H host with one-particle energies ε k < 0. Furthermore, the impurity sites a and b are assumed as fully occupied at time t = 0. Hence, the initial state of the full system is the state As the system is non-interacting, it is sufficient to formulate an equation of motion in terms of the one-particle reduced density matrix ρ.Its elements are defined as The indices I, J run over the L + 2 host and impurity sites: I, J ∈ {1, 2, ..., L, a, b}.Initially, at t = 0, the density matrix ρ(0) has a block-diagonal form with an L × L block representing the host system and two 1 × 1blocks representing the impurities.We are interested in the time evolution of the occupation numbers At t = 0, we have where Θ denotes the Heaviside step function, µ = 0 the chemical potential, and T host the hopping matrix of the host system.Furthermore, ρ aa (0) = ρ bb (0) = 1.
The time dependence of the density matrix ρ(t) is obtained via the von Neumann equation of motion Here, T is the hopping matrix of the full system, Eq. ( 1).The formal solution of Eq. ( 10) is given by where the diagonal matrix of one-particle eigenenergies ε and the unitary matrix U formed by the one-particle eigenstates of T are obtained by solving the eigenvalue problem The t = 0 state, Eq. ( 6), is not an eigenstate of the full Hamiltonian with V > 0. It instead represents a state that is locally excited, in the vicinity of the impurities.One naively expects that the local excess energy and the locally enhanced fermion density at the impurity sites are dissipated to the bulk of the system over time and that the system approaches the fully relaxed state that is locally characterized by the ground-state energy density and the ground-state impurity occupations.
Numerical results for a system with L = 500 host sites and impurities a and b coupling to nearest-neighbor (n.n.) sites i a = 250 and i b = 251 located symmetrically around the chain center are shown in Fig. 3. On a time scale t ∼ 100, the impurity occupations n a (= n b ) relax and approach a value n a ≈ 0.505 close to their ground-state value n (gs) a = 0.5.Similarly the occupations n ia = n i b of the host sites closest to the impurities and also the occupations at more distant sites, e.g., n ia+20 relax.The characteristic time scale for the dissipation of particles (and of energy) can be read off from Fig. 3 by comparing the dynamics of n ia with that of n ia+20 .The total system size (L = 500) is large enough such that reflections of the propagating wave packets at the open system boundaries do not yet interfere with the dynamics in the vicinity of the impurities -on the time interval considered.The time evolution turns out to be completely different, however, when coupling the two impurities to nextnearest-neighbor (n.n.n) host sites.This is demonstrated with Fig. 4.While we do find a fast initial relaxation, the relaxation process stops at t ∼ 2, and the impurity occupations start to oscillate around a value (n a ≈ 0.625) that is considerably larger than the ground-state value.These oscillations are undamped and persistent (until finite-size effects in form of interference with reflections from the boundaries become important).
This qualitative difference between n.n. and n.n.n.impurities, i.e., complete or incomplete relaxation, is likewise observed for impurities at arbitrarily large but odd or even distances, respectively.Key to this effect is the absence or presence of bound single-particle energy eigenstates of the post-quench Hamiltonian.
Two different types of localized eigenstates can be distinguished: (i) For each of the two impurities, there is a pair of bound states that split off from the lower and upper edges of the valence band.These four high-excitation energy bound states are localized near the impurities with a weight that decays exponentially at large distances.The case of strong hybridization V is instructive: For V → ∞, the hopping term, Eq. ( 2), can be ignored, and the Hamiltonian describes a system of two decoupled dimers with two degenerate eigenstates at −V , and two more degenerate states at +V (V > 0).This degeneracy is lifted for finite T ≪ V , and two bonding-antibonding pairs of bound states are formed, one with negative energies below the bottom of the band and one with positive energies.As V decreases, the bound states remain localized and centered around the impurities but their weight is increasingly distributed over the lattice.At V = 1, only a single state from each of the two pairs remains, a spatially symmetric bound state with negative energy and an antisymmetric bound state with positive energy, while the other two states have merged with the bulk continuum.This first type of bound states is generic and thus present for both cases, impurities coupled to n.n. and to n.n.n.host sites.
(ii) A bound state of a different, second type is present in the case of n.n.n.impurities only.It is given by where i a + 1 = i b − 1 denotes the host site between the sites coupled to the impurities.This state has a finite weight on this and on the two impurity sites only, it is "superlocalized".Furthermore, its eigenenergy, ε loc = 0, resides within the continuum of band states.This type of states is well known from flat-band systems [43][44][45][46][47][48]: When coupling a stub impurity to every second host site, the resulting translationally invariant tight-binding lattice model has a unit cell consisting of three sites, and its band structure features a flat band, resulting from superlocalized states, besides two dispersive bands.
To discuss the impact of bound states on the postquench relaxation dynamics, we can straightforwardly adapt some concepts developed in Ref. [10].Accordingly, we decompose the expectation value O(t) ≡ ⟨O(t)⟩ of a local operator O with a Heisenberg time dependence as where the first term is the long-time average, and where the time-dependent second term δO(t) is a fluctuation part with a vanishing long-time average.As a measure for the strength of persistent temporal fluctuations, we consider the long-time average of the absolute square of the fluctuation part: For O = O IJ = c † J c I , see Eq. ( 7), the expectation values O IJ (t) = ρ IJ (t) are given by the elements of the one-particle reduced density matrix.A straightforward computation yields: Here, U Iµ is the I-th component of the µ-th eigen vector of the total hopping matrix, see Eq. ( 12), and is an element of the one-particle reduced density matrix in the basis of the eigenstates of the total hopping matrix.Furthermore, we have assumed that there are no degeneracies and no gap degeneracies in the spectrum of the one-particle eigenenergies.In our specific case, the spectrum is indeed non-degenerate.However, particlehole symmetry implies the presence of gap degeneracies, i.e., This modifies the derivation that led to Eq. ( 17), see also Ref. [10], and we find The contribution to temporal fluctuations, measured with δ 2 ρ IJ in Eq. ( 19), that stems from extended eigenstates must vanish for L → ∞, as the components of the eigenvectors are proportional to L −1/2 and hence On the contrary, the components U iµ of a localized eigenstate µ are independent of L and can be large at some sites i, as compared to the components of delocalized eigenstates, such that their contributions to ρ µν (t = 0) in Eq. ( 19) become significant.19), as a function of the system size L for n.n. and n.n.n.impurities.This is demonstrated with Fig. 5, where the timeaverage of fluctuations of the occupation number is shown for one of the impurity sites.For n.n.impurities, δ 2 na decreases with increasing L and eventually vanishes in the thermodynamical limit L → ∞.Note that the first type of bound states, discussed above under point (i), does not prevent relaxation due to spatial inversion symmetry, as discussed in the Appendix A.
Contrary, in the case of n.n.n.impurities, the strength of the fluctuations is essentially independent of the system size for L ≳ 10 3 .The non-vanishing temporal fluctuations result from a bound state of the second type, see point (ii) above.This explains the observed incomplete relaxation for the n.n.n.case, see Fig. 4.
In the L → ∞ limit, the effect of gap degeneracies is vanishing in the n.n.case.In the n.n.n.case and for large L, the fluctuation δ 2 na , as computed via Eq.( 19), is smaller by 3 • 10 −4 (or by about 20%) compared to the value obtained from Eq. ( 17), i.e., disregarding gap degeneracies.
Note that in the n.n.case and for L = 500, see Fig. 3, the plateau values for the occupation numbers at t ≈ 100, i.e., n a ≈ 0.505 and n ia ≈ 0.525, are close to but different from their L → ∞ expectation values n (gs) a = n (gs) ia = 0.5 in the ground state of the post-quench Hamiltonian.The latter are fixed by particle-hole symmetry.For longer propagation times t max ≫ 100, the plateau in the time evolution (see Fig. 3) is repeatedly interrupted by revivals (not visible on the time scale in Fig. 3).When computing the long-time averages, Eq. ( 15), including the revivals, we find converged t → ∞ values n a ≈ 0.510 and n ia ≈ 0.526.At L = 500, a propagation time t < t max ≈ 0.5 × 10 4 has turned out sufficient.
In fact, the dynamics of a system of non-interacting fermions is constrained by the constants of motion c † µ c µ .Hence, the system will relax to a non-thermal state with long-time averages of n a and n ia equal to the averages in the generalized Gibbs ensemble or, equivalently, in the diagonal ensemble (see Ref. [10]).For an arbitrary oper-ator O, the diagonal average is defined as with C J = ⟨J|Ψ(t = 0)⟩ and with O JJ = ⟨J|O|J⟩.At L = 500, for example, the numerical values n In general, for an integrable system, such as a noninteracting fermion impurity model, the expectation values of local one-particle observables generically do not relax to a thermal state, regardless of the presence or absence of bound states.However, the presence of bound one-particle eigenstates µ of the post-quench Hamiltonian is crucial for the question whether there is any relaxation at all or whether the system is trapped in a metastable state.In the case of n.n.n.impurities, there is a superlocalized bound state of the stub impurity model that prevents relaxation and forces the system into a metastable state without any further dissipation.This explains the qualitatively different relaxation dynamics for n.n. and for n.n.n.impurities.

III. CLASSICAL HEISENBERG IMPURITY MODEL
A similar effect in the relaxation dynamics is found for a purely classical spin model, i.e., for a one-dimensional Heisenberg model of L classical spins s i (i = 1, ..., L) with nearest-neighbor antiferromagnetic exchange coupling J > 0, where in addition two classical impurity spins S m (m = 1, 2) are locally exchange coupled to the host spins at sites i 1 and i 2 .The geometry is the same as for the stub impurity model, see Fig. 6 and compare with Fig. 2. The coupling strength of the local exchange is denoted by K. We assume that K ≪ J and K > 0, i.e., weak antiferromagnetic exchange interaction.The classical Hamilton function is given by where the products of vectors are to be understood as dot products.The length of the host spins and of the impurity spins is set to s ≡ |s i | = 1 2 and S ≡ |S m | = 1 2 , respectively.We consider a lattice with open boundaries.The sites i 1 and i 2 , where the impurity spins are coupled to the host, are assumed to be n.n. or as n.n.n.sites in the center of the lattice.The host nearest-neighbor exchange coupling fixes the energy scale, J = 1, and we assume K = 0.01 unless otherwise stated.
The equations of motion are easily derived within the classical Hamilton formalism by making use of the spin Poisson bracket [49,50].They attain the form of Landau-Lifschitz equations [51].For the impurity spins we have where "×" indicates the cross product, while for the host spins It is immediately apparent that the length of each individual spin represents a constant of motion, such that the spin dynamics is constrained to a configuration space given by the L + 2-fold direct product S ≡ S 2 × • • • S 2 of 2-spheres with radius 1/2.Furthermore, the total energy and, due to the SO(3) spin-rotation symmetry of H, the total spin m S m + i s i are conserved.Moreover, the system has an SO(3)-degenerate ground-state manifold, opposed to non-degenerate singlet state of the quantum-spin model [52].
The spin dynamics is initiated by a parameter quench of the local exchange coupling from zero to a finite value K at time t = 0. We assume that the initial state of the system at t = 0 is given by one of the ground states for K = 0.For the host spins the ground state is given by an antiferromagnetic Néel state with respect to an arbitrary axis.Specifically, we choose The impurity-spin configuration at t = 0 is taken to be non-collinear, i.e., S 1 = Se x and S 2 = Se y .This implies that, after switching on K at t = 0, the spin dynamics is immediately driven by a finite spin torque.Naively, one again expects that the excitation energy stored in the center of the chain is dissipated into the bulk of the system and that locally, in the vicinity of the impurities, the system approaches a ground-state spin configuration after a sufficiently long propagation time, assuming that the host of the system is sufficiently large to avoid unwanted interactions with excitations backscattered from the chain boundaries.
The equations of motion ( 22) and ( 23) are easily solved numerically for systems with L ≃ 10 4 host sites.For this  system size and for J = 1, there are no finite-size effects in the form of reflection of spin excitations from the system boundaries up to a propagation time of t ≃ 10 4 .Fig. (7) displays the time dependence of the x components of two impurity spins and of the cosine of the enclosed angle for the case that the impurity spins are locally coupled to nearest-neighbor (n.n.) host spins at the chain center.We find that after t ≈ 800 the dynamics has stopped and the system has reached one of its local ground states with an antiferromagnetic impurity-spin configuration and with an antiferromagnetic configuration of the host spins (S 1 ↑↓ S 2 ) in the vicinity of the chain center.
For two impurity spins coupled to n.n.n.sites, however, the time evolution is fundamentally different.As can be seen in Fig. 8, the system does not relax to a local ground state, at least not on the numerically accessible time scale.Rather, we find that after a propagation time t ≈ 1000, the system state becomes trapped in a stationary state, in which the impurity spins precess around a common axis.Up to t = 10 4 there is hardly any relaxation to the expected ferromagnetic (S 1 ↑↑ S 2 ) impurity-spin configuration.The angle enclosed by S 1 and S 2 starts to deviate only slightly from its initial zero value (see green curve).We also note that the qualitative difference between the relaxation dynamics for n.n. and for n.n.n. is not due to the larger distance of the impurity spins in the n.n.n.case.While the distance between the impurities does have an effect on the relaxation dynamics since it determines the time it takes for the spins to "communicate" with each other, this distance dependence turns out as negligible compared to the odd-even effect that is seen in calculations with larger inter-impurity distances d.In fact, we find quick relaxation for distances d = 1, 3, 5, ... and trapping in a stationary state for d = 2, 4, ..., very similar to the results shown in Figs.7 and 8 and analogous to the results for the stub impurity model.
The observed odd-even effect is actually related to the different ground-state spin configurations, i.e., an antiferromagnetic and a ferromagnetic impurity-spin configuration for n.n. and for n.n.n.impurities, respectively, and to the small coupling constant K ≪ J.This becomes obvious when looking at the time derivative of the scalar product S 1 S 2 .From Eq. ( 22) we get In the following we argue that, for n.n.n.impurities, the term on the right-hand side is small, as compared to the inverse of the considered propagation time, in contrast to the n.n.case, and that this explains the observed oddeven effect.We consider three lines of argument with increasing level of sophistication.
In the first and crudest approach, we approximate s i2 − s i1 in an ad hoc way by its ground-state value.For n.n.n.impurities, the Néel ground-state value is s i1 = 0, and thus S 1 S 2 is a constant of motion that cannot relax to its ground-state value S (0) 1 S (0) 2 = 1/4 (see Fig. 8).On the other hand, for the n.n.case, the right-hand side is generically finite for the Néel ground state and of the order of K.In fact, there is a nontrivial dynamics of S 1 S 2 , and S 1 S 2 approaches S (0) 1 S (0) 2 = −1/4 on a time scale ∆t ≈ 8 • 10 2 ∼ K −1 (see Fig. 7).
In our second approach we derive an upper bound for |s i1 − s i2 | for the case of n.n.n.impurities.For the initial state considered here, where the host spins form a Néel state, the excitation energy ∆E is solely stored in the interaction term ∝ K, see Eq. (21).It is given by Hence, for an arbitrary initial impurity-spin configuration, the maximum excitation energy is given by ∆E max = K and is realized for a ferromagnetic alignment of S m and s (0) im .Assuming that this excitation energy ∆E is distributed among the bonds between the closest host spins s i1 , s c and s i2 only (c The right-hand side is at a minimum if the central host spin is Using With the above argument, ∆E ≤ ∆E max = K, we find This upper bound is a very conservative estimate as in the course of time the excitation energy will be further dissipated to the bulk of the system, and thus |s i1 − s i2 | will be even smaller.We conclude that for K ≪ J, the small available excitation energy of order K very much restricts the host spin dynamics.Via Eqs. ( 25) and (30) this imply that S 1 S 2 is almost conserved if S 1 and S 2 couple to n.n.n.host spins.
Our third approach is based on a linearization of the equations of motion.We start from Eqs. ( 22) and ( 23 Linearization of Eqs. ( 22) and ( 23) yields: and ṡi = J s where we used to estimate the magnitude of the neglected terms.We see that linearization of the equations of motion is possible although |δS m | = O (1) is not necessarily small.The reasoning is the same that led to Eq. ( 30), i.e., the maximum values for |δS m | and |δs i | are limited by the available initial excitation energy ∆E ≤ K. Impurity spins contribute on the order of K to the total energy, while host spins contribute on the order of J.The estimates (33) are well supported by our numerical results underlying Fig. 8.We proceed by computing the time derivative of S 1 S 2 within the linearized theory.With Eqs. ( 31), ( 32) we find For n.n.n.impurity spins, the first term on the right-hand side vanishes since s i2 .Furthermore, we can write With S i = S (0) + δS i and s m = s m + δs m , exploiting that ground-state spin configurations are collinear, and finally using Eqs.(33), one has This means that, within the linearized theory, (d/dt)(S 1 S 2 ) must be considered as zero and that S 1 S 2 is a constant of motion with a correction of the same order of magnitude as the linearization error only.
In the case of n.n.impurity spins, restarting from Eq. ( 34) with a completely analogous calculation but with an antiferromagnetic ground-state alignment s (0) i.e., there is a nontrivial dynamics on an energy scale that is by an order of magnitude larger than the linearization error, such that, even within the linearized theory, S 1 S 2 cannot be considered as a constant of motion.
We have also studied the dynamics beyond the weakcoupling regime.For K and J of the same order of magnitude, one finds a relaxation of S 1 S 2 already after a very short propagation time of t ≃ 100 for both, the case of n.n. and of n.n.n.impurity spins.
For the weak-coupling regime K ≪ J, we conclude that after an initial local excitation of n.n.n.impurity spins, these show an anomalous relaxation dynamics.There is almost no relaxation of S 1 S 2 , i.e., the enclosed angle is almost a constant of motion, on a time scale of about t ∼ 10 4 .This must be contrasted with the case of n.n.impurity spins, where complete relaxation is reached after a propagation time of t ≃ 800.In contrast to the stub-impurities model discussed above, there is no local symmetry of the (classical) Hamiltonian that would lead to a conserved local observable.(Quasi-)conservation of S 1 S 2 is rather emerging in the course of time.After a certain pre-relaxation process (t ≃ 10 3 ) with a sufficient dissipation of energy and spin, the system state has evolved sufficiently close to one of the ground states locally, i.e., in the vicinity of the impurities, such that the further dynamics is very well captured by linearized equations of motion.Indeed, within the linearlized theory, S 1 S 2 is strictly conserved.Its validity range, however, is not only controlled by the a weak local exchange K ≪ J but also by the propagation time.Residual perturbative deviations from the linear dynamics accumulate over time, such that complete relaxation of the system, also for the n.n.n.case, is expected on a long time scale.In fact, indications for full long-time relaxation are seen at t ≃ 10 4 in Fig. 8.

IV. QUANTUM-CLASSICAL IMPURITY MODEL
In the case of the quantum-classical impurity model, we again find a qualitatively very similar effect in the relaxation dynamics.However, to explain the observed incomplete relaxation, it turns out again that a different methodological approach is necessary.
We consider a spinful single-orbital tight-binding model on a one-dimensional lattice of L sites with hopping between nearest neighbors T , where in addition two classical impurity spins S m (m = 1, 2) are locally exchange coupled to the local electron spins s i = 1 2 σσ ′ c † iσ τ σσ ′ c iσ ′ at sites i = i m of the lattice via an antiferromagnetic exchange interaction K. Here, τ denotes the vector of Pauli matrices, and σ =↑, ↓ is the electron spin projection.A sketch of the system is shown in Fig. 9.The geometry is the same as for the previous models.As for the classical Heisenberg model, we assume that K is weak and can be treated perturbatively.The quantum-classical Hamiltonian is We choose T = 1 to fix the energy (and time) scale and an antiferromagnetic exchange K > 0. As for the stub impurity model, we consider half-filling, i.e., N = L electrons.
The equations of motion [53,54] couple the classical and the quantum sector of the theory.For the classical impurity spins, we obtain Landau-Lifshitz-type equations, similar to Eq. ( 22).Here, ⟨s im ⟩ t is the expectation value in the N -electron state |Ψ(t)⟩ at time t.For a given classical spin configuration at time t, the quantum system is uncorrelated, and hence its real-time dynamics is completely described by the one-particle reduced density matrix ρ(t) with elements Its equation of motion is essentially the same as for the stub-impurity model, see Eq. ( 10), but the hopping matrix T replaced by the time-dependent effective hopping matrix T (eff) (t), which includes the classical impurity spins as time-dependent external parameters: We study the time evolution of the full system starting at t = 0 from an initial state where the two impurity spins are in an excited non-collinear configuration (as in the classical Heisenberg impurities case, S 1 = Se x and S 2 = Se y with S = 1  2 ), while the electron system is in its ground state corresponding to this spin configuration.The excitation energy stored in the vicinity of the impurities is dissipated to the bulk of the electron system on a time scale that, even for K of the order of T , typically exceeds by far the time scale that is numerically accessible when using open boundaries and when reflections of propagating excitations from the boundaries of the system shall be avoided.Since the propagation is essentially ballistic, unphysical reflections from the boundaries that disturb the dynamics near the impurities will occur at time t ∼ L/T , i.e., one would have to work with effective hopping matrices of very high matrix dimension.For this reason and as indicated in Fig. 9, we employ so-called absorbing boundary conditions (ABC), which have been developed and extensively tested previously, see Ref. 55. Apart from the conserving von Neumann term, the resulting equations of motion contain a dissipative term and are given by Here, ρ 0 is the initial ground-state one-particle reduced density matrix, {•, •} denotes the anticommutator, and γ is a diagonal matrix controlling the dissipation rate.It has nonzero entries only for the outermost two "absorbing" sites on both sides of the chain.For the concrete computations, we have fixed the values for γ 1 = γ L and γ 2 = γ L−1 , as in Ref. [55], by comparing the resulting spin dynamics using ABC with the exact spin dynamics, i.e., without the dissipative term in Eq. ( 42).This has been done for a shorter propagation time of t = 5 • 10 2 and a larger system size such that reflections from the boundaries are avoided.We find perfect agreement with γ 1 = 0.230 and γ 2 = 0.115.However, the results are quite insensitive to the precise choice.As compared to the system studied in Ref. [55], the optimal parameters are smaller, because the spin dynamics is much slower.The impurity-spin dynamics for nearest-neighbor spins as obtained by solving the coupled equations of motion (39) and ( 42) is displayed in Fig. 10.At short times t ≲ 10 3 there is a pronounced precession dynamics with a small frequency ω ∼ 0.01.This is explained by the indirect RKKY exchange [56][57][58] mediated by the electron system, which is rather weak, even for an exchange interaction of K = T = 1.On a longer time scale t ∼ 10 4 , the system shows complete relaxation, and the two spins reach their antiferromagnetic ground-state configuration (see green line in the figure).
For impurity spins coupled to n.n.n.sites, see Fig. 11, the same precessional motion is found, but with an even smaller precession frequency.This reflects the smaller RKKY exchange due to the increased distance between the spins.However, the real-time dynamics is qualitatively different, as there is hardly any relaxation to the ferromagnetic ground-state spin configuration visible on the numerically accessible time scale.At time t = 10 4 , the angle enclosed by S 1 and S 2 deviates by less than 1% from its initial value only.We conclude that, as for the other models studied, the system is trapped in an intermediate stationary state and that complete relaxation, if at all, takes place on a still much longer time scale.
For smaller K (not shown here), the time evolution is essentially the same in qualitative terms.The only notable difference is that the dynamics is even slower, i.e. characterized by smaller precession frequencies and longer relaxation times.
None of the explanations for incomplete relaxation used for the previously discussed systems is easily applicable to the quantum-classical model.The linearization of the coupled equations of motion is not helpful to identify possible local conserved observables, and in fact is not informative due to the large number of ∼ 4L 2 of degrees of freedom in the quantum sector, i.e., the density-matrix elements.
With linear-response theory [38,54,[59][60][61], we choose a different approach.Conceptually, this is limited to the weak exchange-coupling regime and directly addresses the dynamics of the classical impurity spins.For weak K, the expectation value of the local spin at site i m can be obtained via the Kubo formula as Here, the response function, the K = 0 retarded magnetic susceptibility of the electron system is isotropic χ mα,m ′ α ′ (t) = δ αα ′ χ mm ′ (t) and independent of the spatial direction α = x, y, z.It is given by where ⟨• • • ⟩ is the K = 0 ground-state expectation value and where η is a positive infinitesimal.A straightforward computation yields Note that the chemical potential µ = 0 at half-filling.
For the evaluation of Eq. ( 45), we consider large systems with up to L = 10 5 sites and choose periodic boundary conditions.Hence, the hopping matrix is diagonalized as T = U εU † , where the unitary matrix U with elements U ik = e ikRi / √ L describes Fourier transformation from lattice sites i to wave "vectors" k in the first Brillouin zone.The entries of the diagonal matrix are given by the tight-binding dispersion for short.This yields: or, after Fourier transformation from time to frequency space, the frequency-dependent susceptibility Note that we have the symmetry χ mm ′ (ω) = χ m ′ m (ω) for the nonlocal elements m ̸ = m ′ , while the local elements are m independent, χ mm (ω) = χ m ′ m ′ (ω), due to translation invariance.The representation ( 47) is well suited to compute the Gilbert damping: and the RKKY indirect magnetic exchange which determine the effective equations of motion for the classical spin dynamics (see Ref. 54): In practice, the results of various calculations for different system sizes L as well as for different η must be extrapolated to obtain physical results in the thermodynamic limit L → ∞ and in the limit η → 0. Here, it is important to take the thermodynamic limit first.This is demonstrated with Fig. 12, where the local, α mm and the nonlocal (n.n.n.) Gilbert damping α mm ′ (m ̸ = m ′ ) are shown as function of η for different L. We start the discussion with the local damping (solid lines).First, we see that for any fixed η ≳ 10 −4 , the values for the local Gilbert damping nicely converge with increasing L. System sizes of about L = 100, 000 are sufficient for numerical convergence unless even smaller values of η are considered.Second, the converged values lim L→∞ α mm become independent of η with decreasing η, once η is sufficiently small.We find a rather precise value lim η→0 lim L→∞ α mm ≈ −0.0398K 2 .Here, we note that taking the limits in the opposite order yields the unphysical result lim L→∞ lim η→0 α mm = 0.This is easily understood: For any finite L, the spectrum of oneparticle energies is gapped.Close to ω = 0, the finite-size gap is δ ≈ 2π/L.This implies that for η ≲ δ ≈ 2π/L, the Gilbert damping must start to deviate from its physical value and approach α mm = 0, as there is no damping in a finite system.For the practical calculations, it has turned out that when fixing the "infinitesimal" at η ≈ 2π/L, it is sufficient to control the convergence with respect to L only.
Our considerations for computing the Gilbert damping likewise apply to the case of n.n.n.spins.There is, however, an important physical result that can be read off from Fig. 12: In the case of n.n.n.spins, the converged value for the nonlocal Gilbert damping (see dashed lines) is exactly the same as the local damping, i.e., lim η→0 lim L→∞ α mm ′ ≈ −0.0398K 2 for both, m = m ′ and m ̸ = m ′ within numerical accuracy.We note that a similar result for the nonlocal Gilbert damping has been found for metallic ferromagnets with quadratic energy-momentum dispersion [62].
The equality between the local and the nonlocal damping has in fact important consequences for the spin dynamics, as can be easily seen when rewriting Eq. ( 50) explicitly for two classical spins but with a single Gilbert damping constant α ≡ α 11 = α 12 = α 12 = α 21 : Note that only the nonlocal RKKY exchange coupling J ≡ J mm ′ = J m ′ m (m ̸ = m ′ ) enters the equations.We immediately see that the total impurity spin S tot = S 1 + S 2 and thus S 1 S 2 are constants of motion, as in the case of the classical Heisenberg impurity model, see Sec.III.This implies that there is no relaxation to the groundstate spin configuration at all.So far we have discussed the case of n.n.n.impurity spins, where as a consequence of α mm = α mm ′ (m ̸ = m ′ ) there is no relaxation to the ground-state spin configuration.While for n.n.impurity spins the local Gilbert damping α mm ≈ −0.398K 2 stays the same, we find, on the other hand, α mm ′ ≈ 0.0021K 2 (m ̸ = m ′ ) for the nonlocal Gilbert damping.The signs are such that a solution of Eq. ( 50) must approach the ground state, i.e., in the n.n.case an antiferromagnetic spin configuration.This is consistent with the computed positive RKKY exchange coupling J 12 ≈ 0.0342K 2 (H RKKY = J 12 S 1 S 2 ) for the n.n.case.On the contrary J 12 ≈ −0.0189K 2 for the n.n.n.case with ferromagnetic ground-state spin configuration.
Returning to the n.n.n.case, the equality of the local and the nonlocal damping can be understood analytically.Using Eqs.(47) and ( 48) one finds Nonzero contributions to the double sum are obtained from wave vectors close to the Fermi points, i.e., for k, k ′ = ±π/2 + O(1/L) only (note that ε(k = ±π/2) = 0 = µ).Contrary, for k, k ′ = ±π/2 + O(1), the imaginary infinitesimal can be disregarded, since we may take η = O(1/L), as argued above, and thus , and the two fractions in Eq. ( 52) cancel exactly in the thermodynamical limit.
It is thus sufficient to analyze the contributions from k = ±π/2 + δk and k ′ = ±π/2 + δk ′ with δk, δk ′ = O(1/L) for L → ∞ and show that these give the same result for i m ′ = i m and for i m ′ = i m + 2. The i m , i m ′ dependence of the damping is due to the weight factor A imi m ′ only.We therefore focus on A imi m ′ .Its imaginary part does not contribute to the double sum in Eq. (52).For the discussion of the real part, we first consider k = π/2 + δk and k ′ = π/2 + δk ′ : and, hence, It is also invalid for n.n.impurities and, more generally, for odd distances between the impurities, because for k = π/2 + δk and k ′ = −π/2 + δk ′ , e.g., we have and for odd distances.
Since our explanation of the incomplete relaxation is based on perturbative-in-K linear-response theory, it is necessary to compare corresponding results with those of the full theory (using absorbing boundary conditions), Eqs.(39) and (42).We choose K = T for this comparison.This provides us with a comparatively fast spin dynamics.Results are displayed in Figs. 13 and 14 for the cases of n.n. and n.n.n.impurity spins.
We find a slight phase offset in the precessional motion for the n.n.n.case (Fig. 14).On the logarithmic time FIG.13.Comparison between the full spin dynamics (solid lines), as obtained from the exact equations of motion ( 39) and ( 42) and absorbing boundary conditions, and linearresponse spin dynamics (dashed lines), as obtained from Eq. ( 50) with numerically determined parameters α11 = α22 = −0.0398,α12 = α21 = 0.0021, and J12 = J21 = 0.0342, see Eqs. ( 48) and ( 49), respectively.Time evolution of the x component of S1 (blue) and cosine of the angle enclosed by S1 and S2 (green) for the case of n.n.impurity spins.All other parameters as in Fig. 10.In particular, K = T = 1.scale, this offset is constant.Furthermore, at late times t ∼ 10 4 , a tiny deviation of the angle enclosed by the two spins from its initial t = 0 value is visible in the results from the full theory, hinting towards complete relaxation on a much longer time scale.This is missing in the linearresponse approach.For the n.n.case, where the spin dynamics is much more complicated, the perturbative method also does an almost perfect job, see Fig. 13.While we observe the same but slightly larger phase shift and a slightly longer relaxation time, all the qualitative features of the spin dynamics are fully captured.
We conclude that linear-response approach itself, i.e., perturbation theory in K is quite reliable even for comparatively strong K ∼ T and errors accumulating up to a time scale t ∼ 10 4 /T do not affect the qualitative trend of the spin dynamics.This also holds for the typical additional approximations that are necessary to arrive at Eqs. ( 48) and (49), i.e., weak retardation effects and time independence of the Gilbert damping, see Refs.[54,60].All in all the numerical results demonstrate that the proposed mechanism based on the analysis of the nonlocal Gilbert-damping term in fact captures the essence of the incomplete relaxation.

V. CONCLUSIONS
Using numerical simulations, we have studied the exact real-time dynamics of three different prototypical onedimensional model systems with two impurities coupled locally to nearest-neighbor or to next-nearest-neighbor sites of the host.In all cases we considered an initial state with a local excitation at or near the impurities.The unifying theme of all three models studied is the conservation (or the approximate conservation) of observables that are localized in the vicinity of the impurities.Furthermore, in all cases, the presence of these (quasi) conserved observables depends on the system geometry, and in all cases this crucial for the relaxation dynamics.
The independent-electron tight-binding quantum model with two stub impurities is conceptually the simplest.Due to the lack of interactions, it is integrable; its real-time dynamics is strongly constrained by a macroscopically large number of conserved observables.This implies that local one-particle observables do not relax to their ground-state values but to a non-thermal GGE-like state respecting the constraints.However, it depends crucially on the geometry, i.e., on the relative position of the impurities, whether or not a complete relaxation to a time-independent state for t → ∞, for an infinitely extended system, actually occurs.In the case of impurities coupled to n.n.n.sites, we find persistent oscillations up to the numerically accessible time scale, i.e., before unwanted reflections from the boundaries set in.
The classical Heisenberg model with two locally exchange coupled classical impurity spins shows very similar behavior.For n.n.impurity positions, there is a fast and complete relaxation to the ground state impurityspin configuration.In contrast, in the n.n.n.case, an (almost) undamped oscillatory spin dynamics is found, when the local exchange coupling K is weak compared to the host exchange J. Complete relaxation to the ground state configuration is not observed on the numerically accessible time scale.However, the numerical data indicate that complete relaxation is possible on a much longer time scale, so that the system actually exhibits pre-relaxation.This is an essential difference from the non-interacting quantum system.
Qualitatively the same results are found for the quantum-classical hybrid model with two classical impurity spins locally exchange coupled to an independentelectron system on the one-dimensional lattice, i.e., fast complete relaxation to the ground-state spin configuration in the n.n.case, while in the n.n.n.case and after a fast pre-relaxation, a metastable intermediate state is formed, in which the impurity spins undergo an (almost) undamped oscillation.This intermediate state is stable up to t ≳ 10 4 in units of the inverse hopping.Again, we assume that the impurity-host coupling, the local exchange K, is sufficiently weak.
An explanation for the observed very different behavior for n.n. vs. n.n.n.geometries, common to all three models, does not seem obvious.In fact, quite different theoretical concepts have been put forward as explanations: The incomplete relaxation of the quantum system with n.n.n.impurities is due to the presence of a superlocalized energy eigenstate bound to the impurities and thus due to a local observable commuting with the Hamiltonian.The superlocalized state is reminiscent of the states forming flat bands in tight-binding models on lattices with characteristic geometries.
The metastability of the classical spin model, on the other hand, could be traced back to an approximately conserved local observable, reminiscent of explanations for the pre-thermalization of interacting lattice models parametrically close to an integrable point.In fact, we had to assume that K ≪ J, which places the model parametrically close to the trivial K = 0 point.Here, the weak-coupling limit has allowed us to linearize the equations of motion and thus to understand the approximate conservation law.This is remarkable because the fluctuations δS m of the impurity spins around their groundstate configuration S (0) m are not at all small, since the metastable state is far from the ground state.Rather, the argument can be based on the fact that the excitation energy is of the order of K ≪ J, and that each impurity spin contributes of the order of K to the total energy, as opposed to the host spins which contribute of the order of J.We note that this reasoning seems possible only for impurity models.
The analysis of the real-time dynamics of the quantumclassical hybrid model is much more complicated, since a simple linearization of the equations of motions is not very feasible and not justified.However, the limit of weak local exchange coupling K ≪ T could be exploited in another way, namely by time-dependent perturbation or linear-response theory.This turns out to be reliable even up to intermediate couplings K ∼ T and propagation times t ≲ 10 4 in units of the inverse hopping parameter.Within the linear-response framework, the stability of persistent oscillations in the spin dynamics in the case of n.n.n.impurities is nicely explained by a perfect cancellation of local with nonlocal Gilbert damping constants.However, the exact dynamics obtained numerically clearly indicates that, beyond the perturbatively accessible time scale, the non-equilibrium steady state is actually metastable and that there is further relaxation on a much longer time scale.
While it seems to make no qualitative difference for the relaxation dynamics if quantum degrees of freedom are replaced by classical ones or vice versa, the geometry is a crucial factor.Common to all three models studied is the bipartite system geometry, i.e., the one-dimensionality of the host lattice with nearest-neighbor couplings between host sites, and the local host-impurity coupling.Admittedly, we expect that next-nearest-neighbor host or nonlocal host-impurity couplings (hopping or spinexchange couplings) break the (meta)stability of the nonequilibrium state in the case of n.n.n.impurities and lead to (faster) complete relaxation.However, parametric proximity to the bipartite geometry, e.g., weak nonlocal host-impurity couplings, should still lead to a significantly different relaxation dynamics between impurities at n.n. and n.n.n.positions.Thus, we believe that our results provide valuable insights for our understanding of metastable states and the control of non-equilibrium dynamics.The study of the relaxation dynamics of impurities coupled to a host system on a two-or threedimensional lattice is one of the promising avenues of further research.
Furthermore, it would be very interesting to study the relaxation dynamics for systems including quantum rather than classical spins or, more generally, for correlated quantum impurity models and to check the robustness of the results against quantum fluctuations.It is quite conceivable that also in such systems the geometry plays a crucial role for the existence of (approximately) conserved local quantities.Of course, it will be technically more challenging to reach the relevant time scales.For one-dimensional systems, however, matrix-productstate approaches [63] seem to be suitable to study relaxation dynamics, see e.g.Ref. [64].
FIG.1.Typical structure of an impurity model.The impurities are coupled to a much larger host system.After an initial local excitation the full system, impurities plus host, is expected to relax locally in its ground state in the vicinity of the impurities after the excitation energy is dissipated into the bulk.

FIG. 4 .
FIG. 4. The same as Fig. 3, but for stub impurities coupled to next-nearest-neighbor (n.n.n.) sites ia = 249 and i b = 251.We also choose L = 499.Therewith, inversion symmetry enforces n b = na and ni b = ni a .
ia ≈ 0.526 perfectly agree with the abovementioned long-time averages.Repeating the computations for larger system sizes (up to L = 5000) and extrapolating to L → ∞ yields slightly smaller values n (D) a ≈ 0.505 and n (D) ia ≈ 0.525.

FIG. 6 .
FIG.6.Sketch of a system consisting of two classical impurity spins (red) locally exchange coupled to a one-dimensional classical Heisenberg model (blue spins) with open boundary conditions.K: weak antiferromagnetic local exchange, J: antiferromagnetic nearest-neighbor exchange interaction of the host spins.

FIG. 7 .
FIG. 7. Time evolution of the x components of two classical impurity spins S1 and S2 and of the cosine of the enclosed angle S1S2/S1S2.Excited state at t = 0: host spins are in an antiferromagnetic Néel state aligned to the z axis, impurity spins point into the x and y direction, i.e., S1 = 1 2 ex and S2 = 1 2 ey.Geometry parameters: L = 10000, i1 = 5000, i2 = 5001.Coupling strengths: J = 1 and K = 0.01.
) and substitute S m = S (0) m + δS m and s i = s state spin orientations while δS m and δs i denote the deviations from the ground state.

FIG. 9 .
FIG.9.Sketch of the quantum-classical hybrid model: Two classical impurity spins (red) are locally exchange coupled to a system of conduction electrons on a one-dimensional lattice (blue).K: antiferromagnetic local exchange, T : nearestneighbor hopping.Green: absorbing boundary conditions (ABC).

FIG. 10 .
FIG. 10.Time evolution of the x components of the two classical impurity spins S1 and S2 (blue and orange) exchange coupled to the conduction-electron system at neighboring sites i1 = 34 and i2 = 35 at the center of a chain with L = 68 sites.Green: cosine of the angle enclosed by S1 and S2.Initial excited state at time t = 0: host electron system in its ground state corresponding to the initial (excited) impurity-spin configuration S1 = 1 2 ex and S2 = 1 2 ey.Further parameters: T = 1, K = 1.Absorbing boundary conditions with nonzero diagonal elements γ2 = γL−1 = 0.115 and γ1 = γL = 0.230 (see text).

FIG. 12 .
FIG. 12.Local (m = m ′ , solid lines) and nonlocal nextnearest-neighbor Gilbert damping (m ̸ = m ′ , dashed lines) α/K 2 as function of η for different system sizes L as indicated.Results for large systems with periodic boundary conditions, T = 1.