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

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

Manas Kulkarni Chemical Physics Theory Group, Department of Chemistry, University of Toronto, 80 Saint George St. Toronto, Ontario, Canada M5S 3H6 Department of Physics, University of Toronto, 60 Saint George St. Toronto, Ontario, Canada M5S 1A7    Kunal L. Tiwari Chemical Physics Theory Group, Department of Chemistry, University of Toronto, 80 Saint George St. Toronto, Ontario, Canada M5S 3H6 Department of Physics, University of Toronto, 60 Saint George St. Toronto, Ontario, Canada M5S 1A7    Dvira Segal Chemical Physics Theory Group, Department of Chemistry, University of Toronto, 80 Saint George St. Toronto, Ontario, Canada M5S 3H6
July 5, 2019
###### Abstract

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.

###### pacs:
05.30.-d, 03.65.Aa, 03.65.Yz, 72.10.-d

## 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 Rev (). These investigations have been spurred by recent experimental breakthroughs in the ability to watch out-of-equilibrium dynamics, for example, in cold atomic gases bloch (), or in on-chip superconducting circuits houck (). This endeavor is fundamentally important for resolving basic issues in quantum dynamics, and in particular, for understanding equilibration and thermalization in quantum systems Rev (). The nonequilibrium dynamics of the eminent Anderson model anderson (), 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 rmk ().

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 infinite-dissipative 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 anders (), real-time diagrammatic Monte Carlo techniques marco () and path integral approaches egger1 (); egger2 (). 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 electronic states and 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 ford (); sen (); DharDM (), 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 ford (); sen (); DharDM (), adopted here both in the noninteracting limit and in the mean-field (MF) regime, (ii) fermionic trace formula israel (), used here for simulating the exact dynamics of the noninteracting model, and (iii) an influence functional path integral method dvira_pccp (); dvira_prb (), 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 Eran (); Thoss (); Thorwart (); Andrei (), phenomenological techniques have been developed in their stand, Buttiker1 (); Buttiker2 (); Beenakker (). 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 -dot- 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 qle (). 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 h1 (); h2 (); barra (); qle (). 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.

The paper is organized as follows. In Sec. II we present the finite-size closed Anderson model and outline the implementation of the probes. In Section III we present our developed numerical and analytic treatments of the density matrix dynamics: First, we extend the standard quantum Langevin equation approach to include reservoirs dynamics. The method can treat both the noninteracting model (III.1.1), the case with interactions, only at the level of mean-field (Hartree) theory (III.1.2), and the probe model (III.1.3). Second, we present the fermionic trace formula, useful for studying the Anderson model without interactions and with implemented dephasing and inelastic effects in Sec. III.2. The third method, presented in Sec. III.3, is an influence-functional path integral approach dvira_prb (); dvira_pccp (). This non-perturbative tool can treat the model with interactions in a numerically exact manner. Applications are included in Sec. IV. The effects of Coulomb repulsion effects on the dot are studied using mean-field QLE and the path integral technique in Sec. IV.1. In Sec. IV.2 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, , including (each) a collection of dense electronic levels initially populated by noninteracting electrons up to the chemical potential , at temperature . 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.2, where it is applied in the context of quantum equilibration.

### ii.1 The interacting model

In the absence of decoherence and dissipation effects the interacting Hamiltonian takes the form

 H=HL+HR+HW+VL+VR, (1)

where represents the Hamiltonian for the left reservoir, right reservoir and the dot, respectively. The term denotes the coupling of the dot to the reservoir,

 HL = ∑l,σϵlc†l,σcl,σ,HR=∑r,σϵrc†r,σcr,σ VL = ∑lσvlc†d,σcl,σ+h.c.VR=∑rσvrc†d,σcr,σ+h.c. HW = ∑σϵdc†d,σcd,σ+Und,↑nd,↓. (2)

Here, () are fermionic operators of the left reservoir, , right reservoir, and the dot (). The symbol stands for the spin state ( or ) and accounts for the onsite repulsion energy. We assume that and are real numbers and that the Hamiltonians of the leads are diagonal in momentum basis and define the hybridization , 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 .

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(t)≡⟨c†k(t)cj(t)⟩≡Tr[ρ(t0)c†k(t)cj(t)], (3)

written here in the Heisenberg representation with 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’ DM to be diagonal,

 ⟨c†d(t0)cd(t0)⟩=0,⟨c†l(t0)cl(t0)⟩=fL(ϵl)≡fl,⟨c†r(t0)cr(t0)⟩ = fR(ϵr)≡fr, (4)

with the population . As a convention, we use the symmetric chemical potential bias .

### ii.2 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 Buttiker1 (); Buttiker2 (); Beenakker (). The probe technique has been of extensive use in mesoscopic physics, for describing the disappearance of quantum effects in transport Buttiker2 (), dissipation Buttiker1 (), and equilibration dynamics Buttiker2b (). Recent advances include a full-counting statistics analysis of the probe model Buttiker3 (), and an extension of the probe technique to the AC regime Buttiker4 (). 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 ExpVolt (). 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 , to the Hamiltonian (2), again discarding the spin degree of freedom,

 HP=H+HG+VG, (5)

where

 HG=∑gϵgc†gcg,VG=∑gvgc†gcd+h.c. (6)

We naturally define the hybridization , and take it as a constant. As always, our objective here is the resolution of all system expectation values of two-body operators (), . As initial conditions we assume Eq. (4), where the 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, , and its chemical potential is set such that the net charge current from the dot to the unit vanishes for all times

 iG≡ddt∑g⟨c†gcg⟩=0. (7)

With the motivation to explore situations beyond the linear response regime DharSC (), we retrieve numerically, by employing the Newton-Raphson method NR (),

 μ(m+1)G=μ(m)G−iG(μ(m)G)/i′G(μ(m)G). (8)

is the initial guess, denotes the first derivative with respect to . In principle, one should adjust throughout the simulation, to eliminate population leakage from the -dot- system into . However, we have found in our simulations that the bath has lawfully behaved as a probe once we determined from the steady-state limit using the following analytic expression for the charge current

 iG(ϵ) = 2ΓGπΓR[fG(ϵ)−fR(ϵ)]+ΓL[fG(ϵ)−fL(ϵ)](ϵ−ϵd)2+Γ2, iG = ∫iG(ϵ)dϵ, (9)

with comm (). 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 to fulfill the relation

 ∫dϵfG(ϵ)(ϵ−ϵd)2+Γ2=1ΓL+ΓR∫dϵfL(ϵ)ΓL+fR(ϵ)ΓR(ϵ−ϵd)2+Γ2. (10)

Dephasing probe. Implementation of the dephasing probe, fabricating elastic decoherence, necessitates the stronger requirement , i.e., the charge current at a given energy should vanish. Using the steady-state behavior (9), we obtain a non-Fermi distribution

 fG(ϵ)=ΓRfR(ϵ)+ΓLfL(ϵ)ΓR+ΓL. (11)

We emphasize that or have been determined here in the steady-state limit, assuming fixed chemical potentials for the and baths. Indeed, at short time, , before a (quasi) steady-state sets in, we find that . However, we have confirmed numerically that beyond this time throughout all our simulations , thus the reservoir plays the role of a proper probe.

Three different approaches for the calculation of the full DM are described in Sections III.1, III.2 and III.3. Applications are included in Sec. IV.

## Iii Methods

### iii.1 Quantum Langevin Equation

The dynamics of the Anderson model in the absence of interactions (), with interactions at the mean-field level, or with a probe, is described here within a quantum Langevin equation framework ford (). 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 sen (); DharDM (). 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.

#### iii.1.1 Noninteracting case (U=0)

In the Heisenberg representation the fermionic operators satisfy the following equations of motion (EOM),

 ˙cd = −iϵdcd−i∑lvlcl−i∑rvrcr ˙cl = −iϵlcl−ivlcd ˙cr = −iϵrcr−ivrcd (12)

Formal integration of the reservoirs EOM yields, e.g. at the end,

 cl(t) = e−iϵl(t−t0)cl(t0)−ivl∫tt0dτe−iϵl(t−τ)cd(τ)dτ. (13)

We substitute Eq. (13), and the analogous expression for , into the dot EOM [Eq. (12)], and retrieve

 ˙cd = −iϵdcd−i∑lvle−iϵl(t−t0)cl(t0)−i∑rvre−iϵr(t−t0)cr(t0) (14) − ∫tt0dτ∑lv2le−iϵl(t−τ)cd(τ)−∫tt0dτ∑rv2re−iϵr(t−τ)cd(τ).

In this (exact) equation the second and third terms are interpreted as “noise” ford (),

 ηL(t) ≡ ∑lvle−iϵl(t−t0)cl(t0) ηR(t) ≡ ∑rvre−iϵr(t−t0)cr(t0). (15)

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 . This is justified by following two assumptions: (i) The hybridization may be taken as a positive real self-energy function ford (), 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 oscillatory behavior, . Its EOM is

 ˙~cd = −i[ηL(t)+ηR(t)]eiϵdt−∫tt0dτ∑lv2le−iϵld(t−τ)~cd(τ)−∫tt0dτ∑rv2re−iϵrd(t−τ)~cd(τ) (16)

where . We then change variables, , and make the assumption that the dot evolution (now missing the fast phase oscillation) is slow with respect to other time scales in the system, . This results in

 ∫tt0dτ∑lv2le−iϵld(t−τ)~cd(τ)≈~cd(t)∑lv2l∫t−t00e−iϵldxdx

If the time is long, , integration gives

 ∑lv2l∫t→∞0e−iϵldxdx=ΓL(ϵd)−2iv2llimt→∞[∫+∞−∞DL(ϵ)sin2(ϵ−ϵd)tϵ−ϵddϵ]

with as the density of states of the metal, taken as flat here. We also take the interaction parameters to be independent of the index. The imaginary term introduces an energy shift, which can be absorbed into the definition of . 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 , introducing a small correction to . We return to Eq. (14) and conclude that it obeys the quantum Langevin equation

 ˙cd=−iϵdcd−iηL(t)−iηR(t)−Γ(ϵd)cd(t), (17)

with . The dynamics of the dot occupation, , can be reached by a formal integration of Eq. (17),

 cd(t)=cd(t0)e(−iϵd−Γ)(t−t0)−i∫tt0e(−iϵd−Γ)(t−τ)[ηL(τ)+ηR(τ)]dτ, (18)

to provide the standard expression Komnik ()

 ⟨nd(t)⟩ ≡ ⟨c†d(t)cd(t)⟩=∑k=l,r|vk|2fkϵ2dk+Γ2[1+e−2Γ(t−t0)−2e−Γ(t−t0)cos[ϵdk(t−t0)]]. (19)

This derivation relied on the initial conditions (4). The summation runs over all the reservoirs degrees of freedom foot_dense (). We now use Eq. (13) and its analogous expression for , together with Eq. (18), and derive analytical expressions for other quadratic expectation values, , . These results are valid as long as one can faithfully rely on Eq. (17). In what follows we take to simplify our notation. The reservoir-dot coherence can be obtained analytically,

 (20)

Here, includes contributions from the side only,

 B1=ivlflΓ−iϵdl[1−e−t(Γ−iϵdl)]. (21)

includes electron transmission pathways from the side, through the dot, to the grain,

 B2=−ivl∑k∈L,Rv2kfkΓ2+ϵ2dk ×{1−eiϵkltiϵlk+e−2Γt−e−t(Γ−iϵdl)iϵld−Γ−e−t(Γ−iϵkd)−e−iϵlktiϵld−Γ−e−t(Γ−iϵdk)−e−t(Γ+iϵld)iϵlk}. (22)

Using Eq. (20), we derive an expression for the charge current at the contact,

 iL(t) ≡ ddt∑l⟨c†l(t)cl(t)⟩=−2I∑lvl⟨c†d(t)cl(t)⟩ (23) = 2ΓLΓRπ∑lfL−fRϵld2+Γ2−2ΓLπe−Γt∑l1ϵld2+Γ2 × {e−Γt(flΓL+frΓR)−(2flΓL+2frΓR−Γfl)cos(ϵldt)−flϵldsin(ϵldt)}.

Here stands for the imaginary part. An analogous expression can be written for . Ref. Komnik () 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,

 p(ϵl) ≡ ⟨c†l(t)cl(t)⟩=⟨c†l(t0)cl(t0)⟩ (24) + ivle−iϵl(t−t0)∫tt0eiϵl(t−τ)⟨c†d(τ)cl(t0)⟩dτ+c.c. + v2l∫tt0∫tt0dτ1dτ2⟨c†d(τ1)cd(τ2)⟩eiϵl(t−τ1)e−iϵl(t−τ2).

The two-times correlation functions can be obtained from Eqs. (13) and (18) without additional approximations, and the explicit expressions are given in Appendix A. Similarly, closed analytic expressions can be written for inter and intra-reservoir coherences, e.g. , , see Appendix A. Eqs. (19), (20), (24), (A3) and (A4), and the analogous -bath expressions, form the time-dependent full density matrix of the system.

Timescales. We now comment on the applicability of the QLE approach. Given infinite reservoirs, a current carrying steady-state behavior develops, and the dot occupation, as well as the charge current, reach a fixed value after a short time, (see for example Fig. 8). However, since the reservoirs are finite in the present treatment, recurrence effects should eventually manifest in our system. These effects cannot be handled by the QLE technique since an irreversible behavior has been assumed for the dot, as Eq. (17) breaks unitary evolution. The technique can still excellently reproduce the exact dynamics in the so called “quasi steady-state” (QSS) region, up to yuka (). Here is the number of electronic states in each bath and is the band cutoff. Within this time, the dot occupation and the charge current are constant, similar to a real steady-state situation. Around the time the dot occupation should begin to vary, showing (partial) recurrence behavior, and the QSS limit breaks down.

Fig. 2 clarifies this timescale issue. The left panel displays the dot occupation as a function of time using either the QLE approach (full line) or an exact method (dashed line), described in Sec. III.2. Results agree up to , and deviations before this time are due to the finite band used in QLE, while neglecting the energy correction in Eq. (17). Around the time exact simulations show a partial reversal of the dot occupation, while QLE still produces the QSS value. At a later time, , the QLE data diverges. Panel (b) presents the bath occupation for two selected energies, and we find that nonphysical values, such as a population exceeding unity, can be obtained with QLE when . Thus, the QLE method can be used within the interval 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.

#### iii.1.2 Mean-field theory

The QLE description of Sec. III.1.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 pathologies, for the study of dot occupation or charge current in the steady-state limit KomnikMF (); BerMF (), 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 produce physical results up to . 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 KomnikMF (); BerMF ()

 Und,↑nd,↓→U[⟨nd,↑⟩nd,↓+nd,↑⟨nd,↓⟩]. (25)

This assumption reduces the Hamiltonian to an effectively noninteracting one, with a renormalized dot energy

 ~ϵd(t)=ϵd+U⟨nd(t)⟩. (26)

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 . In such situations the validity of MF equations is governed by another energy scale besides and , namely . The dot occupation is determined in a self-consistent manner at every instant by modifying Eq. (19) to contain the dot renormalized energy,

 (27)

The solution provides the renormalized dot energy , which is then used to replace in Eqs. (20), (24), (A3) and (A4), to provide the full DM at the MF level.

#### iii.1.3 Probe model

To implement elastic dephasing or inelastic effects with a probe, the set of equations (12) is augmented by an additional equation for . The EOM for must be modified to include its coupling to the bath,

 ˙cg = −iϵgcg−ivgcd ˙cd = −iϵdcd−i∑lvlcl−i∑rvrcr−i∑gvgcg. (28)

It can be easily shown that under the QLE basic assumptions, as discussed in Sec. III.1.1, the dot still satisfies Eq. (17) with an additional noise term and with a re-defined total hybridization, . The noise obeys a relation analogous to Eq. (15), and Eq. (18) is generalized to

 cd(t)=cd(t0)e(−iϵd−Γ)(t−t0)−i∫tt0e(−iϵd−Γ)(t−τ)[ηL(τ)+ηR(τ)+ηG(τ)]dτ. (29)

We can now recognize that, in the presence of the probe, the expressions for the DM elements (19), (20), (24), (A3) and (A4) stay formally intact. The technical adjustments are as follows: (i) We re-define the total hybridization, . (ii) We augment summations that run over both and baths by terms. For example, the summation in [Eq. (A2)] should include such terms. (iii) We set the bath distribution to satisfy the probe conditions, explained in Sec. II.

### iii.2 Fermionic trace formula (U=0)

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 . 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 israel ()

 Tr[eM1eM2...eMp]=det[1+em1em2...emp], (30)

where is a single-particle operator corresponding to a quadratic operator . () are fermionic creation (annihilation) operators. Our objective is the dynamics of a quadratic operator , either given by system or bath degrees of freedom, , ,

 ⟨A(t)⟩=Tr[ρ(t0)eiHtAe−iHt]=limλ→0∂∂λTr[ρLρRρdeiHteλAe−iHt]. (31)

We introduce the parameter, taken to vanish at the end of the calculation. The initial condition is taken to be factorized, , , is the partition function, describes the dot initial density matrix. These density operators follow an exponential form, , with a quadratic operator. The application of the trace formula leads to

 ⟨eλA(t)⟩=det{[IL−fL]⊗[IR−fR]⊗[Id−fd]+eihteλae−ihtfL⊗fR⊗fd}. (32)

Here, and are single-body matrices of the and operators, respectively. The matrices and are the identity matrices for the space and for the dot. The functions and are the band electrons occupancy . Here they are written in matrix form and in the energy representation. represents the dot initial occupation, again written in a matrix form. Since we are working with finite-size reservoirs, Eq. (32) can be readily simulated numerically-exactly.

The fermionic trace formula can be trivially generalized to include a probe. We add the bath into the expectation value expression, ,

 ⟨A(t)⟩=Tr[ρ(t0)eiHPtAe−iHPt]=limλ→0∂∂λTr[ρLρRρGρdeiHPteλAe−iHPt], (33)

where as before , is the partition function, . stands for the dot initial density matrix, and we trace over all DOF, the two reservoirs, the probe, and the dot.

Timescales. Simulations with the trace formula are not restricted to a certain time scale. The method is unitary, providing (physical) recurrence behavior due to finite size effects. Since the time evolution scheme is not iterative, the accuracy of results does not deteriorate in time.

### iii.3 Numerically exact path integral simulations, U≠0

The time evolution of the closed and interacting Anderson model can be simulated by employing a numerically-exact iterative influence-functional path integral (INFPI) approach dvira_prb (); dvira_pccp (). This method relies on the fact that in out-of-equilibrium (and nonzero temperature) cases bath correlations have a finite range, allowing for their truncation beyond a memory time dictated by the voltage-bias and the temperature. Based on this finite-memory assumption, an iterative-deterministic time-evolution scheme has been developed, where convergence with respect to the memory length can, in principle, be reached. The principles of the INFPI approach have been detailed in Refs. dvira_prb (); dvira_pccp (), where it has been developed to investigate dissipation effects in the nonequilibrium spin-fermion model, and the population and current dynamics in correlated quantum dots. Recently, it has been used to examine the effects of a magnetic flux on the intrinsic coherence dynamics in an Aharonov-Bohm quantum dot interferometer AB (). The INFPI method relies on the existence of a finite decorrelation time, thus it is suited for simulating the dynamics of an impurity coupled to a bath. Here we show that it can be used to retrieve the total DM in a system-bath setup. While in principle the method could encompass both interactions and probe, we focus exclusively on the first element.

The method is based on the fermionic trace formula (30), incorporating many-body effects within a path integral expression. Our work starts with the time evolution expression (31) under the Hamiltonian (2). We factorize the time evolution operator, , , and adopt the Trotter decomposition , where with

 H0 = ∑ν=L,R(Hν+Vν)+∑σ(ϵd+U2)c†d,σcd,σ H1 = U[nd,↑nd,↓−12(nd,↑+nd,↓)]. (34)

extracts many-body interactions on the dot, and it is eliminated by introducing auxiliary Ising variables via the Hubbard-Stratonovich (HS) transformation Hubb-Strat (),

 e±iH1δt=12∑seH±(s), eH±(s)≡e−sκ±(nd,↑−nd,↓). (35)

Here, , , . The uniqueness of this transformation requires that . Incorporating the Trotter decomposition and the HS transformation into Eq. (31), we find that the time evolution of is dictated by

 ⟨A(t)⟩=limλ→0∂∂λ{∫ds±1ds±2...ds±NtI(s±1,s±2,...,s±Nt)}. (36)

The integrand, referred to as as the “Influence Functional” (IF), is given by (, )

 I(s±q,...,s±q+p)=122(p+1)Tr[ρ(t0)G+(s+q+p)...G+(s+q)eiH0(q−1)δteλAe−iH0(q−1)δtG−(s−q)...G−(s−q+p)],

where and . Eq. (36) is exact in the limit. Practically, it is evaluated by truncating the IF beyond a memory time , corresponding to the time beyond which bath correlations may be ignored dvira_prb (), is an integer. The following (non-unique) breakup has been suggested by dvira_prb (),

 I(s±1,s±2,...s±Nt)≃I(s±1,s±2,...,s±Ns)Is(s±2,s±3,...,s±Ns+1)...Is(s±Nt−Ns+1,s±Nt−Ns+2,...,s±Nt), (38)

where each element in the product, besides the first one, is given by a ratio between truncated IFs,

 Is(sq,sq+1,...,sq+Ns−1)=I(s±q,s±q+1,...,s±q+Ns−1)I(s±q,s±q+1,...,s±q+Ns−2). (39)

We now define a multi-time object,

 R(s±q+1,s±q+2,...,s±q+Ns−1) ≡∑s±1,s±2,...,s±qI(s±1,s±2,...,s±Ns)Is(s±2,s±3,...,s±Ns+1)...×Is(s±q,s±q+1,...,s±q+Ns−1), (40)

and evolve it iteratively by multiplication with the subsequent truncated IF, followed by summation over the time variables at the head,

 R(s±q+2,s±q+3,...,s±q+Ns)=∑s±q+1R(s±q+1,s±q+2,...,s±q+Ns−1)Is(s±q+1,s±q+2,...,s±q+Ns). (41)

The behavior at a particular time is reached by summation over the internal variables,

 ⟨eλA(tq)⟩=∑s±q+2−Ns,...,s±qR(s±q+2−Ns,s±q+3−Ns,...,s±q). (42)

This procedure is repeated for several (small) values of , and the expectation value is retrieved by numerical differentiation in . The truncated IF, Eq. (LABEL:eq:IF), is the core of this calculation. 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 dvira_prb (); dvira_pccp (). 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 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 . In practice, since these calculations are intensive, we have computed dynamics within a relatively short interval, , 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.

### iv.1 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 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.1 perfectly agree with numerically-exact fermionic trace formula simulations. The three panels present results using different values for the dot energy. (a) When the dot energy is placed within the bias-window () a resonance feature develops around the position of the dot level, with a dip (peak) showing in the () bath. In contrast, if the dot energy is positioned either above the bias window (b) 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 . The dynamics shown in Fig. 3 is reversible, with a characteristic time , is the mean spacing between energy levels and is the number of states in the baths yuka (). 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 BerMF (), we stress that this calculation provides an intuitive understanding of the role of interactions: Within MF, the effect of finite is to shift features in concert with the renormalized dot energy, [Eq. (26)]. Interestingly, panel (c) demonstrates a change in transport mechanism, from a dot-assisted tunneling at small , to resonance transmission at large , 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.

Mean-field results are compared to numerically-exact INFPI simulations in Fig. 5 for and , with the bare dot energy centered within the bias window. Data was produced by time evolution of all , , expectation values up to . In agreement with MF QLE results, the basic effect of e-e interactions observed here is a shift in the resonance position. Overall, we conclude that MF simulations can reproduce the dynamics for this set of parameters, up to . Qualitative features are correct through .

Convergence of INFPI is verified with respect to the time step adopted, , and the memory time accounted for, . Representative convergence curves for 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 . 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 (for all ) with . This is because the resonant feature may change its magnitude with , as we see here, as well as its absolute position. Fig. 6 only analyzes the effect of on the peak magnitude.