A Nonequilibrium DMFT for the Holstein model

# Coexistence of excited polarons and metastable delocalized states in photo-induced metals

## Abstract

We study how polaronic states form as a function of time due to strong electron-phonon coupling, starting from a hot electron distribution which is representative of a photo-induced metallic state immediately after laser excitation. For this purpose we provide the exact solution of the single-electron Holstein model within nonequilibrium dynamical mean-field theory. In particular, this allows us to reveal key features of the transient metallic state in the numerically most challenging regime, the adiabatic regime, in which phonon frequencies are smaller than the electronic bandwidth: Initial coherent phonon oscillations are strongly damped, leaving the system in a mixture of excited polaron states and metastable delocalized states. We compute the time-resolved photoemission spectrum, which allows to disentangle two contributions. The existence of long-lived delocalized states suggest ways to externally control transient properties of photo-doped metals.

71.10.Fd

## I Introduction

Ultra-short laser pump-probe techniques in condensed-matter systems have opened the possibility to generate correlated nonequilibrium phases, such as photo-induced metallic states in Mott insulators (1), and to study their dynamics on femtosecond timescales. On a fundamental level, seeing how correlations evolve in time can shed new light on many-body effects which have been investigated for decades under equilibrium conditions. A paradigm example is the formation of polaronic quasiparticles, i.e., the self-trapping of an electron in a lattice distortion, or “phonon cloud”. This phenomenon was predicted in the early days of quantum mechanics (2) and has been thoroughly investigated for a large class of systems (3); (4); (5), more recently also for ultra-cold gases and trapped ions (6); (7). In nonequilibrium, however, many questions related to the dynamics of systems with strong electron-phonon coupling remain only partially understood.

Signatures of strong electron-phonon coupling and polaronic effects in photo-excited systems have been found for the self-trapping of excitons (8); (9); (10); (11), in Mott insulators (12); (13), and organic materials (14); (15); (16). A direct observation of the self-localization process was achieved by two-photon photoemission for electrons in surface states which couple to adsorbate layers (17); (18); (19); (20). While polaronic effects can be visible already within one phonon period after photo-excitation, it is not entirely clear how, and how fast, the actual polaron ground state can be reached. The presence of non-equilibrated polarons, on the other hand, would determine carrier mobilities in transient metallic states and can thus be of importance also for potential technological applications like ultra-fast switches. It is therefore important to pinpoint signatures of excited polarons, to understand their properties, and whether these can be controlled, e.g., by the photo-excitation process.

These questions have motivated considerable effort to understand the nonequilibrium polaron dynamics from a theoretical perspective. A large body of work has been performed for the Holstein model (21), which describes tight-binding electrons coupled to an optical phonon with frequency . The physical picture for the polaron formation process which has emerged from these studies suggests two important bottlenecks: For large , one finds long-lived beating between well-separated polaron sub-bands in the many-body spectrum (22); (23), while in the opposite and experimentally very relevant adiabatic regime, in which is smaller than the electron hopping, a semiclassical argument (24); (25) predicts an energy barrier between delocalized and localized states. In the present work we solve the model exactly in the large coordination limit to see how relaxation of high energy electrons by emission of phonons, strongly damped coherent oscillations, long-lived delocalized states, and trapping in excited polaron states come together in particular in the adiabatic limit and how they are reflected in characteristic signatures of the photoemission spectrum.

Even for a single electron (the relevant limit to describe diluted polarons), the Holstein model is difficult to be solved in nonequilibrium, because established approaches like Quantum Monte Carlo (26) cannot be used. Variants of exact diagonalization (22); (23); (27); (28); (29); (30); (31) provide an accurate and versatile description of the electron-phonon dynamics in many regimes, but they rely on a cutoff of the phonon Hilbert space and become challenging in the adiabatic regime in which the phonon cloud involves a large number of oscillator quanta. Our work is based on nonequilibrium dynamical mean-field theory (DMFT) (32), which is exact in the limit of large lattice coordination numbers (33). In DMFT, a lattice model is mapped onto a single impurity coupled to a self-consistent bath. While the real-time dynamics of this impurity problem can usually be solved only approximately (see, e.g., Ref. (34) for the Holstein model), the limit of low electron density in the Holstein model provides a remarkable exception. In equilibrium, the DMFT equations for this case can be written exactly in terms of a continued fraction for the electron Green’s function (35). Technically, this solution is similar to earlier diagrammatic approaches (36); (37), and also to the momentum averaged technique (38); (39), which have provided a solution throughout all regimes of the single-electron Holstein model in equilibrium (adiabatic and non-adiabatic, weak and strong coupling). These diagrammatic techniques rely on a momentum-independent self-energy which is approximate in finite dimensions, but shows good agreement with Monte Carlo particularly in the strong-coupling regime (40). Here we generalize the exact DMFT solution of Ref. (35) to the case of nonequilibrium DMFT.

## Ii Model and Methods

The Holstein model (21) is defined by the Hamiltonian

 H=−J∑⟨ij⟩(c†icj+h.c.)+∑iH(i)loc, (1) H(i)loc=ω0b†ibi+gni(bi+b†i)+ϵfni. (2)

The first term in Eq. (1) describes tight-binding electrons with nearest-neighbor hopping on a lattice; and are the creation and annihilation operator of an electron on lattice site , respectively. The local part (2) of the Hamiltonian represents one harmonic oscillator with frequency at each lattice site, i.e., a dispersion-less optical phonon mode, whose coordinate is linearly coupled to the electron occupancy ; defines the zero of the energy. We focus on the dilute limit, where correlations between electrons are negligible, so that expectation values of observables are proportional to the density and can obtained from the solution of the model with only one electron. The hopping is taken as a unit of energy, and times are measured in terms of . The results are obtained for a Bethe lattice with a semi-elliptic density of states .

To get an understanding of polaron formation in the Holstein model, the limit of isolated lattice sites (atomic limit) is a convenient starting point (41). In this limit, the presence of an electron on the site shifts the equilibrium position of the oscillator: omitting site-indices for convenience, the local part of the Hamiltonian can be rewritten as

 Hloc=ω02[(X+X0)2+P2]+(ϵf−EP)^n, (3)

where and are coordinate and momentum of the oscillator, respectively, , and

 EP=g2ω0, (4)

is the lowering of the ground state energy which defines the bare polaron binding energy. In the lattice model, the energy ratio distinguishes the regimes of weak-coupling () and strong coupling (). For strong coupling, self-localized electron states at energy at different sites are coupled by the hopping and form a band of delocalized polaronic states; it’s bandwidth is reduced with respect to the free bandwidth by the Frank-Condon factor, which takes into account the coherent motion of the lattice distortion with the electron, i.e., the overlap between the ground states and of the oscillator and the displaced oscillator (3) respectively. A second important scale for the Holstein model is the ratio , which distinguishes the adiabatic behavior (), in which the phonon is slow compared to the electron, from the non-adiabatic behavior (). In the adiabatic strong-coupling regime, the number of oscillator quanta in the phonon cloud proliferates (in the atomic limit, ), which makes the dynamics in this regime qualitatively distinct from the non-adiabatic regime.

To study polaron formation in time, we start the simulations from an initial state in which electrons and lattice are decoupled, and the mean kinetic energy of the electron is comparable to the free bandwidth, whereas the lattice temperature is low (). This initial state may be taken as a simple model for the situation immediately after electrons have been promoted into an empty valence band by photo-excitation from a conduction band, because the process of rapid inter-band excitation leaves the lattice unaffected up to a good approximation. The precise form of the initial electron energy distribution is not important for the subsequent dynamics as long as it is broad on the scale of the bandwidth, and we take it to be a hot electron distribution with electron temperature .

To monitor the dynamics of the model we compute the time-resolved photoemission spectrum, which can be obtained from the electronic Green’s function. In the low density limit, the relevant propagators for adding an electron () and removing an electron () are given by

 ˜G>(t,t′) =−iZ0TrN=0[e−βHci(t)c†i(t′)], (5) ˜G<(t,t′) =iZ1nelTrN=1[e−βHc†i(t′)ci(t)], (6)

where is the trace over the -electron sector, and . (Note that we have assumed translational invariance and normalized by the electron density , so that .) The photoemission spectrum as a function of probe time and energy is obtained from by partial Fourier-transform and convolution with the envelope of the probe pulse (42),

 I(ω,t)=∫dt1dt22πiS(t1)S(t2)eiω(t1−t2)˜G<(t+t1,t+t2). (7)

In equilibrium, is translationally invariant in time, so that is given by the convolution

 I(ω)=∫dω′A<(ω−ω′)|~S(ω′)|2, (8)

of the power spectrum of the probe pulse with the occupied density of states, . In addition to the photoemission spectrum, we will compute time-local observables, i.e., the kinetic energy per site, , as well as the average number of oscillation quanta in the phonon cloud (i.e., at a site occupied by an electron), (the expectation values are translationally invariant and normalized by the electron density).

We compute the dynamics of the Holstein model using the nonequilibrium generalization of DMFT (32). In the limit of low density, the solution can be made exact, yielding both Green’s functions (5) and (6). In equilibrium (35), computing the propagator is sufficient, because and are related by a fluctuation dissipation relation. For the nonequilibrium case, we thus have to reformulate the equations of Ref. (35) in real-time and provide additional equations for (or equivalently, one set of equations on the Keldysh contour). The resulting equations are Volterra integral equations whose numerical solution is controlled by the maximum number of phonons on each site; the computational effort increases however only linearly with , so that we can obtain converged results with for several tens of hopping times. To keep the presentation concise, the detailed formalism is explained in the appendix A.

## Iii Results

### iii.1 Weak coupling regime

The weak-coupling regime is rather well described by rate equations (see below), which can capture the cooling of the initial hot electron state by emission of phonons. Nevertheless it is illustrative to look at the corresponding DMFT solution, to contrast the behavior for strong-coupling below. Figure 1a and b show the relaxation of the kinetic energy and the phonon number for various coupling strength . After a short transient, the time-evolution of both quantities follows a monotonous relaxation, which becomes faster with increasing coupling strength. Similarly, the relaxation can be seen in the time-resolved photoemission spectrum (Fig. 1c). At early times, the occupied density of states reflects the initial hot electron state and is smeared over the full band. (In the uncorrelated equilibrium state, the occupied density of states is .) Subsequently, electrons reduce their kinetic energy by the emission of phonons, and spectral weight is concentrated closer to the lower band edge.

For weak electron-phonon coupling, relaxation phenomena at long times are captured by a kinetic equation (43), which is also in agreement with exact diagonalization studies (22); (31). For low lattice temperature (), an electron with band energy can only emit phonons, at a rate determined by the coupling and the density of (final) states,

 1τ(ϵ)=g2D(ϵ−ω0). (9)

This result is obtained from Fermi’s golden rule, or equivalently, the imaginary part of the equilibrium self-energy . The -dependence of the relaxation time is indeed confirmed by the DMFT results when one fits the time-dependence of the photoemission spectrum in a certain energy window with a simple exponential function (this will be analyzed further below, see the curve in Fig. 5d). Furthermore, from Eq. (9) one sees that a thermal equilibrium state can never be reached, because the density of states vanishes if the final energy is below the lower band edge. This phase-space effect can be seen explicitly in our data: At long times, the time-resolved photoemission spectrum remains shifted with respect to the spectrum of the equilibrium state at temperature (see dotted horizontal lines in Figs. 1d).

Finally, we note that due to the energy-dependent relaxation time, the long-time asymptotic behavior of averaged quantities is not necessarily exponential. This can be seen for the kinetic energy: For a density of states with a van-Hove singularity at the lower band edge (as for a three-dimensional lattice, or the semi-elliptic density of states used here), the rate Eq. (9) implies a power-law long-time asymptotic behavior of with . (For a one-dimensional density of states, one would expect an exponential decay (31).) This behavior is observed in the numerical data (see Fig. 1a, inset), which is a nice confirmation of the rate equation analysis.

### iii.2 Strong coupling regime: Overview

In the remainder of this paper we focus on the intermediate and strong coupling regime, where small polarons are formed in equilibrium. Figure 2 shows the relaxation of and for couplings to , and phonon frequencies and in the adiabatic and non-adiabatic regime, respectively. The sudden coupling of the electron and phonons leads to coherent oscillations, which are more pronounced for large . Furthermore, the absolute value of the kinetic energy becomes smaller with increasing , indicating a stronger localization of the carriers, and shows a pronounced enhancement of the phonon cloud. These effects provide a first glance at the crossover from intermediate to strong coupling. A further analysis of the photoemission spectrum (Fig. 3) will show that the observed dynamical results from a mixture of two different relaxation path, involving either delocalized and localized states.

In the adiabatic case, (Figs. 3a and b), we can distinguish several characteristic features in the photoemission spectrum: (i) A rapid decay of the weight at high energies (, ), starting from the broad distribution of the initial hot electron state. (ii) Buildup of spectral weight far below the lower edge of the free band (around ) within less than one period , and a beating of weight between this region and at the frequency . Finally, (iii), even though the oscillations are damped, the spectrum is still different from the spectrum in the thermal state at temperature (dashed line in panel b), and displays two peaks instead of a single polaron band. Other than at weak-coupling, the differences between transient and equilibrium spectra occur on energy scales considerably larger than . Spectra for the non-adiabatic regime () are shown in Figs. 3c and d: Coherent oscillations are reflected in a rigid-like shift of the occupied density of states, and a two-peak structure of the transient state is not observed.

To develop a physical understanding of these observations, we will perform an analysis in two directions: a comparison to the spectrum of an isolated site will allow us to single out characteristic spectral signatures of (excited) polaron states and show how they reflect the structure of the phonon cloud, and a momentum-resolved spectrum will distinguish contributions from polarons and delocalized electrons.

### iii.3 Atomic limit and spectroscopic signatures of excited polarons

In the atomic limit, the Holstein model can be solved analytically, both in and out of equilibrium, using a Lang-Firsov transformation (41) or it’s time-dependent generalization (34). Details of the solution are summarized in Appendix B. In the ground state, the polaron corresponds to the displaced oscillator [Eq. (3) with ], and the occupied density of states is given by a set of delta-peaks,

 A<(ω)=∞∑n=0P(n)δ(ω−EP−nω0), (10)

where the weights are given by the phonon number distribution in the polaron state. This result has an intuitive understanding: photoemission removes an electron from the bound state at energy and transfers the oscillator into it’s excited state with a probability which is given by the overlap of and the oscillator state before removing the electron, . At zero temperature, is the ground state of the oscillator (3) with , and is a Poisson distribution with mean . The corresponding photoemission spectrum, Eq. (8), already matches the lattice result quite accurately for the parameters of Fig. 3, as shown by the curves labelled in Figs. 4a and b. It is thus worthwhile to take the isolated site also as a starting point to analyze the peculiar double peak spectra of the non-thermal state after dephasing of oscillations transient state at . (The dephasing of oscillations is studies in more detail in Sec. III.5 below.)

At first sight, one may assume that a peak in which is shifted several multiples of with respect to the ground state polaron implies a highly excited state. We will now argue, however, that the two-peak structure of the spectrum in the adiabatic case can be taken as the characteristic signature of a low lying excited polaron state. For this purpose we compute the photoemission spectrum for an isolated site, assuming that the displaced oscillator is initially in it’s th excited eigenstate. In this case Eq. (10) still holds, with the phonon excitation energy in the delta function replaced by . The phonon distribution function of the exited state, , is given by

 Pm(n+m)=P0(n)n!m!(n+m)!L(n)m(γ2)2, (11)

where is the Poisson distribution of the ground state (), and is a generalized Laguerre polynomial (see Appendix B). In particular, we have , i.e., the distribution function is suppressed at (close to the maximum of ), which implies a double peak. In general the th polynomial has zeros, reflecting the probability distribution function of the oscillator coordinate. A comparison of these excited state spectra with the time-dependent spectra of the lattice model shows that the splitting of the two peaks in (Fig. 4a) or the width of the distribution (Fig. 4b) after the decay of the oscillations is well in agreement with the fact that a low lying excited polaron state () is reached. The main difference to the lattice result is a strong enhancement of the peak around in the adiabatic case, which will be analyzed in Sec. III.4 below.

Because in the atomic limit the photoemission spectrum reflects the number distribution function in the phonon cloud, it is interesting to analyze directly in the lattice model and see whether a similar relation can be established. The phonon-number distribution in the lattice, which is defined by the translation-invariant correlation function

 Pph(n,t)=1Lnel∑i⟨niδb†ibi,n(t)⟩, (12)

is plotted in Fig. 5 for various coupling strength in the adiabatic limit. Initially (at time zero, not shown), the distribution is a Boltzmann distribution . In the equilibrium state at coupling (dashed lines), the formation of a polaron is indicated by a peak at finite , which approaches the Poisson result for large , see Fig. 5d. The real-time data (solid lines in Fig. 5a-c) show an initial increase of phonon numbers (phonon states up to must be kept to simulate the dynamics in this regime). For the weaker coupling case (Fig. 5a), then evolves towards the equilibrium distribution. For couplings beyond a crossover scale (, ), where the polaron peak forms in equilibrium, a maximum which is shifted with respect to appears in addition to the zero-centered distribution (Figs. 5b,c). Comparison of with the position of the maximum of the distribution of the excited polaron states [Eq. (11)] for also confirms the previous finding that the polaron is transferred into an a low-lying excited state. A similar characterization of excited polaron states by their number distribution has also been discussed for an isolated Holstein impurity (23).

The relation (10) in the atomic limit would imply that the separation of the two maxima and in is related to the separation of two peaks in the photoemission spectrum by (up to the energy resolution of the probe pulse). This relation indeed holds quite accurately in the lattice, see Fig. 5d: the position of the maxima at large time () depends on coupling and time, but it quite accurately matches the value (open and filled blue circles in symbols in Fig. 5d). Hence the photoemission spectrum is a good measure for the phonon cloud also in the lattice model. In particular we note that in the adiabatic case excited polarons appear generically for couplings beyond crossover scale for polaron formation in equilibrium, and since the splitting is of the order of rather than the small scale , this feature could be taken to monitor the polaron crossover in experiment. On the other hand, it is interesting to see that no signature of the crossover is seen in the behavior of high-energy electrons. For this we integrate the spectrum over the high-energy part ( in Fig. 3a) and fit the result with an exponential function . The relaxation rate is a smooth function and almost linear with of over the whole crossover regime (see red line in Fig. 5d).

### iii.4 Disentangling free and bound states

We now focus on the marked asymmetry of the two peaks in the transient spectra, Fig. 3b. Because the peak at higher energy also roughly coincides with the energy of the lower band edge in the free band, one may assume that the additional weight of the peak at higher energy is due to a contribution from delocalized states. To confirm this picture, we look at the momentum-resolved photoemission spectrum , to show that the asymmetric contribution is localized in . is obtained from Eq. (7) by replacing the local Green’s function with the momentum-resolved Green’s function . With a momentum-independent self-energy, dependence on appears only via the electron dispersion , which extends from to for the semi-elliptic density of states. The local spectrum is simply .

In Fig. 6, is plotted for two different times. At early time one observes one maximum in for each (see white dotted line in Fig. 6a). The linear relation still reflects the behavior of free electrons. At later times, a flat band with two maxima and appears which reflects the polaron states (white dotted lines in Fig. 6b). The ratio of the two maxima, , is however strongly enhanced at ; it is , , and for , , , respectively. This confirms that the asymmetry of the two peaks in the -integrated spectrum indeed comes mainly from the region , and thus may be assigned to an additional contribution from delocalized states, which could not be disentangled from the upper polaron peak by the energy-resolved spectrum alone.

The presence of metastable delocalized states has long been predicted from semiclassical arguments (24); (25) from the existence of a potential energy barrier between delocalized and polaron states in the adiabatic potential . In high-dimensions (35), the latter is given by the sum of the classical energy cost for displacing the oscillator at one lattice site, and the corresponding lowering of the ground state due to the impurity with potential . Since the electronic ground state energy is not lowered if the impurity potential lies within the bandwidth, there is always an energy cost for creating small distortions, and thus an energy barrier for bringing the system into a self-trapped state. In infinite-dimensions, can be computed analytically (35). In weak coupling, slightly deviates from the zero-centered harmonic oscillator. A second minimum in appears for , and becomes the global minimum for , see inset Fig. 6b. Note that the scale is nicely in agreement with the crossover scale in Fig. 5, beyond which we observe the formation of excited polarons. The global minimum describes the ground state properties of the localized state, and the local minimum at corresponds to a delocalized state in the semiclassical picture.

### iii.5 Coherent oscillations

In this section we will finally discuss the initial coherent oscillations which follow the coupling of the electrons to the lattice and the resulting sudden displacement of the oscillator zero. In the non-adiabatic regime, oscillations are reflected in a rigid-like shift of the band (Fig. 3c). One can see that this is the behavior expected for a single oscillator: In the atomic limit, the Green’s function for a sudden switch-on of the coupling can be obtained exactly; it is related to the time-translationally invariant equilibrium one [, with Eq. (10)] by a simple time-dependent factor (see Appendix B),

 ˜G<(t,t′) =˜G

In the photoemission spectrum, Eq. (7), the oscillating factor roughly acts like a shift of the probing frequency by when the probe pulse is shorter than , so that the in . The resulting photoemission spectrum is shown in Figs. 4c and d. (Longer pulses, which average over many cycles, would lead to time-independent bands split by .)

From the comparison of Fig. 3 with Figs. 4c and d it is apparent that only in the non-adiabatic regime does the lattice result reflect the coherent oscillations found in the atomic limit. This shows a qualitative difference between the two regimes. In the adiabatic regime, the same bare polaron binding corresponds to a larger number of phonon energy quanta. An electron can thus easily emit several phonons to neighboring sites, so that vibrational dephasing occurs already on the timescale of one phonon-period. In the non-adiabatic regime, in contrast, the total excitation energy corresponds to very few oscillator quanta right from the beginning, so that emission of phonons is restricted by phase space effects and the system remains in long-lived beating oscillations, which is in agreement with results from exact diagonalization (22); (23).

## Iv Conclusion

In conclusion, we have obtained the numerically exact solution of the single-electron Holstein model within nonequilibrium DMFT. The results provide a comprehensive picture how an excited “hot” electron distribution relaxes due to optical phonons, both at weak and strong coupling, and in the adiabatic and non-adiabatic regimes. Most important are the results for small phonon frequencies (adiabatic regime) and strong coupling, where polaronic states are expected in equilibrium. After a quick dephasing of initial coherent oscillations, the system reaches a state in which excited polarons coexist with metastable delocalized states. While we cannot resolve the final relaxation to the ground state (the time range of our simulations extends to several phonon periods), the observed transient features are expected to be important for a photo-induced metallic state at strong-electron-phonon coupling. (In fact, in real systems the lifetime of the entire photo-induced state may be shorter than the final equilibration time.)

Moreover, we discuss how the photoemission spectrum reflects properties of the phonon cloud and can thus be used to characterize the transient state: Excited polarons lead to a characteristic double-peak structure of the almost flat (i.e., weakly momentum dependent) polaron band. Delocalized states, on the other hand, can be identified because their distribution is peaked in momentum space. Nonequilibrium polarons and metastable delocalized states appear beyond a well-defined polaron crossover scale. At the same time, no signature of the crossover is seen in the relaxation behavior of high-energy electrons. This suggest that the high-energy relaxation rates can be used in experiment to estimate the coupling by a analysis in terms of Fermi Golden rule (44) even in the regime where small polarons are formed.

As far as a comparison is possible, our results are in qualitative agreement with earlier predictions, and with results for low-dimensional systems: A beating between excited polaron states in the non-adiabatic case is in agreement with exact diagonalization results for one dimension (22); (23). The dynamics of the strong coupling adiabatic regime most difficult to describe in a quantum mechanical lattice calculation. A barrier for relaxation from delocalized states to self-trapped states was predicted by semiclassical arguments (24); (25), and it is in agreement with the occurrence of a level anti-crossing between localized and delocalized ground states in the energy spectrum (22).

Even though the simple Holstein model is not directly applicable to many experiments, the coexistence of long-lived polarons and metastable delocalized states may be qualitatively correct for systems which at the moment do not allow for a simple modeling. In fact, the coexistence of a Drude peak and polaronic features in photo-excited states has been observed in optical experiments on TaS (12). If delocalized states are stabilized by an energy barrier, this suggests unique ways to control the properties of photo-excited states: The number of mobile carriers may be modified by second pulse that helps to bring electrons over the barrier, either by field-localization of the electrons, which can transiently increase the electron-lattice effects (45), or by exciting the delocalized carriers. In this way the carrier mobility could be lowered by a pulse, allowing for a controlled switch-on of a metallic state (by photo-exciting carriers), followed by a switch-off (by localizing carriers). Such possibilities will be investigated in future work. From a technical perspective, we note that the structure of the DMFT equations in equilibrium (a continued fraction) is similar to the momentum averaged technique (39). Hence the Green’s function formalism presented in our work can be directly applied to extend the latter approach to nonequilibrium, which would be a promising way to study the time-resolved optical conductivity of the transient state in finite dimensions (40).

###### Acknowledgements.
We thank J. Bonča, U. Bovensiepen, D. Golež, Z. Lenarčič, P. Prelovśek, Ph. Werner, and L. Vidmar for fruitful discussions. ME acknowledges the Aspen Center for Physics and the NSF Grant No. 1066293 for hospitality during writing of the manuscript. The calculations were run in part on the supercomputer HLRN of the North-German Supercomputing Alliance.

## Appendix A Nonequilibrium DMFT for the Holstein model

### a.1 Nonequilibrium DMFT for the Holstein model

In this appendix we present the exact solution of the nonequilibrium DMFT equations for the Holstein model in the low density limit. For an introduction to the Keldysh formalism, as well as the notation for Keldysh Green’s functions, self-energies and Dyson equations, we refer to Ref. (32). Nonequilibrium DMFT for the Holstein model has been discussed in Ref. (34), so we will only briefly outline the general formalism and then focus on the low-density limit. In DMFT, the lattice model Eq. (1) is mapped onto a single impurity model with action

 S=−i∫Cdt[Hloc(t)−μc†c]−i∫Cdtdt′c†(t)Δ(t,t′)c(t′) (15)

on the Keldysh contour (see Fig. 7). The action describes coupling of one “Holstein atom” to a bath of non-interacting electrons by the hybridization function ( is the chemical potential). From the impurity problem one obtains the local contour-ordered Green’s functions

 G(t,t′)=−i1ZTr[TCeSc(t)c†(t′)]. (16)

The hybridization function is determined self-consistently; we will use the closed form self-consistency relation , which corresponds to a Bethe lattice with a semi-elliptic density of states . In general, the action with the hybridization function is equivalent to an Anderson impurity Hamiltonian in which the impurity site is coupled to a number of bath orbitals with a certain choice of the parameters and , i.e., may also be computed with action . (For the mapping in non-equilibrium, see Ref. (46)). In the following discussion we will switch between the action and the Hamiltonian notation as appropriate. The final result can always be written in terms of the hybridization function, so the precise time-dependence of the parameters is not needed.

### a.2 Green’s functions in the low density limit

To implement the low density limit, one takes the limit , keeping only leading terms in the fugacity . With the equivalent impurity action, the contour-ordered Green’s function (16) can be written as

 G(t,t′)= Extra open brace or missing close brace (17)

by taking the term , which commutes with , out of the integral. One can then perform an expansion in powers of by separating the trace in contributions from particles, . The result is written as

 G(t,t′)= e(it−it′)μ× [ΘC(t,t′)+ξΘC(t′,t)][˜G(t,t′)+O(ξ)], (18)

where the factor takes care of the fact that the leading contribution is if is later than on and otherwise, and is given by

 ˜G(t,t′) =−iZ0TrN=0[U(−iβ,t)cU(t,t′)c†U(t′,0)], (19) ˜G(t,t′) =iZ0TrN=1[U(−iβ,t′)c†U(t′,t)cU(t,0)], (20)

for and , respectively; is the standard time-evolution operator which is given by the equation

 i∂tU(t,t′)=Himp(t)U(t,t′);U(t,t)=1. (21)

Because the normalization factor is the average particle number to leading order in , one can see that the propagators (5) and (6) are given by Eqs. (19) and (20), respectively.

In the following we refer to the function which contains the leading terms of in the low-density limit as the projected Green’s function. Before discussing the interacting case, it is useful to have a brief look at properties of the projected Green’s function in the noninteracting case (). For the action (15), is given by the standard Dyson equation on the contour , where is the convolution . By using the ansatz (18) for and one can show that the projected functions and satisfy the integral-differential equation

 (i∂t−ϵf)˜G0−[˜Δ↻˜G0](t,t′)=0, (22)

to be solved for with an initial condition , where is the cyclic convolution, i.e., the convolution integral is restricted to the range in which the variables appear along in clock-wise order,

 ∫t′C,tdsf(s) =∫t′≻s≻tdsf(s) t′≻t, (23) ∫t′C,tdsf(s) =∫t′≻s≻−iβdsf(s)+∫0+≻s≻tdsf(s) t≻t′. (24)

The Dyson equation for the projected functions has been derived and discussed in great detail in relation to the nonequilibrium generalization of the non-crossing approximation Ref. (47); the latter can also be obtained as the low-density limit of a (pseudo-particle) theory, and hence the mathematical structure of the Green’s function theory is the same. The numerical solution of the integral equation is also discussed in Ref. (47).

### a.3 The interacting case

To obtain a solution for the interacting projected Green’s function we insert an identity for the electron sector in Eqs. (19) and (20) to the left of the annihilation operator ; in the sector, a full basis is given by the phonon-number states . To re-group the resulting terms, it is convenient to introduce a cyclic propagator,

 Uc(t′,t) Unknown environment '% (25) =TCexp(−i∫t′C,tdsHimp(s)), (26)

and auxiliary propagators

 Gpp′(t,t′) =−i[ΘC(t,t′)−ΘC(t′,t)]⟨p|cUc(t,t′)c†|p′⟩. (27)

With these definitions it is straightforward to re-group Eqs. (19) and (20) into

 ˜G(t,t′)=1Z0∞∑p=0⟨p|Uc(t′,t)|p⟩Gpp(t,t′). (28)

The factor satisfies

 ⟨p|Uc(t,t′)|p′⟩=δpp′e−i(t−t′)pω0[ΘC(t,t′)+ΘC(t′,t)e−βpω0]. (29)

The next step is to derive the Dyson equation for . For this purpose, we evaluate by expanding in terms of the electron-phonon interaction ,

 Gpp′(t,t′)=∞∑n=0(−i)n+1∫tC,t′dt1∫t1C,t′dt2⋯∫tn−1C,t′dtn××⟨p|TCe−i∫CdsH0(s)cHep(t1)⋯Hep(tn)c†|p′⟩. (30)

Here is the noninteracting part of , and the operator acts on the one-electron sector, in which it can be expressed in terms of phonon number states,

 Hep =∑mm′c†|m⟩Xmm′(t)⟨m′|c, (31)

where is given by

 Xmm′ =g(t)√m+1δm′,m+1+g(t)√mδm′,m−1. (32)

Hence we have

 Gpp′(t,t′)=∞∑n=0(−i)n+1∫tC,t′dt1∫t1C,t′dt2⋯∫tn−1C,t′dtn∑m1,m′1,⋯,mn,m′nG0,pm1(t,t1)Xm1m′1G0,m′1m2(t1,t2)⋯×Xmnm′n(tn)G0,m′np′(tn,t′), (33)

where is the noninteracting resolvents. Eq.(33) could be shortened into

 G =G0+G0↻X↻G0+G0↻X↻G0↻X↻G0⋯ =G0+G0↻X↻G, (34)

where all objects are matrices, and is the cyclic convolution defined above.

This matrix-integral equation has to be solved for the diagonal elements . Before doing this, we look at the noninteracting resolvent . Since electrons and phonons decouple, one can see from the definition that is the product , where is the projected Green’s function for the pure electrons, which satisfies the projected Dyson equation (22). Because is just an exponential factor [cf. Eq. (29)], it is easy to show that

 (i∂t−ϵf−pω0)G0,pp′−[˜Δp↻G0,pp′](t,t′)=0 (35)

with the initial condition , where

 ˜Δp(t,t′)=˜Δ(t,t′)⟨p|Uc(t,t′)|p⟩. (36)

Hence we can bring Eq.(34) to differential form by acting with from the left,

 (i∂t−ϵf−pω0)Gpp′− [˜Δp↻Gpp′](t,t′) − ∑a=±1Xp,p+a(t)Gp+a,p(t,t′)=0, (37)

with the boundary condition . This equation has a tridiagonal structure. Like any matrix equation of that form, the equations for diagonal elements can be derived in recursive form. (For example, a similar recursive structure is solved in the non-equilibrium variant of inhomogeneous DMFT (48).) We summarize the results:

 (i∂t−ϵf−pω0)Gpp −[(˜Δp+˜Ap+˜Bp)↻Gpp](t,t′)=0, (38)
 ˜Ap(t,t′)=pg(t)˜G[p]p−1(t,t′)g(t′), (39) (i∂t−ϵf−(p−1)ω0)˜G[p]p−1(t,t′) −[(˜Δp+˜Ap−1)↻˜G[p]p−1](t,t′)=0, (40)
 ˜Bp(t,t′)=(p+1)g(t)˜G[p]p+1(t,t′)g(t′), (41) (i∂t−ϵf−(p+1)ω0)˜G[p]p+1(t,t′) −[(˜Δp+1+˜Bp+1)↻˜G[p]p+1](t,t′)=0, (42)

where initial conditions are . Solving Eqs.(38, 39, 40, 41, 42), consistently and plugging the solution into Eq.(28) by tuning such that the , enable us to come up with an exact numerical solution of a single-polaron problem.

Phonon-occupation numbers [Eq. (12)] can be read off directly from the Green’s function,