Full density matrix dynamics for large quantum systems: Interactions, Decoherence and Inelastic effects

We develop analytical tools and numerical methods for time evolving the total density matrix of the finite-size Anderson model. The model is composed of two finite metal grains, each prepared in canonical states of differing chemical potential and connected through a single electronic level (quantum dot or impurity). Coulomb interactions are either excluded all together, or allowed on the dot only. We extend this basic model to emulate decoherring and inelastic scattering processes for the dot electrons with the probe technique. Three methods, originally developed to treat impurity dynamics, are augmented to yield global system dynamics: the quantum Langevin equation method, the well known fermionic trace formula, and an iterative path integral approach. The latter accommodates interactions on the dot in a numerically exact fashion. We apply the developed techniques to two open topics in nonequilibrium many-body physics: (i) We explore the role of many-body electron-electron repulsion effects on the dynamics of the system. Results, obtained using exact path integral simulations, are compared to mean-field quantum Langevin equation predictions. (ii) We analyze aspects of quantum equilibration and thermalization in large quantum systems using the probe technique, mimicking elastic-dephasing effects and inelastic interactions on the dot. Here, unitary simulations based on the fermionic trace formula are accompanied by quantum Langevin equation calculations.


I. INTRODUCTION
There has been recently a great deal of interest in simulating the real-time dynamics of quantum systems, open or closed, prepared in a nonequilibrium state [1]. These investigations have been spurred by recent experimental breakthroughs in the ability to watch out-of-equilibrium dynamics, for example, in cold atomic gases [2], or in on-chip superconducting circuits [3]. This endeavor is fundamentally important for resolving basic issues in quantum dynamics, and in particular, for understanding equilibration and thermalization in quantum systems [1]. The nonequilibrium dynamics of the eminent Anderson model [4], composed of a single electronic level (quantum dot) coupled to two metals, has in particular been of great interest. This is because it is perhaps the simplest platform for probing both equilibrium and out-of-equilibrium physics in a many-body system. The model is integrable, even when electron-electron (e-e) repulsion effects are accounted for on the dot, and its integrability has been exploited for resolving its transport behavior [5].
A central strategy in most analytic and numerical tools devoted to the Anderson model, and impurity systems at large, is the separation of the total system into a subsystem (dot) and the environment (metals, referred to as reservoirs). The latter are typically assumed to be infinitedissipative and are maintained in one of the canonical ensembles of statistical mechanics. This assumption allows one to treat the effect of the reservoirs on the subsystem within a self-energy term. However, once the reservoirs are traced out, one cannot describe their explicit dynamics.
Among the numerical approaches developed along these lines we list the time-dependent numerical renormalization-group method [6], real-time diagrammatic Monte Carlo techniques [7] and path integral approaches [8,9]. These methods place focus on quantities such as the dot occupancy, transmission probability, conductance, current, noise, and correlations on the impurity. The dynamics of the total system, including the electron reservoirs, has not yet been explored since general tools for simulating the overall dynamics in a system-bath scenario are still missing.
The current work develops analytical and numerical treatments of global system evolution based on established impurity dynamics techniques. These tools allow investigation of the roles of e-e interactions and decoherence and dissipation effects on nonequilibrium reservoirs dynamics. We focus on the finite-size Anderson model composed of two metallic grains weakly coupled through a single electronic level. We refer to the metal grains, each composed of N ∼ 100 − 500 electronic states and n ∼ 50 − 200 electrons as "reservoirs" alluding to their high density of states (DOS).
This large DOS allows the reservoirs' effect on the dot (subsystem) to be absorbed into a positive real self-energy function lending to a quantum Langevin equation (QLE) description [10][11][12], as we explain below. A schematic representation is presented in Fig. 1. We are interested in following the real-time dynamics of both the dot and the reservoirs degrees of freedom. As an initial condition, we assume that each reservoir is prepared in a distinct Gibbs-like grand canonical state at a different chemical potential but at the same temperature.
In the absence of dephasing and inelastic effects, the dynamics of the total density matrix is followed by extending three approaches: (i) the quantum Langevin equation method [10][11][12], adopted here both in the noninteracting limit and in the mean-field (MF) regime, (ii) fermionic trace formula [13], used here for simulating the exact dynamics of the noninteracting model, and (iii) an influence functional path integral method [14,15], employed to treat interactions beyond the perturbative regime. In the latter half of the paper, these techniques are used to study reservoir population evolution both without and with Coulomb repulsion effects on the dot, and in the presence of emulated dephasing and inelastic scattering effects.
While Coulomb interactions are explicitly introduced here, the inclusion of dephasing and inelastic effects warrants further discussion. The origin of such processes are many-body interactions in the system, e.g., electron phonon coupling. Since an explicit and exact inclusion of these interactions is extremely challenging [16][17][18][19], phenomenological techniques have been developed in their stand, [20][21][22]. In the case of elastic decoherring processes the technique is referred to as a "dephasing probe". In the case of inelastic scattering processes, it is referred to as a "voltage probe". These probes are electron reservoirs, prepared such that, there is no either energy resolved or total net electron flow from the L-dot-R system towards these probes. For a scheme of this model, see Fig. 1(b). It should be noted that elastic-decoherring processes or inelastic effects are only emulated here by the probes. The overall dynamics can be still simulated using the unitary trace formula technique [23]. We also extend the QLE method to include a probe, and time evolve the system. Since our calculations provide the real-time dynamics of the full density matrix (DM), the process of equilibration and thermalization in a finite quantum system can now be studied [23][24][25][26]. Particularly, we find that when only decoherence effects are allowed, the system approaches a non-canonical equilibrium state. In contrast, when inelastic processes are included, the reservoirs relax towards a common Gibbs-like state. In Sec. IV B quantum equilibration and thermalization is investigated using the probe technique.
In this case, the total density matrix is resolved using the QLE method and the fermionic trace formula. The paper is summarized, along with an outlook, in Sec. V.

II. MODEL
The closed-system Anderson model consists two metal grains, ν = L, R, including (each) a collection of N ν dense electronic levels initially populated by noninteracting electrons up to the chemical potential µ ν , at temperature T = β −1 . The two baths couple only through their (weak) hybridization with a single level quantum dot. Work presented in this study concerns three variants of the model. The simplest version is the "noninteracting case", where electron-electron repulsion effects and any decoherring and relaxation mechanisms are excluded. The second case, the "interacting model", allows for e-e interactions on the dot only. The third model variant, the "probe model", phenomenologically contains elastic decoherring and inelastic scattering processes on the dot, using the probe technique and excluding e-e repulsion effects. This model is discussed in detail in Sec. IV B, where it is applied in the context of quantum equilibration.
A. The interacting model In the absence of decoherence and dissipation effects the interacting Hamiltonian takes the form where H L,R,W represents the Hamiltonian for the left reservoir, right reservoir and the dot, respectively. The term V ν denotes the coupling of the dot to the ν reservoir, Here, c k,σ (k = l, r, d) are fermionic operators of the left reservoir, l ∈ L, right reservoir, r ∈ R and the dot (d). The symbol σ stands for the spin state (↑ or ↓) and U accounts for the onsite repulsion energy. We assume that v l and v r are real numbers and that the Hamiltonians of the leads are diagonal in momentum basis and define the hybridization Γ ν (ǫ) = π k∈ν v 2 k δ(ǫ − ǫ k ), taken in practice to be energy independent. The Hamiltonian (2) disregards magnetic fields, yielding spin-degenerate energy levels, thus it is sufficient to consider observables for one spin species. We note that the noninteracting case arises simply from the suppression of U .
Our objective in this paper is to calculate the time evolution of the expectation values of all two-body operators in the system (k, j = l, r, d) written here in the Heisenberg representation with ρ(t 0 ) = ρ d ⊗ ρ L ⊗ ρ R representing the factorized time-zero density matrix of the system, and with the trace performed over all degrees of freedom.
We suppress the spin degree of freedom in the density matrix since its elements are identical for the two spin configurations. As an initial condition, we take the dot to be empty and the reservoirs' with the population f L,R (ǫ) = [e β(ǫ−µ L,R ) + 1] −1 . As a convention, we use the symmetric chemical potential bias µ L = −µ R > 0.
B. The probe model The Anderson probe model, a variant of the basic model, Eq. (2), can emulate memory loss and energy redistribution in a quantum system without explicitly introducing many-body interactions [20][21][22]. The probe technique has been of extensive use in mesoscopic physics, for describing the disappearance of quantum effects in transport [21], dissipation [20], and equilibration dynamics [27]. Recent advances include a full-counting statistics analysis of the probe model [28], and an extension of the probe technique to the AC regime [29]. We model here either a dephasing probe, allowing for quasi-elastic decoherence processes, or a voltage probe, where inelastic effects are further mimicked. In both cases we suppress electron-electron interaction effects in the system.
As we explain next, in our study the "probe" terminology refers to a setup slightly different from the conventional one. The standard construction refers to an open system scenario, where the probe practically performs which-path experiments through repetitive measurements of the system [30]. In contrast, in our picture the probe is a finite-closed quantum system, only initialized with a certain-special distribution. After its preparation, the probe, similarly to other parts of the system, is left undisturbed. Thus, we can use exact unitary approaches and simulate the dynamics of the total system. While this picture abuses to some extent the standard notion of a "probe", we maintain this terminology here since practically our implemented probe acts like a proper one, inducing phase loss or/and energy reorganization in the system.
We introduce a probe into the model by adding an additional Fermi-sea reservoir, denoted by the letter G, to the Hamiltonian (2), again discarding the spin degree of freedom, where We naturally define the hybridization Γ G (ǫ) = π g v 2 g δ(ǫ − ǫ g ), and take it as a constant. As always, our objective here is the resolution of all system expectation values of two-body operators (k, j = l, r, d), ρ k,j (t). As initial conditions we assume Eq. (4), where the G bath initial condition is set according to the particular probe condition, explained below.
Voltage probe. Inelastic scattering effects of electrons on the dot are effectively included by implementation of a voltage probe. The probe has a canonical distribution, f G (ǫ) = [e β(ǫ−µ G ) +1] −1 , and its chemical potential µ G is set such that the net charge current from the dot to the G unit vanishes for all times With the motivation to explore situations beyond the linear response regime [31], we retrieve µ G numerically, by employing the Newton-Raphson method [32], G is the initial guess, i ′ G denotes the first derivative with respect to µ G . In principle, one should adjust µ G throughout the simulation, to eliminate population leakage from the L-dot-R system into G. However, we have found in our simulations that the G bath has lawfully behaved as a probe once we determined µ G from the steady-state limit using the following analytic expression for the charge current with Γ = Γ L + Γ R + Γ G [33]. The lower and upper integration limits are determined by the band simulated. Substituting Eq. (9) in Eq. (7), a voltage probe condition is set by demanding f G (ǫ) to fulfill the relation Dephasing probe. Implementation of the dephasing probe, fabricating elastic decoherence, necessitates the stronger requirement i G (ǫ) = 0, i.e., the charge current at a given energy should vanish.
Using the steady-state behavior (9), we obtain a non-Fermi distribution We emphasize that µ G or f G have been determined here in the steady-state limit, assuming fixed chemical potentials for the L and R baths. Indeed, at short time, Γ L,R t 2, before a (quasi) steady-state sets in, we find that i G = 0. However, we have confirmed numerically that beyond this time throughout all our simulations |i G /i L,R | < 10 −4 , thus the G reservoir plays the role of a proper probe.
Three different approaches for the calculation of the full DM are described in Sections III A, III B and III C. Applications are included in Sec. IV.

A. Quantum Langevin Equation
The dynamics of the Anderson model in the absence of interactions (U = 0), with interactions at the mean-field level, or with a probe, is described here within a quantum Langevin equation framework [10]. The basis of our method has been used in the past to follow the dot evolution or the charge and energy currents in the system [11,12]. Here, we show results for the full DM. We begin our analysis with the trivial treatment of the impurity (dot) and review the steps involved. This review helps highlight underling approximations and establishes limits for the method's applicability.

Noninteracting case (U = 0)
In the Heisenberg representation the fermionic operators satisfy the following equations of mo- Formal integration of the reservoirs EOM yields, e.g. at the L end, We substitute Eq. (13), and the analogous expression for c r (t), into the dot EOM [Eq. (12)], and In this (exact) equation the second and third terms are interpreted as "noise" [10], The last two terms in Eq. (14) can be reduced, each, into decay terms, further inducing an energy shift of the dot energy, absorbed into the definition of ǫ d . This is justified by following two assumptions: (i) The hybridization Γ ν (ǫ) = π k∈ν v 2 k δ(ǫ − ǫ k ) may be taken as a positive real self-energy function [10], and (ii) the dot dynamics is slow relative to the reservoirs' evolution. We now explain the steps involved. First, we define a new operator for the dot, by absorbing its fast where ǫ kj ≡ ǫ k − ǫ j . We then change variables, x ≡ t − τ , and make the assumption that the dot evolution (now missing the fast phase oscillation) is slow with respect to other time scales in the with D L (ǫ) = l δ(ǫ − ǫ l ) as the density of states of the L metal, taken as flat here. We also take the interaction parameters v l to be independent of the l index. The imaginary term introduces an energy shift, which can be absorbed into the definition of ǫ d . It diminishes when the density of states does not depend on energy (the case used later), and when the bandwidth is large enough.
In our numerical calculations we have used a finite bandwidth with a cutoff D = ±1, introducing a small correction to ǫ d . We return to Eq. (14) and conclude that it obeys the quantum Langevin with Γ(ǫ) = Γ L (ǫ) + Γ R (ǫ). The dynamics of the dot occupation, n d (t) , can be reached by a formal integration of Eq. (17), to provide the standard expression [34] This derivation relied on the initial conditions (4). The summation runs over all the reservoirs degrees of freedom [35]. We now use Eq. (13) and its analogous expression for c r (t), together with Eq. (18), and derive analytical expressions for other quadratic expectation values, c † k (t)c j (t) , k, j = l, r, d. These results are valid as long as one can faithfully rely on Eq. (17). In what follows we take t 0 = 0 to simplify our notation. The reservoir-dot coherence can be obtained analytically, Here, B 1 includes contributions from the L side only, B 2 includes electron transmission pathways from the L side, through the dot, to the K = L, R grain, Using Eq. (20), we derive an expression for the charge current at the L contact, Here ℑ stands for the imaginary part. An analogous expression can be written for i R (t). Ref. [34] includes a Green's function based derivation for the time dependent current in the symmetric limit . This derivation results in a surplus nonphysical term at the initial time. We now turn our attention to the reservoirs' states population. Using Eq. (13), we find that it is given by three contributions, The two-times correlation functions can be obtained from Eqs. (13) and (18)   can be obtained with QLE when t > τ d . Thus, the QLE method can be used within the interval t < τ d only, to be consistent with its underlying assumptions. However, interesting nonequilibrium physics takes place within this window, thereby making this approach valuable considering that exact computational schemes are very expensive.

Mean-field theory
The QLE description of Sec. III A 1 can be generalized to accommodate electron-electron repulsion effects on the dot at the mean-field level. We refer to this extension as a MF QLE treatment, and note that it is not trivial: While a MF theory has been developed, suffering from some patholo- gies, for the study of dot occupation or charge current in the steady-state limit [37,38], here we present a MF scheme to describe the real-time dynamics of the full density matrix. By comparing MF results to exact numerical simulations, see Sec. IV, we conclude that a MF description can The effectiveness of the method also delicately depends on the dot level position, see for example Fig. 7.
The MF prescription, treating Coulombic repulsion, takes us back to Eq. (2). We now assume that the many-body interaction term can be factorized [37,38] This assumption reduces the Hamiltonian to an effectively noninteracting one, with a renormalized dot energyǫ The spin-index has been dropped here, as we choose not to study magnetic effects. The formalism could be feasibly generalized to include magnetic fields, resulting in ǫ d,↑ = ǫ d,↓ . In such situations the validity of MF equations is governed by another energy scale besides U/Γ and U/∆µ, namely The dot occupation is determined in a self-consistent manner at every instant by modifying Eq. (19) to contain the dot renormalized energy, The solution provides the renormalized dot energyǫ d (t), which is then used to replace ǫ d in Eqs. (20), (24), (A3) and (A4), to provide the full DM at the MF level.

Probe model
To implement elastic dephasing or inelastic effects with a probe, the set of equations (12) is augmented by an additional equation for c g . The EOM for c d must be modified to include its coupling to the G bath,ċ It can be easily shown that under the QLE basic assumptions, as discussed in Sec. III A 1, the dot still satisfies Eq. (17) with an additional noise term η G and with a re-defined total hybridization, Γ = Γ L +Γ R +Γ G . The noise η G obeys a relation analogous to Eq. (15), and Eq. (18) We can now recognize that, in the presence of the probe, the expressions for the DM elements (19), We describe here an exact brute force calculation that can provide numerically all the elements of the density matrix in the noninteracting case. We begin without the presence of a probe, opting to include its effects later. This unitary method complements the QLE description, whose validity is governed by τ d . Since the method is unitary, a recurrences behavior is expected to manifest at long enough time. The core of the method is the trace formula for fermions [13] Tr where m p is a single-particle operator corresponding to a quadratic operator M p = i,j (m p ) i,j c † i c j . c † i (c j ) are fermionic creation (annihilation) operators. Our objective is the dynamics of a quadratic operator A, either given by system or bath degrees of freedom, A ≡ c † j c k , j, k = l, r, d, We introduce the λ parameter, taken to vanish at the end of the calculation. The initial condition is taken to be factorized, ρ(t 0 ) = ρ d ⊗ ρ L ⊗ ρ R , ρ ν = e −β(Hν −µν Nν ) /Z ν , Z ν is the partition function, ρ d describes the dot initial density matrix. These density operators follow an exponential form, e M , with M a quadratic operator. The application of the trace formula leads to Here The fermionic trace formula can be trivially generalized to include a probe. We add the G bath into the expectation value expression, A ≡ c † j c k , where as before ρ ν = e −β(Hν −µν Nν ) /Z ν , Z ν is the partition function, ν = L, R, G. ρ d stands for the dot initial density matrix, and we trace over all DOF, the two reservoirs, the probe, and the dot.
It is achieved numerically-exactly using the fermionic trace formula (30).
Timescale. Previous studies for dense reservoirs have confirmed that INFPI can provide accurate results in both short time and in the quasi steady-state region [14,15]. However, the method is not restricted to such dense-reservoirs situations, and it can describe the dynamics of small metallic grains since it handles all states explicitly. It should be still noted that the basic working assumption behind INFPI is the existence of a finite bath-induced decorrelation time. If the metal grains are very small, including few discrete states, this memory time τ c does not exist or it becomes large, hindering convergence. Roughly, one could expect that a decorrelation time can be identified when a system-bath picture still holds, in the sense that a QLE description can be written i.e., Eq. (17) is valid. In such situations, INFPI simulations should converge and generally hold beyond τ d . In practice, since these calculations are intensive, we have computed dynamics within a relatively short interval, Γt < 5, where the QSS description is still valid.

IV. APPLICATIONS
We now turn our attention to applications of the preceding methods. We first study the effects of Coulombic interactions on the reservoirs' DOF evolution. We later investigate the equilibration process in the system, through the implementation of probes.

A. Anderson model with electron-electron interactions
In this section we study the evolution of the finite-size Anderson model with or without interactions, based on the three methods described earlier in Sec. III. As mentioned above, these techniques provide the dynamics of the total DM. While the fingerprints of many-body effects are disguised in the time evolution of conventional quantities, e.g., in the dot occupation and the charge current, they are well manifested in the reservoirs' population dynamics, allowing us to discern microscopic many-body scattering processes from single-particle events.
The population p(ǫ k ) = c † k c k of both reservoirs, in the noninteracting case, is displayed in Fig.  3 at different times. We note that results obtained using the QLE framework of Sec. III A perfectly agree with numerically-exact fermionic trace formula simulations. The three panels present results or below it (c), a dot-assisted tunneling feature develops, with population transfer taking place around available states that are the nearest in energy to ǫ d . The dynamics shown in Fig. 3 is reversible, with a characteristic time τ d ∼ 2π/∆E, ∆E = 2D/N is the mean spacing between energy levels and N = N L,R is the number of states in the L, R baths [36]. At this characteristic time the dot population begins to vary from its QSS value due to finite size effects. This behavior can be captured with trace formula simulations, but not within the QLE approach.
In Fig. 4 we display the dynamics under a mean-field QLE treatment with parameters corresponding to Fig. 3. While we are mindful of the technique's known pathologies [38], we stress that this calculation provides an intuitive understanding of the role of interactions: Within MF, the effect of finite U is to shift features in concert with the renormalized dot energy,ǫ d (t) = ǫ d +U n d (t) [Eq. (26)]. Interestingly, panel (c) demonstrates a change in transport mechanism, from a dotassisted tunneling at small U , to resonance transmission at large U , since the renormalized dot energy enters the bias window at a large enough interaction strength. Therefore, e-e interactions can enhance or suppress electronic transport, depending on the dot bare energy position. Overall, we conclude that MF simulations can reproduce the dynamics for this set of parameters, up to U/Γ 2. Qualitative features are correct through U/Γ ∼ 6.
Convergence of INFPI is verified with respect to the time step adopted, δt, and the memory time accounted for, τ c . Representative convergence curves for p(ǫ l = ǫ d ) are depicted in Fig.   6. While the time-step used does not affect our results, we note that the data is not yet fully converged with respect to τ c . This slow convergence could be attributed to the long decorrelation time experienced by an electron residing on any particular bath level, since its decorrelation process should take place by following a two-step procedure: the electron should first leave the particular bath state and populate the quantum dot. From the dot, it may subsequently transfer to any other bath state. One should also note that we display here a convergence curve for a particular level.
It is more accurate, but computationally demanding, to look at the overall evolution of p(ǫ) (for all ǫ) with τ c . This is because the resonant feature may change its magnitude with τ c , as we see here, as well as its absolute position. Fig. 6 only analyzes the effect of τ c on the peak magnitude.  In Fig. 7 we study the dynamics of the reservoirs' levels population when the renormalized dot energy sits above the bias window,ǫ d > µ L . The bare energy is taken at ǫ d = 0. 25  In this case the population shows a high-energy tail that is further enhanced with respect to MF data, representing significant deviations from a single particle description. However, as we did not manage to fully converge these results, they are not included here.
The breakdown of the single-particle picture is difficult to discern in cumulative quantities such as the charge current and the dot occupation, the latter is presented in Fig. 8, where we follow the time evolution of this quantity using four different techniques: MF QLE equations, first-order perturbation theory method [34], Monte-Carlo simulations [41][42][43][44] and INFPI [15]. The comparison shows a very good agreement between the latter two exact methods up to Γt ∼ 1 [15]. Before addressing the equilibration problem in detail we present in Fig. 9 a more standard quantity, the steady-state dot occupation. This will serve us for motivating the study of the total DM for resolving transport mechanisms. The dot occupation is displayed here as a function of Γ G , the dot-probe hybridization strength, and we show results using either a voltage probe or a dephasing probe, for different values of the dot energy position. We find that the dot occupation is insensitive to the probe condition, an observation that can be proved analytically by studying the long time behavior of Eq. (19) under either probes, Eqs. (10) or (11). It is interesting to note the crossover behavior of the dot occupation when its energy is placed above the bias window: When Γ G < Γ L,R occupation grows linearly with Γ G . However, it decays as Γ −1 G at large values, when effective dephasing and inelastic effects are strong. This behavior is similar to the thermally assisted tunneling behavior observed when using more detailed modeling [46]. It can similarly be shown that the steady-state current in the system is identical irrespective of the probe condition.
Note that we restrict ourselves to the quasi steady-state region since the G bath does not serve as a proper probe before quasi steady-state sets in.
The underlying transport mechanisms are therefore obscured in cumulative quantities (current, occupation) in this probe model, as well as in more explicit electron-phonon modeling [47]. Details about the involved mechanisms can be resolved by studying, e.g. the current noise, inelastic electron tunneling spectra [47], and the evolution of the reservoirs population [48], as we show below in the context of quantum equilibration and thermalization. The equilibration problem in quantum mechanics could be considered within different setups: a closed system [49][50][51], a system-bath scenario [52,53], or taking peer quantum systems [25]. Here, we consider the Anderson model with a probe, excluding e-e interactions, and simulate the following setup: At time t = t 0 we put into contact through a quantum dot two reservoirs each separately prepared in grand canonical states at chemical potentials µ ν and temperature β −1 . Electrons on the dot are susceptible to either elastic-decoherring processes or inelastic effects, mimicked by coupling to the relevant probe. For a schematic representation, see Fig. 1(b). Given this scenario, we investigate whether the two reservoirs can equilibrate or even thermalize in time, and furthermore, the nature of the equilibrium state. As we show below, when only elastic dephasing effects are mimicked, the system approaches a non-canonical equilibrium state. When inelastic processes are emulated, the two reservoirs relax towards a common canonical state. It should be noted that these results can be obtained for a finite and closed system, under a unitary evolution [23].
We identify thermal equilibration in our peer quantum system setup, by adjusting the conditions of Refs. [52,54], demanding that: (i) The system should equilibrate, i.e., evolve towards some particular state, and stay close to it for almost all time. Furthermore, the equilibrium state should be (ii) independent of the dot properties-energetics and initial state, (iii) insensitive to the precise initial state of each reservoir, (iv) close to diagonal in the energy basis of its eigen-Hamiltonian, and (v) a canonical state.
We use the trace formula, an exact unitary method, and follow the reservoirs' mutual equilibration process. We evolve the system using either a dephasing probe or a voltage probe, see Fig. 10, up to the time where recurrence features start to manifest, found here to scale as τ rec ∝ i=L,R,G N i .
As a reference, panel (a) displays results for the model without a probe, showing the development of a resonance feature around the dot energy position at ǫ d = 0. A clear evolution towards an equilibrium state is demonstrated when a probe is presented. With a dephasing probe, (b)-(c), the population of the two reservoirs relax to a two-step function with p(µ R < ǫ < µ L ) ∼ 0.5.
Because electrons from the L grain loose their phase memory on the dot, half populate the R side, on average, in the long time limit. This equilibrium state is sensitive to the precise details of the initial electron distribution, as energy redistribution is not allowed. We build a large G to delay recurrence behavior, but note that results at earlier times do not depend on the size of G, reinforcing the observation that G acts as an agent in driving the L-R mutual equilibration. When inelastic effects are mimicked with a voltage probe, and Γ G is large enough, panel (f), the system approaches a Gibbs-like thermal state-a step function at zero temperature. Results are shown up to the time τ rec at which recurrence features develop, which emerges here before full thermalization takes place. In order to achieve full thermalization one should further increase the size of the G bath, so as to delay recurrences. Alternatively, a dissipative mechanism could be introduced into G, e.g. by building a hierarchy of its interactions with the L-R system. Using a smaller value for Γ G , a non-canonical equilibrium distribution develops (d)-(e), reflecting the contribution of coherent and (effectively) incoherent electrons in the dynamics. It is also interesting to compare panels (c) and (e), displaying the equilibration progress for a dephasing probe and a voltage probe, respectively, while maintaining the value of Γ G . We find that the characteristic timescale to reach equilibrium is very similar in both cases. Thus, while the probe type dictates the structure of the equilibrium state, it does not affect the equilibration timescale. Fig. 11 shows that while under coherent evolution the resonance peak emerges around the energy ǫ d , in the presence of a voltage probe with (large enough) Γ G , the buildup of the equilibrium state systematically occurs around the equilibrium Fermi energy. This holds even when the dot is placed outside the bias window (not shown). Analogous trends take place when allowing for dephasing only.
A thermal equilibrium state should be diagonal in the energy eigenbasis of its Hamiltonian [52].
In Fig. 12 we display the density matrix ρ k,k ′ = c † k c k ′ , k, k ′ = l, r, excluding diagonal elements ρ k,k , without a voltage probe (a)-(b), and with a one (c)-(d), using the QLE technique. This quantity is expected to oscillate in the long time limit since the Hamiltonian is not diagonal in the are scattered, yet they appear more prominently around the equilibrium Fermi energy, ǫ = 0. These three features should become more pronounced at longer times, which can be simulated using the trace formula approach.

V. CONCLUSION
We have extended analytical and numerical methods, developed for simulating the dynamics of impurities, i.e., subsystems attached to large reservoirs, to reveal the dynamics of the total system.
As an example, we have focused on the Anderson model, a quantum dot coupled to two metal grains, and obtained the evolution of the total density matrix, focusing on the reservoirs' evolution from an initial nonequilibrium state. We have studied the noninteracting model, as well as a model with interactions and a probe model, emulating elastic dephasing and dissipation effects.
The three methods presented are the analytic quantum Langevin equation approach, a simulation based on a trace formula, and an exact numerical path integral scheme that can accommodate e-e repulsion effects. Notably, the extension of the QLE treatment to provide the total DM is of general importance as it can be used in multitude of other systems, as long as one can identify a "subsystem" within the total system. Making use of the methods developed, we have investigated the total system dynamics in the presence of distinct effects: (i) e-e interactions on the impurity, and (ii) dephasing and inelastic scattering effects. Addressing the prior, our calculations allow us to energy resolve the effect of e-e interactions on electron transfer in the Anderson dot model. In the resonant regime we found that the dynamics observed for noninteracting electrons is largely preserved up to U/Γ 2, and the main effect of interactions on the reservoirs' occupation is apparently a simple shift in the position of features affected by the renormalization of the dot energy. Away from resonance, in the tunneling domain, the presence of weak interactions already manifested itself in scattering electrons to high energy levels, an effect that is not captured within a mean-field treatment. In the case of the later effect, we found numerically that the presence of dephasing and inelastic effects on a weak link only can lead to global system equilibration and even thermalization. It is important to note that no restrictions were enforced on the metals' band structure and the dot energy. This is significant in light of many other studies in which equilibration requires the "nondegenerate energy gap" condition to be satisfied [25,49,52].
Future directions include the study of finite temperature and electron-electron interaction effects in the equilibration process [15], and the behavior given a quantum dot chain between the two metal grains. In a linear chain of impurities we expect that the coherent-diffusive crossover in the charge current behavior [55] would similarly manifest itself in the energy reorganization process of the reservoirs. The methods developed here could also be adopted for the study of bosonic systems, e.g., to describe the dynamics of bosonic degrees of freedom interacting with harmonic baths.