{'fraction_non_alphanumeric': 0.04148167493688518, 'fraction_numerical': 0.014581701053364673, 'mean_word_length': 4.231102003642987, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonequilibrium dynamics of the Holstein polaron driven by external electric field

This work represents a fundamental study of a Holstein polaron in one dimension driven away from the ground state by a constant electric field. Taking fully into account quantum effects we follow the time-evolution of the system from its ground state as the constant electric field is switched on at t = 0, until it reaches a steady state. At weak electron phonon coupling (EP) the system experiences damped Bloch oscillations (BO) characteristic for noninteracting electron band. An analytic expression of the steady state current is proposed in terms of weak EP coupling and large electric field. For moderate values of EP coupling the oscillations are almost critically damped and the system reaches the steady state after a short time. In the strong coupling limit weakly damped BO, consistent with nearly adiabatic evolution within the polaron band, persist up to extremely large electric fields. A traveling polaron under the influence of the electric field leaves behind a trail of phonon excitations absorbing the excess energy gained from the electric field. The shape of the traveling polaron is investigated in details.


I. INTRODUCTION
Research in the field of non-equilibrium dynamics of complex quantum systems constitutes a formidable theoretical challenge. Many advanced numerical techniques, ranging from exact-diagonalization 1 , expansion using Chebyshev polynomials 2 , time-dependent density matrix renormalization group 3 , to non-equilibrium dynamical mean field techniques 4 have been developed to tackle this complex problem.
More than forty years ago using a path-integral approach Thornber and Feynman 5 discovered that an electron in a parabolic band, driven by the electric field, acquires a constant velocity due to emission of phonons. Later approaches to polaron motion in high electric field used Boltzmann equations 6 , the high-field drift velocity was estimated via phonon-assisted hopping between different rungs of Wannier-Stark states using rate equations 7,8 . In ref. 9 the one-dimensional Holstein polaron problem in strong electric field has been mapped on a nonstandard Bethe lattice. It has been realized that keeping full quantum coherence between many-body states is crucial to obtain finite drift velocity for dispersionless optical phonons. Extensive research of polaron dynamics has been conducted within the semiclassical Su-Shrieffer-Heeger model to describe properties of conjugated polymers that may be used in a variety of applications like molecular electronics or light-emitting diodes [10][11][12][13][14][15] . Polaron formation and its influence on transport properties has also been investigated in the context of DNA molecules within the semiclassical Peyrard-Bishop-Holstein model [16][17][18][19][20][21][22] and other polaron-like models [23][24][25] .
Bloch oscillations (BO) represent a fundamental phenomenon in quantum mechanics where a charged particle in a periodic potential exhibits a periodic motion when exposed to an uniform external electric field. Since the electrons in solids can dissipate energy due to scattering from inelastic degrees of freedom on a time scale usually shorter than a typical Bloch time t B , it took a long time until the first experimental observation of BO was carried out on semiconducting superlattices [26][27][28][29][30] and later in optical potentials [31][32][33] . Nowadays the concept of BO is frequently present in a variety of different fields, for instance atomic Bose-Einstein condensates in optical lattices [34][35][36][37] , interacting quantum few-body systems [38][39][40] or organic molecules 14,21,22,41 . However, the description of damping of BO in dissipative medium remains a challenging task.
By choosing the Holstein Hamiltonian as one of the simplest model systems describing the interaction between a fermion and phonons, we are able to investigate the field-induced acceleration of the polaron, which simultaneously dissipates the energy by inelastic scattering on optical phonons while maintaining the full quantum nature of the problem. Following the time evolution of the ground state when the electric field is switched on at time t = 0, we show how the polaron reaches the steady state and consequently develops a constant non-zero velocity. In particular, we calculate the steady-state current vs. voltage characteristics of the Holstein polaron for different regimes of electron-phonon couplings.
We discuss the Holstein model in one spacial dimension and give a brief overview of the numerical method in the second paragraph. In the third we discuss numerical results.
Here we give special emphasis on the time evolution from the ground state towards the steady state by presenting various correlation functions in different EP coupling regimes. We compare our results with a simple Landau Zener model and follow the time evolution of the polaron as it starts propagating after switching on the electric field. As a focal point of this work we discuss the dependence of the steady state current on the EP coupling and electric field. In the last paragraph we give conclusions.

II. MODEL AND NUMERICAL METHOD
We analyze the one-dimensional Holstein model with a single electron, threaded by a time-dependent flux: where c † i and a † i are electron and phonon creation operators at site i, respectively, and n i = c † i c i is electron density. ω 0 denotes a dispersionless optical phonon frequency and t 0 nearest-neighbor hopping amplitude. The dimensionless EP coupling strength is λ = g 2 /2t 0 ω 0 . The constant electric field F that is switched on at time t = 0 enters the Hamiltonian in Eq. 1 through the time-dependent phase θ(t) = −F t for t ≥ 0. We measure electric field F in units of [t 0 /e 0 a] where e 0 is the unit charge and a is the lattice distance. We furthermore measure time in units of [ /t 0 ]. Unless otherwise specified, from here on we set a = e 0 = = t 0 = 1.
To solve the time-dependent Hamiltonian for a single electron coupled to phonon degrees of freedom we use an improved numerical method, originally introduced in Ref. 42 , that led to numerically exact solutions of the polaron ground and low-lying excited state properties. The method constructs the variational Hilbert space (VHS) starting from the singleelectron Bloch state c † k |∅ with no phonons on an infinite lattice. The VHS is then generated by applying the off-diagonal terms of Hamiltonian (1) where H kin and H g correspond to the first and the second term of the Hamiltonian in Eq. 1, respectively. Parameters N h and M determine the size of the VHS. In addition, N h − 1 represents the maximum distance between the electron and the phonon quanta and N h * M is the maximum number of phonon quanta contained in the Hilbert space. The parameter M > 1 (Ref. 43 ) ensures good convergence in the strong EP coupling regime that contains multiple phonon excitations.
To reach weak coupling regime, λ << 1, we introduce an additional parameter N phmax limiting the maximum number of phonon quanta, which enables construction of VHS as large as N h = 40.
We first solve the Hamiltonian in Eq. 1 for F = 0, i.e., we calculate the (zero temperature) polaron ground state 42,44-53 . Then we switch on the uniform electric field and start the time propagation from the initial state using the time-dependent Lanczos technique 54 . We manage to find numerically accurate results of the model away from equilibrium while maintaining full quantum description of phonons. Since we are dealing with a single particle in an infinite system, we compute time-dependent average of the current operator j(t) = Î (t) , whereÎ In the case of a time-independent field F , the time-integral of the current is directly related to a change of the total energy where ∆h(t) = H(t) − H(t = 0) and x(t) represents the travelled distance 55 .

A. Time evolution from the ground state towards the steady state
We first present results obtained near the noninteracting limit, i.e. at λ = 0.01, where j(t) displays damped BO around j(t → ∞) > 0, see Fig. 1(a). The period and the initial amplitude of BO at small t are consistent with BO of a free electron, denoted with thin dashed line in Fig. 1(a). Damping is due to inelastic scattering on phonons that is in turn reflected in a monotonic increase of the average phonon number n ph with time, as depicted in Fig. 1(c). Damping is, however, not the most important consequence of inelastic scattering. Notably, j(t) approaches a positive steady state currentj for t > t s ∼ 4t B and 6t B at F = 1/5 and F = 1/2, respectively, where t B denotes the Bloch oscillation period t B = 2π/F . Note thatj << j max = 2. The dependence of the steady state currentj on F will be discussed further in the text. The steady state current as well emerges as a linear dependence of the total energy on time: ∆h(t) = Fjt + ∆h 0 , see results in Fig. 1(e), where with increasing t, ∆h(t) approaches a straight line. In the steady state we as well observe a linear increase of n ph vs. t. When comparing ∆h(t) and n ph in the linear regime we find that ∆ḣ(t) = ω 0 d n ph /dt. This equality confirms an intuitive expectation that in the steady state the total energy gain is entirely absorbed by the lattice.
On a more technical side we note that to reach a steady state, the Hilbert space used in our calculation must contain large enough set of excited states that in turn represent the reservoir for the absorption of energy. For this reason, different Hilbert spaces were used, depending on the strength of EP coupling and the size of F , see as well Caption of Fig. 1.
At a larger value of EP coupling, λ = 0.2, a somewhat different physical picture emerges, as shown in Figs. 1(b,d) and (f). The main differences can be summarized as: (i) BO become overdamped, (ii) j(t) remains positive at all t, and (iii) j(t) reaches a steady state after a short time t s t B . Characteristic for a steady state are linear t dependencies of n ph and ∆h(t) in Figs. 1(d) and (f), respectively. Common to all cases presented in Fig. 1, is the emergence of a constant steady-state current for t > t s .
In Fig. 2  wavevectors k = F t = 2πt/t B . The response of the system to external field is nearly elastic, since ∆h(t = l * t B ) ∼ 0 for any integer value of l. The average current remains indistinguishable from zero in the largest time interval tested with our calculation, i.e. t ≤ 20t B . In order to illuminate this behavior we note that in the strong coupling limit a large gap ∆ exists in the polaron excitation spectrum being of the order of ω 0 . The low-energy polaron excitation spectrum is presented for ω 0 = 1 and λ = 2 in the inset of Fig. 2(a) where a gap ∆ ∼ 0.64 separates the polaron band from the excited polaron band 42,56-58 , located just below the continuum denoted by the grey area. At small F << ∆ there exist exponentially small probability for a nonadiabatic transition from the polaron band to the excited polaron band or/and into the continuum.
For sufficiently large F (see Figs. 2(b,d,f) for F = 1/3 and 1), BO in j(t) lose periodicity even though remnants of BO remain clearly visible, and the time averaged current becomes finite (nonzero). Additional frequencies appear in n ph that indicate multiple phonon excitations due to polaron transitions to excited polaron bands. Moreover, the average value of n ph between successive t B intervals increases. The total energy ∆h(t) as well increases in time. At large field,  The observed behavior in the strong coupling regime due to a large gap in the spectrum resembles Landau-Zener (LZ) transition 59,60 , where the probability for tunneling between bands in a two level system, is given by where ∆ is the energy gap between the two levels and v is the velocity. Using Eq. 6 we estimate the threshold electric field F th using Such an estimate has been used to determine the dielectric breakdown of the insulating half-filled Hubbard model 1,61 . Applying Eq. 7 to the specific case of λ = 2 presented in Fig. 2, using ∆ ∼ 0.64 and v = j max ∼ 0.1, we obtain F th ∼ 1.0. LZ formalism gives roughly the correct order of magnitude of F th , since a noticeable current appears around F = 1/3, as seen in Fig. 2(b) and (f) as well as from steady state current, presented in Fig. 4(c). One should however be mindful when considering the LZ formalism. In the polaron case the band structure deviates significantly from the ideal LZ model with two hyperbolic bands. In the realistic case, multiple transitions occur from the polaron band to a continuum of excited states composed from a polaron and additional phonon degrees of freedom.

C. Time evolution of polaron
In Fig. 3 we follow the time evolution of the polaron towards the steady state at F = 1/2. We compute the average number of phonon quanta located at a given distance r from the electron fulfilling the following sum-rule n ph = r γ(r). At t = 0, γ(r) displays a pronounced peak at the position of the electron, i.e. at r = 0, consistent with the shape of the polaron in its k = 0 ground state. After the electric field is switched on, γ(r) experiences a compelling time evolution with three outstanding characteristics: (i) the overall increase of γ(r) with time, (ii) development of pronounced asymmetry of γ(r) with respect to the electron position at r = 0, and (iii) increased amount of polaron excitations in the forward direction. The overall increase of γ(r) is consistent with the absorption of energy that is deposited in increasing number of phonon excitations. The asymmetry is a result of a growing phonon tail, extending behind the moving polaron. Note that the polaron is moving from left to right. In the long time limit, γ is expected to be approximately constant, independent of r and t, for r sufficiently negative. The average height of the polaron tailγ is due to energy conservation requirement independent of λ: thereforeγ ∼ F/ω 0 , compare Figs. 3(d) and (h). Note that this relation holds only when the system has reached a steady state. The length of the polaron tail is given by the expression for the travelled distance x(t) = ∆h(t)/F . At t = 8t B and λ = 0.01 we obtain from Fig. 1(e) x t = x(8t B ) ∼ 14.8, that fits well with the length of the phonon tail in Fig. 3(d).
Rather unexpected is the pronounced increase of phonon excitations in the forward direction where γ(r) 0 up to r ≤ r f ∼ 5 − 7 for all t > 0 presented in Fig. 3. Since time evolution starts from the ground state at zero temperature, there are no phonon excitations present far ahead from the moving electron. A substantial forward tail of phonon excitations is a consequence of damped BO. Indeed, r f compares well with the Stark localization length, i.e. r f ∼ L S = 4/F = 8. Yet another intriguing feature in γ(r) emerges as regular oscillations in the polaron tail with a period K = ω 0 /F = 2, clearly seen in the small λ = 0.01 limit, see Figs. 3(b) through (d). At larger λ = 0.2 these oscillations become overdamped.

D. Steady state current
The focal point of this work is the calculation of the steady state currentj and analysis of its dependence on F and λ. In Fig. 4 we present the current-voltage characteristics, i.e.j vs. F for different values of λ. Note that the upper limit of j is given by the current amplitude j max = 2 in the noninteracting system. We have limited our calculations to commensurate values F = ω 0 /K with K being integer, with two exceptions: (i) large F > ω 0 where we have chosen F = 2ω 0 , 3ω 0 , . . . and (ii) results presented with disconnected triangles in Fig. 4(a), with details given in the figure caption. In the regime λ ≤ 0.1, presented in Fig. 4(a),j decreases with increasing F for F 0.1. Our method does not yield steady state results in the regime F 0.1 due to large Stark localization length L S = 4t 0 /F . Sincej = 0 for F = 0 as well as in the opposite limit, when F → ∞, there must exist a global maximum valuej max that depends on λ. For λ = 0.1,j max ∼ 0.82, while for λ < 0.1,j max is reached somewhere in the interval 0 < F < 0.1, not accessible by the present numerical method. Choosing rational or even irrational values of ω 0 /F leads to a decrease ofj that nevertheless remains non-zero even in the latter case. A sweep over continuous values of F would lead to spikes inj located at integer values of ω 0 /F , as consistent with observations in previous works [7][8][9]62 .
To gain further insight into the decrease ofj with F , we plot in Fig. 4(b)j/ √ λ vs. 1/ √ F and realize that curves approximately collapse onto a straight line. The revealed scaling with 1/ √ F is a clear signature that we are dealing with a coherent propagation between Stark states with identical total energy that are spaced by K = ω 0 /F . This is in contrast with the assumption of an incoherent hopping between localized states 7,8 , which would predict a dependence j ∝ 1/F . In turn our derivation, as presented in the Appendix, leads for integer K > 1 as well as for ω 0 < W = 4t 0 to a scaling of the maximum steady current While j 0 cannot be directly compared to the average current j, the functional dependence on λ and F is in good agreement with scaling in Fig. 4(b) that leads to α ∼ 0.89 (fit is represented by a dashed line). The expression in Eq. 10 is valid in the small λ and large F regime (however, F < ω 0 ) wherej decreases due to decreasing overlap between Stark states. The scaling breaks down when with decreasing F ,j approaches the maximumj max . In Fig. 4(c) we present results for larger λ ∈ [0.2, 2.0] which enables us to observe the evolution ofj vs. F as the system evolves from the weak EP coupling (λ < 1) towards the strong EP coupling (λ > 1) regime. With increasing λ, the position ofj max shifts towards larger values of F while it decreases in its magnitude. The main difference between the weak and strong EP coupling regime emerges due to increasing energy gap ∆ in the polaron excitation spectrum 56-58 that for λ >> 1 approaches ∆ ∼ ω 0 . Due to large ∆ at large λ > 1,j remains zero until F ∼ F th .
We summarize the overview of numerical results with a diagram describing different regimes characterized by distinct short-time behaviors (after switching on F ), presented in Fig. 4(d). We distinguish four different regimes: (i) The regime of damped free particle BO (DBO) for small values of g. (ii) Almost critically damped (CD) regime where steady state current is reached in a time shorter than or of the order of t B , and the oscillations in the current are still visible, however j(t) > 0 for any t > 0. (iii) Polaron BO regime (PBO) where system evolves nearly adiabatically. Polaron Bloch oscillates within the polaron band and damping is exponentially small (numerically undetectable). In PBO average current remains zero and total energy remains periodic within numerical accuracy and up to the largest measured time t ≤ 20t B . (iv) Damped polaron BO regime (DPBO) where remnants of PBO are seen in j(t) while there exists a measurable average currentj > 0 within t ≤ 20t B .

IV. CONCLUSIONS
In summary we list our main results. Using a time dependent Lanczos method we have followed the time evolution of the polaron from its ground state towards the steady state after the electric field has been switched on. Different sizes of VHS have been used to ensure that presented results are valid in the thermodynamic limit. Steady state conditions have been reached at intermediate to high electric fields and the current vs. voltage characteristics has been plotted for different regimes of EP couplings. By calculating the electron-phonon correlation function representing the time evolution of the polaron, we show that the absorbed energy in the steady state is deposited as an increasing number of phonon excitations arranged as a growing tail behind the moving polaron.
The damped BO can be observed in the extremely weak EP coupling limit. In the former case, period of BO t B = 2π/F should be less than the relaxation time t 0 /g 2 related with the emission of phonons. A large gap in the spectrum in the strong coupling regime is responsible for observation of nearly perfect BO arising from the polaron motion within the polaron band. The breakdown of this quasiadiabatic regime qualitatively resembles the Landau-Zener transition from the polaron band to higher excited states. Analytical estimate for the steady state current on the electric field and EP coupling constant at large fields is proposed and numerically tested. The unusual steady state current vs. electric field dependence, j ∝ λ/F , valid at large F and small λ, reflects the significance of coherent processes for proper description of polaron motion. In contrast, approaches calculating the steady state current relying on probabilities for transitions between neighboring Wannier-Stark states, mediated by the EP coupling, yieldj ∝ λ/F 7,8 .
connects the average displacement in the Stark basis ∆ǫ l = F ∆l to phonon generation ω 0 ∆N ph . In the following we consider only the simple commensurate case where ω 0 /F = K is integer (K > 1, i.e. F < ω 0 ), where electrons can perform coherent hopping between Stark states with ∆l = K keeping E 0 constant by emitting (or absorbing) single phonon via the coupling term H ′ = g i n i (a † i + a i ). Also restricting phonon frequencies to ω 0 < W = 4t 0 we remain in the regime K < L S . We now construct the basis of possible coherent states having the same E 0 ∼ 0 starting with a bare electron state at l = 0 and generating novel states by application of H ′ whereby j m denote location of phonons. The matrix elements between subsequent states can be evaluated explicitly by neglecting multiple occupations of sites (being rare for K ≫ 1), i.e. j 1 = j 2 . . . or equivalently simplifying boson factor for multiply occupied sites, i.e., which leads to group velocities in the tight-binding form v p ∝ v p = v 0 sin(pK) with the maximum The derivation can be made more rigorous taking into account the actual Stark wavefunctions w j and matrix elements Eq.(A.5). Still it is not expected to change qualitatively the scaling of coherent group velocities v p with the maximum v 0 , Eq.(A.11). It should be, however, reminded that we did not yet match the actual solution Eq.(A.10) with the boundary condition, as determined, e.g., with the first equation in the system Eq. (A.7). Anyhow, it is expected that an eigenstate of stationary Hamiltonian, as in Eq.(A.1), cannot posess a finite steady current (solution being a superposition of ±p eigenstates). On the other hand, the driven system and the time-dependent model, Eq.(1), clearly can generate the current j(t) and in this sense induce solutions with the steady currentj ∝ j 0 = v 0 following Eq.(A.11). Evidently, more rigorous relation between the eigenstates of the stationary case and the driven problem is still desired.