Real-time decay of a highly excited charge carrier in the one-dimensional Holstein model
F. Dorfner, L. Vidmar, C. Brockt, E. Jeckelmann, F. Heidrich-Meisner
RReal-time decay of a highly excited charge carrier in the one-dimensional Holsteinmodel
F. Dorfner, L. Vidmar, ∗ C. Brockt, E. Jeckelmann, and F. Heidrich-Meisner Department of Physics and Arnold Sommerfeld Center for Theoretical Physics,Ludwig-Maximilians-Universit¨at M¨unchen, D-80333 M¨unchen, Germany Institut f¨ur Theoretische Physik, Leibniz Universit¨at Hannover, Appelstrasse 2, D-30167 Hannover, Germany (Dated: Draft of June 30, 2018)We study the real-time dynamics of a highly excited charge carrier coupled to quantum phononsvia a Holstein-type electron-phonon coupling. This is a prototypical example for the non-equilibriumdynamics in an interacting many-body system where excess energy is transferred from electronic tophononic degrees of freedom. We use diagonalization in a limited functional space (LFS) to studythe non-equilibrium dynamics on a finite one-dimensional chain. This method agrees with exactdiagonalization and the time-evolving block-decimation method, in both the relaxation regime andthe long-time stationary state, and among these three methods it is the most efficient and versatileone for this problem. We perform a comprehensive analysis of the time evolution by calculating theelectron, phonon and electron-phonon coupling energies, and the electronic momentum distributionfunction. The numerical results are compared to analytical solutions for short times, for a smallhopping amplitude and for a weak electron-phonon coupling. In the latter case, the relaxationdynamics obtained from the Boltzmann equation agrees very well with the LFS data. We also studythe time dependence of the eigenstates of the single-site reduced density matrix, which defines so-called optimal phonon modes. We discuss their structure in non-equilibrium and the distribution oftheir weights. Our analysis shows that the structure of optimal phonon modes contains very usefulinformation for the interpretation of the numerical data.
PACS numbers: 71.38.-k, 71.38.Ht, 05.70.Ln, 05.10.-a
I. INTRODUCTION
The energy transfer from charge carriers to their en-vironment represents one of the focal points of researchin the dynamics of condensed-matter systems. The exci-tation of the electronic subsystem and the measurementof its relaxation in real time may provide valuable in-sight into fundamental processes in strongly correlatedsystems [1], and it may also be exploited in energy con-version devices such as solar cells [2]. The interest in non-equilibrium systems with strong correlations has recentlybeen intensified by a rapid development of time-resolvedexperiments, which enable one to follow the evolution ofthe excited system from the early stage (of the order of afew electronic time units [3–5]) up to the final, thermal-ized state. One of the central mechanisms of relaxation isthe energy transfer from charge carriers to phonons sincethe latter type of excitations is ubiquitous in solids.The precise role of phonons in dynamics is neverthe-less a very subtle issue that depends on many details ofthe underlying many-body system. In some cases, thecoupling of phonons to the charge carriers is so strongthat polaronic effects represent the major many-body ef-fect in the dynamics. These effects were studied exper-imentally in the context of polaron formation (or moregenerally, self-trapping of charge carriers) [6–13], and po-laron transport [14–16]. Phonons also represent a keyingredient to understand the dynamics of photo-excited ∗ E-mail: [email protected] charge carriers in semiconductors [17]. In addition, evenif electronic correlations strongly influence the many-body spectrum, phonons may still provide the dominant(and the fastest) relaxation channel. For example, thekey role of phonons for relaxation has been conjecturedfor one-dimensional (1D) Mott insulators [18, 19] (wherethe spin relaxation channel is inefficient [20] due to spin-charge separation), and has been discussed in the contextof the mechanism of ultrafast demagnetization of ferro-magnets [21, 22]. On the contrary, recent results for two-dimensional systems with antiferromagnetic correlationshave shown that relaxation due to coupling to antiferro-magnetic spin excitations can be very fast [23–26], andon a short time scale, these excitations can absorb moreenergy than phonons [27].To clarify the role of phonons in non-equilibrium sys-tems in more detail, the main theoretical questions thatmotivate our investigation are: (i)
How efficient is theenergy transfer to phonons, depending on the character-istic energy scales of the electrons and phonons and theelectron-phonon coupling strength? (ii)
What is the rel-evant time scale for the energy transfer to phonons? (iii)
What is the relation between closed and open systemswith respect to dissipation and what is the dynamicsbeyond the weak-coupling regime? In connection withthat, to which extent is the knowledge of the unitarytime evolution required to describe the dynamics of aquantum many-body system, or in which cases are semi-classical approaches sufficient? (iv)
How to numericallytreat quantum many-body systems with bosonic degreesof freedom far away from equilibrium? a r X i v : . [ c ond - m a t . s t r- e l ] M a r Various directions of the recent developments of non-equilibrium techniques offer a rich perspective to an-swer the last question. However, before the era of high-precision time-resolved experiments, most of the methodswere developed to address transport properties withinsemi-classical approaches (see, e.g., Ref. [28]). Methodsretaining full quantum coherence have recently been de-veloped to study transport [29–31] and relaxation [32–35]of isolated carriers coupled to phonons. (Note that ex-act solutions can be obtained in the 1D adiabatic regimeand for the linear electronic dispersion [36, 37].) Sev-eral techniques are also now available to study many-electron systems beyond the mean-field [38] and Boltz-mann approaches [39, 40]. For example, dynamics ofelectrons coupled to phonons has been studied withinthe Holstein model using dynamical mean-field theory(DMFT) [41], continuous-time quantum Monte Carlo [42]and Keldysh Green functions within the Migdal approxi-mation [43, 44]. In addition, by extending the problem toelectronic correlations, studies of the Hubbard-Holsteinand the t - J Holstein model have been performed usingmethods based on exact diagonalization [27, 45–47], thedensity-matrix renormalization group [18] and dynamicalmean-field theory [48, 49].In this work we apply wavefunction-based methodsto study the non-equilibrium dynamics of a coupledelectron-phonon system. Already on the level of ground-state calculations, the maximal number of local phonons N max has been identified as the bottleneck for efficientsimulations. An efficient truncation of the phonon ba-sis therefore represents a crucial step to overcome therapid growth of the Hilbert space [50–52]. Density-matrix renormalization group (DMRG) algorithms cantreat much larger system sizes than exact diagonaliza-tion, but they also scale unfavorably in N max (see, e.g.,Refs. [53, 54] for the Holstein model and Refs. [55–58] forthe Hubbard-Holstein model). Here we focus on a sin-gle charge carrier, which, as the main advantage, allowsfor the numerically reliable treatment of the time evolu-tion in a broad parameter regime. We study the HolsteinHamiltonian H = − t (cid:88) j (cid:16) c † j c j +1 + c † j +1 c j (cid:17) + (cid:126) ω (cid:88) j b † j b j (1) − γ (cid:88) j (cid:16) b † j + b j (cid:17) n j , where b j ( c j ) annihilates a phonon (electron) at site j and the local electronic density on site j is n j = c † j c j .We set (cid:126) ≡ et al. [59] to describe the FIG. 1. (Color online) Sketch of the initial condition, Eq. (2),and the time evolution: We start from the state with oneelectron at momentum k = π and no phonon. The totalenergy is thus equal to the initial electronic kinetic energy E kin ( t = 0) = (cid:15) π = 2 t , where (cid:15) k = − t cos k . Due tothe coupling to phonons, the electron loses energy by excitingphonons while moving through the lattice, which also resultsin the redistribution of its momentum occupations. Holstein polaron ground state. We apply the method ona finite lattice and show that it allows for the efficientsimulation of dynamics in both the relaxation regime aswell as in the long-time stationary regime. This comple-ments a previous work using the same method [34], wherethe relaxation regime on an infinite lattice was studied.Second, we want to analytically and numerically de-scribe the non-equilibrium dynamics in limiting cases andin the crossover from weak to strong electron-phonon cou-pling, as well as in the crossover from adiabatic to anti-adiabatic regime. To be specific, we are interested in thedynamics emerging from an initial state with all excessenergy contained in the electronic sector | ψ (cid:105) = c † K |∅(cid:105) ele ⊗ |∅(cid:105) ph , (2)where c † k represents the creation operator for an elec-tronic state with momentum k (see also Fig. 1). Wechoose k = π , which also sets the total crystal momen-tum K of the coupled electron-phonon system. The twostates in Eq. (2) represent the electron and the phononvacuum, respectively. This state is an eigenstate of thesystem when the electron-phonon coupling is zero and isone of the simplest initial states where the charge carrieris highly excited. The physical motivation for choosingthe state in Eq. (2) is to model the dynamics of a free-electron wave-packet [32, 33, 35] emerging after a suddenexternal perturbation of a many-body system. Prevent-ing density fluctuations (i.e., by using the fully delocal-ized initial state with a sharp momentum) simplifies theanalysis in the stationary regime, discussed in Sec. V B 2,since our state is homogeneous for all times. We definetwo distinct regimes of the time evolution: (i) The re-laxation regime , which is characterized by a net transferof energy between the electronic and phononic system.For some parameter regimes, coherent oscillations withthe period 2 π/ω emerge in the dynamics. In such cases,the relaxation regime is defined by a nonzero net energytransfer averaged over one oscillation period. However,to have a meaningful definition of relaxation and energytransfer, we need to require that the amplitude of oscil-lations is sufficiently small compared to the mean valueof energy transfer. (ii) The stationary regime , in whichno net energy transfer is observed. This definition coversboth the case where the observables indeed become timeindependent, as well as the case with persisting coherentoscillations (i.e., when the amplitude of oscillations doesnot decay to zero). Such terminology differs from otherthermalization studies where temporal fluctuations aboutthe time average are usually required to be arbitrarilysmall to obtain a stationary state (see, e.g., Refs. [61–64]and references therein).As a third goal, we extract the single-site reduced den-sity matrix from the time-dependent total wavefunction.Motivated by the idea of the DMRG algorithms [65–67],we calculate the von Neumann entropy and the eigen-system of the single-site density matrix. In the electron-phonon coupled systems, the corresponding eigenstateswith the largest eigenvalues represent the optimal phononmodes . Zhang et al. [50, 51] first used the optimalphonon modes to truncate the phonon Hilbert spaceself-consistently in a study of the Holstein model inequilibrium. Later, this idea was applied to the spin-Peierls model [68] as well as to the Holstein-Hubbardmodel [52, 69]. Here we analyze the optimal phononmodes in a non-equilibrium set-up. We do not use themas a tool for truncating the Hilbert space, but ratherto gain additional insight into the dynamics. We showthat the structure of the optimal phonon modes carriesvaluable information about physical processes, and maycomplement the analysis based on observables only.
Summary of the results.
One of the central measuresof whether the energy transfer between the electron andphonons is substantial or not is the adiabaticity parame-ter, i.e., the ratio between the phonon energy ω and theelectronic bandwidth, set by the electronic hopping am-plitude t . If this ratio is sufficiently small, relaxation dy-namics can be observed. We construct a set of Boltzmannequations to describe the relaxation and compare the re-sults to the numerical data at weak electron-phonon cou-pling, obtaining a very good agreement. The Boltzmannequation also provides a convenient framework to ana-lytically derive the characteristic relaxation time τ for aconstant density of states. Its quantitative value is givenby τ ω = (16 /π )( γ/t ) − , and interestingly, it also rea-sonably well describes the numerical results obtained forthe density of states of a one-dimensional tight-bindingsystem.On the contrary, in the anti-adiabatic limit ω (cid:29) t thedynamics in both weak- and strong-coupling regimes isgoverned by coherent oscillations with the period 2 π/ω and a negligible energy transfer. We have applied per-turbation theory in both limits γ (cid:28) t and γ (cid:29) t andobtained very similar results for the time dependence ofobservables. The paradigmatic case to describe the oscil-latory dynamics is the single-site problem t = 0: here, the time evolution of observables can be obtained ana-lytically for all times. While the frequency of oscillationsis given by ω , the amplitude of the oscillations is set bythe ratio γ/ω .We then numerically calculate the time-dependent ex-pectation values of observables in the crossover from adi-abatic to anti-adiabatic regime as well as from weak tostrong electron-phonon coupling. The generic evolutionof the electronic momentum-distribution function is con-sistent with the sketch presented in Fig. 1. It undergoesstrong redistributions for the initial state in Eq. (2): ini-tially, all the weight is located at k = π , while at latertimes a maximum at k = 0 develops, accompanied by areduction of the electronic kinetic energy. The station-ary regime is typically characterized by the maximumof the momentum-distribution function being at k = 0and persistent coherent oscillations of the observablesabout their average; however, we show that the oscil-lations vanish with increasing lattice size for sufficientlyweak electron-phonon coupling and away from the anti-adiabatic limit. Note that both the electronic kinetic en-ergy and the electron-phonon coupling energy are lowerin the stationary regime in comparison to its initial value,hence both contributions should be accounted for whenquantifying the energy redistribution. We construct aheuristic estimate, supported by the numerical data, ofthe average number of emitted phonons in the stationaryregime: they equal the total energy minus the sum ofthe kinetic and the coupling energy in the correspondingground state. The time evolution can therefore simplybe viewed as the energy transfer from the subsystem,containing both the electronic and the electron-phononcoupling energy, to the phononic subsystem.We finally calculate the entanglement entropy and theeigensystem of the single-site reduced density matrix toinvestigate how the optimal phonon modes evolve in time.For the optimal phonon modes with the largest weights,the weight distribution remains approximately exponen-tial in non-equilibrium as found for ground states, how-ever, with a much slower decay. The structure of themost relevant optimal phonon modes is very interestingat strong electron-phonon coupling. In particular, theoptimal phonon mode with the largest weight exhibitstime-dependent oscillations between the phonon-vacuumand a state with a Poisson-like distribution over phononoccupation numbers. This structure strongly resemblesthe results obtained in the t = 0 limit, suggesting thatthe single-site coherent oscillations govern the dynamicsalso when none of the model parameters is vanishinglysmall.Note that the same problem, albeit with slightly dif-ferent initial states and in the limit of an infinite co-ordination number, has very recently been studied bySayyad and Eckstein using a non-equilibrium DMFTmethod [70]. Their main interest is in the adiabaticregime of phonon frequencies and in an extended analy-sis of photoemission spectra. Their results for the weak-coupling regime are in agreement with our discussion. Outline.
The plan of the paper is the following. InSec. II we introduce the model and discuss its limitingcases. A comparison of the numerical methods, i.e., exactdiagonalization, diagonalization in a limited functionalspace and the time-evolving block-decimation method,is performed in Sec. III. We then discuss the dynam-ics in limiting cases by applying perturbative techniquesin Sec. IV, and we investigate the relaxation dynamicswithin the Boltzmann approach. In Sec. V, we studythe time evolution and steady-state properties of sev-eral observables around the crossover from adiabatic toanti-adiabatic regime as well as from weak to strongelectron-phonon coupling. Finally, we focus on optimalphonon modes and the single-site entanglement entropyin Sec. VI. We present our conclusions in Sec. VII.
II. THE HOLSTEIN MODEL
The Holstein model [71] is a widely studied model incondensed matter theory, mimicking the interaction ofcharge carriers with local lattice vibrations. Recently,the possibility of the quantum simulation of the Hol-stein model has been discussed in the context of ul-tracold bosons [72], polar molecules [73, 74], trappedions [75, 76], Rydberg atoms [77] and superconductingcircuits [78]. The Holstein model consists of three terms H = H kin + H ph + H coup , (3)which represent the electronic kinetic energy H kin , thephonon energy H ph and the electron-phonon coupling en-ergy H coup . We provide these terms both in real- andmomentum-space representations H kin = − t (cid:88) j (cid:16) c † j c j +1 + c † j +1 c j (cid:17) = (cid:88) k (cid:15) k c † k c k , (4) H ph = ω (cid:88) j b † j b j = ω (cid:88) q b † q b q , (5) H coup = − γ (cid:88) j (cid:16) b † j + b j (cid:17) n j (6)= − γ √ L (cid:88) q,k (cid:16) b † q c † k − q c k + h.c. (cid:17) . We are interested in a single-electron system only, hencethe filling is kept at n = 1 /L where L is the num-ber of sites on a 1D lattice. Unless stated other-wise, periodic boundary conditions are assumed. In themomentum-space representation, the annihilation oper-ator of a particle with momentum k is defined by adiscrete Fourier transformation d k = (cid:80) j e − ijk d j / √ L ,where d j ∈ { b j , c j } . The electron in Eq. (4) has thestandard tight-binding dispersion (cid:15) k = − t cos k, (7)while the phonons in Eq. (5) are dispersionless. We alsostudy the electronic momentum distribution function n k , which is defined as n k = (cid:104) c † k c k (cid:105) = 1 L (cid:88) j,l e i ( j − l ) k (cid:104) c † j c l (cid:105) . (8)Throughout the manuscript we will denote the expecta-tion values of H , H ph , H coup and H kin as E total , E ph , E coup and E kin , respectively.In the ground state of the model at non-zero electron-phonon coupling γ , the electron and a cloud of phononsare spatially correlated and they form a composite par-ticle with a renormalized effective mass. This quasi-particle is called a polaron [79, 80]. There is no ana-lytical solution for the Holstein polaron in the full pa-rameter regime, however, results from numerical simula-tions [81] and perturbative limits may be very instruc-tive [82]. The limit of weak electron-phonon coupling γ (cid:28) t leads to a ground state with a highly mobileelectron and very few phonons. Because of the largespatial extent of electron-phonon correlations, this caseis called the large-polaron limit. On the contrary, thestrong-coupling anti-adiabatic limit γ, ω (cid:29) t is char-acterized by a reduced mobility and a rapid decay ofelectron-phonon correlations with distance. Hence it asalso referred to as the small-polaron limit. In the ex-treme limit of zero hopping amplitude t = 0 one is leftwith a fully local Hamiltonian. This limit is particularlyinstructive since the Hamiltonian H = ω (cid:88) j b † j b j − γ (cid:88) j (cid:16) b † j + b j (cid:17) n j (9)can be diagonalized analytically by shifting the phononoperators at the site j of the electron via b j = a j + g, (10)which leads to the diagonal Hamiltonian H j = ω a † j a j − ε b . (11)Here, ε b represents the polaron binding energy ε b = γ ω (12)and we have introduced the shift given by g = γω , (13)which emerges as the only relevant model parameter inthe t → | ψ GS (cid:105) = e − g / √ L (cid:88) j (cid:104) e gb † j |∅(cid:105) ph ⊗ c † j |∅(cid:105) ele (cid:105) . (14)The average number of phonons in this state is N ph = (cid:80) j (cid:104) b † j b j (cid:105) = g . The transformation of operators definedby Eq. (10) corresponds to a unitary transformation tothe new orthogonal basis, known also as the coherent-state basis. We will elaborate in more detail on thistransformation in Sec. IV E where we derive an analyti-cal expression for the real-time dynamics in the limit ofsmall t . In addition, the coherent-state basis also repre-sents a very instructive framework for the discussion ofoptimal modes in Sec. VI.The intermediate regime of parameters is characterizedby a crossover from the large to the small polaron casewith increasing electron-phonon coupling. The proper-ties are usually characterized as a function of the dimen-sionless electron-phonon coupling parameter λ = γ t ω = ε b t . (15)This is the ratio between the ground-state energies in twoextreme cases: the free fermion energy − t , which rep-resents the ground-state energy at zero electron-phononcoupling, and the polaron binding energy − ε b , whichrepresents the ground-state energy at infinitely strongelectron-phonon coupling. Hence the crossover (i.e., thepoint where most observables exhibit the most rapid vari-ation) emerges roughly at λ ∗ ≈
1. The precise crossovervalue also depends on another dimensionless parameter α = ω /t [83], called the adiabaticity parameter, whichis used to distinguish the adiabatic limit ( α (cid:28)
1) fromthe anti-adiabatic limit ( α (cid:29) III. NUMERICAL METHODS
Electron-phonon lattice models such as the Holsteinmodel (3) do not conserve the number of phonons in thesystem and the Hilbert space is infinite even for a systemon a finite lattice. Wavefunction based methods usedto study the Holstein polaron (e.g., exact diagonaliza-tion [84–87] and DMRG [88]) are flexible and can yieldessentially exact results, however, they face an obviouslimitation since they require a finite Hilbert-space ba-sis. The central task is hence to efficiently construct atruncated basis that is able to represent the states oneis interested in. Typically, one truncates the total num-ber of phonons, however, the efficiency of the trunca-tion strongly depends on the regime of model parame-ters considered. In Sec. III A we introduce the conceptof diagonalization in a limited functional space, while inSec. III B we describe the TEBD method. Other methods widely used to study properties of the Holstein polaronare dynamical mean-field theory [89], different versionsof quantum Monte Carlo (e.g., diagrammatic [90, 91],continuous-time [92, 93], Lang-Firsov based [94] and vari-ational Monte Carlo [95]), and a momentum-average ap-proach [91, 96–99].
A. Diagonalization in a limited functional space
The main idea in this truncation scheme comes fromthe structure of the polaron: Since the electron andphonons are correlated in space, the phonons in the vicin-ity of the electron are more relevant than the phononsfar away from the electron. One should therefore find away to keep a large amount of distinct phonon configu-rations around the electron while neglecting some stateswith less important phonon configurations. This is ef-ficiently achieved by constructing a limited functionalspace (LFS), first introduced by Bonˇca et al. [59]. In-stead of the full basis, it only picks up a limited setof wavefunctions, i.e., configurations in the occupation-number basis, represented in a translationally invariantform. For a given size of the functional space, the sys-tem is then diagonalized exactly. The method has beenshown to be both very accurate and efficient, and it wasapplied to many different problems [27, 30, 34, 100–105](see also [35, 83, 106–109]).One of the main advantages of the method is that itprovides a systematic way of generating a limited Hilbertspace. The generation procedure is initiated by a simplestarting state, which is generally the bare electron in agiven momentum eigenstate (i.e., in a translationally in-variant state). Then, the off-diagonal elements of theHamiltonian are applied to the initial state. The max-imal number of generations of new states is labeled bythe parameter N h . We represent the entire set of states (cid:110)(cid:12)(cid:12)(cid:12) φ ( N h ,M ) k (cid:69)(cid:111) forming the LFS by a sum N h (cid:88) n h =0 (cid:32) H kin + M (cid:88) m =1 ( H coup ) m (cid:33) n h c † k |∅(cid:105) . (16)Beside N h , there is another parameter M that determinesthe size of the LFS: it counts the number of applicationsof the electron-phonon coupling term within a single gen-eration, i.e., it may increase the amount of phonons thatare created in each generation. In this work we definestates on a finite lattice, which is in contrast to the initialapplication of the method [59] where an infinite latticewas used. The spatial extent of the phonon cloud aroundthe electron is governed by N h . For the states formingthe LFS on a large lattice and at small N h , the largestrelative distance between the electron and a phonon is N h − N h > L , hence the states with all thepossible relative distances between the electron and a sin-gle phonon are incorporated in the LFS. The generationscheme is illustrated in Fig. 2 for N h = 3 and M = 1. FIG. 2. (Color online) Sketch of the states in the limitedfunctional space that are generated by Eq. (16) during thefirst three iterations (we take M = 1 for simplicity). The ringrepresents a lattice with L = 6 sites and periodic boundaryconditions. Each state in the figure represents the parent stateof a translationally invariant state. The site of the electron isrepresented by an arrow for clarity, while the local phononsare represented by harmonic oscillators. Starting from thefree electron state ( n h = 0), a single phonon (since M = 1) iscreated on the site of the electron. In the next generation, the n h = 1 state is taken as the starting state: the electron hopsto the left or right of the phonon (leading to two new states)and a state with an additional phonon is created. Next, thesecond generation ( n h = 2) states serve as starting statesand the new unique states of the third generation ( n h = 3)are added to the basis. This procedure is repeated until thedesired number of generations are set up (when a state isgenerated that already exists in an earlier generation, it isomitted from the generation). The method provides numerically exact results for theground-state properties of the Holstein polaron [59, 110],and can be applied to different parameter regimes. In theweak-coupling regime, the spatial extent of the ground-state electron-phonon correlations may be large. In ad-dition, in most excited states there are unbound phononsthat can be arbitrarily far away from the electron. Nev-ertheless, by systematically increasing N h it has beenshown that an accurate description of low-lying excitedstates is possible [104]. In the strong-coupling regime, thepolaron is very small and has a huge amount of phononsin its vicinity. This property can be captured by tun-ing the parameter M [102] since the maximal number ofphonons on the site of the electron is N h × M . Note that t.t E ph ( t ) /t D = 1711247, N h = 20D = 144898, N h = 16D = 19547, N h = 13D = 2392, N h = 100123456 LFS, D = 89992ED, D = 125970L = 12(a)(b) L = 8LFS FIG. 3. (Color online) Time evolution of the phonon energy E ph ( t ) for λ = 0 . ω = t . (a) Comparison of the numer-ical methods for L = 8. The solid line represents results ob-tained by exact diagonalization (ED) in the complete Hilbertspace, using the global phonon cutoff N max = 12. The dashedline represents results obtained by diagonalization in a lim-ited functional space (LFS) by using the parameters N h = 17and M = 1 [see Eq. (16)]. (b) Convergence of E ph ( t ) by usingLFS for L = 12. For the parameters of the basis generation,we keep M = 1 while N h is varied. The legends display thecorresponding basis dimension D . the generator of states in Eq. (16) does not represent theonly possibility to efficiently truncate the phonon num-ber; for alternative methods of creating distinct phononconfigurations see, e.g., Refs. [111, 112].Recently, the LFS has also been applied to study timeevolution under non-equilibrium conditions. Two spe-cific cases were studied: the Holstein polaron [30] orbipolaron [105] driven by a constant electric field, andthe Holstein polaron excited by a short pump-pulse [34].The LFS was constructed in both cases on an infinitelattice and hence in the limit of a large time, the quasi-stationary state may include states with a diverging rela-tive distance between the electron and the phonons. Thiseffect limits the largest times available by using the LFS.Here we generalize the method to deal with a finite L , anduse the L → ∞ limit only in a specific case in Sec. V B 2.The benefit of this generalization is many-fold: one can,e.g., time evolve the system for a longer time interval,efficiently calculate the electronic momentum distribu-tion function n k , and analyze the L -dependent scaling ofthe amplitude of coherent oscillations in the stationaryregime.We demonstrate the efficiency of the LFS method inFig. 3 for the weak-coupling regime at λ = 0 . ω = t by comparing to exact diagonalization in afull Hilbert space, subjected by a global phonon cutoff N max . For both methods, we time-evolve the wavefunc-tion | ψ ( t ) (cid:105) = e − iHt | ψ (cid:105) by using the iterative Lanczosalgorithm [113]. We use a time step of ∆ t = 0 . N max = 12 for a system size L = 8. Both data are in excellent agreement. However,when the electron-phonon coupling is further increased,a larger phonon cutoff is required, and consequently thefull Hilbert space exceeds the available computational re-sources while the structure of the LFS can be tuned byincreasing the parameter M in Eq. (16). Figure 3(b)shows the convergence of the time-dependent phonon en-ergy E ph ( t ) at λ = 0 . N h and M of the basis generator in Eq. (16). With in-creasing N h , the data converge fast. For the system size L = 12 shown in Fig. 3(b), the Hilbert space of dimen-sion D ∼ is sufficient for an accurate description ofdynamics in all time regimes. In all subsequent figureswhere the LFS results are shown, we set the parameters N h and M large enough such that the results are con-verged for all times. B. Time-evolving block decimation
Another method that we use is the time-evolving block-decimation (TEBD) algorithm [60], closely related to thetime-dependent DMRG method [114, 115]. The advan-tage of this method is that it can also treat the caseof large distances between the electron and (unbound)phonons, because the computational cost increases lin-early with the system size for the polaron problem. As itis based on the matrix-product-state (MPS) formalism,this algorithm is efficient only for representing states withsmall entanglement [116]. In this work we always startwith a slightly entangled state as the matrix dimensionis m = 2 for the initial state (2), but the entanglementincreases rapidly over time. Thus we can simulate thepolaron dynamics accurately only for a limited period oftime, which varies with the model parameters. In addi-tion, we use open boundary conditions because this yieldsa better performance with TEBD.The TEBD method can readily be applied to systemswith bosonic degrees of freedom represented by a bareboson basis with a cutoff N max (cid:28) ∞ as in the exactdiagonalization methods [117]. However, since the com-putational cost (time) of TEBD scales as N , we haveused only up to N max = 17 states per boson site and thuswe are limited to moderately large values of the coupling γ (cid:46) t =0 . ω . In addition to the left and right normaliza-tion conditions, we also use a representation of the MPSthat conserves the particle number explicitly. During thetime evolution we keep all eigenstates of the reduced den-sity matrices that correspond to eigenvalues greater than10 − until the maximal matrix dimension ( m = 50 or t.t E ph ( t ) /t TEBDLFSL = 12
FIG. 4. (Color online) Comparison of the phonon energy E ph ( t ) obtained by diagonalization in LFS and the TEBDmethod. We set λ = 0 . ω = t and L = 12. Open boundaryconditions are used in this figure (we use the free-electroneigenstate with the highest energy as the initial state). m (cid:48) = N max m up to m (cid:48) = 1700in TEBD simulations of XY Z spin-1 / − . For each timestep, the resulting truncation error is smaller than theexpected Trotter error. We estimate that all TEBD re-sults presented here have negligible errors on the figurescales.The result of the TEBD method is compared to diago-nalization within LFS in Fig. 4 for λ = 0 . ω = t ona open chain with L = 12 sites. The perfect agreementconfirms the accuracy of both methods for long times t , even for system sizes L exceeding the maximal sizethat can be treated with exact diagonalization methods.Since the entanglement growth makes the TEBD methodcomputationally more expensive as time evolves and thepresence of chain edges complicates the interpretation ofresults, we primarily use diagonalization in LFS in thefollowing to comprehensively analyze the dynamics forvarious parameter regimes. IV. PERTURBATIVE RESULTS
The goal of this section is to understand the dynam-ics in limiting cases, in which analytical expressions canbe obtained. There are two main benefits of this treat-ment: (i)
It provides a distinction between parameterregimes where relaxation dynamics sets in and regimeswhere the response is governed by coherent oscillations; (ii)
It clarifies the role of model parameters on dynamics,in particular, it determines the functional dependence ofthe relaxation time.We apply time-dependent perturbation theory to ob-tain the time evolution of observables and compare themto the numerically exact results. In Sec. IV A we applythe general procedure to extract the lowest orders of thetime-evolution operator, while in Sec. IV B we carry out aperturbative expansion for the case when one parameterof the Hamiltonian is much smaller than the others. Theresults of Sec. IV B are applied in Secs. IV C and IV Eto the cases of a weak electron-phonon coupling γ anda small hopping amplitude t , respectively. In addition,for weak electron-phonon coupling we construct the setof Boltzmann equations and compare them to the numer-ical data in Sec. IV D. In the analytical calculations, weassume periodic boundary conditions. A. Short-time dynamics
The time dependence of an operator O in the Heisen-berg picture can be obtained for short times by using theBaker-Campbell-Hausdorff (BCH) formulae X O e − X = O + [ X, O ] + 12! [ X, [ X, O ]] + · · · , (17)where X = iHt , and the n -th order expansion yieldsresults up to O ( t n ). For the expectation value of thephonon-energy operator, one gets E ph ( t ) = (cid:104) H ph (cid:105) − ω (cid:68) H (1)coup (cid:69) t (18)+ (cid:18) ω γ + ω (cid:104) H coup (cid:105) − ω (cid:68) H (2)coup (cid:69)(cid:19) t + O ( t ) , where (cid:104) . . . (cid:105) denotes expectation values in the initialstate. Here, we have introduced expectation values ofthe generalized electron-phonon coupling term H ( ϑ )coup = − γ √ L (cid:88) q,k (cid:16) M ( ϑ ) k,q b † q c † k − q c k + h . c . (cid:17) , (19)where M (1) k,q = i and M (2) k,q = (cid:15) k − (cid:15) k − q . Even though H (1)coup and H (2)coup do not appear in the Holstein Hamil-tonian (3), which governs the time evolution, their ex-pectation values may become relevant for the short-timedynamics. For our initial state introduced in Eq. (2),however, most of the expectation values in Eq. (18) van-ish and we get E ph ( t ) = ω ( γt ) + O ( t ) . (20)Similar derivations can be carried out for the otherterms in the Hamiltonian. For the kinetic energy, onegets E kin ( t ) = E kin (0) (cid:2) − ( γt ) (cid:3) + O ( t ) , (21)while for the coupling energy, the result is E coup ( t ) = [ E kin (0) − ω ] ( γt ) + O ( t ) . (22)This clearly conserves the total energy during the timeevolution, E total = E kin (0) = E ph ( t ) + E coup ( t ) + E kin ( t ).For the initial state considered in our work, the short-time dynamics is therefore controlled by 1 /γ . In Fig. 5 wecompare the short-time evolution of E ph ( t ) from Eq. (20)with the numerical results using LFS. The figure showsthat the short-time expansion describes the numericaldata well for γt (cid:46) .
5, while the agreement is the betterthe larger λ . t. γ E ph ( t ) /t Pert. theory λ = 0.5 λ = 1 λ = 2 λ = 4.5 FIG. 5. (Color online) Short-time dynamics of the phononenergy E ph ( t ). The bold-dashed line represents the resultsfrom perturbation theory, given by Eq. (20). All other curvesare numerical results using LFS at ω = t and different λ (we use L = 12 for λ = 0 . L = 8 for λ = 2 and4 . /γ . B. Perturbation theory in the interaction picture
If one model parameter (for the sake of generalitycalled η ) is much smaller than other ones, we split theHamiltonian into two terms H = H + ηV. (23)At short times, ηt represents a small parameter inthe argument of the time-evolution operator U ( t ) =exp ( − iHt ). It is convenient to expand U ( t ) in the in-teraction picture U ( t ) = U ( t ) + ηU ( t ) + η U ( t ) + O ( η ) , (24)where U ( t ) = e − iH t , (25) U ( t ) = − i (cid:90) t dt e − iH t V e − iH ( t − t ) (26)and U ( t ) = − (cid:90) t dt (cid:90) t − t dt e − iH t V e − iH t × (27) V e − iH ( t − t − t ) . Hence, the time-dependent expectation values can be ex-pressed in the initial state as [we omit the time variablein U i ( t )] O ( t ) = (cid:68) U † OU (cid:69) + η (cid:68) U † OU + U † OU (cid:69) (28)+ η (cid:68) U † OU + U † OU + U † OU (cid:69) + O ( η ) , which in many cases leads to accurate results in a longertime interval compared to the short-time expansion per-formed in Sec. IV A. C. Time evolution for weak electron-phononcoupling
Here we explicitly derive O ( t ) from Eq. (28) for thecase η = γ (cid:28) ω , t and our initial state (2). We performthe derivation up to O ( γ ). We then distinguish betweenthe anti-adiabatic and adiabatic regime, i.e., Sec. IV C 1refers to the case ω > t and Sec. IV C 2 to the case ω < t .For the first-order contribution we get (cid:68) U † OU + U † OU (cid:69) =2 √ L (cid:88) q R (1) K,q (cid:18) − cos ( δE K,q t ) δE K,q (cid:19) , (29)where δE K,q = (cid:15) K − q + ω − (cid:15) K (30)is the energy difference between one-phonon states (con-sisting of the electron with momentum K − q and aphonon with momentum q ) and the initial state (the elec-tron with momentum K and no phonon). Electronic en-ergies are given by the dispersion relation in Eq. (7). Inaddition, R (1) K,q = (cid:68) ∅ (cid:12)(cid:12)(cid:12) c K O (ele) c † K − q (cid:12)(cid:12)(cid:12) ∅ (cid:69) ele (cid:68) ∅ (cid:12)(cid:12)(cid:12) O (ph) b † q (cid:12)(cid:12)(cid:12) ∅ (cid:69) ph (31)represents the matrix element decomposed into the elec-tronic and the phononic part, and O ≡ O (ele) ⊗ O (ph) .The second-order contributions consist of two terms, (cid:68) U † OU + U † OU (cid:69) = − R (2 a ) K (cid:18) − cos ( δE K,q t )( δE K,q ) (cid:19) , (32)and (cid:68) U † OU (cid:69) = 2 L (cid:88) q R (2 b ) K,q (cid:18) − cos ( δE K,q t )( δE K,q ) (cid:19) . (33)The corresponding matrix elements are R (2 a ) K = (cid:68) ∅ (cid:12)(cid:12)(cid:12) c K O (ele) c † K (cid:12)(cid:12)(cid:12) ∅ (cid:69) ele (cid:68) ∅ (cid:12)(cid:12)(cid:12) O (ph) (cid:12)(cid:12)(cid:12) ∅ (cid:69) ph (34)and R (2 b ) K,q,q (cid:48) = (cid:68) ∅ (cid:12)(cid:12)(cid:12) c K − q (cid:48) O (ele) c † K − q (cid:12)(cid:12)(cid:12) ∅ (cid:69) ele (cid:68) ∅ (cid:12)(cid:12)(cid:12) b q (cid:48) O (ph) b † q (cid:12)(cid:12)(cid:12) ∅ (cid:69) ph , (35)where the identity R (2 b ) K,q,q (cid:48) = R (2 b ) K,q δ ( q, q (cid:48) ) has been as-sumed to derive Eq. (33) (this is the case for all operatorsconsidered in this work), and O (ph) was assumed to beat most linear in b in Eq. (34).The expressions for the matrix elements R (1) K,q , R (2 a ) K and R (2 b ) K,q are listed in Appendix A for the observables t. ω × -4 × -4 × -4 E ph ( t ) / ω Pert. theoryLFS
FIG. 6. (Color online) Time evolution of the phonon energy E ph ( t ) in the weak-coupling anti-adiabatic limit for L = 12.We set t /ω = 0 . γ/ω = 0 .
01, which corresponds to λ = 0 . × − . The solid line is the numerical result usingLFS, while the dashed line is the result from perturbationtheory, given by Eq. (37). of interest. For the expectation values of the terms thatappear in the Holstein Hamiltonian (3), we then get E kin ( t ) = (cid:15) K (cid:32) − γ L (cid:88) q − cos ( δE K,q t )( δE K,q ) (cid:33) (36)+ 2 γ L (cid:88) q (cid:15) K − q − cos ( δE K,q t )( δE K,q ) E ph ( t ) = 2 ω γ L (cid:88) q − cos ( δE K,q t )( δE K,q ) (37) E coup ( t ) = − γ L (cid:88) q − cos ( δE K,q t ) δE K,q . (38)Note that the expansion of the above expressions upto t matches the short-time results from Eqs. (20), (21)and (22). Nevertheless, the results from Eqs. (36)-(38)can be applied to a much larger time domain, as we aregoing to demonstrate in the following.
1. Weak-coupling anti-adiabatic regime
We first discuss the regime when the phonon energyis larger than the electronic bandwidth, ω > t [moregenerally, if the initial K (cid:54) = π , it is sufficient to require ω > t (1 − cos K )]. Hence the energy of a single quan-tum phonon exceeds the maximal electronic energy dif-ference, which already indicates the inefficiency of en-ergy transfer. To understand the relevant energy scalesin this regime, we first consider the anti-adiabatic limit ω (cid:29) t where we can replace δE K,q by ω in Eqs. (36),(37) and (38) and thus obtain simple formulas for thetime evolution of these energies. For the time-dependentphonon energy, this results in E ph ( t ) = 2 ω (cid:18) γω (cid:19) [1 − cos ( ω t )] . (39)There are two main observations from the above equa-tion: (i) The time evolution is governed by oscillations0with frequency ω ; (ii) The maximal number of emit-ted phonons approaches 4( γ/ω ) . This is qualitativelythe same result as in the strong-coupling anti-adiabaticregime ( ω , γ (cid:29) t ) to be discussed in Sec. IV E, eventhough here we assume a weak coupling γ (cid:28) t .Nevertheless, when 1 < ω / (4 t ) (cid:28) ∞ , Eq. (39) repre-sents only a poor approximation, and one should calcu-late the sum over q in Eq. (37) explicitly. In Fig. 6, wecompare the resulting phonon energy with the numeri-cally exact LFS data for t /ω = 0 . γ/ω = 0 . tω <
30. Since Eq. (37) has been derivedby taking into account transitions between the phononvacuum state and one-phonon states only, this impliesthat no higher-order processes are relevant for the dy-namics in this parameter regime.
2. Relaxation in the weak-coupling regime
One of the central goals of our study is to investi-gate the relaxation dynamics, i.e., the situation wherethe majority of the electronic energy is, at sufficientlylarge times, transferred to phonons. It is convenient toaddress this goal in the weak-coupling regime and for ω < t (1 − cos K ). We study the electronic relaxationby measuring the temporal decrease of the electronic ki-netic energy. One can notice from Eq. (36) that the ki-netic energy is reduced with respect to the initial valueby ∆ E kin ( t ) = (cid:80) q ( (cid:15) K − q − (cid:15) K ) n K − q ( t ), where n K − q ( t ) = 4 γ L sin (cid:16) δE K,q t (cid:17) ( δE K,q ) . (40)To address the relaxation, therefore, it is important tounderstand the time evolution of the momentum distri-bution function n k .We define the transition rate W = t (cid:80) k (cid:54) = K n k (∆ t ),which is the probability per time interval ∆ t for the tran-sition to electronic states different from the initial state.In principle, this rate is time dependent, W = W (∆ t ).For our problem, following the notation of Eq. (40), itcan be expressed as W (∆ t ) = 4 γ L (cid:88) q (cid:54) =0 t sin (cid:16) δE K,q ∆ t (cid:17) ( δE K,q ) . (41)In the limit of large L , one can replace the sum by anintegral and rewrite W in the energy representation byintroducing the dimensionless electronic density of states D ele ( E + (cid:15) − ω ) = 2 t d q d (cid:15) K − q = 1 (cid:114) − (cid:16) E + (cid:15) − ω t (cid:17) , (42)where E represents the electronic energy after the transi-tion and (cid:15) is the electronic energy before the transition (i.e., if we start from the initial state, (cid:15) = (cid:15) K ). Thetransition rate then equals W (∆ t ) = γ t (cid:90) t − ( (cid:15) − ω ) − t − ( (cid:15) − ω ) d E × (43) D ele ( E + (cid:15) − ω ) 2 π t sin (cid:0) E ∆ t (cid:1) E . In the last part of the expression, one can recognize thedelta function δ ( E ) = 2 sin (cid:0) E ∆ t (cid:1) / ( πE ∆ t ) for suffi-ciently large time ∆ t . Hence for short times, the electroncan make transitions to various different energy levelswith δE K,q (cid:54) = 0 , while for longer times, it transfers onlyto states with a well-defined energy satisfying the totalenergy conservation δE K,q = 0. The time scale is givenby the energy uncertainty principle ∆ E ∆ t > π .Provided that D ele ( E + (cid:15) − ω ) is a smooth functionaround E = 0, this results in a time-independent transi-tion rate W = γ t D ele ( (cid:15) − ω ) = 2 λω D ele ( (cid:15) − ω ) . (44)The transition rate is therefore proportional to the elec-tronic density of states after the transition, known asFermi’s golden rule [118]. Equation (44) suggests express-ing W in units of ω since λ and D ele are dimensionlessby definition. We will elaborate more on this issue inSec. IV D 2 where we show that it is also convenient touse the same unit to measure the characteristic time ofthe entire relaxation process.In the adiabatic limit ω → K is W K = γ t | sin( K ) | (45)in leading order. This result agrees up to a factor oftwo with the exact rate obtained in a 1D polaron modelwith a linear electronic dispersion [36]. The factor twois due to the presence of two electronic branches (andthus two relaxation paths) in the model (3) as opposedto the single branch of the linear-dispersion model. Wecan also understand why both results agree in the adia-batic limit only. On the one hand, our model, Eq. (3),and the linear-dispersion model of Ref. 36 are equivalentonly if the transfer of electronic momentum k is small.On the other hand, Eq. (40) shows that the transitionprobability is maximal for δE K,q = 0. Both conditionscan be fulfilled simultaneously only in the adiabatic limit ω (cid:28) t .However, in 1D systems the electronic density of statesdiverges at the band edges and a more careful analysisis required for transitions to the bottom of the electronicband, i.e., for (cid:15) − ω ≈ − t . In that case, Eq. (43)should be calculated explicitly. By expanding D ele ( E − t ) around E = 0, we get W (∆ t ) = γ √ t (cid:90) t d E π t sin (cid:0) E ∆ t (cid:1) E / . (46)1For 4 t ∆ t (cid:29) W (∆ t ) ≈ γ √ t √ √ π √ ∆ t. (47)This shows that the transition rate to the lowestelectronic state cannot be approximated by a time-independent quantity.In the special case when the single phonon energy ex-actly matches the maximal electronic energy difference, ω = 2 t (1 − cos K ), the result of Eq. (47) can be ap-plied directly to get the time evolution of observables.Consequently, we find an anomalous power-law behavior∆ E kin ( t ) ≈ − t (1 − cos K ) (cid:114) γt √ √ π ( γt ) / (48)for times t , ω (cid:29) t − (cid:29) γ .To summarize the discussion, the above equationswere derived within the second-order perturbation the-ory describing the single-phonon emission only. Hence, if ω (cid:28) t (1 − cos K ), the two-level transition described by Eq. (44) is followed by a cascade of transitions form-ing the complete quasi-particle relaxation process (seealso Fig. 1). A possible way to describe this relaxationis to discretize time in small time steps and assume thatat each step, the transition rates between the allowedelectronic levels are given by Eq. (44) [or, in case of tran-sitions to the bottom of the band, by Eq. (47)]. We aregoing to pursue this idea in the context of the Boltzmannequation to calculate the time evolution of the electronicmomentum distribution n k . Such dynamics is Marko-vian, but may still be a good approximation for weakenough electron-phonon couplings (for studies on the in-fluence of Markovian effects on polaron dynamics see,e.g., Refs. [119, 120]). D. Boltzmann equation
The Boltzmann equation is a semi-classical approachfor the non-equilibrium dynamics of electronic distribu-tion functions in the thermodynamic limit. In a systemwithout any external force or inhomogeneity, the changeof the momentum distribution comes solely from colli-sions. The corresponding set of equations for the 1Delectron-phonon system is˙ n k = (cid:88) q W k,q [( n k − q (1 − n k ) N q − n k (1 − n k − q )( N q + 1)) δ ( (cid:15) k − (cid:15) k − q − (cid:126) ω q ) (49)+ ( n k + q (1 − n k )( N q + 1) − n k (1 − n k + q ) N q ) δ ( (cid:15) k − (cid:15) k + q + (cid:126) ω q )] . The first half of the right-hand-side term describes tran-sitions between the electronic state with momentum k accompanied by N q phonons with momentum q and the(lower-energy) electronic state with momentum k - q andone more phonon. The second half of the right-hand-sideterm describes transitions between the electronic statewith momentum k accompanied by N q + 1 phonons withmomentum q and the (higher-energy) electronic statewith momentum k + q and one less phonon. The matrixelement W k,q represents the transition rate for these pro-cesses. In the following, we will rewrite Eq. (49) in the en-ergy representation n k → n (cid:15) and use the transition ratesfrom Eqs. (44) and (47). In a 1D tight-binding model,there are two distinct k -points at the same energy (theexceptions being the top and the bottom of the band).This was already taken into account in the derivation ofthe transition rates in Sec. IV C 2, and the correspondingrenormalization should be (1 − n k ) → (1 − n (cid:15) / s = 4 t /ω transitions between s + 1 electronic lev-els. We assume, without loss of generality, that s is aninteger.
1. Numerical solution
We now solve the Boltzmann equation (49) for the Hol-stein polaron problem (3) with our initial condition (2).We set ω q = ω since our phonons are dispersionless. Ourinitial state is a vacuum for phonons, and we thereforeset N q = 0 at t = 0 and assume the phonons to remainin that equilibrium state throughout the time evolution.This implies that the electron never scatters again offphonons that it has excited, which is a reasonable as-sumption when dealing with a single electron in a largeempty lattice. In total, there are s + 1 equations for s transitions, where in each transition, the electronic en-ergy is lowered by ω ,˙ n (cid:15) = − W n (cid:15) (1 − n (cid:15) / n (cid:15) = − W n (cid:15) (1 − n (cid:15) /
2) + W n (cid:15) (1 − n (cid:15) / . . . (50)˙ n (cid:15) s − = − W s n (cid:15) s − (1 − n (cid:15) s ) + W s − n (cid:15) s − (1 − n (cid:15) s − / n (cid:15) s = W s n (cid:15) s − (1 − n (cid:15) s ) . The nearest energy levels are hence related by (cid:15) m +1 = (cid:15) m − ω . For 1 ≤ m < s , the transition rate is W m = ( γ /t ) D ele ( (cid:15) m ). In the last equation, W s is time-2 t/ τ -2-1012 E k i n ( t ) /t LFSBoltzmann γ /t = √ ω = t FIG. 7. (Color online) Relaxation of the electronic kinetic en-ergy E kin ( t ) at ω = t in the weak-coupling regime ( γ/t ) =0 .
4, which corresponds to λ = 0 .
2. The solid line is the nu-merical result using LFS (for L = 12), the dashed line isthe result from the numerically integrated Boltzmann equa-tions (50). Time is measured in units of τ , defined in Eq. (53). dependent and given by Eq. (47).When solving Eq. (50) numerically, some care isrequired when choosing the appropriate discrete timestep. On the one hand, in the derivation of the time-independent values of W m [see Eq. (43)] we require thatthe time step is not too small; on the other hand, theequations must necessarily fulfill W m ∆ t <
1, which im-plies ∆ tω (cid:46) (2 λ ) − (we assumed D ele ∼ tω =1. In Fig. 7 we compare the electronic kinetic energy E kin ( t ) = (cid:80) k (cid:15) k n k ( t ) obtained by the Boltzmann ki-netic equation, with the numerical solution using LFS.Remarkably, the results show perfect agreement in therelaxation regime. This indicates that, as long as theelectron-phonon coupling is sufficiently small, the pro-cesses described by Eq. (50) provide the most relevantcontribution to the dynamics. Deviations set in whenthe steady state is approached: there, the dynamics ona finite lattice is governed by the interplay between theelectron and the phonons emitted during the relaxation.Clearly, such processes are not included in the Boltzmannequations (50), and thus an agreement between the solu-tion of the Holstein model and the Boltzmann equationis not expected in the steady state. We mention thatrecently, by using the Boltzmann approach, the mobilityof the Holstein polaron was accurately reproduced in theweak-coupling regime [121].
2. Relaxation for a constant density of states
The occupation of different electronic k -states can alsobe obtained analytically under some approximations. Wecan linearize the Boltzmann equations (50) because wehave only one electron in the system and we can assumea constant energy density D ele = π/
4, which gives theenergy-independent transition rate Γ = ( π/ λω (thiswould be exact for a linear electronic dispersion [36]). t/ τ -2-1012 E k i n ( t ) /t ω /t = 0.5 ω /t = 0.8 ω /t = 1 ω /t = 1.2 ω /t = 1.5 γ /t = √ FIG. 8. (Color online) Relaxation of the electronic kinetic en-ergy E kin ( t ) in the weak-coupling regime. Results are shownfor a fixed γ and t [setting ( γ/t ) = 0 .
4] and different val-ues of ω /t = 0 . , . , . , . , .
5, which correspond to λ =0 . , . , . , . , .
13, respectively. We use LFS for L = 12.The dashed line is the function E kin ( t ) /t = 2[1 − t/τ )] withthe relaxation time τ defined in Eq. (53). Then one needs to solve equations of the form˙ n (cid:15) = − Γ n (cid:15) ˙ n (cid:15) = − Γ n (cid:15) + Γ n (cid:15) . . . (51)˙ n (cid:15) s − = − Γ n (cid:15) s − + Γ n (cid:15) s − ˙ n (cid:15) s = Γ n (cid:15) s − . The solution is, for a fixed 0 ≤ m < s , n (cid:15) m ( t ) = (Γ t ) m m ! e − Γ t , (52)which is a Poisson distribution and has a maximum at m ≈ Γ t . To obtain the characteristic relaxation time τ ,we simply require τ ≡ s/ Γ. This definition provides areasonable estimate for the time when the lowest single-electron state (i.e., the one at k = 0) becomes dominantlyoccupied. It gives τ ω = 16 π (cid:18) γt (cid:19) − = 8 π (cid:18) λ ω t (cid:19) − . (53)The main observation from this result is that in gen-eral, the quantitative value of the characteristic relax-ation time depends on all the three parameters of theHolstein Hamiltonian. Nevertheless, an important in-sight from Eq. (53) is that 1 /ω is a very convenient timeunit to measure relaxation. The advantage of measuringtime as tω is twofold: first, since the derivation stemsfrom the weak-coupling regime, τ ω is always larger thanone, and second, the quantitative value of the relaxationtime in these units is then only a function of the electron-phonon coupling energy in units of the electron delocal-ization energy. In addition, the result (53) also affirmsthat the relaxation time can not be deduced from theshort-time expansion in Eq. (21).We test the prediction of the relaxation time (53) us-ing the numerical results for the Holstein Hamiltonian.In Fig. 8 we plot E kin ( t ) for different ω , but for fixed3values of γ and t . When time is measured in units of ω , the curves fairly well collapse on the same line in therelaxation regime. Remarkably, the simple estimate (53),using a constant density of states, provides a reasonablemeasure for the relaxation time. Physically, the reasonfor this agreement is that the relaxation is dominatedby the slowest transitions W m in the Boltzmann equa-tions (50), i.e., when the electron is away from the bandedges and their divergent density of states. Hence eventhough both initial and final electronic states are at theband edges, a decent estimate for the relaxation time canalready be obtained from the assumption of a constantdensity of states. E. Time evolution for small hopping amplitude
We now consider the case t (cid:28) γ, ω . This is thelimit where the net energy transfer between the electronand phonons is expected to be small or even does not oc-cur. Instead, the time dependence of observables exhibitsstrong oscillations. Analytical solutions for the lowestnon-zero orders in η = t allow us to show explicitly thatthe period of oscillations is given by ω , while the ampli-tude of oscillation is governed by the ratio g = γ/ω .
1. Single-site dynamics
For t = 0 the unperturbed Hamiltonian H is givenby Eq. (9). As all lattice sites are decoupled, it can bediagonalized by means of a rotation to the coherent-statebasis [122] as shown in Sec. II. Therefore, the exact an-alytical solution can be obtained for all times. To carryout a perturbative expansion in t , however, it is moreconvenient to rotate the Hamiltonian operator instead ofthe basis. Explicitly, any operator O can be rotated bya unitary Lang-Firsov transformation (cid:101) O = e iS Oe − iS (54)with the hermitian operator S = g (cid:88) j p j n j , (55)where p j = i ( b † j − b j ) represents the local oscillator mo-mentum operator. The Hamiltonian is then diagonal (cid:101) H = ω (cid:88) j b † j b j − ε b (56)with the polaron binding energy ε b defined in Eq. (12)and the initial state (2) is now represented by | (cid:101) ψ (cid:105) = e iS | ψ (cid:105) (57)= e − g / √ L (cid:88) j e ijK (cid:104) e − gb † j |∅(cid:105) ph ⊗ c † j |∅(cid:105) ele (cid:105) . The time evolution of this state can be calculated exactly | (cid:101) ψ ( t ) (cid:105) = e − i (cid:101) Ht | (cid:101) ψ (cid:105) (58)= e iε b t e − g / √ L (cid:88) j e ijK (cid:104) e − g ( t ) b † j |∅(cid:105) ph ⊗ c † j |∅(cid:105) ele (cid:105) with g ( t ) = ge − iω t . Clearly, this state describes a de-localized composite quasi-particle made of the electrondressed by a phonon cloud. This is qualitatively similarto the small polaron found in the ground state of theHolstein model. However, here the phonon cloud fluc-tuates with time in contrast to the static cloud of theground-state polaron.Finally, time-dependent expectation values can be cal-culated directly in this representation using O ( t ) = (cid:104) (cid:101) ψ ( t ) | (cid:101) O | (cid:101) ψ ( t ) (cid:105) . (59)Using the BCH formula (17) we easily get (cid:101) x j = b † j + b j + 2 gn j (60)for the local oscillator displacement operator x j = b † j + b j and (cid:101) p j = p j for the corresponding momentum operator.Thus their expectation values are x j ( t ) = 2 gL [1 − cos ( ω t )] (61) p j ( t ) = 2 gL sin ( ω t ) . (62)This result can be easily understood. Each site has aprobability 1 /L to be occupied by the electron. If thesite is occupied, the boson degree of freedom represents aharmonic oscillator with frequency ω and an equilibriumposition (cid:104) x (cid:105) eq = 2 g . Our initial condition (2) correspondsto x j ( t = 0) = p j ( t = 0) = 0. Thus, the oscillator swingsbetween x = 0 and x = 4 g like a classical harmonic oscil-lator. Similarly, for the phonon energy H ph = ω (cid:80) j b † j b j we obtain (cid:101) H ph = ω (cid:88) j b † j b j + gω (cid:88) j (cid:16) b † j + b j (cid:17) n j + ε b (63)and then the expectation value E ph ( t ) = 2 g ω [1 − cos ( ω t )] . (64)This is identical to the result (39) obtained in the anti-adiabatic limit of the weak-coupling perturbation expan-sion, i.e., for γ (cid:28) t (cid:28) ω , although we have assumed t (cid:28) γ, ω to derive the present result. Thus in the anti-adiabatic limit, Eq. (64) seems to hold in both weak- andstrong-coupling regime. This is also confirmed by ourTEBD simulations for t /ω = 10 − and finite γ . We alsonote that the number of phonons N ph ( t ) = E ph ( t ) /ω reaches a maximum value 4 g as a function of time, i.e.,four times larger than the polaron ground-state value N ph = g in the anti-adiabatic strong-coupling regime.In all observables of interest, the oscillations never decay4 t. ω -0.500.511.52 E k i n ( t ) /t Pert. theoryt / ω = 0.001t / ω = 0.1 FIG. 9. (Color online) The electronic kinetic energy E kin ( t )in the strong-coupling anti-adiabatic regime at γ/ω = 1 and L = 12. The dashed line is the result from perturbationtheory, given by Eq. (68), while the solid and thin dashedline represent the TEBD result for t /ω = 10 − and 10 − (corresponding to λ = 500 and 5, respectively). in time and their amplitude is determined by the bareelectron-phonon coupling γ expressed in units of ω . Allresults also unambiguously show that 2 π/ω or integerfractions thereof set the period of these oscillations. Thetime evolution therefore imposes harmonic oscillations ofthe system around its polaron ground state.
2. Finite but small hopping amplitude
If the hopping amplitude is finite but still t t (cid:28) η = t (cid:28) t − . The timeevolution of the first-order term of the kinetic energy E kin ( t ) = (cid:88) k (cid:15) k n k ( t ) (65)can also be obtained from zeroth-order expectation val-ues (59). Using (cid:101) c † k (cid:101) c k = 1 L (cid:88) j,l e ik ( j − l ) e ig ( p j − p l ) c † j c l (66)we first calculate the electronic momentum distributionin zeroth order resulting in n k ( t ) = 1 L + e g [cos ( ω t ) − (cid:18) δ k,K − L (cid:19) (67)and then obtain E kin ( t ) = − t cos( K ) e g [cos( ω t ) − . (68)This result is plotted in Fig. 9. We see that the kinetic en-ergy is also a periodic function with period 2 π/ω that os-cillates between its initial value 2 t cos( K ) and an expo-nentially reduced value exp( − g )2 t cos( K ). This con-firms that the system oscillates without any relaxationup to leading orders in t in the limit of small hoppingterms t . Since the result (68) is valid for both open and E ph ( t ) / ( t ) ω /t = 10 ω /t = 6 ω /t = 5 ω /t = 2 ω /t = 1 t. ω -2-1012 E k i n ( t ) /t ω /t = 10 ω /t = 6 ω /t = 5 λ = 0.2(a)(b) FIG. 10. (Color online) Time evolution for the two partsof the Hamiltonian: (a) the phonon energy divided by thebandwidth E ph ( t ) / (4 t ) and (b) the electron kinetic energy E kin ( t ), plotted for several values of the adiabaticity ratio ω /t . We use LFS at λ = 0 . L = 12. periodic boundary conditions, we can compare it directlyto TEBD simulations using an initial standing wave withthe wave number K = πL/ ( L + 1). We observe a per-fect agreement between Eq. (68) and the TEBD resultfor t /ω = 10 − , as shown in Fig. 9. We note thatthe kinetic energy scale is a tiny fraction of the otherenergy scales, i.e., E kin ( t ) (cid:46) − ω , ε b , but it is nev-ertheless perfectly reproduced by the TEBD data. Thisshows that TEBD simulations are very accurate in thisregime. When the ratio t /ω increases, the numericalresults and perturbation theory cease to perfectly agreewith each other. However, we show in Fig. 9 that at t /ω = 0 .
1, Eq. (68) still provides the correct qualita-tive behavior of the dynamics.
V. NUMERICAL RESULTS FORINTERMEDIATE PARAMETER REGIMES
In this section we focus on non-perturbative results byapplying diagonalization in the LFS. This complementsSec. IV that has addressed the regimes where one pa-rameter of the Hamiltonian is vanishingly small. Thecentral question here is to which degree the properties ofthe weak- and strong-coupling limits as well as adiabaticand anti-adiabatic limits persist at intermediate valuesof parameters.5
A. Crossover from adiabatic to anti-adiabaticregime
The solution of the second-order perturbation theoryin the weak-coupling anti-adiabatic regime discussed inSec. IV C 1 yields oscillatory behavior of all parts of theHamiltonian [see Fig. 6 for E ph ( t )]. In this case relax-ation barely takes place since the net energy transferfrom electron to phonons is negligible and comparableto the amplitude of oscillations. On the other hand,the adiabatic regime at weak electron-phonon couplingis reasonably well described by the Boltzmann equation,see Sec. IV D. This is the paradigmatic case for relax-ation since the majority of the electronic kinetic energyis transferred to phonons. The question remains howthese two regimes evolve into each other when the ratio ω /t is varied.A convenient way of addressing this is to study theweak-coupling regime where γ is smaller than ω and4 t . In this case, the amount of emitted phonons inthe adiabatic regime roughly compensates the reduc-tion of the electronic kinetic energy. In Figs. 10(a)and 10(b) we plot, for λ = 0 .
20, the phonon energy di-vided by the bandwidth E ph ( t ) / (4 t ) and the electron ki-netic energy E kin ( t ), respectively. At ω = t we observe E ph ( t ) / (4 t ) ≈ ω (cid:29) t , E ph ( t ) / (4 t ) < ω (cid:28) t , dom-inated by relaxation of the electronic kinetic energy, to-wards a different type of behavior at ω > t governedby coherent oscillations and only weak redistribution ofenergies. B. Crossover from weak to strong coupling at ω = t In the following we focus on the role of the electron-phonon coupling at ω = t and study the dynamicsat the crossover from weak to strong coupling. We de-scribe the features in two different time regimes, the re-laxation regime and the stationary regime, in Secs. V B 1and V B 2, respectively.
1. Relaxation regime
The relaxation regime is the regime in which excessenergy is transferred from the electron to the phonon de-grees of freedom. This occurs in the first stage of the timeevolution, i.e., in the time interval roughly given by thecharacteristic relaxation time. In Sec. IV, perturbativearguments allowed us to clearly distinguish between thelimiting cases where relaxation is expected to take placeand where it is inefficient: a sufficient condition for the
FIG. 11. (Color online) Density plot showing the electronicmomentum distribution function n k ( t ) using LFS at ω = t and L = 12. (a) λ = 0 .
5, (b) λ = 2. first scenario is that the phonon energy is much smallerthan the electronic bandwidth. The decrease of kineticenergy is a consequence of the redistribution of the elec-tronic momentum distribution, discussed in the contextof the Boltzmann equation in Sec. IV D. The simplestpicture of momentum redistribution (relevant, in partic-ular, in the weak-coupling regime) has already been dis-cussed in the context of Fig. 1: the electron reduces itskinetic energy by emitting phonons, where each of theprocesses conserves the total energy and momentum. InFig. 11(a) we show the evolution of the electronic mo-mentum distribution function n k as a function of time atmoderate coupling λ = 0 . ω = t . The electronstarts off with a momentum k = π due to our initial con-dition (2). The momentum redistribution can be seen atshort times in the density plot of Fig. 11(a) as a narrowbright region expanding roughly linearly towards k = 0.When the maximum occupation is at k = 0 [this occursat tt ≈ λ = 0 . ω = t in Fig. 11(a)], therelaxation is completed and the system enters into thestationary regime. Note that in Sec. IV D we have de-fined the characteristic relaxation time τ as the time atwhich n k develops a maximum at k = 0 [see the discus-sion below Eq. (52)].When the electron-phonon coupling is increased, acrossover from the large to the small polaron regimetakes place in the ground state of the Holstein model [81].The dimensionless electron-phonon coupling λ , Eq. (15),represents a natural measure of the crossover: λ ∗ ≈ E ph ( t ), E coup ( t )and E kin ( t ) for λ = 0 .
5, 1 and 2 at fixed ω /t = 1. Interms of the energy transfer from the electron to phonons,no drastic features occur at λ ≈
1, and the relaxationtime of the electronic energy E kin ( t ) in Fig. 12(c) con-tinuously decreases with increasing λ . Larger differencescan be observed in the stationary regime, where oscil-lations of E ph ( t ) and E coup ( t ) become noticeable in thestrong-coupling regime (this will be discussed in moredetail in Sec. V B 2). The momentum distribution func-tion n k at λ = 2 [see Fig. 11(b)] also exhibits oscillations6 E ph ( t ) /t λ = 0.5 λ = 1-15-10-5 E c oup ( t ) /t λ = 2 t.t -2-101 E k i n ( t ) /t (a)(b)(c) FIG. 12. (Color online) Time evolution for the three compet-ing parts of the Hamiltonian: (a) phonon energy E ph ( t ), (b)electron-phonon coupling energy E coup ( t ), and (c) electronickinetic energy E kin ( t ). We set ω = t and show results usingLFS for the three parameter values λ = 0 .
5, 1 and 2 (we use L = 12 for the first two systems and L = 8 for the last one). in the stationary regime, however, it remains peaked at k = 0 for most of the time.Figure 12 also addresses a parameter regime that hasnot been covered by our perturbative analysis in Sec. IV:this is the regime where the electron-phonon coupling islarge while the phonon energy is smaller than the elec-tronic bandwidth (c.f. λ = 2 and ω = t ). In such a case,one can still observe the relaxation regime at short timesin the time evolution of the kinetic energy, Fig. 12(c),which starts at E kin ( t = 0) = 2 t and becomes negativeat sufficiently large time. In contrast to the kinetic en-ergy, E ph ( t ) and E coup ( t ) clearly exhibit oscillations withthe period 2 π/ω and the entire energy transfer takesplace already within a single period.
2. Stationary regime
In the introductory part, we defined the stationaryregime as the regime where no net energy transfer takesplace between the electron and the phonon system. If thetime dependence of observables exhibits persistent coher-ent oscillations, the latter statement corresponds to theaverage over the oscillation period (in most cases, thisequals 2 π/ω ).One of the main questions about the oscillations inthe stationary regime is whether they persist in the limit L → ∞ . To investigate this, we define the variance oftemporal fluctuations of an observable A ( t ) about the σ ∆ E ph ( L ) Linear fitLFS f ( L ) = 1.45 (1/ L ) f ( L ) = 1.10 + 1.89 (1/ L )(a) λ = 0.5(b) λ = 2 FIG. 13. (Color online) Standard deviation σ ∆ E ph [seeEq. (69)] of the phonon energy E ph ( t ) in the stationaryregime, as a function of the inverse lattice size. Circles repre-sent results using LFS, while the dashed lines in (a) and (b)are the linear fitting functions f ( L ) with one and two free pa-rameters, respectively. (a) λ = 0 .
5, where oscillations vanishin the thermodynamic limit (the two points for the smallest L were not included in the fit). (b) λ = 2 .
0, where oscillationspersist for all L . average in the stationary regime σ A = 1 t − t t (cid:90) t (cid:0) A ( t ) − ¯ A (cid:1) dt, (69)where t and t are chosen to be at the minimum andmaximum of the oscillations, and ¯ A represents the timeaverage. Figures 13(a) and 13(b) show the standard de-viation of temporal fluctuations of the phonon energy at λ = 0 . λ = 2, respectively. In the first case, theoscillations vanish with the system size, while at λ = 2,the oscillations seem to persist for any lattice size. Sinceanalytical arguments in Sec. IV E yield the time depen-dence with undamped oscillations in the strong-couplinganti-adiabatic limit t →
0, it is not surprising that awayfrom this limit, the oscillations still govern the dynamicsfor all times.It is nevertheless still an intriguing question how muchinformation about the dynamics is already included inthe solution of the strong-coupling anti-adiabatic limitin Sec. IV E. In Fig. 12, we have shown that the ki-netic energy at λ = 2 is already close to zero in thestationary state. Because of energy conservation, thisimplies that the phonon and the coupling energies os-cillate with the same phase and a similar amplitude, E ph ( t ) = − E coup ( t ) + E total , consistent with the resultfrom Eq. (64). In Figs. 14(a) and 14(b) we compare nu-merical results for E ph ( t ) at t = ω to the analyticalresult at t = 0 [see Eq. (64)] at the same ratio γ/ω .We take λ = 2 and 4 . t. ω E ph ( t ) / ω Single siteLFS(a) λ = 2(b) λ = 4.5 FIG. 14. (Color online) Time evolution of the phonon energy E ph ( t ) in the strong-coupling regime. The solid lines representnumerical results using LFS for t = ω and L = 8, whilethe dashed lines represent the single-site dynamics ( t = 0),given by Eq. (64). (a) λ = 2, which corresponds to g = 2. (b) λ = 4 .
5, which corresponds to g = 3. correspond to g = 2 and 3, respectively. The oscilla-tion period 2 π/ω matches almost perfectly in all cases,suggesting its universality for phonon-related quantities.However, the amplitude of oscillations at t = 0 (givenby 2 g ) is larger than at t = ω . There are, in fact,two effects when λ increases: both the amplitude of os-cillations and the average value approach the result (64).In Sec. VI we will discuss the optimal modes emergingfrom the single-site reduced density matrix, which repre-sent a complementary approach to study the structure ofphonon modes in non-equilibrium. We will further con-firm that the oscillations at finite λ are indeed relatedto the single-site coherent oscillations characterizing the t → E ph ( t ), E coup ( t ) and E kin ( t ) in the stationary regime as afunction of electron-phonon coupling λ (if the observableoscillates in the stationary regime, we take the averageover a sufficiently large time interval). We also com-pare these values to their respective ground-state valuesfor the Holstein polaron with the same total momentum, K = π . Results are shown in Fig. 15. While kineticand coupling energies in general relax to values close totheir respective ground-state values, this is not the casefor the phonon energy [see Fig. 15(b)]. Obviously, sincethe total energy of our closed system exceeds the ground-state energy of the Holstein polaron, it can not happenthat all three contributions to the total energy relax totheir ground-state values. The important message com-ing from Fig. 15 is therefore that practically all the ex-cess energy is stored in the form of excess phonons in thestationary state. The precise quantitative analysis of de- -10-50 ( E c oup + E k i n ) /t λ E ph /t t-AVGGSE kin (t=0)-E ( λ )-E ( λ ) 0 1 2 λ -8-4 E c oup /t λ -2-10 E k i n /t (a)(b) FIG. 15. (Color online) Time-averaged values (t-AVG) of thethree competing parts of the Hamiltonian in the stationarystate as a function of electron-phonon coupling λ at L = 12and ω = t , obtained by using LFS. The values are comparedto the polaron ground state (GS) energies E kin , E coup and E ph at K = π . (a) Main panel: E kin + E coup versus λ inthe ground state (diamonds) compared to the correspondingtime-averaged values in the stationary state (circles). Theinsets display the comparison of the individual terms. Thedata at λ → L -dependence due to the largespatial extent of the polaron. (b) Phonon energy E ph versus λ . Since the averages in the stationary state clearly exceedthe ground-state values, we also compare to E kin ( t = 0) − ( E , coup + E , kin ), where E kin ( t = 0) = E total , and E , coup , E , kin are measured in the ground state. viations of the kinetic and the electron-phonon couplingenergy with respect to their ground-state values also rep-resents an interesting aspect, which is, however, beyondthe scope of this paper.The latter result allows us to address a more ambitiousquestion: If we know the energy of the initial state, andif we know the ground-state energies for the system forgiven model parameters, is it possible to estimate thenumber of phonons at asymptotically long times withoutperforming the time evolution? Or, in other words, doesthe numerical data support a simple rule of thumb topredict ¯ N ph ? Figure 15(b) shows that such a heuristicestimate is indeed possible. The number of phonons is,to a good approximation, given by¯ N ph = E total − ( E , kin + E , coup ) ω , (70)where E , kin and E , coup are the ground-state energyof the kinetic and electron-phonon coupling part of theHamiltonian. Even though the latter result appears tobe very simple, it contains a potentially non-intuitive as-pect: If the electron-phonon coupling is not small, it isnot correct to view the time evolution solely as the energytransfer from the electron (i.e., the electronic kinetic en-8ergy) to phonons (i.e., the phonon energy E ph ). Instead,a considerable part of energy is stored in the couplingenergy as well, and the appropriate number of excessphonons can only be obtained if the electronic and thecoupling energy are taken together. Moreover, Eq. (70)is also helpful for setting up parameters for numericalsimulations for which the maximal number of phononsrepresents the bottleneck for the efficiency. C. Crossover from weak to strong coupling in theanti-adiabatic regime
We complete our investigation by studying thecrossover from weak to strong electron-phonon couplingat ω /t = 10. The evolution of the phonon energy E ph ( t ) for weak electron-phonon coupling has alreadybeen analyzed in Figs. 6 and 10(a). The response isgoverned by coherent oscillations without any significantenergy redistribution. This is corroborated by showing E kin ( t ) in the weak-coupling regime in Fig. 16 (see thecurves at g = γ/ω = 0 .
01 and 0 . E kin ( t ) is only atiny fraction of the electronic bandwidth. In contrastto the case ω = t studied in Sec. V B, therefore, theweak-coupling anti-adiabatic regime is characterized bynegligible energy transfer between different parts of theHamiltonian.Nevertheless, weak- and strong-coupling anti-adiabaticregimes are different with respect to the temporal evo-lution since in the latter case the amplitude of oscil-lations becomes very large. We show in Fig. 16 how E kin ( t ) changes when g is varied from small ( g (cid:28)
1) tolarge values ( g > E kin = 2 t and E kin ≈
0. This range of oscillationsis consistent with the result from perturbation theory,Eq. (68), and exhibits the revivals of the initial state atthe integer multiples of the phonon period 2 π/ω . Thesetwo regimes of oscillations are connected by a crossoverregime where the amplitude of oscillations decreases (seethe g = 0 . t . However, there is noindication of the amplitude of oscillations decay to zerofor large L based on the available system sizes. VI. OPTIMAL PHONON MODES ANDSINGLE-SITE ENTANGLEMENT ENTROPYA. Optimal phonon modes
As mentioned in the introductory section, the numberof phonon states needed to correctly describe an electron-phonon coupled system can be quite large. For the single-electron problem studied in this paper, we have shownthat the limited functional space basis works very effi-ciently. However, it is typically restricted to systems that t. ω E k i n ( t ) /t g = 0.01g = 0.2 g = 0.5g = 5 FIG. 16. (Color online) The electron kinetic energy E kin ( t ) at ω /t = 10 for several values of the electron-phonon coupling g = γ/ω . We use LFS and set L = 12. contain only a single (or a few) charge carriers. Zhang et al. first introduced the optimal mode basis for exactdiagonalization [50] with the goal to obtain an efficienttruncation scheme for the on-site Hilbert space of bosonsthat is able to represent the target state with fewer basisstates than in the occupation number basis.A physical example for such an optimal basis canbe given for the strong-coupling, anti-adiabatic limit γ, ω (cid:29) t . In this case, the Hamiltonian is well ap-proximated by Eq. (9) and only two local phonon statesare needed to construct the ground state of Eq. (14): thefirst one is the phonon vacuum state if the site is not oc-cupied by the electron, and the second one is a coherentstate of phonons if the site is occupied by the electron.Thus two phonon states (for each site) build an optimalbasis for the ground state. Any coherent state | β (cid:105) with b j | β (cid:105) = β | β (cid:105) (for any site j ) can be expressed in the localphonon occupation number basis | β (cid:105) = (cid:80) ∞ n =0 (cid:104) n | β (cid:105) | n (cid:105) .The weight of each occupation number state is given bya Poisson distribution P ( n ; | β | ) = |(cid:104) n | β (cid:105)| = | β | n e −| β | n ! (71)with n = 0 , , , . . . . The average number of localphonons is equal to the variance of this distribution N ( j )ph = (cid:104) b † j b j (cid:105) = Var[ P ( n ; | β | )] = | β | . (72)For the ground-state Holstein polaron (14) we have β = g and thus N ph = g . In other words, this ground state canbe constructed in a local basis of dimension d = 2 consist-ing of the phonon vacuum and a coherent phonon state,while one needs a much larger number ∝ g of states inthe occupation number basis. Similarly, we can deter-mine the optimal states for the time-dependent solution | ψ ( t ) (cid:105) = e − iS | (cid:101) ψ ( t ) (cid:105) with | (cid:101) ψ ( t ) (cid:105) given by Eq. (58). Wefind that for all times t only two optimal phonon statesare required to describe the state | ψ ( t ) (cid:105) : the phonon vac-uum state and a coherent state with a complex eigenvalue β = g (1 − e − iω t ). This corresponds to a number of barephonon states N ph = 2 g [1 − cos( ω t )] in agreement withthe phonon energy (64).9We next describe how to extract this basis from anarbitrary state. This is closely related to the procedurein DMRG methods [50, 66, 67]. The lattice is split intotwo parts, say S and E , and the reduced density matrix ρ S of a quantum state | ψ (cid:105) for the subsystem S is definedby ρ S = tr E ( | ψ (cid:105) (cid:104) ψ | ) = (cid:88) α w α | α (cid:105) (cid:104) α | . (73)We get an optimal basis by selecting the density-matrixeigenstates | α (cid:105) with the highest eigenvalues w α . If thedistribution of weights w α is sharp enough, the proper-ties of the quantum state | ψ (cid:105) are almost completely de-termined by the highest optimal states only. The optimalphonon basis is a special case where the subsystem S isa single site.Since the Hamiltonian conserves the number of elec-trons, the reduced density matrix ρ S is block-diagonalin the number of electrons N e . In the case of the Hol-stein model with one electron, there are only two distinctblocks that we denote by N e = 0 and N e = 1. Since we re-strict our subsystem S to a single site only, each block hasthe dimension N max + 1 of the phonon occupation num-ber basis (so-called bare basis). Thus the eigenstates ofthe reduced density matrix ρ S can be grouped into twosets of N max + 1 states according to the correspondingelectron occupation number. Thereafter we denote theeigenstates in the block N e = 0 by | α (cid:48) (cid:105) and those in theblock N e = 1 by | α (cid:48)(cid:48) (cid:105) . In an algorithm that uses thisapproach for truncating the local basis, the number N opt of kept states is smaller than the number of bare modesin the basis, i.e, N opt < N max . Here we keep the fulloptimal basis, i.e., N opt = N max .The sum rule for the weights ω α coming from the twoblocks can then be reorganized as1 = N max +1 (cid:88) α =0 w α = N max (cid:88) α (cid:48) =0 w α (cid:48) + N max (cid:88) α (cid:48)(cid:48) =0 w α (cid:48)(cid:48) (74)with the partial sum rules (cid:88) α (cid:48) w α (cid:48) = L − L (75) (cid:88) α (cid:48)(cid:48) w α (cid:48)(cid:48) = 1 L (76)representing the probability of finding no electron or oneelectron on a given site, respectively. B. Von Neumann entropy
A convenient measure of the entanglement between asingle site and the rest of the system is the von Neumannentropy S vN . It is defined as S vN = − N max +1 (cid:88) α =0 w α ln ( w α ) , (77) t.t S v N ( t ) λ = 0.5 λ = 2 λ = 4.5 Time evolutionGS FIG. 17. (Color online) Time evolution of the von Neumannentropy S vN ( t ) for ω = t and different values of the electron-phonon coupling λ = 0 .
5, 2 and 4 .
5. The values at t = 0 aredetermined by Eq. (78). The horizontal lines represent thecorresponding ground-state (GS) values. We use LFS with L = 8. where the w α are the eigenvalues of the single-site re-duced density matrix ρ S , see Eq. (73). In our work weuse the von Neumann entropy to get insight into the in-crease of the number of relevant optimal modes: a largerentanglement entropy indicates that more optimal modesare relevant for the dynamics. In a more general con-text, since the excess energy is an intensive quantity,this classifies our set-up as a so-called local quench prob-lem [123, 124].For our initial state (2) there is no phonon in the sys-tem, hence the reduced density matrix contains only twonon-zero entries, w = L − L and w = L , which corre-spond to the zero-phonon state on a site without andwith the electron, respectively. The initial entropy istherefore S vN = 1 L ln( L ) + L − L ln (cid:18) LL − (cid:19) . (78)The entropy is hence finite because of the block structureof the reduced density matrix. At the same time, thisentropy also represents the minimal value of S vN in adelocalized electronic system. C. Results
We now turn to the discussion of numerical resultsfor the optimal modes and the von Neumann entropy,obtained by time evolving the initial wavefunction (2)within LFS. We first show the results for the time de-pendence of the von Neumann entropy S vN ( t ) in Fig. 17for three different values of the electron-phonon coupling λ = 0 .
5, 2 and 4 .
5. It starts from the minimum valueat t = 0, given by Eq. (78), and it monotonically in-creases until it reaches the stationary regime. The ini-tial slope of S vN ( t ) increases with λ , in apparent correla-tion with the initial slope of E kin ( t ) shown in Fig. 12(c).Interestingly, even though in the stationary regime at λ = 2 the phononic energy E ph ( t ) and the electron-phonon coupling energy E coup ( t ) exhibit coherent oscil-0 α -12 -9 -6 -3 w ~ α λ = 0.5 λ = 2 λ = 4.5Open symbols: GSFilled symbols: t = t max FIG. 18. (Color online) Weights (cid:101) w α of the optimal modesin descending order for ω = t and different values of theelectron-phonon coupling λ = 0 .
5, 2 and 4 . L = 12 for λ = 0 . L = 8 for the latter twosystems). Filled symbols represent weights in the station-ary regime while open symbols represent the correspondingground-state (GS) values. We choose a time t max in the sta-tionary regime at which the phonon energy E ph ( t ) has a localmaximum (c.f. Figs. 12(a) and 14), i.e., t max t = 15 . λ = 0 . t max t = 9 . λ = 2 and 4 . lations [c.f. Figs. 12(a) and 12(b)] with a non-vanishingamplitude in the thermodynamic limit [c.f. Fig. 13(b)],the von Neumann entropy barely exhibits any tempo-ral dependence. Since S vN ( t ) measures the entanglemententropy between a single site and the rest of the system,this result supports the view that the coherent oscilla-tions indeed stem from the single-site dynamics discussedin Sec. IV E. The time evolution at very strong coupling λ = 4 . tt < t .For the remainder of the study, we focus on the sta-tionary regime only and analyze the eigensystem of thesingle-site reduced density matrix ρ S defined in Eq. (73).In particular, we study the eigensystem on the site occu-pied by the electron, i.e., the N e = 1 block of ρ S . Sincethe corresponding eigenvalues w α (cid:48)(cid:48) do not sum up to one[see Eq. (76)], we normalize them by a factor L , whichis the inverse of the probability to find the electron on agiven site. To simplify the presentation of our results inFigs. 18-20 we use the notation (cid:101) w α = L · w α (cid:48)(cid:48) , | α (cid:105) = | α (cid:48)(cid:48) (cid:105) N e =1 . (79)Figure 18 shows the weights (cid:101) w α in decreasing orderfor three different values of the coupling λ = 0 .
5, 2 and4 .
5. For each λ , we choose a time t max in the station-ary regime where the oscillating phonon energy E ph ( t )has a local maximum [c.f. Figs. 12(a) and 14]. We alsocompare the weights to the corresponding ground-statevalues. In general, both in the ground state and in thestationary state, the weights (cid:101) w α decay roughly exponen-tially (in fact, in Fig. 18 that shows that data at t max ,there is a crossover from a slower to a faster decay).This suggests that the concept of optimal modes couldbe used to reduce the computational cost in correlated | < α | n > | t min .t = 18.4t max .t = 15.80 5 10 n α = 0 (b) α = 1(c) α = 2 (d) α = 3w~ = 0.672 w~ = 0.168w~ = 0.108 w~ = 0.036w~ = 0.648 w~ = 0.192w~ = 0.096 w~ = 0.048 FIG. 19. (Color online) Components |(cid:104) n | α (cid:105)| of the four high-est optimal modes | α (cid:105) in the phonon occupation number basis {| n (cid:105)} . Results are shown in the stationary regime using LFSat λ = 0 . ω = t and L = 12 for two different times: at t min (squares) the phonon energy E ph ( t ) has a local minimum,while at t max (circles) it has a local maximum [see the timeevolution of E ph ( t ) in Fig. 12(a)]. (cid:101) w min α and (cid:101) w max α representthe optimal mode weights (79) at the given times. electron-phonon systems in non-equilibrium as alreadydone for ground states and dynamical correlation func-tions [50, 51]. Nevertheless, the decay of the weights (cid:101) w α is much more rapid in the ground state, consistent witha much smaller entanglement entropy. In addition, inthe stationary regime, the weights decay much slower forlarger electron-phonon coupling λ . Since the weights arenormalized, the smaller decay constant for larger λ in-dicates that more eigenstates (optimal modes) need tobe taken into account to get the same truncation error,which in turn would make the numerical computationmore demanding.Finally, we focus on the structure of the optimal modes | α (cid:105) in the stationary regime. We express their compo-nents |(cid:104) n | α (cid:105)| in the local phonon occupation numberbasis {| n (cid:105)} and plot them at two different times t max and t min corresponding to a maximum and a minimum of thephonon energy E ph ( t ), respectively. Results for the fourhighest weighted optimal modes are shown in Fig. 19 for λ = 0 . λ = 4 . t = t max and t = t min look very much alike, which can be re-lated to the fact that the oscillations are very weak in thisregime and even fully vanish in the thermodynamic limit(c.f. Fig. 13). In addition, the modes are clearly peakedand formed only by a few bare modes (i.e., occupationnumber states). This is expected in the weak-couplingregime where only a relatively small amount of phononsis excited. The highest weighted mode | α = 0 (cid:105) , shownin Fig. 19(a), has its peak at zero phonons, but the sec-ond and third highest optimal modes | α = 1 (cid:105) and | α = 2 (cid:105) shown in Figs. 19(b) and 19(c) seem to have swappedrank because their components peak at the positions of1 | < α | n > | t min .t = 12.1t max .t = 9.30 10 20 30 40 n α = 0 (b) α = 1(c) α = 2 (d) α = 3w~ = 0.368 w~ = 0.128w~ = 0.120 w~ = 0.104w~ = 0.360 w~ = 0.160w~ = 0.128 w~ = 0.088 FIG. 20. (Color online) Components |(cid:104) n | α (cid:105)| of the four high-est optimal modes | α (cid:105) in the phonon occupation number basis {| n (cid:105)} . Results are shown in the stationary regime using LFSat λ = 4 . ω = t and L = 8 for two different times: at t min (squares) the phonon energy E ph ( t ) has a local minimum,while at t max (circles) it has a local maximum [see the timeevolution of E ph ( t ) in Fig. 14(b)]. (cid:101) w min α and (cid:101) w max α representthe optimal mode weights (79) at the given times. Dashedlines in (a) and (b) represent Poisson distributions (71). the two- and one-phonon states, respectively.On the other hand, in the strong-coupling case, the op-timal modes at t = t max and t = t min differ significantly.Figure 20 shows results for λ = 4 . g = 3). Themost striking difference can be observed for the highestweighted mode | α = 0 (cid:105) in Fig. 20(a). At t = t min it isessentially equal to the bare mode with zero phonon oc-cupation (i.e., the phonon vacuum), while at t = t max it strongly resembles a Poisson distribution P ( n ; | β | )with a width set by the actual number of phonons in thesystem, i.e., | β | = N ph = E ph ( t max ) /ω [compare thedashed line with circles in Fig. 20(a)]. This actual valueof N ph differs significantly from the value | β | = 4 g inthe limit t = 0 as seen in Fig. 14(b). Nevertheless, thisresult agrees qualitatively with the picture of a coherentphonon state obtained in the t = 0 limit in Sec. VI A.Therefore, a similar dynamics can be observed for the op-timal mode with the largest weight in a system with finiteparameters λ = 4 . ω = t and in the anti-adiabaticstrong-coupling limit.Apart from the highest optimal mode, also the sec-ond highest optimal mode | α = 1 (cid:105) shown in Fig. 20(b)exhibits intriguing properties. At t max , it represents adistribution peaked at some finite phonon number. Weplot the Poisson distribution P ( n ; g = 9) as a dashedline in Fig. 20(b), which corresponds to the ground-statephonon distribution at t = 0. Interestingly, the latterfunction strongly resembles the numerically calculatedoptimal mode. It suggests that at the local maximum of E ph ( t ), there is a coexistence of distinct coherent states.On the other hand, the modes with lower weights shownin Figs. 20(c) and 20(d) apparently do not exhibit any simple structure. VII. SUMMARY, DISCUSSION, CONCLUSION
In summary, we performed a comprehensive analysisof the non-equilibrium dynamics in the Holstein polaronmodel. We considered the initial state where the elec-tron is highly excited while the phonon energy is zero.Two qualitatively different limiting scenarios exist de-pending on whether the phonon energy is much smalleror much larger than the electronic bandwidth: in the firstone, relaxation sets in, i.e., there is a net energy trans-fer from the electron to the phonon system, while in thesecond case, the response of the system is usually gov-erned by persistent coherent oscillations. The analysisof the real-time dynamics in the entire parameter regimewas done using advanced numerical methods: by com-paring with exact diagonalization and analytical resultsin limiting cases, we benchmarked the TEBD methodand diagonalization in LFS. The increase of the phononnumber as a function of time, which typically happens inthe transient dynamics, is the main problem preventingthe simulation of arbitrarily long times for all parameterregimes. Nonetheless, if this initial increase can be cor-rectly captured, then the LFS method is very efficient indescribing the dynamics for very long times on relativelylarge lattices and in different parameter regimes.The dynamics at weak electron-phonon coupling andsmall phonon energy compared to the electronic band-width is very instructive: the electron dissipates its en-ergy in the relaxation regime, and then enters the sta-tionary regime where temporal fluctuations about theaverage value vanish in the thermodynamic limit. Therelaxation is characterized by the redistribution of elec-tronic momenta from the top ( k = π ) towards the bottom( k = 0) of the electronic band. This represents a proto-typical example of dissipative dynamics in a closed quan-tum system. As one of the main advantages, it allowsfor a direct comparison of the dynamics emerging fromunitary time evolution (calculated using numerical algo-rithms) with semi-classical approaches. We calculatedthe relaxation dynamics from the Boltzmann equationand obtained good agreement with the numerical data.The results show that the relaxation of the many-bodysystem can well be described by a linear decrease of theelectronic kinetic energy with time, and the characteristicrelaxation time is determined by the interplay of all threeparameters of the Holstein model. Assuming a constantelectronic density of states, a compact expression for therelaxation time τ follows from the Boltzmann equation, τ ω = (16 /π )( t /γ ). Our numerical results for relax-ation, obtained by using the actual density of states ofa 1D tight-binding system, are interestingly consistentwith this prediction. As a simple consequence, we ob-serve that there is a lower bound for the relaxation timein this parameter regime: it is set by the inverse phononenergy, i.e., τ ω >
1. Here, the relaxation time τ rep-2resents the complete quasi-particle relaxation time. Itdiffers qualitatively from the characteristic time (cid:101) τ of asingle transition, i.e., a process where a single phononis emitted: in the latter case, the characteristic transi-tion time for the constant density of states is given by (cid:101) τ ω = (4 /π )( t ω /γ ).When the electron-phonon coupling or the phonon fre-quency are increased, coherent temporal oscillations areenhanced and it becomes more difficult to disentanglethe relaxation regime from the stationary regime. Wediscussed two well-defined regimes of model parameterswhere persistent coherent oscillations govern the dynam-ics: when ω is much larger than the electronic band-width and when the electron-phonon coupling γ getsmuch larger than the hopping amplitude t . The sim-plest model which captures the essence of both scenariosis the single-site Holstein model (i.e., when t = 0). Itsground state is a shifted harmonic oscillator and ana-lytical expressions for the time evolution of observablescan be obtained for all times. This limiting case unveilstwo important scales for dynamics: the period of oscilla-tions, which equals 2 π/ω , and the amplitude of oscilla-tions, given by γ/ω . For our initial state without anyphonon, the interpretation of the single-site dynamics isparticularly simple: the phonon state of an occupied siteis a coherent state with the time-dependent eigenvalue β = g (1 − e − iω t ). As a consequence, one can set anupper bound for the number of emitted phonons in thetime evolution: it can not exceed four times the valuefound in the ground state. Our numerical results suggestthat this criterion represents a reasonable estimate forthe dynamics at finite t as well.We also calculated the time dependence of the single-site reduced density matrix: it allows for the investigationof the optimal phonon modes and the entanglement en-tropy during the time evolution. We showed that thestructure of the optimal modes carries valuable informa-tion about the dynamics of the system, and hence com-plements the investigation of non-equilibrium dynam-ics based on observables only. In particular, for strongelectron-phonon coupling, the structure of the highest-weighted optimal mode resembles the phonon distribu-tion in the t = 0 limit. This is consistent with thenumerical observation that the single-site coherent oscil-lations, which can be described analytically for t = 0,persist for finite electron-phonon coupling and finite hop-ping t . On the level of observables, the persistent oscil-lations can be seen, e.g., when measuring the phononenergy: our numerical results suggest that the amplitudeof these oscillations does not vanish in the thermody-namic limit. The structure of other optimal modes alsoexhibit interesting, yet less well understood properties.In any case, non-equilibrium optimal modes representa promising concept for characterizing the dynamics aswell as for designing novel numerical algorithms [50] totreat strongly correlated electron-phonon systems out ofequilibrium, and hence deserve more attention in futurestudies. ACKNOWLEDGEMENT
We acknowledge stimulating discussions with JanezBonˇca, Martin Eckstein, Marcin Mierzejewski and Ines de Vega. We thank Volker Meden for comments on theprevious version of a manuscript. C.B., F.D., F.H.-M.,and E.J. acknowledge support from the DFG (DeutscheForschungsgemeinschaft) through grants Nos. JE 261/2-1and HE 5242/3-1 in the Research Unit
Advanced Compu-tational Methods for Strongly Correlated Quantum Sys-tems (FOR 1807). L.V. is supported by the Alexandervon Humboldt Foundation.
Appendix A: Matrix elements for time evolution inthe weak-coupling regime
In Sec. IV C, we obtained analytical expressions forthe time dependence of observables in the limit of weakelectron-phonon interaction. Contributions up to O ( γ )consist of the three terms given by Eq. (29), (32) and (33).Here we provide the corresponding matrix elements forthe individual parts of the Holstein Hamiltonian, Eq. (3),as well as for the momentum distribution function n k .For the kinetic energy E kin ( t ), only the second-orderterms are non-zero, R (1) K,q = 0 (A1) R (2 a ) K = (cid:15) K = − t cos K (A2) R (2 b ) K,q = (cid:15) K − q = − t cos ( K − q ) , (A3)and the same holds for the momentum distribution func-tion n k R (1) K,q = 0 (A4) R (2 a ) K = δ k,K (A5) R (2 b ) K,q = δ k,K − q . (A6)For the phonon energy E ph ( t ), only a single second-orderterm is non-zero R (1) K,q = 0 (A7) R (2 a ) K = 0 (A8) R (2 b ) K,q = ω , (A9)while for the electron-phonon coupling energy E coup ( t ),which is the only off-diagonal operator in the momentum-space basis, there is a non-zero term already in the firstorder, R (1) K,q = − γ √ L (A10) R (2 a ) K = 0 (A11) R (2 b ) K,q = 0 . (A12)We apply these values to study dynamics in Sec. IV Band IV D.3 [1] J. Orenstein, Phys. Today , 44 (2012).[2] Z. G. Sheng, M. Nakamura, W. Koshibae, T. Makino,Y. Tokura, and M. Kawasaki, Nat. Commun. , 4584(2014).[3] S. Dal Conte, L. Vidmar, D. Goleˇz, M. Mierzejew-ski, G. Soavi, S. Peli, F. Banfi, G. Ferrini, R. Comin,B. M. Ludbrook, L. Chauviere, N. D. Zhigadlo,H. Eisaki, M. Greven, S. Lupi, A. Damascelli, D. Brida,M. Capone, J. Bonˇca, G. Cerullo, and C. Giannetti,arXiv:1501.03833 (to appear in Nature Physics) (2015).[4] C. Gadermaier, A. S. Alexandrov, V. V. Kabanov,P. Kusar, T. Mertelj, X. Yao, C. Manzoni, D. Brida,G. Cerullo, and D. Mihailovic, Phys. Rev. Lett. ,257001 (2010).[5] F. Novelli, G. De Filippis, V. Cataudella, M. Espos-ito, I. Vergara, F. Cilento, E. Sindici, A. Amaricci,C. Giannetti, D. Prabhakaran, S. Wall, A. Perucchi,S. Dal Conte, G. Cerullo, M. Capone, A. Mishchenko,M. Gr¨uninger, N. Nagaosa, F. Parmigiani, andD. Fausti, Nat. Commun. , 5112 (2014).[6] N.-H. Ge, C. M. Wong, R. L. Lingle, J. D. McNeill, K. J.Gaffney, and C. B. Harris, Science , 202 (1998).[7] S. Tomimoto, H. Nansei, S. Saito, T. Suemoto,J. Takeda, and S. Kurita, Phys. Rev. Lett. , 417(1998).[8] S. L. Dexheimer, A. D. Van Pelt, J. A. Brozik, and B. I.Swanson, Phys. Rev. Lett. , 4425 (2000).[9] A. Sugita, T. Saito, H. Kano, M. Yamashita, andT. Kobayashi, Phys. Rev. Lett. , 2158 (2001).[10] A. D. Miller, I. Bezel, K. J. Gaffney, S. Garrett-Roe,S. H. Liu, P. Szymanski, and C. B. Harris, Science , 1163 (2002).[11] C. Gahl, U. Bovensiepen, C. Frischkorn, and M. Wolf,Phys. Rev. Lett. , 107402 (2002).[12] N. Dean, J. C. Petersen, D. Fausti, R. I. Tobey,S. Kaiser, L. V. Gasparov, H. Berger, and A. Caval-leri, Phys. Rev. Lett. , 016401 (2011).[13] F. X. Morrissey, J. G. Mance, A. D. Van Pelt, andS. L. Dexheimer, J. Phys.: Condens. Matter , 144204(2013).[14] M. Scherff, J. Hoffmann, B. Meyer, T. Danz, andC. Jooss, New J. Phys. , 103008 (2013).[15] D. G. Ouellette, P. Moetakef, T. A. Cain, J. Y. Zhang,S. Stemmer, D. Emin, and S. J. Allen, Sci. Rep. , 3284(2013).[16] Y. M. Sheu, S. A. Trugman, L. Yan, J. Qi, Q. X. Jia,A. J. Taylor, and R. P. Prasankumar, Phys. Rev. X ,021001 (2014).[17] F. Rossi and T. Kuhn, Rev. Mod. Phys. , 895 (2002).[18] H. Matsueda, S. Sota, T. Tohyama, and S. Maekawa,J. Phys. Soc. Jpn. , 013701 (2012).[19] M. Mitrano, G. Cotugno, S. R. Clark, R. Singla,S. Kaiser, J. St¨ahler, R. Beyer, M. Dressel, L. Bal-dassarre, D. Nicoletti, A. Perucchi, T. Hasegawa,H. Okamoto, D. Jaksch, and A. Cavalleri, Phys. Rev.Lett. , 117801 (2014).[20] K. A. Al-Hassanieh, F. A. Reboredo, A. E. Feiguin,I. Gonz´alez, and E. Dagotto, Phys. Rev. Lett. ,166403 (2008).[21] B. Koopmans, G. Malinowski, F. Dalla Longa,D. Steiauf, M. F¨ahnle, T. Roth, M. Cinchetti, and M. Aeschlimann, Nature Mater. , 259 (2010).[22] K. Carva, M. Battiato, D. Legut, and P. M. Oppeneer,Phys. Rev. B , 184425 (2013).[23] D. Goleˇz, J. Bonˇca, M. Mierzejewski, and L. Vidmar,Phys. Rev. B , 165118 (2014).[24] E. Iyoda and S. Ishihara, Phys. Rev. B , 125126(2014).[25] M. Eckstein and P. Werner, Phys. Rev. Lett. ,076405 (2014).[26] M. Eckstein and P. Werner, arXiv:1410.3956v1 (2014).[27] L. Vidmar, J. Bonˇca, T. Tohyama, and S. Maekawa,Phys. Rev. Lett. , 246404 (2011).[28] D. Emin and C. F. Hart, Phys. Rev. B , 2530 (1987).[29] J. Bonˇca and S. A. Trugman, Phys. Rev. Lett. , 4874(1997).[30] L. Vidmar, J. Bonˇca, M. Mierzejewski, P. Prelovˇsek,and S. A. Trugman, Phys. Rev. B , 134301 (2011).[31] A. K. C. Cheung and M. Berciu, Phys. Rev. B ,035132 (2013).[32] L.-C. Ku and S. A. Trugman, Phys. Rev. B , 014307(2007).[33] H. Fehske, G. Wellein, and A. R. Bishop, Phys. Rev. B , 075104 (2011).[34] D. Goleˇz, J. Bonˇca, L. Vidmar, and S. A. Trugman,Phys. Rev. Lett. , 236402 (2012).[35] G. Li, B. Movaghar, A. Nitzan, and M. A. Ratner, J.Chem. Phys. , 044112 (2013).[36] V. Meden, J. Fricke, C. W¨ohler, and K. Sch¨onhammer,Z. Phys. B , 357 (1996).[37] J. Fricke, V. Meden, C. W¨ohler, and K. Sch¨onhammer,Annals of Physics , 177 (1997).[38] H. Krull, D. Manske, G. S. Uhrig, and A. P. Schnyder,Phys. Rev. B , 014515 (2014).[39] V. V. Kabanov and A. S. Alexandrov, Phys. Rev. B ,174514 (2008).[40] V. V. Baranov and V. V. Kabanov, Phys. Rev. B ,125102 (2014).[41] Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Phys.Rev. B , 045128 (2015).[42] M. Hohenadler, Phys. Rev. B , 064303 (2013).[43] A. F. Kemper, M. Sentef, B. Moritz, C. C. Kao, Z. X.Shen, J. K. Freericks, and T. P. Devereaux, Phys. Rev.B , 235139 (2013).[44] M. Sentef, A. F. Kemper, B. Moritz, J. K. Freericks, Z.-X. Shen, and T. P. Devereaux, Phys. Rev. X , 041033(2013).[45] K. Yonemitsu and N. Maeshima, Phys. Rev. B ,125118 (2009).[46] G. De Filippis, V. Cataudella, E. A. Nowadnick, T. P.Devereaux, A. S. Mishchenko, and N. Nagaosa, Phys.Rev. Lett. , 176402 (2012).[47] J. Kogoj, Z. Lenarˇciˇc, D. Goleˇz, M. Mierzejewski,P. Prelovˇsek, and J. Bonˇca, Phys. Rev. B , 125104(2014).[48] P. Werner and M. Eckstein, Phys. Rev. B , 165108(2013).[49] P. Werner and M. Eckstein, arXiv:1403.7376 (2014).[50] C. Zhang, E. Jeckelmann, and S. R. White, Phys. Rev.Lett. , 2661 (1998).[51] C. Zhang, E. Jeckelmann, and S. R. White, Phys. Rev.B , 14092 (1999). [52] A. Weiße, H. Fehske, G. Wellein, and A. R. Bishop,Phys. Rev. B , 747(R) (2000).[53] E. Jeckelmann, C. Zhang, and S. R. White, Phys. Rev.B , 7950 (1999).[54] S. Ejima and H. Fehske, EPL , 27001 (2009).[55] H. Fehske, G. Wellein, G. Hager, A. Weiße, and A. R.Bishop, Phys. Rev. B , 165115 (2004).[56] M. Tezuka, R. Arita, and H. Aoki, Phys. Rev. Lett. ,226401 (2005).[57] H. Fehske, G. Hager, and E. Jeckelmann, EPL ,57001 (2008).[58] A. Nocera, M. Soltanieh-ha, C. A. Perroni,V. Cataudella, and A. E. Feiguin, Phys. Rev. B , 195134 (2014).[59] J. Bonˇca, S. A. Trugman, and I. Batisti´c, Phys. Rev.B , 1633 (1999).[60] G. Vidal, Phys. Rev. Lett. , 040502 (2004).[61] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854(2008).[62] A. Polkovnikov, K. Sengupta, A. Silva, and M. Ven-galattore, Rev. Mod. Phys , 863 (2011).[63] T. Langen, R. Geiger, and J. Schmiedmayer,arXiv:1408.6377 (2014).[64] J. Eisert, M. Friesdorf, and C. Gogolin, arXiv:1408.5148(2014).[65] S. R. White, Phys. Rev. Lett. , 2863 (1992).[66] U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005).[67] U. Schollw¨ock, Annals of Physics , 96 (2011).[68] B. Friedman, Phys. Rev. B , 6701 (2000).[69] A. Weiße, G. Wellein, and H. Fehske, Density-Matrix Algorithm for Phonon Hilbert Space Reduction inthe Numerical Diagonalization of Quantum Many-BodySystems, in High Performance Computing in Scienceand Engineering ’01 , edited by E. Krause and W. J¨ager,pp. 131 (Springer Berlin Heidelberg, 2002).[70] S. Sayyad and M. Eckstein, arXiv:1410.4298v1 (2014).[71] T. Holstein, Annals of Physics , 325 (1959).[72] M. Bruderer, A. Klein, S. R. Clark, and D. Jaksch,Phys. Rev. A , 011605 (2007).[73] F. Herrera and R. V. Krems, Phys. Rev. A , 051401(2011).[74] F. Herrera, K. W. Madison, R. V. Krems, andM. Berciu, Phys. Rev. Lett. , 223002 (2013).[75] V. M. Stojanovi´c, T. Shi, C. Bruder, and J. I. Cirac,Phys. Rev. Lett. , 250501 (2012).[76] A. Mezzacapo, J. Casanova, L. Lamata, and E. Solano,Phys. Rev. Lett. , 200501 (2012).[77] J. P. Hague and C. MacCormick, Phys. Rev. Lett. ,223001 (2012).[78] F. Mei, V. M. Stojanovi´c, I. Siddiqi, and L. Tian, Phys.Rev. B , 224502 (2013).[79] L. D. Landau, Phys. Z. Sowjetunion , 644 (1933).[80] J. T. Devreese and A. S. Alexandrov, Rep. Prog. Phys. , 066501 (2009).[81] H. Fehske and S. A. Trugman, in Polarons in AdvancedMaterials , Springer Series in Materials Science, Vol. 103(Springer Netherlands, 2007) pp. 393–461.[82] O. S. Bariˇsi´c and S. Bariˇsi´c, Eur. Phys. J. B , 1 (2008).[83] A. Alvermann, H. Fehske, and S. A. Trugman, Phys.Rev. B , 165113 (2010).[84] F. Marsiglio, Physica C , 21 (1995).[85] G. Wellein, H. R¨oder, and H. Fehske, Phys. Rev. B ,9666 (1996).[86] G. Wellein and H. Fehske, Phys. Rev. B , 4513 (1997). [87] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske,Rev. Mod. Phys. , 275 (2006).[88] E. Jeckelmann and S. R. White, Phys. Rev. B , 6376(1998).[89] S. Ciuchi, F. de Pasquale, S. Fratini, and D. Feinberg,Phys. Rev. B , 4494 (1997).[90] N. V. Prokof’ev and B. V. Svistunov, Phys. Rev. Lett. , 2514 (1998).[91] G. L. Goodvin, A. S. Mishchenko, and M. Berciu, Phys.Rev. Lett. , 076403 (2011).[92] P. E. Kornilovitch, Phys. Rev. Lett. , 5382 (1998).[93] J. P. Hague, P. E. Kornilovitch, A. S. Alexandrov, andJ. H. Samson, Phys. Rev. B , 054303 (2006).[94] M. Hohenadler, H. G. Evertz, and W. von der Linden,Phys. Rev. B , 024301 (2004).[95] T. Ohgoe and M. Imada, Phys. Rev. B , 195139(2014).[96] M. Berciu, Phys. Rev. Lett. , 036402 (2006).[97] G. L. Goodvin, M. Berciu, and G. A. Sawatzky, Phys.Rev. B , 245104 (2006).[98] O. S. Bariˇsi´c, Phys. Rev. Lett. , 209701 (2007).[99] M. Berciu, Phys. Rev. Lett. , 209702 (2007).[100] J. Bonˇca, T. Katraˇsnik, and S. A. Trugman, Phys. Rev.Lett. , 3153 (2000).[101] J. Bonˇca, S. Maekawa, and T. Tohyama, Phys. Rev. B , 035121 (2007).[102] J. Bonˇca, S. Maekawa, T. Tohyama, and P. Prelovˇsek,Phys. Rev. B , 054519 (2008).[103] L. Vidmar, J. Bonˇca, S. Maekawa, and T. Tohyama,Phys. Rev. Lett. , 186401 (2009).[104] L. Vidmar, J. Bonˇca, and S. A. Trugman, Phys. Rev.B , 104304 (2010).[105] D. Goleˇz, J. Bonˇca, and L. Vidmar, Phys. Rev. B ,144304 (2012).[106] O. S. Bariˇsi´c, Phys. Rev. B , 144301 (2002).[107] O. S. Bariˇsi´c, Phys. Rev. B , 064302 (2004).[108] Z. Li, D. Baillie, C. Blois, and F. Marsiglio, Phys. Rev.B , 115114 (2010).[109] M. Chakraborty and B. I. Min, Phys. Rev. B , 024302(2013).[110] L.-C. Ku, S. A. Trugman, and J. Bonˇca, Phys. Rev. B , 174306 (2002).[111] V. Cataudella, G. De Filippis, F. Martone, and C. A.Perroni, Phys. Rev. B , 193105 (2004).[112] G. De Filippis, V. Cataudella, A. S. Mishchenko, andN. Nagaosa, Phys. Rev. B , 094302 (2012).[113] T. J. Park and J. C. Light, The Journal of ChemicalPhysics , 5870 (1986).[114] A. Daley, C. Kollath, U. Schollw¨ock, and G. Vidal, J.Stat. Mech. (2004), P04005.[115] S. R. White and A. E. Feiguin, Phys. Rev. Lett. ,076401 (2004).[116] G. Vidal, Phys. Rev. Lett. , 147902 (2003).[117] M. Einhellinger, Transport simulations in nanosystemsand low-dimensional systems , Ph.D. thesis, Leibniz Uni-versit¨at Hannover (2013).[118] J. M. Ziman,
Electrons and phonons: the theory oftransport phenomena in solids (Oxford, ClarendonPress, 1960).[119] A. Gelzinis, D. Abramavicius, and L. Valkunas, Phys.Rev. B , 245430 (2011).[120] A. Dey and S. Yarlagadda, Phys. Rev. B , 064311(2014).[121] A. S. Mishchenko, N. Nagaosa, G. De Filippis, A. de Candia, and V. Cataudella, arXiv:1406.6486v1(2014).[122] I. G. Lang and Y. A. Firsov, Zh. Eksp. Teor. Fiz.43