Iterative path integral summation for nonequilibrium quantum transport

# Iterative path integral summation for nonequilibrium quantum transport

S. Weiss\Ast1 1 Theoretische Physik, Universität Duisburg-Essen and CENIDE, Lotharstr.1, 47048 Duisburg, Germany
2 Institut für Theoretische Physik IV, Heinrich-Heine-Universität Düsseldorf, Universitätsstr.1, 40225 Düsseldorf, Germany
3 Departement Physik, Universität Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland
4 I. Institut für Theoretische Physik, Universität Hamburg, Jungiusstraße 9, 20355 Hamburg, Germany
R. Hützen2 1 Theoretische Physik, Universität Duisburg-Essen and CENIDE, Lotharstr.1, 47048 Duisburg, Germany
2 Institut für Theoretische Physik IV, Heinrich-Heine-Universität Düsseldorf, Universitätsstr.1, 40225 Düsseldorf, Germany
3 Departement Physik, Universität Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland
4 I. Institut für Theoretische Physik, Universität Hamburg, Jungiusstraße 9, 20355 Hamburg, Germany
D. Becker3 1 Theoretische Physik, Universität Duisburg-Essen and CENIDE, Lotharstr.1, 47048 Duisburg, Germany
2 Institut für Theoretische Physik IV, Heinrich-Heine-Universität Düsseldorf, Universitätsstr.1, 40225 Düsseldorf, Germany
3 Departement Physik, Universität Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland
4 I. Institut für Theoretische Physik, Universität Hamburg, Jungiusstraße 9, 20355 Hamburg, Germany
J. Eckel2 1 Theoretische Physik, Universität Duisburg-Essen and CENIDE, Lotharstr.1, 47048 Duisburg, Germany
2 Institut für Theoretische Physik IV, Heinrich-Heine-Universität Düsseldorf, Universitätsstr.1, 40225 Düsseldorf, Germany
3 Departement Physik, Universität Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland
4 I. Institut für Theoretische Physik, Universität Hamburg, Jungiusstraße 9, 20355 Hamburg, Germany
R. Egger2 1 Theoretische Physik, Universität Duisburg-Essen and CENIDE, Lotharstr.1, 47048 Duisburg, Germany
2 Institut für Theoretische Physik IV, Heinrich-Heine-Universität Düsseldorf, Universitätsstr.1, 40225 Düsseldorf, Germany
3 Departement Physik, Universität Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland
4 I. Institut für Theoretische Physik, Universität Hamburg, Jungiusstraße 9, 20355 Hamburg, Germany
M. Thorwart4 1 Theoretische Physik, Universität Duisburg-Essen and CENIDE, Lotharstr.1, 47048 Duisburg, Germany
2 Institut für Theoretische Physik IV, Heinrich-Heine-Universität Düsseldorf, Universitätsstr.1, 40225 Düsseldorf, Germany
3 Departement Physik, Universität Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland
4 I. Institut für Theoretische Physik, Universität Hamburg, Jungiusstraße 9, 20355 Hamburg, Germany
XXXX, revised XXXX, accepted XXXX
###### Abstract
\abstcol

We have developed a numerically exact approach to compute real-time path integral expressions for quantum transport problems out of equilibrium. The scheme is based on a deterministic iterative summation of the path integral (ISPI) for the generating function of nonequilibrium observables of interest, e.g., the charge current or dynamical quantities of the central part. Self-energies due to the leads, being nonlocal in time, are fully taken into account within a finite memory time, thereby including non-Markovian effects. Numerical results are extrapolated first to vanishing (Trotter) time discretization and, second, to infinite memory time. The method is applied to nonequilibrium transport through a single-impurity Anderson dot in the first place. We benchmark our results in various regimes of the rich parameter space. In the respective regime of validity, ISPI results are shown to match those of other state-of-the art methods. Especially, we have chosen the mixed valence regime of the Anderson model to compare ISPI to time dependent density matrix renormalization group (tDMRG) and functional RG calculations. Secondly, we determine the nonequilibrium current through a molecular junction in presence of a vibrational mode. We have found an exact mapping of the single impurity Anderson-Holstein model to an effective spin-1 problem. In analytically tractable regimes, as the adiabatic phonon or weak molecule-lead coupling regime, we reproduce known perturbative results. Studying the crossover regime between those limits shows that the Franck-Condon blockade persists in the quantum limit. At low temperature, the Franck-Condon steps in the characteristics are smeared due to nonequilibrium conditions. The third system under investigation here is the magnetic Anderson model which consists of a spinful single-orbital quantum dot with an incorporated quantum mechanical spin-1/2 magnetic impurity. Coulomb interaction together with the exchange coupling of the magnetic impurity with the electron spins strongly influence the dynamics. We investigate the nonequilibrium tunneling current through the system as a function of exchange and Coulomb interaction as well as the real-time impurity polarization. From the real-time evolution of physical observables, we are able to determine characteristics of the time-dependent nonequilibrium current and the relaxation dynamics of the impurity. These examples illustrate that the ISPI technique is particularly well suited for the deep quantum regime, when all time and energy scales are of the same order of magnitude.

molecular quantum transport, molecular junctions, quantum dots, nonequilibrium path integrals, exchange coupling, electron-phonon coupling, Keldysh fomalism.
\mail

e-mail weiss@thp.uni-due.de, Phone:+49-203-379-2969

\published

XXXX

## 1 Introduction

Small condensed matter systems that behave as transistors have attracted numerous research activities in recent years. When placed in a transport setup, current and noise measurements are used to determine the properties of the usually strongly correlated quantum mechanical system. In as well as out of thermal equilibrium, several quantum many-body properties of systems such as quantum dots or single molecule junctions are accessible in experiments [1, 2]. By down-scaling electronic transistors to the nanometer scale, quantum dots (QD’s) allow to control electronic and/or spin properties of single electrons and might be used to impress and read out information. Single molecules, inheriting also vibrational degrees of freedom, are ideal candidates for such devices since they can rather easily be functionalized such that current switches are implemented in a nanoscale environment [3, 4, 5, 6, 7, 8]. Furthermore, a controlled design of functional groups attached to molecules seems in reach. The field of molecular electronics is at the intersection of interdisciplinary research areas, combining chemical and physical properties of molecules. Under nonequilibrium conditions transport through molecules challenges experimentalists as well as theoreticians. Nonequilibrium in this review is referred to as the situation that the Fermi energies, provided by source and drain electrodes and located ’left’ and ’right’ of the structure, are not equal such that electrons are transferred through the system by tunneling. The voltage is assumed to drop at the nanostructure. In such finite voltage bias situations, many interesting physical effects arise due to the quantum nature of the electrons and their strong Coulomb interaction among each other and/or by their interaction with molecular vibrations or local magnetic moments. Prominent examples are resonant tunneling or the Kondo effect [9, 10, 11]. Different techniques are used to approach the strong coupling limit. For instance, the regime of low energy, low temperature/bias has been approached by Fermi liquid theory [12], via interpolative schemes [13], using integrability concepts [14], or by the perturbative renormalization group (RG) [15, 16]. A perturbative RG analysis has been performed by Schoeller and König [17]. The nonequilibrium generalization of Wilson’s numerical RG approach [18] and of the fRG method [19] have been discussed. Transport features have also been discussed by perturbation theory in the interaction strength [20]. On the other hand, perturbatively treating the tunneling matrix elements is a powerful way of describing the incoherent regime, see Ref. [9]. Density matrix renormalization group techniques have been extended to the nonequilibrium regime [21, 22]. In addition the flow equation method allows to study the Toulouse point of the Kondo model [23].

In case that a molecule is placed between electrodes, an appropriate theoretical description is needed to deal with its characteristic vibrational (phonon) degrees of freedom [2, 24]. The interplay between mechanical and electronic degrees of freedom is of interest in other areas of physics as well, e.g., for inelastic tunneling spectroscopy [25], nanoelectromechanical systems [26], break junctions [27], and suspended semiconductor or carbon-based nanostructures [28, 29, 30, 31]. Including vibrational degrees of freedom via a simple Anderson Holstein (AH) model, where a spinless electronic level is coupled to a single oscillator mode, shows several effects as the Franck-Condon blockade, negative differential conductance, or current induced heating or cooling [32, 33]; for a review, see Ref. [25]. As for the interacting Anderson model, analytical approaches typically address different corners of parameter space, a full theory that connects those corners seems not in reach at present.

A third topic adressed in this review is the investigation of a magnetic QD. Those setups have been studied experimentally in ensembles which are particularly suited for the investigation by laser and electromagnetic fields [34, 35, 36, 37, 38, 39, 40]. They are designed with standard lithographic methods and are technologically well established. Moreover, embedding individual magnetic Mn ions into quantum dots and studying the electrical properties is possible [41, 42, 43]. Small quantum dots with few charge carriers and a single magnetic impurity may become important candidates for efficient high density spintronic devices.

There is a considerable need for numerical methods which describe small quantum systems out of equilibrium accurately and, ideally, treat simple model Hamiltonians numerically exactly. Numerical renormalization group [44] or quantum Monte Carlo (QMC) calculations [45, 46, 47, 48, 49, 50] provide a possible line of attack to those problems. Due to the dynamical sign problem, these calculations become increasingly difficult at low temperatures, but in several parameter regions, the stationary steady-state regime seems accessible. Based on non-standard ensembles, the steady state is described in a recent work by Han [51], using an imaginary-time QMC approach followed by a double analytical continuation scheme. This last step is numerically the most difficult part [52].

We here review the development of a novel numerical scheme denoted as iterative summation of real-time path integrals (ISPI), in order to address quantum transport problems out of equilibrium[53]. Many-body systems driven out of equilibrium are known [15, 54, 55] to acquire a steady state that may be quite different in character from their ground-state properties. Details of the steady state may depend on the nature of the correlations, as well as on the way in which the system is driven out of equilibrium. Our ISPI approach, described in detail below, provides an alternative and numerically exact method to tackle out-of-equilibrium transport in correlated quantum dots. Based on the evaluation of the full nonequilibrium Keldysh generating function, along with the inclusion of suitable source terms, observables of interest are computed. It builds on the fact that nonlocal in time correlations, induced by the fermionic leads, decay exponentially in the long-time limit at any finite temperature. Within a characteristic time , all correlations are taken into account, while for larger times, the correlations are dropped due to their exponentially small contributions. This allows us to construct an iterative scheme to evaluate the generating function. An appropriate extrapolation procedure allows to eliminate the Trotter time discretization error (the Hubbard-Stratonovich (HS) transformation below requires to discretize time) as well as the finite memory-time error. This yields the desired numerically exact value for the observables of interest. Note that the need for a finite memory time makes our approach difficult to apply at very low energies (). Fortunately, other methods are available in this regime. At finite or , the requirement of not too long memory times is exploited and the spin path summation remains tractable. Recently, also Segal et al., see Ref. [56] have provided an alternative formulation of the ISPI approach in terms of Feynman Vernon like influence functionals.

The ISPI scheme is implemented here for three characteristic impurity models. First, the single-level Anderson impurity model [57, 85, 59, 60, 61], second the spinless Anderson Holstein model [25], which mimics the behavior of a molecular quantum dot. Finally, when magnetic molecules are investigated, additional couplings to localized spin impurities are important. Results for the relaxation dynamics of such a system are presented in this article as well.

The present paper is organized as follows. In Sec. 2, we introduce the model for the quantum dot coupled to normal leads. We present the computation of the generating function for the nonequilibrium Anderson model there as well. The presence of an external source term allows to calculate the current as a functional derivative. In Sec. 3 we introduce the numerical iterative path integral summation method, from which we obtain observables of interest. We give a detailed discussion of the convergence properties of our method and describe the extrapolation scheme. For several sets of parameters, we will present results and benchmark checks in Sec. 4. In addition, results for the mixed valence regime of the Anderson model are presented. A good agreement with tDMRG and fRG is found. The Anderson Holstein model is studied by means of the ISPI method in Sec. 5, the main finding of a sustained Franck-Condon blockade at low temperatures is discussed upon validating our spin-1 mapping of the AH Hamiltonian. By far the most complex model of this work, a magnetic interacting QD is studied in Sec.6. We briefly review the necessary changes in the appearing Keldysh generating function and discuss our results on the relaxation time and impurity polarization as well as the tunneling current. A summary is given in Sec. 7.

## 2 Keldysh generating function for non-interacting impurity models

We consider the generic Hamiltonian for a quantum dot that is coupled to metallic leads to the left and right side ()

 H = Hdot+Hleads+HT (1) = ∑σE0σnσ+∑kpσ(ϵkp−μp)c†kpσckpσ −∑kpσ[tpc†kpσdσ+h.c.].

Here, with is the energy of a single electron with spin on the isolated dot. Tuning a back gate voltage or a Zeeman magnetic field term changes the value of . The latter is assumed not to affect the electron dispersion in the leads. The corresponding dot electron annihilation/creation operator is , with the density operator . Possible eigenvalues of are , corresponding to the empty or occupied electronic state with spin . Interactions on the quantum dot are treated in Sec. 2.1 and not considered for the moment. In Eq. (1), denotes the energies of the noninteracting electrons (operators ) in lead , with chemical potential . Quantum dot and leads are connected by the tunnel couplings . The observable of interest is the (symmetrized) tunneling current ,

 I(t)=−ie2∑kpσ[ptp⟨d†σckpσ⟩t−pt∗p⟨c†kpσdσ⟩t], (2)

where with . The stationary steady-state dc current follows as the asymptotic long-time limit, . We have explicitly confirmed that current conservation, , is numerically fulfilled for the ISPI scheme.

In the presence of a finite bias voltage, , the Keldysh technique [62, 63, 81] provides a way to study nonequilibrium transport. In this formalism, the time axis is extended to a contour with branches, see Ref. [81], along with an effective doubling of fields. The Keldysh Green function (GF) exhibits a matrix structure where denotes the contour ordering of times along the Keldysh contour, and correspond to fields representing lead or dot fermions, respectively. We omit the spin indices here, remembering that each entry still is a diagonal matrix in spin space. The Keldysh partition function contains all relevant information about the physics of the system. In order to obtain it, we first integrate over the noninteracting lead fermion fields. Subsequently, we integrate over the dot fields as well. In a fermion coherent state basis, the generating function is

 Z[η]=∫D[∏σ¯dσ,dσ,¯ckpσ,ckpσ]eiS[¯dσ,dσ,¯ckpσ,ckpσ], (3)

with Grassmann fields . The external source term, which allows to compute the current at measurement time , is chosen such that

 I(tm)=−i∂∂ηlnZ[η]∣∣∣η=0. (4)

Correspondingly, it is also possible to evaluate other observables, e.g., the zero-frequency shot noise, by introducing appropriate source terms and performing the corresponding derivatives. The action is , see Ref. [53] for the explicit expressions. After integrating over the leads’ degrees of freedom, the effective action for the dot becomes nonlocal in time. The generating function for the noninteracting system reads

 Zni[η]=∫D[∏σ¯dσdσ]ei(Sdot,0+Senv) (5)

with

 Senv = ∫Cdt∫Cdt′∑σ¯dσ(t){γL(t,t′)+γR(t,t′) (6) + ieη2[γL(t,t′)−γR(t,t′)] × [δ(t−tm)+δ(t′−tm)]}dσ(t′).

For the source term, the physical measurement time is fixed on the upper branch, hence the Keldysh element of the source term self-energy vanishes. The matrices in Eq. (6) represent the leads, their Fourier transforms in frequency space are explicitly given as Keldysh matrices

 γp(ω)=iΓp(2f(ω−μp)−1−2f(ω−μp)2−2f(ω−μp)2f(ω−μp)−1). (7)

With the usual assumption that the leads are in thermal equilibrium, . Taking the wide-band limit with a constant density of states per spin channel around the Fermi energy, the hybridization of the dot level with lead enters. We focus on symmetric contacts and on symmetrically applied bias voltages as well. The generalization to asymmetric contacts is straightforward.

In the next step (still for vanishing on-dot interactions), we integrate over the dot degrees of freedom. This yields the noninteracting generating function

 Zni[η]=∏σdet[iG−10σ(t,t′)+ηΣJ(t,t′)]. (8)

The function follows from

 G0σ(ω) = [(ω−ϵ0σ)τz−γL(ω)−γR(ω)]−1 = 11+[(ω−ϵ0σ)/Γ]2 × (ω−ϵ0σ+iΓ(1−F)iΓFiΓ(F−2)−ω+ϵ0σ+iΓ(1−F)),

where is the standard Pauli matrix in Keldysh space, and . Moreover, the self-energy for the source term is obtained as

 ΣJ(t,t′) = e2[γL(t,t′)−γR(t,t′)] (9) × [δ(t−tm)+δ(t′−tm)].

Up to this point, we have discussed the noninteracting case. The next section describes how to include the electron-electron interaction.

### 2.1 Hubbard-Stratonovich transformation

In the presence of on-dot Coulomb interactions, we add the Coulomb term

 Hint=Un↑n↓ (10)

to the Hamiltonian in Eq. (1). For our purpose, it is convenient to use the operator identity , which results in a shift of the single-particle energies , and we may rewrite . Now, the action in Eq. (3) in real time contains quartic terms of dot Grassmann fields and a Gaussian integration is not possible. We will use a time-discrete path integral for the following discussion [65]. In order to decouple the quartic term, we discretize the full time interval, , with the time increment . On each time slice, we perform a Trotter breakup of the dot propagator according to , where . According to Refs. [66, 67, 68] the emerging Trotter error can be systematically eliminated from the results [69, 70], see below. On a single Trotter slice, a discrete Hubbard-Stratonovich transformation [67, 71, 72, 73] allows to decouple the interaction.This locally introduces Ising-like discrete spin fields on the branches of the Keldysh contour with on the -th Trotter slice. For a given Trotter slice, we define

 e±iδtU(n↑−n↓)2/2=12∑s±=±e−δtλ±s±(n↑−n↓). (11)

The HS parameter is obtained from the equation,

 cosh(δtλ±)=cos(δtU/2)±isin(δtU/2),

under the condition that , see Ref. [53]. Note that the arbitrarily chosen overall sign of does not influence the physical result. Uniqueness of this HS transformation requires . To ensure sufficiently small time discretizations, we meet the condition in all calculations in general.

After the HS transformation, the remaining fermionic Grassmann variables appear quadratically and are integrated out at the cost of the full path summation over the discrete HS Ising spins according to

 Z[η]=∑{s}∏σdetG−1σ[{s},η]. (12)

The full Keldysh GF written in time-discretized () form is

 (G−1σ)αβkl[{s},η]=(G−10σ)αβkl+iηΣJ,αβkl−iδtδklλαsαkδαβ, (13)

where labels the Keldysh branches, and the noninteracting GF is

 G0σ,kl=∫∞−∞dω2πeiδt(k−l)ωG0σ(ω). (14)

Note that depends only on time differences due to time-translational invariance of the noninteracting part which holds at thermal equilibrium of the leads. The respective time discrete version of the self-energy kernel for the source term is, cf. Eqs. (6) and (9),

 ΣJ,αβkl=e2[γαβL,kl−γαβR,kl]δmkδα,++δmlδβ,+δt, (15)

where and the measurement time is .

## 3 Formulation of the iterative scheme

In this section we review the construction for the iterative expressions for the Keldysh partition function which form the ISPI scheme. For this, we exploit the property (see, e.g., Ref. [53, 73, 74, 75]) that each Keldysh component of in Eq. (14) decays exponentially at long time differences () for any finite temperature, see Eq. (14). The related time scale is denoted as correlation or memory time . In Fig. 1, we show typical examples of for different bias voltages . The exponential decrease with time is presented in the inset of Fig. 1(a), where the absolute value is plotted in log-linear representation for the same parameters. For large bias voltages and low enough temperatures, e.g., at and , the decay is superposed by an oscillatory behavior. Since the lead-induced correlation function decays as , the respective correlations decay on a time scale given by . Thus, the exponential decay suggests to neglect lead-induced correlations beyond . This motivates an iterative scheme which exactly accounts for correlations within an interval , while neglecting them outside. Notice that the exponential decay is only present for finite and/or , whereas at correlations die out only algebraically, and our approach is not applicable.

Let us then face the remaining path sum in Eq. (12). In the discrete time representation, we denote by the initial and final times and . The discretized GF and self-energy kernels for given spin are then represented as matrices of dimension . For explicit calculations, we arrange the matrix elements related to Keldysh space (characterized by the Pauli matrices ) and to physical times as . In particular, the ordering of the matrix elements from left to right (and from top to bottom) follows increasing times. The lead-induced correlations thus decrease exponentially with growing distance from the diagonal of the matrix.

For numerical convenience, we evaluate the generating function in the equivalent form

 Z[η]=N∑{s}∏σdetDσ[{s},η], (16)

with , see Ref. [53] for details, and thus get

 Dαβσ,kl[{s},η]=δαβδkl+iδtλαGαβ0σ,klsαk−iη∑j,α′Gαα′0σ,kjΣJ,α′βjl. (17)

By construction, we have to sum over auxiliary Ising spins, which appear line-wise. The total number of possible spin configurations is .

Next, we exploit the truncation of the GF by setting for , where is the correlation time, with the respective number of Trotter time slices. All GF matrices have a band structure whose band width is given by . Equivalently, we can use for the spin-spin correlation the definition

 sαi⋅sβj={sαi⋅sβj,if|i−j|≤K,0,else. (18)

We note that in the continuum limit, i.e., and , the approach is formally exact.

To proceed, we exploit that the determinant of a quadratic block matrix is given by . {widetext} In time space, we obtain the -Keldysh GF band matrix

 D≡D(1,NK)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝D11D1200…0D21D22D230…⋮0D32D33D34…⋮00D43D44…0⋮⋮⋮⋮⋱DNK−1NK0……0DNKNK−1DNKNK⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

where the single blocks are block matrices defined as ()

 Dll′=⎛⎜ ⎜ ⎜⎝D(l−1)K+1,(l′−1)K+1…D(l−1)K+1,l′K⋮⋱⋮DlK,(l′−1)K+1…DlK,l′K⎞⎟ ⎟ ⎟⎠.

Without loss of generality, the number of Trotter slices is chosen as integer such that and the matrix elements follow from Eq. (17). We keep their dependence on the Ising spins implicit. Each still has a Keldysh structure and a spin structure. If we apply the formula for the determinant from above, the generating function (16) is represented as

 Z[η] = ∑s±1,…,s±Ndet{D11[s±1,…,s±K]} (19) × det{D(2,NK)[s±K+1,…,s±N]−D21[s±K+1,…,s±2K][D11[s±1,…,s±K]]−1D12[s±K+1,…,s±2K]},

where the matrix is obtained from by removing the first line and the first column.

In order to set up an iterative scheme, we use the following observation: to be consistent with the truncation of the correlations after a memory time , we have to neglect terms that directly couple Ising spins at time differences larger than , see also Eq. (18), consequently matrix products of the form

 Dl+2,l+1[Dl+1,l+1]−1Dl+1,l[Dl,l]−1Dl,l+1[Dl+1,l+1]−1Dl+1,l+2=0 (20)

within the Schur complement in each further iteration step. We do not neglect the full Schur complement but only those parts which are generated in the second-next iteration step. With this, we rewrite the generating function as

 Z[η] = ∑s±1,…,s±Ndet{D11[s±1,…,s±K]}NK−1∏l=1det{Dl+1,l+1[s±lK+1,…,s±(l+1)K] (21) −Dl+1,l[s±lK+1,…,s±(l+1)K][Dl,l[s±(l−1)K+1,…,s±lK]]−1Dl,l+1[s±lK+1,…,s±(l+1)K]}.

Exchanging the sum and the product, and reordering the sum over all Ising spins, we obtain

 Z[η]=∑s±N−K+1,…,s±NZNK[s±N−K+1,…,s±N], (22)

where is the last element obtained from the iterative procedure defined by ()

 Zl+1[s±lK+1,…,s±(l+1)K]=∑s±(l−1)K+1,…,s±lKΛl[s±(l−1)K+1,…,s±lK,s±lK+1,…,s±(l+1)K]Zl[s±(l−1)K+1,…,s±lK]. (23)

The propagating tensor is read off from Eq. (21) as

 Λl = det{Dl+1,l+1[s±lK+1,…,s±(l+1)K] (24) −Dl+1,l[s±lK+1,…,s±(l+1)K][Dl,l[s±(l−1)K+1,…,s±lK]]−1Dl,l+1[s±lK+1,…,s±(l+1)K]}.

The iteration starts with .

The current is numerically obtained by evaluating Eq. (4) for a small but fixed value of ; we have taken for all results shown. The current as a function of time shows a transient oscillatory or relaxation behavior at short times. We extract the stationary value when the current has saturated to a plateau.

### 3.1 Convergence and extrapolation procedure

By construction there are two systematic errors in the scheme: (i) the Trotter error due to finite time discretization , and (ii) the memory error due to a finite memory time . The scheme becomes exact in the limit and . We can eliminate both errors from the numerical data in the following systematic way: Step 1: We choose a fixed time discretization and a memory time . A reasonable estimate for is the minimum of and (see above). With that, we calculate the current , and, if desired, the differential conductance (the derivative is performed numerically for a small ). The calculation is then repeated for different choices of and . Step 2: Next, the Trotter error can be eliminated by exploiting the fact that it vanishes quadratically for [66, 67, 68]. For a fixed memory time , we can thus extrapolate and obtain , which still depends on the finite memory time . The quadratic dependence on is illustrated in Fig. 2 (a) for different values of . Note that each line corresponds to the same fixed memory time . Step 3: In a last step, we eliminate the memory error by extrapolating to , and obtain the final numerically exact value . For the dependence on , we empirically find a regular and systematic behavior as shown in Fig. 2(b). The value is approached with corrections of the order of , see Fig. 2(b).

We have implemented the iterative scheme together with the convergence procedure on standard Xeon 2GHz machines. Computations are then only possible for due to the limited memory resources available. Typical running times for the shown simulation data are approximately hours for .

## 4 Benchmarking the approach: comparison with exact and perturbative results

In this section, we discuss the results obtained for the Anderson model. We measure energies in units of . Unless noted otherwise, all error bars for the shown data points, which are due to the Trotter and memory extrapolation scheme, are of the order of the symbol sizes in the figures.

We stop short on reporting that we have recovered the exact current for the solvable case, as a function of the bias and the gate voltage. By construction, the ISPI method includes all tunneling processes exactly to arbitrary orders.

In order to benchmark our code for finite , we compare the numerical results to a perturbative calculation at the charge degeneracy point , where the interaction self-energy can be computed up to second order in [20, 76]. For a detailed comparison, we plot mostly the interaction corrections, , with being the current , the linear conductance , or the nonlinear conductance , respectively. Figure 3 shows the results for as a function of the bias voltage for , and . For , we perfectly recover the perturbative results, which confirms the reliability of our code even in the regime of nonlinear transport. Clearly, the corrections are small and negative, which can be rationalized in terms of indications of Coulomb blockade physics, as transport is suppressed by a finite on-dot interaction. For , the current decreases even further, and the deviations between the ISPI and perturbative results increase. The relative deviation for is already , illustrating that perturbation theory is already of limited accuracy in this regime. Although it well reproduces the overall tendency, there is a significant quantitative difference. It is even more pronounced for , as shown in the figure. Here, second-order perturbation theory does not even reproduce qualitative features.

### 4.1 Comparison with a master equation approach

Next, we compare our numerically exact results with the outcome of a standard classical rate equation calculation [9, 85]. The rate equation is expected to yield reliable results in the incoherent (sequential) tunneling regime, when . Then, a description in terms of occupation probabilities for the isolated many-body dot states is appropriate. Results for are shown in Fig. 4, both for the interaction corrections to the nonlinear and the linear conductance. When the temperature is lowered, , quantum coherent interaction effects become more important, as seen from the exact ISPI results. They are clearly not captured by the master equation in the sequential tunneling approximation. However, for , interaction corrections are washed out, and the master equation becomes accurate, cf. Fig. 4. Similarly, from our ISPI results, we have found (data not shown) that interaction corrections are suppressed by an increasing bias voltage as well.

### 4.2 ISPI vs. tDMRG and fRG

The mixed valence regime is characterized by the fact that the QD’s occupation fluctuates between the different charge states, namely empty, singly and doubly occupied, see Ref. [77] for details. Furthermore the coupling to the leads is strong, i.e. the Kondo temperature , see below in Sec. 4.3. When controlling the gate voltage by tuning , it is possible to tune the nanostructure into this limit. From the theoretical point of view another energy scale enters, again, not as a small parameter and simulation methods are rare in this regime. For the case of intermediate Coulomb repulsion we have studied [75] the mixed valence regime and compare our results to other state-of the art methods, functional renormalization group (fRG) [19] and time dependent density matrix renormalization group (tDMRG)[21, 22]. The results are shown in Fig. 5 for the steady state current for , i.e., away from the charge degeneracy point. The results from fRG and ISPI match perfectly from small bias voltages up to the strong non-equilibrium regime. We note that by construction, it becomes increasingly cumbersome (and finally impossible) to obtain converged ISPI results in the limit of vanishing bias voltages and low temperatures [53], since then the correlations do not decay sufficiently to be truncated. In the present setup, tDMRG has a tendency to overestimate the currents in the mixed valence regime, see Ref. [78] for details, hence we find slight deviations between tDMRG and the two other methods (see , lower left panel). Overall agreement away from the symmetric point between the three methods is very good. Furthermore, curves show agreement for a wide range of magnetic fields, applied to the QD, see Ref. [75] for details.

### 4.3 Small bias regime eV≪Γ:

For sufficiently small bias voltage, the current is linear in , and we can focus on the linear conductance . Figure 6 shows for different magnetic fields , taking and (for ). For , two spin-degenerate transport channels contribute, and a single resonant-tunneling peak at results. For , the spin-dependent channels are split by , resulting in a double-peak structure. Furthermore, since also lifts the degeneracy, the spin-resolved levels are now located at due to the Zeeman splitting. We find an interaction-induced broadening, cf. Fig. 6, of the resonant-tunneling peak as compared to the noninteracting case. The width of the Lorentzian peak profile for is determined by at sufficiently low , and broadens as increases. Here, the double-peak structure, with two clearly separated peaks for finite , is not yet fully developed. The two peaks largely overlap, and the distance of the peaks is below the expected , since tunneling significantly broadens the dot levels.

Next, we address the temperature dependence of the linear conductance (numerically evaluated for ). In Fig. 7, we show for different values of (up to ) at . For , the deviations from the -result is small. For larger , deviations become more pronounced at low temperatures where the on-dot interaction becomes increasingly relevant. Up to present, we have obtained converged results in the regime of small bias voltages for interaction strengths for temperatures above or close to the Kondo temperature, . The corresponding Kondo temperatures are for and for . In the regime , we compare our results to the result of Hamann [77, 79],

 G(T)=e2h(1−ln(T/TKH)[ln2(T/TKH)+3π2/4]1/2), (25)

for the linear conductance, where , see Fig. 7 (solid lines). In Ref. [77], it has been shown that the results of the numerical RG coincide with those of Eq. (25) in this regime. Fig. 7 illustrates that the agreement between the two approaches is satisfactory and shows that the ISPI provides reliable results in the linear regime above or close to the Kondo temperature. Due to the construction of the approach, the situation is more favorable for large bias voltages, where short to intermediate memory times allow us to obtain convergent results.

### 4.4 Large bias regime eV≥Γ:

Let us then turn to nonequilibrium transport for voltages . Here, the transport window is given by , and a double-peak structure for emerges even for , see Fig. 8 with distance between the peaks. We show results for but otherwise the same parameters as in Fig. 6. For an additional finite magnetic field, each peak of the double-peak structure itself experiences an additional Zeeman splitting, resulting in an overall four-peak structure. For and the depicted values of , the two innermost peaks (closest to ) overlap so strongly that they effectively form a single peak at again. The two outermost peaks are due to the combination of the finite magnetic field and the bias voltage.

Increasing the on-dot interaction downsizes the differential conductance peaks as compared to the noninteracting case, i.e., the interaction corrections are again largest when the level energy matches the chemical potential in the leads. Note that the four-peak structure is already present in the noninteracting case (with ) and, hence, is not modified qualitatively by a finite .

Finally, we address the temperature dependence of the differential conductance . ISPI results for are shown in Fig. 9. Again, as in the linear regime, the conductance increases with lower temperatures, and finally saturates, e.g., at for and . Clearly the conductance decreases when the bias voltage is raised. Increasing renders this suppression yet more pronounced, see also inset of Fig. 9 for the corresponding corrections. At high temperatures, thermal fluctuations wash out the interaction effects, and the interaction corrections die out.

## 5 Sustained Franck-Condon blockade in molecular quantum dots

As a second example, we apply the ISPI scheme to investigate vibrational effects in the nonequilibrium tunneling current through a molecular quantum dot, see Ref. [73]. We extend the single-impurity Anderson model by a linear phonon of frequency (annihilation operator ) which couples to a single spinless electronic level with energy (operators ). Hence, we have the molecular Hamiltonian

 Hm=Ω b†b+[E0+λ(b+b†)]nd (26)

with electron-phonon coupling strength and .

The short-time propagator on the forward/backward branch of the Keldysh contour, , then allows for a Trotter breakup, , with , where the auxiliary relation

 e∓iδtH1=1−nd+nde−λ2δ2t/2e∓iδtE0e∓iδtλb†e∓iδtλb (27)

holds. This effectively decouples the electron-phonon interaction in terms of a three-state variable defined at each (discretized) time step along the forward/backward () part of the Keldysh contour, where . Below, we also use the notation with periodic boundary conditions on the Keldysh contour. The “Ising spin” variable picks up the three terms in Eq. (27) and acts like a Hubbard-Stratonovich auxiliary field, similar to the Ising field employed in the Hirsch-Fye formulation of the Anderson model [53, 66, 67, 68]. The bosonic (phonon) scalar field and the fermionic (dot and lead electrons) Grassmann fields appearing in the Keldysh path integral are noninteracting but couple to the time-dependent auxiliary spin variable. Hence, those fields can be integrated out analytically and the time-dependent current follows from a path summation as above in Sec. 2. The resulting matrix (in time and Keldysh space) depends on the complete spin path . Specifically, we obtain , where has spin-dependent matrix elements . We find only when , where it coincides with the usual (wide-band limit) expression [80]. Finally, the diagonal matrix (quoted here for ) with

 Bηη=Asηe−λ2δ2t∑η′αα′[iGph]η,η′+1|sηsη′| (28)

encapsulates all phonon effects, where is the discretized phonon Green’s function, see Ref. [81], and we have used the notation and . We comment shortly on the peculiar convergence properties of the present model. Convergence of the extrapolation requires intermediate or -values, for otherwise the necessary memory times become exceedingly long. For the results below, we have used and . The shown current follows by averaging over the -window, with error bars indicating the mean variance. Additional ISPI runs for and were consistent with these results, and we conclude that small error bars indicate that convergence has been reached. When dealing with summands in the evaluation of the generating function in the AH case as compared to summands for the Anderson model, our ISPI code to calculate runs for  CPU hours on a 2.93 GHz Xeon processor.

Here, we have once more convinced ourselves that the numerical ISPI results for the curves are consistent with known analytical theory in the respective parameter limits. We employ the following perturbative methods: (i) For , perturbation theory in the electron-phonon coupling applies and yields a closed expression for arbitrary values of all other parameters [82]. We note that the solution of the AH model with a very broad dot level [83, 84] corresponds to this small- regime. (ii) For high temperatures, , a description in terms of a rate equation is possible [80]. We here use the sequential tunneling approximation with golden rule rates [85]. For small , the corresponding results match those of perturbation theory, while in the opposite strong-coupling limit, the Franck-Condon blockade occurs and implies a drastic current suppression at low bias voltage [30, 86]. (iii) For small oscillator frequency, , the nonequilibrium Born-Oppenheimer (NEBO) approximation is appropriate and allows us to obtain from a Langevin equation for the oscillator [87, 88]. For small , this approach is also consistent with perturbative theory, while for high , NEBO and rate equation results are found to agree. For clarity, we focus on a resonant level with here. The case of weak electron-phonon coupling, is shown Fig. 10. We compare our ISPI data for to the results of perturbation theory in and of the rate equation. Perturbation theory essentially reproduces the ISPI data. The rate equation is quite accurate for high temperatures, but quantitative agreement with ISPI was obtained only for . We note that the ISPI error bars increase when lowering due to the growing memory time () demands.

Next, Fig. 11 shows ISPI results for a slow phonon mode, , with larger electron-phonon coupling . In that case, perturbation theory in is not reliable and likewise, the rate equation is only accurate at the highest temperature () studied, cf. the upper left inset of Fig. 11. However, we observe from Fig. 11 that for such a slow phonon mode, NEBO provides a good approximation for all temperatures and/or voltages of interest. We conclude that the ISPI technique is capable of accurately describing three different analytically tractable parameter regimes.

In the limit of strong electron-phonon coupling , the classical rate equation predicts a Franck-Condon blockade of the current for low bias and [86]. Sufficiently large can be realized experimentally, and the Franck-Condon blockade has indeed been observed in suspended carbon nanotube quantum dots [30]. For a nonequilibrated phonon with intermediate-to-large , understanding the Franck-Condon blockade in the quantum coherent regime of low temperature, , is an open theoretical problem. Here, multiple phonon excitation and deexcitation effects generate a complicated (unknown) nonequilibrium phonon distribution function, and the one-step tunneling interpretation in terms of Franck-Condon factors between shifted oscillator parabolas [86] is no longer applicable. We here study this question using ISPI simulations, which automatically take into account quantum coherence effects.

In Fig. 12, the crossover from weak to strong electron-phonon coupling is considered. The inset shows curves for , where we observe a current blockade for low voltages once . The blockade becomes more pronounced for increasing and is lifted for voltages above the polaron energy [86]. Remarkably, the Franck-Condon blockade persists and becomes even sharper as one enters the quantum-coherent regime (here, ), despite of the breakdown of the sequential tunneling picture. We also observe a nonequilibrium smearing of phonon step-like features in the curves in Fig. 12, cf. also Refs. [30, 86].

## 6 Nonequilibrium quantum dynamics in the magnetic Anderson model

The third example to which we have applied the ISPI scheme [74] is the magnetic Anderson model. We extend the Hamiltonian Eq. (1) in the presence of Coulomb interactions, see Eq. (10), by a magnetic impurity localized on the quantum dot which interacts via an exchange interaction with the spins of the confined electrons on the QD. The magnetic part of the Hamiltonian reads

 Himp+HJint= Δimp2τz+Jτz(d†↑d↑−d†↓d↓)H∥int+J2(τ+d†↓d↑+τ−d†↑d↓)H⊥int.

The generating function is obtained by integrating again over the corresponding Grassmann fields for dot and lead operators as well as the discrete paths and for the real spin variables and the HS Ising fields, respectively, and

 Z[η]=∑{τ,ζ}∫D[¯ckpσ¯dσckpσdσ](−1)ℓ(−iJδt2)mP[{τ}]eiS. (30)

The path sums over impurity and HS spin-fields are performed over the -tuples and with . Within an impurity path , flip-flop transitions occur on the Keldysh contour, where of them lie on the lower branch. The action includes tunneling and lead effects as in Eq. (6). Correspondingly, the magnetic part of the action is

 Simp=−Δimpδt2N∑k=2(τk−τ2N−k+1)=−Δimp2∫Kdtτ(t). (31)

The polynomial in Eq. (30) depends on the impurity path for the forward (backward) branch of the contour. Then, we collect all indices of the flips into the tuple (sorted in ascending order) along the forward path with for all . Accordingly, is the tuple of ascending flip indices along the backward path with for all . 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 electronic Grassmann fields as

 P[{τ}]:=∏j∈T−flip¯dj+1τjdj−τj∏k∈T+flip¯dk−τkd