Time evolution of the Kondo resonance in response to a quench

TL;DR

This paper studies how the Kondo resonance evolves over time after a sudden change using the TDNRG method, revealing the different time scales for the development of the resonance and satellite peaks.

Contribution

It introduces a new TDNRG-based approach to compute two-time nonequilibrium Green functions and spectral functions in the Anderson impurity model after a quench.

Findings

Kondo resonance fully develops at long times $t \,\gtrsim\, 1/T_K$

Satellite peaks develop rapidly within $|t| \,\lesssim\, 1/\Gamma$

Initial and final spectral functions are recovered at $t\to -\infty$ and $t\to +\infty$

Abstract

We investigate the time evolution of the Kondo resonance in response to a quench by applying the time-dependent numerical renormalization group (TDNRG) approach to the Anderson impurity model in the strong correlation limit. For this purpose, we derive within TDNRG a numerically tractable expression for the retarded two-time nonequilibrium Green function $G(t+t',t)$, and its associated time-dependent spectral function, $A(\omega,t)$, for times $t$ both before and after the quench. Quenches from both mixed valence and Kondo correlated initial states to Kondo correlated final states are considered. For both cases, we find that the Kondo resonance in the zero temperature spectral function, a preformed version of which is evident at very short times $t\to 0^{+}$, only fully develops at very long times $t\gtrsim 1/T_{\rm K}$, where $T_{\rm K}$ is the Kondo temperature of the final state. In contrast, the final state satellite peaks develop on a fast time scale $1/\Gamma$ during the time interval $-1/\Gamma \lesssim t \lesssim +1/\Gamma$, where $\Gamma$ is the hybridization strength. Initial and final state spectral functions are recovered in the limits $t\rightarrow -\infty$ and $t\rightarrow +\infty$, respectively. Our formulation of two-time nonequilibrium Green functions within TDNRG provides a first step towards using this method as an impurity solver within nonequilibrium dynamical mean field theory.

Figures

Equations

$$G(\omega +i\eta ,t)=\sum _{m=m_0}^N\sum _{rsq}^{\notin K{K}^{\prime }{K}^{\prime \prime }}{\rho }_{sr}^{i\to f}(m)e^{-i(E_s^m-E_r^m)t}$$

$$\displaystyle\times\Big{(}\frac{B^{m}_{rq}C^{m}_{qs}}{\omega+E^{m}_{r}-E^{m}_{q}+i\eta}+\frac{C^{m}_{rq}B^{m}_{qs}}{\omega+E^{m}_{q}-E^{m}_{s}+i\eta}\Big{)},$$

$$G_{BC}(t+t^{\prime },t)$$

$$=-i\theta (t^{\prime })\mathrm{Tr}\{\hat{\rho }[e^{i{H}_f(t+t^{\prime })}\hat{B}{e}^{-i{H}_f(t+t^{\prime })},e^{i{H}_{f}t}\hat{C}{e}^{-i{H}_{f}t}{]}_s\}$$

$$=-i\theta (t^{\prime })\mathrm{Tr}\{e^{-i{H}_{f}t}\hat{\rho }{e}^{i{H}_{f}t}[e^{i{H}_f{t}^{\prime }}\hat{B}{e}^{-i{H}_f{t}^{\prime }},C{]}_s\}$$

$$=-i\theta (t^{\prime })\mathrm{Tr}\{\hat{\rho }(t)[\hat{B}(t^{\prime }),\hat{C}{]}_s\},$$

$$I_1(t+t^{\prime },t)=-$$

$$=-$$

$${}_f⟨r{e}^{\prime }m|\hat{B}(t^{\prime })\hat{C}|sem{⟩}_f$$

$$=\sum _{q{e}^{\prime \prime }}{}_f⟨r{e}^{\prime }m|\hat{B}(t^{\prime })|q{e}^{\prime \prime }m{⟩}_f{}_f⟨q{e}^{\prime \prime }m|\hat{C}|sem{⟩}_f$$

$$+\sum _{m^{\prime }=m_0}^{m-1}\sum _{l{e}^{\prime \prime }}{}_f⟨r{e}^{\prime }m|\hat{B}(t^{\prime })|l{e}^{\prime \prime }{m}^{\prime }{⟩}_f{}_f⟨l{e}^{\prime \prime }{m}^{\prime }|\hat{C}|sem{⟩}_f.$$

$$\sum _{m^{\prime }=m_0}^{m-1}\sum _{l{e}^{\prime \prime }}{}_f⟨r{e}^{\prime }m|(1_{m^{\prime }}^{+}+1_{m^{\prime }}^{-})\hat{B}(t^{\prime })|l{e}^{\prime \prime }{m}^{\prime }{⟩}_f{}_f⟨l{e}^{\prime \prime }{m}^{\prime }|\hat{C}(1_{m^{\prime }}^{+}+1_{m^{\prime }}^{-})|sem{⟩}_f$$

$$=$$

$${}_f⟨r{e}^{\prime }m|\hat{B}(t^{\prime })\hat{C}|sem{⟩}_f=$$

$$+$$

$$I_1(t+t^{\prime },t)=$$

$$\displaystyle+\sum_{m^{\prime}=m_{0}}^{m-1}\sum_{lk_{1}k_{2}}\sum_{\alpha_{m}...\alpha_{m^{\prime}+1}}[A^{\alpha_{m}\dagger}_{XK}...A^{\alpha_{m^{\prime}+1}\dagger}_{KK}]_{rk_{1}}B^{m^{\prime}}_{k_{1}l}e^{i(E^{m^{\prime}}_{k_{1}}-E^{m^{\prime}}_{l})t^{\prime}}C^{m^{\prime}}_{lk_{2}}[A^{\alpha_{m^{\prime}+1}}_{KK}...A^{\alpha_{m}}_{KX^{\prime}}]_{k_{2}s}\Big{\}},$$

$$I_1(\omega +i\eta ,t)=$$

$$\displaystyle+\sum_{m^{\prime}=m_{0}}^{m-1}\sum_{lk_{1}k_{2}}\sum_{\alpha_{m}...\alpha_{m^{\prime}+1}}[A^{\alpha_{m}\dagger}_{XK}...A^{\alpha_{m^{\prime}+1}\dagger}_{KK}]_{rk_{1}}\frac{B^{m^{\prime}}_{k_{1}l}C^{m^{\prime}}_{lk_{2}}}{\omega+E^{m^{\prime}}_{k_{1}}-E^{m^{\prime}}_{l}+i\eta}[A^{\alpha_{m^{\prime}+1}}_{KK}...A^{\alpha_{m}}_{KX^{\prime}}]_{k_{2}s}\Big{\}}.$$

$$I_2(\omega +i\eta ,t)=$$

$$\displaystyle+\sum_{m^{\prime}=m_{0}}^{m-1}\sum_{lk_{1}k_{2}}\sum_{\alpha_{m}...\alpha_{m^{\prime}+1}}[A^{\alpha_{m}\dagger}_{XK}...A^{\alpha_{m^{\prime}+1}\dagger}_{KK}]_{rk_{1}}\frac{C^{m^{\prime}}_{k_{1}l}B^{m^{\prime}}_{lk_{2}}}{\omega+E^{m^{\prime}}_{l}-E^{m^{\prime}}_{k_{2}}+i\eta}[A^{\alpha_{m^{\prime}+1}}_{KK}...A^{\alpha_{m}}_{KX^{\prime}}]_{k_{2}s}\Big{\}}.$$

$${\tilde{\rho }}_{k_2{k}_1}^{i\to f}(m^{\prime },t)=\sum _{m=m^{\prime }+1}^N\sum _{rs}^{\notin K{K}^{\prime }}\sum _{{\alpha }_m...{\alpha }_{m^{\prime }+1}}[A_{KK}^{{\alpha }_{m^{\prime }+1}}...A_{K{X}^{\prime }}^{{\alpha }_m}{]}_{k_{2}s}{e}^{-i(E_s^m-E_r^m)t}{\rho }_{s,r}^{i\to f}(m)[A_{XK}^{{\alpha }_m†}...A_{KK}^{{\alpha }_{m^{\prime }+1}†}{]}_{r{k}_1}$$

$$\displaystyle=\begin{dcases}0&\text{if }m^{\prime}=N;\\ \sum_{\alpha_{m^{\prime}+1}}\Big{\{}\sum^{\notin KK^{\prime}}_{rs}A^{\alpha_{m^{\prime}+1}}_{k_{2}s}e^{-i(E^{m^{\prime}+1}_{s}-E^{m^{\prime}+1}_{r})t}\rho^{i\to f}_{sr}(m^{\prime}+1)A^{\alpha_{m^{\prime}+1}\dagger}_{rk_{1}}+\sum_{kk^{\prime}}A^{\alpha_{m^{\prime}+1}}_{k_{2}k}\tilde{\rho}^{i\to f}_{kk^{\prime}}(m^{\prime}+1,t)A^{\alpha_{m^{\prime}+1}\dagger}_{k^{\prime}k_{1}}\Big{\}}&\text{otherwise}\end{dcases}.$$

$$G(\omega ,t)=$$

$$\displaystyle+\sum_{m^{\prime}=m_{0}}^{N-1}\sum_{k_{1}k_{2}}\tilde{\rho}^{i\to f}_{k_{2}k_{1}}(m^{\prime},t)\sum_{l}\Big{\{}\frac{B^{m^{\prime}}_{k_{1}l}C^{m^{\prime}}_{lk_{2}}}{\omega+E^{m^{\prime}}_{k_{1}}-E^{m^{\prime}}_{l}+i\eta}+\frac{C^{m^{\prime}}_{k_{1}l}B^{m^{\prime}}_{lk_{2}}}{\omega+E^{m^{\prime}}_{l}-E^{m^{\prime}}_{k_{2}}+i\eta}\Big{\}},$$

$$I_1(t+t^{\prime },t)=-$$

$$+$$

$$-i\sum _{m=m_0}^N\sum _{rs}^{\notin K{K}^{\prime }}\sum _e{}_f⟨sem|\hat{\rho }(t)|rem{⟩}_f\sum _q{e}^{i(E_r^m-E_q^m)t^{\prime }}{B}_{rq}^m{C}_{qs}^m,$$

$$-i\sum _{m=m_0}^N\sum _{rs}^{\notin K{K}^{\prime }}\sum _{e{e}^{\prime }}{}_f⟨sem|\hat{\rho }(t)|r{e}^{\prime }m{⟩}_f\sum _{m^{\prime }=m_0}^{m-1}\sum _{l{e}^{\prime \prime }}\sum _{k_1{k}_2}{}_f⟨r{e}^{\prime }m|k_1{e}^{\prime \prime }{m}^{\prime }{⟩}_f{B}_{k_{1}l}^{m^{\prime }}{e}^{i(E_{k_1}^{m^{\prime }}-E_l^{m^{\prime }})t^{\prime }}{C}_{l{k}_2}^{m^{\prime }}{}_f⟨k_2{e}^{\prime \prime }{m}^{\prime }|sem{⟩}_f$$

$$=$$

$$I_1(t+t^{\prime },t)=$$

$$-i\sum _{m=m_0}^{N-1}\sum _{k{k}^{\prime }}\sum _e{}_f⟨kem|\hat{\rho }(t)|k^{\prime }em{⟩}_f\sum _l{B}_{k^{\prime }l}^m{e}^{i(E_{k^{\prime }}^m-E_l^m)t^{\prime }}{C}_{lk}^m$$

$$=$$

$$I_2(t+t^{\prime },t)=$$

$$G(t+t^{\prime },t)=$$

$$G(\omega +i\eta ,t)=$$

$$\hat{\rho }$$

$$\sum _e{}_f⟨sem|\hat{\rho }(t)|rem{⟩}_f$$

Full text

Introduction.— The nonequilibrium properties of strongly correlated quantum impurity models continue to pose a major theoretical challenge. This contrasts with their equilibrium properties, which are largely well understood Hewson (1993), or can be investigated within a number of highly accurate methods, such as the numerical renormalization group method (NRG) Wilson (1975); Krishna-murthy et al. (1980); Bulla et al. (2008); Gonzalez-Buxton and Ingersent (1998), the continuous time quantum Monte Carlo (CTQMC) approach Gull et al. (2011), the density matrix renormalization group White (1992), or the Bethe ansatz method Tsvelick and Wiegmann (1983); Andrei (2013). Quantum impurity models far from equilibrium are of direct relevance to several fields of research, including charge transfer effects in low-energy ion-surface scattering Brako and Newns (1981); Kasai and Okiji (1987); Merino and Marston (1998); Langreth and Nordlander (1991); Shao et al. (1994a, b); Pamperin et al. (2015); He and Yarmoff (2010), transient and steady state effects in molecular and semiconductor quantum dots Hershfield et al. (1992); Hershfield (1993); Meir et al. (1993); Bruder and Schoeller (1994); Kretinin et al. (2011, 2012); Pletyukhov and Schoeller (2012); Scott et al. (2013); Nordlander et al. (1999); Park et al. (2002); Kogan et al. (2004); Hemingway et al. (2014); Jauho et al. (1994); Kennes et al. (2012); Cohen et al. (2014); Schmidt et al. (2008); Dorda et al. (2014); Rosch et al. (2003); Antipov et al. (2016), and also in the context of dynamical mean field theory (DMFT) of strongly correlated lattice models Metzner and Vollhardt (1989); Georges et al. (1996); Kotliar and Vollhardt (2004), as generalized to nonequilibrium Schmidt and Monien (2002); Freericks et al. (2006); Aoki et al. (2014). In the latter, further progress hinges on an accurate non-perturbative solution for the nonequilibrium Green functions of an effective quantum impurity model. Such a solution, beyond allowing time-resolved spectroscopies of correlated lattice systems within DMFT to be addressed Eckstein and Kollar (2008); Freericks et al. (2009); Perfetti et al. (2006); Loukakos et al. (2007); Iyoda and Ishihara (2014), would also be useful in understanding time-resolved scanning tunnelling microscopy of nanoscale systems Loth et al. (2010) and proposed cold atom realizations of Kondo correlated states Nishida (2013); Bauer et al. (2013); Nishida (2016); Riegger et al. (2017), which could be probed with real-time radio-frequency spectroscopy Goold et al. (2011); Knap et al. (2012); Cetina et al. (2016).

In this Letter, we use the time-dependent numerical renormalization group (TDNRG) approach Anders and Schiller (2005, 2006); Anders (2008a, b); Güttge et al. (2013); Nghiem and Costi (2014a, b) to calculate the retarded two-time Green function, $G(t_{1}=t+t^{\prime},t_{2}=t)$, and associated spectral function, $A(\omega,t)$, of the Anderson impurity model in response to a quench at time $t=0$, and apply this to investigate in detail the time evolution of the Kondo resonance. This topic has been addressed before within several approaches, including the non-crossing approximation Nordlander et al. (1999); Randi et al. (2017), conserving approximations Bock et al. (2016) and within CTQMC for quantum dots out of equilibrium Cohen et al. (2014). Related work on the temporal evolution of the spin-spin correlation function in the Kondo model and thermalization in the Anderson impurity model following initial state preparations has also been carried out Lobaskin and Kehrein (2005); Heyl and Kehrein (2010). Formulations of the time-dependent spectral function within TDNRG are also available Anders (2008b); Weymann et al. (2015), but only for positive times. Here, we derive expressions for the two-time Green function and spectral function $A(\omega,t)$ which are numerically tractable at arbitrary times, including negative times. The main advantages of the TDNRG over other approaches for calculating time-dependent spectral functions is that it can access arbitrary long times ($t\to\pm\infty$) and arbitrary low temperatures and frequencies, is non-perturbative and calculates spectral functions directly on the real frequency axis. It is therefore well suited for investigating the formation in time of the exponentially narrow and low temperature Kondo resonance 111The use of a discretized Wilson chain within TDNRG results, to a small degree, in incomplete thermalization at $t\to\infty$ in thermodynamic observables Anders and Schiller (2006); Rosch (2012); Nghiem and Costi (2014a); Nghiem et al. (2016), and in spectral functions Anders (2008b); Weymann et al. (2015), but it suffices for a consistent nonperturbative picture of the overall time evolution of $A(\omega,t)$.

*Model and quenches.—*We consider the time-dependent Anderson impurity model, $H=\sum_{\sigma}\varepsilon_{d}(t)n_{d\sigma}+U(t)n_{d\uparrow}n_{d\downarrow}+\sum_{k\sigma}\epsilon_{k}c^{\dagger}_{k\sigma}c_{k\sigma}+\sum_{k\sigma}V(c^{\dagger}_{k\sigma}d_{\sigma}+d^{\dagger}_{\sigma}c_{k\sigma})$, where $\varepsilon_{d}(t)=\theta(-t)\varepsilon_{i}+\theta(t)\varepsilon_{f}$ is the energy of the local level, $U(t)=\theta(-t)U_{i}+\theta(t)U_{f}$ is the local Coulomb interaction, $\sigma$ labels the spin, $n_{d\sigma}=d_{\sigma}^{\dagger}d_{\sigma}$ is the number operator for local electrons with spin $\sigma$, and $\varepsilon_{k}$ is the kinetic energy of the conduction electrons with constant density of states $\rho(\omega)=\sum_{k}\delta(\omega-\varepsilon_{k})=1/2D$ with $D=1$ the half-bandwidth. We take $\Gamma\equiv\pi\rho(0)V^{2}=0.001$ throughout and consider two types of quench [referred to subsequently as quench (A) or quench (B)]: (A), from a symmetric Kondo regime with $\varepsilon_{i}=-15\Gamma,U_{i}=30\Gamma$ and a vanishingly small Kondo scale $T^{i}_{\rm K}=3\times 10^{-8}$ 222We use the Bethe ansatz expression $T_{\rm K}=\sqrt{\Gamma U/2}e^{-\pi U/8\Gamma+\pi\Gamma/2U}$ valid in the symmetric Kondo limit $U/\pi\Gamma\gg 1$ Zlatić and Horvatić (1983); Hewson (1993). to a symmetric Kondo regime with $\varepsilon_{f}=-6\Gamma,U_{f}=12\Gamma$ and a larger Kondo scale $T_{\rm K}=2.5\times 10^{-5}$, and, (B), from a mixed valence regime with $\varepsilon_{i}=-\Gamma,U_{i}=8\Gamma$ to a symmetric Kondo regime with $\varepsilon_{f}=-4\Gamma,U_{f}=8\Gamma$ and a Kondo scale $T_{\rm K}=1.0\times 10^{-4}$.

Spectral function $A(\omega,t)$.— We obtain the time-dependent spectral function via $A(\omega,t)=-\frac{1}{\pi}{\rm Im}[G(\omega+i\eta,t)]$, where $G(\omega+i\eta,t)$, with infinitesimal $\eta>0$, is the Fourier transform of $G(t+t^{\prime},t)\equiv-i\theta(t^{\prime})\langle[d_{\sigma}(t+t^{\prime}),d_{\sigma}^{\dagger}(t)]_{+}\rangle_{\hat{\rho}}$ with respect to the relative time $t^{\prime}$ and $\hat{\rho}$ denotes the full density matrix of the initial state Weichselbaum and von Delft (2007); Peters et al. (2006); Costi and Zlatić (2010). In the notation of Ref. Nghiem and Costi, 2014a, we find for the case of positive times 333See Supplementary Material [URL] for derivations and additional results, including Ref. Bulla et al., 2001.

[TABLE]

where $B=d_{\sigma}$, $C=d_{\sigma}^{\dagger}$, and $\rho_{sr}^{i\to f}(m)=\sum_{e}{{}_{f}}\langle sem|\hat{\rho}|rem\rangle_{f}$ is the full reduced density matrix projected onto the final states Nghiem and Costi (2014a). A somewhat more complicated expression can be derived for negative times Note (3). From Eq. (1), we see that the spectral function can be calculated highly efficiently at all times and frequencies from a knowledge of $\rho_{sr}^{i\to f}(m)$, the final state matrix elements, and excitations at each shell $m$. Our expressions for $A(\omega,t)$ in the two time domains $t<0$ and $t>0$ recover the initial and final state spectral functions for $t\to-\infty$ and $t\to+\infty$, respectively and satisfy the spectral sum rule $\int_{-\infty}^{+\infty}d\omega A(\omega,t)=1$ exactly Note (3). Below, we shall first focus on positive times, where the main time evolution of the Kondo resonance occurs, then on negative to positive times, showing how the high energy final state features in $A(\omega,t)$ evolve from their initial state counterparts already at negative times.

Results for positive times.— Consider quench (A), i.e., switching between symmetric Kondo regimes with $T_{\rm K}^{i}\ll T_{\rm K}$. Figure 1(a) shows the overall time-dependence of the spectral function $A(\omega>0,t>0)=A(-\omega,t>0)$. Two structures, associated with two energy scales, are visible at all times $t>0$ : the satellite peak at $\omega=\varepsilon_{f}+U_{f}\approx 240T_{\rm K}$ and a structure on the scale of $T_{\rm K}$ around the Fermi level. The former has negligible time-dependence, indicating that the satellite peak in the spectral function has already formed by time $t=0$ (its evolution at negative times from the initial state satellite peak at $\omega=\varepsilon_{i}+U_{i}>\varepsilon_{f}+U_{f}$ is discussed below). In contrast to this, the structure around the Fermi level has significant time-dependence at $t>0$ and evolves into the fully formed final state Kondo resonance only on time scales $t\gtrsim 1/T_{\rm K}$ [Figs. 1(c) and 1(d)] in agreement with Ref. Nordlander et al., 1999 for the $U=\infty$ Anderson model. For $tT_{\rm K}\gg 1$, the height of the Kondo resonance at the Fermi level approaches its unitary value given by the Friedel sum rule $\pi\Gamma A(\omega=0,t\to\infty)=1$ to within $15\%$ [Fig. 1(d)]. The small deviation from the expected value is a result of incomplete thermalization due to the discretized Wilson chain used within TDNRG Rosch (2012); Weymann et al. (2015); Note (3). Consequently, evaluating $A(\omega,t\to\infty)$ via the self-energy Bulla et al. (1998) does not improve the Friedel sum rule further in this limit Anders (2008b). In the opposite limit, $t\to-\infty$, where thermalization is not an issue, we recover the Friedel sum rule to within $3\%$ (discussed below). The use of a discrete Wilson chain is also the origin of the small substructures at $|\omega|\lesssim T_{\rm K}$ in Figs. 1(b)-1(d), effects seen in the time evolution of other quantities, such as the local occupation, and explained in terms of the discrete Wilson chain Eidelstein et al. (2012). On shorter time scales, $tT_{\rm K}\lesssim 1$, states in the region $T_{\rm K}^{i}\ll|\omega|<T_{\rm K}$, initially missing [Fig. 1(b)], are gradually filled in by a transfer of spectral weight from higher energies [Fig. 1(c)] to form the final state Kondo resonance at long times [Fig. 1(d)]. The presence of a structure on the final state Kondo scale $T_{\rm K}$ at short times $t\to 0^{+}$ is understood as follows: the Fourier transform with respect to $t^{\prime}=t_{1}-t_{2}$ necessarily convolutes information about the final state at large $t_{1},t_{2}$ into the spectral function at short-times $t$ Turkowski and Freericks (2005). Hence, the gross features of the spectral function, even at short times $t\to 0^{+}$, are close to those of the final state spectral function $A(\omega,t\to\infty)$, and far from those of the initial state spectral function. Clear signatures of the latter, such as the much narrower initial state Kondo peak, only appear at negative times.

Consider now quench (B), in which the system, is switched from the mixed valence to the symmetric Kondo regime. Figures 2(a)-(b) show the overall time-dependence of the spectral function for $\omega<0$ [Figure 2(a)] and $\omega>0$ [Figure 2(b)]. As for quench (A), two structures associated with two energy scales are again visible at all times $t>0$: the satellite peaks at $\omega=\varepsilon_{f}\approx-40T_{\rm K}$ [Figure 2(a)] and $\omega=\varepsilon_{f}+U\approx+40T_{\rm K}$ [Figure 2(b)] and a structure on the scale of $T_{\rm K}$ around the Fermi level [Figs. 2(a) and 2(b)]. In contrast to quench (A), the former have some non-negligible time-dependence at short positive times as can be seen in Fig. 2(c) for $tT_{\rm K}=10^{-4}$ ($t\Gamma=10^{-3}$), where the weight of the satellite peaks has still not equalized. This asymmetry vanishes on time scales exceeding $1/\Gamma$ [Figs. 2(d) and 2(e) for $tT_{\rm K}=1$ ($t\Gamma=10$) and $tT_{\rm K}=10^{4}$ ($t\Gamma=10^{3}$), respectively]. The low energy structure of width $T_{\rm K}$, initially asymmetric and exceeding the unitary height $1/\pi\Gamma$, has significant time-dependence for $t>0$ and evolves into the fully developed Kondo resonance at $t\gtrsim 1/T_{\rm K}$ [Figs. 2(d) and 2(e)]. The deviation from the Friedel sum rule $\pi\Gamma A(\omega=0,t\to\infty)=1$ is comparable to that for quench (A) and reflects the incomplete thermalization due to the discrete Wilson chain used within TDNRG. The discrete Wilson chain also results in the substructures at $|\omega|\lesssim T_{\rm K}$ in Figs. 2(c) and 2(d) and in the small remaining asymmetry of the fully developed Kondo resonance in Fig. 2(e).

From negative to positive times.— Figures 3(a) and 3(b) show the overall time-dependence of the spectral function for negative and positive times, respectively, for quench (A), on a linear frequency scale. As for positive times [Fig. 1(a) and Fig. 3(b)], low and high energy structures are visible also for negative times [Fig. 3(a)]. Moreover, it is clear from Figs. 3(a) and 3(b) that the transition from the initial to the final state spectral function occurs on different time scales for the different structures. Consider first the high energy structures, which carry essentially all the spectral weight. Initially, these are located at $\omega=\pm\varepsilon_{i}\approx\pm 600T_{\rm K}$ as is clearly visible in Fig. 3(a) or in Fig. 3(c) for $tT_{\rm K}=-10^{3}$ ($t\Gamma=-4\times 10^{4}\ll-1$). They cross over to their final state positions at $\omega=\pm\varepsilon_{f}=\pm 240T_{\rm K}$ when $tT_{\rm K}\gtrsim-10^{-2}$ ($t\Gamma\gtrsim-0.4$) [Figs. 3(a) and 3(e)], i.e., on the charge fluctuation time scale $1/\Gamma$. This can also be seen in Figs. 3(d) and 3(e). This large shift in spectral weight from $\omega=\pm\varepsilon_{i}$ to $\omega=\pm\varepsilon_{f}$ in the time-range $-10^{-2}\lesssim tT_{\rm K}\lesssim-10^{-3}$ ($-0.4\lesssim t\Gamma\lesssim-0.04$), clearly seen in Fig. 3(a), is accompanied by small regions of negative spectral weight in this transient time range Note (3). This does not violate any exact results for time-dependent, as opposed to steady-state, spectral functions, and is observed in other systems Dirks et al. (2013); Freericks and Turkowski (2009); Jauho et al. (1994). The spectral sum rule is satisfied analytically exactly at all times and numerically within $1\%$ at all negative times and to higher accuracy at positive times for all quench protocols Note (3). Turning now to the low energy structure, i.e., the Kondo resonance, the use of a linear frequency scale now allows the initial state Kondo resonance at $\omega=0$ to be clearly seen in Fig. 3(a) [see also Fig. 3(c)]. This structure, of width $T_{\rm K}^{i}\ll T_{\rm K}$ at $t\to-\infty$ and satisfying the Friedel sum rule $\pi\Gamma A(\omega=0,t\to-\infty)=1$, gradually broadens and acquires a width of $T_{\rm K}$ at short negative times Note (3), and then evolves into the fully developed Kondo resonance on positive time scales $tT_{\rm K}\gtrsim 1$ [Fig. 3(e)].

Even more interesting is the negative [Fig. 4(a)] to positive [Fig. 4(b)] time evolution of the spectral function upon quenching from the mixed valence to the symmetric Kondo regime [quench (B)]. At large negative times [Fig. 4(c)], one recovers the initial state spectral function of the mixed valence regime (with $\varepsilon_{i}=-\Gamma$) showing a mixed valence resonance, renormalized by many-body effects to lie close to, but just above the Fermi level $\varepsilon_{i}\to\tilde{\varepsilon}_{i}\gtrsim 0$ and satisfying the Friedel sum rule $A(0,t\to-\infty)=\sin^{2}(\pi n_{d}/2)/\pi\Gamma$ to within $3\%$ Costi et al. (1996) [Figs. 4(a) and 4(c), $n_{d}=0.675$]. The upper satellite peak at $\omega=\varepsilon_{i}+U_{i}=7\Gamma\approx 70T_{\rm K}$ is more clearly visible in Fig. 4(c). These peaks give rise to the final state satellite peaks at $\omega=\pm\varepsilon_{f}=\pm 4\Gamma\approx\pm 40T_{\rm K}$ which start to form already at negative times $tT_{\rm K}\gtrsim-10^{-1}$ ($t\Gamma\gtrsim-1$), i.e., on the charge fluctuation time scale $1/\Gamma$, as for quench (A). While the positions of these peaks start to shift to their final state values at negative times $tT_{\rm K}\gtrsim-10^{-1}$ ($t\Gamma\gtrsim-1$), their weights remain disparate [see Fig. 4(e)] and only equalize at $tT_{\rm K}\gtrsim+10^{-1}$ ($t\Gamma\gtrsim+1$) as clearly seen in Fig. 4(b), i.e., the formation of the high energy final state satellite peaks occurs on a fast time scale $t\approx 1/\Gamma$ in the interval $-1/\Gamma\lesssim t\lesssim+1/\Gamma$ (dashed lines in Fig.4). Going into more details, we see in Figs. 4(a) and 4(c)-(e) the deconstruction of the mixed valence resonance in the time range $-1/\Gamma<t<0$. While this resonance carries essentially all the spectral weight at $t\ll-1/\Gamma$, weight is gradually transferred to $\omega<0$, with precursor oscillations starting at $tT_{\rm K}=-1$ ($t\Gamma=-10$) [Fig. 4(d)], to form the lower final state satellite peak at $\omega=\varepsilon_{f}$ for $-1/\Gamma<t<0$ [Fig. 4(e)] . Simultaneously, the mixed valence resonance narrows from its original width $\Gamma\approx 10T_{\rm K}$ and shifts towards the Fermi level to form a low energy structure on the scale of $T_{\rm K}$ [Fig. 4(e)]. The latter eventually evolves into the final state Kondo resonance at $tT_{\rm K}\gtrsim 1$. The final state spectral function is recovered in the long-time limit $tT_{\rm K}\gg 1$ [Fig. 4(f)].

Conclusions.— In summary, we investigated within the TDNRG the time evolution of the spectral function of the Anderson impurity model in the strong correlation limit. Quenching into a Kondo correlated final state, we showed that the Kondo resonance in the zero temperature spectral function only fully develops at very long times $t\gtrsim 1/T_{\rm K}$, although a preformed version of it is evident even at very short times $t\to 0^{+}$. The latter can be used as a smoking gun signature of the transient build up of the Kondo resonance in future cold atom realizations of the Anderson impurity model Bauer et al. (2013). The satellite peaks evolve from their initial state values at negative times on a much faster time scale $t\approx 1/\Gamma$ in the time-interval $-1/\Gamma\lesssim t\lesssim 1/\Gamma$. Our formulation of sum rule conserving two-time nonequilibrium Green functions within TDNRG, including lesser Green functions, and their explicit dependence on both times Note (3), yields the basic information required for applications to time-dependent quantum transport Jauho et al. (1994) and constitutes a first step towards using TDNRG as an impurity solver within nonequilibrium DMFT Freericks et al. (2006); Gramsch et al. (2013); Aoki et al. (2014).

I Retarded two-time nonequilibrium Green function in TDNRG

We consider the retarded two-time Green function $G_{BC}(t+t^{\prime},t)=-i\theta(t^{\prime})\operatorname{Tr}\{\hat{\rho}[\hat{B}(t+t^{\prime}),\hat{C}(t)]_{s}\}$ for a system undergoing a quantum quench at $t=0$ as described by the time-dependent Hamiltonian $H(t)=(1-\theta(t))H_{i}+\theta(t)H_{f}$, with $\hat{\rho}=e^{-\beta H_{i}}$ the density matrix of the initial state, represented by the full density matrix in Eq. (S18) Weichselbaum and von Delft (2007). Since $t^{\prime}>0$, we have two cases to consider, (i), $t>0$, in which case both operators $B$ and $C$ evolve with respect to $H_{f}$, and, (ii), $t<0$, in which case, either both operators evolve with respect to $H_{i}$ if $t<t+t^{\prime}<0$, or if $t<0<t+t^{\prime}$ operator $C$ evolves with respect to $H_{i}$, while operator $B$ evolves with respect to $H_{f}$. While the Fourier transform with respect to the time difference $t^{\prime}$ yielding $G(\omega+i\eta,t)$ and hence $A(\omega,t)$ is straightforward in case (i), in case (ii), expressions for $G(t+t^{\prime},t)$ are needed from both time domains $t<t+t^{\prime}<0$ and $t<0<t+t^{\prime}$ in order to construct $G(\omega+i\eta,t<0)$ and hence $A(\omega,t<0)$. Depending on the physical system considered, both positive and negative times may be of interest. Thus, in problems where the quench represents an initial state preparation, the main interest is in the evolution at $t>0$ following this preparation Lobaskin and Kehrein (2005). However, if the quench is considered to be a perturbation applied to the system at time $t=0$, the full time evolution is of interest. This case is also required for applications to nonequilibrium dynamical field theory (DMFT)Freericks et al. (2006); Turkowski and Freericks (2005); Aoki et al. (2014). From a theoretical point of view, $t<0$ is also of interest to fully describe the evolution of the spectral function $A(\omega,t)$ from its initial state value at $t=-\infty$ to its final state value at $t=+\infty$.

I.1 Positive time-dependence $t>0$

We first consider the case $t>0$, treating $t<0$ in the next subsection. We have for the retarded Green function,

[TABLE]

where $s=\pm 1$ for fermionic/bosonic Green functions, respectively. In the notation of Ref. Nghiem and Costi, 2014a and following the approach of Anders Anders (2008), we have for the first ($BC$) term of the anticommutator with $t>0$

[TABLE]

where $|lem\rangle,m=m_{0},\dots,N$ is the complete set of eliminated states in the NRG diagonalization procedure, with $m_{0}$ the first iteration at which states are eliminated and $N$ is the last NRG iteration (see Refs. Anders and Schiller (2006); Nghiem and Costi (2014a) for details). We evaluate ${{}_{f}}\langle re^{\prime}m|\hat{B}(t^{\prime})\hat{C}|sem\rangle_{f}$ by inserting the decomposition of unity Anders and Schiller (2006) $I=I_{m}^{+}+I_{m}^{-}$ between $\hat{B}(t^{\prime})$ and $\hat{C}$,

[TABLE]

The first term in the above expression is diagonal in the environment variables $e$, $e^{\prime}$ , and $q$ runs over all states (kept and discarded) at the shell $m$. We put $1=1_{m^{\prime}}^{+}+1_{m^{\prime}}^{-}$ in the last term to obtain

[TABLE]

Substituting (S4) into (S3) and using the NRG approximation, we have

[TABLE]

Substituting $\quad{{}_{f}}\langle re^{\prime}m|k_{1}e^{\prime\prime}m^{\prime}\rangle_{f}=\delta_{e^{\prime}_{m}e^{\prime\prime}_{m}}[A^{\alpha_{m}\dagger}_{XK}...A^{\alpha_{m^{\prime}+1}\dagger}_{KK}]_{rk_{1}}$ and $\quad{{}_{f}}\langle k_{2}e^{\prime\prime}m^{\prime}|sem\rangle_{f}=\delta_{e^{\prime\prime}_{m}e_{m}}[A^{\alpha_{m^{\prime}+1}}_{KK}...A^{\alpha_{m}}_{KX^{\prime}}]_{k_{2}s}$ into Eq. (S5), results in the following expression for Eq. (S2)

[TABLE]

in which $\rho^{i\to f}_{s,r}(m)=\sum_{e}{{}_{f}}\langle sem|{\rho}|rem\rangle_{f}$ is the projected full reduced density matrix known from Ref. Nghiem and Costi, 2014a. Fourier transforming the above equation with respect to $t^{\prime}$ gives

[TABLE]

with $\eta$ is a positive infinitesimal. Similarly, the second ($CB$) term of the anticommutator in (S1) gives us

[TABLE]

Hence, for positive times we obtain $G(\omega,t)=I_{1}(\omega,t)+I_{2}(\omega,t)$.

In order to calculate the time-dependent spectral function from the above, one can follow Anders Anders (2008) by defining the following time-dependent density matrix

[TABLE]

Then we have the following expression for the Green’s function

[TABLE]

from which the time-dependent spectral function can be calculated. This expression of Anders Anders (2008), generalized here within the full density matrix approach, and hence valid at arbitrary temperature Nghiem and Costi (2014a), requires, however, the time-dependent reduced density matrix $\tilde{\rho}^{i\to f}_{k_{2}k_{1}}(m^{\prime},t)$ at each point in time, and the latter is in turn obtained via the recursion relation in Eq. (S9), resulting in a numerically highly time consuming calculation. This motivated us to develop an alternative and numerically more tractable expression for the retarded two-time Green function, to which we now turn.

A different and numerically more feasible expression for $G(\omega,t>0)$ than the expression above, can be obtained if we go back to Eq. (S5) and substitute this into Eq. (S2),

[TABLE]

The first term in this expression is simply

[TABLE]

and the second term can be rearranged as follows

[TABLE]

Combining the above two terms, we have

[TABLE]

in which $\sum_{e}{{}_{f}}\langle sem|\hat{\rho}(t)|rem\rangle_{f}=e^{i(E^{m}_{r}-E^{m}_{s})t}\rho^{i\to f}_{sr}(m)$ and all the other matrix elements are known quantities. Together with the second term ($CB$) in the anticommutator,

[TABLE]

we have the retarded Green’s function as follows

[TABLE]

This expression is useful for the non-equilibrium DMFT, which requires the dynamical fields expressed in terms of two-time Green functions Freericks et al. (2006); Gramsch et al. (2013); Aoki et al. (2014). Fourier transforming with respect to the time difference $t^{\prime}$ gives

[TABLE]

which together with $\sum_{e}{{}_{f}}\langle sem|\hat{\rho}(t)|rem\rangle_{f}=\rho^{i\to f}_{sr}(m)e^{-i(E_{s}^{m}-E_{r}^{m})t}$ results in a time-dependent spectral function $\displaystyle A(\omega,t)=-{\rm Im}[G(\omega,t)]/\pi$ that can be evaluated at all times in a numerically highly efficient manner: only the time-independent projected density matrix $\rho^{i\to f}_{sr}(m)$ together with the NRG excitations and matrix elements are required to evaluate $A(\omega,t)$ at all positive times. While the above expression looks deceptively similar to the first term in Anders expression for the positive time Green function in Eq. (S10), this is not the case [note the different sum $\sum_{rsq}^{\notin KK^{\prime}K^{\prime\prime}}$ in Eq. (S17) as compared to the sum $\sum_{rs}^{\notin KK^{\prime}}\dots\times\sum_{q\in(K^{\prime\prime}D)}$ in the first term of Eq. (S10)]. It includes also the second term in Eq. (S10), involving $\tilde{\rho}^{i\to f}_{k_{2}k_{1}}(m^{\prime},t)$, but within a different approximation that follows from the second line of Eq. (S14). In this way, the recursive evaluation of a reduced density matrix depending explicitly on time is circumvented in our expression, making it numerically tractable. A similar expression to Eq. (S17) has been derived for initial state density matrices corresponding to either pure states, such as the ground state, or to decoupled initial states ($\Gamma=0$ in the Anderson model) with or without excitations in the bath and used to study thermalization in the Anderson model following such an initial state preparation Weymann et al. (2015).

As a check on Eq. (S17), we can verify that it reduces to the equilibrium retarded Green function in the case that $H^{f}=H^{i}$ (vanishing quench size). In this limit, using the definition of the full density matrix of the initial stateWeichselbaum and von Delft (2007),

[TABLE]

and of the reduced full density matrix $R^{m}_{red}(k,k^{\prime})$ in Refs. Costi and Zlatić, 2010; Nghiem and Costi, 2014a, we have

[TABLE]

Substituting this expression into Eq. (S17) for $G(\omega,t)$, we obtain the following time-independent expression

[TABLE]

which is identical to the expression in the equilibrium case Weichselbaum and von Delft (2007); Costi and Zlatić (2010); Nghiem and Costi (2014a).

It is interesting to compare the spectral function obtained from our Eq. (S17) with that of Anders obtained from Eq. (S10). For this purpose, we consider quench (A) of the main text and three representative times: $t=0$, $t=1/T_{\rm K}$, and $t\to+\infty$. From Figs. S1(a) and S1(c) we see that both expressions give the same results at $t=0$ and $t\to+\infty$, while at finite times small differences arise for $|\omega|\lesssim T_{\rm K}$ [see Fig. S1(b)]. As we discussed in connection with Eq. (S10) above the main advantage of our new expression Eq. (S17), is that it can be evaluated numerically very efficiently. By contrast, to date, no results at finite times have been published using Eq. (S10) and only results for infinite times have been published Anders (2008).

I.2 Negative time-dependence ($t<0$)

At negative times $t<0$, we can have $t+t^{\prime}$ before or after the quench. When $t+t^{\prime}<0$ or $t<-t^{\prime}$, we have

[TABLE]

which is $t$-independent, and just corresponds to the equilibrium propagator of the initial state Hamiltonian (as long as $t<-t^{\prime}$). While for $t+t^{\prime}>0$ or $t>-t^{\prime}$, we have

[TABLE]

In general, then, we have for the retarded Green function at $t<0$

[TABLE]

For the first part ($BC$ term) of the anticommutator at $t+t^{\prime}<0$ we have

[TABLE]

Next, we consider $\sum_{e}[C\rho]^{m}_{sr}$, which with Eq. (S18) becomes

[TABLE]

In this sum, only the parts with $m_{1}\geq m$ are finite, and we obtain

[TABLE]

with $R_{kk^{\prime}}^{m},\tilde{R}^{m}_{qr}$ as in Ref. Nghiem and Costi, 2014a. Substituting the above into Eq. (S23) we obtain

[TABLE]

Similarly, the first part ($BC$ term) of the anticommutator for $t+t^{\prime}>0$ is given by

[TABLE]

where $\widetilde{C\rho}(t)\equiv e^{iH_{i}t}C\rho e^{-iH_{i}t}$ (i.e., the tilde signifies that the time evolution operators apply to the composite operator $C\rho$). We decompose the above sum into three parts corresponding to $m_{1}>m$, $m_{1}=m$, and $m_{1}<m$, then simplify them as follows

[TABLE]

Substituting $\sum_{e}[C\rho]^{m}_{sr}$ into the above equation, we have

[TABLE]

The second part ($CB$ term) of the anticommutators for the two parts at $t+t^{\prime}<0$ and $t+t^{\prime}>0$, denoted by $I^{-}_{2}(t+t^{\prime},t^{\prime})$ and $I^{+}_{2}(t+t^{\prime},t^{\prime})$, respectively, can be derived in a similar way and combining all terms gives the expression for the two-time Green function $G(t^{\prime}+t,t)$ at $t<0$. Fourier transforming with respect to the time difference $t^{\prime}$ results in the following expression for the spectral function at $t<0$

[TABLE]

From this expression, we easily see that $G(\omega,t\to-\infty)$ recovers the initial state Green function. While from the starting definition, we have $G(\omega,t\to 0^{-})$ equal to $G(\omega,t\to 0^{+})$, in the approximate expressions above this is no longer guaranteed as a consequence of the NRG approximation and the different derivations for $t\to 0^{\pm}$. Figures S2(a)-S2(b) show the spectral functions at $t\to 0^{+}$ and $t\to 0^{-}$ and quantify the size of the discontinuity at $t=0$. While the spectral functions for $t\to 0^{\pm}$ match to high accuracy at high frequencies $|\omega|\gg T_{\rm K}$ for both quench (A) [Fig. S2(a)] and quench (B) [Fig. S2(b)] of the main text, there is a mismatch at low frequencies.

In Fig. S3 we show the negative time spectral function for quench (A) of the main text using logarithmic axes for both time $t$ (in units of the initial state Kondo scale $T_{\rm K}^{i}=1.2\times 10^{-3}T_{\rm K}$) and frequency, $\omega$ (in units of the final state Kondo scale $T_{\rm K}$). The data is the same as the positive frequency data in Fig. 3(a) of the main text, but the use of a logarithmic frequency axis now makes clearer the statement made there concerning the structure around the Fermi level, that “this structure, of width $T_{\rm K}^{i}\ll T_{\rm K}$ at $t\to-\infty$ and satisfying the Friedel sum rule $\pi\Gamma A(\omega=0,t\to-\infty)=1$, gradually broadens and acquires a width of $T_{\rm K}$ at short negative times …”. In addition, Fig. S3 shows that the low energy structure on a scale $T_{\rm K}$ at short negative times is formed by drawing spectral weight from both the initial state Kondo resonance (diagonal stripes at $t\gtrsim-1/T_{\rm K}^{i}$), and also from the satellite peaks (starting at $t\gtrsim-1/\Gamma$). As for short positive times [Fig. 1(a) of the main text], the “preformed” Kondo resonance at short negative times ($t\to 0^{-}$) is seen to have missing states in the vicinity ($\omega\ll T_{\rm K}$) of the Fermi level. Finally, notice that the relevant time scale for the “devolution” of the initial state Kondo resonance at $t=-\infty$ is $1/T_{\rm K}^{i}=1.2\times 10^{-3}T_{\rm K}$.

I.3 Numerical evaluation of $A(\omega,t)$

I.3.1 Positive time spectral function

We first consider the evaluation of the spectral function for positive times. From (S17), we have for $A(\omega,t)=-\frac{1}{\pi}{\rm Im}G(\omega+i\eta,t)$

[TABLE]

where $E_{sr}=E_{s}-E_{r}$. The first contribution,

[TABLE]

is evaluated in the usual way by replacing $\delta(\omega-E)$ by the logarithmic Gaussian $\frac{1}{\sqrt{\pi}b|E|}e^{-b^{2}/4}e^{-(\frac{\ln(|w/E|)}{b})^{2}}$ and summing over the excitations $E$ Bulla et al. (2001, 2008).

In order to evaluate the second contribution above, we define the auxiliary function $F^{\prime\prime}(\omega,t)$ via

[TABLE]

and evaluate this in the usual way. Taking it’s principle value integral then gives the second contribution to the spectral function:

[TABLE]

To summarize,

[TABLE]

I.3.2 Negative time spectral function

We now consider the evaluation of $A(\omega,t)$ for negative times starting from the expression for the retarded Green function in Eq. (S30). This is a sum of two terms $G(\omega+i\eta,t)=G_{1}(\omega+i\eta,t)+G_{2}(\omega+i\eta,t)$, with

[TABLE]

where $\{C\tilde{R}\}_{sr}^{m}\equiv\sum_{q}(C^{m}_{sq}\tilde{R}^{m}_{qr}+\tilde{R}^{m}_{sq}C^{m}_{qr})$. Correspondingly, the spectral function is also written as a sum of two parts $A(\omega,t)=A_{1}(\omega,t)+A_{2}(\omega,t)$. Consider first $A_{1}(\omega,t)=-\frac{1}{\pi}{\rm Im}G_{1}(\omega,t)$, for finite $0>t>-\infty$, $G_{1}(\omega+i\eta,t)$ is regular, having no poles on the real axis, then $A_{1}(\omega,t)$ is evaluated directly

[TABLE]

with a finite $\eta=b|E_{sr}^{m}|$ and $b\geq 1/N_{z}$ where $N_{z}$ is the number of bath realizations in the $z$ averaging procedure.

For consistency, we also evaluate $G_{2}(\omega,t)$ in Eq. (S36) with the same Lorentzian broadening and thereby obtain $A_{2}(\omega,t)=-\frac{1}{\pi}{\rm Im}G_{2}(\omega,t)$ and thus $A(\omega,t)=A_{1}(\omega,t)+A_{2}(\omega,t)$.

I.4 Spectral sum rule

The spectral weight sum rule $\int_{-\infty}^{+\infty}d\omega A(\omega,t)=1$ with $A(\omega,t)=-\frac{1}{\pi}{\rm Im}G(\omega,t)$ is exactly satisfied for all positive times within our expression Eq. (S17) and the same holds true for Eq. (S10) (see Ref. Anders, 2008). Here, we prove that it is exactly satisfied also for all negative times. With $G(\omega,t)$ defined by Eq. (S30), we need to evaluate the following terms

[TABLE]

with $\eta$ a positive infinitesimal. Defining $E^{m}_{r}-E^{m}_{s}=E^{m}_{rs}$ and $E^{m}_{r_{1}}-E^{m}_{s_{1}}=E^{m}_{r_{1}s_{1}}$, we have

[TABLE]

where use was made of $\int_{-\infty}^{+\infty}\cos(tx)/(x^{2}+a^{2})=\pi e^{ta}/a$ and $\int_{-\infty}^{+\infty}x\sin(tx)/(x^{2}+a^{2})=-\pi e^{ta}$ ($t<0$). Using Eqs. (S40-S41) and Eq. (S30), we have

[TABLE]

Therefore the sum rule is proved at $t<0$. The sum rule also holds for $t=0^{-}$, as can be seen by noting that the last integrals contributing to $I_{1}(t)$ and $I_{2}(t)$ return [math] in this case.

The numerically evaluated negative-time spectral function also satisfies this sum rule to within $\approx 1\%$ for all quench protocols and all negative times, as shown in Fig. S4(a). In Fig. S4(b) we show $w_{-}(t)=\int_{-\infty}^{+\infty}d\omega F(\omega,t)$ with $F(\omega,t)=A(\omega,t)$ if $A(\omega,t)<0$ and $F(\omega,t)=0$ if $A(\omega,t)>0$, i.e., the contribution to the total weight coming from regions of negative spectral density. Regions of negative spectral weight appear for transient times in many other systems Dirks et al. (2013); Jauho et al. (1994); Freericks and Turkowski (2009), while in steady state limits $t\to\pm\infty$ the spectral function is positive definite within canonical density matrix approaches. We see that both $w_{-}(t)$ and the error in the sum rule are both largest in the time region $t\Gamma\gtrsim-1$ where the major part of the spectral weight (located in the satellite peaks) is being rearranged from $\omega=\varepsilon_{i},\varepsilon_{i}+U_{i}$ to $\omega=\varepsilon_{f},\varepsilon_{f}+U_{f}$. The maximum size of $w_{-}(t)$ correlates with the quench size $\Delta\epsilon_{d}=|\epsilon_{f}-\epsilon_{i}|$, and reaches up to $15\%$ for the largest quench (A) in Fig. S4(b).

I.5 Friedel sum rule, thermalization and discretization effects

The TDNRG replaces the continuum conduction electron bath $H=\sum_{k\sigma}\epsilon_{k}c^{\dagger}_{k\sigma}c_{k\sigma}$ in the Anderson model by a logarithmically discretized bath $\epsilon_{k}\to\pm\Lambda^{-n}$ whose tight binding representation in energy space (the so called Wilson chain) is given by $H=\sum_{n=o,\sigma}^{N}t_{n}(f_{n\sigma}^{\dagger}f_{n+1\sigma}+H.c.)$ with $t_{n}\approx\Lambda^{-(n-1)/2}$ for $n\gg 1$ and where $\Lambda>1$ is the discretization parameter. The continuum limit corresponds to $\Lambda\to 1^{+}$, which is not possible to take numerically within the iterative NRG diagonalization scheme due to the increasing slow convergence in this limit. The effect of using such a discrete Wilson chain on the time evolution of physical quantities in response to a quench is twofold: (i), incomplete thermalization at long times after the quench due to the fact that a Wilson chain (even in the limit $N\to\infty$) cannot act as a proper heat bath Rosch (2012) (see also Refs. Anders and Schiller, 2006; Nghiem and Costi, 2014a, b; Weymann et al., 2015) and ,(ii), additional real features appear in the time evolution, due to the logarithmic discretization of the bath. We address these two effects in turn.

Incomplete thermalization in the long time limit is reflected in deviations of observables from their expected values in the final state. These deviations can be reduced by decreasing $\Lambda$, as shown in Ref. Nghiem and Costi, 2014a for the case of the local occupation $n_{d}(t)$ in the Anderson model, where $n_{d}(t\to+\infty)$ was found to approach its expected value more closely upon decreasing $\Lambda$ towards $1$. In the present context of spectral functions, incomplete thermalization is reflected in a $\approx 15\%$ deviation of $\pi\Gamma A(\omega=0,t\to\infty)$ from the continuum result $\sin^{2}(\pi n_{d}/2)$, where $n_{d}$ is the occupation number in the final state. However, this deviation is not an error of the TDNRG but is expected because of (i). The Friedel sum rule (FSR) only holds for equilibrium states. In the infinite past, where thermalization issues play no role, and one achieves the equilibrium state at $t=-\infty$, the FSR is satisfied to within a $1-3\%$ for all quenches studied, which is comparable to the accuracy achievable in equilibrium NRG calculations Costi et al. (1996) (the spectra for large negative $t$ in Fig.3(c) and Fig.4(c) of the main text). This demonstrates that the $15\%$ deviation in the value of $\pi\Gamma A(\omega=0,t\to\infty)$ is not an error of the TDNRG but is the correct result for the Wilson chain used.

We note, furthermore, that other approaches to time dependent spectral densities, such as the non-crossing approximation Nordlander et al. (1999) suffer, even for equilibrium spectral densities, from much larger errors in the Friedel sum rule (see Table I in Ref. Costi et al., 1996), and still other methods Bock et al. (2016) can result in errors which exceed $50\%$ despite the use of a continuum bath. In the latter approaches the deviations from the Friedel sum rule represent real errors in the underlying methods, whereas in the TDNRG, the deviation observed is that expected from using a logarithmically discretized chain.

The second effect of using a Wilson chain is that additional features appear in the time evolution of physical observables that would be absent for a continuum bath. Examples are the small oscillations seen in Figs .1 and 2 of the main text at low energies $|\omega|<T_{\rm K}$. Physically, these oscillations, or ”substructures” result from the highly nonequilibrium situation created by the quench: following the quench, the local change in energy has to be transported by electrons propagating outwards in the process of thermalization. These electrons are reflected off the different sites ($n=0,1,...,N$) of the inhomogeneous Wilson chain (with hoppings $t_{n}\sim\Lambda^{-(n-1)/2}$) arriving back at the impurity site where they interfere at specific times to give additional features in the time evolution of local quantities (such as in $A(\omega,t)$). This was originally explained in great detail by Eidelstein et al. in Ref. Eidelstein et al., 2012 for observables like the occupation number $n_{d}(t)$. These authors also showed a comparison between $n_{d}(t)$ using TDNRG and $n_{d}(t)$ using exact diagonalization for a noninteracting resonant level model, both calculated using the same Wilson chain. The additional features in the time evolution of $n_{d}(t)$, not present in the continuum model, were identified as real effects and not as being due to errors or unphysical features of the TDNRG method. Thus, the features seen in Figs.1 and 2 of the main text at low energies $|\omega|<T_{\rm K}$ have their origin in the logarithmically discretized bath. In the remote past, before the quench is able to act, such features are expected to be absent. In support of this interpretation, we therefore show on a logarithmic scale the spectral function in the distant past $A(\omega,t=-\infty)$ in Figs. S5(a)-S5b, which indeed show the absence of all “substructures”. This is to be compared with the presence of such ”substructures” at finite and long positive times [Figs. 1(b)-1(d) and Fig. 2(c)-2(e) in the main text]. Finally, we mention that since the limit $\Lambda\to 1^{+}$ is not feasible within NRG, a simpler approach that is useful to obtain results closer to those of the continuum limit is the averaging of time-dependent quantities over $N_{z}\gg 1$ realizations of the bath Anders and Schiller (2006).

II Additional results

Our derivation for the non-equilibrium two-time retarded Green function at both positive and negative times can be used to derive expressions for the other commonly used Green functions in many-body theory, including lesser (and greater) Green functions, and can be used for arbitrary quenches. In the next subsections, we show the results for the cases of the reverse of quench (B) in the main text, i.e. quenching from a symmetric Kondo into a mixed valence regime, a finite temperature quench as in Nordlander et al. Nordlander et al. (1999), and a hybridization quench in which the coupling $\Gamma$ is switched on at time $t=0$. We also show results for the lesser Green function, which together with the retarded Green function constitute the basic ingredients for many applications, e.g., to transient and non-equilibrium transport through correlated quantum dots Hershfield et al. (1992); Meir et al. (1993); Jauho et al. (1994). In Sec. II.5, we show the explicit two-time dependence of the retarded Green function, a basic ingredient in nonequilibrium DMFT applications Freericks et al. (2006).

II.1 Symmetric Kondo to mixed valence [reverse of quench (B)]

In this subsection, we show the time-dependent spectral functions for the case of the reverse of quench (B) in the main text. i.e., from the symmetric Kondo to the mixed valence regime. We see in Fig. S6 how the spectral function evolves from the spectral function of the initial state at $t\to-\infty$ [Fig. S6(c) for $tT_{\rm K}=-10^{3}$ ($t\Gamma=-10^{4}$)], with the Friedel sum rule satisfied to within a few $\%$ in this limit, to its value in the long-time limit with a mixed valence peak close to the Fermi level and a satellite peak at $\omega=\varepsilon_{f}+U$ above the Fermi level [Fig. S6(f) for $tT_{\rm K}=+10^{3}$ ($t\Gamma=+10^{4}$)]. Similar to the other cases in the main text, the initial state satellite peaks at $\omega=\pm\varepsilon_{i}$ rapidly relocate at $t=-1/\Gamma$ [dashed line in Fig. S6(a)] with spectral weight being shifted to form the upper satellite peak and the mixed valence resonance, a process which essentially has completed by $t=+1/\Gamma$ [dashed line in Fig. S6(b)]. The weights of these peaks weights are also close to those of the final state for $t\Gamma\gtrsim+1$. The central peak which represents the Kondo resonance at $t\to-\infty$ also varies strongly at $t\gtrsim-1/\Gamma$, and evolves into the mixed valence peak by time $t\gtrsim+1/\Gamma$.

II.2 Nordlander quench (finite temperature)

In Ref. Nordlander et al., 1999 a quench is made within the $U=\infty$ Anderson model via a level shift $\varepsilon_{d}(t)=\varepsilon_{i}\theta(-t)+\varepsilon_{f}\theta(t)$ with $\epsilon_{i}=-10\Gamma$ ($T_{\rm K}^{i}\approx 10^{-7}$) and $\epsilon_{f}=-4\Gamma$ ($T_{\rm K}\approx 1.8\times 10^{-3}$). This corresponds to a quench from one asymmetric Kondo regime to another with disparate Kondo scales; in contrast we previously investigated quenches in which one of the states was in a symmetric Kondo regime whereas the other was in an mixed valence regime. The quench in Ref. Nordlander et al., 1999 also differs from those studied so far since it is at a finite temperature $T=2.5\times 10^{-3}$ such that $T_{\rm K}^{i}\ll T\approx T_{\rm K}$. Thus, initially the Kondo resonance is strongly temperature suppressed whereas in the final state it is only moderately suppressed by temperature. This quench can therefore serve to illustrate the application of our TDNRG formalism for time dependent spectral functions to finite temperatures.

In figure S7, we show the time-dependent spectral function from negative to positive times. The calculations were carried out for $U\gg D$ to simulate the $U=\infty$ case. We therefore observe only the satellite peak below the Fermi level in the negative frequency range, both in the initial and final states. Similar to the other calculations, the satellite peak rapidly relocates at $t\approx-1/\Gamma$ from $\varepsilon_{i}$ to $\varepsilon_{f}$ as shown in Fig. S7 (a). At the same time, the spectral function develops small regions of negative spectral weight, with the total sum-rule remaining satisfied to within $1\%$ as shown in Fig. S4 (a). The central peak at $\omega=0$ is absent at $t\to-\infty$ since the calculation is at finite temperature $T\gg T_{\rm K}^{i}$. Since the temperature $T\approx T_{\rm K}$ is finite and comparable to the final state Kondo scale, the Kondo resonance does not fully develop at long times [Fig. S7(b) and S7(c)] with $\pi\Gamma A(\omega=0,t\to\infty)$ reaching only about $59\%$ of its $T=0$ value. This is better seen in Fig. S7(d), which shows a close up of the low frequency region around the Fermi level. Nevertheless, despite the finite temperature, one sees the build up of the Kondo resonance at $t\gtrsim 1/T_{\rm K}$.

II.3 Hybridization quench

In this subsection, we show the time-dependent spectral functions for the case of a hybridization quench as in Ref. Weymann et al., 2015, where the hybridization between the impurity and the conduction electrons, initially turned off at $t<0$, is suddenly turned on at $t=0$. In figure S8 (a)-(b), we see that for this quench also, low and high energy features are present at all times. The high energy features correspond to the final state satellite peaks at $\varepsilon_{f}$ and $\varepsilon_{f}+U$, whereas the low energy feature of width on the scale of the final state Kondo temperature $T_{\rm K}$ represents the Kondo resonance. While the former have little temperature dependence at all $t>0$, as in Weymann et al. Weymann et al. (2015), the latter has significant time dependence, developing fully only at $tT_{\rm K}\gtrsim 1$ [Fig. S8(a)] with weight drawn in from higher energies in the process. Notice that this low energy peak appears even at $t=0$, which is different from Weymann et al,Weymann et al. (2015) since the broadening parameter is set to be time-independent in our calculation, while it is time-dependent (and large of order $\Gamma$ at $t=0$) in Weymann et al. Weymann et al. (2015) . While the strong time dependence of the Kondo resonance can be seen on a logarithmic frequency scale from Fig. S8(a), it is barely discernible on the linear frequency scale of Fig. S8(b).

II.4 Lesser Green functions

We consider explicitly the lesser Green function for the local level in the Anderson impurity model, defined by

[TABLE]

For equal times ($t^{\prime}=0$),

[TABLE]

i.e., $\text{Im}[G^{<}(t,t)]=n_{d\sigma}(t)$, so the lesser Green function at equal times gives the time evolution of the local occupation number. Following the derivation for the retarded Green function in Sec. I, we similarly obtain the following expression for the lesser Green function

[TABLE]

with $B\equiv d^{\dagger}_{\sigma}$ and $C\equiv d_{\sigma}$. The time evolution of the occupation number calculated from this expression by setting $t^{\prime}=0^{+}$ can be compared with that calculated directly from the thermodynamic observable $n_{d\sigma}(t)$ Nghiem and Costi (2014a). The two results, shown in Fig. S9, match perfectly at short times and differ slightly on longer time scales ($t\Gamma\gtrsim 1$). This small difference arises because the NRG approximation enters differently in the expressions for thermodynamic and dynamic quantities.

II.5 Retarded Green function: explicit dependence on times

From Eq. (S17), we can directly evaluate the dependence of the retarded Green function on its two time arguments. This is shown for the imaginary and real parts in Figs. S10(a)-S10(b) versus the time difference $t^{\prime}>0$ and time $t>0$. At equal times we see from Fig. S10(a) that $-{\rm Im[G(t,t))]}=1$ for all $t>0$, recovering the canonical anticommutation relation for fermions, ande hence the spectral sum rule for $t>0$. Non-equilibrium DMFT Freericks et al. (2006); Schmidt and Monien (2002); Gramsch et al. (2013); Aoki et al. (2014) requires impurity Green functions in real time, and the ability to calculate these within TDNRG, which we here demonstrated, is a useful first step for future applications to the former.