Nonequilibrium quantum dynamics of the magnetic Anderson model

We study the non-equilibrium dynamics of a spinful single-orbital quantum dot with an incorporated quantum mechanical spin-1/2 magnetic impurity. Due to the spin degeneracy, double occupancy is allowed, and Coulomb interaction together with the exchange coupling of the magnetic impurity influence the dynamics. By extending the iterative summation of real-time path integrals (ISPI) to this coupled system, we monitor the time-dependent non-equilibrium current and the impurity spin polarization to determine features of the time-dependent non-equilibrium dynamics. We particulary focus on the deep quantum regime, where all time and energy scales are of the same order of magnitude and no small parameter is available. We observe a significant influence of the non-equilibrium decay of the impurity spin polarization both in the presence and in the absence of Coulomb interaction. The exponential relaxation is faster for larger bias voltages, electron-impurity interactions and temperatures. We show that the exact relaxation rate deviates from the corresponding perturbative result. In addition, we study in detail the impurity's back action on the charge current and find a reduction of the stationary current for increasing coupling to the impurity. Moreover, our approach allows us to systematically distinguish mean-field Coulomb and impurity effects from the influence of quantum fluctuations and flip-flop scattering, respectively. In fact, we find a local maximum of the current for a finite Coulomb interaction due to the presence of the impurity.


Introduction
Diluted magnetic semiconductors [1,2,3,4] are an important class of materials for the spintronics community since they combine properties of (ferro-)magnets and semiconductors. Thus they allow for the all-electrical control of the magnetic degrees of freedom of spintronic devices. The magnetisation of semiconducting devices is mainly caused by the interaction of magnetic impurities in the sample with itinerant charge carriers. To understand the microscopic details of these magnetic materials, reduced model systems are required in which the individual constituents are well under control. This allows to study fundamental questions concerning the interplay of coherent quantum dynamics and dissipation under nonequilibrium conditions. An ideal candidate for such a model system is a small quantum dot which connects metallic leads and also carries a magnetic degree of freedom described by a quantum mechanical spin (magnetic Anderson model).
Magnetic quantum dots have been studied experimentally in ensembles which are particularly suited for the investigation by laser and electromagnetic fields [5,6,7,8,9,10,11]. They are designed with standard lithographic methods and are technologically well established. Moreover, embedding individual Mn ions into quantum dots and studying the electrical properties is possible. [12,13,14]. Small quantum dots with few charge carriers and a single magnetic impurity may become important candidates for efficient high density spintronic devices. Although Mn ions in quantum dots have usually spin larger than S = 1/2, there are good reasons to study the low spin case first. From a methodological point of view, this simple model is an ideal basis to develop the numerical tools necessary to treat the real-time dynamics of more general coupled electron-impurity spin systems as well. In addition, the magnetic Anderson model serves as a generic model to study electronic transport through magnetic atoms or molecules placed between the tip of a scanning tunnelling microscope and a substrate [15,16,17]. Likewise, our model is a phenomenological basis for microscopic studies of molecular junctions based on organic radicals [18]. We mention the switching of the spin state of the central iron ion in a single molecular complex in a double layer on gold by a low temperature scanning tunnelling microscope [19] for which our model also is applicable. Finally, it also mimics features of the dynamics of electrons in quantum dots that are subject to hyperfine interactions with the nuclei of the host material. In that respect, the understanding of the electron-impurity coupling and its influence on the dephasing and relaxation times of qubits is essential for the experiments realized in Refs. [20,21].
Nonequilibrium quantum transport in strongly correlated systems continues to remain a challenging problem. Especially, reliable theoretical treatments prove to be difficult when no small parameter is present in the system, i.e., in the deep quantum regime. Approximate methods are often based on advanced perturbative treatments such as quantum master equations [22,23,24], which base on Markovian approximations and weak tunnelling coupling. The renormalisation group (RG) approach [25,26,27,28,29,30] as well as the functional RG method are able to capture essential nonequilibrium features [31,32,33,34]. Real time density matrix RG's [35,36,37,38] require to represent the system plus the macroscopic reservoirs by a large but finite lattice. Due to the possible appearance of finite size effects, the maximal propagation times are limited.
Quantum Monte-Carlo (QMC) concepts are a priori numerically exact and have been adopted to nonequilibrium quantum transport [39,40,41,42,43]. As opposed to the fermionic sign problem, which shows up in equilibrium simulations of quantum many-body systems, the real-time (nonequilibrium) QMC weight function is itself highly oscillatory and causes the dynamical sign problem. Reaching the stationary nonequilibrium state at asymptotic times remains notoriously difficult [44,39,40,41,42,43]. An analytic continuation to the "complex voltage plane" via imaginary "Matsubara voltages" attempts to circumvent the dynamical sign problem. However, the postprocessed back transformation to real times is plagued by numerical instabilities [45,46]. Recently, long-time and steady state results have been obtained by a QMC sampling of the diagrammatic corrections to the non-crossing approximation [47], stating the reduction of the sign problem.
In this work, we adopt the method of the iterative summation of path integrals (ISPI), see Ref. [48] to the magnetic Anderson model. This approach evaluates a real time path integral expression for the Keldysh partition function in a numerically exact manner and is particularly suitable for nonlinear transport. Recently, this scheme has been carefully verified by a comparison to existing approximations in the appropriate parameter regimes for the single-impurity Anderson model [48]. Furthermore, the prediction of a sustained Franck-Condon blockade in the deep quantum regime has been reported based on ISPI data as well [49]. In contrast to the stochastic evaluation of the real time path integral in the QMC approach, the ISPI scheme calculates the real time path-integral deterministically. Hence the sign problem is avoided. The fact that metallic leads at finite temperature suppress all electronic correlations exponentially beyond a finite memory time window is exploited by the ISPI method. This allows for an iterative propagation of an augmented reduced density operator of the system to arbitrary long times. By construction, the technique is limited to finite temperatures and/or finite bias voltages and not too large system sizes. Whenever numerical results are converged with respect to the memory time, they are numerically exact. Recently, Segal et al. [50] have provided an alternative formulation of this approach in terms of Feynman-Vernon-like influence functionals.
We theoretically investigate the real time dynamics of a single-level quantum dot containing a magnetic impurity with spin-1/2. Exchange interaction with on-dot electrons results in dissipation, induced by the metallic leads to the localised impurity in the quantum dot. The deep quantum regime, characterised by the absence of a small parameter is in the focus of the present article, i.e., all interaction strengths are of the same order of magnitude. In particular, we are interested in the nonlinear transport regime, where a finite bias voltage is applied and linear response theory is no longer applicable. The real time relaxation of the impurity spin as a function of various system parameters is investigated. Due to the additional degree of freedom of the magnetic impurity, an extension of the ISPI scheme [48] is necessary. Although the inclusion of a magnetic impurity affects the single particle dynamics in the first place, coupling of the impurity to the electronic density on the dot renders this extension nontrivial. Our accurate results show that in the considered cross-over regime, the nonequilibrium charge current significantly influences the quantum relaxation dynamics of the impurity. Likewise, the back action of the impurity dynamics on the nonequilibrium current becomes significant. Most importantly, we clearly show that this crossover regime is not accessible by perturbative means.
The structure of the article is as follows. After introducing the model in Sec. 2, we express the Keldysh partition function as a real time path integral in Sec. 3 and show how to evaluate this path sum by a deterministic iteration scheme. We calculate in Sec. 4 the stationary transport current and the impurity spin dynamics. Some of the results are compared to outcomes of perturbative methods for the appropriate parameter combinations. We analyse the influence of the nonequilibrium current on the impurity relaxation. Furthermore, we present results of the stationary tunnelling current in the deep quantum regime, where rate equation results are not reliable. The dependence of the relaxation rate on various model parameters is presented.

Model system
We extend the single-impurity Anderson model for the electronic degree of freedom of a quantum dot (QD) coupled to two metallic leads (L/R) via tunnelling barriers in order to study magnetic quantum dots, see Fig. 1 for a sketch. We assume equal tunnelling barriers at the left and right side, the generalisation to the asymmetric case is straightforward. The total Hamiltonian is given by H = H dot + H leads + H T where we use units such that = 1. The Hamiltonian of the magnetic dot includes the electronic part H el dot and the part H imp to describe the spatially fixed magnetic impurity as well as the coupling term between electron and impurity denoted as H int .
We write the electronic part of the QD as where the operator d σ annihilates a dot electron with spin σ. The dot is formed by a single spin degenerate level, which can be controlled by a gate electrode that shifts the electrostatic potential of the dot Φ D . Hence, the electronic degree of freedom can assume four values {0, ↑, ↓, d}, indicating whether the dot is empty (0), contains one electron with spin σ =↑, ↓= ±1 and energy ǫ σ = Φ D + σ∆/2, or is in a spin singlet state with double occupation (d). The Coulomb repulsion is modelled via U > 0 when the Figure 1. (a) Two metallic leads are coupled to a quantum dot (QD) via tunnelling barriers. The gate electrode allows to tune the electrostatic potential Φ D and the bias voltage V induces an electron current. The QD has a single orbital electronic level which is empty or occupied either by one electron in spin up or down state (b) or by two electrons in a Coulomb-interacting singlet state (c). Moreover, the QD carries a localised magnetic impurity M with spin 1/2. In (b), a single dot electron interacts with M via an exchange interaction J, while in (c) only Coulomb interaction is present for double occupation.
dot is in the doubly occupied state d. The Zeeman level splitting ∆ might be present due to possible external or internal crystallographic magnetic fields. The magnetic impurity (quantum spin) is included via the Hamiltonian with the Zeeman energy ∆ imp . Spin operators of the impurity are given in terms of the Pauli matrices τ x,y,z with τ ± = τ x ± iτ y . The impurity spin τ and the dot electron spins σ are coupled by the exchange interaction of strength J which is captured by the Hamiltonian While the first term shifts the single-particle energies, the second term induces mutual flips of an electron-and the impurity spin, which we denote as flip-flop processes. The interaction vanishes for double occupation of the dot, then two electrons form a spinsinglet state.
As usual, we model the noninteracting leads by H leads = kσp ǫ k c † kσp c kσp , where c kσp annihilates an electron with spin σ and wave vector k in lead p = L, R = ±1. A bias voltage V is symmetrically applied between the (thermal) leads and shifts their electrochemical potentials such that µ p = peV /2. Finally, the tunnelling Hamiltonian H T = kσp γd † σ c kσp + γ * c † kσp d σ describes the energy-and spin-independent tunnelling of electrons with the amplitude γ. We assume a constant density of states ̺(ǫ F ) around the Fermi energy and work in the wide-band limit. Then, the tunnelling is parametrised by the parameter Γ := 2π |γ| 2 ̺(ǫ F ).

The ISPI approach
In order to determine the nonequilibrium electron current and to study the relaxation properties of the impurity, we generalise the approach of the iterative summation of path integrals (ISPI) of Ref. [48]. This is a nontrivial step and requires to incorporate the impurity quantum dynamics consistently within the quantum many-body formalism.

Path integral, generating function and short-time propagator
The ISPI approach is based on the Keldysh partition function The time ordering operator is T K and ρ(−∞) is the density matrix of the system's initial state [51,52]. A generating function Z[η(t)] is constructed to calculate the expectation value of an Operator O(t) via is represented by a (fermionic) path integral. Throughout this work, we assume an initially (t i = −∞) empty quantum dot prepared in the spin-up impurity state. The full Keldysh time interval ∆ t := t f − t i = Nδ t is discretised into time steps δ t . A short-time propagator U δt is introduced such that it is related to the exact propagator U(t + δ t , t) = U δt + O(δ 2 t ). Subsequently, a complete basis of fermionic coherent states is inserted between every two short-time propagators, accounting for the Coulomb and the flip-flop interactions.
To construct a particular U δt , we separate the total Hamiltonian H into the diagonal part H 0 and a nondiagonal part H 1 with respect to appearing dot operators d σ and d † σ . This yields H = H 0 + H 1 with In the interaction picture, the full real-time propagator from the initial to the final time can be written as where the Keldysh contour is divided into N −1 pieces of free propagation by U 0 (t k+1 , t k ) that are connected by N − 2 interaction vertices −iH 1 dt k acting during the transition time dt k . Here, we define t 1 := t i and t N := t f . When t k+1 = t k + δ t and in the limit of very small δ t , the system can either propagate freely during δ t via the free propagator U 0 (δ t ) = U 0 (t k+1 , t k ) = U 0 (t k+1 − t k ) or undergo at most one transition induced by H 1 within the interval 0 < t ′ < δ t . Hence, the short time propagator takes the form where the commutator [U 0 (t), H 1 ] = O(t) is used to bring H 1 to the right of U 0 (t). The remaining time integral is evaluated exactly up to order O(δ 2 t ). A comparable error is obtained by normal ordering, denoted by colons in Eq. (8).

Discrete Hubbard-Stratonovich Transformation
Next, we treat the Coulomb term H U dot in a similar way as in Ref. [48]. Since [H 0 , H U dot ] = 0, the propagator for the system with Coulomb interaction factorises into a free part U 0 (δ t ) and the interacting part exp{−iH U dot δ t }. We apply the Hubbard-Stratonovich (HS) transformation and obtain with the HS parameter λ δt defined via The variable ζ = ±1 is interpreted as a fluctuating Ising-like spin field. Note that the solution given in Eq. (10) is uniquely defined as long as 0 ≤ Uδ t ≤ π. By this step, the interacting problem with local-in-time Coulomb repulsions is effectively mapped to a noninteracting problem, at the price that the appearing spin fields interact over a finite time range with each other. In general, the condition δ t Γ ≪ 1 is needed to minimise the systematic error of order δ 2 t . If δ t is bounded from below, however, U should be adjusted in agreement with the condition in Eq. (10). The exponential in Eq. (9) commutes with H 0 and we can write the full short-time propagator as exp{−i(H 0 + H U dot )δ t } = 1/2 ζ=±1 exp{−i H ζ 0 δ t }, thereby absorbing the classical part of the Coulomb interaction into the single-particle dot energies according to In passing, we note that due to the imaginary energy component, H ζ 0 should not be considered as a Hamiltonian. Instead, we obtain a short-time propagator by enforcing normal ordering, i.e., Combining the short-time propagators again into the full path integral, the path sum extends over all tuples {ζ} := (ζ N , . . . , ζ 1 ) of the HS fields.

The remaining interaction terms
The tunnelling term H T is quadratic in the number of fermions but contains both dot and lead operators. Therefore, the stationary state of the isolated system is in general not an eigenstate of the system with tunnelling. However, for arbitrary electronic coherent and impurity states |Ψ τ ≡ |Ψ |τ , the matrix elements can be written as Hence, the short time propagator is obtained by adding the coupling term H T to H ζ 0 in Eq. (12). To be specific, we introduce the total electronic coherent state as where d σ and c kσp are Grassmann numbers for dot and lead electrons. Flip-flops of the electron-and impurity spin are mediated by the non-diagonal exchange coupling term H ⊥ int in Eq. (6) . According to Eq. (8) it enters the short time propagator as In fact, this structure of the short-time propagator motivates our choice of a "mixed" basis of electron coherent states |Ψ and impurity states |τ . Most importantly, a straightforward derivation of a path integral in this basis becomes possible as the paths are separated into parts that contribute to the matrix element Ψ τ ′ |U δt |Ψ τ either for aligned (∝ ½) or opposite (∝ H ⊥ int ) impurity spin orientations. This form is particularly useful with respect to a numerical summation over discrete impurity paths. An equivalent short-time propagator in form of a single exponential : exp{−iH ζ δ t }: with H ζ := H ζ 0 + H T + H ⊥ int is much more cumbersome.
In contrast, the short-time propagator ½ − iH ζ δ t , though exact up to O(δ 2 t ), is not convenient to construct a path integral. Consider, the case of opposite impurity spin states τ = τ ′ . The phase accumulated by the free propagation between two instantaneous flip-flop events is missing, i.e., Such a propagator would only be correct if the transition process lasted the whole time span δ t instead of being instantaneous. In the resulting path integral, the system does not propagate freely between consecutive flip-flop processes with the resulting continuous limit δ t → 0 being unphysical.

Constructing the full path integral
The path integral for the generating function Z[η] is obtained by using Eq. (15) and the corresponding Grassmann fields Ψ as well as the discrete paths {τ } and {ζ}. This yields where the path sums over impurity and HS spin-fields are performed over the 2N-tuples {τ j } = (τ 2N , . . . , τ 1 ) and {ζ j } = (ζ 2N , . . . , ζ 1 ) with τ j , ζ j = ±1. Within an impurity path {τ }, m flip-flop transitions occur on the Keldysh contour, where ℓ of them lie on the lower branch. We discuss the building blocks of Eq. (17) in the following. In Eq. (17), the action S = S imp + S el + S O [η] is given by the sum of the free impurity action and the electronic action. Due to possible source terms, depending on the observable of interest O(t), S O [η] is added when necessary (see below). The electronic action contains all coupling terms of the electrons to the impurity, to the leads (via the tunnelling) and to the HS fields. In particular, we have the action of the free impurity with the discrete and continuous representation given as the first and the second expression, respectively. This equation also illustrates how a well-defined continuous limit of a discrete path can be obtained (cf. Fig. 2). The electronic action S el can be represented as with the inverse electronic Keldysh Green's function (G el ) −1 t,t ′ naturally given in terms of the action S el with the three contributions We note that one of the time integrations in Eq.
Since the bias voltage enters through the respective lead equilibrium density matrix, see Eq. (23), electronic energies are the same in both leads. We have used the effective dot energy in Eq. (20) defined as on the forward branch, while E − σ (t) = E + σ (t) * on the backward branch.   . The polynomial written in (b) follows upon using Eq. (22). We note that the electron's spin flips in the opposite direction as compared to the impurity (not shown).
Note that a flip index on the backward path is labelled according to the smaller step index of the flipping spins corresponding to the later time. The impurity polynomial can be expressed in terms of the Grassmann fields as Figure 2 illustrates an example of an impurity path. The fields of the forward and backward parts of the continuous inverse Green's function (G el ) −1 are not coupled since the corresponding Hamiltonian is diagonal. However, the associated discrete version connects fields from both branches via the upper right element (1, 2N), which is due to the system's initial state ρ(t i ), see Ref. [52] for details. Throughout this work, we assume factorising initial conditions where When constructing the action S O [η] for an observable O(t m ), evaluated at the measurement time t m , we replace every instance of an impurity-or electron operator in the observable by the corresponding spin-and Grassmann fields from the forward branch at the time step closest to t m . This choice is arbitrary and a replacement on the backward branch would not change physical results. Let us assume that t m − t i = (k m − 1)δ t . The electron operators are replaced by their Grassmann field with time index k m . Correspondingly, we substitute the Pauli matrices according to Since the matrix elements τ ′ |τ x,y |τ are non-zero only for τ = τ ′ , the Pauli matrices should be replaced by field expressions that include neighbouring spins. In other words, only if a flip-flop occurs at time t m , the fields τ A generalisation to observables with two and more time parameters, e.g., correlation functions, is possible via higher-order derivatives of the generating function. We calculate the charge current I(t m ) via the source term and the expectation value of the impurity τ z (t m ) from S τz = ητ km .

Tracing Out Electron Degrees of Freedom
Next, we perform traces over the lead degrees of freedom and perform the path integral over all Grassmann fields c kσp , c kσp . Contour time integrations are transferred to their respective real time counterparts. This results in two integrals over real time for the (+) and (−) branch and generates the 2 × 2-matrix structure for the Keldysh Green's function. The resulting generating function is written as a path integral with the effective electronic dot action S el = S dot el + S env with the (effective) environmental action Here, we have introduced the time non-local Keldysh matrix This equation only holds for t−t ′ = 0. The singularity for t = t ′ will be addressed below. The exponential decay of γ for |t−t ′ | → ∞ is the cornerstone of the ISPI method, which allows us to truncate certain long-time correlations (see below). Equation (26) can also be given in terms of an environmental Green's function This expression is well defined for all values of t and t ′ . An explicit expression for (G 0 env ) t,t ′ is given below in Eq. (31). Before we proceed, we address the (t − t ′ ) −1 singularity of (G 0 env ) −1 . In contrast to its inverse, the Green's function G 0 env itself is finite at t = t ′ . Still, matrix elements decay exponentially with growing time differences. The calculation of observables does not suffer from this singularity either, since it cancels in the fraction which is used to numerically evaluate O . Since the divergence does neither depend on the paths of the impurity-and the HS-spins nor on η, we may rescale Z[η] by the singular factor without affecting O (t m ). For clarity, we keep the same notation for the rescaled generating function. The singularity originates from (G 0 env ) −1 . We collect all non-interacting contributions, which do not depend on η into the sum with ω U σ := (ǫ σ + U/2). The Fourier transform of (G el 0,σ ) −1 t,t ′ is obtained as . (32) where we have introduced F (ω) = f L (ω) + f R (ω) as the sum of the two lead Fermi distributions. This 2 × 2-matrix is inverted algebraically and transformed back into time space by complex contour integration [53] or numerical integration. The function Fig. 3. Whenever this divergence arises, we will remove it, see for instance in Eq. (42), by multiplication of −i det G el 0,σ , or equivalently, we replace with the self-energy The particular form of Σ η σ depends on the observable O.
for two combinations of ω U σ and eV for βΓ = 1/5. Except for t = t ′ , the function is smooth for all real times. The imaginary part is discontinuous at t = t ′ , which turns out to be harmless for our numerical scheme. A possible way to cure it, is to introduce a high frequency cutoff exp{− |ω| /ω c } in Eq. (32), see dotted lines in the inset for the same vertical scale for ω c = 100Γ.
The discrete form of the matrix D σ follows via the relation To remove the singularity at k = l, we choose the regularisation Alternatively, introducing a high frequency cut off exp{− |ω| /ω c } in the Green's function [48] yields a consistent result, see the inset of Fig. 3(b). For ω U σ = 3Γ and eV = 4Γ, Eq. (36) yields (G el 0,σ ) l,l ≈ 0.3384i and with the cutoff method with ω c = 100Γ, we obtain 0.3245i (the difference of ∼ 4% decreases for larger ω c ). However, using Eq. (36) has two advantages: first, it does not modify off-diagonal Green's matrix elements, and, second, we do not need an additional parameter ω c .
We emphasize that both methods of regularisation obey the necessary causality relation, This follows from the causality structure of the Green's matrix and the self-energies Σ [52]. Here, the diagonal elements of both G el 0,σ and the Σ matrices have to be understood as the average, i.e., the integral of the time non-local matrix elements in an interval δ t around the point t − t ′ = 0. The discrete version of Σ η σ follows accordingly. For the current, e.g., using Eq. (29) yields We note that also source terms for observables O containing dot fields may be added to the action. The effective full inverse Green's function (G eff ) −1 is given either by Plugging in those pieces, the remaining path integral is recast as a discrete sum, i.e., Hence, we have obtained the Keldysh partition function as a sum over expectation values of the polynomial P of Grassmann numbers in a system with Green's function G eff σ . The next step is to derive an expression for P [{τ }] that refers to G eff σ only by applying Wick's theorem [52]. With Eq. (22), we have For an odd number m of flip-flops the expectation value vanishes, since each process contributes a creator and an annihilator for electrons with opposite spins. Odd m implies an odd number of alternating products of creators and annihilators, the number of d is different from the number of d. When applied to any state |ψ in the trace, this changes the particle count by 1 and the projection with ψ| vanishes. Therefore we have to consider paths with even m only. As an example, we evaluate Eq. (39) for the impurity path shown in Fig. 4 before we turn to the general formalism and arbitrary paths {τ }.
The exemplary path features m = 4 flip-flops and we find S eff and the expectation value of the mixed operator product factorise with respect to the spin degree of freedom. Even in the presence of spin-mixing terms this factorisation remains valid. Applying Wick's theorem to Eq. (40) yields P σ ≡ − det Ξ σ with Figure 5. Dependence of γ ++ on t − t ′ for α = Γ/(2β) and eV β = 4π. We depict the absolute value and the real as well as the imaginary part. For π|∆t|/β 2, the absolute value decays exponentially with a decay time proportional to T −1 . The period of the oscillations in the real and imaginary parts is determined by V .
The general procedure to construct Ξ σ for an arbitrary path {τ } having an even number m of flip-flops is as follows, Using matrix element (Ξ σ ) k,l = −i d q k σ d r l σ = (G eff σ ) q k r l , the final expression for the generating function follows as where the summation over impurity paths is restricted to tuples {τ } with τ 1 = τ 2N = τ i , i.e., correct boundary conditions along the Keldysh contour are imposed. The limit δ t → 0 appears explicitly here, since there is no continuous measure used for the discrete spin paths, neither for the HS-nor for the impurity spins.

Iterative summation of the path integral
The exact generating function in Eq. (42) is intractable due to the exponentially growing size of the matrices for long propagation times and an adept numerical treatment is necessary to proceed. The off-diagonal elements of (G eff ) −1 , given in terms of γ(p, t − t ′ ) in Eq. (27) decay exponentially with increasing distance t−t ′ from the diagonal, at finite T and/or V [48]. This fact is illustrated in Fig. 5. In addition, the bias voltage induces oscillations. The correlation or memory time τ c =: (K −1)δ t determines the range of the lead-induced correlations, they are exponentially suppressed for time differences larger than τ c . This allows us to restrict the path summation to a finite memory time window and thus, the number of paths that contribute to Z[η], originating from the magnetic impurity as well as from the Coulomb interaction, is drastically reduced. We use K as a memory time parameter in the ISPI scheme in an extrapolation procedure for τ c → ∞. The latter eventually gives converged results independent of τ c .
To proceed, we rotate the basis unitarily, thereby rearranging the matrix elements of (G eff σ ) −1 such that increasing distances with respect to the diagonal correspond to increasing time differences. For a given K, we obtain the block band-matrix where N C = N/K and the blocks ✄ ✂ ✁ A σ k,l are K−dimensional matrices, whose entries are given by those of (G eff σ ) −1 taken from the rows and columns in the range of {(j − 1)K + 1, . . . , jK} with j = k, l. Since we neglect exponentially small components, the To illustrate the scheme, a particular impurity path is shown in Fig. 6, which is of length 8δ t , consisting of 2N = 18 vertices having 12 flip-flops. The calculation of the determinant in our scheme needs quadratic block matrices on the diagonal, which we obtain again by a unitary transformation. The off-diagonal blocks are reshaped accordingly. In Fig. 6 this procedure is illustrated for Ξ ↓ with τ c = 2δ t (K = 3), obtained after the rearrangement. The hatched matrix elements are disregarded. The path with N = 9 is divided into N C = 3 segments with K vertices on each branch. In analogy to the blocks and d − ↓,6 from the second segment to the first and third, respectively (blue arrows). Such a reordering is always possible and renders all diagonal blocks of Ξ σ quadratic. We define the recursive notation X = A, B to compactify the computation of the determinant for the blocked Keldysh partition function as The double line denotes matrices which themselves consists of blocks with the subscript D giving their dimension in blocks. The determinant of X is calculated iteratively [48] in D − 1 steps of which each performs the following manipulations: (i) Perform a Gaussian elimination of the block in the second row, first column.
(ii) Neglect in the resulting element of second row, second column products like , which connect segments beyond the nearest neighbour.
(iii) Expand the determinant after the first column, thus reducing the problem by one in block dimension.
While step (i) and (iii) are exact algebraic operations, step (ii) is the second building block of the ISPI method, necessary for the scheme to remain consistent with neglecting correlations beyond τ c . A step k → k + 1 is performed solely based on the determinant after step k and the spin orientations in segments k and k +1. Here we stress that within the time window τ c , the ISPI scheme takes into account important non-Markovian effects in a natural way. In the present work, these correlations are lead-induced. Within a typical Markovian approximation, the real-time dependence of the GF in Fig. 6 is replaced by δ(t − t ′ ). Including terms in the iteration that connect segments beyond nearest neighbouring K-blocks would require information about impurity spins in "earlier" segments < k, which are beyond τ c . After step (iii), we arrive at Next, we construct the Ξ σ -matrices for finite correlation times starting from Eq. (42) by filling the entries with the elements of G eff σ [{τ, ζ}]. The latter depends on the entire spin path {τ, ζ} between t i and t f . In principle, inversion of the full inverse Green's function is possible but out of reach for practical applications since the numerical effort grows exponentially with propagation time. To remain consistent with the finite correlation time approach, we have to find approximations of ✄ ✂ ✁ B σ i,i(±1) that depend on spins of the nearest neighbour segments i(±1). We observe that blocks on the secondary diagonals contribute much less than those on the main diagonal. Hence, we expand (G eff σ ) −1 in powers of the off-diagonal blocks, yielding for all 1 ≤ k, l ≤ N C and |k − l| = 1 (the index σ was omitted). The blocks B -blocks are filled as described above. Next, we use from Eq. (33) the relation G eff σ = D −1 σ G el 0,σ , see Sec. 3.5. The Ξ-matrices are free of any singularity as well and from Eq. (47), we obtain the approximate inverse of D σ . We multiply the result with the free Green's matrix in block form. For step Exemplary paths for the impurity (a) and the HS-field (b). For τ c = 2δ t (K = 3) both are decomposed into segments of length K = 3 in real time, which contain 2K = 6 spins. For example, path segment . Accordingly, we find that Collecting all parts, we can finally express Z[η] iteratively as The real time propagator of the ISPI scheme is introduced as with the chosen initial configuration Furthermore, we use the definition {τ, ζ} j = (τ ∓ jK , . . . , τ ∓ (j−1)K+1 , ζ ∓ jK , . . . , ζ ∓ (j−1)K+1 ) as the tuple of those impurity-and HS-spins that lie in the j-th path segment of length K, see Fig. 7. The propagator Λ j,j−1 (and matrix blocks) depends on all HS-and impurity spins in segments j − 1 and j. The prefactor is related to the number and the position of the flip-flops. m j is the number of flip-flops in segment j, out of which ℓ j lie on the backward branch. The phase is defined as

Extrapolation procedure
By construction, the numerical value of Z[η] contains two systematic errors, namely the finite discretisation time step δ t and the finite correlation time τ c = (K − 1)δ t . In the limit δ t → 0 and K → ∞, however, the iterative procedure yields an exact result. The major benefit of this iterative scheme is that the numerical costs scale linearly with evolution time t m − t i . The systematic errors are eliminated [48] by an extrapolation of the numerical results to vanishing Trotter increment δ t and infinite memory time τ c → ∞.
Due to the Trotter time discretisation, all expressions are by construction exact up to order δ 2 t terms. For a fixed τ c and small enough values of δ t , we extrapolate the numerical values of some observable to δ t → 0 and thus completely eliminate the Trotter error. An example is shown in Fig. 8 (a) and (b) for the current and the impurity orientation, respectively. Depending on the observable, Trotter convergence may be achieved on different scales [57]. Note that one source of errors is the numerical derivative in Eq. (30) which results in tiny imaginary parts of observables, ranging between 10 −5 to 10 −3 . For typical parameters, it is at least one order of magnitude smaller than the numerical error from the linear extrapolation.
For the memory extrapolation τ c → ∞, we do not have a strict mathematical argument at hand, in contrast to the Trotter extrapolation. Whenever results are convergent, however, we find empirically two typical behaviours: (i) either the numerical results for O (τ c ) depend linearly on 1/τ c with small deviations, (ii) or their dependence on 1/τ c is reasonably smooth and exhibits a local extremum. The latter case is consistent with the principle of least dependence, see Refs. [48,49,54,55] for a verification when compared to analytical results. An example of the linear scaling of the numerical results with τ −1 c is illustrated in Fig. 8 (c) and (d). When O (τ −1 c ) shows a weak dependence on τ −1 c in a certain corner of parameter space, we still try to apply criterion (ii). Such a behaviour results from a trade-off between accuracy and computational costs. A minimal Trotter error requires minimal δ t , while, at the same time, a maximally large correlation time τ c is desirable. Naturally, these requirements are limited by the exponentially increasing numerical costs. This is illustrated in Fig. 9.

Restricting the number of flip-flops in the memory window
In order to reduce the exponentially growing number of contributing paths (∼ 4 K ) without affecting the accuracy, we may exploit that F j ∼ (Jδ t /2) m j in Eq. (52). The number m j of flip-flops in path segment j is 0 ≤ m j ≤ 2(K − 1). We observe that the smaller the weight of each segment is, the more flip-flops it contains. On the other hand, the number of path segments {τ } j with m j flip-flops (given by 4C Rapidly decreasing weights of the paths may not be (over-)compensated by increasing weights for 0 ≤ m j ≤ K − 1, since each contribution is small and the number of paths decreases again for larger m j ≥ K. The behaviour of the impurity weights is illustrated as follows. Consider the case when m j is close to the maximum 2(K − 1). Both path classes with m j = 0 and m j = 2(K − 1) contain the same number of elements (four), while each path contribution in the second class is weighted by (Jδ t /2) 2(K−1) . , For typical values of K = 4, δ t Γ = 1/2, and J = Γ, the weight is ∼ 2.5 × 10 −4 . This also holds for all K ≤ m j ≤ 2(K − 1). Since m max j is unknown a priori , we include it into our code as an additional parameter. Then, we perform a numerical estimate by a spot sample of the parameter space. It turns out that for many cases, it is sufficient already to choose m max j = 2. Fig. 10 shows an example where both I and τ z converge quickly for increasing m max j . This drastically reduces the CPU running times, e.g., for parameters chosen as in Fig. 10, from more than one month to typically three to five days.

Charge Current and Impurity Dynamics
The ISPI scheme has originally been developed for the Coulomb-interacting singlelevel quantum dot (SLQD). By reproducing established analytical and experimentally confirmed results, the general validity of the approach has been shown in Refs. [48,49].
Here, we focus on novel transport features caused by the magnetic impurity and its interaction with dot electrons. We emphasize that novel dynamical and transport features are mediated by the transverse or flip-flop interactionĤ ⊥ int , given by Eq. (4). Without the possibility of flip-flops the orientation of the impurity spin and its quantum state could not change and not participate in the dynamics. The remaining longitudinal part of the interactionĤ int causes a renormalisation of rates and energies which adds to the effect stemming from a magnetic field. Necessarily, flip-flop processes are involved from the beginning to investigate the non-trivial impurity dynamics by considering the time dependence of the impurity orientation τ z .
In all presented results below, the impurity is initially polarized and then the coupling to the leads is switched on. We observe that the time-decay of the polarisation is well described by a single exponential with a constant decay rate. In order to single out the relevant physical processes, we compare our numerical ISPI results to a diagrammatic perturbation theory in the weak-to intermediate exchange interaction regime.
We show that for the appropriate parameter regime, the exact numerical results are in accordance with the perturbative result and, by this, we obtain a first intuitive explanation of the impurity dynamics. A next step is the transfer of its plausibility to the ISPI results. However, interaction-induced deviations from the perturbative theory are large enough to clearly illustrate the need for a non-perturbative theoretical description, provided by the ISPI results.

Real-time decay of the impurity polarisation
In Fig. 11 we present the time evolution of the impurity polarisation τ z for different values of the exchange interaction J. The remaining model parameters are chosen as Φ D = ∆ = ∆ imp = U = 0, and βΓ = 1 as well as eV = 0.6Γ. As a function of time, the impurity polarisation shows a decaying behaviour, well described by an exponential relaxation for intermediate to long propagation times. A faster decay of the polarisation is observed as the impurity interacts stronger with the electron spins. In passing, we note that a rate equation ansatz with constant rate τ −1 R , where τ R is the relaxation time, results in an exponential decay as well, i.e., τ z (t) ∝ e −(t−t i )τ −1 R . The parameters are chosen to yield an isotropic (symmetric with respect to [relative] spin orientations) model system. In this case the anti-ferromagnetic interaction favours anti-parallel orientation of electron-and impurity spin. Over long propagation times, the coupling to the unpolarised leads then destroys any polarisation of the impurity. It is therefore reasonable to assume, that the rates for up-and down flips are equal. In the chosen parameter range, the polarisation of the impurity in contact with the leads is well described by a Markovian dynamics, i.e., solely by the time dependent probabilities P τ (t) of finding the impurity in state |τ at time t. Apparently, this simple theoretical prediction agrees well with the numerical results, see Fig. 11. While the impurity interaction energy is comparable to the tunnel coupling and considerably affects the transport behaviour as we show below (see Fig. 13), the rather high temperature and bias voltage nevertheless reduce the relevance of coherent dynamics due to on-dot interactions to a secondary role.
We then turn to the investigation of the inverse relaxation time τ −1 R . In Fig. 12(a), we present results for varying J and U = 0, and three different bias voltages. These show a nearly quadratic behaviour growing from zero (no relaxation) in the sense that for a fit of the results for 0 ≤ J ≤ Γ/2 to a polynomial function aJ b the exponent b lies between ∼ 1.8 and ∼ 1.9. An exact quadratic dependence of τ −1 R on J is obtained only in cases where the dynamics is strongly dominated by sequential (incoherent) flip-flop processes. This is only realized when J ≪ Γ. The corresponding rates may be obtained based on the real-time diagrammatic technique developed by Schoeller and Schön [24] yielding Details of the derivation can be found in Ref. [53]. It reveals the physical structure and allows for the intuitive interpretation of the processes contributing to sequential flipflops. In the numerator of the integrand, we have the sum of all four possible ways to multiply one of the lead's occupations (f + ) with another or the same lead's probability to find an empty state (f − ) at some energy. Each of these four combinations is then multiplied by the Lorentzian spectral density for the two different spin states each shifted by ±J (longitudinal interaction energy). This suggests the following interpretation: A  Fig. 12(b) shows, how the ISPI result (blue crosses) compares to the sequential relaxation time (blue solid line). Although the latter is of the correct order of magnitude, it is systematically larger than the exact value by 10%. Since J is not a small parameter of the system, we can presume that the deviations are mostly coming from coherent higher-order flip-flop processes, which are neglected in Eq. (54). Another source of those deviations may be that free Green's functions are used for the derivation of the rates in Eq. (54). In their qualitative features, however, both results agree. From their finite value at zero bias voltage they grow monotonically. While for small voltages the relaxation shows a nearly quadratic functional form (power-law), for larger bias voltages eV 1.25Γ, it exhibits a linear behaviour.

Impact of the impurity interaction on the current
As opposed to the slow impurity dynamics, measured in terms of Γ −1 , the current is relaxing fast into the stationary state. This behaviour is caused by the strong coupling between the leads and the dot. For the parameters considered here, the upper limit for reaching stationarity is about t ST Γ −1 . Therefore, we consider the stationary current. Fig. 13 depicts the current as a function of J with V = 0.6Γ and βΓ = 1 both for vanishing Coulomb interaction and for U = Γ/2. The current decreases with stronger The similar characteristics of the LB curve and the ISPI data points suggest that the current is mainly affected by the longitudinal part of the electron-impurity interaction. This is probably due to the relatively high temperature and, consequently, a short coherence time, which strongly limits the influence of coherent dynamics. While the LB theory and ISPI agree for J = Γ/8, growing differences for increasing J show the effect of flip-flop processes. The current for small finite Coulomb interaction (no error bars given), though consistently smaller than the LB values and also decreasing with J, drops slower than for vanishing U .
impurity interaction. To distinguish the influence of the longitudinal (single-particle) and transversal (spin-scattering) part of the interaction, we compare the ISPI results with the Landauer-Büttiker (LB) current I LB (see [56]), which can be written here as For J = 0 this is an approximate expression, as it only includes the effect of the longitudinal impurity interaction, which acts as an effective magnetic field in the sense of a mean field. Similar to the sequential relaxation rates of Eq. (54), the current formula has a simple physical interpretation. The joint density of dot-electron states is given by a Breit-Wigner function, whose width equals the tunnel coupling strength and whose resonance lies at the single-electron energy ω σ ± J. Hence, the (non-interacting) current is given by the integral over the energy-dependent difference of the left and right lead's occupation multiplied by the density of available dot states at that same energy. The difference in occupation of the lead electronic states is largest around the Fermi level, where it has the biggest overlap with the density of states for J = 0. With increasing J, the density resonances "move away" from the Fermi level, where f + L (ω)−f + R (ω) decreases and the current drops. This effect is explained in terms of the single-particle energy shift due to the longitudinal component of the impurity interaction only. Fig. 13 shows that the flip-flop termĤ ⊥ int has a considerably smaller influence on the charge current at this rather large temperature (incoherent regime) than the longitudinal part of the interaction. Despite the qualitatively similar behaviour of the LB current and the exact data, the flip-flop scattering causes an additional significant current drop that grows for growing J. A finite Coulomb interaction of U = Γ/2 increases the resistivity of the dot and the ISPI points are consistently lower than the LB values. The decreasing effect of the longitudinal part of the impurity, however, is partially compensated by a broadening of the joint density of states due to Coulomb fluctuations.

Finite Impurity Interaction and Coulomb Repulsion
In the deep quantum regime, where no small parameter exists, ISPI is certainly applicable and able to describe physical properties not predictable by perturbative methods. In this section we study, how the relaxation rate and the current behave as functions of bias voltage, Coulomb interaction and temperature, respectively. Figure 14 presents results in a voltage range 0 ≤ eV ≤ 3Γ, with J = Γ, U = Γ/2 and for temperatures βΓ = 1, 2 and 5. The ISPI data of τ −1 R are indicated by the symbols, while the solid lines mark the corresponding perturbative rates. The latter exhibit the same features (power-law growth, followed by a [quasi-]linear behaviour, which finally saturates) as in Fig. 12(b), which are more pronounced for lower temperatures. As expected, the ISPI data points deviate considerably from this lowest-order approximation. Quantum coherent effects are increasingly relevant since, all energies are of the same order. Both, the degree of quantitative differences and the deviations in the qualitative behaviour increase with lower temperatures.
In Fig. 15, for each of two different temperatures four different current curves are shown-one for each possibility to either have (i) only mean field dynamics (LB), (ii) the full Coulomb interaction without flip-flop processes ("no flips"), (iii) flip-flop dynamics without Coulomb fluctuations ("mean-field U"), or (iv) the fully interacting dot ("full int."). For J = Γ and V = 2Γ, the Coulomb energy is varied between 0 ≤ U ≤ Γ. The  situation "mean-field U" is implemented by setting Φ D = U/2 and the HS-parameter λ = 0 to illustrate the effect of the "classical" or mean-field part of the Coulomb interaction. This is tantamount to setting the third term in Eq. (11) to zero thereby neglecting the Coulomb interaction induced fluctuations, which results in a shift of each single electron energy by U/2. Since this shift tends to move the transport channels away from the Fermi level, i.e., the region with the highest density of state in the leads, the current drops. By its nature, this decrease is equivalent to the one observed for the LB current. Only for the "single-interaction" currents ("no flips"), we show the error bars. The reason why no margin of confidence is given for the fully interacting case, regards the comparability of the error data. Calculating the "full int." current is a very time consuming task and thus, the extrapolation involves considerably fewer data points. Nevertheless, this does not render these values unreliable. We still see a compelling linear behaviour of the 1/τ c extrapolation with errors of the order of 1% based on the sample standard deviation. Notice that with about 10%, the relative error of the current values is rather small, the small variations of the "full int." data are solely due to the extrapolation errors and have no physical meaning.
For both temperatures, both the LB current and the current without Coulomb scattering show only a weak dependence on U due to the single-particle energy shift. The current with full Coulomb interaction but fixed impurity shows a local maximum for U ∼ Γ/2, which is more pronounced for βΓ = 5. In this case, the fixed impurity acts as an effective static magnetic field. The ISPI values for the fully interacting dot vary strongly over the considered U interval, but are scattered around the "no flips" and "no U" curves.
As long as the Coulomb interaction is small, all current values in both figures lie close. The difference between the respective current values is given by the inclusion or exclusion of flip-flop processes only. Hence the rather good agreement of the U = 0 values suggests, that even at this temperature flip-flop processes alone affect the current only weakly. Nevertheless, for lower temperatures, the flip-flop processes start to influence the current more strongly, which results in an increased resistivity. The case of the "no flip" current (fixed impurity) is equivalent to a Coulomb-interacting SLQD in a magnetic field. Both curves in (a) and (b) show a very similar dependence on U, featuring a local maximum at around U = Γ/2. For the lower temperature, the relative height of the broad current peak is twice as big as for βΓ = 1. This effect is also caused by the broadening of the dot's joint density of states due to the Coulomb fluctuations.

Conclusions
We have studied the real-time nonequilibrium dynamics of a single-orbital magnetic quantum dot including Coulomb interactions.
In order to obtain stationary nonequilibrium states at asymptotic times, the ISPI scheme is employed and extended to cases, when an additional magnetic degree of freedom is present. Besides the appearance of a Hubbard-Stratonovich field, which decouples the Anderson repulsion term in the Hamiltonian, we have to include the impurity interaction on the same level. This nontrivial task requires an additional summation over paths of the impurity spin degree of freedom. The resulting action in the path integral formalism involves the Green's function of the quantum dot as well as its inverse. Inversion of the Green's matrix enlarges the numerical effort tremendously. From the technical point of view, appearing matrices, dealing with the impurity dynamics may violate the necessary block diagonal structure. However, a unitary transformation helps to build up proper ingredients for the algorithm. Then, also impurity induced correlations become tractable and do not violate the exponential decay of quantum many-body correlations. We have presented how an efficient truncation scheme provides accurate results for the coupled spin dynamics. Results are given for a quantum spin-1/2 impurity on the dot, whereas the generalization to an impurity with a larger spin is possible. This would be necessary when, e.g. a Mn system is under investigation.
For the same kind of exchange coupling, the implementation of a Mn impurity essentially increases the dimension of the impurity path sum. Instead of summing over all step-like paths of a spin-1/2, it involves the sub-class of step-like paths (steps between orientations differing by one) in the space of the six possible orientations of spin 5/2. Of course, the numerical efforts also increase but still are within reach of the ISPI method. Further work will be dedicated to this goal.
The exponential drop of time correlations due to the leads' coupling allows for an efficient truncation of the appearing sums-the main building block of the ISPI scheme. Its application yields high quality numerical data for the impurity relaxation time and the tunnelling current as a function of the bias voltage in the presence of (Coulomb) and electron-impurity interactions. We have performed necessary checks to compare our findings to established results. In the regime of small impurity interaction, where sequential flip-flops dominate the impurity dynamics, we have found good agreement with a classical rate equation. This is a useful tool to gain insight into the dominating processes in the incoherent regime. Relaxation is described reasonably well by a rate equation when lead-induced coherences are absent.
In the deep quantum regime, however, we find that the ISPI method is the only tool to obtain both the correct order of magnitude and the qualitative features of the relaxation rate as it depends on the system parameters U and J. The same holds for the influence of J and U on the current in the deep quantum regime. Most importantly, the ISPI scheme proves to be useful to cover the full cross-over regime where no small parameter can be identified and thus any perturbative approach becomes invalid.
Furthermore, Kondo physics in such a single spin system under nonequilibrium conditions is of course an interesting subject to study. It emerges when the Coulomb interaction is large compared to the tunnel coupling and the temperature is sufficiently low (also in comparison to the tunnel coupling). In the present work, the two interaction terms in the Hamiltonian and the strong tunnel coupling presently limit the application of ISPI to intermediate temperatures. Therefore, Kondo features have not yet been obtained. In that regard, the further development of the method is still demanding, see also the discussion in Ref. [48]. Nevertheless, due to the suppression of longtime correlations at finite voltages, the regime of nonlinear transport where the Kondo correlations compete with the finite bias is accessible and will be treated in future work.
We have provided a first glimpse on the interesting new physics that comes into reach with the ISPI scheme. A generalisation of the model to several localised magnetic impurities, with electrons mediating a finite magnetisation between them should be possible. The real-time dynamics and all-electrical control of such devices could be simulated. The presented scheme is also applicable to provide the x-and y-components of the impurity spin, thus yielding the complete spin dynamics and the real time dephasing on the Bloch sphere.