Dynamics of spin relaxation in nonequilibrium magnetic nanojunctions

We investigate nonequilibrium phenomena in magnetic nano-junctions using a numerical approach that combines classical spin dynamics with the hierarchical equations of motion technique for quantum dynamics of conduction electrons. Our focus lies on the spin dynamics, where we observe non-monotonic behavior in the spin relaxation rates as a function of the coupling strength between the localized spin and conduction electrons. Notably, we identify a distinct maximum at intermediate coupling strength, which we attribute to a competition that involves the increasing influence of the coupling between the classical spin and electrons, as well as the influence of decreasing local density of states at the Fermi level. Furthermore, we demonstrate that the spin dynamics of a large open system can be accurately simulated by a short chain coupled to semi-infinite metallic leads. In the case of a magnetic junction subjected to an external DC voltage, we observe resonant features in the spin relaxation, reflecting the electronic spectrum of the system. The precession of classical spin gives rise to additional side energies in the electronic spectrum, which in turn leads to a broadened range of enhanced damping in the voltage.

Both fully quantum-mechanical and fully classical approaches possess strengths as well as limitations.For example, quantum methods are typically restricted to analyzing small systems or short time scales.On the other hand, the LLG equation offers a wider range of applicability.However, when using the LLG equation to describe realistic setups, it often necessitates the inclusion of phenomenological terms and the need to fit its parameters to experimental results or extract them from ab initio calculations [85][86][87].
In this work, we utilize a variation of these techniques, namely a quantum-classical equations of motion (QC-EOM) approach for open quantum systems [15,102].QC-EOM is a QC method which combines the equations of motion for the classical spins with the hierarchical equations of motion technique for the conduction electrons [103,104].Its advantage is that, in the case of noninteracting fermions and when focusing on singleparticle properties [105], the hierarchy terminates exactly at the second tier, respectively, first tier in the case of wide-band limit [106][107][108][109][110][111].‡ The method is therefore numerically exact even far away from equilibrium, allows one to reach long simulation times, and does not contain additional approximations or phenomenological terms that could distort analysis of the relaxation processes.
We use QC-EOM to study the spin dynamics, in particular relaxation processes, of a single classical spin embedded in a chain of conduction electrons controlled by an external voltage bias.We show that the effective damping and, therefore, also the relaxation rates of the classical spin are nonmonotonic functions of spin-electron coupling and that they are strongly affected by external voltage with clear resonant features reflecting the static as well as dynamic electronic spectrum of the system.We analyze and discuss some important differences between our results and the results obtained for an equivalent setup in the previous work [6].Furthermore, we demonstrate that the transient dynamics of a long isolated central system can be modeled by a short open one, i.e., a short chain of several sites coupled to reservoirs.We also discuss the usability of the simpler LLG approach in different regimes of the system.This paper is organized as follows.In section 2 we define the model of a quantum-classical magnetic nano-junction.We describe the QC-EOM formalism for open quantum-classical systems in section 3. The results for the spin dynamics are presented in section 4. First, we investigate in section 4.1 the spin relaxation dynamics in an isolated system, described by a single-impurity Kondo chain with a classical spin.We then in section 4.2 show that the dynamics of a long isolated system can be simulated by a short system coupled to semi-infinite leads.In section 4.3 we introduce a finite voltage drop to the system and investigate the current-driven spin dynamics of the classical spins.We then further elaborate on some of the findings in section 4.4 where we focus on a system weakly coupled to the leads which models a simplified magnetic nanojunction.Section 5 summarizes our findings.Some more technical aspects are outlined in the appendices.

Model of Hybrid Magnetic Nano-Junction
We consider a magnetic nano-junction consisting of a central hybrid quantum-classical system [6,94] coupled to reservoirs of non-interacting electrons.The central part is modeled by a one-dimensional electronic tight-binding chain, in which electrons interact locally with a classical spin positioned at the center of the chain (see figure 1).The Hamiltonian of the quantum sector comprises three contributions and reads where the central electronic chain is described by H C , the fermionic reservoirs ℓ = {L, R} by H ℓ , and the coupling between the chain and reservoir ℓ by H Cℓ .The tight-binding ‡ Note, that the expectation values of many-particle operators might require a higher truncation tier as discussed in detail Refs.[105,112] Figure 1.Schematic setup of a quantum-classical magnetic nano-junction.The central system is modeled as a one-dimensional electronic tight-binding chain with a classical spin, and it is coupled to reservoirs L and R of non-interacting electrons.
chain coupled to a classical spin is modeled by the (classical) single-impurity Kondo model (with s-d-type interaction) The first term in equation ( 2) is the kinetic energy associated with electron hopping between neighboring sites.Here, c † jσ and c jσ denote the fermionic creation and annihilation operators with spin σ = {↑, ↓} at site j and γ is the hopping integral.In the analysis below, we fix γ = 1 meaning that we are using γ as the energy unit and time is in units of γ −1 .All relevant physical constants are absorbed into the model parameters.
The second term controls the electronic filling via electrochemical potential µ.Here, we assume that the system is small enough that its equilibrium electrochemical potential can be set externally, e.g., by an auxiliary gate electrode, which does not contribute to the charge and spin transport.If not stated otherwise, µ is taken to zero which sets the half-filling condition.
The last term describes the local coupling of a localized spin S m (t) ∈ R 3 to the spin polarization of the electrons at the chain site m, where σ denotes the vector of Pauli matrices.Note that H C (t) acquires a time dependence through S m (t), as discussed in detail below.
The electronic reservoirs and their coupling to the central chain is modeled by respectively.Here, d † ℓασ and d ℓασ denote the fermionic creation and annihilation operators of the reservoirs.Further, ε ℓ ασ is the single-particle energy of state α with spin σ in reservoir ℓ, and γ ℓ αjσ is the coupling constant, connecting site j at the edge of the central chain to state α in the reservoir ℓ (see figure 1).This constraint implies that the system-reservoir coupling matrix is nonzero only at the interface sites.For a linear chain, these are j = 1 for ℓ = L and j = N for ℓ = R. Here, the spin-flip processes between the reservoirs and the central system are forbidden γ ℓ αjσσ ′ = 0 (where σ ̸ = σ ′ ).Therefore, the coupling is diagonal in the basis of σ z .The reservoirs are assumed to be in chemical and thermal equilibrium, with chemical potential µ ℓ and inverse temperature β ℓ .

Quantum-Classical Equations of Motion
We employ the QC-EOM technique, which consists of a set of quantum equations of motion for the conduction electrons coupled with classical equations of motion for the classical spin.The time-evolution of the quantum sector of the magnetic junction is governed by the Liouville-von Neumann equation for the electronic density matrix When considering an isolated system, that is in the absence of fermionic reservoirs, H(t) in equation ( 5) is identical to H C (t) from equation (2).
In the presence of fermionic reservoirs, H(t) is identified with H O (t) from equation (1).A system of equations of motion for the reduced density matrix ρ C (t) = tr L+R {ρ(t)} is obtained by tracing out the reservoir degrees of freedom from ρ(t) (where tr X {•} denotes the partial trace over the X sub-system).Therefore, the dynamics of the magnetic junction is reduced to the central part only.To describe the dynamics of the magnetic junction, we use a hierarchical equations of motion approach [102][103][104].For the case of non-interacting fermions and single-particle properties studied in this paper, the hierarchy of equations of motion for the auxiliary density matrices terminates at the second tier exactly [105][106][107][108][109][110][111][112].The equation of motion for the zeroth-tier auxiliary density matrix ρ C , that is the reduced density matrix of the central system contains a unitary time-evolution under the (time-dependent) Hamiltonian H C (t) of the central system.The second term on the right hand side of equation ( 6) generates dissipation, a non-unitary time-evolution due to the coupling of the central system to the fermionic reservoirs.The dissipation operator is defined in terms of first-tier auxiliary density matrices Π ℓ (called current matrices for brevity), which can be expressed via time-dependent nonequilibrium Green's functions with lesser and greater Green's functions G ≶ and Σ ≶ ℓ being the lesser and greater tunneling self-energies due to presence of reservoir ℓ [67,110,111,113].For further details, see Appendix A. We use the wide-band approximation, where the line-width functions are energy independent and read Given the one-dimensional geometry employed throughout this work, the coupling is finite at the interface sites j = 1, N only.The previously stated constraints on γ ℓ αjσ due to the considered geometry and spin conservation at the interface lead to an analogous matrix-structure of Γ ℓ jk (ε) as that of γ ℓ αjσ .
By utilizing the Padé-decomposition of the reservoir Fermi-Dirac distribution [114] where χ ± pℓ (t) = µ ℓ (t) ± iξ p β −1 and η p are Padé coefficients of the N p poles of Padéexpansion, the current matrices take the following form with Padé-resolved auxiliary current matrices Π ℓ,p following the equation of motion [108,110,111] Here, Γ = ℓ Γ ℓ denotes the total line-width function.From the current matrices Π ℓ , the charge I ℓ and spin currents Q α ℓ flowing through the interface between the reservoir ℓ and the system can be obtained Note, that neither the Padé expansion [108] nor quantum-classical methods in general are in anyway limited to the here used wide-band limit [15,37,51,92].We use it for convenience, in particular to reduce the parameter space and the number of tiers in the hierarchy.In particular, the hierarchy for single-particle operators terminates at the first tier exactly within the wide-band limit [109,115].The wide-band limit also allows for a more straightforward analysis of the results.The classical sector contains a single localized spin S m (t), and its dynamics is generated by the classical Hamiltonian The classical spin couples to the expectation value of the local conduction electron spin polarization defined in terms of the reduced density matrix ρ m = tr Λ⧹m {ρ} at site m, obtained by tracing out all chain degrees of freedom Λ except site m.Furthermore, we assume an external magnetic field B(t) acting on the classical spin only, which gives rise to a Zeeman contribution in equation (14).
Using the extension of classical Poisson-brackets to spin systems [116], the classical spin equation of motion is Herein, the local effective field B eff is obtained from the gradient ∇ Sm of the classical Hamiltonian H with respect to the classical spin.In the following, we assume |S m (t)| = S = const.with S = 1.
The QC-EOM method thus consists of solving the coupled set of equations of motion equations ( 6), ( 10) and (11) in the quantum sector and simultaneously the classical equation of motion (16), coupled by the s-d term in equation ( 2) and the spin polarization expectation value (15).We note that the quantum-classical approach employed here is an Ehrenfest-type method [94][95][96].The Ehrenfest approach has been used to study nuclear dynamics in quantum transport in Refs.[117,118], and, in particular, has been combined with the hierarchical equations of motion approach in Ref. [115].

Results
Employing QC-EOM, we investigate the spin relaxation dynamics.We start our analysis with long isolated central system and then compare the dynamics to results obtained for the open system, i.e., short chain coupled to reservoirs.After demonstrating that shorttime dynamics of long isolated systems and an open system are equivalent, we investigate the influence of external driving on the dynamics of the central spin.By addressing a junction with single and five chain points, we discuss the role of the electron states and their polarization on the relaxation.

Isolated Central System
We first analyze an isolated tight-binding chain with a single classical spin adsorbed at its center.An analogous system was addressed by Sayad and Potthoff [6] who investigated the relaxation dynamics by focusing on the switching time of the classical spin after reversing the external magnetic field.In contrast, our emphasis lies in extracting the effective (time-independent) Gilbert damping α.The switching and relaxation time is proportional to α/(1 + α 2 ) [6,119], allowing us to reconstruct its profile from the fitted damping coefficient.Fitting the latter has two main advantages.First, it can be extracted at much shorter simulation times, often way before the classical spin can be considered numerically relaxed.Consequently, we can analyze shorter chains in our study.Second, with the introduction of a suitable fitting model, determining the effective Gilbert damping does not require an arbitrary criterion for identifying the precise moment when the spin is considered fully relaxed in a numerical simulation.This aspect proves particularly significant in regimes characterized by long-lived nutations.
In order to analyze the dependence of spin relaxation on spin-electron coupling J sd , we apply a two-stage switching protocol.In the first stage, the localized spin is set to S m0 = e x (e.g., prepared by a strong external field pointing to the x-direction) and electrons are in the ground state with density matrix Here, U describes the unitary transformation to the eigenstates of the initial quantum Hamiltonian [equation ( 2), at t = 0] and j, k are indexing system sites as well as spinprojections.We consider the electronic system at half-filling set by electrochemical potential µ = 0.The second stage is initiated at time t = t 0 = 0 by a sudden switch of the external magnetic field B → B ′ = Be z , which drives the classical spin out of its initial orientation towards a new steady state."In the following analysis we assume  For further analysis, we also consider the approximation of these results using the LLG equation [120] for the single spin in a constant magnetic field b pointing to the z direction (its components are (0, 0, b)) which has, for our initial conditions, the exact solution: where both the effective (time independent) Gilbert damping α and the effective magnetic field b are fitting parameters.Note that b cannot be identified with the real magnetic field B z .The reason is that the dominant effective precession frequency ω p = b/(1 + α 2 ) of the full system depends on J sd due to the geometrical torque as was already discussed in Ref. [38].
We have performed a nonlinear least-squares analysis of S x (and independently S y ) for various J sd .Representative fits (red lines) for J sd = 1.8 and J sd = 5 are shown in figure 2. Figure 3 shows the fitted effective Gilbert damping α as red empty circles in panel (a) and the fitted effective magnetic field b (black line) as well as the relaxation rate From figure 3(a) it is clear that the spin-electron coupling significantly influences the damping and, therefore, also the relaxation times.The overall shape of the α-J sd dependence, which contains a broad maximum at J sd ≈ 4, can be qualitatively understood following the approximate formula for Gilbert damping § which was derived several times in the literature by different methods, see, e.g., Refs.[6,[121][122][123][124][125].Here D ε F is the local density of the states at the Fermi level (ε F = 0) evaluated at site m.This quantity can be obtained by taking advantage of the Green's function formalism [67,126].Assuming a fixed classical spin, we can calculate D ε F as the density of states of a non-interacting quantum dot in a magnetic field J sd S/2 coupled to two semi-infinite chains of non-interacting electrons.Here, the hopping between the dot and the chain is γ = 1 and chains can be fully represented by their edge (surface) density of states per spin [126,127] By using the spin-resolved retarded G r (ε) and advanced G a (ε) Green functions of the coupled dot, we can express the on-dot density of states at the Fermi level as For further details see Refs.[99][100][101]127].As it is shown in figure 3(a), D ε F is a decreasing monotonic function of J sd (solid green line).Its fast decrease reflects the fact that the classical spin acts as an external magnetic field, which splits the, otherwise degenerated, energy level of the dot symmetrically around ε = 0 by ±J sd /2.The finite D ε F results from the overlap of the these splitted and broadened levels.The broadening results from the coupling to the rest of the chain acting as two semi-infinite leads.
Because the overlap quickly decays with increasing J sd so does the density of the states at the Fermi level.For example, if we would assume a simple Lorentzian broadening functions (i.e., constant DOS e ) the D ε F and, therefore, also α F would for large J sd decrease as ∼ J −2 sd .Consequently, the resulting α F dependence, plotted by dashed blue line in figure 3(a), contains a broad maximum for J sd ≈ 4, because the increase of J sd cannot compensate the decrease of D ε F above this value.In addition, because the relaxation rate is linear in b which rises from b/B = 1 at J sd ≪ 1 to b/B ≈ 2 for J sd ≫ 1, the maximum of T −1 is shifted to higher J sd ≈ 5 when compared to the position of the maximum of effective Gilbert damping (J sd ≈ 4).
Similar behavior was also observed for the switching time in Ref. [6].However, there are some differences.First, the authors of the cited work stated that their switching time does not scale as 1/J 2 sd , as expected for weak J sd .This was attributed to finite-size effects (backscattering) and related numerical instabilities due to very long spin-reversal times for small J sd presented in their work.In contrast, because relaxation, respectively, switching times, are proportional to (1 + α 2 )/α we can conclude that in our case the relaxation times scale with 1/J 2 sd for J sd ≪ 1 where D ε F is approximately constant.This confirms the stability of our method in the weak J sd regime.Second, in Ref. [6] the extremal point of the switching time was placed at J sd ≈ 30 which is much higher than our result J sd ≈ 5.This difference cannot be accounted for by the difference in the used magnitude of the classical spin (|S| = 1/2 in Ref. [6] and |S| = 1 here) or external magnetic field.The latter does not have a significant effect on the position of this maximum, as we discuss in Appendix C.
A plausible explanation for this shift of the maximum in Ref. [6] is the influence of high-frequency oscillations imposed on top of the dominant precession, e.g., nutations [28,36,128].These higher-order oscillations emerge in the dynamics for intermediate and strong J sd ≳ 4 [14,36], as discussed in Appendix B, and are longlived [14].Therefore, depending on the criterion chosen for the fully relaxed (switched) classical spin in a numerical simulation, they can significantly influence the extracted switching time for strong J sd .Our fitting model (19) does not take into account these higher-order terms.A clear difference between the fit and the full dynamics at long simulation times can be seen already for J sd = 5 in figure 3(b).The spin dynamics of the fit (red line) is quickly damped, but small steady oscillations (of different frequency than ω p ) can be seen in the full dynamics (black line) for a very long time in figure 2(c).Therefore, the results in figure 3 for J sd > 4 should be understood as the analysis of the relaxation of the dominant spin dynamics.Nevertheless, this seems to be the correct approach, as it was already shown in Ref. [14] that the quantum-classical models highly overestimate the relaxation time of high-order spin processes, such as nutations, when compared with the results of fully quantum methods.

Short chain coupled to leads
Qualitatively, spin relaxation can be understood as a dissipation of the local nonequilibrium electronic spin excitation into the remaining chain in form of spinpolarized electronic wave-packets [6,7].Due to the finite size of the system, these traveling spin wave-packets reflect at the boundaries and after time τ = N/(2γ) interact with the classical spin, leading to recurrences.In the above analysis we have always used long enough chains to make the results free of any finite-size effects.A counterexample demonstrating recurrences in the spin dynamics for different chain lengths is shown in figure 4(a).
Finite-size effects can be mitigated even for short chains by coupling these to a reservoir [39,129].In our case we couple the system to semi-infinite fermionic leads.approximate the broadening function of the semi-infinite tight-binding chain with γ = 1 around the Fermi level, for which, using the semicircle DOS in equation ( 23), we get Γ ℓ = 2πDOS e (ε = 0) = 2 (i.e., the flat DOS e is set to 1/π).Therefore, spin recurrences are strongly suppressed as only a small fraction of electrons is reflected at the systemreservoir interface.Also note, that we use a so-called partition-free method [109,130], meaning that the leads had been coupled to the system in a very distant past.The nonequilibrium situation at time zero is introduced only by the switching of the external magnetic field.The comparison in figure 4(b) confirms that the dynamics in a long isolated system can be to a large extent reproduced by a much shorter magnetic junction.We also show this in figure 3(a), where the black circles represent effective Gilbert damping obtained from fitting the dynamics of a central system with only N = 11 points coupled to leads.The results are in very good agreement with the data for long chains.In addition we show in Section Appendix F, where we compare the dynamics of quantum and classical spins, that this technique can be used also to address dynamics of a large quantum spin coupled to a long chain.

Influence of Nonequilibrium Electron Transport on Spin Dynamics
The introduction of reservoirs allows us, besides the possibility to mitigate finite-size effects, to investigate the influence of nonequilibrium electron transport, driven by a finite voltage bias between the leads, on the classical spin dynamics.
We first consider a relatively long central system analogous to the above case, where we set N = 21 and Γ ℓ = 2. Afterwards, to clarify the role of particular system states on the spin relaxation, we address short systems N = 1 and N = 5 with a much smaller broadening Γ ℓ = 0.1.In all scenarios, we initialize S 0 = e x and employ the partitionfree method [109,130] where the finite voltage drop was introduced in a distant past.In practice, this means that we evolve the system with a fixed initial classical spin orientation until the steady state is reached.Only then, at t = t 0 , the external field B = Be z , perpendicular to the initial orientation of the spin, is switched on.As in the above analysis of the closed system, we use equations ( 20)-( 21) to extract the effective damping and effective field, respectively, the relaxation rate.However, here this analysis has some important limitations.A typical comparison of the classical spin evolution and the corresponding LLG dynamics with fitted α and b is illustrated in figure 5(b,d), where we set J sd = 3, N = 21 and Γ ℓ = 2.In general, the fit of S x (or S y ) spin component dynamics is stable up to V ≈ 4. Below this value, the voltage probes the energy window of the broadened states of the central system.Yet, the comparison of the time evolution of S z from the model equation ( 21), with α and b fitted from the x-component of full evolution, with the full S z dynamics reveals that the simple LLG model is not correctly capturing the short time dynamics for large enough V .For J sd = 3 a clear deviation can already be seen at V ≈ 3.This is related to the fact, that although the frequency of the dominant precession is stable, the real damping is actually time dependent.In addition, the LLG with a constant effective α correctly describes the long-time spin dynamics only if we avoid the regime of long-living nutations, i.e., J sd > 5.
Nevertheless, the fitting procedure can be used to capture the general qualitative behavior of relaxation represented by the effective relaxation rate T −1 for small J sd and V .Figure 5(a) shows dependencies of fitted rates T −1 on the voltage for three values of J sd .The relaxation rate is relatively stable for small voltages (V < 2) with only small deviations from its equilibrium value.It rapidly declines above V ≈ 2 and cannot be reliably extracted for V ≳ 4.This can be partially explained by realizing that damping is determined by the dissipative processes of electrons leaving or entering the central system.These processes are proportional to the systems density of states at the Fermi energy of the reservoirs [122].Therefore, the damping decreases with vanishing DOS.A crucial role for the damping plays the spin polarization of the relevant states and also the precession of the spin.These result in spin-polarized currents Q ℓ and spin-torques Q L + Q R acting on the classical spin.Although the overall charge current increases monotonically with V (not shown here), the spin torque maxima diminish as illustrated in figure 5(c).
The slightly non-monotonic dependence of the relaxation rate on the voltage in figure 5(a) reflects the details of the DOS of the central system.These effects of particular states are to a large extent overshadowed by the natural broadening from the leads.To elevate these details and to study them in a more controlled way, in the next section, we significantly lower the coupling to the leads and also reduce the central chain to one and then five sites.In this regime, the setup can be understood, e.g., as a  simple model of a molecular magnet weakly coupled to metallic leads.

Magnetic junction
Figure 6 shows the time evolution of the classical spin S z (t) (a), the charge current I L (t) (c) and the spin current Q z L (t) (e) measured at the left system-reservoir interface for various DC voltages V in a single-site magnetic junction.Time slices (marked by vertical cuts of respective colors) taken for times t ∈ {25, 100, 250, 500}γ −1 are shown for S z , I L and Q z L in figure 6 (b), (d), and (f), respectively.To emphasize the role of particular states, we consider weak symmetric broadening functions Γ ℓ ≡ Γ L = Γ R = 0.1 and spin-electron coupling J sd = 2.This coupling is sufficient to significantly split the energy eigenvalues (ε ± = ±J sd /2) and simultaneously small enough to suppress the conditions which can therefore be neglected in the analysis of the spin dynamics.The external magnetic field was set to B = 0.1.
The results in figure 6 (a),(b) reveal that spin relaxation is enhanced when the chemical potential of the reservoir ℓ matches one of the energy levels ε ± of the central system.As before, this can be partially attributed to the elevated local density of states at these energies.The enhancement is caused by resonant tunneling of electrons between (e) the reservoirs and the central energy level.To be more specific, the spin relaxation in this voltage regime is a two-step process: First, precession of the classical spin in the x−y plane induces electronic spin excitations.Second, these excitations are transmitted into the reservoirs via polarized electron hopping, as seen, for example, from the enhanced transient spin current Q z L at V = 2, shown in figure 6 (e),(f).Equivalent relaxation mechanisms have already been described in a number of studies addressing different structures, e.g., thin metal films, and are known as spin-pumping [60,[131][132][133][134].
On the other hand, spin relaxation is strongly suppressed for |µ ℓ | ≪ |ε ± | and even more so for |µ ℓ | ≫ |ε ± |.Because of the weak coupling to the leads, there is only a small charge current for small voltages, i.e., in the |µ ℓ | ≪ |ε ± | regime.Consequently, the dynamic excitations generated by the classical spin in the electronic part are mostly localized in the central system as the tunneling into the reservoirs is limited.On the contrary, for |µ ℓ | ≫ |ε ± | both levels contribute to electronic transport [figure 6(c)] and we see large charge transport.However, both levels are populated almost equally P + (t) ≈ P − (t) which leads to a vanishing spin polarization s z (t) ∼ P + (t) − P − (t) ≈ 0 and consequently to a strong suppression of spin relaxation.This is clear from the vanishing transient spin current Q z L in this voltage-regime as shown in figure 6(e),(f).Nevertheless, the aforementioned qualitative analysis neglects crucial mechanisms.It relies on the assumption that the system consists only of static energy levels ε ± , which holds true only for small magnetic fields B. However, for stronger magnetic fields, one cannot neglect the fact that spin dynamics influences the electronic states [25].Specifically, H C (t) varies with time due to the precession of the classical spin at a frequency ω p , which differs from the Larmor frequency as discussed in Refs.[36,37,135,136] and also in Appendix B. The main resulting effect can be understood by assuming a (nearly) time-periodic classical spin S(t + 2π/ω p ) ≈ S(t) (which is fulfilled in the high-voltage regime).This introduces an essential effect whereby the oscillations generate frequency-dependent side bands, denoted by ε ± = ε ± ∓|ω p | and ε± = ε ± ± |ω p |, to the central levels ε ± [25].Two of these additional channels (in our case ε± ) are available for transport via inelastic tunnel processes.For a more detailed analysis of such processes, we refer the interested reader to the work of Filipović et al. in Ref. [25].It is important to emphasize that in the intermediate voltage regime, where the aforementioned assumption does not hold, a feedback mechanism between the conduction electrons and the classical spin leads to damping of spin dynamics.In this regime, the additional electronic states ε± should be understood as time-dependent.Technically, their time-dependence arises due to the non-perturbative character of our approach.Although these additional channels are of a transient nature, they have a significant influence on the voltage-dependence of spin-relaxation rates.They result in a magnetic field-dependent broadening of the voltage range in which the spin relaxation rates are enhanced.This effect is illustrated in figure 7, where panel (a) depicts the influence of the external magnetic field strength (B = 0.5, 1, 2) on the voltage-dependent relaxation rate T −1 .Consistent with previous settings, we fix J sd = 2 and Γ ℓ = 0.1.To quantify the precession frequency at different voltages, we analyze the spectral density |S x (ω)| 2 obtained through the expression where in practice we, however, use a finite upper limit of the integration time τ set to 100γ −1 , since on this timescale the classical spin is typically fully relaxed.This is exemplified in figure 7(b) for B = 2. Further details can be found in the discussion in Appendix B.Moreover, we determine the dominant precession frequency ω p by fitting it using equation (20) as shown in figure 7(c).The extracted frequency aligns well with the frequency obtained at the maximum of The vertical green dashed line in figure 7 signals the voltage at which the chemical potential of the left lead is aligned with ε + .The other three vertical lines indicate positions of 2ε + .The results shown in figure 7 (a) affirm an enlargement of the voltage window in which spin relaxation is enhanced due to the side bands.As shown by the vertical dashed lines in figure 7, the precession frequency ω B p = 1.2B induced by the external field (for J sd = 3, obtained from the analysis shown in Appendix C) provides a quantitative bound on this voltage-window as 2ε + ≤ V ≤ 2(ε + + ω B p ) (shown in figure 7 by dashed lines, colored according to the respective B). Figure 7 (c) shows that as the voltage is increased and reaches the regime |µ ℓ | ≈ |ε ± |, the dominant precession frequency gets strongly shifted towards the Larmor precession frequency ω L but stays well above it until |µ ℓ | ≲ 2ε + .Finally, in the high-voltage regime |µ ℓ | ≫ |ε ± | the precession frequency is identical to ω L .This result further accentuates the previously mentioned decoupling of electrons from the classical spin in the large voltage regime.
Similar observations hold for central systems consisting of multiple sites with the classical spin at the center of the chain.Figure 8 shows the differential conductance dI L /dV (a) and classical spin z-projection S z (b) both at time t = 200γ −1 and the time evolution of S z (c), all three as functions of voltage drop V for N = 5, Γ ℓ = 0.1, J sd = 5 and B = 0.1.We use these directly obtained quantities to illustrate the resonance features instead of fitted rates T −1 as the fitting with model equations ( 20), (21) proved to be unstable.Panels (b) and (c) in figure 8 reveal that the relaxation is maximal when the chemical potential of one of the reservoirs is in resonance with some of the 2N energy levels ε n of H C [plotted as color bars in the bottom part of panels figure 8(a)-(c)].The analysis of a multi-site central system also shows that not all states contribute equally to spin relaxation.In particular, despite high conductance, the degenerate states (ε n = ε m for n ̸ = m, black bars) contribute only weakly to spin relaxation.This can be deduced from the local minima in S z for V = ±2 although dI/ dV is maximal there.The degeneracy of these states leads to a vanishing local spin polarization.Therefore, these states do not couple to the classical spin and hence leave spin dynamics unaffected.Moreover, we observe that longer chains exhibit significantly enhanced suppression of relaxation between resonant features when compared to single-site systems.This suppression becomes apparent at energy values close to zero (|V | ≈ 0) or at energy values approximately equal to ±1 (|V | ≈ 2).The underlying cause can be attributed to a substantial decrease in the local density of states at corresponding energies, which stems from the rapid decline in broadening as the distance from the leads increases, particularly for strong J sd coupling [99,137].Additionally, as discussed in Section Appendix D, the localization of electrons also contributes to these processes.

Summary
We have investigated the relaxation dynamics of a hybrid quantum-classical spin system using the QC-EOM method.In a first step, we have analyzed spin dynamics in an isolated hybrid system, where we have extracted the relaxation rates and Gilbert damping by fitting numerical solutions to the classical LLG dynamics.The results for the damping as well as relaxation rates of the classical spin show a pronounced maximum at finite electron-spin coupling, which can be qualitatively understood by investigating the local density of states on the Fermi level.We compare our results with the prior work conducted by Sayad and Potthoff [6].We show that our results follow the quadratic scaling of the relaxation rate with J sd as predicted by approximate analytical solutions.This confirms the stability of our method even in the regime of small relaxation times, which proved to be difficult in the previous work.We also argue that the discrepancy in the position of its maxima is due to the differences in the extraction of the relaxation time between the two works.We have also shown that the dynamics of a closed system with a very long chain can be reproduced by a short system coupled to infinite leads with a flat band.
Extending the study to a magnetic nanojunction coupled to fermionic reservoirs and subjected to an external dc voltage reveals a strong influence of the bias voltage on spin relaxation.We have observed that a clear signature of the electronic spectrum is imprinted in the voltage-dependent spin dynamics in an non-trivial way where the spin relaxation is boosted only by the magnetically polarized single-particle eigenstates.The dynamical side bands induced to the electronic spectrum due to the precession boost the damping in a wide range of the voltage drop.
We note that despite the fact that it is in good quantitative agreement with the full quantum approach if the spin is large, as briefly discussed in Sec Appendix F, the used quantum-classical mode has its limitations [138].The classical representation of the localized spin cannot account for some inherently quantum phenomena.These include the Kondo effect [6,139], damping of nutations [14] or possible torques caused by the many-body character of quantum states [62,140].However, the formulation of the QC-EOM as a hierarchical master equation invites its future generalizations to a fully quantum system by going beyond the second tier in the expansion [141,142].
where G < (t, τ ) is the lesser nonequilibrium Green function, G r (t, τ ) is the retarded nonequilibrium Green function and Σ a/< ℓ (τ, t) is the advanced/lesser self-energy due to the coupling to the leads.Following Croy and Saalmann [108] this expression can be rewritten using the general relation valid for both Green functions and self-energies to the current matrix form as stated in equation (7) where The related equation of motion, given in equation (11), can be derived using the Kadanoff-Baym relations by using again the expressions in equation (A.2) together with the identity ρ αβ (t) = −iG < αβ (t, t).In our case we further simplify the expression by assuming the wideband limit and time independent coupling of the system to the leads as stated in equation ( 8), meaning that the leads have been always coupled to the central system (partition-free method).In addition, the finite voltage drop is introduced in distant past t v ≪ t 0 and we initially evolve the system with fixed classical spin till the steady state is reached.In principle any other protocol including time-dependent voltage where χ + pℓ (t) = µ ℓ (t) + iξ p β −1 , can be treated by the method.A detailed derivation of such and similar cases (with different Padé coefficients as used in the present work) can be found in Ref. [108].

Appendix B. Influence of s-d Coupling on Spin Precession
Here we address some additional details of the influence of J sd on the spin-dynamics, in particular, on the spin precession frequency ω p .In accordance with the case presented in the main text in Sec.4.1 we set the chain length to N = 151 and B = 1.We show the spectral density |S x (ω)| 2 (a) for various J sd and the dominant frequency as a function of J sd (b) in figure B1.Finite J sd shifts the dominant precession frequency away from the Larmor value ω p = B [panel (b)] which was already discussed by Stahl and Potthoff [36] and others, e.g., [37,135,136].From our analysis, we obtain the shifted precession frequency due to spin-electron coupling as ω B p = 1.2B for J sd = 3 as employed in Sec.4.3.
In addition, strong J sd , e.g.J sd = 15 in panel (a), also induces distinct highfrequency oscillations.The high-frequency oscillations can be attributed to higher-order terms in spin dynamics, e.g., by assuming a Taylor-series expansion and inertia ∼ S × S giving rise to nutation on a short time-scale [28,36,128].The observed fast oscillations, and, hence, the high-frequency peaks in figure B1 (a), can be attributed to nutation or even higher-order effects.Due to the increasing significance of these contributions with increasing J sd , we assess electronic dynamics as the root cause for these fast oscillations [27].statement, we show in figure D1 the site-resolved amplitude |ϕ n | 2 of single-particle eigenstates |ϕ n ⟩ corresponding to eigenenergies in the lower band ε n < 0, in a closed chain with N = 5 with J sd = 5.The state with energy ε = −3.19(figure D1 black line) is clearly localized around site j = 3 and its probability density rapidly decays with increasing distance from this site.All other states, on the other hand, show a less pronounced spatial distribution.Degenerate states (e.g., states corresponding to eigenenergy ε = −1, blue and green line), have vanishing probability density at the site of the classical spin.Thus, these states have a vanishing contribution to nonequilibrium spin dynamics when tuning the reservoir chemical potentials in resonance with the corresponding eigenenergies of these states although the electronic transmission function for these energies is finite.

Appendix E. Derivation of Hybrid Spin Equation of Motion
In this appendix we derive the effective spin equation of motion within the QC-EOM approach.We start with the differentiation of s(t) (equation 15) with respect to time and employ the hierarchical equations of motion (equation 6), resulting in  Note, that due to the integral over the full history of the system, contained in the timedependence of H C (τ ) ≡ H C [S(τ ); τ ], the resulting equation of motion for S takes that of a non-Markovian Master equation.

Appendix F. Classical versus quantum spin
Here we briefly address the difference between the dynamics of a classical spin and that of a quantum spin coupled to a fermionic chain.To this end, we focus on a system previously investigated by Sayad et al. in Ref. [14], namely a classical or quantum spin coupled sideways (at the left edge) to the chain of otherwise non-interacting electrons.This system can be described by the Hamiltonian (2) with m = 1 and where S 1 is either a quantum spin characterized by the quantum number S q or a classical spin with fixed length |S 1 | = S q (S q + 1) ≡ S. We investigated dynamics where the spin was initially prepared in the x direction and the dynamics is triggered at t = 0 by introducing a magnetic field B that points along the z axes with B z = 2.In Fig. F1 we show two examples for S q = 1/2 (S = 3/4) in panels (a), (b); and for a large spin S q = 5 (S = √ 30) in (c), (d).Black circles represent the result for the localized quantum spin obtained by the time-dependent density matrix renormalization group technique (tDMRG) taken graphically from Fig. 1 in Ref. [14].The lines show our results for the classical spin in three different setups.The red lines represent the time evolution of the x (a,c) and z (b,d) components for the closed system where the classical spin is coupled to a tight binding chain of N = 51 sites.The dashed blue line describes the time evolution of a classical spin coupled to chain of N = 11 sites, which is on the right edge coupled to a WBL electronic reservoir with Γ R = 2. Finally, the turquoise line shows the system with N = 3 coupled to the same type of lead as in the previous case.
Fig. F1 shows that the spin dynamics of the quantum spin with a large quantum number, represented here by S q = 5, is approachable by that of a classical spin of length S = S q (S q + 1).Moreover, one can address this dynamics even with a short chain if coupled to an electronic reservoir in a way that mitigates the finite-size effects.On the other hand, as expected, the agreement between the QC and quantum approach decreases with the decreasing spin-quantum number.The S q = 1/2 case is in this respect the most extreme case, as it represents in a sense the "most quantum" case.The quantum case shows stronger damping at S q = 1/2 than the quantum-classical one.This is related to a more pronounced polarization of the spin density of the electrons at the site 1, as discussed in more detail in Ref. [14].Here, the dots represent tDMRG results for a quantum spin with S q = 1/2 coupled to a chain with 50 sites taken graphically from Ref. [14].The lines show our QC-EOM results for S = S q (S q + 1) = 3

Figure 3 .
Figure 3. (a) Left axis: Comparison of the effective Gilbert damping α fitted to the dynamics of long closed system (red empty circles) and the short open system using Γ ℓ = 2(black circles) with approximation α F from equation (22) (dashed blue line).Right axis: Local (on-dot) density of states at the Fermi level from equation (24) (green solid line).(b) Left axis: Relaxation rate T −1 (blue line) as functions of J sd .Right axis: The fitted effective field b (black line).

Figure 5 .
Figure 5. (a) Extracted relaxation rate T −1 = αb/(1 + α 2 ) as a function of voltage drop for N = 21, Γ ℓ = 2, β ℓ = 40 and three values of J sd .(b),(d) Examples of S x and S z dynamics for small and large voltages together with their LLG fits (red dashed lines).(c) Sums of spin currents which equal the spin torque experienced by the central spin.

Figure 6 .
Figure 6.Time-evolution of S z (t, V ) (a) and S z (V ) at times t ∈ {25, 100, 250, 500}γ −1 (b), charge current I L (t, V ) (c) and I L (V ) (d), and spin current Q z L (t, V ) (e) and Q z L (V ) (f) measured at the left system-reservoir interface.The vertical lines in the first-row panels signal and color-code the time positions of V dependencies in the second row.All panels show data for model parameters Γ ℓ = 0.1, J sd = 2.0, β ℓ = 40 and B = 0.1.

Figure 8 .
Figure 8. Charge-current conductance dI/dV calculated at time t = 200γ −1 (a), snapshot of S z at time t = 200γ −1 (b) and the time-evolution of S z (c), all three as functions of voltage V .Parameters are N = 5, Γ ℓ = 0.1, J sd = 5, β ℓ = 40, and B = 0.1.Eigenenergies of an isolated system are depicted by vertical bars, colored according to their averaged spin polarization ⟨s z n ⟩ = k ⟨ϕ n (r k )|σ z |ϕ n (r k )⟩ of eigenstate |ϕ n (r k )⟩ over all sites k: Red for ⟨s z n ⟩ = + 1 2 , blue for ⟨s z n ⟩ = − 1 2 , and degenerate states are depicted in black in panels (a),(b) and white in (c).

Figure B1 .
Figure B1.Spectral density of the classical spin |S x (ω)| 2 obtained from the Fouriertransform of S x (t) for selected J sd (a) and therefrom obtained dominant precession frequency ω p for various J sd (b).

Figure C1 .
Figure C1.Effective Gilbert damping α (a) and relaxation rates T −1 (b), both as functions of J sd and for different field strengths B in a classical single-impurity Kondo chain with N = 151.The full line in (b) shows the relaxation rate obtained from S z , whereas the dashed line the one obtained from S x .Variance due to the fits in (b) is smaller than the line widths.

Figure D1 .
Figure D1.Space-resolved single-particle eigenstates |ϕ n ⟩ of Hamiltonian given in equation 2 with energy ε n with same parameters as in figure 8, N = 5, J sd = 5 and µ = 0, for a classical spin S = e z positioned at the center of the chain.

5 Figure F1 .
Figure F1.Comparison of the transient spin dynamics of a single spin coupled sideways (left edge) to a fermionic chain with J sd = 1.The spin was initially prepared in the x direction, and at time zero a magnetic field B z = 1 in the z direction was switched on.(a), (b) Time evolution of the S x and S z components of the spin.Here, the dots represent tDMRG results for a quantum spin with S q = 1/2 coupled to a chain with 50 sites taken graphically from Ref.[14].The lines show our QC-EOM results for S = S q (S q + 1) = 3/4.Here the red solid line represents a closed system with 51 sites, the dashed blue lines open system with N = 11 sites and the turquoise dashed line represents an open system with three sites.In open systems, the chains are coupled to a lead, treated in the WBL with Γ R = 2, at the opposite end of the chain to where the spin is located.Panels (c) a (d) show the same as (a) and (b) but for S q = 5 and S = 5(5 + 1) = √ 30.

/ 4 .
Figure F1.Comparison of the transient spin dynamics of a single spin coupled sideways (left edge) to a fermionic chain with J sd = 1.The spin was initially prepared in the x direction, and at time zero a magnetic field B z = 1 in the z direction was switched on.(a), (b) Time evolution of the S x and S z components of the spin.Here, the dots represent tDMRG results for a quantum spin with S q = 1/2 coupled to a chain with 50 sites taken graphically from Ref.[14].The lines show our QC-EOM results for S = S q (S q + 1) = 3/4.Here the red solid line represents a closed system with 51 sites, the dashed blue lines open system with N = 11 sites and the turquoise dashed line represents an open system with three sites.In open systems, the chains are coupled to a lead, treated in the WBL with Γ R = 2, at the opposite end of the chain to where the spin is located.Panels (c) a (d) show the same as (a) and (b) but for S q = 5 and S = 5(5 + 1) = √ 30.