Transient dynamics of a one-dimensional Holstein polaron under the influence of an external electric field

TL;DR

This study investigates the transient behavior of a one-dimensional Holstein polaron under an external electric field using an advanced variational approach, revealing detailed dynamics, effects of exciton-phonon coupling, and differences in finite systems.

Contribution

It introduces a multi-D₂ Ansatz method for accurately simulating Holstein polaron dynamics, aligning well with hierarchy equations and exploring effects of coupling and system size.

Findings

Bloch oscillations observed from broad wave packets.

Weak exciton-phonon coupling reduces current amplitude.

Strong diagonal coupling causes exciton localization and decay of oscillations.

Abstract

Following the Dirac-Frenkel time-dependent variational principle, transient dynamics of a one-dimensional Holstein polaron with diagonal and off-diagonal exciton-phonon coupling in an external electric field is studied by employing the multi-D$_2$ {\it Ansatz}, also known as a superposition of the usual Davydov D$_2$ trial states. Resultant polaron dynamics has significantly enhanced accuracy, and is in perfect agreement with that derived from the hierarchy equations of motion method. Starting from an initial broad wave packet, the exciton undergoes typical Bloch oscillations. Adding weak exciton-phonon coupling leads to a broadened exciton wave packet and a reduced current amplitude. Using a narrow wave packet as the initial state, the bare exciton oscillates in a symmetric breathing mode, but the symmetry is easily broken by weak coupling to phonons, resulting in a non-zero exciton current. For both scenarios, temporal periodicity is unchanged by exciton-phonon coupling. In particular, at variance with the case of an infinite linear chain, no steady state is found in a finite-sized ring within the anti-adiabatic regime. For strong diagonal coupling, the multi-$\rm D_2$ {\it Anstaz} is found to be highly accurate, and the phonon confinement gives rise to exciton localization and decay of the Bloch oscillations.

Equations

$$\hat{H}={\hat{H}}_{\mathrm{ex}}+{\hat{H}}_{\mathrm{ph}}+{\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{\mathrm{diag}}+{\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{o.d.}$$

$${\hat{H}}_{\mathrm{ex}}$$

$${\hat{H}}_{\mathrm{ph}}$$

$${\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{\mathrm{diag}}$$

$${\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{o.d}$$

$${\hat{H}}_{\mathrm{ph}}$$

$${\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{\mathrm{diag}}$$

$${\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{o.d.}$$

$${\hat{b}}_q^{†}=N^{-1/2}\sum _n{e}^{iqn}{\hat{b}}_n^{†},{\hat{b}}_n^{†}=N^{-1/2}\sum _q{e}^{-iqn}{\hat{b}}_q^{†}.$$

$$|{\mathrm{D}}_2^{\mathrm{M}}\left(t\right)⟩=\sum _i^M\sum _n^N{\psi }_{in}|n⟩|{\lambda }_i⟩,$$

$$=\sum _i^M\sum _n^N{\psi }_{in}{\hat{a}}_n^{†}{|0⟩}_{\mathrm{ex}}\exp \left\{\sum _q\left[{\lambda }_{iq}{\hat{b}}_q^{†}-{\lambda }_{iq}^{\ast }{\hat{b}}_q\right]\right\}{|0⟩}_{\mathrm{ph}},$$

$$\frac{d}{dt}\left(\frac{\partial L_2}{\partial \dot{{\psi }_{in}^{\ast }}}\right)-\frac{\partial L_2}{\partial {\psi }_{in}^{\ast }}$$

$$\frac{d}{dt}\left(\frac{\partial L_2}{\partial \dot{{\lambda }_{iq}^{\ast }}}\right)-\frac{\partial L_2}{\partial {\lambda }_{iq}^{\ast }}$$

$$L_2$$

$$\frac{i}{2}\left[⟨{\mathrm{D}}_2^M(t)|\frac{\overrightarrow{\partial }}{\partial t}|{\mathrm{D}}_2^M(t)⟩-⟨{\mathrm{D}}_2^M(t)|\frac{\overleftarrow{\partial }}{\partial t}|{\mathrm{D}}_2^M(t)⟩\right]$$

$$=\frac{i}{2}\sum _{i,j}^M\sum _n\left({\psi }_{jn}^{\ast }{\dot{\psi }}_{in}-{\dot{\psi }}_{jn}^{\ast }{\psi }_{in}\right)S_{ji}$$

$$+\frac{i}{2}\sum _{i,j}^M\sum _n{\psi }_{jn}^{\ast }{\psi }_{in}{S}_{ji}\sum _q[\frac{{\dot{\lambda }}_{jq}^{\ast }{\lambda }_{jq}+{\lambda }_{jq}^{\ast }{\dot{\lambda }}_{jnq}}{2}$$

$$-\frac{{\dot{\lambda }}_{iq}{\lambda }_{iq}^{\ast }+{\lambda }_{iq}{\dot{\lambda }}_{iq}^{\ast }}{2}+{\lambda }_{jq}^{\ast }{\dot{\lambda }}_{iq}-{\lambda }_{iq}{\dot{\lambda }}_{jq}^{\ast }],$$

$$⟨{\mathrm{D}}_2^M\left(t\right)|\hat{H}|{\mathrm{D}}_2^M\left(t\right)⟩$$

$$=⟨{\mathrm{D}}_2^M\left(t\right)|{\hat{H}}_{\mathrm{ex}}|{\mathrm{D}}_2^M\left(t\right)⟩+⟨{\mathrm{D}}_2^M\left(t\right)|{\hat{H}}_{\mathrm{ex}}|{\mathrm{D}}_2^M\left(t\right)⟩$$

$$+⟨{\mathrm{D}}_2^M\left(t\right)|{\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{\mathrm{diag}}|{\mathrm{D}}_2^M\left(t\right)⟩+⟨{\mathrm{D}}_2^M\left(t\right)|{\hat{H}}_{\mathrm{ex}-\mathrm{ph}}^{o.d.}|{\mathrm{D}}_2^M\left(t\right)⟩,$$

$$S_{ij}=\exp \left\{\sum _q{\lambda }_{iq}^{\ast }{\lambda }_{jq}-\frac{1}{2}\left({\left|{\lambda }_{iq}\right|}^2+{\left|{\lambda }_{jq}\right|}^2\right)\right\}.$$

$$P_{\mathrm{ex}}(t,n)$$

$$X_{\mathrm{ph}}(t,n)$$

$$c(t)$$

$${\sigma (t)}^2$$

$$\rho \left(n,t=0\right)={\left[{\left(2\pi \right)}^{1/2}{\sigma }_0\right]}^{-1}\exp \left(-n^2/2{\sigma }_0^2\right)$$

$$\psi \left(n,t=0\right)={\left[{\left(2\pi \right)}^{1/2}{\sigma }_0\right]}^{-1/2}\exp \left(-n^2/4{\sigma }_0^2\right)$$

$$H_0=-J\sum _n\left(a_{n+1}^{†}{a}_n+a_n^{†}{a}_{n+1}\right),$$

$$H_{int}=F\sum _n{a}_n^{†}n{a}_n,$$

$$\tilde{\psi }\left(t\right)=e^{-iA\left(t\right)x}\psi \left(t\right)$$

$$i\frac{\partial \psi \left(t\right)}{\partial t}=\left(H_0+H_{int}\right)\psi \left(t\right)$$

$$i\frac{\partial \tilde{\psi }\left(t\right)}{\partial t}=\tilde{H}\tilde{\psi }\left(t\right)$$

$$i\frac{\partial \tilde{\psi }\left(t\right)}{\partial t}$$

$$=i\frac{\partial \left(-iA\left(t\right)x\right)}{\partial t}{e}^{-iA\left(t\right)x}\psi \left(t\right)+e^{-iA\left(t\right)x}i\frac{\partial \psi \left(t\right)}{\partial t}$$

$$=-Fx{e}^{-iA\left(t\right)x}\psi \left(t\right)+e^{-iA\left(t\right)x}i\frac{\partial \psi \left(t\right)}{\partial t}$$

$$=-H_{int}{e}^{-iA\left(t\right)x}\psi \left(t\right)+e^{-iA\left(t\right)x}\left(H_0+H_{int}\right)\psi \left(t\right)$$

$$=e^{-iA\left(t\right)x}{H}_0\psi \left(t\right)$$

$$\tilde{H}\tilde{\psi }\left(t\right)=\tilde{H}{e}^{-iA\left(t\right)x}\psi \left(t\right)$$

$$e^{-iA\left(t\right)x}{H}_0\psi \left(t\right)=\tilde{H}{e}^{-iA\left(t\right)x}\psi \left(t\right)$$

Full text

I Introduction

Early evidence of time-domain Bloch oscillations (BOs) in superlattices was revealed by transient degenerate four-wave mixing (DFWM) fe_92 ; me_93 , yielding results in agreement with the work of Esaki and Tsu es_70 . BOs emerge when an electron subject to a perfectly periodic lattice potential executes periodic motion in the presence of an external electric field bl_28 ; ze_34 ; ch_83 ; wa_59 ; ho_40 ; Bloch1 ; Bloch2 . In recent decades, BOs were found in various physical systems. For example, with the progress of laser cooling and manipulation, cold atoms in an optical lattice were found to exhibit BOs and Wannier-Stark ladders da_87 ; ch_91 , in accordance with theoretical predictions da_96 ; wi_96 ; mo_01 . Further experimental evidence of BOs was found in atomic Bose-Einstein condensates (BEC) in optical lattices an_98 ; ko_10 ; cr_02 , in a semiconductor waveguide mo_99 , and in a thermo-optic polymer array subjected to a temperature gradient pe_99 . It was also suggested that waveguide arrays with a changeable effective index of the individual guides would be an optimal system to detect optical BOs in the space domain pe_98 . In addition, BOs have received much attention in the past decades for its potential application in THz generation and negative differential conductance es_70 . Investigations of BOs have been carried out in THz emission wa_93 , electro-optic detection fo_2000 , and DFWM experiments fe_92 . Recently, a new mesoscopic amplifier by the name of Bloch Oscillating Transistor has been proposed based on BOs de_03 .

Following the remarkable experimental progress, recent theoretical work explores the presence of BOs in a variety of contexts. For example, formation of photonic BOs was investigated in an exponentially chirped one-dimensional Bragg grating using Hamiltonian optics, where paths of geometrical rays are determined from Hamilton’s equations wi_02 . BOs were also theoretically predicted to exist in magnetic systems, such as soliton-like domain walls in anisotropic spin $1/2$ chains under magnetic fields kr_98 . Furthermore, BOs in interacting quantum few body systems have been modeled with the Bose-Hubbard model bu_03 .

One important issue for BOs dynamics is the effect of carrier-phonon interactions, which is of essential importance for systems such as semiconductor superlattices dekorsy_2000 ; ghosh_2000 , organic materials jo_01 ; ba_02 ; jo_04 ; qu_09 ; pe_89 ; di_08 ; Rakhmanova_1999 and quantum dots dmitriev_2001 ; petroff_2001 . The phonon bath is treated classically and the associated carrier-phonon coupling semi-classically in the Su-Schrieffer-Heeger model for polymers jo_01 ; ba_02 ; jo_04 ; qu_09 and in the Peyrard-Bishop-Holstein model for DNA pe_89 ; di_08 . More than forty years ago, Thornber and Feynman studied the motion of an electron in a polar crystal in a strong electric field using Fröhlich’s model of polaron, and found that the electron acquires a constant velocity due to the emission of phonons th_70 . Much theoretical work based on the rate equations or the Boltzmann equation later focused on calculating transition rate probabilities, rather than complex quantum amplitudes kummel_1978 ; Emin1987 . However, full quantum coherence was revealed to be retained in an inelastic quantum transport process and a steady state was found to be reached subject to an infinite lattice, leaving the dispersionless optical phonon absorbing the excess energy from the external field bo_97 ; vi_11 .

Despite dedicated studies of BOs dynamics, the effect of complex interplay between the electron and its accompanying phonon cloud in a one-dimensional lattice remains an open question zhang_02 ; cheung_13 . Spatial dynamics influenced by exciton-phonon coupling is also inadequately studied given initial broad and narrow exciton wave packets ly_1997 ; su_98 . Moreover, a unified treatment on various types of exciton-phonon coupling remains elusive.

The intramolecular (diagonal) and intermolecular (off-diagonal) exciton-phonon coupling has been demonstrated to coexist in organic materials coropceanu_07 , and it is shown that off-diagonal coupling plays a crucial role in polaronic diffusion tamura_2012 . In the presence of an external electric field, polaron motion with off-diagonal coupling in polymer chains has been simulated with a semiclassical method Rakhmanova_1999 ; Johansson_2002 , neglecting the quantum nature of phonons and exciton-phonon coupling. Fully quantum mechanical treatments are few in the literature with little attention paid to Hamiltonians with off-diagonal exciton-phonon coupling and an external field due to inherent difficulties to obtain reliable solutions vi_11 .

Based on the Holstein molecular crystal model, which describes the motion of an exciton enveloped by a cloud of phonons ho_59 , we include both diagonal and off-diagonal exciton-phonon coupling. Despite apparent simplicity of its Hamiltonian, the Holstein model never ceases to surprise us with rich physics related to exciton-phonon correlations. The off-diagonal coupling was investigated earlier by the Munn-Silbey theory mu_85 and a variational method using the Davydov D2 Ansatz zh_12 and the global-local Ansatz zh_97 . Recently, Zhao et al. developed a refined trial state, the multiple Davydov D2 Ansatz, to accurately treat the dynamics of the Holstein model with simultaneous diagonal and off-diagonal coupling zh_15 ; zhou_2016 . In this work we will apply this method to polaron dynamics in an external field.

In addition, only interactions between the exciton and optical phonons were considered previously th_70 ; vi_11 ; zh_15 . Quantum-mechanical calculations indicated that excitons are coupled to both acoustic and optical phonons bredas_2012 ; li_2013 . Discussions on the role played by acoustic modes are inadequate in the literature, especially in the simultaneous presence of acoustic and optical phonons accompanied by an external field davydov_1967 ; la_1997 .

In this work, we investigate the impact of diagonal and off-diagonal exciton-phonon coupling on polaron dynamics in a ring system in the presence of a constant electric field using the multiple Davydov D2 trial state with the Dirac-Frenkel variational principle. Accuracy of the method is verified by the comparison with benmark calculations obtained from the numerically exact hierarchical equations of motion (HEOM) method Tanimura1 ; ch_15 which is also firstly used in this work to deal with field-driven cases to the best of our knowledge. The validity of our variational method is also carefully examined by quantifying how faithfully our result follows the Schrödinger equation in balance with the computational efficiency. For weak exciton-phonon coupling in the anti-adiabatic regime, we consider both acoustic and optical phonons, and calculate the time evolution of various quantities with special attention paid to the effect of weak dissipation on the spatial amplitudes of BOs using different initial conditions. Treating strong diagonal coupling is a formidable challenge for the HEOM method, but our variational method using the multi-$\rm D_{2}$ Ansatz remains numerically affordable (in fact, with a higher precision with stronger coupling).

The rest of the paper is structured as follows. In Sec. II, we present the Hamiltonian in the gauge transformed form and our trial wave function, the multi-$\rm D_{2}$ Ansatz. In Sec. III.1, our variational results are compared to those of the benchmark HEOM calculations. In Sec. III.2, the influence of weak diagonal coupling in the anti-adiabatic regime on the exciton wave packet is investigated using initial Gaussian wave packets with varying widths with particular attention paid to the possible existence of a steady state. Finally, polaron dynamics is examined in the strong diagonal coupling regime in Sec. III.3. Conclusions are drawn in Sec. IV.

II methodology

II.1 Model

The Hamiltonian of the one-dimensional Holstein polaron reads

[TABLE]

where $\hat{H}_{\rm ex},\hat{H}_{\rm ph}$,$\hat{H}_{\rm ex-ph}^{\rm diag}$ and $\hat{H}_{\rm ex-ph}^{\rm o.d.}$ denote the exciton Hamiltonian, the bath (phonon) Hamiltonian, the diagonal and off-diagonal exciton-phonon coupling Hamiltonian, respectively. An extra term $F\sum_{n}na_{n}^{\dagger}a_{n}$ is added to represent the scalar potential induced by a constant external electric field $F$. After bending a linear chain of atoms into a ring, however, potential of the field would become discontinuous at end points, leaving the use of periodic boundary conditions questionable za_68 . A later treatment kr_86 using a gauge transformed vector potential ki_63 avoids this problem and is thus suitable for ring systems. We will use the gauge transformed Hamiltonian for the ring system under a constant electric field (see Supporting Information for details), which can be defined as

[TABLE]

where $\hat{a}_{n}^{\dagger}$ ($\hat{a}_{n}$) and $\hat{b}_{n}^{\dagger}$ ($\hat{b}_{n}$) are the exciton and phonon creation (annihilation) operators for the $n$-th site, respectively. Written in the phonon momentum space,

[TABLE]

where $\omega_{q}$ is the phonon frequency with momentum $q$, and $\hat{b}_{q}^{\dagger}$ ($\hat{b}_{q}$) is the creation (annihilation) operator of a phonon with momentum $q$,

[TABLE]

The parameters $J$, $g$ and $\phi$ represent the transfer integral, the diagonal coupling strength and off-diagonal coupling strength, respectively. In this work, we take $\omega_{q}=\omega_{0}\left|\sin\left(q/2\right)\right|$ as the dispersion relation for acoustic phonons, where $q=2\pi l/N$ represents the momentum index with $l=-\frac{N}{2}+1,\ldots,\frac{N}{2}$; in the case of optical phonons, we consider the Einstein dispersionless model, i.e., $\omega_{q}=\omega_{0}$. For simplicity, the Debye frequency is taken to be equal to the Einstein frequency $\omega_{0}$. In the remainder of the paper, $\omega_{0}$ is set to unity as the energy unit.

II.2 Multiple Davydov trial states

The multiple Davydov trial states with multiplicity $M$, which are essentially $M$ copies of the corresponding single Davydov Ansatz, have been developed to investigate the time evolution of the Holstein polaron following the Dirac-Frenkel variational principle. The multi-${\rm D}_{2}$ Ansatz has less variational variables than the multi-${\rm D}_{1}$ Ansatz when same $M$ is used, nonetheless it can more accurately describe polaron dynamics in the presence of off-diagonal coupling zh_15 . In this work, we employ the multi-${\rm D}_{2}$ Ansatz. It can be constructed as

[TABLE]

where $\psi_{in}$ and $\lambda_{iq}$ are time-dependent variational parameters for the exciton probability and phonon displacement, respectively, $n$ represents the site number, and $i$ labels the coherent superposition state. If $M=1$, the $|\rm D_{2}^{M}(t)\rangle$ Ansatz is restored to the usual Davydov $\rm D_{2}$ trial state. The equations of motion of the variational parameters $\psi_{in}$ and $\lambda_{iq}$ are then derived by adopting the Dirac-Frenkel variational principle,

[TABLE]

for the multi-$\rm D_{2}$ Ansatz, and the Lagrangian $L_{2}$ is formulated as

[TABLE]

where the first term yields

[TABLE]

and the second term takes the form

[TABLE]

where the Debye-Waller factor is formulated as

[TABLE]

Detailed derivations of the equations of motion for the variational parameters are given in Supporting Information. In the numerical calculations, the fourth-order Runge-Kutta method is used to integrate the equations of motion.

With the wave function $|{\rm D}^{M}_{2}(t)\rangle$ available, the total energy $E_{\rm total}=E_{\rm ex}+E_{\rm ph}+E_{\rm diag}+E_{\rm off}$ is calculated, where $E_{\rm ex}=\langle{\rm D}^{M}_{2}|\hat{H}_{\rm ex}|{\rm D}^{M}_{2}\rangle,\leavevmode\nobreak\ E_{\rm ph}=\langle{\rm D}^{M}_{2}|\hat{H}_{\rm ph}|{\rm D}^{M}_{2}\rangle,\leavevmode\nobreak\ E_{\rm diag}=\langle{\rm D}^{M}_{2}|\hat{H}_{\rm ex-ph}^{\rm diag}|{\rm D}^{M}_{2}\rangle$, and $E_{\rm off}=\langle{\rm D}^{M}_{2}|\hat{H}_{\rm ex-ph}^{\rm o.d.}|{\rm D}^{M}_{2}\rangle$. Additionally, time evolution of the exciton probability $P_{\rm ex}(t,n)$ and the phonon displacement $X_{\rm ph}(t,n)$ is also calculated

[TABLE]

We further use a standard deviation $\sigma(t)$ and a mean value $c(t)$ of the exciton wave packet to characterize the motion of the exciton wave, which are defined as follows

[TABLE]

where the mean value $c(t)$ describes the centroid of the exciton wave packet, and the standard deviation is used to measure how far the exciton wave spreads out from the mean position. It is noted that the standard deviation $\sigma(t)$ and the mean value $c(t)$ are sensitive to the initial standard deviation $\sigma_{0}$ of the exciton wave packet.

III Numerical results and discussions

III.1 Validity of variation for transient dynamics

In this subsection, we show that the multi-$\rm D_{2}$ Anstaz with sufficiently large multiplicity $M$ yields quantitatively accurate solutions to the dynamics of the Holstein polaron with both diagonal and off-diagonal coupling, in perfect agreement with the benchmark calculations of the numerically exact HOEM method (see Sec. $4$ of Supporting Information) Tanimura1 ; ch_15 . For simplicity, only optical phonons are used in this subsection.

III.1.1 Diagonal Coupling

First, we study the case of diagonal exciton-phonon coupling. Dynamics of the Holstein polaron under a constant external field is investigated by using the multi-$\rm D_{2}$ Anstaz, in comparison with those obtained with the single Davydov D2 Ansatz and the HEOM method. Using these approaches, the time evolution of the exciton probability $P_{\rm ex}(t,n)$ as shown in Fig. 1 is simulated in the case of $J=0.1,g=0.28$ and $F=0.1$ in a small ring with $N=8$ sites for simplicity. The exciton is created on two nearest neighboring sites $\psi_{n}=(\delta_{n,N/2}+\delta_{n,N/2+1})/\sqrt{2}$, and the phonon displacement $\lambda_{i,q}(t=0)=0$ is set. As depicted in Figs. 1(a) and (b), distinguishable deviations in $P_{\rm ex}(t,n)$ can be found between the variational results from the $\rm D^{M=1}_{2}$ and $\rm D^{M=16}_{2}$ Ansätze. Interestingly, $P_{\rm ex}(t,n)$ obtained from the HEOM method in Fig. 1(c) is almost identical to that by the $\rm D_{2}^{M=16}$ Ansatz in Fig. 1(b). Furthermore, the difference in the time evolution of the exciton probability between the variational and HEOM methods, $\Delta P_{\rm ex}(t,n)$, as displayed in Fig. 1(d), is two orders of magnitude smaller than the value of $P_{\rm ex}(t,n)$, indicating that the variational dynamics of the Holstein polaron under the external field can be numerically exact if the multiplicity $M$ of the multi-$\rm D_{2}$ Ansatz is sufficiently large.

It should be noted that the HEOM method is numerically expensive and thus impractical when the system size is large, while time dependent variational approaches are still valid to treat the polaron dynamics for large systems once a proper trial wave function is adopted.

III.1.2 Off-diagonal Coupling

We extend our discussion to off-diagonal coupling, which is a formidable problem due to the intrinsic difficulties in achieving reliable results. The off-diagonal coupling was emphasized as modulations of electron-electron interactions by ion vibrations in Mahan’s celebrated textbook on many-particle physics mahan . Due to a lack of dependable solutions, a complete understanding of out-of-equilibrium dynamics for off-diagonal coupling remains elusive. Below we start with a validity check of our variational method.

Dynamics of the Holstein polaron with off-diagonal coupling under a constant external field is examined by using the multi-$\rm D_{2}$ Ansatz with different multiplicity $M$. As shown in Figs. 2(a)-(c) and (e), the time evolution of the exciton probability $P_{\rm ex}(t,n)$ is obtained by the $\rm D_{2}^{M=1}$, $\rm D_{2}^{M=8}$, $\rm D_{2}^{M=20}$ Ansätze and the HEOM method, respectively. Figs. 2(b),(c) and (e) display quite similar patterns and all three are largely different from Fig. 2(a). Difference between $P_{\rm ex}(t,n)$ obtained by the $\rm D_{2}^{M=8}$ and $\rm D_{2}^{M=20}$ Ansatz as shown in Fig. 2(d) is two orders of magnitude smaller than the value of $P_{\rm ex}(t,n)$, pointing to the nearly converged results already obtained by $M=8$. Moreover, the difference in $P_{\rm ex}(t,n)$ between the $\rm D_{2}^{M=20}$ Ansatz and the HEOM method, $\Delta P_{\rm ex}(t,n)$, as depicted in Fig. 1(f), is also two orders of magnitude smaller than the value of $P_{\rm ex}(t,n)$, showing the superior accuracy of the multi-D2 Ansatz.

In addition, a quantity named the relative deviation $\Sigma$ (see Supporting Information) is also used to test the validity of our time-dependent variational approach by quantifying how closely the trial state follows the Schrödinger equation, as depicted in Fig. 3. As the multiplicity $M$ increases, the relative deviation $\Sigma$ decreases and approaches zero as $M$ goes to infinity. This is supported by the relationship $\Sigma\sim M^{-\mu}$ with an exponent of $\mu=3.15(1)$ in the inset of Fig. 3. Therefore, in the limit of large $M$, our variational method using the multi-$\rm D_{2}$ Ansatz provides a numerically exact solution to the Schrödinger equation in the presence of the external field.

III.2 Anti-adiabatic regime

The effect of a bosonic environment on BOs, i.e., BOs in a polaron framework, is further investigated in this subsection. As usually envisioned for electronic transport in crystals, motion of electrons is occasionally scattered by lattice vibrations th_70 . Pronounced modulations of the four-wave mixing signal with characteristics of the temporal periodicity of BOs has been detected experimentally in an optical investigation of BOs in a semiconductor superlattice, and can be attributed to lattice scattering fe_92 . Bouchard et al. excluded interband transitions from being responsible for the signal modulations and confirmed the single-band model as an appropriate approximation bo_95 . Several pioneering theoretical studies also demonstrated the modulations of BOs by the electron-phonon interaction th_70 ; vi_11 ; lakhno_2004 . Despite tremendous experimental and theoretical efforts dedicated to the phonon modulated BOs, the underlying mechanisms are still not well understood.

To be specific, both the acoustic and optical phonons have been experimentally revealed to coexist in the semiconductor superlattice kelvin_1986 , and amplitudes of exciton wave packets has been confirmed to be very sensitive to the precise excitation conditions in a weak external field ly_1997 . However, there is a lack of investigation on BOs dynamics influenced by the two phonon branches when the exciton starts from a Gaussian wave packet with varying widths.

In this subsection we will focus on dynamic properties of the Holstein polaron including the time evolution of the exciton probability and the phonon displacement for weak exciton-phonon coupling subject to a constant external electric field. Two typical scenarios are discussed by considering initial Gaussian wave packets with initial widths $\sigma_{0}=1$ and $0.2$, as shown in Fig. 4 and Fig. 5, respectively. The initial Gaussian distribution st_2015 is described as follows

[TABLE]

which is centered at site $n=0$, and the associated initial wave function is taken as the square root of the distribution:

[TABLE]

In particular, the roles played by the acoustic and optical phonons on the exciton transport are explicitly compared under the same external field. It is difficult to partition these branches rigorously bredas_2012 , and in order to avoid the complexity induced by partitioning these branches, we include only one phonon branch in each separate calculation and then compare the results of two calculations.

We start from using an initial broad Gaussian wave packet of $\sigma_{0}=1$. Transfer integral $J=0.1$ and an external field of $F=0.1$ is set in the calculations of this subsection. In the absence of exciton-phonon coupling, as shown in Fig. 4(a), the exciton exhibits typical BOs, with the center of mass of the wave packet oscillates while its shape is essentially unchanged, in agreement with earlier theoretical work based on an ideal GaAs/AlxGa1-xAs superlattice in a uniform electric field as described by a conventional flat-band picture st_2015 ; bo_95 . In the presence of the weak diagonal coupling ($g=0.4$), we focus on the anti-adiabatic regime where the phonon frequencies are larger than the transfer integral J. Only $P_{\rm ex}(t,n)$ of optical phonons is shown in Fig. 4(b) because the effect of the two phonon branches turned out to be similar. By comparing Figs. 4(a) and (b), the time periods of the exciton transport are the same because the temporal periodicity is determined by the external field. This agrees with the experimental observation that the detected signal is found to be modulated over time but the time period of the signal is found to be equal to the temporal periodicity of BOs fe_92 . The largest oscillation amplitude of the center of mass of the exciton wave packet in Fig. 4(a) is in accordance with the theoretical value of $4J/F$ ha_04 .

Even though the temporal periodicity is preserved after the exciton-phonon coupling is turned on, the spatial periodicity is changed over time. The addition of the exciton-phonon coupling moves the center of mass of the exciton wave packet of $g=0$ closer to the initial location, which is in line with the contrast in the mean value $c(t)$ between $g=0$ and $0.4$ in Fig. 6(a). Not only the motion of the center of mass of the exciton wave packet is changed, but also the width of the wave packet is enlarged as another characteristics of an effect of the weak coupling, giving rise to a broadened wave packet. The width $\sigma(t)$ in Fig. 6(c) is increased due to smearing of the wave packet.

The phonon displacement $X_{\rm ph}(t,n)$ is shown in Figs. 4(c) and (d) for the acoustic and optical phonons, respectively. Before the exciton creation at $t=0$, the phonons are in their vacuum states for both scenarios. By comparing the two figures, parts of $X_{\rm ph}(t,n)$ (red and yellow) are found to propagate with BOs frequency $\omega_{B}$ because they are generated by the moving exciton wave packet and would in turn smear out the exciton wave fronts. As for the remainder of $X_{\rm ph}(t,n)$ (blue), a V-shaped feature and an oscillatory component with the phonon frequency $\omega_{0}$ are respectively found in Figs. 4(c) and (d). Moreover, the existence of ten peaks in one Bloch period $t_{B}$ can be attributed to the ratio $\omega_{0}/\omega_{B}$ (phonon frequency $\omega_{0}$ over BOs frequency). In the presented cases, $\omega_{0}$ is ten times of $\omega_{B}=de|F|/\hbar$, where $d$ is the lattice constant, $e$ is the charge of the exciton, and $d=e=\hbar=1$ is set. Meanwhile, the weak coupling affects the transport and the exciton current in Fig. 6(e). At zero coupling, the exciton current $j(t)$ (see Supporting Information) in the case of an initial broad Gaussian wave packet is $j(t)=2\sin Ft$ with the largest amplitude among all cases studied. The amplitude of the exciton current is decreased after the coupling is switched on, mitigating unidirectional exciton transport.

Next, we consider the effect the weak exciton-phonon coupling on the dynamics with an initial narrow Gaussian wave packet of $\sigma_{0}=0.2$. As shown in Fig. 5(a), the exciton undergoes a symmetric breathing mode at zero coupling. The exciton wave packet propagates with its center of mass fixed at the original location and its width oscillates with the Bloch period, in accord with previous studies on the breathing mode ly_1997 ; bo_95 ; do_10 . After the addition of the weak diagonal coupling $g=0.4$, the exciton amplitude becomes larger in the left branch than that in the right branch (Fig. 5(b)) and the center of mass of the exciton wave packet is moved away from the original location. These features of the exciton indicate that the symmetry of the breathing mode is broken by even minute exciton-phonon interactions, and are further reflected by the mean value $c(t)$ plotted in Fig. 6(b), which oscillates for the weak coupling case of $g=0.4$ in contrast to being zero at all times for $g=0$. Recently, the breakdown of the breathing mode has been confirmed in semiconductor superlattices di_94 and cold atoms trapped in optical lattice Trombettoni2001 ; ga_09 ; ga_11 , making our results here highly relevant in the study of breathing modes in those systems. Furthermore, the width of the exciton wave packet for the weak coupling case of $g=0.4$ is reduced in comparison to that for $g=0$, as depicted in Fig. 6(d). A comparison between Figs. 6(c) and (d) reveals that weak diagonal coupling has opposite influences on the width of the exciton wave packet: the exciton wave packet is broadened in the case of $\sigma_{0}=1$ while that in the case of $\sigma_{0}=0.2$ is suppressed after the exciton-phonon coupling is switched on.

Similar to the scenario of $\sigma_{0}=1$, the phonon displacement $X_{\rm ph}(t,n)$ of $\sigma_{0}=0.2$ propagates along with the movement of the exciton. As shown in Fig. 5(c), $X_{\rm ph}(t,n)$ maintains its V-shaped feature in the case of acoustic phonons. As for Fig. 5(d), characteristic oscillations with frequency $\omega_{0}$ is found in the case of optical phonons. Reciprocally, the generated $X_{\rm ph}(t,n)$ induces smearing of the exciton in its center. A brief comparison between Fig. 4 and Fig. 5 reveals that the exciton triggered phonon displacement leads to the weakening of the exciton wave packet in its center and edge, respectively. As shown in Fig. 6(f), the current in the zero coupling case is zero due to the spatial symmetry of the exciton wave packet. The exciton current in the weak coupling exhibits fast beating with the characteristic frequency $\omega_{0}$ superimposed by slower BOs.

So far, we have focused on the impact of weak exciton-phonon coupling on the quantum transport with respect to the two typical scenarios. In general, BOs are very sensitive to any kind of dephasing generated by the electron-hole Coulomb interaction effects or lattice imperfections, since they rely on the coherent reflection of waves ga_11 . The electron-hole Coulomb interaction is pointed out to destroy the breathing mode by Dignam et al. after they examined the nature of the exciton wave packets in undoped semiconductor superlattices in a uniform along-axis electric field di_94 . The effects of atom-atom interactions in the dilute BEC trapped in a periodic potential were actively studied using a discrete nonlinear Schrödinger equation (DNLSE), and the atom-atom interactions can result in the decoherence of atomic BOs Trombettoni2001 ; kolovsky_09 ; ko_10 and enhancement or suppression of the breathing width ga_09 ; ga_11 . To the best of our knowledge, ours is the first investigation on the influence of the exciton-phonon interaction on both the typical BOs and the breathing mode.

Finally, we analyze the characteristics of the nonequilibrium dynamics of the ring system, including a comparison with the steady state of an infinite linear chain. First of all, in the absence of an external electric field, in the anti-adiabatic limit of $\omega_{0}\gg J$, exciton dynamics in both weak and strong coupling regimes is dominated by coherent oscillations and negligible energy transfer, leaving a conserved total energy do_15 . After the an external field is applied, the total energy is not conserved as plotted in Fig. 7 because of the addition of the extra energy acquired from the external field. Secondly, a steady state with a saturated current can be obtained in a linear chain that extends to infinity in both directions, where the exciton acquires kinetic energy from the applied field while dissipating its energy by a net emission of optical phonons th_70 . Being a steady state situation, the expectation value of the rate of change of electron momentum is zero, and the steady state is left with an equation which balances the applied field, or rate of increase of electron momentum, against the rate of loss of momentum due to the lattice scattering. This result is an explicit field dependence in terms of the steady-state velocity of the electron. The key feature of the steady state is also explored in Ref. [43; 81; 82]. Equivalently, the condition of the steady state can be that there exist linear regimes of the total energy and phonon energy following the equality of $dE_{\rm total}/dt=dE_{\rm ph}/dt$ vi_11 , meaning that the total energy gain is entirely absorbed by the lattice. Specifically, with increasing time $t$, the total energy vs time approaches a straight line, and the phonon energy has a linear time dependence as well. But in our ring system, the condition of the steady state can not be met as this equation can not be satisfied from the curves of $E_{\rm total}$ and $E_{\rm ph}$ as shown in Fig. 7. For example, the value of $dE_{\rm ph}/dt$ at $t=2t_{B}$ is not equal to but $9$ times of the value of $dE_{\rm total}/dt$. The energy imported by the external field can not cancel out the energy emitted to the phonons, leaving the exciton energy oscillates over time. Thus, the electron’s acceleration in the ring can no longer be balanced out by the lattice’s deceleration and the electron can not drift with a saturated velocity. Consequently, no saturated constant current is found in the anti-adiabatic regime for the small ring systems, and steady states in the adiabatic regime will be investigated elsewhere.

To sum up, weak exciton-phonon coupling in the anti-adiabatic regime breaks the spatial periodicity but retains the temporal periodicity of both the typical BOs and the breathing mode. The exciton movement shows similar features in spite of the difference between the acoustic and optical phonon branches. In fact, in a weak external field, when diagonal coupling is strong between the exciton and both the acoustic and optical phonons, the exciton become localized, leading to the diminishing of the typical BOs and the breathing mode. The localization is similar to that in the moderate external field cases in the next subsection, therefore the details are not shown. We also show that no saturated constant current is found in the ring system for the duration of our simulation.

III.3 Strong diagonal coupling

In this subsection we seek to present a general picture of the exciton wave packet evolution in the regime of the strong diagonal coupling. We first investigate the exciton wave packets as well as the exciton currents under a moderate external field of $F=1$. The time evolution of the exciton probability $P_{\rm ex}(t,n)$ for different diagonal coupling strengths $g=0$, $3$ and $4$ is illustrated in Figs. 8(a), (b) and (c), respectively. The exciton is created on two nearest neighbouring sites. We set $J=1$ and $N=16$. As shown in Figs. 8(a)-(c), the bare exciton exhibits a partially BOs pattern, while the strong exciton-phonon coupling is found to localize the exciton to the initial excitation sites. BOs are largely quenched due to strong exciton-phonon coupling lakhno_2007 . This behaviour is also demonstrated in Fig. 8(d), which shows that the amplitude of the exciton current decreases with increasing exciton-phonon coupling and fluctuates around zero.

We then discuss the case of a strong external field. Using the same parameters as the case of $F=1$, the time evolution of the exciton probability $P_{\rm ex}(t,n)$ and exciton current $j(t)$ under a strong field of $F=70$ are presented in Figs. 9(a)-(c) and Fig. 9(d), respectively. Due to the strong external field, the exciton is found to be localized at the initial site of excitation irrespective of the exciton-phonon coupling strength. A bare exciton case is shown in Fig. 9(a) for comparison. The localization of the exciton wave packet can be rationalized in terms of Wannier-Stark states wa_59 . In the Wannier representation, the Wannier-Stark states are denoted as $|\Psi_{n}\rangle=\sum_{m}J_{m-n}(\gamma)|m\rangle$, where $J_{m-n}(\gamma)$ is the Bessel function of order $m-n$ with $\gamma=2J/Fd$ ha_04 . Due to properties of the Bessel functions, the Wannnier-Stark states extend over the interval $L\simeq 2J/F$, which leads to well localized states in the limit of a strong external field Emin1987 . As shown in Figs. 9(b) and (c), strong exciton-phonon coupling leads to the decay of BOs. The exciton currents in the presence of strong exciton-phonon coupling exhibits damped oscillations with a time scale of $8t_{B}$ ($6t_{B}$) for $g=3$ ($4$).

IV Conclusion

We have studied transient dynamics of the Holstein polaron in a one-dimensional ring under a constant external field using the Dirac-Frenkel time-dependent variational principle and the novel multi-D2 Ansatz, which is a linear combination of the usual Davydov $\rm D_{2}$ trial states. In both the diagonal and off-diagonal coupling cases, our efficient variational calculations are in perfect agreement with those obtained from the numerically exact HEOM method. Moreover, the relative deviation is found to decay with the increasing multiplicity of the multi-D2 Ansatz, which vanishes in the limit of $M\to\infty$, inferring that our approach is numerically exact in that limit.

Firstly, the influence of the initial condition is studied in the absence of exciton-phonon coupling. For an initial broad Gaussian wave packet, typical BOs are found with the center of the wave packet oscillating but its shape essentially unchanged. Starting from an initial narrow Gaussian wave packet, the exciton wave exhibits a symmetric breathing mode with its width oscillating with the Bloch period and its center of mass fixed at the original location.

The effect of the exciton-phonon coupling is the focus of our investigation. In general, weak diagonal coupling breaks the spatial periodicity while keeping the temporal periodicity of the exciton wave. For an initial broad Gaussian wave packet, the application of weak diagonal coupling modifies BOs with the exciton wave packet broadened and the exciton current reduced. For an initial narrow Gaussian wave packet, after the addition of weak coupling, the spatial symmetry of the exciton wave is broken and the center of mass of the exciton wave packet oscillates away from the original position, leading to a non-zero exciton current. In particular, a saturated current which exists in an infinite linear lattice is not found in a finite-sized ring within the anti-adiabatic regime, since the energy imported by the external field is not entirely absorbed by the lattice, leaving the steady state unreachable.

For strong exciton-phonon coupling, the variational method using the multi-D2 Ansatz is found to be highly accurate while it is prohibitively expensive for the HEOM method to tackle higher phonon excited states. The exciton wave packet is found to be localized due to either strong diagonal coupling or a strong external field. Finally, strong diagonal coupling gives rise to the decay of BOs under a strong external field, leading to damped oscillations of the exciton current.

Acknowledgments

Support from the Singapore National Research Foundation through the Competitive Research Programme (CRP) under Project No. NRF-CRP5-2009-04 is gratefully acknowledged.

Appendix A A gauge transformed Hamiltonian

The Hamiltonian without the phonon part in infinite lattice can be written in two parts: the tight-binding Hamiltonian

[TABLE]

scalar potential interaction

[TABLE]

we can introduce the vector potential $A\left(t\right)=-Ft$ and employ the gauge transformation

[TABLE]

where $x=\sum_{n}a_{n}^{\dagger}na_{n}$ is the position operator.

[TABLE]

can be satisfied, and spatial translational symmetry can be restored in $\tilde{H}$.

The form of $\tilde{H}$ is shown in the following.

[TABLE]

[TABLE]

[TABLE]

both side left multiply by $e^{iA\left(t\right)x}$

[TABLE]

using the following rule

[TABLE]

make

[TABLE]

we can obtain

[TABLE]

[TABLE]

[TABLE]

thus

[TABLE]

therefore

[TABLE]

Appendix B The Multi-${\rm D}_{2}$ trial state

The energies of the system is shown in the following,

[TABLE]

The Dirac-Frenkel variational principle results in equations of motion including

[TABLE]

and

[TABLE]

the current of the system under external field is described as,

[TABLE]

[TABLE]

[TABLE]

Appendix C Relative deviation

Assuming the trial wave function $|{\rm D}^{M}_{2}(t)\rangle=|\Psi(t)\rangle$ ($|\Psi(t)\rangle$ is the real wave function) at the time $t$, we introduce a deviation vector $\vec{\delta}(t)$ to quantify the accuracy of the variational dynamics based on the multiple Davydov trial states,

[TABLE]

where the vectors $\vec{\chi}(t)$ and $\vec{\gamma}(t)$ obey the Schrödinger equation $\vec{\chi}(t)=\partial|\Psi(t)\rangle/\partial t=\frac{1}{i}\hat{H}|\Psi(t)\rangle$ and the Dirac-Frenkel variational dynamics $\vec{\gamma}(t)=\partial|{\rm D}^{M}_{2}\rangle/\partial t$, respectively. Using the Schrödinger equation and the relationship $|\Psi(t)\rangle=|{\rm D}^{M}_{1,2}(t)\rangle$ at the moment $t$, the deviation vector $\vec{\delta}(t)$ can be calculated as

[TABLE]

Thus, the deviation from the exact Schrödinger dynamics can be indicated by the amplitude of the deviation vector $\Delta(t)=||\vec{\delta}(t)||$. In order to view the deviation in the parameter space $(W,J,g,\phi)$, a dimensionless relative deviation $\Sigma$ is calculated as

[TABLE]

where $N_{\rm err}(t)=||\vec{\chi}(t)||$ is the amplitude of the time derivative of the wave function,

[TABLE]

Appendix D Hierarchy equation of motion

The reduced density matrix element for the exciton system in the site basis can be written in the path integral form with the factorized initial condition as

[TABLE]

where $S[n]$ is the action of the exciton system, and $F[n,n^{\prime}]$ is the Feynman-Vernon influence functional given by

[TABLE]

where $\beta$ is the inverse of temperature $(\beta=1/k_{B}T)$, and we define following abbreviations

[TABLE]

with $V_{q}^{\ast}$ denoting the matrix representation of the operator $\hat{V}_{q}^{\dagger}=g\sum_{n}\hat{a}_{n}^{\dagger}\hat{a}_{n}e^{iqn}$

Taking derivative of Eq. (40),we have

[TABLE]

with the following super-operators defined as

[TABLE]

We then introduce the auxiliary operator $\rho_{m_{1\pm},m_{2\pm},\cdots,m_{N\pm}}(n,n^{\prime};t)$ by its matrix element as

[TABLE]

for non-negative integers $m_{1\pm},m_{2\pm},...,m_{N\pm}$. Note that $\hat{\rho}_{0......0}(t)=\hat{\rho}(t)$ as well as other auxiliary density matrices contain the complete information on the Liouville space wave packets. Differentiating $\rho_{m_{1\pm},m_{2\pm},...,m_{N\pm}}(n,n^{\prime};t)$ with respect to $t$, we can obtain a set of equations

[TABLE]

It should be noted that the HEOM consists of an infinite number of equations, but they can be truncated at finite number of hierarchy elements by the terminator for the practical calculation as

[TABLE]

—