Transient dynamics in interacting nanojunctions within self-consistent perturbation theory

We present an analysis of the transient electronic and transport properties of a nanojunction in the presence of electron-electron and electron-phonon interactions. We introduce a novel numerical approach which allows for an efficient evaluation of the non-equilibrium Green functions in the time domain. Within this approach we implement different self-consistent diagrammatic approximations in order to analyze the system evolution after a sudden connection to the leads and its convergence to the steady state. These approximations are tested by comparison with available numerically exact results, showing good agreement even for the case of large interaction strength. In addition to its methodological advantages, this approach allows us to study several issues of broad current interest like the build up in time of Kondo correlations and the presence or absence of bistability associated with electron-phonon interactions. We find that, in general, correlation effects tend to remove the possible appearance of charge bistability.


Introduction
For decades, studies of quantum transport in nanoscale devices have mainly focused on steady state properties [1]. While the potential interest of transient dynamics was pointed out long ago [2,3] such studies have recently received an increasing attention in connection with advances in experimental techniques for time-resolved measurements [4][5][6][7][8][9][10][11][12][13]. These studies are also motivated by the important technological goal of increasing the operation speed of devices while reducing their energy consumption. Moreover, studies of the transient dynamics after a quench of a given parameter are currently undertaken in many fields of physics ranging from cold atoms [14,15], correlated materials [16], dynamical phase transitions [17] and, more generally, in connection to the question on the existence of a well defined stationary state for any given model of interacting particles [18].
In this work we present an efficient algorithm for the integration of the time-dependent Dyson equation for the non-equilibrium Green functions applied to different models of correlated nanoscale systems, including electron-electron and electron-phonon interactions. To deal with these correlations we use a diagrammatic expansion of the system self-energies at different levels of approximation including self-consistency effects. In the case of electron-phonon interactions we introduce novel theoretical tools for solving the Dyson equations associated with the phonon propagator in order to account properly for the build up of a non-equilibrium phonon population. As a check of these approximations we study the convergence of the system properties like mean charge, current and spectral density to their stationary values and also compare them to available numerically exact results. When not available we have implemented our own NRG calculations. We show how this time-dependent approach is quite convenient for including self-consistency in a straightforward way. We exemplify the use of this methodology to investigate the issue of bistability for the molecular junction, demonstrating how the inclusion of correlation effects beyond the mean-field approximation tends to eliminate the bistable behavior of charge and current for certain parameter regimes.
The paper is organized as follows: in section 2 we introduce the formalism and the numerical techniques used for computing the transient electronic and transport properties; in section 3 we analyze the dynamics of a system with strong electron-electron interactions taking the non-equilibrium Anderson model as a paradigmatic example. Section 4 is devoted to the study of the transient properties in the presence of electronphonon interactions by means of the spinless Anderson-Holstein model. In section 5 we consider a situation where both electron-electron and electron-phonon interactions are present using the spin-degenerate Anderson-Holstein model. Finally we present the conclusions and provide a brief overlook of our main results in section 6.  where c kσ,ν , with ν=L(R) labeling the left (right) electrode, and c 0σ are annihilation operators for electrons in the leads and in the central region respectively and v kσ,ν (t) is the tunneling amplitude which will depend on time.

Keldysh formalism for the transient regime
The two electrodes can be kept at different chemical potentials via a constant bias voltage eV=μ L −μ R . For simplicity the central region will consist of a single quantum level denoted by ò 0 . The last term, H int , in Ĥ describes the many body interactions in the central region, which we shall specify later. Hereafter we assume In what follows we will consider the wide-band approximation for the electrodes. Within this approximation the tunneling rates can be taken as constants where ρ F is the density of states at the Fermi edge, the resonant level width being Γ=Γ L +Γ R . Our aim is to analyze the transient dynamics of such a correlated system after a sudden quench of the coupling to the electrodes at an initial time that we take at t=0. Thus, v kσ,ν (t)=θ(t)v kσ,ν , which allows us to define a time-dependent tunneling rate Γ(t)=θ(t)Γ. Although this work is focused on this sudden connection case, more general time-dependent Hamiltonians could be considered within the formalism presented below.
The dynamical electronic and transport properties can be obtained from the central level Green functions in where T ĉ is the chronological time-ordering operator along the Keldysh contour [76] (see figure 1(a)). In the absence of interactions the problem is exactly solvable even in the presence of an arbitrary time dependent potential [2,3]. However, in the presence of interactions the problem of obtaining the dynamical behavior of the system usually becomes extraordinarily demanding. On the one side, there is the problem of finding an appropriate treatment of correlation effects by means of an adequate selfenergy. This is not always a simple task in the dynamical problem. On the other hand, even if an appropriate selfenergy is found, the numerical solution of the Dyson equation for the Keldysh propagators (which in the time domain becomes an integral equation) is a formidable numerical problem.
In this section we present an efficient numerical procedure for the calculation of the Keldysh propagators in the transient regime. It allows us to obtain accurate results for the electronic and transport properties such as the central region charge and current. The power of the method is additionally checked by analyzing the convergence of these quantities (together with the central region spectral density) to their expected stationary values.
We start from the Dyson equation for the central level Green function in Keldysh space, which can be formally inverted is the inverse free electron propagator of the uncoupled central level, T , S ŝ the tunneling self-energy and ,int S ŝ the interaction self-energy. Interactions mixing the spin degree of freedom could be also included in the equation as discussed in [77,78]. Equation (2) can be numerically solved by discretizing time in the Keldysh contour (see figure 1(a)). From now on the discretized matrix propagators and self-energies will be denoted in bold type. The inverse free level Green function discretized on the contour is then given by [79] h h h h g i where h ± =1miò 0 Δt, Δt indicates the time step in the discretization with N=t/Δt. In this expression the initial level charge is determined by n σ (0)=ρ σ /(1+ρ σ ). Note that the discretization over the contour is made starting from t=0 to the final time through the positive Keldysh branch and returning to t=0 through the negative one. The time-dependent tunneling self-energies can be evaluated straightforwardly and at zero temperature have the simple form [44] t t t t t t t t , e e , , e e , 4 2 being the leads bandwidth. Alternatively, it is possible to take the limit D  ¥ provided that a finite temperature, taken as the smallest energy parameter, is introduced (see [44]). In all the results given below we consider this infinite bandwidth limit except when comparing with numerically exact methods where an energy cutoff with a precise value is used. The other Keldysh self-energy components are then given by T  T   T  T  T   , , , where θ(t) is the Heaviside step function. Notice that there is an ambiguity in the definition of these self-energies at equal times. It turns out that the different possible choices in the definition of t t , T, S s ++ ( )and t t , T, We have checked that this choice appropriately recovers the correct stationary limit and perfectly reproduces the transient behavior in the cases where an analytic expression is available (see section 3.1). The evaluation of the interaction self-energy will be discussed in sections 3-5 for the cases of electronelectron and electron-phonon interactions. For computing the correlation part of the interaction self-energy we also find that the most stable algorithm consists on the calculation of the non-diagonal Keldysh components ( t t , int S ¢ +-( )and t t , int S ¢ -+ ( )) and then using the relations of equations (5), (6) for the diagonal ones. The self-energies are then evaluated in the discrete time mesh ( figure 1(a)). The propagators in Keldysh space can now be obtained by numerically inverting the matrix Notice the factor (Δt) 2 introduced by the discretization procedure. The knowledge of G t t , ¢ ŝ ( ) enable us to calculate the evolution with time of the electronic and transport properties of the system such as the central level charge, the spectral density and the current. Thus, the level charge can be calculated as n t while the current through the interface between the central region an the electrodes is given by Finally, following [54,80], it is possible to define a time dependent auxiliary spectral density function per spin A σ (ω, t) by calculating the current to weakly coupled probes and which tends to the correct stationary value . For the present system we have

Electron-electron interaction: the Anderson model
In this section we will consider the Anderson model and U is the local Coulomb repulsion.

Hartree-Fock (HF) solution
The dynamical HF solution of the Anderson model provides an ideal test of the accuracy of the numerical method presented in section 2 as in this case the time-dependent problem can be exactly solved [2,3]. Thus, within this approximation, the model becomes an effective single electron problem with a spin and timedependent central level where n σ (t) is the central level occupation per spin. As commented in the previous section, the problem of an impurity level in a time-dependent potential coupled to metallic leads is exactly solvable using the Keldysh method. For the HF case addressed in this paper, the Keldysh Green function can be written in a very compact way as The time evolution of the central level occupation is then obtained as n t G t t i , ,HF = s s +-( ) ( ) and has the form One can compare the result of the numerical method proposed in section 2 with equation (13). In the HF approximation the self-energy is given by the left diagram of figure 1 (b) and has the form HF, where α, β=± are the Keldysh branch indexes. Notice that the Dirac delta in the previous equation is converted to a Kronecker δ function, including an additional 1/Δt factor when discretizing in the time mesh. We can now obtain the propagators in the HF approximation by inverting and following the numerical procedure presented in the previous section. The dynamical properties of the system can be now calculated from G HF,ŝ . It is worth remarking that the self-consistency condition on the charge in this approximation is particularly straightforward as it is simply achieved by storing the charge values obtained in the discrete mesh by inverting equation (15) at each time step, starting from the initial one n σ (0). The undefined components of the self-energy at each final time can be accurately approximated as the self-energy one time step before i.e. t t t t , , ). The error introduced by this approximation becomes negligible for a sufficiently small Δt. In the finite bandwidth situation this means Δt  1/D and in the wide band limit Δt has to be taken smaller than the inverse of the greatest energy scale. It can be checked that this procedure leads to the proper stationary values of n σ (t) in the unrestricted self-consistent HF approximation.
In figure 2 we show the time evolution of the central level charge per spin at different levels of approximation. In figure 2(a) we compare the exact MC results from [37] with the ones obtained for the selfconsistent and the non self-consistent (first order) HF approximation for the electron-hole symmetric situation (ò 0 =−U/2) and the n n 0 , 0 0, 0 =   ( ( ) ( )) ( ) initial configuration. As can be observed, the non self-consistent approximation tends to deviate from the exact results, leading to a stationary charge overpassing the electronhole symmetric stationary value. This result is in agreement with [37], where the authors analyzed the level population by means of a first order perturbation theory in the Coulomb interaction parameter U/Γ. Although a good agreement is found for small values of U/Γ, the charge progressively deviates from the exact results for increasing U/Γ. This pathological behavior is corrected within a fully self-consistent HF treatment, where the average charge per spin n σ (t) tends to the correct singlet state for all U/Γ values. As shown below, inclusion of correlations provided by the second order diagrams further improve the agreement with the numerically exact results.
In figure 2(b) we show the level population evolution for an initially trapped spin, n n 0 , 0 We have chosen a case with electron-hole symmetry (ò 0 =−U/2) and with parameters such that U/πΓ>1, which leads to a magnetic solution in the stationary case within the HF approximation [81]. As can be observed, the numerical solution is in remarkable agreement with the exact expression of equation (13). Let us comment that for initial conditions with unbroken spin symmetry, i.e. n n 0 , 0 0, 0 , 1, 1 =   ( ( ) ( )) ( ) ( ), the system always tends to a non-magnetic solution for all values of U/Γ.
Finally, it is worth remarking that the prediction of a magnetic solution within the HF approximation at zero-temperature is well known to be unphysical as the ground state of the system should be always a singlet [82][83][84]. This behavior should be corrected when including electronic correlations in an appropriate way. In the next section we will analyze the effect of correlations beyond the HF approximation in the transient regime.  (13) for U/Γ=4, ò 0 =−U/2 and in the infinite bandwidth limit. The arrows denote the stationary solution. The red lines correspond to the second-order self-energy case. In (c) the level charge for the same parameters and three different Coulomb interactions U/Γ=4 (red), 6 (green) and 8 (black) is shown for the second order approximation. In a continuous line we show the evolution of the spin up and in dashed the evolution of the down spin.

Effects of correlation beyond the Hartree-Fock approximation
Within a Green functions approach, correlation effects are included in the electron self-energy. In a stationary situation an appropriate second-order self-energy in the interaction parameter U/Γ can include these effects in a rather satisfactory way. Indeed it can be shown that the exact self-energy in the limit U G  ¥ has a functional form close to the second order one and is in fact proportional to U 2 [85]. This fact has been used in different interpolative approaches based on the second-order self-energy giving a reasonable approximation for the Anderson model between the weak and strong coupling limits [85][86][87][88].
We will concentrate in the symmetric case, ò 0 =−U/2, where correlations effects are more important. It can be shown that the inclusion of the second-order self-energy yields a spectral density in the equilibrium stationary case in rather good agreement with numerical renormalization group (NRG) calculations [89]. Indeed in this approximation the charge peaks in the spectral density are well described, fulfilling the Friedel sum rule at zero energy, although somewhat overestimating the width of the Kondo resonance at very large U/Γ. It is important to notice that the second-order self-energy diagram has to be calculated with propagators including the HF correction to the energy level (i.e. the HF propagators) in order to ensure electron-hole symmetry. On the other hand, it can be shown that a fully self-consistent calculation of the diagrams (i.e. using fully dressed propagators) yields instead a poor description of the spectral density [90].
In a general time-dependent non-equilibrium situation the self-energy diagrams must be calculated in Keldysh space. The +-(-+) components of the second order self-energy have the simple expressions here the HF propagators are calculated as indicated in equation (15). The other Keldysh components are then given by the usual Keldysh relations, making the same choice for equal times as in equation (6). The propagators in Keldysh space can now be evaluated inverting equation (7) with int, ) .
We will now analyze the effect of correlations on the electronic and transport properties of the system. In figure 2(a) we show the population evolution for the case discussed in the previous section and an initial configuration n n 0 , 0 0, 0 =   ( ( ) ( )) ( ). As can be observed, the inclusion of electron correlation effects improve the agreement with the exact MC calculations.
In figure 2(b) we show the evolution of n σ (t) with an initial configuration n n 0 , 0 1, 0 =   ( ( ) ( )) ( ) in which a magnetic solution was predicted by the HF approximation. As it can be observed, when including correlations (electron-hole pair creation) the system evolves to a non-magnetic solution corresponding to a singlet state in the stationary limit. In figure 2(c) we analyze the evolution of n σ (t) for the same initial magnetic configuration for increasing values of the electron-electron interaction parameter. It is found that for U 8  G the initial localized spin is no longer screened by the electrodes, tending to a magnetic solution. This indicates a shortcoming of the approximate self-energy for sufficiently large interaction strength. The singlet stationary state is, however, always reached when starting from a configuration without spin-symmetry breaking.
In figure 3(a) we analyze now the long time evolution of the DOS, A t , w  ¥ ( ). These results demonstrate that the second order self-energy provides a good approximation to the problem [91], leading to a remarkable agreement with results from NRG calculations for moderate U/Γ values [89]. The inset in this panel shows a blow up of the Kondo resonance, where it can be observed that the second order self-energy tends to overestimate its width for large U/Γ values.
It should be remarked that the convergence time increases with U/Γ. In this respect it is interesting to analyze the convergence in time of the Kondo resonance, an issue that has been addressed in previous works [63,92]. One would expect this convergence time to be of the order of T T , ] , for these cases we have the ratio T K (U/Γ=4)/T K (U/Γ=8) ; 3.4. Thus, one would expect a Kondo resonance formation time for the U/Γ=8 case roughly 3.4 times larger than for U/Γ=4. The ratio of formation times that can be estimated from figures 3(b) and (c) is somewhat smaller due to the slight overestimation of the width of the Kondo peak by the second order diagrammatic approximation for the larger U/Γ value. On the other hand, figure 3(d) shows the evolution of the height of the central peak, A(ω=0, t), to its stationary value fixed by the Friedel sum rule A t 0, . A kink in the evolution is observed at times ∼1/U, mainly visible for large U/Γ values, due to the appearance of the charge bands.
Let us discuss now the voltage biased situation. In figure 4(a) we show the evolution of the current for the second order perturbation expansion together with results from the MC simulations finding also a good quantitative agreement. For very large interaction strengths the agreement becomes somewhat poorer although still capturing the general trend.
Finally in figure 4(b) we show the asymptotic I(V ) characteristic for increasing U/Γ values compared to the MC results of [93]. As can be observed, there is an overall good agreement specially for V>Γ. However, for V<Γ the second order self-energy tends to slightly overestimate the current due to the already mentioned shortcoming in the description of the Kondo resonance. In fact, this approximation is unable to describe the splitting of this resonance for V<T K . This shortcoming would be removed in this electron-hole symmetric case by including the fourth order diagrams, as shown in [94] in the stationary limit.

Electron-phonon interaction: spinless Anderson-Holstein model
In order to analyze the transient regime in the presence of electron-phonon interactions we will first consider the spinless Anderson-Holstein model [95]. In this model an electron in the central level is coupled to a single vibrational mode. The Hamiltonian of the system is given by

Hartree solution
As in the previous section, we begin our analysis with the self-consistent mean-field approximation in which the self-energy is approximated by the 'tadpole' diagram of figure 5 (Hartree approximation). Within this approximation, the self-energy in Keldysh space can be evaluated as where n(t) is the self-consistent central level charge and dˆis the free phonon propagator in Keldysh space given by where n e 1 ) is the free phonon population, described in a thermal equilibrium situation by the Bose-Einstein distribution. Most of the calculations are performed at zero or very small temperature, considering n p =0. Using the Keldysh relations, equations (18) can then be written as It is worth noticing that, at variance with the case of the electron-electron interaction discussed in the previous section, the electron-phonon interaction introduces retardation effects even in the Hartree approximation. These effects will be important in the transient regime except in the limit of a sufficiently fast phonon (ω 0 ? ò 0 , Γ) [33] with a central charge evolving adiabatically. In this limit equation (20) tends to We can now follow a similar procedure to the one used in the previous section to calculate G Ĥ and the central level self-consistent charge. Figures 6(a) and (b) show the evolution of the level charge in the transient regime. As in the case of electron-electron interactions, the charge evolves to the stationary value, indicated by the arrows in the figure. Figures 6(a) and (b) also illustrate how the solution progressively deviates from the adiabatic approximation given by equation (22) when reducing the value of ω 0 . The full self-consistent solution as given by the self-energy in equation (20), describes the time-dependent modification of the central level charge at time t induced by its past history at time τ<t. Retardation effects of the phonon dynamics results in a coherent superposition of oscillations with period 2π/ω 0 but with different amplitudes (∝n(τ)). In the intermediate regime where the electron and the phonon dynamics are equally fast (ω 0 ≈Γ), the coherence between those oscillations is lost at long times (t ? 1/Γ,2π/ω 0 ), thus resulting in an effective damping of the central level charge, see figures 6(a) and (b). However, in the adiabatic regime (ω 0 ? Γ) the dynamics of the electrons and phonons decouple, and small charge oscillations persist on time, mostly in the n(0)=1 case (black curve in figure 6(a)). A natural lifetime describing the decay of those oscillations could be included by dressing the phonon line in the Hartree term depicted in figure 5(a).
Finally, one can observe in figures 6(a) and (b) that for the smallest values of ω 0 two different asymptotic charge values are reached depending on the initial level population. This is the charge bistable behavior predicted by the self-consistent Hartree approximation in the strong-coupling limit [33,97]. For the case of electron-hole symmetry considered in figure 6 and at zero temperature and bias voltage, the condition for the appearance of bistability is 2λ 2 /πΓω 0 >1. The possibility of a bistable regime for a molecular quantum dot with strong electron-phonon interaction was suggested some time ago [98][99][100]. The interest in investigating such a Figure 5. Second order self-energy diagrams for the spinless Anderson-Holstein interaction. (a) Diagrams for the Born approximation, using the bare phonon propagator (wavy line). In (b) similar approximations are shown with two different schemes for dressing the phonon propagator: RPA [96], where the electronic propagators are considered to be undressed, and the selfconsistent MIGDAL [27], where the electronic propagators are fully dressed. phenomenon has experienced a recent revival. For instance, it has been shown that the displacement fluctuation spectrum of a nanomechanical oscillator strongly coupled to electronic transport, either in the regime of semiclassical phonons [101,102], or for a quantum nanomechanical oscillator entering the Franck-Condon regime [103] bears clear signatures of a transition to a bistable regime. Moreover, by making a mapping to the Kondo problem, the bistability was shown to be destroyed in equilibrium conditions by quantum fluctuations if the temperature is lower than a phonon mediated Kondo temperature [53]. Notice, however, that this phonon displacement bistability does not correspond necessarily to a bistable behavior for the charge or the current, as predicted by the mean field approximation. As even this simple spinless Anderson-Holstein model is not exactly solvable, this issue is still under debate [49,104,105]. It seems to us plausible that, at least for equilibrium conditions and T=0 correlation effects destroy the charge and current bistability predicted by the mean field solution. We address this issue in the following section.

Effects of correlation beyond Hartree approximation
We will go beyond the mean-field solution by analyzing three different approximations for the self-energy. We first consider the self-consistent Born approximation given by the diagrams in figure 5(a). This is a conserving approximation in which the diagrams are calculated from the fully dressed electron propagators. The phonon propagator is however not renormalized. Within this approximation both diagrams appearing in figure 5(a) have the expression This fully self-consistent approximation can be straightforwardly implemented within the numerical procedure of section 2. For each time in the discretized mesh, the self-energies of equations (23) are calculated from the final Green functions and then stored. As in the case of the Hartree solution previously discussed, when inverting equation (7) for each time the self-energies at the final time in each iteration are not well defined but its value can be extrapolated from the ones calculated at the previous mesh point in the time grid. For sufficiently small Δt the error introduced by this approximation becomes negligible. We have checked the accuracy of this procedure by verifying that the solution tends to the proper stationary one.
In figure 6(c) we show the evolution of the central level charge for a choice of parameters in which the Hartree approximation predicts a bistable behavior. As can be observed, the inclusion of correlations eliminates the charge bistability appearing in the Hartree approximation, tending to the same asymptotic value for the initially empty and full cases. We have checked that this behavior is maintained up to quite large values of λ 2 /ω 0 Γ, although eventually the self-consistent Born approximation breaks down in the strong polaronic regime. This indicates that another kind of approximation has to be used to explore this parameter regime, like for instance in the lines of the ones discussed in [105][106][107][108][109]. These results suggest that the bistable behavior of the . The dotted lines represent the evolution using an instantaneous Hartree term equation (22), while the solid ones correspond to the full Hartree self-energy equation (20). The dependence on phonon frequency is also illustrated for the values: ω 0 =8Γ (black), 2Γ(green) and Γ (blue). (c) Charge evolution for an initially empty (dashed line) and initially full level (solid line) for the ω 0 =Γ case. The blue and red lines correspond to the Hartree and self-consistent Born approximation, respectively. The remaining parameters are λ=1.5Γ, V=0 and the central level is set to ò 0 =λ 2 /ω 0 , thus preserving electron-hole symmetry.
central level charge predicted in [33,97] is a spurious feature of the mean field approximation which disappears when correlation effects are included. This is in agreement with the predictions of exact numerical calculations of [53], at least for the equilibrium case and at sufficiently low temperatures. It does not imply that an apparent bistability might not be observed for a continuous bath model or adopting more general initial conditions for the phonon mode density matrix [110,111].
So far, the renormalization of the phonon propagator has been neglected. The simplest way to include this effect is by means of an RPA-like approximation [96]. The phonon propagator will satisfy a Dyson equation in Keldysh space similar to the electronic one given in equation is the inverse free phonon propagator and P is the phonon self-energy given by i , , .
As in the electronic case, equation (24) can be discretized in a time mesh along the Keldysh contour. In order to solve numerically the corresponding matrix equation, an expression for the inverse free phonon propagator discretized on the contour must be obtained. This is a task which, to best of our knowledge, has not been achieved in the literature, the mathematical difficulty lying in the fact that the inverse phonon propagator becomes singular in the free limit. This singularity must be then somehow regularized. To obtain an expression of d  . The parameter η can be interpreted as a small phonon relaxation rate which has to be taken such as 1/η ? t,1/ω 0 for a good convergence to the expected free propagator when inverting equation (26).
It should be noticed that this problem with the inversion of the free phonon propagator has been avoided in the literature by neglecting fast oscillating terms of the type in the diagrammatic expansion of D . This corresponds to the so-called rotating wave approximation, which describes the regime where the phonon timescale is much faster than the electron dynamics ( , ) [112,113]. For the calculation of the phonon self-energy, P , we will analyze two different approximations. In the first one (that will be denoted as RPA) the propagators in the electron 'bubble' are the non-interacting ones, whereas the fully dressed propagators will be used in the second one (denoted as MIGDAL), see figure 5(b).
In figure 7 we show the long-time DOS at the central level for the three approximations considered in this section using the same parameters as in figure 6 with ω 0 =2Γ; a case with a rather strong electron-phonon coupling although still far from the polaronic limit (λ 2 /(ω 0 Γ) ? 1). Notice the dip in the DOS at ω≈ω 0 in the self-consistent Born approximation, which is a feature due to the logarithmic divergence of the second order self-energy X w Ŝ ( )at ω=ω 0 [106,114]. As it can be observed, both RPA and MIGDAL approximations, which include phonon renormalization eliminate this pathological divergence. A slight shift of the resonance around ω 0 due to the renormalization of the phonon mode in both approximations can be observed. Notice also that all these approximations lead to an additional feature at ∼2ω 0 , associated to the appearance of a second phonon sideband. As an additional remark, in all cases the zero energy spectral density tends to reach the expected value predicted by the Friedel sum rule [115].
A further check of these approximations can be made by comparing their long-time DOS with the one predicted by a NRG calculation. To this end we have performed a NRG calculation of the stationary DOS for the parameters of figure 7. As can be observed the agreement with the results of both RPA and MIGDAL is quite good for this parameter range. It should be commented that neither of these approximations are expected to be valid in the strong polaronic limit. Thus, features like the exponential decrease of the central resonance together with the appearance of a multiphonon structure in the DOS [105][106][107][108][109] would require an approximation for the self-energy valid in the polaronic regime, as commented above.
Finally, in figure 8 we show results from the three approximations for the transient left, right and average currents compared to results obtained using MC simulations in [48]. Both cases correspond to a rather strong interaction (λ=8, ω 0 =10 and Γ=1) but two different bias voltages. Strikingly, as can be observed, RPA captures remarkably well the quantitative behavior of the numerically exact results in the small voltage case, see figures 8(a)-(c), whereas for very large bias it is the MIGDAL approximation that gives a better quantitative agreement with the MC numerical results, see figures 8(d)-(f). This would indicate that the inclusion of phonon renormalization and non-equilibrium effects (like heating of the local vibrational mode under increasing bias voltage) are essential for a good description of this regime. Furthermore, the higher the bias voltage the better this effects are included in the fully self-consistent approach given by MIGDAL.

Electron-electron and electron-phonon interactions
In this section we study the transient regime in the presence of both electron-electron and electron-phonon interactions. We consider the spin-degenerate Anderson-Holstein model defined as   figure 9(a) we show the long time spectral density compared to the exact NRG results from [116] using the RPA for e ph S -. Similar results are obtain for the MIGDAL approximation. As can be observed, for the smaller λ case the RPA exhibits an overall agreement with the exact results. However, for larger values of the electronphonon interaction the agreement becomes poorer (blue curve). In fact, this diagrammatic self-consistent approximations would not describe properly the transition to an insulating phase which is expected when increasing the electron-phonon interaction for U 2 2 0  l w [106,116,117]. To explore this parameter regime, one would need to develop an approximation correctly interpolating between the perturbative regime and the strong polaronic limit.
Finally, in figures 9(b) and (c) we show the time evolution of the spectral density for λ/Γ=0 in figure 9(b) and λ/Γ=2 in figure 9(c), with U/Γ=8 for the RPA. We show that, even in the Kondo dominated regime, the electron-phonon interaction modifies significantly the system dynamics, leading to longer convergence times. This is illustrated in c) where the height of the central resonance, A(ω=0, t), is represented. We show that, although the central resonance width in the long time regime is not significantly modified with respect to the pure Kondo case, it exhibits different dynamical properties like oscillations with a period ∼2π/ω 0 . Furthermore, the decay time of these oscillations is considerably longer with respect to the U=0 case (not shown), indicating that the electron-electron interaction increases phonon retardation effects.

Conclusions
We have presented an accurate and stable algorithm to calculate the transient transport properties of interacting nanojunctions. We have shown how different self-consistent diagrammatic approximations can be implemented within this framework, yielding accurate results for both the transient and the steady state regimes. The method has allowed us to address several issues of great current interest in the condensed matter community like the dynamical build up of Kondo correlations and the possible existence of bistability in the presence of strong electron-phonon interactions.
For the Anderson model we have analyzed the evolution of the spectral density explicitly exhibiting the formation of the Kondo resonance. In both cases of zero and finite voltage bias, the results are in good agreement with available numerically exact calculations. For the electron-phonon case we have implemented two different schemes for dressing the phonon propagator (denoted as RPA and MIGDAL), showing the importance of a good description of the phonon dynamics to obtain accurate results. As a technical requirement for this implementation we have derived an expression for the inverse of the time-discretized Keldysh free phonon propagator, allowing us to go beyond previous approaches to the problem based on a rotating-wave like approximation. Comparison with numerically exact results shows that the RPA and the MIGDAL approximation can provide accurate results for the transient currents up to rather strong coupling values in the  [116], for two different electron-phonon coupling parameters: λ=1.89 (red) and λ=3.14 (blue). The remaining parameters are U=6.3, ω 0 =3.14, ò 0 =λ 2 /ω 0 −U/2 and Γ R =Γ L =0.5. (b) and (c) time evolution of the density of states for λ=0 and λ=2. In (d) we represent the central peak height evolution, showing in red the λ=0 and in blue the λ=2 case. The remaining parameters are U=8, ω 0 =2, ò 0 =λ 2 /ω 0 −U/2, Γ=1 and V=0. low and high voltage regimes, respectively. Regarding the possible bistable behavior, we have found that electron correlation effects beyond the mean-field approximation tend to suppress its appearance, in agreement with recent numerically exact results [53]. However, this does not imply that upon choosing a different initial condition for the vibron density matrix in a model including low frequency modes, one should not observe an apparent bistability, as indicated in [104,110].
Finally, we have analyzed the situation where both interactions are present showing a reasonable agreement with the available numerically exact results for moderate electron-phonon coupling. We have also shown that the presence of electron-phonon interactions in the Kondo dominated regime introduces additional dynamical features in the evolution of this resonance. We notice, however, that addressing the strong polaronic limit would require the implementation, within the present framework, of non-perturbative approximations for the selfenergy in the spirit of [107][108][109].