# Transient dynamics in interacting nanojunctions within self-consistent perturbation theory

R Seoane Souto, Rémi Avriller, a Levy Yeyati, A Martín-Rodero

► To cite this version:

R Seoane Souto, Rémi Avriller, a Levy Yeyati, A Martín-Rodero. Transient dynamics in interacting nanojunctions within self-consistent perturbation theory. New Journal of Physics, 2018, 20 (8), pp.083039 (1-16). 10.1088/1367-2630/aad99d. hal-02014525

HAL Id: hal-02014525

https://hal.science/hal-02014525

Submitted on 13 Feb 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Transient dynamics in interacting nanojunctions within self-consistent perturbation theory

To cite this article: R Seoane Souto et al 2018 *New J. Phys.* **20** 083039

View the [article online](#) for updates and enhancements.

## Related content

- Full-counting statistics of energy transport of molecular junctions in the polaronic regime

Gaomin Tang, Zhizhou Yu and Jian Wang

- Functional renormalization group study of the Anderson-Holstein model

M A Laakso, D M Kennes, S G Jakobs et al.

- NRG study of the Anderson-Holstein model

A C Hewson and D Meyer

<figure><img src="image_2.png" /></figure>

**IOP** | **ebooks™**

Bringing you innovative digital publishing with leading voices to create your essential collection of books in STEM research.

Start exploring the collection - download the first chapter of every title for free.
<figure><img src="image_3.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

20 March 2018

REVISED

8 August 2018

ACCEPTED FOR PUBLICATION

10 August 2018

PUBLISHED

28 August 2018

Original content from this
work may be used under
the terms of the *Creative
Commons Attribution 3.0
licence*.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_4.png" /></figure>

PAPER

Transient dynamics in interacting nanojunctions within self-
consistent perturbation theory

R Seoane Souto¹, , R Avriller², A Levy Yeyati¹ and A Martín-Rodero¹

¹ Departamento de Física Teórica de la Materia Condensada, Condensed Matter Physics Center (IFIMAC) and Instituto Nicolás Cabrera, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

² Univ. Bordeaux, CNRS, LOMA, UMR 5798, F-33405 Talence, France

**Keywords:** nanoelectronics, molecular electronics, interactions at the nanoscale, non-equilibrium phenomena, dynamical properties, Kondo effect

Abstract

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.

**1. 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–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].

On the theoretical side transport transient dynamics has been addressed using different methods valid for different regimes. Thus, the scattering approach or the non-equilibrium Green function formalism have been used for describing the dynamics in the coherent non-interacting regime [19–30]. However, the inclusion of interactions is essential to analyze the transport dynamics through localized states, as is the case of molecular junctions or semiconducting quantum dots. For these cases, rate equations approaches, adequate for the sequential tunneling regime, have been extensively used [31, 32]. The most interesting and general coherent-interacting regime constitutes a great theoretical challenge. This regime has been addressed using several complementary approaches: diagrammatic techniques [33–47], quantum Monte-Carlo (MC) [48–55], time-dependent NRG [56–63], time-dependent DFT [64–70] among others [71–75]. However, all of these techniques as they are actually implemented have some limitations. For instance, numerically exact methods like quantum MC are strongly time-consuming, require finite temperature and typically do not allow to reach long time scales. Similar concerns can be applied to the case of time-dependent NRG.
This situation suggests the convenience of revisiting perturbative diagrammatic methods for analyzing transport transient dynamics in interacting nanoscale devices. Although these methods have been partially explored in previous works [33, 37], these implementations did not, in general, include self-consistency which can become of essence in order to increase the accuracy and range of validity of these methods. Moreover, in the case of models including electron–phonon interactions further methodological developments are needed in order to take into account properly the dynamical build up of a non-equilibrium phonon distribution.

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 electron–phonon 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.

## 2. Keldysh formalism for the transient regime

For describing a nanoscale central region coupled to metallic electrodes we consider a model Hamiltonian of the form $\hat{H} = \hat{H}_{\text{leads}} + \hat{H}_{c} + \hat{H}_{T} + \hat{H}_{\text{int}}$, where

$$ \hat{H}_{\text{leads}} = \sum_{k\sigma,\nu} \epsilon_{k\sigma,\nu} c_{k\sigma,\nu}^{\dagger} c_{k\sigma,\nu}, \quad \hat{H}_{c} = \sum_{\sigma} \epsilon_{0} c_{0\sigma}^{\dagger} c_{0\sigma}, \quad \hat{H}_{T} = \sum_{k\sigma,\nu} [v_{k\sigma,\nu}(t) c_{k\sigma,\nu}^{\dagger} c_{0\sigma} + \text{h.c.}], \quad (1) $$

where $c_{k\sigma,\nu}$ with $\nu = L(R)$ labeling the left (right) electrode, and $c_{0\sigma}$ are annihilation operators for electrons in the leads and in the central region respectively and $v_{k\sigma,\nu}(t)$ is the tunneling amplitude which will depend on time. The two electrodes can be kept at different chemical potentials via a constant bias voltage $eV = \mu_{L} - \mu_{R}$. For simplicity the central region will consist of a single quantum level denoted by $\epsilon_{0}$. The last term, $\hat{H}_{\text{int}}$, in $\hat{H}$ describes the many body interactions in the central region, which we shall specify later. Hereafter we assume $e = \hbar = k_{B} = 1$.

In what follows we will consider the wide-band approximation for the electrodes. Within this approximation the tunneling rates can be taken as constants, $\Gamma_{\nu} = \pi \sum_{k} |v_{k\sigma,\nu}|^2 \delta(w - \epsilon_{k\sigma,\nu}) \sim \pi |v|^2 \rho_{\text{F}}$, where $\rho_{\text{F}}$ is the density of states at the Fermi edge, the resonant level width being $\Gamma = \Gamma_{L} + \Gamma_{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\sigma,\nu}(t) = \theta(t)v_{k\sigma,\nu}$, which allows us to define a time-dependent tunneling rate $\Gamma(t) = \theta(t)\Gamma$. 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 Keldysh space, $\hat{G}_{\sigma}(t, t') = -i \langle \hat{T}_{c} c_{0\sigma}(t) c_{0\sigma}^{\dagger}(t') \rangle$, where $\hat{T}_{c}$ 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 self-energy. This is not always a simple task in the dynamical problem. On the other hand, even if an appropriate self-energy 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
<figure><img src="image_5.png" /><figcaption>Figure 1. (a) Keldysh contour considered to analyze the transient regime, Δt being the time step in the discretized calculation of the time dependent Green function. (b) Self-energy diagrams for the Coulomb interaction up to second order. The solid line represents the electron propagator and the wavy line the interaction. In the HF diagram the double line indicates the charge calculated using the dressed propagators.</figcaption></figure>

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

$$ \hat{G}_{\sigma} = [\hat{g}_{\sigma}^{-1} - \hat{\Sigma}_{\sigma,T} - \hat{\Sigma}_{\sigma,\text{int}}]^{-1}, \quad (2) $$

where $\hat{g}_{\sigma}^{-1}$ is the inverse free electron propagator of the uncoupled central level, $\hat{\Sigma}_{\sigma,T}$ the tunneling self-energy and $\hat{\Sigma}_{\sigma,\text{int}}$ 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]

$$ ig_{\sigma}^{-1} = \begin{pmatrix} -1 & & & & & -\rho_{\sigma} \\ h_{-} & -1 & & & & \\ & h_{-} & -1 & & & \\ & & \ddots & \ddots & & \\ & & & h_{+} & -1 & \\ & & & & h_{+} & -1 \end{pmatrix}_{2N \times 2N}, \quad (3) $$

where $h_{\pm} = 1 \mp i\epsilon_0\Delta t$, $\Delta t$ indicates the time step in the discretization with $N = t/\Delta t$. In this expression the initial level charge is determined by $n_{\sigma}(0) = \rho_{\sigma}/(1 + \rho_{\sigma})$. 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]

$$ \Sigma_{T,\sigma}^{+-}(t, t') = \frac{\Gamma}{\pi} \sum_{\nu} \frac{e^{-i\mu_{\nu}(t-t')} - e^{iD(t-t')}}{(t-t')}, \quad \Sigma_{T,\sigma}^{-+}(t, t') = \frac{\Gamma}{\pi} \sum_{\nu} \frac{e^{-i\mu_{\nu}(t-t')} - e^{-iD(t-t')}}{(t-t')}, \quad (4) $$

2D being the leads bandwidth. Alternatively, it is possible to take the limit $D \to \infty$ 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

$$ \Sigma_{T,\sigma}^{++}(t, t') = -\theta(t-t')\Sigma_{T,\sigma}^{+}(t, t') - \theta(t'-t)\Sigma_{T,\sigma}^{-+}(t, t'), \\ \Sigma_{T,\sigma}^{-+}(t, t') = -\theta(t-t')\Sigma_{T,\sigma}^{-+}(t, t') - \theta(t'-t)\Sigma_{T,\sigma}^{+}(t, t'), \quad (5) $$

where $\theta(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 $\Sigma_{T,\sigma}^{+}(t, t)$ and $\Sigma_{T,\sigma}^{-}(t, t)$ can significantly affect the convergence and stability of the system properties with time. Although the precise value of $\Sigma_{T,\sigma}^{+}(t, t)$ and $\Sigma_{T,\sigma}^{-}(t, t)$ depends on the whole energy range of the leads density of states, if one is not interested in the dynamics on time scales smaller than $1/D$ there is freedom to choose this value. We have found that the most stable algorithm corresponds to the choice

$$ \Sigma_{T,\sigma}^{++}(t, t) = \Sigma_{T,\sigma}^{-+}(t, t) = -\frac{\Sigma_{T,\sigma}^{+}(t, t) + \Sigma_{T,\sigma}^{-}(t, t)}{2}. \quad (6) $$

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 electron–electron 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 ($\hat{\Sigma}_{\text{int}}^{+}(t, t')$ and $\hat{\Sigma}_{\text{int}}^{-}(t, t')$) 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

$$ \hat{G}_{\sigma}^{-1} = \hat{g}_{\sigma}^{-1} - (\Delta t)^2 (\hat{\Sigma}_{T,\sigma} + \hat{\Sigma}_{\text{int},\sigma}). \quad (7) $$

Notice the factor $(\Delta t)^2$ introduced by the discretization procedure.

The knowledge of $\hat{G}_{\sigma}(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_{\sigma}(t) = iG_{\sigma}^{+}(t, t)$, while the current through the interface between the central region and the electrodes is given by

$$ I_{\nu} = \sum_{\sigma} \int_{0}^{t} [G_{\sigma}^{+}(t, t_{1}) \Sigma_{T,\sigma\nu}^{-}(t_{1}, t) - G_{\sigma}^{-}(t, t_{1}) \Sigma_{T,\sigma\nu}^{+}(t_{1}, t)] dt_{1}. \quad (8) $$

Finally, following [54, 80], it is possible to define a time dependent auxiliary spectral density function per spin $A_{\sigma}(\omega, t)$ by calculating the current to weakly coupled probes and which tends to the correct stationary value at large times $A_{\sigma}(\omega) = \text{Im}[G_{\sigma}^{A}(\omega) - G_{\sigma}^{R}(\omega)]/2\pi$. For the present system we have

$$ A_{\sigma}(\omega, t) = \text{Im} \int_{0}^{t} dt' \frac{e^{-i\omega(t-t')}}{2\pi} [G_{\sigma}^{+}(t', t) - G_{\sigma}^{-}(t', t)], \quad (9) $$

and the spin averaged spectral density as $A(\omega, t) = \sum_{\sigma} A_{\sigma}(\omega, t)/2$.

### 3. Electron–electron interaction: the Anderson model

In this section we will consider the Anderson model [81] consisting of a single spin degenerate level with on-site electron–electron repulsion, coupled to metallic electrodes. The interaction term in the Hamiltonian of section 2 is given by $\hat{H}_{e-e} = U\hat{n}_{1}\hat{n}_{1}$, where $\hat{n}_{\sigma} = c_{0\sigma}^{\dagger}c_{0\sigma}$ and $U$ is the local Coulomb repulsion.

#### 3.1. 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 time-dependent central level

$$ \epsilon_{\sigma}(t) = \epsilon_{0} + U n_{\bar{\sigma}}(t), \quad (10) $$

where $n_{\sigma}(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

$$ G_{\text{HF},\sigma}^{+}(t, t') = \theta(t)\theta(t')ie^{-i[t\bar{\epsilon}(t)-t'\bar{\epsilon}(t')]}e^{-\Gamma(t+t')} \left\{ n_{\sigma}(0) + \frac{1}{\pi} \int d\omega \left[ \sum_{\nu=L,R} \Gamma_{\nu} f_{\nu}(\omega) \right] g_{\sigma}(\omega, t) g_{\sigma}^{*}(\omega, t') \right\}, \quad (11) $$

where

$$ \bar{\epsilon}_{\sigma}(t) = \frac{1}{t} \int_{0}^{t} d\tau \epsilon_{\sigma}(\tau), \quad g_{\sigma}(\omega, t) = \int_{0}^{t} d\tau e^{-i[\omega + i\Gamma - \bar{\epsilon}_{\sigma}(\tau)]\tau}. \quad (12) $$

The time evolution of the central level occupation is then obtained as $n_{\sigma}(t) = iG_{\sigma,\text{HF}}^{+}(t, t)$ and has the form

$$ n_{\sigma}(t) = e^{-2\Gamma t} \left\{ n_{\sigma}(0) + \int \frac{d\omega}{\pi} \left[ \sum_{\nu=L,R} \Gamma_{\nu} f_{\nu}(\omega) \right] |g_{\sigma}(\omega, t)|^2 \right\}. \quad (13) $$

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

$$ \sum_{\text{HF},\sigma}^{\alpha\beta}(t, t') = \alpha U n_{\bar{\sigma}}(t)\delta(t-t')\delta_{\alpha\beta}, \quad (14) $$

where $\alpha, \beta = \pm$ are the Keldysh branch indexes. Notice that the Dirac delta in the previous equation is converted to a Kronecker $\delta$ function, including an additional $1/\Delta t$ factor when discretizing in the time mesh. We can now obtain the propagators in the HF approximation by inverting the matrix
<figure><img src="image_6.png" /><figcaption>Figure 2. Time evolution of the level charge for the Anderson model. (a) Average population per spin at different levels of approximation compared to MC simulations from [37], with $U/\Gamma = 8$, $V = 0$, $\epsilon_0 = -U/2$, $T = 0.2\Gamma$, $D = 10\Gamma$ and $(n_\uparrow(0), n_\downarrow(0)) = (0, 0)$ initial configuration. (b) and (c) results for the up (solid lines) and down (dashed lines) spin and the initial configuration (1,0). In (b) the HF approximation (blue lines) is compared to the analytic expression (black points), given by equation (13) for $U/\Gamma = 4$, $\epsilon_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/\Gamma = 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.</figcaption></figure>

$$ \hat{G}_{\text{HF},\sigma}^{-1} = \hat{g}_{\sigma}^{-1} - (\Delta t)^2 (\hat{\Sigma}_{T,\sigma} + \hat{\Sigma}_{\text{HF},\sigma}), \quad (15) $$

and following the numerical procedure presented in the previous section. The dynamical properties of the system can be now calculated from $\hat{G}_{\text{HF},\sigma}$.

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_\sigma(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.

$\Sigma_{\text{HF},\sigma}^{\alpha\alpha}(t_N, t_N) \approx \Sigma_{\text{HF},\sigma}^{\alpha\alpha}(t_{N-1}, t_{N-1})$. The error introduced by this approximation becomes negligible for a sufficiently small $\Delta t$. In the finite bandwidth situation this means $\Delta t \lesssim 1/D$ and in the wide band limit $\Delta 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_\sigma(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 self-consistent and the non self-consistent (first order) HF approximation for the electron-hole symmetric situation ($\epsilon_0 = -U/2$) and the $(n_\uparrow(0), n_\downarrow(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 electron-hole 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/\Gamma$. Although a good agreement is found for small values of $U/\Gamma$, the charge progressively deviates from the exact results for increasing $U/\Gamma$. This pathological behavior is corrected within a fully self-consistent HF treatment, where the average charge per spin $n_\sigma(t)$ tends to the correct singlet state for all $U/\Gamma$ 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_\uparrow(0), n_\downarrow(0)) = (1, 0)$. We have chosen a case with electron-hole symmetry ($\epsilon_0 = -U/2$) and with parameters such that $U/\pi\Gamma > 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_\uparrow(0), n_\downarrow(0)) = (0, 0)$, $(1, 1)$, the system always tends to a non-magnetic solution for all values of $U/\Gamma$.

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–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.
## 3.2. 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/\Gamma$ can include these effects in a rather satisfactory way. Indeed it can be shown that the exact self-energy in the limit $U/\Gamma \to \infty$ 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–88].

We will concentrate in the symmetric case, $\epsilon_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/\Gamma$. 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

$$ \begin{aligned} \hat{\Sigma}_{\sigma}^{(2)+-}(t, t') &= -U^2 \hat{G}_{\text{HF},\sigma}^{+-}(t, t') \hat{G}_{\text{HF},\bar{\sigma}}^{+-}(t, t') \hat{G}_{\text{HF},\bar{\sigma}}^{-+}(t', t), \\ \hat{\Sigma}_{\sigma}^{(2)-+}(t, t') &= -U^2 \hat{G}_{\text{HF},\sigma}^{-+}(t, t') \hat{G}_{\text{HF},\bar{\sigma}}^{-+}(t, t') \hat{G}_{\text{HF},\bar{\sigma}}^{+-}(t', t), \end{aligned} \quad (16) $$

where 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 $\hat{\Sigma}_{\text{int},\sigma} = \hat{\Sigma}_{\text{HF},\sigma} + \hat{\Sigma}_{\sigma}^{(2)}$.

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_{\uparrow}(0), n_{\downarrow}(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_{\sigma}(t)$ with an initial configuration ($n_{\uparrow}(0), n_{\downarrow}(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_{\sigma}(t)$ for the same initial magnetic configuration for increasing values of the electron–electron interaction parameter. It is found that for $U/\Gamma \gtrsim 8$ 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(\omega, t \to \infty)$. 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/\Gamma$ 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/\Gamma$ values.

It should be remarked that the convergence time increases with $U/\Gamma$. 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_K^{-1}$, $T_K$ being the Kondo temperature. In figures 3(b) and (c) we show the time evolution of the spectral density for two values of the interaction strength, $U/\Gamma = 4$ and $8$. The formation in time of the Kondo resonance is illustrated, showing a longer time for the larger interaction. Considering the expression for the Kondo temperature in the electron–hole symmetric Anderson model, i.e. $T_K = \sqrt{U\Gamma/2} \exp[-\pi U/8\Gamma]$, for these cases we have the ratio $T_K(U/\Gamma = 4)/T_K(U/\Gamma = 8) \simeq 3.4$. Thus, one would expect a Kondo resonance formation time for the $U/\Gamma = 8$ case roughly 3.4 times larger than for $U/\Gamma = 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/\Gamma$ value. On the other hand, figure 3(d) shows the evolution of the height of the central peak, $A(\omega = 0, t)$, to its stationary value fixed by the Friedel sum rule $A(\omega = 0, t \to \infty) = 1/\pi\Gamma$. A kink in the evolution is observed at times $\sim 1/U$, mainly visible for large $U/\Gamma$ 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
<figure><img src="image_7.png" /><figcaption>Figure 3. (a) Long time DOS for different electron–electron interactions, U/Γ = 4 (red curves) and 8 (blue curves) and V = 0. We compare the results of the perturbation expansion up to second order (full lines) with those from the exact NRG calculation of [89] (dashed lines). The inset show the convergence of the Kondo resonance. (b) and (c) time evolution of the density of states, showing the formation of the Kondo peak for U/Γ = 4 and U/Γ = 8, respectively. In (d) the height of the Kondo peak is represented as a function of time for these two cases.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. (a) Short time symmetrized current (I) = (I_L - I_R)/2, comparing the results for the perturbation expansion up to second order (solid lines) with the ones obtained using MC in [37] (symbols) for U/Γ = 0, 4 and 8, with V = 10Γ, ϵ_0 = -U/2, T = 0.2Γ and D = 10Γ. The asymptotic current as a function of the voltage is shown in (b) for increasing values of the electron–electron interaction, compared with the exact MC results from [93] (symbols).</figcaption></figure>

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.

## 4. 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

$$ \hat{H} = \hat{H}_0 + \hat{H}_{\text{ph}} + \hat{H}_{\text{e-ph}}, \qquad (17) $$

where $\hat{H}_0$ is the non-interacting part in the Hamiltonian of section 2, $\hat{H}_{\text{ph}} = \omega_0 b^\dagger b$, $\omega_0$ being the frequency of the local phonon mode and b (b$^\dagger$) the phonon annihilation (creation) operator. The electron–phonon interaction at the central region is described by the term $\hat{H}_{\text{e-ph}} = \lambda (b^\dagger + b) d^\dagger d$, where λ measures the electron–phonon coupling strength.
a)

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

b)

<figure><img src="image_11.png" /><figcaption>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 self-consistent MIGDAL [27], where the electronic propagators are fully dressed.</figcaption></figure>

## 4.1. 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

$$ \Sigma_H^{\alpha\alpha}(t, t') = \alpha\delta(t-t')\lambda^2 \int d\tau [d^{++}(t, \tau) - d^{+-}(t, \tau)]n(\tau), \quad \Sigma_H^{+-}(t, t') = \Sigma_H^{-+}(t, t') = 0, \quad (18) $$

where $n(t)$ is the self-consistent central level charge and $\hat{d}$ is the free phonon propagator in Keldysh space given by

$$ \hat{d}(t, t') = -i \begin{pmatrix} 2n_p \cos \omega_0 (t-t') + e^{-i\omega_0 |t-t'|} & n_p e^{-i\omega_0 (t-t')} + (n_p+1)e^{i\omega_0 (t-t')} \\ n_p e^{i\omega_0 (t-t')} + (n_p+1)e^{-i\omega_0 (t-t')} & 2n_p \cos \omega_0 (t-t') + e^{i\omega_0 |t-t'|} \end{pmatrix}, \quad (19) $$

where $n_p = (e^{\omega_0/T} - 1)^{-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

$$ \Sigma_H^{\alpha\beta}(t, t') = \alpha\lambda^2\delta(t-t')\delta_{\alpha\beta} \int_0^t d\tau d^R(t, \tau)n(\tau), \quad (20) $$

where $d^R(t, t')$ is the retarded free phonon propagator

$$ d^R(t, t') = -2\theta(t)\theta(t-t')\sin[\omega_0(t-t')]. \quad (21) $$

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 ($\omega_0 \gg \epsilon_0, \Gamma$) [33] with a central charge evolving adiabatically. In this limit equation (20) tends to

$$ \Sigma_H^{\alpha\beta}(t, t') \approx -\alpha\delta(t-t')\delta_{\alpha\beta} \frac{2\lambda^2}{\omega_0} n(t). \quad (22) $$

We can now follow a similar procedure to the one used in the previous section to calculate $\hat{G}_H$ 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 $\omega_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 $\tau < t$. Retardation effects of the phonon dynamics results in a coherent superposition of oscillations with period $2\pi/\omega_0$ but with different amplitudes ($\propto n(\tau)$). In the intermediate regime where the electron and the phonon dynamics are equally fast ($\omega_0 \approx \Gamma$), the coherence between those oscillations is lost at long times ($t \gg 1/\Gamma, 2\pi/\omega_0$), thus resulting in an effective damping of the central level charge, see figures 6(a) and (b). However, in the adiabatic regime ($\omega_0 \gg \Gamma$) 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 $\omega_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\lambda^2/\pi\Gamma\omega_0 > 1$. The possibility of a bistable regime for a molecular quantum dot with strong electron–phonon interaction was suggested some time ago [98–100]. The interest in investigating such a
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Time evolution of the central level charge for an initially full in (a), and empty level in (b). 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: $\omega_0 = 8\Gamma$ (black), $2\Gamma$(green) and $\Gamma$ (blue). (c) Charge evolution for an initially empty (dashed line) and initially full level (solid line) for the $\omega_0 = \Gamma$ case. The blue and red lines correspond to the Hartree and self-consistent Born approximation, respectively. The remaining parameters are $\lambda = 1.5\Gamma$, $V = 0$ and the central level is set to $\epsilon_0 = \lambda^2/\omega_0$, thus preserving electron-hole symmetry.</figcaption></figure>

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.

## 4.2. 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

$$ \hat{\Sigma}_{H}^{\alpha\beta}(t, t') = -2\lambda^2\alpha\delta_{\alpha\beta}\delta(t-t') \int_0^t d\tau \sin[\omega_0(t-\tau)]n(\tau), \qquad \hat{\Sigma}_{X}^{\alpha\beta}(t, t') = i\alpha\beta\lambda^2\hat{G}^{\alpha\beta}(t, t')\hat{d}^{\alpha\beta}(t, t'), \qquad (23) $$

where $\hat{G}^{\alpha\beta}$ denotes the Keldysh components of the fully dressed electron propagators and $n(t)$ is the final self-consistent charge.

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 $\Delta 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 $\lambda^2/\omega_0\Gamma$, 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–109]. These results suggest that the bistable behavior of the
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 (2)

$$ \hat{D} = (\hat{d}^{-1} - \hat{\Pi})^{-1}, \qquad (24) $$

where $\hat{D}(t, t') = -i\langle \hat{T}_c[\hat{\phi}(t)\hat{\phi}^{\dagger}(t')] \rangle$, with $\hat{\phi} = b + b^{\dagger}$. $\hat{d}^{-1}$ is the inverse free phonon propagator and $\hat{\Pi}$ is the phonon self-energy given by

$$ \hat{\Pi}^{\alpha\beta}(t, t') = -i\alpha\beta\lambda^2 G^{\alpha\beta}(t, t') G^{\beta\alpha}(t', t). \qquad (25) $$

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 $\hat{d}^{-1}$ we have developed a regularization procedure which is discussed in the appendix, finding

$$ \mathbf{d}^{-1} = \frac{1}{2\delta} \begin{bmatrix} \begin{pmatrix} h_0^+ & -1 & & & & & h_{0N} \\ -1 & h & -1 & & & & \\ & \ddots & \ddots & & & & \\ & & h & -1 & & & \\ & & & -1 & h_N^+ & & \\ & & & & -1 & h_0^+ & \\ & & & & & c & \\ & & & & & & c \end{pmatrix} & h_{0N} \\ \begin{pmatrix} h_0^+ & -1 & & & & & h_{0N} \\ -1 & h & -1 & & & & \\ & \ddots & \ddots & & & & \\ & & h & -1 & & & \\ & & & -1 & h_N^+ & & \\ & & & & -1 & h_0^+ & \\ & & & & & c & \\ & & & & & & c \end{pmatrix} & h_N^- \\ \begin{pmatrix} h_0^+ & -1 & & & & & h_{0N} \\ -1 & h & -1 & & & & \\ & \ddots & \ddots & & & & \\ & & h & -1 & & & \\ & & & -1 & h_N^- & & \\ & & & & -1 & h_0^- & \\ & & & & & c & \\ & & & & & & c \end{pmatrix} & h_0^- \end{bmatrix}_{2N \times 2N}, \qquad (26) $$

where $\delta = \Delta t \omega_0$ and $h = 2(1 - \delta^2/2)$. The information about the initial phonon state is encoded in the components $h_0^\pm = \pm h/2 + i\delta(1 + \rho_0^2)/(1 - \rho_0^2)$ and $h_{0N} = -2i\delta\rho_0/(1 - \rho_0^2)$, where $\rho_0 = n_p(0)/[n_p(0) + 1]$ and $n_p(0)$ is the initial phonon population. We will consider that phonons are initially in thermal equilibrium and thus $\rho_0 = e^{-\omega_0/T}$. The regularization procedure requires introducing an infinitesimal quantity $\eta$ which enters in the matrix elements connecting both Keldysh branches: $c = -2i\delta/\eta$ and $h_N^\pm = \pm h/2 - c$. The parameter $\eta$ can be interpreted as a small phonon relaxation rate which has to be taken such as $1/\eta \gg t, 1/\omega_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 $\langle \hat{T}_c[\hat{b}(t)\hat{b}^\dagger(t')] \rangle$ and $\langle \hat{T}_c[\hat{b}^\dagger(t)\hat{b}^\dagger(t')] \rangle$ in the diagrammatic expansion of $\hat{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 ($\omega_0 \gg \Gamma, \lambda$) [112, 113]. For the calculation of the phonon self-energy, $\hat{\Pi}$, 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 $\omega_0 = 2\Gamma$; a case with a rather strong electron-phonon coupling although still far from the polaronic limit ($\lambda^2/(\omega_0\Gamma) \gg 1$). Notice the dip in the DOS at $\omega \approx \omega_0$ in the self-consistent Born approximation, which is a feature due to the logarithmic divergence of the second order self-energy $\hat{\Sigma}_X(\omega)$ at $\omega = \omega_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 $\omega_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 $\sim 2\omega_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
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Long time spectral density for the self-consistent MIGDAL (solid line), RPA (dashed) and Born (dotted) approximations, compared to NRG calculations (yellow dots). The inset shows the convergence of the central peak to the expected stationary value for the Born and RPA approximations. Parameters: $\lambda = 1.5$, $\omega_0 = 2$, $V = 0$, $\epsilon_0 = \lambda^2/\omega_0$, $\Gamma = 1$ and $D = 30$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Comparison of the left (a) and (d), right (b) and (e) and symmetrized (c) and (f) current results with the MC simulations of [48]. Results for $V = 4$, (a)–(c) and 32 (d)–(f) cases are shown, and for the three approximations described in the text: self-consistent Born (dotted line), RPA (dashed line) and MIGDAL (solid line) approximations. The remaining parameters are: $\lambda = 8$, $\omega_0 = 10$, $D = 20$, $T = 0.2$ and $\Gamma = 1$.</figcaption></figure>

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–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 ($\lambda = 8$, $\omega_0 = 10$ and $\Gamma = 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.

## 5. 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><img src="image_15.png" /><figcaption>Figure 9. Spectral density in the spin-degenerate Anderson-Holstein model. (a) Long time values for RPA (lines) compared with equilibrium NRG results (symbols) from [116], for two different electron-phonon coupling parameters: $\lambda = 1.89$ (red) and $\lambda = 3.14$ (blue). The remaining parameters are $U = 6.3$, $\omega_0 = 3.14$, $\epsilon_0 = \lambda^2/\omega_0 - U/2$ and $\Gamma_R = \Gamma_L = 0.5$. (b) and (c) time evolution of the density of states for $\lambda = 0$ and $\lambda = 2$. In (d) we represent the central peak height evolution, showing in red the $\lambda = 0$ and in blue the $\lambda = 2$ case. The remaining parameters are $U = 8$, $\omega_0 = 2$, $\epsilon_0 = \lambda^2/\omega_0 - U/2$, $\Gamma = 1$ and $V = 0$.</figcaption></figure>

$$ \hat{H} = \sum_{\sigma=\uparrow,\downarrow} \hat{H}_{0,\sigma} + \hat{H}_{e-e} + \hat{H}_{ph} + \hat{H}_{e-ph}, \quad (27) $$

where $\hat{H}_{0,\sigma}$ is the non-interacting part of the Hamiltonian given in section 2, $\hat{H}_{e-e} = U\hat{n}_{\uparrow}\hat{n}_{\downarrow}$, $\hat{H}_{ph} = \omega_0 b^{\dagger}b$ and $\hat{H}_{e-ph} = \lambda(b + b^{\dagger})\sum_{\sigma} \hat{n}_{\sigma}$. In this case we combine the approximations used in section 4 for the electron- electron self-energies with the ones in the previous section for the electron-phonon case, i.e. $\Sigma_{int} = \Sigma_{e-e} + \Sigma_{e-ph}$ (see figures 1 and 5).

In figure 9(a) we show the long time spectral density compared to the exact NRG results from [116] using the RPA for $\Sigma_{e-ph}$. Similar results are obtain for the MIGDAL approximation. As can be observed, for the smaller $\lambda$ case the RPA exhibits an overall agreement with the exact results. However, for larger values of the electron-phonon 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 $\lambda^2/\omega_0 \gtrsim U/2$ [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 $\lambda/\Gamma = 0$ in figure 9(b) and $\lambda/\Gamma = 2$ in figure 9(c), with $U/\Gamma = 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(\omega = 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 $\sim 2\pi/\omega_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.

## 6. 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
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 self-energy in the spirit of [107–109].

## Acknowledgments

RSS, ALY and AMR acknowledge financial support by Spanish MINECO through Grants No. FIS2014-55486-P and FIS2017-84860-R, and the María de Maeztu Programme for Units of Excellence in R&D (Grant No. MDM-2014-0377). RA acknowledges support of Conseil Regional de la Nouvelle Aquitaine.

## Appendix. Inverse free boson propagator

In this appendix we discuss the problem of obtaining the inverse of the free phonon propagator discretized along the Keldysh contour. This problem has already been discussed by Kamenev in [79], where the author considers the problem of bosonic particles occupying a single level of energy $\omega_0$

$$ \hat{H}_{\text{ph}} = \omega_0 b^{\dagger} b, \qquad (\text{A1}) $$

with the free phonon propagator defined as $\hat{d}_0(t, t') = -i\langle T_c b(t) b^{\dagger}(t') \rangle$. The inverse propagator in this case is formally similar to the electronic one (3), finding

$$ \text{id}_0^{-1} = \begin{pmatrix} -1 & & & & \rho(\omega_0) \\ h_- & -1 & & & \\ h_- & -1 & \ddots & & \\ & & \ddots & \ddots & \\ & & & h_+ & -1 \end{pmatrix}_{2N \times 2N}, \qquad (\text{A2}) $$

with $h_{\pm} = 1 \pm i\Delta t \omega_0$. This expression constitutes a discretized version of the $i\partial_t - \omega_0$ operator on the time contour with an initial condition $\rho(\omega_0) = n_p(0)/[1 + n_p(0)]$, which depends on the initial phonon population $n_p(0)$. The obtention of the inverse free phonon propagator defined as $\hat{d}(t, t') = -i\langle \hat{T}_c[\hat{\varphi}(t) \hat{\varphi}^{\dagger}(t')] \rangle$, with $\hat{\varphi} = b + b^{\dagger}$ becomes more demanding since it involves the discretization of the second order differential operator $\hat{H} = p^2/2 + \omega_0^2 x^2/2$ with $p = -i\partial_x$ and $x = \sqrt{1/2\omega_0} \hat{\varphi}$. Moreover, it can be checked that the discretized version of the free phonon propagator given in equation (19) is not invertible as it becomes singular. In this section we discuss the way to obtain this inverse propagator by including a regularization procedure. By definition, the system partition function is given by [79]

$$ Z = \frac{\text{Tr}[\mathcal{U}_c \rho]}{\text{Tr}[\rho]}, \qquad (\text{A3}) $$

where $\mathcal{U}_c = \mathcal{U}^{\dagger}(t_{2N}, t_{N+1})\mathcal{U}^{-}(t_N, t_1)$ is the contour evolution operator and $\rho = e^{-H/T}$ is the initial density matrix. Expanding $Z$ in coordinate space and for $N = 3$ we find

$$ \text{Tr}[\mathcal{U}_c \rho] = \int dx_1 \dots dx_6 \langle x_6 | \mathcal{U}_{-\Delta t} | x_5 \rangle \langle x_5 | \mathcal{U}_{-\Delta t} | x_4 \rangle \langle x_4 | \mathbf{1} | x_3 \rangle \langle x_3 | \mathcal{U}_{\Delta t} | x_2 \rangle \langle x_2 | \mathcal{U}_{\Delta t} | x_1 \rangle \langle x_1 | \rho | x_6 \rangle, \qquad (\text{A4}) $$

where $\mathcal{U}_{\Delta t} = e^{-iH\Delta t}$. It is worth noticing that the last term in the integrand correspond to the contour closing and the third one is the branch changing in the Keldysh contour at the final time. The relevant matrix components are given by so-called Mehler kernel [118]

$$ \langle x | e^{-iHt} | y \rangle = \frac{\exp\{i[(x^2 + y^2)\cos(\omega_0 t) - 2xy]/2 \sin(\omega_0 t)\}}{\sqrt{2\pi i \sin(\omega_0 t)}}. \qquad (\text{A5}) $$
Discretizing the expression and considering the time step $\Delta t$ as the smallest timescale we find

$$ \langle x|e^{\mp iH\Delta t}|y\rangle = \frac{\exp\{\pm i[(x^2 + y^2)(1 - \delta^2/2) - 2xy]/2\delta\}}{\sqrt{2\pi i\delta}} \quad (A6) $$

with $\delta = \omega_0\Delta t$. A similar expression can be found for the contour closing term

$$ \langle x|\rho|y\rangle = \sqrt{\frac{\rho_0}{\pi(1 - \rho_0^2)}} \exp\left[-\frac{(1 + \rho_0^2)(x^2 + y^2)}{2(1 - \rho_0^2)} + \frac{2xy\rho_0}{1 - \rho_0^2}\right], \quad (A7) $$

where $\rho_0 = n_p(0)/[n_p(0) + 1]$ contains information about the initial phonon population, $n_p(0)$. The final step for obtaining the inverse is to regularize the delta function, i.e. we should take

$$ \langle x|\mathbf{1}|y\rangle \approx \frac{1}{\sqrt{2\pi\eta}} e^{-(x-y)^2/2\eta}, \quad (A8) $$

being $\eta$ an infinitesimum. Finally, the inverse of the free phonon propagator can be obtained identifying the components of

$$ Z = \int dx_1 \dots dx_{2N} e^{ix^T d^{-1} x}, \quad (A9) $$

finding the expression of equation (26). It is worth commenting that all the prefactors in the Mehler kernel expression normalize the partition function, without affecting the phonon dynamics. The particular case for $N = 2$ can be written as

$$ \mathrm{id}_{N=2} = \begin{pmatrix} 1 - \frac{\delta^2}{2} + i\delta \frac{1 + \rho_0^2}{1 - \rho_0^2} & -1 & 0 & -2i \frac{\delta \rho_0}{1 - \rho_0^2} \\ -1 & 1 - \frac{\delta^2}{2} + 2i \frac{\delta}{\eta} & -2i \frac{\delta}{\eta} & 0 \\ 0 & -2i \frac{\delta}{\eta} & -1 + \frac{\delta^2}{2} + 2i \frac{\delta}{\eta} & 1 \\ -2i \frac{\delta \rho_0}{1 - \rho_0^2} & 0 & 1 & -1 + \frac{\delta^2}{2} + i \delta \frac{1 + \rho_0^2}{1 - \rho_0^2} \end{pmatrix}. \quad (A10) $$

## ORCID iDs

R Seoane Souto https://orcid.org/0000-0002-2978-3534

## References

[1] Nazarov Y and Blanter Y 2009 *Quantum Transport: Introduction to Nanoscience* (Cambridge: Cambridge University Press)

[2] Blandin A, Nourtier A and Hone D W 1976 *J. Phys.* **37** 369

[3] Jauho A-P, Wingreen N S and Meir Y 1994 *Phys. Rev.* B **50** 5528

[4] Fève G, Mahé A, Berroir J-M, Kontos T, Plaçais B, Glattli D C, Cavanna A, Etienne B and Jin Y 2007 *Science* **316** 1169

[5] Vink I T, Nooitgedagt T, Schouten R N, Vandersypen L M K and Wegscheider W 2007 *Appl. Phys. Lett.* **91** 123512

[6] Flindt C, Fricke C, Hohls F, Novotný T, Netočný K, Brandes T and Haug R J 2009 *Proc. Natl Acad. Sci.* **106** 10116

[7] Terada Y, Yoshida S, Takeuchi O and Shigekawa H 2010 *Nat. Photon.* **4** 869

[8] Loth S, Etzkorn M, Lutz C P, Eigler D M and Heinrich A J 2010 *Science* **329** 1628

[9] Latta C et al 2011 *Nature* **474** 627

[10] Yoshida S, Aizawa Y, Wang Z-H, Oshima R, Mera Y, Matsuyama E, Oigawa H, Takeuchi O and Shigekawa H 2014 *Nat. Nanotechnol.* **9** 588

[11] Otsuka T et al 2017 *Sci. Rep.* **7** 12201

[12] Karnetzky C, Zimmermann P, Trummer C, Duque-Sierra C, Wörle M, Kienberger R and Holleitner A 2018 *Nat. Commun.* **9** 2471

[13] Du S, Yoshida K, Zhang Y, Hamada I and Hirakawa K 2018 arXiv:1712.07339

[14] Hung C-L, Gurarie V and Chin C 2013 *Science* **341** 1213

[15] Schreiber M, Hodgman S S, Bordia P, Lüschen H P, Fischer M H, Vosk R, Altman E, Schneider U and Bloch I 2015 *Science* **349** 842

[16] Cazalilla M A 2006 *Phys. Rev. Lett.* **97** 156403

[17] Heyl M, Polkovnikov A and Kehrein S 2013 *Phys. Rev. Lett.* **110** 135704

[18] Eisert J, Friesdorf M and Gogolin C 2015 *Nat. Phys.* **11** 124

[19] Cini M 1980 *Phys. Rev.* B **22** 5887

[20] Stefanucci G and Almbladh C-O 2004 *Phys. Rev.* B **69** 195318

[21] Albert M, Flindt C and Büttiker M 2011 *Phys. Rev. Lett.* **107** 086805

[22] Albert M, Haack G, Flindt C and Büttiker M 2012 *Phys. Rev. Lett.* **108** 186806

[23] Dasenbrook D, Flindt C and Büttiker M 2014 *Phys. Rev. Lett.* **112** 146801

[24] Tang G-M, Xu F and Wang J 2014 *Phys. Rev.* B **89** 205310

[25] Tang G-M and Wang J 2014 *Phys. Rev.* B **90** 195422

[26] Ridley M, MacKinnon A and Kantorovich L 2015 *Phys. Rev.* B **91** 125433

[27] Murakami Y, Werner P, Tsuji N and Aoki H 2015 *Phys. Rev.* B **91** 045128

[28] Seoane Souto R, Martín-Rodero A and Yeyati A L 2017 *Fortschr. Phys.* **65** 1600062
[29] Odashima M M and Lewenkopf C H 2017 Phys. Rev. B 95 104301

[30] Covito F, Eich F G, Tuovinen R, Sentef M A and Rubio A 2018 J. Chem. Theory Comput. 14 2495–504

[31] Kambly D, Flindt C and Büttiker M 2011 Phys. Rev. B 83 075432

[32] Stegmann P and König J 2017 Phys. Status Solidi b 254 1600507

[33] Riwar R-P and Schmidt T L 2009 Phys. Rev. B 80 125109

[34] Nordlander P, Pustilnik M, Meir Y, Wingreen N S and Langreth D C 1999 Phys. Rev. Lett. 83 808

[35] Plihal M, Langreth D C and Nordlander P 2005 Phys. Rev. B 71 165321

[36] Goker A, Friedman B A and Nordlander P 2007 J. Phys.: Condens. Matter 19 376206

[37] Schmidt T L, Werner P, Mühlbacher L and Komnik A 2008 Phys. Rev. B 78 235110

[38] Komnik A 2009 Phys. Rev. B 79 245102

[39] Myöhänen P, Stan A, Stefanucci G and van Leeuwen R 2009 Phys. Rev. B 80 115107

[40] Myöhänen P, Tuovinen R, Korhonen T, Stefanucci G and van Leeuwen R 2012 Phys. Rev. B 85 075105

[41] Perfetto E and Stefanucci G 2013 Phys. Rev. B 88 245437

[42] Latini S, Perfetto E, Uimonen A-M, van Leeuwen R and Stefanucci G 2014 Phys. Rev. B 89 075306

[43] Vinkler-Aviv Y, Schiller A and Anders F B 2014 Phys. Rev. B 90 155110

[44] Seoane Souto R, Avriller R, Monreal R C, Martin-Rodero A and Yeyati A L 2015 Phys. Rev. B 92 125435

[45] Chen H-T, Cohen G, Millis A J and Reichman D R 2016 Phys. Rev. B 93 174309

[46] Tang G, Yu Z and Wang J 2017 New J. Phys. 19 083007

[47] Tang G, Xing Y and Wang J 2017 Phys. Rev. B 96 075417

[48] Mühlbacher L and Rabani E 2008 Phys. Rev. Lett. 100 176403

[49] Albrecht K F, Wang H, Mühlbacher L, Thoss M and Komnik A 2012 Phys. Rev. B 86 081412

[50] Cohen G, Gull E, Reichman D R, Millis A J and Rabani E 2013 Phys. Rev. B 87 195108

[51] Cohen G, Gull E, Reichman D R and Millis A J 2014 Phys. Rev. Lett. 112 146802

[52] Härtle R, Cohen G, Reichman D R and Millis A J 2015 Phys. Rev. B 92 085430

[53] Klatt J, Mühlbacher L and Komnik A 2015 Phys. Rev. B 91 155306

[54] Albrecht K F, Martin-Rodero A, Schachenmayer J and Mühlbacher L 2015 Phys. Rev. B 91 064305

[55] Ridley M, Singh V N, Gull E and Cohen G 2018 Phys. Rev. B 97 115109

[56] Anders F B and Schiller A 2005 Phys. Rev. Lett. 95 196801

[57] Anders F B and Schiller A 2006 Phys. Rev. B 74 245113

[58] Heidrich-Meisner F, Feiguin A E and Dagotto E 2009 Phys. Rev. B 79 235336

[59] Eckel J, Heidrich-Meisner F, Jakobs S G, Thorwart M, Pletyukhov M and Egger R 2010 New J. Phys. 12 043042

[60] Eidelstein E, Schiller A, Güttge F and Anders F B 2012 Phys. Rev. B 85 075118

[61] Güttge F, Anders F B, Schollwöck U, Eidelstein E and Schiller A 2013 Phys. Rev. B 87 115115

[62] Nghiem H T M and Costi T A 2014 Phys. Rev. B 89 075118

[63] Nghiem H T M and Costi T A 2017 Phys. Rev. Lett. 119 156601

[64] Zheng X, Wang F, Yam C Y, Mo Y and Chen G 2007 Phys. Rev. B 75 195127

[65] Kurth S, Stefanucci G, Khosravi E, Verdozzi C and Gross E K U 2010 Phys. Rev. Lett. 104 236801

[66] Uimonen A-M, Khosravi E, Stan A, Stefanucci G, Kurth S, van Leeuwen R and Gross E K U 2011 Phys. Rev. B 84 115103

[67] Khosravi E, Uimonen A-M, Stan A, Stefanucci G, Kurth S, van Leeuwen R and Gross E K U 2012 Phys. Rev. B 85 075103

[68] Kwok Y, Zhang Y and Chen G 2014 Frontiers Phys. 9 698

[69] Dittmann N, Splettstoesser J and Helbig N 2018 Phys. Rev. Lett. 120 157701

[70] Kurth S and Stefanucci G 2018 Eur. Phys. J. B 91 118

[71] Wang H and Thoss M 2003 J. Chem. Phys. 119 1289

[72] Weiss S, Eckel J, Thorwart M and Egger R 2008 Phys. Rev. B 77 195316

[73] Kennes D M, Jakobs S G, Karrasch C and Meden V 2012 Phys. Rev. B 85 085113

[74] Kennes D M and Meden V 2012 Phys. Rev. B 85 245101

[75] Perfetto E and Stefanucci G 2015 J. Comput. Electron. 14 352

[76] Keldysh L V 1965 Sov. Phys.—JETP 20 1018

[77] Souto R S, Martín-Rodero A and Yeyati A L 2016 Phys. Rev. Lett. 117 267701

[78] Souto R S, Martín-Rodero A and Yeyati A L 2017 Phys. Rev. B 96 165444

[79] Kamenev A 2011 Field Theory of Non-Equilibrium Systems (Cambridge: Cambridge University Press)

[80] Chen H-T, Cohen G, Millis A J and Reichman D R 2016 Phys. Rev. B 93 174309

[81] Anderson P W 1961 Phys. Rev. 124 41

[82] Wiegmann P B 1980 Phys. Lett. 80A 163

[83] Kawakami N and Okiji A 1981 Phys. Lett. A 86 483

[84] Andrei N, Furuya K and Lowenstein J H 1983 Rev. Mod. Phys. 55 331

[85] Martin-Rodero A, Flores F, Baldo M and Pucci R 1982 Solid State Commun. 44 911

[86] Martín-Rodero A, Louis E, Flores F and Tejedor C 1986 Phys. Rev. B 33 1814

[87] Yeyati A L, Martín-Rodero A and Flores F 1993 Phys. Rev. Lett. 71 2991

[88] Kajueter H and Kotliar G 1996 Phys. Rev. Lett. 77 131

[89] Anders F B 2008 J. Phys.: Condens. Matter 20 195216

[90] White J A 1992 Phys. Rev. B 45 1100

[91] Han J E and Heary R J 2007 Phys. Rev. Lett. 99 236808

[92] Nordlander P, Pustilnik M, Meir Y, Wingreen N S and Langreth D C 1999 Phys. Rev. Lett. 83 808

[93] Werner P, Oka T and Millis A J 2009 Phys. Rev. B 79 035320

[94] Fujii T and Ueda K 2003 Phys. Rev. B 68 155310

[95] Holstein T 1959 Ann. Phys., NY 8 325

[96] Utsumi Y, Entin-Wohlman O, Ueda A and Aharony A 2013 Phys. Rev. B 87 115407

[97] D'Amico P, Ryndyk D A, Cuniberti G and Richter K 2008 New J. Phys. 10 085002

[98] Gogolin A and Komnik A 2002 arXiv:cond-mat/0207513

[99] Alexandrov A S and Bratkovsky A M 2003 Phys. Rev. B 69 235312

[100] Galperin M, Ratner M A and Nitzan A 2005 Nano Lett. 5 125

[101] Micchi G, Avriller R and Pistolesi F 2015 Phys. Rev. Lett. 115 206802
[102] Micchi G, Avriller R and Pistolesi F 2016 *Phys. Rev. B* **94** 125417

[103] Avriller R, Murr B and Pistolesi F 2018 *Phys. Rev. B* **97** 155414

[104] Wilner E Y, Wang H, Cohen G, Thoss M and Rabani E 2013 *Phys. Rev. B* **88** 045137

[105] Albrecht K F, Martin-Rodero A, Monreal R C, Mühlbacher L and Yeyati A L 2013 *Phys. Rev. B* **87** 085127

[106] Martin-Rodero A, Yeyati A L, Flores F and Monreal R C 2008 *Phys. Rev. B* **78** 235112

[107] Maier S, Schmidt T L and Komnik A 2011 *Phys. Rev. B* **83** 085401

[108] Dong B, Ding G H and Lei X L 2013 *Phys. Rev. B* **88** 075414

[109] Seoane Souto R, Yeyati A L, Martin-Rodero A and Monreal R C 2014 *Phys. Rev. B* **89** 085412

[110] Wilner E Y, Wang H, Thoss M and Rabani E 2014 *Phys. Rev. B* **90** 115145

[111] Wilner E Y, Wang H, Thoss M and Rabani E 2015 *Phys. Rev. B* **92** 195143

[112] Meyer D, Hewson A C and Bulla R 2002 *Phys. Rev. Lett.* **89** 196401

[113] Agarwalla B K, Kulkarni M, Mukamel S and Segal D 2016 *Phys. Rev. B* **94** 035434

[114] Laakso M A, Kennes D M, Jakobs S G and Meden V 2014 *New J. Phys.* **16** 023007

[115] Jovchev A and Anders F B 2013 *Phys. Rev. B* **87** 195112

[116] Jeon G S, Park T-H and Choi H-Y 2003 *Phys. Rev. B* **68** 045106

[117] Hewson A C and Meyer D 2002 *J. Phys.: Condens. Matter* **14** 427

[118] Erdélyi A 1953 *Higher Transcendental Functions* vol 2 (New York: McGraw-Hill) Bateman Project