Dissociation of a Hubbard--Holstein bipolaron driven away from equilibrium by a constant electric field
aa r X i v : . [ c ond - m a t . s t r- e l ] M a r Dissociation of a Hubbard–Holstein bipolaron driven away from equilibrium by aconstant electric field
D. Goleˇz, J. Bonˇca,
1, 2 and L. Vidmar J. Stefan Institute, SI-1000 Ljubljana, Slovenia Department of Physics, FMF, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia
Using a variational numerical method we compute the time-evolution of the Hubbard-Holsteinbipolaron from its ground state when at t = 0 the constant electric field is switched on. The systemis evolved taking into account full quantum effects until it reaches a quasistationary state. In thezero-field limit the current shows Bloch oscillations characteristic for the adiabatic regime wherethe electric field causes the bipolaron to evolve along the quasiparticle band. Bipolaron remainsbound and the net current remains zero in this regime. At larger electric fields the system enters thedissipative regime with a finite quasistationary current. Concomitantly, the bipolaron dissociatesinto two separate polarons. By examining different parameter regimes we show that the appearanceof a finite quasistationary current is inevitably followed by the dissociation of the bipolaron. PACS numbers: 71.38.Mx,63.20.kd,72.10.Di,72.20.Ht
I. INTRODUCTION
Nonequilibrium phenomena in interacting many–bodysystems have recently drawn significant attention. One ofthe most intriguing challenges is to understand how a sys-tem responds to an external electric field and to unravelthe microscopic mechanism that gives rise to a constantcurrent and a steady increase of the total energy of thesystem. From a theoretical point of view, a non–linearcurrent–voltage characteristics with a threshold behav-ior is well–established, and recently the appearance ofa regime of negative differential resistivity has been ob-served in several interacting systems.
Moreover, theonset of Bloch oscillations (well–known for the nonin-teracting systems) has been observed in different corre-lated systems at very large electric fields and evenin the integrable models. Oscillations in the current re-sponse can also exhibit a more complex pattern leadingto beats in a certain regime of model parameters.
Nev-ertheless, the heating arising from the energy flow to thesystem represents a serious obstacle in calculation of aconstant nonzero current for systems at half–filling, therefore most of the calculations investigating current–voltage characteristics have been limited to one dimen-sion.
One of the widely used concepts is to couple thesystem to a thermal bath, which at a loss of quantum co-herence enables calculations of a quasistationary current.Recently, the influence of coupling to a thermal bath onthe current–voltage characteristics of the Hubbard modelhas been investigated. The reported influence is similarto the role of quantum phonons on dissipation of theexcess energy and consequently on the current–voltagecharacteristics of the Holstein model. A complementary approach to investigation of quan-tum systems out of equilibrium is to study dynamics of asingle carrier propagating in dissipative medium, inparticular a charge carrier driven by a constant electric–field. Such problems have been recently solved for adriven carrier doped into Mott insulator and a driven carrier coupled to Holstein phonons, where the nonzerocurrent arises due to constant emission of magnons andphonons, respectively. Recently, using a state–of–the–arttruncated Lanczos method the quasistationary cur-rent in a doped 2D Mott insulator coupled to phononswas calculated. The advantage of such approaches isreflected in calculation of a stable well-defined currentwhile preserving the full quantum nature of the prob-lem. In addition, the total energy gained by hopping of acharge carrier along the direction of the electric field is en-tirely absorbed within the microscopic model. This ideaenables the formulation of the time–evolution of differentsubsystems without coupling to an external reservoir.The aim of our study is to extend a recent study of adriven charge carrier coupled to Holstein phonons to atwo–particle problem as described within the Hubbard–Holstein Hamiltonian. There are two sources of theparticle-particle interaction in this problem: a) an indi-rect interaction, mediated by the electron–phonon cou-pling and b) the direct on–site Coulomb interaction. Thismotivates us to address the following questions: (i) Is itpossible to drive a bound pair of electrons out of equilib-rium in such way that a finite current is obtained withoutbreaking the pair? We show that a Hubbard–Holsteinbipolaron at initial time t = 0, driven by a constantelectric field, begins to dissociate as soon as a nonzeroquasistationary current is reached. (ii) How is the qua-sistationary current calculated per particle changed in atwo–electron system with respect to the single–electroncase? (iii)
We analyze beats in the transient current re-sponse. We show that the beats in our model emergeif the energy scale associated with the Bloch oscillationswithin the quasiparticle (QP) band coincides with theenergy gap between the QP band and the low–energycontinuum of excited states.The paper is organized as follows. We introduce themodel and numerical method in Sec. II. In Sec. III weshow numerical results for the short–time behavior of adriven Hubbard–Holstein bipolaron with a focus on thedissociation of a bound state due to the constant externalelectric field. In Sec. IV we address properties of thecurrent–field characteristics in the quasistationary statewhile in Sec. V we analyze the Fourier spectrum of thereal–time current. We give conclusions in Sec. VI.
II. MODEL AND NUMERICAL METHOD
We define the one-dimensional time-dependentHubbard–Holstein Hamiltonian, threaded by an externalflux H = − t X l,σ h e iφ ( t ) ˜ c † l,σ ˜ c l +1 ,σ + H .c. i + g X j n j ( a † j + a j )+ ω X j a † j a j + U X j n j ↑ n j ↓ , (1)where c † j and a † j are local electron and phonon creationoperators on site j , respectively, and n j = P σ c † j,σ c jσ is electron density. ω denotes the dispersionless opti-cal phonon frequency, t is the nearest–neighbor hoppingamplitude and U represent Hubbard repulsion betweenon-site electrons. We set ω /t = 1 . λ = g / t ω . The electric field isapplied through a time–dependent Aharonov-Bohm flux φ ( t ) = − F t , which generates force on the charged parti-cle due to the Faraday’s law. The constant electric field F is switched on at time t = 0 and is measured in unitsof [ t /e a ], where e is the unit charge and a is latticeconstant. The time is furthermore measured in units of[ t / ~ ], and we set t = ~ = e = a = 1.In our calculations we obtain reliable numerical resultsfor electric fields in the range between 10 − V /cm ,if the model parameters are set to t = 0 . eV and a = 10 − m . These fields become relevant for compar-ison with experiments as presented in Ref. if the ex-trapolation to zero temperature T → In the case of the Hubbard–Holsteinmodel the method generates the VHS starting with thetranslationally invariant initial state | ϕ i where bothelectrons are located on the same site with no phonondegrees of freedom. The VHS is then generated by re-peatedly applying the off-diagonal terms of Hamiltonianin Eq. (1), n | ϕ N h ,Ml i o = [ H kin + M X m =1 ( H g ) m ] N h | ϕ i , (2)where H kin and H g are the first and the second term ofthe Hamiltonian in Eq. (1). Parameters N h and M deter-mine the size of the VHS. The parameter M > λ <
1, weintroduced an additional parameter N phmax which limitsthe maximal number of phonon quanta, enabling largervalues of N h .We first calculate the ground state at F = 0 of theHubbard–Holstein Hamiltonian in Eq. (1). Equilibriumproperties of bipolarons have been extensively studiedin the literature, with particular emphasis on theconditions for formation of a bound bipolaron state. Weswitch on the uniform electric field F at t = 0 and startthe time propagation using the time-dependent Lanczostechnique. The time step is taken to be small enough( dt/t B = 0 . × − where t B = 2 π/F denotes the Blochoscillation period) to ensure good numerical accuracy.We investigate only the S ztot = 0 subspace. An importantobservable is the time–dependent average of the currentoperator per particle, j ( t ) = h ˆ I ( t ) i /
2, whereˆ I ( t ) = i X l,σ e − iF t c † l,σ c l +1 ,σ − H . c . . (3)Since we are dealing with the time–independent field F ,the time integral of the current is directly related to achange of the total energy Z t j ( t ′ ) dt ′ = ∆ h ( t ) / F, (4)where ∆ h ( t ) = h H ( t ) i − h H ( t = 0) i . We also calcu-late the average number of phonon quanta at time t , h n ph i ( t ) = h P j a † j a j i .Due to the finite size of the VHS we are able torun the propagation only until the system reachesthe boundary of the variational space. There are twotypical boundaries the system can reach: (i) the sizeof the bipolaron becomes larger than N h , or (ii) thetotal number of phonons in the system approachesthe maximal allowed value N phmax . We checked thefinite size effects by comparing the expectation valuesof different observables for different system’s sizes. InFig. 1 we represent the time dependence of the relativedifference between energies for different systems sizes,namely | E i − E N h =15 | /E N h =15 , where i represents theparameter N h of the functional generator in Eq. (2)and the maximum number of phonons were limited to N phmax = 15. Our simulations were stopped when therelative energy difference between largest and secondlargest system was more that 5%. Using the sum rule,Eq. (4), we were able to independently control theprecision of the time-evolution, where errors occur dueto time discretization. III. REAL–TIME BEHAVIOR
Previous studies have shown the importance ofthe low-energy spectrum for the short–time behavior of | d E i | / E t/t b i=13i=12i=11 FIG. 1: (Color online) Time dependence of the relative dif-ference between energies for different systems sizes | E N h − E | /E , where E N h is energy of the system generated by thefunctional generator, see Eq. (2), using N h as parameter andmaximal number of phonons were limited to N phmax = 15.We set λ = 0 . F = 1 . the interacting systems after switching on the static elec-tric field. In terms of short–time dynamics there existat least three distinct energy scales. First is the bind-ing energy of the bipolaron δ representing the energy re-quired to dissociate a bound bipolaron into two separatepolarons, δ = E bi − E pol , where E bi and E pol are theenergies of a bipolaron and a polaron, respectively. Thesecond is the phonon frequency ω , the third is the gapin the excitation spectrum ∆ separating the top of thebipolaron band from either the continuum of states or theexcited bipolaron band. Due to the existence of a finitegap, there exists a threshold field F th separating the adi-abatic regime for F < F th and a dissipative regime witha finite current at F > F th , see also Refs. . Character-istic features of the adiabatic regime are periodic Blochoscillations as the bipolaron moves along the bipolaronband. As a consequence, a time–averaged current in theadiabatic regime is zero. On the other hand, in the dis-sipative regime a quasistationary state is reached aftera short transient time t > t tr . A quasistationary stateis characterized by a constant time–independent current j ( t ) = ¯ and a linear growth of the total energy of thesystem, ∆ ˙ h ( t ) = ¯ F. (5)In the case of a driven polaron it has been shown that thesteady increase of ∆ h ( t ) is entirely due to the increaseof phonon excitations left in the wake of the travelingpolaron. A reasonable conjecture in the case of the bipolaronis that for a large binding energy, i.e. , for δ > ω , thesystem on short time scale propagates as a bound bipo-laron, depositing the excess energy in the form of excitedphonons left in its wake. Our calculations show that inthe weak coupling regime ( λ ≪
1) the bipolaron dissoci-ates on a time scale much shorter than t B . We therefore focus in the following on the intermediate and strongcoupling regimes ( λ = 0 . .
9, respectively), wherethe propagation of a bound pair under the electric field ismore likely. However, as we shall show, such propagationis not observed in our calculations, instead, the bipolaronbegins to dissociate as soon as the system enters the dis-sipative regime with a finite electric current.In this section we focus on U = 0. The motivation forchoosing λ = 0 . λ = 0 . ω ( δ ≈ . < ω ) for λ = 0 . λ = 0 . ω ( δ ≈ . > ω ).We also note that in the case of λ = 0 . ∼ .
85 above the groundstate at k = 0. This state represents a bound state ofthe bipolaron with excited phonon cloud. Similar boundexcited states have been well-studied for the case of theHolstein polaron. They form dispersive bands justbelow the one-phonon continuum that starts at ω abovethe ground state energy.For λ = 0 . F = 1 / j ( t ) as shown in Fig. 2(a), which are consistentwith the bipolaron ground state dispersion. Weak damp-ing due to inelastic scattering on phonons is observed,that is in turn reflected in slow increase in the total sys-tem energy ∆ h ( t ) in Fig. 2(c). For larger electric field F = 1 / t >
0. From the available data we estimate the thresh-old field F th ∼ .
2. Although for λ = 0 . h ( t ) shows increase in time which is approximately lin-ear. The increase of ∆ h ( t ) allows us to use a linear fitaccording to Eq. (5) and to calculate the value of qua-sistationary current ¯ . The corresponding current–fieldcharacteristics and the influence of model parameters on¯ will be discussed in Sec. IV.We also observe a quasi–linear increase of the totalphonon quanta in the system h n ph i ( t ), see Fig. 2(b), andwe find that ∆ ˙ h ( t ) > ∼ ω d h n ph i /dt . While the equalitywould mean that the entire energy absorbed from theelectric field is absorbed by the lattice, the inequalitysignifies that there exists excess energy that is responsiblefor the dissociation of the bipolaron. Further insight intothis process can be obtained by calculating the averageparticle distance d ( t ) defined as d ( t ) = sX i,j h n i n i + j i j , (6)shown in Fig. 2(d). We detect an overall increase of d ( t )for both values of F , signaling the dissociation of thebipolaron. We discuss the process of dissociation in moredetail in Sec. III A and Sec. III B.In the strong–coupling regime at λ = 0 . F = 0 . ≈ F th ,shown in Figs. 3(a), 3(c), 3(e) and 3(g). Remarkably, FIG. 2: (Color online) j ( t ), h n ph i ( t ), ∆ h ( t ) and d ( t ) vs t/t B for λ = 0 .
5. We represent characteristic fields in nearly adia-batic regime and when
F > ∼ F th , i.e. , F = 1 / F = 1 / λ = 0 . N h = 17, M = 2, N phmax = 15. beside the previously mentioned effects, we observe beatson a longer time–scale. We discuss their origin in Sec. V.If the values of electric field are increased above F ≈ F th , the beats become less pronounced. For F = 1 and F = 2 we obtain, after initial oscillations, a well definedquasistationary current j ( t ) and a nearly linear time–dependence of h n ph i ( t ) and ∆ h ( t ) characteristic for thedissipative regime, see Figs. 3(b), 3(d) and 3(f). A. Bipolaron dissociation
The central goal of our work is to understand whetherin the long–time limit the bipolaron remains bound (insome parameters regimes) as it travels under the influ-ence of the constant electric field, or it dissociates intotwo separate propagating polarons. The most interestingis the regime of strong EP interaction where the bipo-laron binding energy δ is larger than the phonon energy ω . Such case is shown in Figs. 3(g) and 3(h) for λ = 0 . d ( t ) is steadily increasing even inthe near-adiabatic regime, and the slope of the linear in-crease of d ( t ) is enhanced in the dissipative regime, asseen for F = 1 and 2 in Fig. 3(h).We should stress that the overall increase of d ( t = t max ) after the maximal propagation time relative to theinitial ground state distance d ( t = 0) is more than 6–foldin both cases, i.e. , at F = 1 and F = 2. Even though themaximal propagation time is in our time-evolution lim- j ( t ) D h ( t ) < n ph > j ( t ) < n ph > D h ( t ) d ( t ) d ( t ) t /t B t /t B (a)(g) (h)F=2(c)(e) (f)(b)F=1F=0.5(d)F=2 F=0.5F=1 FIG. 3: (Color online) j ( t ), h n ph i ( t ), ∆ h ( t ) and d ( t ) vs t/t B for λ = 0 .
9. Left column: F = 0 . t/t B = 15.Right column: F = 0 .
5, 1 and 2 with a maximal t/t B = 8.Thin horizontal line in (b) indicate the quasistationary cur-rent ¯ . Thin lines in (e) and (f) represent a linear depen-dence of ∆ h ( t ), which indicates the quasistationary state.Parameters defining the functional generator of Eq. (2) for λ = 0 . N h = 13, M = 3 and N phmax = 15. ited with the maximal allowed size of the Hilbert space,our results strongly suggest that in the limit when t → ∞ the bipolaron dissociates into two separate polarons.Additional information about the bipolaron dissocia-tion can be obtained if d ( t ) is compared to the averagedistance traveled by the center of mass of two particles∆ x ( t ), where ∆ x ( t ) = ∆ h ( t )2 F . (7)Using ∆ x ( t ) it is convenient to define a ratio η ( t ) = ∆ x ( t )∆ d ( t ) , (8)which represents the ratio between the average distance∆ x ( t ) traveled by two particles along the field direction,and the relative increase of the distance d ( t ) between theparticles defined as ∆ d ( t ) = d ( t ) − d ( t = 0).If the bipolaron propagated without dissociating intoseparate polarons, we would expect a linear increase of η ( t ) in the quasistationary state. Instead, results forboth λ = 0 . . F . For small F ∼ F th , η ( t ) shows damped Bloch oscillations while for larger F the approach toward a constant is rather monotonous.This results suggest that the dissociation of the bipo-laron and the appearance of the quasistationary current h ( t ) t/t b (a) l =0.5 F=0.14F=0.20F=0.25F=0.33F=0.50 0 1 2 3 4 5 2 4 6 8 10 12 14 16 h ( t ) t/t B (b) l =0.9 F=0.33F=0.5F=1.0F=2.0 FIG. 4: The ratio η ( t ) as defined in Eq. (8), vs t/t B . Up-per panel and lower panel correspond to λ = 0 . . emerge simultaneously as the system evolves from a tran-sient regime to the quasistationary state. B. Phonon correlation function
We also compute the average number of phonon quantalocated at a given distance r from (both) electrons. Tostudy this effect we define γ ( r ) = 12 h n ph i h X i,σ n i,σ a † i + r a i + r i , (9)fulfilling the sum rule P r γ ( r ) = 1 . Correlations func-tion γ ( r ) for λ = 0 . .
9) and F = 1 / γ ( r ) with respect to theelectron positions at r = 0 that grows with time. Theasymmetry emerges due to phonon excitations extendingbehind the moving particles, which contain the excessenergy absorbed from the external electric field. Here we FIG. 5: (Color online) log ( γ ( r )) for λ = 0 . .
9) and F = 1 / d at a certain time. The arrow represents the directionof the moving bipolaron. Green dashed–dotted lines representthe corresponding log ( γ ( r )) of the Holstein polaron model. We show log ( γ ( r )) for | r | ≪ N h where results for differentsystem sizes show good convergency. note that the electron is moving in the direction of r > The second fea-ture is the increased amount of phonon excitations in theforward direction, which arises due to two distinct con-tributions. While the first contribution is due to dampedBloch oscillations analogous to the polaron propagation, the second contribution is entirely due to a two–particleeffect, in which case there is always one particle travel-ing ahead in the electric field. It is the second particlewhich follows the first one that in turn detects phononexcitations left by the first particle, thus generating extraweight of γ ( r ) in the region r >
0. Therefore, to obtainthe complete explanation of the increasing γ ( r ) in theforward direction one has to take also into account theincreasing average distance between two particles.For comparison we added results for the polaron case,presented by dashed–dotted lines in Fig. 5, which showa farther extend of γ ( r ) behind the moving polaron anda smaller reach of γ ( r ) in the forward direction, the lat-ter being consistent with the argument given above. Ashorter phonon disturbance behind the traveling bipo-laron indicates a smaller center of mass velocity in com-parison to the polaron. We further elaborate on this issuein the following section. IV. QUASISTATIONARY CURRENT
We next present results for the quasistationary current¯ , its dependence on F and U , and comparison with thepolaron case. In Fig. 6(a) we display the current–fieldcharacteristics for λ = 0 .
5. We have limited our calcula-tions to commensurate values of F = ω /n for F < ω and F = nω for F > ω and n is integer. The cur-rent in the bipolaron case is in this case smaller thanthe current in the polaron case. This is a consequence ofthe fact that in the former case one of the two electronsis traveling in the wake of the other electron where thewake consists of the finite number of phonon excitations.The intuitive explanation is that in this region the ef-fective temperature is elevated that in turn leads to thelowering of the current. To check this assumption we per-formed a numerical test, where we have propagated thepolaron from the ground state until it has reached thequasistationary state. Then we changed the direction ofthe electric field, which prompted the polaron to propa-gate backwards into the phonon rich area. After initialoscillations the polaron has reached the quasistationarystate with a lower current, i.e. , for λ = 0 . is reduced in a two–particle case dueto the arguments given above, and setting η ( t ) ≈ λ = 0 . – F characteristics of a drivenbipolaron from a polaron ¯ – F characteristics. If ∆ ˙ x ( t ) = v (1 + a ) /
2, where v is the polaron velocity in the qua-sistationary state and a < d ( t ) = v (1 − a ). Usinga rough estimate η ≈ a ≈ /
3. This enables us to estimate the two–particlecurrent as ¯ ( F ) ≈ ¯ ( F ). Indeed, these results shownin Fig. 6(a) are in reasonable agreement for F > ∼ . d ( t ) in Fig. 2(d) and 3(h), since the particle the followsthe first one moves with a renormalized but a constantvelocity.At larger lambda λ = 0 . F th ≈ . F th ∼ . F th on modelparameters is better understood when loosely applying asimple Landau-Zener (LZ) formalism to a multi-bandenergy spectrum of the bipolaron, which gives F LZth = (∆ / v , (10)where v is the group velocity of the quasiparticle. Muchlarger F th in the case of bipolaron is attributed tothe much narrower bipolaron bandwidth W bi in com-parison to the polaron W pol , which leads to a largergap ∆ and a lower group velocity v of the quasipar-ticle. We can estimate the scaling of F LZth for bipo-laron in the strong coupling limit. The gap is ∆ ≈ ω (1 − t exp( − λt /ω )) and v is approximated as max-imal value of dE bi /dk , which in the asymptotic expansionscales as v ≈ t exp( − λt /ω ). Therefore to the lead-ing order in λ , F th scales as F th ≈ ( ω /t ) exp(8 λt /ω ). j F (a) l =0.5j , l =0.52h, l =0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.2 0.3 0.5 1 1.5 2 3 j F (b) l=0.9 U=0.0U=1.0U=2.0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.2 0.3 0.5 1 1.5 2 3 j F (b) l=0.9 U=0.0U=1.0U=2.0 0.01 0.1 1 10 0.4 0.6 0.8 1 1.2 1.4 1.6 F t h l (c) F t h l (c) F t h U FIG. 6: (Color online) (a) Quasistationary current ¯ vs elec-tric field F at λ = 0 . (see textfor details). The uncertainty of ¯ is represented with errorbars and is a consequence of transient region; see discussionin text. (b) Quasistationary current ¯ vs electric field F for λ = 0 . U = 0, U = 1, and U = 2. The uncertainty for ¯ is representedwith error bars and is again a consequence of transient re-gion. (c) The threshold electric field F th as a function of theelectron-phonon coupling λ for U = 0. The inset representsdependence of the threshold electric field F th as a function ofelectron-electron repulsion U at fixed λ = 0 . To extract F th from ¯ we used¯ ( F ) = σ F e − πF th /F , (11)where σ and F th are fitting parameters.Nearly linear scaling of log( F th ) with λ is representedin Fig. 6(c). Here we should emphasize that in the leadingorder of λ , F th is increased predominantly due to thenarrowing of the bandwidth and consequently reducingof the group velocity v in Eq. (10).In Fig. 6(b) we also present the influence of the Hub-bard repulsion U on the quasistationary current. Themost important effect is the increase of ¯ with increas-ing U for F < ∼
2. In several recent studies of interactingsystems driven out of equilibrium, an increase of a qua-sistationary current due to stronger correlations has beenobserved.
In the case of the Hubbard–Holstein bipo-laron the effect of U is to reduce the double occupancyand consequently the binding energy δ . Moreover, in theframework of the strong EP coupling limit, U reducesthe effective mass and thus increases the bandwidth, which in turn leads to a decrease of the gap in the bipo-laron excitation spectrum. These combined effects leadto a predominantly linear decrease of F th with increas-ing U in the regime U < ∼
1, as represented in the insetof Fig. 6(c). Nevertheless, F th tends to a constant value F th ∼ . U ∼ . U c ∼
3, see Ref. .We also note that a similar threshold value F th ∼ . − . The bipolaronthus reaches threshold field values of a single polaron evenin the regime when
U < U c .Here we should note that in case of λ = 0 . i.e. , linear de-pendence of energy with time, see Fig. 3(f). ¯ – F depen-dence for λ = 0 . λ = 0 . V. FOURIER ANALYSIS
Following the idea of Refs. we perform the Fourieranalysis of j ( t ) in order to extract additional informationabout the dynamics of the driven system. In doing so wefirst subtract the linear and quadratic terms from j ( t )data in all cases where applicable. In Fig. 7 we repre-sent the normalized Fourier transform (the integral overwhole spectrum is 1) of the real–time current for λ = 0 . π/F and the angular frequency isaccordingly given by ˜ ω = ω ∗ t B .In the (nearly) adiabatic regime, for F = 1 /
3, the bipo-laron exhibits Bloch oscillations, with the doubled Blochfrequency ˜ ω B = 2˜ ω B = 4 π . This is reflected in a wellpronounced peak in the Fourier spectrum of j ( t ) at ˜ ω B .We as well observe well separated peaks at higher har-monics of ˜ ω B , however, with decreasing amplitude.In the strong–coupling picture processes responsiblefor such peaks are diagonal transitions, i.e. , processeswhere bipolaron exhibits coherent propagation throughthe lattice. For example, if the bipolaron jumps as abound pair for one site along the field direction, the en-ergy difference between the initial and the final state is2 F . The related frequency for such a process is ω = 2 F or ˜ ω = 4 π . Jumps for more sites represent higher orderprocesses, which are responsible for generation of higherharmonics. Other peaks are related to different energyscales in the system, such as ω and ∆ as indicated inFig. 7.For F = 1 / F = 1 / ω = 4 π , however, it this case the bipolaron Bloch fre-quency ˜ ω B matches the phonon frequency ω . This is incontrast to the case at F = 1 / j ( t ). An- F ou r w ~ w D w B D w B w F=1/3j 0.001 0.01 0.1 0 10 20 30 40 50 F ou r w ~ w D w B D w B w F=1/3j 0.001 0.01 0.1 F ou r w = w D w B w B F=1/2j 0.001 0.01 0.1 F ou r w = w D w B w B F=1/2j 0.001 0.01 0.1 F ou r w D w B w F=1j
FIG. 7: (Color online) Normalized Fourier transform of thereal–time current j ( t ), where linear and quadratic terms havebeen extracted. We set electron-phonon coupling λ = 0 . F = 1 / , / ,
1. Arrows mark the angularfrequencies corresponding to the distinct energy scales of thesystem. The horizontal dashed line represent the frequencycorresponding to the bipolaron binding energy. other condition for the appearance of the beats is that theelectric field is in the vicinity of the threshold field, i.e. , F ∼ F th , since at larger F increased energy dissipationto phonons overdamps Bloch oscillations. Pronouncedbeats in j ( t ) were observed as well by other authors ex-ploring different driven models. For F = 1 the system reaches the quasistationary statein a time t < ∼ t B with a nearly linear total energy in-crease ∆ h ( t ), see Fig. 3(f), and a rapid increase of d ( t ),see Fig. 3(h). For this reason the Fourier spectrum inFig. 7 of j ( t ) shows less pronounced, broader peaks. Wenote that the frequency corresponding to the bipolaronbinding energy δ as marked with the dashed line is com-parable to ˜ ω B .The influence of the Hubbard repulsion on the dynam-ics of j ( t ) is presented via the Fourier transform in Fig. 8for the case of λ = 0 . F = 1 /
2. The spectrum atfinite U = 1 . U = 0 results. Theeffect of increasing U is in agreement with the consider-ations given is Sec. IV. We again note that the spectrumis broadened when the frequency corresponding to δ be-comes comparable to ω B . Despite much broader spec-trum at U = 1, ¯ remains small, comparable to the U = 0case. At larger U = 2, the peaks in the spectrum becomeindistinguishable, while simultaneously a significant qua-sistationary current appears, as seen in Fig. 6(b). VI. CONCLUSIONS
We have studied the time-evolution of the Holstein-Hubbard bipolaron from its ground state at zero temper- F ou r w ~ w = w D w B w B U=0.0j 0.001 0.01 0.1 0 10 20 30 40 50 F ou r w ~ w = w D w B w B U=0.0j 0.001 0.01 0.1 F ou r w = w D w B w B w B U=1.0j 0.001 0.01 0.1 F ou r w = w D w B w B w B U=1.0j 0.001 0.01 0.1 F ou r w = w D w B w B w B U=2.0j
FIG. 8: (Color online) Normalized Fourier transform of thereal–time current j ( t ), where linear and quadratic terms havebeen extracted, for different values of Hubbard repulsion U .We set λ = 0 . F = 0 .
5. Arrows mark the angularfrequencies corresponding to the distinct energy scales of thesystem. The horizontal dashed line represent the frequencycorresponding to the bipolaron binding energy. ature after a constant electric field has been switched onat t = 0. Using an efficient variational method, definedon an infinite one-dimensional chain, we time-evolved themany–body Schr¨odinger equation while preserving thefull quantum nature of the problem, until the systemhas reached a quasistationary state. In the limit of smallelectric field the bipolaron evolves along the quasiparticleband while the current shows characteristic Bloch oscil-lations. In this regime the propagation is adiabatic, thereis no increase of the total energy, the net current remainszero and the bipolaron remains bound. At larger electricfields we detect an overall increase of the total energy,a finite net current appears and the system enters thedissipative regime. Due to dissipation of the gained po-tential energy into lattice vibrations, Bloch oscillations inthis regime become damped. After a transient time thesystem enters a quasistationary state with a constant cur-rent. We note that there is no clear distinction betweenthe adiabatic and dissipative regime. While a true adi-abatic regime exists only in the limit when F →
0, asizable net current signaling the dissipative regime ap-pears for
F > ∼ F th .The focal result of our work is that in the dissipa-tive regime the bipolaron dissociates into two separatepolarons. By examining different parameter regimes werealized that the appearance of a finite quasistationarycurrent is closely connected with the dissociation of thebipolaron. Even though our calculations were limited tolarge F and short propagation times, our results stronglysuggest that bipolaron has a finite lifetime as long asthe system displays a non-zero ¯ . This hypothesis issupported by the calculation of η ( t ) that in the large-time limit tends to a constant, largely independent of F . This result demonstrates that the dissociation and theappearance of the quasistationary current emerge simul-taneously as the system evolves from a transient regimeto the quasistationary state.Here we should comment that due to dispersionlessphonons the bipolaron can gain energy when two par-ticles occupy the same site, but to achieve a non-zeroquasistationary current charged particles have to starthopping along the direction of electric field, and in doingso breaking the bound pair. In Ref. the Green functionstudy at finite temperature (without electric field) calcu-lated the bipolaron dissociation time, which goes to zerofor dispersionless phonons, while this is not the case forphonons with dispersion. For the problem studied herethis may indicate that phonons with dispersion would notalways lead to dissociation of bipolaron as soon as a finitequasistationary current is obtained. Further research inthis field is necessary to resolve the issue.A linear time–dependence of a real–space expansionof particle densities or related quantities has been ob-served in several recent nonequilibrium problems. Even though we study a two–particle problem on an infi-nite lattice and we do not deal with a well–defined ther-malization, we can consider a cloud of emitted phononexcitations as a subsystem with an elevated effective tem-perature. The quantum Monte Carlo study of the ther-mal dissociation of the Hubbard–Holstein bipolaron pro-vides evidence that a bipolaron dissociates as the tem-perature becomes comparable to the binding energy δ .In comparison, our results indicate that the bipolarondissociates as soon as a finite net current appears around F ∼ F th . This occurs even in the case of δ > ω , whichis realized for λ = 0 . λ ≤ . λ the qualitative behavior remains thesame, the bipolaron as well dissociates with the disso-ciaton time that is increasing with λ and the thresholdfield is also roughly exponentially increased as shown inFig. 6(c).Yet another noteworthy finding reveals that the qua-sistationary current per particle is in case of the disso-ciated bipolaron smaller than in the single polaron casesince the velocity of one particle is diminished due to thescattering on phonons generated by the motion of theother particle. Comparison of the bipolaron ¯ – F charac-teristics at λ = 0 . Acknowledgments
L.V., J.B. and D.G. acknowledge stimulating discus-sions with M. Mierzejewski, V. Zlatiˇc and T. Tohyama. We thank also P. Prelovˇsek for his inspiring remark. Thiswork has been support by the Program P1-0044 of theSlovenian Research Agency (ARRS), REIMEI project,JAEA, Japan, and CINT user program, Los Alamos Na-tional Laboratory, NM USA. T. Oka, R. Arita, and H. Aoki, Phys. Rev. Lett. , 066406(2003). S. Kirino and K. Ueda, J. Phys. Soc. Jpn , 093710(2010). M. Eckstein, T. Oka, and P. Werner, Phys. Rev. Lett. ,146404 (2010). M. Mierzejewski, L. Vidmar, J. Bonˇca, and P. Prelovˇsek,Phys. Rev. Lett. , 196401 (2011). L. Vidmar, J. Bonˇca, M. Mierzejewski, P. Prelovˇsek, andS. A. Trugman, Phys. Rev. B , 134301 (2011). C. Aron, G. Kotliar, and C. Weber, Phys. Rev. Lett. ,086401 (2012). A. Amaricci, C. Weber, M. Capone, and G. Kotliar,arXiv:1106.3483v2 (2011). J. K. Freericks, Phys. Rev. B , 075109 (2008). M. Mierzejewski and P. Prelovˇsek, Phys. Rev. Lett. ,186405 (2010). M. Mierzejewski, J. Bonˇca, and P. Prelovˇsek, Phys. Rev.Lett. , 126601 (2011). D. Karlsson, A. Privitera, and C. Verdozzi, Phys. Rev.Lett. , 116401 (2011). F. Heidrich–Meisner, I. Gonzalez, K. A. Al–Hassanieh, A.E. Feiguin, M. J. Rozenberg, and E. Dagotto, Phys. Rev.B , 205110 (2010). L. Vidmar, J. Bonˇca, T. Tohyama, and S. Maekawa, Phys.Rev. Lett. , 246404 (2011). L.–C. Ku and S. A. Trugman, Phys. Rev. B , 014307(2007). H. Fehske, G. Wellein, and A. R. Bishop, Phys. Rev. B ,075104 (2011). K. K. Thornber and R. P. Feynman, Phys. Rev. B , 4099(1970). J. Bonˇca, S. Maekawa, and T. Tohyama, Rev. B , 035121(2007). J. Bonˇca, S. Maekawa, T. Tohyama, and P. Prelovˇsek,Phys. Rev. B , 054519 (2008). L. Vidmar, J. Bonˇca, S. Maekawa, and T. Tohyama, Phys.Rev. Lett. , 186401 (2009). J. Bonˇca, S. A. Trugman, and I. Batisti´c, Phys. Rev. B ,1633 (1999). J. Bonˇca, T. Katraˇsnik, and S. A. Trugman, Phys. Rev.Lett. , 3153 (2000). J. Bonˇca and S. A. Trugman, Journal of Superconductivity 13, 999 (2000). A. Macridin, G. A. Sawatzky, and M. Jarrell, Phys. Rev.B , 245111 (2004). J. P. Hague, P. E. Kornilovitch, J. H. Samson, and A. S.Alexandrov, Phys. Rev. Lett. , 037002 (2007). M. Berciu, Phys. Rev. B , 081101(R) (2007). D. M. Eagles, R. M. Quick, and B. Schauer, Phys. Rev. B , 054305 (2007). J. T. Devreese and A. S. Alexandrov, Rep. Prog. Phys. ,066501 (2009). O. S. Bariˇsi´c and S. Bariˇsi´c, arXiv:1006.3160v1 (2010). M. Hohenadler, M. Aichhorn, and W. von der Linden,Phys. Rev. B , 014302 (2005). M. Hohenadler and W. von der Linden, Phys. Rev. B ,184309 (2005). T. J. Park and J. C. Light, J. Chem. Phys. , 5870 (1986). S. Fishman, K. Mullen, and E. Ben-Jacob, Phys. Rev. A , 5181 (1990). O. S. Bariˇsi´c, Phys. Rev. B , 214304 (2006). L. Vidmar, J. Bonˇca, and S. A. Trugman, Phys. Rev. B , 104304 (2010). L. Landau, Phys. Z. Sowjetunion 2,46 (1932). C. Zener, Proc. R. Soc. London, Ser. A 137, 696 (1932). A. S. Alexandrov and V. V. Kabanov, Sov. Phys. SolidState , 631 (1986). F. Claro, J. F. Weisz, and S. Curilef, Phys. Rev. B ,193101 (2003). R. Khomeriki, D. O. Krimer, M. Haque, and S. Flach,Phys. Rev. A , 065601 (2010). G. D. Mahan,
Many-particle physics (Plenum Press, NewYork, 1990). Y. Kuroda and D. L. Mills, Phys. Rev. B , 7624 (1985). S. Langer, M. J. A. Schuetz, I. P. McCulloch, U.Schollw¨ock, and F. Heidrich–Meisner, arXiv:1109.4364v1(2011). S. Langer, M. Heyl, I. P. McCulloch, and F. Heidrich–Meisner, Phys. Rev. B , 205115 (2011). T. Oka, H. Aoki,