Ultrafast molecular dynamics in terahertz-STM experiments: Theoretical analysis using the Anderson-Holstein model

We analyze ultrafast tunneling experiments in which electron transport through a localized orbital is induced by a single-cycle terahertz (THz) pulse. We include both electron-electron and electron-phonon interactions on the localized orbital using the Anderson-Holstein model and consider two possible ﬁlling factors, the singly occupied Kondo regime and the doubly occupied regime relevant to recent experiments with a pentacene molecule. Our analysis is based on variational non-Gaussian states and provides the accurate description of the degrees of freedom at very different energies, from the high microscopic energy scales to the Kondo temperature T K . To establish the validity of this method we apply this formalism to study the Anderson model in the Kondo regime in the absence of coupling to phonons. We demonstrate that it correctly reproduces key properties of the model, including the screening of the impurity spin, formation of the resonance at the Fermi energy, and a linear conductance of 2 e 2 / h . We discuss the suppression of the Kondo resonance by the electron-phonon interaction on the impurity site. When analyzing THz-STM experiments we compute the time dependence of the key physical quantities, including current, the number of electrons on the localized orbital, and the number of excited phonons. We ﬁnd long-lived oscillations of the phonon that persist long after the end of the pulse. We compare the results for the interacting system to the noninteracting resonant level model.


A. Motivation
Ultrafast experiments constitute a new approach to exploring quantum many-body systems and provide a platform for developing new types of solid state devices for nanotechnology and quantum information processing (for review see Refs.[1][2][3][4]).One of the promising techniques is a recently developed terahertz STM (THz-STM) that integrates femtosecond lasers with scanning tunneling microscopes (STM) [5][6][7][8].This technique allows to combine atomic spatial resolution of STM with subpicosecond coherent temporal control of electron currents.Such experiments pose a new challenge to manybody theory to develop methods for analyzing the far out of equilibrium quantum dynamics of interacting manybody systems.Motivated by recent experiments in this paper we provide a theoretical analysis of THz-STM experiments of tunneling through a single localized orbital, such as a HOMO orbital in a pentacene molecule used by Cocker et al. [7] (see Fig. 1).Our analysis extends earlier theoretical studies of such systems (see [9,10] and references therein) by including a non-perturbative treatment of electron-phonon and electron-electron interactions.
In writing this paper we set ourselves two goals.Firstly, we present a new theoretical approach for analyzing non-equilibrium dynamics of the Anderson-Holstein impurity model.This method is based on the variational non-Gaussian wavefunctions introduced in Ref. [11] and further extended in the context of quantum impurity models in Refs.[12,13].This technique is versatile and can be applied to a broad class of nonequilibrium problems, including quenches, pump and probe experiments, and the analysis of AC and DC transport.Most of the previously developed approaches for the Anderson impurity problem focus either on high energy degrees of freedom on the scale of the electron-electron repulsion U [14], or the low energy sector with the scale set by the Kondo temperature T K [15][16][17][18].The advantage of our approach is that within the same framework it describes both high and low energy degrees of freedom without requiring numerical resources of NRG [19][20][21] or DMRG [22][23][24].To establish the validity of the new method we verify that it correctly reproduces basic results of the canonical Anderson impurity model: it captures the Kondo resonance in the spectral function at the Fermi energy and gives a conductance G = 2e 2 /h in the linear response regime.
We extend the analysis of the Anderson model to include the electron-phonon interaction on the impurity site and demonstrate that coupling to the phonons can strongly suppress the Kondo peak.Our second goal in this paper is to analyze a specific type of nonequilibrium dynamics: THz STM experiments through a single molecule.
For such experiments we compute the time dependence of the key observables: current through the system, number of electrons on the molecule, and number of excited phonons.To illustrate the important role of interactions in the system, we also discuss a non-interacting resonant level model (RLM)in the infinite bandwidth approximation.We show that, in this case, the transient current dynamics in the THz-STM experiemnts can be computed analytically.We discuss the difference in the results between the non-interacting RLM system and the Anderson-Holstein model.

B. Anderson-Holstein Model
The system we consider is shown schematically in Figs. 1 and 2. In the absence of laser light the chemical potential of the organic molecule is in the middle of the gap between the HOMO and LUMO orbitals of the molecule and there is no current through the system [7].The light pulse changes the relative energies of the molecular level and the reservoirs (see Fig. 2b) and allows for an ultrafast electron burst as a co-tunneling process between the STM tip, the molecule, and the metal.Due to the asymmetry of the electric field in the pulse [5,6] it is suffi-cient to consider only one of the orbitals in the molecule, which we take to be the HOMO orbital.We model the interaction between HOMO electrons using the Anderson type local repulsion U and include interaction of electrons with vibrations of the molecule [7] using Holstein type coupling of the phonon displacement to the number of electrons in the localized orbital [25].Thus, molecular degrees of freedom are given by the Hamiltonian Here, d σ are the annihilation operators of electrons in the HOMO, We model the STM tip and the metallic electrode as one dimensional chains of non-interacting electrons because we do not expect appreciable dependence on the nature of electron reservoirs.Our analysis can be easily extended to other type of reservoirs.
Here, c j,σ,a denote annihilation operators of electrons with spin σ in the reservoir a = L, R at site j.Coupling between the molecule and reservoirs is included via the hybridization term Different chemical potentials of the reservoirs account for a finite time-dependent bias voltage, where we assume that all the time dependent potential is applied between the molecule and the right reservoir (our analysis can be easily generalized to the arbitrary time dependent potentials for both reservoirs) The main effect of the femtosecond laser pulse is to change the bias voltage as shown in Fig. 2b.The full Hamiltonian of the system is then given by

C. Review of the theoretical formalism
The theoretical modeling of light assisted tunneling of electrons through a single molecule requires the analysis of nonequilibrium dynamics of the Anderson impurity model with an added complication, that of the electronphonon coupling.The Anderson impurity problem is one of the most fundamental models in the field of interacting electron systems.It provides the foundation for our understanding of a broad range of physical phenomena, including the Kondo effect in metals [25,26], the electron transport in mesoscopic structures [27][28][29][30], the formation of heavy fermion electron systems [31][32][33][34].A broad range of analytical tools has been developed to study equilibrium properties of this model including a Bethe ansatz solution [35,36], the renormalization group approach [37], the slave-particle method [38,39], and dynamical mean-field theory (DMFT) [40][41][42][43][44].However, the non-equilibrium dynamics of the model remains poorly understood.While most of the theoretical work on quantum dynamics of the Anderson model relied on the non-crossing approximation [16,45], other promising new approaches utilized real-time DMFT [46] combining the bold-line quantum Monte Carlo technique with a memory function formalism [47].Calculations using either of these approaches are demanding in terms of numerical resources.
In this paper we propose a new approach to analyzing the Anderson-Holstein impurity model both in and out of equilibrium.This approach is based on combining a unitary transformation of the generalized-polaron type with a family of non-Gaussian variational states for quantum impurity problems.The electron-phonon part of the transformation entangles the electron and phonon degrees of freedom and allows us to use factorized wavefunctions, in which bosonic and electronic parts are described using a Gaussian ansatz and a family of non-Gaussian variational wavefunctions, respectively.The role of the unitary transformation in the non-Gaussian ansatz for the Anderson model is to utilize exact conservation laws to reduce the number of impurity degrees of freedom at the cost of introducing additional interactions and correlations between the bath particles.The key difference between our method and the traditional approaches based on the polaron transformations [25] is that we allow parameters of the transformation to vary in time.The general philosophy of this method has been introduced earlier in Ref. [11], which demonstrated that this technique successfully describes equilibrium and dynamical properties of many important solid state systems including polarons in the Su-Schrieffer-Heegger and Holstein models, spin-bath models, and superconductivity in the Holstein model.Recently, Ashida et al. [12,13] demonstrated that this method efficiently captures the complicated dynamics of the anisotropic Kondo model, including the finite time crossover between the ferromagnetic and antoferromagnetic couplings [48].Since this is the first time that this method is used to study the Anderson model, we include a separate analysis of the equilibrium properties including the electron spectral function.We demonstrate that our approach accurately captures the formation of the Kondo resonance at the Fermi energy with the width set by the Kondo energy scale.This demonstrates that our method is versatile and captures both short time/high energy and long time/low en-ergy aspects of the Anderson model.

D. Organization of the paper
This paper is organized as follows: In Section II we introduce the general formalism of non-Gaussian variational wavefunctions.In Section III we use the imaginary time flow approach to analyze the ground state properties of the Anderson and Anderson-Holstein models focusing on the regimes of single and double occupancy of the localized orbital.We compute electron spectral functions and demonstrate that in the Kondo regime of the Anderson model our approach correctly reproduces a resonance at the Fermi level.We show how the spectral functions get modified upon adding the interaction with phonons.Two prominent features are the suppression of the Kondo resonance and the appearance of the phonon shake-off peaks.As additional check of the validity of our method we compute spin correlation functions between electrons on the localized orbital and in the electron reservoirs.We demonstrate that in the Kondo regime we get the expected oscillating correlations decaying as a power law of the distance to the localized orbital.In Section IV we introduce a real time analysis of the Anderson and Anderson-Holstein models.For the Anderson model we analyze the DC transport and compute the full non-linear conductance.We demonstrate that in the Kondo regime our method correctly gives a linear conductance equal to 2e 2 /h.For the Anderson-Holstein model, we focus on the analysis of photocurrent induced by the THz pulze in the single and double occupancy regimes.We present results for the time evolution of the current, the occupation number of the localized orbital, and the phonon displacement.We show the dependence of the total transferred charge on the amplitude of the electromagnetic pulse.Section V contains a summary of our results and a discussion of interesting open problems.Many technical details of the calculations are not presented in the main text of the paper but delegated to Appendices.This includes a derivation of the analytical results for photocurrent in the non-interacting resonant level model.

II. NON-GAUSSIAN VARIATIONAL ANSATZ
In this section, we introduce the variational ansatz which we will use to study both the ground state and non-equilibrium real-time dynamics.The unitary transformation U ph entangles the phonon mode with the electronic degrees of freedom in HOMO and in the reservoirs.The unitary transformation U A uses parity conservation to partially decouple the impurity degrees of freedom.
Throughout the paper we use the terminology of the Anderson impurity model and refer to the electronic degrees of freedom in the molecule as impurity degrees of freedom.Finally |Ψ GS f,b are Gaussian states for electrons and the phonon mode.

A. Electron-phonon polaron transformation
The generalized polaron transformation U ph = e S is defined by the generating operator with the time-dependent variational parameters λ = (λ x , λ p ) T .Here R = (x, p) T is the quadrature of the phonon mode.When the polaron transformation ( 7) is applied to the system Hamiltonians we find that the lead Hamiltonians H L,R do not change, whereas the molecular part of the Hamiltonian H mol = H A + H phon becomes and the hybridization term changes to The polaron transformation reduces the electronphonon interaction to G λ = (g − ω b λ p , ω b λ x ) T , renormalizes the single particle energy level εd = ε d − 3(ω b λ T λ − 2gλ p ) and the on-site interaction Ũ = U + 2(ω b λ T λ − 2gλ p ), whereas the hybridization term V e −iR T λ between the phonon-dressed electrons in HOMO and reservoirs aquires polaronic dressing.

B. Parity operator and impurity decoupling transformation
Refs [11,12] introduced an efficient way of solving quantum impurity problems by utilizing parity conservation.In particular, in the Kondo impurity problem one can construct an exact unitary transformation that maps the conserved parity operator to one of the components of the impurity spin.After the transformation, the impurity spin is decoupled from the reservoir at the cost of introducing interactions between bath fermions, which, however, can be effectively handled using the Gaussian part of the wavefunction.In this paper we generalize such transformation to the Anderson impurity problem.
In the Anderson model one can not completely decouple the impurity from the reservoir, but it is possible to reduce the number of degrees of freedom so that the Gaussian ansatz for the wavefunction can be utilized efficiently.
To simplify notations we introduce a four component representation of the electronic Hilbert space on the molecule: We define four component matrices in this space so that the first Pauli matrix acts in the space (1,2) ↔ (3,4).For example, When the Pauli matrix is in the second position in the tensor product, it acts simultaneously in sub-spaces 1 ↔ 2, 3 ↔ 4, so that We define the parity operator for spin ↑ electrons on the molecule The choice of the overall sign is a matter of convenience.In the tensor notations introduced earlier Note that the operator Σ z performs a π rotation in the subspace (| ↑ , | ↓ ) as well as a π rotation in the subspace (|0 , | ↑↓ ).The parity operator for spin ↑ fermions in the bath is defined as The parity operator for the system as a whole is conserved since the Hamiltonian preserves the number of electrons with a given spin.Mathematically this means that [P, H 1 ] = 0.
In the spirit of Refs [12,13] we introduce a transformation that maps the conserved operator P entirely into the impurity degrees of freedom.Following this transformation, the original conservation low for P becomes a conservation law for the imourity electrons.We consider a unitary operator with Σ y = −σ x ⊗ σ y Direct calculation shows that We observe that the operator X does not contain any degrees of freedom of the reservoir electrons.Another important feature of the operator X is that it only connects states with the same electron parity, i.e. it does not mix between subspaces (| ↑ , | ↓ ) and (|0 , | ↑↓ ).Hence the operator X is bosonic and its eigenstates are physically well defined states.After performing the transformation U A on the Hamiltonian we are guaranteed that the operator X commutes with the transformed Hamiltonian H 2 , which means that we reduced the number of degrees of freedom corresponding to the molecule.Let us discuss the implications of this fact in more detail.Operator X has eigenvalues ±1.The eigenstates corresponding to the eigenvalue +1 make a two dimensional Hilbert space with basis vectors The eigenstates corresponding to the eigenvalue −1 make an orthogonal two dimensional Hilbert space with basis vectors Conservation of X means that the dynamics described by the Hamiltonian H 2 can not mix between subspaces (17) and (18).Hence for the dynamics in the +1 subspace we need to consider only two electronic states in the molecule: |+ s and |+ c .We introduce the fermion operator where F = j,σ,a c † j,σ,a c j,σ,a is the total number operator of electrons in both reservoirs.In the +1 subspace Hamiltonian H 2 can be expressed entirely in terms of the reservoir operators c i,σ,a , for the molecule f ± + , and phonon operators.Analogously, in the −1 subspace (18) we introduce the fermion operator that connects two states with different electron parity The explicit expression for H 2 can be obtained from a straightforward but somewhat lengthy calculation.The part of the Hamiltonian corresponding to the left and right reservoirs does not change.The part of H 2 corresponding to the molecule, H mol = H A + H phon , becomes The transformed hybridization term is given by

C. Gaussian part of the wavefunction
A convenient way of defining Gaussian wavefunctions in equation ( 6) is to consider a unitary Gaussian transformation acting on the electron and phonon vacuum GS for bosons and fermions respectively, where δR = R − ∆ R is the quadrature fluctuation around its mean value, and Majorana operators are given by linear combinations of the form [11,49].
In the next section we will obtain the variational ground state of the Anderson-Holstein model by analyzing the flow of variational parameters λ, ∆ R , and Γ b,f (m) in imaginary time.We will then derive the real time equations of motion (EOM) which we will apply to study the dynamics.

III. GROUND STATE PROPERTIES
In this section, we use imaginary time evolution to approximate the ground state of the Anderson-Holstein model.We will demonstrate that in the case of the canonical Anderson model our method correctly captures the non-perturbative Kondo effect, including the formation of a resonance at the Fermi energy and presence of characteristic spin correlations between the impurity and reservoir electrons indicating the presence of Kondo spin screening clouds.We project the equations of motion (EOM) onto the tangential plane of the variational manifold (6), and obtain the EOM for the variational parameters λ, ∆ R , and Γ b,f (m) (see Ref. [11] for details).As τ → ∞ the system reaches a fixed point which approximates to the ground state of the system.The steady state solution of the EOM in the limit τ → ∞ determines the ground state properties, including the occupation number, the magnetization, and the spectral function of electrons in HOMO, as well as correlation functions between HOMO and reservoirs.

A. Equations of motion in imaginary time
We first analyze the structure of the tangential vectors of the variational manifold (6).To this purpose, we introduce the unitary operator that generates the Gaussian state and transforms R and f /2 are related to the symplectic and unitary transformations via S = e iσ y ξ b and U f = e iξ f .Let us consider the tangential vector for the variational wavefunction (6).It is defined as the derivative of |Ψ NGS with respect to τ We separated the V 1 , V 2 , V h terms in equation ( 24) based on the number of creation and annihilation operators.The V 1 term in ( 24) is linear in phonon operators It is determined by the expectation value P z f † γ f γ GS in the Gaussian state.
The V 2 term in (24) is quadratic in phonon and electron creation/annihilation operators Note that since the operator V 2 acts on the vacuum state we should bring it to normal ordered form.This is indicated by "::" in equation (26).Due to the presence of matrices performing rotations of creation/annihilation operators, S for phonons and U f for fermions, this normal ordering is with respect to the instantaneous Gaussian state.The matrix O f is defined from where the diagonal matrix I f has only one non-zero element (I f ) 11 = 1, and The V h term in equation ( 24) contains higher order terms defined by the expansion Equation ( 30) should be understood as the definition of P h .The projection on the tangential vector where and the effective hybridization strength V eff = V e λ is reduced by the factor e λ ≡ e iR T λ By projecting on the tangential vector The explicit form of the Gram matrix G and the vector ξ τ is given in Appendix A. By solving Eqs. ( 31), ( 33), ( 34) and ( 35) numerically, we obtain the ground state configuration in the limit τ → ∞.We note that the EOM ( 35) is quadratic in ∂ τ λ, but it can be reduced to the linear ordinary differential equation (ODE) (see Appendix B).

B. Physical observables
In this subsection we will discuss several physical quantities that can be used to characterize the ground state of the Anderson-Holstein model (here the bias voltage V e = 0).
In the sector γ, the ground state energy

GS
, the molecular energy E mol = E A + E phon : and the hybridization energy E V =Re I dc GS , where h a=L,R = I 2 ⊗ (−t 0 δ i,j±1 − µ a δ ij ).As pointed out in Ref [11] the energy E γ monotonically decreases during the imaginary time evolution.
The occupation number and the magnetization of electrons in HOMO are given by Correlation between electrons in the HOMO and reservoirs are characterized by the correlation functions C α=x,y,z = , where τ α is the Pauli matrix.In the transformed frame, the correlation functions can be expressed as expectation values computed in the Gaussian state The hole excitations in HOMO interact with the charge in the substrate via the Coulomb interaction, which results in the displacement of phonon.The spectral function A(ω) = −ImG R (ω)/π characterizes the properties of excitations above the ground state.More specifically, for the Anderson-Holstein system, it is determined by the Fourier transform G R (ω) = dte iωt G R (t) of the retarded Green function Following the unitary transformation given by U ph U A , the Green's function becomes where the evolution d↓ (t) = e iH2t d↓ e −iH2t of the fermionic operator d↓ = e −iR T λ F is governed by the Hamiltonian H 2 , and We The Green function G R (t) contains the average values on the Gaussian state, which can be obtained analytically in terms of the covariance matrices and ∆ R (see Appendix C).For instance, where α = λ T Γ b λ and ω re is the simplectic eigenvalue of Ω re .Eventually, the imaginary part of the Fourier transform gives the spectral function where we add a small imaginary part iδ in the numerical calculation of Fourier transforms G >(<) .
In the next two subsections, we show the occupation number, the correlation functions, the displacement, and the spectral function for the Anderson-Holstein model.

C. Anderson model in equilibrium
In this subsection, we present results for the ground state properties of the Anderson model, i.e. when the electron-phonon coupling g = 0. We expect that in the Kondo regime ε d < 0, ε d + U > 0, and Γ = V 2 < (U, ε d ), HOMO is singly occupied, i.e., n d ∼ 1 and its magnetization m z = 0, which reflects screening of the impurity spin by electrons in reservoirs.An important signature of the Kondo regime of the Anderson model is the formation of anti-ferromagnetic correlations between spins of electrons in HOMO and on the adjacent sites in the reservoirs.These correlations will be absent when the impurity is in the doubly occupied regime with and n d ∼ 2, which is relevant to the HOMO in the THz STM experiments.In Fig. 3, we plot the correlation functions C α=x,y,z of the electrons in HOMO and one of the reservoirs in the Kondo regime (U, ε d , Γ) = (1, −0.5, 0.16) and (1, −0.5, 0.04), and the double occupation regime (0.25, −0.5, 0.04) and (0.05, −0.5, 0.04).In the Kondo regime (Figs.3a and 3b), the SU (2)-symmetric C α shows significant anti-ferromagnetic correlations of electrons in the HOMO and the reservoirs.In the double occupation regime (Figs.3c and 3d) we observe small values of the correlations C α ∼ 10 −3 , which indicates weak correlations between the singlet electron pair in the HOMO and the reservoirs.In Fig. 4, we show the spectral function A(ω) in the Kondo and the doubly occupied regimes, where we take δ = 0.01.In Fig. 4a, the spectral function in the Kondo regime (U, ε d , Γ) = (1, −0.5, 0.04) displays both the Kondo resonance peak around ω = 0 and the resonance with energy levels ε d and ε d + U .In the the double occupation regime (U, ε d , Γ) = (0.05, −0.5, 0.04), as shown in Fig. 4b, the Kondo resonance peak vanishes.

D. Anderson-Holstein model
In this subsection, we present results for the equilibrium properties of the full Anderson-Holstein model with finite electron-phonon interaction.Consequencies of the electron-phonon interaction g include a change of the effective single particle energy level, softening of the onsite repulsive interaction, and polaronic dressing of hybridization.Note, for example, that when the occupation number on the molecular orbital n d is different from two, we get a finite displacement of the phonon operator x 0 = x , which favors partial occupation of the HOMO.We observe that in the Kondo regime including the electron-phonon intraction tends to suppress the formation of the singlet cloud, while in the doubly occupied regime it favors the creation of the hole excitations in HOMO.
In Fig. 5a, we show the occupation number n d and the displacement x 0 as a function of the coupling strength g in the Kondo regime (U, ε d , Γ) = (1, −0.5, 0.16), where the magnetization m z = 0.In Fig. 5b, we show the correlation function C α for different g = 0.1, 0.5, and 0.9 in the Kondo regime.As g increases, the displacement x 0 becomes larger, which lifts the single particle energy level of HOMO, therefore the the occupation number n d decreases from 1 to 0. The finite electron-phonon interaction also softens the on-site interaction and the effective hybridization strength V eff < V , thus, as shown in Fig. 5b, the larger the coupling strength g, the weaker the correlation C α between HOMO and leads.Eventually, the Kondo singlet is destroyed by the finite electron-phonon coupling.
In Fig. 5c, we show the occupation number n d and the displacement x 0 as a function of g in the double occupation regime (U, ε d , Γ) = (0.05, −0.5, 0.04), where the magnetization m z = 0.In Fig. 5d, we show the correlation functions C α for different g = 0.1, 0.5, and 0.9 in the double occupation regime.Since the single particle energy level is lifted, the particle number n d decreases from 2 to 0 as g increases.For small electron-phonon couplings, e.g., g = 0.1, HOMO is mostly doubly occupied, and C α displays a weak correlation of HOMO and reservoirs.In the intermediate coupling regime, e.g., g = 0.5, the occupation number n d decreases to 1.5, meaning that electrons in HOMO are in the superposition of doubly and singly occupied states.Since the total spin in HOMO is totally screened, i.e., m z = 0, the singly occupied electron forms the singlet state with lead electrons.This explains somewhat the presence of a stronger correlation C α in the intermediate regime than that in the weakly coupling regime.When g increases further, e.g. to g = 0.9, the occupation number n d is reduced to 0, and the correlation C α becomes vanishingly small.
In Fig. 6, we show the spectral function A(ω) in the Kondo and doubly occupied regimes.In Fig. 6a, the spectral function A(ω) in the Kondo regime (U, ε d , Γ) = (1, −0.5, 0.04) is shown, where the coupling strength g = 0.2 and 0.4.For g = 0.2, the renormalized single particle energy level is slightly lifted to εd ∼ −0.40 and the onsite interaction Ũ ∼ 0.93 is softened due to the electronphonon interaction.Since εd < 0 and εd + Ũ > 0 still in the Kondo regime, the Kondo peak survives in the spectral function.However, the system deviates from the symmetric Anderson model, i.e., |ε d | = εd + Ũ , therefore, in the spectral function, the two peaks around εd and εd + Ũ become asymmetric.The small α ∼ 0.02 shows that the electron excitation in HOMO is weakly dressed by the phonon, and the spectral function is dominated by the pure electron excitation with the phonon in the vacuum state, i.e., n = 0 in Eq. (45).For the larger g = 0.4, εd is shifted to −0.06 and the on-site interaction is reduced to Ũ ∼ 0.71.The spectral function shows that the Kondo peak is destroyed, and two peaks around εd and εd + Ũ still survive.Since the electron excitation in HOMO is dressed by the phonon with larger average phonon number α ∼ 0.08, two peaks corresponding to the single phonon excitation appear in the spectral function, as shown in the (black) box of Fig. 6a.
In Fig. 6b we show the spectral function A(ω) in the doubly occupied regime (U, ε d , Γ) = (0.05, −0.5, 0.04) for two values of the electron-phonon coupling strength g = 0.2 and 0.4.We point out that for g = 0.2 the renormalized single particle energy level is shifted to εd ∼ −0.41 and the on-site electron-electron interaction becomes effectively attractive Ũ ∼ −0.01 due to the phonon mediated attraction.The small α ∼ 0.01 implies that the electron in HOMO is weakly dressed by the phonon mode, as a result, the contribution of the phonon excitation to the spectral function is negligible.Comparing with the spectral function without electronphonon interaction, we find that the peak is shifted to the larger frequency around the lifted single particle level εd .When the coupling constant is further increasing, e.g., g = 0.4, εd is shifted to −0.1 and the attractive interaction Ũ ∼ −0.22 becomes stronger.The electron excitation in HOMO is dressed by phonons with larger average number α ∼ 0.06, which results in the visible peak corresponding to the single phonon exitation, as shown in the (black) box of Fig. 6b.

IV. REAL TIME DYNAMICS
In this section, we apply the non-Gaussian ansatz (6) to study real-time dynamics for the Anderson-Holstein system.Similar to the procedure used in the imaginary time evolution, we project the Schrödinger equation to the tangential space of the variational manifold ( 6) and obtain EOM for variational parameters λ, ∆ R , and Γ b,f (m) .Analysis of these differential equations allows us to explore a large variety of non-equilibrium phenomena We now present results for transport through the molecule in the two cases: with and without electronphonon interaction.To find DC conductance we perform a quench-type protocol, when we start with the reservoirs at different chemical potentials but disconnected from the molecule, and hence from each other.We switch on the molecule-reservoirs coupling and then evolve the system in real time until it reaches the steady state but before electron wavepackets reflected from the outer ends arrive back at the molecule.The steady state current I as a function of (finite) V e gives us the full nonlinear current voltage characteristic of the system.Differential conductance is defined through the relation σ c = ∂ Ve I.In the case of ultrafast THz STM experiments, bias voltage is applied as a time dependent pulse given in equation ( 54).This voltage pulse results in a transfer of a finite number of electrons between the two reservoirs.We compute the full time dependent evolution of the current in the system.We demonstrate that in the case of finite electron-phonon interaction a burst of current gives rise to vibrations of the molecules which persist well after the duration of the pulse.This was observed in experiments by Cocker et al. [7].

A. Equations of motion
Projection of Eq. ( 46) to the tangential vectors gives EOM for variational parameters in the real time evolution.When we perform projection on the vectors and To obtain the time evolution of λ, we perform the projection on the tangential vector U ph U A U GS V h |0 and obtain where ξ t = iξ τ .While Eq. ( 50) is non-linear in λ, it can be reduced to a linear ODE.Details are presented in Appendix B. Thus, the description of a broad range of nonequilibrium phenomena in the Anderson-Holstein model can be reduced to solving Eqs. ( 47)-( 50) for the time evolution of variational parameters.Before presenting results of our analysis we comment on one important technical aspect of the calculations.The total electron number operators should be conserved throughout the real time evolution.This fundamental conservation law should not be affected by the transformations of the Hamiltonian U ph and U A .However, when computing long time evolution needed for finding steady states the numerical errors may accumulate and lead to the violation of particle number conservation.To circumvent this numerical problem, we introduce a penalty term in the Hamiltonian, where Λ is chosen to be much larger than all the energy scales in the system, and N↑(↓) is the average number of spin-up (down) electrons.Note that specific value of Λ turns out to be unimportant for all results that we discuss in this paper.The mean-field Hamiltonian H Λ,m in the Majorana basis for the penalty term H Λ is derived in Appendix D, which modifies the mean-field Hamiltonian (A1) as H m → Hm = H m + H Λ,m .In Eq. ( 49), the substitution of H m by Hm leads to the conserved particle numbers N ↑(↓) = N 0 ↑(↓) in the real time evolution.
For the system with bias time (in)dependent V e (t), we can study the occupation number n d , the displacement R , and the correlation function C α .In the non-equilibrium state, the current is defined as I = j,σ ∂ t c † j,σ,L c j,σ,L .The Heisenberg equation of motion leads to where the average values depend on the covariance matrix Γ m (see Appendix A).The derivatives σ c = ∂ Ve I give the conductance.

B. Transport in the Anderson model
In this subsection, we present results for electron transport in the Anderson model without the electron-phonon interaction.We compute the non-linear conductance in both the Kondo and doubly occupied regimes.
For the Kondo regime we set parameters (U, ε d , Γ) = (1, −0.5, 0.16) and for the doubly occupied regime we choose (U, ε d , Γ) = (0.05, −0.5, 0.04), Fig. 7 shows the occupation number n d and the conductance σ c in the steady state as the function of V e .Fig. 7a shows the conductance σ c around zero bias in the Kondo regime.The peak value σ c = 2 originates from the Kondo resonance and agrees with the Friedel sum rule.In the doubly occupied regime, when the bias V e crosses ∼ ε d + U , the energy level of the doubly occupied state in HOMO is higher than the Fermi surface of the right reservoir.As a result, the transport channel is turned on and the electron in HOMO tunnels to the reservoir, resulting in the decay of n d and the appearance of a conductance peak around ε d + U in Fig. 7b.
In the experiment, an ultrafast pulse is applied to shift the chemical potential of the right lead, where the HOMO is doubly occupied.The effect of the ultrafast pulse centered at the instant t c is described by the time-dependent bias where V e,0 is the intensity, α determines the width of the pulse, and ω d is the frequency.
For different pulse intensities V e,0 , the transient current I and the occupation number n d as a function of time t are shown in the left and right panels of Fig. 8.The two rows in Fig. 8 correspond to system parameters (U, ε d ) = (0.05, −0.5) and (1, −0.5) in the double occupation and Kondo regimes, where Γ = 0.04, t c = 5, α = ω d = 1, and t 0 is taken as the unit.The number N tran of electrons transferred from the left lead to the right lead as the function of pulse intensity V e,0 is displayed in Fig. 9.
The pulse has a single period Sine-type shape, which first shifts the Fermi level of the right reservoir downwardly and then lifts it above the Fermi level of the left reservoir after t c .In the double occupation regimes, as shown by the first row of Fig. 8, when the Fermi energy is resonant with ε d + U < 0 at the instant t res < t c , the transport channel is turned on and the transient current flowing to the right reservoir establishes.However, after the instant t c when the Fermi level of the right reservoir is higher than that of the left reservoir, there is no resonant energy level.Eventually, due to the asymmetric spectral structure around the Fermi level ε F = 0, the light pulse induces the non-zero net charge transported ∼ 10 −2 from the left to the right, as shown in Fig. 9a.In the Kondo regime, the energy spectrum is symmetric around the Fermi level, as a result, the net charge transported is highly reduced.We note that a finite value of the net charge ∼ 10 −3 in Fig. 9b most likely arises from the nonlinear electronic dispersion relation in reservoirs.

C. Transport in the Anderson-Holstein model
In this subsection, we concentrate on the phonon excitations generated by the ultrafast pulse in the doubly occupied regime.Since the light pulse induces a transient current, i.e., the transport of electrons between HOMO and leads, the holes are created in the HOMO.The Coulomb interaction between the charge in the substrate and the hole excitation in HOMO induces the molecular vibration after interacting with the ultrafast pulse, which is described by the time-dependent average value R = (x 0 , p 0 ) T of the quadrature.In the first and second rows of Fig. 10, we show the occupation number n d , the average value R = (x 0 , p 0 ) T of the quadrature, and the transient current I in the double occupation and Kondo regimes (U, ε d ) = (0.05, −0.5) and (1, −0.5), respectively.Here, Γ = 0.04, g = 0.2, t c = 5, α = ω d = 1, and t 0 is taken as the unit.The electrons are transferred from the left reservoir to the right one when the terahertz pulse is applied, and the long-lived oscillation of net charge transport around a center value is observed, as shown in the left panel of Fig. 11 for the double occupation and Kondo regimes.Here, the center value of the oscillation as a function of pulse intensity is also shown in the right panel of Fig. 11.
In the initial stage, i.e., t < t c , the ultrafast pulse shifts the Fermi level of the right reservoir downwardly, thus electrons in HOMO flow to the right reservoir, and the occupation number n d decreases, as shown in the first column of Fig. 10.Since more holes are generated in the HOMO, the molecule is driven away from the equilibrium position along the negative direction by the electron-phonon interaction, as shown in the second column of Fig. 10.When the Fermi level of the right reservoir is resonant with εd , the transport channel is turned on, and the significant transient current is formed, as shown in the third column of Fig. 10 and the left panel of Fig. 11.
In the intermediate stage, i.e., t c < t < t c +π/ω d , since the Fermi level of the right reservoir is above that of the left one, the HOMO is repumped, as shown in the first column of Fig. 10.However, the resonant transport is absent since the Fermi level of the right reservoir is detuned from LUMO.Eventually, less electrons flow back to the left reservoir, and the finite net charge is transported to the right reservoir, as shown in the Fig. 11.In the final stage t > t c + π/ω d , when the light pulse is turned off, the long-lived phonon mode with frequency ∼ ω b survives, as shown in Fig. 10b.Due to the indirect coupling to the reservoir through electrons in HOMO, the phonon mode has a finite life-time, which is revealed by the slowly decaying amplitude of the oscillation in 10b.The molecular vibration induces the long-lived oscillation of current and net charge transport around their values in the steady state, as shown in the third column of Fig. 10 and the left panel of Fig. 11.

D. Discussion of results
Before concluding this section we would like to highlight several results of our analysis.
The main difference between the regimes of single and double occupancy is the dependence of the net transferred charge on the amplitude of the THz light pulse V e,0 .In the former case we observe a linear dependence of the transferred charge on V e,0 (see Fig. 8b), whereas in the latter case we find a non-linear dependence (see Fig. 8a).A special feature of the single occupancy regime is the existence of the Kondo resonance at the Fermi energy, which gives rise to finite DC conductance.This however is not the whole story.In Appendix E we present the analysis of the photoinduced current in the RLM with the energy of the localized state set exactly at the Fermi energy, d = 0, which would also allow for finite DC conductance.In the RLM we also find a non-linear dependence of the photocurrent on V e,0 (see fig. 12d).We remind the readers that that the integral of the time dependent electric field in the THz pulse is zero.Hence for a completely linear system the net transferred charge should vanish.The linear dependence of photocurrent on V e,0 in the Kondo regime is thus a surprising feature of the system.
Another interesting result of our analysis is the existence of different time scales characterizing the transient dynamics.The occupation number of electrons on the localized orbital exhibits two distinct scales in its dynamics.There is a relatively short timescale over which the strong amplitude deviation from equilibrium configuration decays.It is followed by small amplitude oscillations that match the frequency of the photo-excited phonon mode, with both oscillations of the phonon amplitude and the electron occupation number exhibiting slow decay.The existence of different time scales is a common feature of pump and probe experiments in strongly correlated elec-tron systems (see e.g.Ref. [50]).Slow relaxation of phonon excitations was one of the key features observed in experiments by Cocker et al. [7].Our analysis suggests that the origin of this relaxation is due to phonon displacements modifying virtual tunneling processes of reservoir electrons into the localized orbital, which ultimately allows the phonon energy to be converted into particle-hole excitations.

V. SUMMARY AND OUTLOOK
We introduced a new method for analyzing the Anderson-Holstein impurity model.Key ingredients of the approach are two unitary transformations that use the parity conservation to partially decouple the impurity degrees of freedom and a generalized polaron transformation to entangle phonons and electrons.An appealing aspect of the new method is that degrees of freedom at very different energy scales, from the local repulsion U to the Kondo temperature T K , are described within the same framework and without the numerical demands of the NRG and DMRG calcualtions.To verify the accuracy of the new approach we computed the properties of the Anderson model, including the equilibrium spectral function, linear and non-linear DC conductance.We demonstrated that we correctly reproduce known results for this model.We extended the analysis of the Anderson model to include electron-phonon interactions and showed that it leads to a suppression of the Kondo resonance.We used our approach to analyze THz-STM experiments of tunneling through a single molecule.We found that a picosecond light pulse that induces the current flow gives rise to strong oscillations of the molecular phonon, which persist long after the end of the pulse.We analyzed the time dependence of the current induced by the THz pulse and showed a strong difference between the Kondo and doubly occupied regimes of the molecule.
Our work can be extended in several directions.One interesting question is developing a deeper understanding of the interplay of electron-electron and electron-phonon interactions in pump and probe experiments [51][52][53][54][55]. Recent time-resolved ARPES/Xray experiments by Gerber et al. measured changes of the electron energy relative to the phonon displacement in a photoexcited iron selenide [56].They observed a strong renormalization of the electron-phonon coupling relative to the value predicted from band structure calculations and attributed it to the effects of electron-electron interactions.Extension of the formalism presented in this paper can be used to study the time dependent spectral function of electrons following a THz light pulse.This analysis will provide a direct comparison between changes of the electron energy and the phonon amplitude.
Another interesting question is the analysis of shot to shot fluctuations in the number of electrons transferred through the junction during the THz pulse.In the case of DC transport the study of shot noise has been a valu-able tool for understanding underlying many-body states (see Refs. [57][58][59] for a review), including the demonstration of fractional statistics in the fractional quantum Hall states [60,61] and understanding the back-scattering mechanism in the Kondo system [62][63][64] as a probe of the low temperature fixed point.
The formalism developed in our paper introduces a powerful new technique that can be used for the theoretical analysis of many different types of nonequilibrium problems.For example, it should provide useful insights into the analysis of resonant XRay scattering experiments [65] in materials with mobile electrons.The essence of these experiments is that an incident photon excites an electron from the core orbital into one of the unoccupied bands, which results in a sudden introduction of a positively charged hole in the sea of conduction band electrons [66].The attractive potential of the hole is strong enough to allow the formation of the excitonic bound state, however, such bound state can be occupied by only one of the electrons.Hence, an analysis of the many-body dynamics following absorption of the incident photon should include the local repulsion on the hole site.The formalism developed in this paper can be used for computing RIXS processes involving the excitation of bosonic modes, such as phonons and magnons, as well as a continuum of electron-hole pairs.Another interesting direction for applying the ideas presented in this paper is to use variational non-Gaussian states as a solver for non-equilibrium DMFT calculations [46,67]. in the Majorana basis, where the matrices of the tangential vector U ph U A U GS V h |0 is determined by The vector is the overlap of the tangential vector and the right hide side of Eq. ( 22), where

FIG. 2 :
FIG. 2: (a) Anderson-Holstein microscopic model for tunneling through a molecule.Only the HOMO orbital on the molecule is considered.Coulomb repulsion between electrons on the molecular orbital is included through the Anderson model.Interaction of electrons with the optical vibrational mode is included using the Holstein model.(b) Light assisted tunneling through a molecule in THz-STM experiments.Electric field from the femtosecond pulse modifies energy differences between the reservoirs and the molecular level.
the total number of electrons in the molecule n σ=↑↓ = d † σ d σ , b is the annihilation operator for the phonon with frequency ω b , and ε d is the energy of the HOMO orbital.
for the covariance matrices Γ b,m , where Ω re = ω b I 2 − 2Re I dc GS λλ T and H m are the mean-field single particle Hamiltonians of phonon and electrons, and

FIG. 5 :
FIG. 5: The occupation number n d , the displacement x0, and the correlation function Cα in the Kondo regime (a)-(b) and the double occupation regime (c)-(d), where ε d = −0.5, ω b = 1, and the hopping amplitude t0 is taken as the unit.(a) The occupation number n d and the displacement x0 for U = 1 and Γ = 0.16; (b) The correlation function Cα for U = 1 and Γ = 0.16; (c) The occupation number n d and the displacement x0 for U = 0.05 and Γ = 0.04; (d) The correlation function Cα for U = 0.05 and Γ = 0.04.

FIG. 7 :
FIG. 7: The occupation number n d and the finite bias conductance σc in the Kondo regime (a) U = 1 and Γ = 0.16 and the double occupation regime (b) U = 0.05 and Γ = 0.04, where ε d = −0.5 and the hopping amplitude t0 is taken as the unit.
approximate d↓ (t) ∼ e iHMFt d↓ e −iHMFt by the mean-field Hamiltonian H MF = H p MF +H e MF with H p MF = δR T Ω re δR/4 and H e MF = iA T H m A/4, where Ω re and H m are determined by the average value of quadrature and covariance matrices in the ground state.Since the bosonic and electronic parts in the Gaussian state are factorized, the retarded Green function reads