Polaron dynamics with off-diagonal coupling: beyond the Ehrenfest approximation

TL;DR

This paper advances the simulation of exciton-phonon dynamics in molecular crystals by employing a variational method with the multi-D2 Ansatz, surpassing the Ehrenfest approximation in accuracy and capturing detailed exciton momentum and temperature effects.

Contribution

The study introduces a multi-D2 variational approach that improves upon the Ehrenfest approximation for modeling off-diagonal exciton-phonon coupling dynamics.

Findings

Multi-D2 Ansatz yields more accurate exciton momentum distributions.

Exciton steady-state distributions depend on transfer integral and off-diagonal coupling.

The variational method is effective across temperature regimes.

Abstract

Treated traditionally by the Ehrenfest approximation, dynamics of a one-dimensional molecular crystal model with off-diagonal exciton-phonon coupling is investigated in this work using the Dirac-Frenkel time-dependent variational principle with the multi-D$_2$ {\it Ansatz}. It is shown that the Ehrenfest method is equivalent to our variational method with the single D$_2$ {\it Ansatz}, and with the multi-D$_2$ {\it Ansatz}, the accuracy of our simulated dynamics is significantly enhanced in comparison with the semi-classical Ehrenfest dynamics. The multi-D$_2$ {\it Ansatz} is able to capture numerically accurate exciton momentum probability and help clarify the relation between the exciton momentum redistribution and the exciton energy relaxation. The results demonstrate that the exciton momentum distributions in the steady state are determined by a combination of the transfer integral and the off-diagonal coupling strength, independent of the excitonic initial conditions. We also probe the effect of the transfer integral and the off-diagonal coupling on exciton transport in both real and reciprocal space representations. Finally, the variational method with importance sampling is employed to investigate temperature effects on exciton transport using the multi-$\rm D_2$ {\it Ansatz}, and it is demonstrated that the variational approach is valid in both low and high temperature regimes.

Equations

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

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

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

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

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

$${\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^{†}.$$

$${\omega }_q={\omega }_0\left[1+(\frac{2\left|q\right|}{\pi }-1)W\right],$$

$$|{\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}{\partial \dot{{\psi }_{in}^{\ast }}}\right)-\frac{\partial L}{\partial {\psi }_{in}^{\ast }}$$

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

$$L$$

$$\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}$$

$$\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 }}_{jq}}{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{ph}}|{\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)⟩.$$

$$\vec{\delta }(t)$$

$$\vec{\delta }(t)=-i\hat{H}|{\mathrm{D}}_2^M(t)⟩-\frac{\partial }{\partial t}|{\mathrm{D}}_2^M(t)⟩.$$

$$\sigma =\frac{\max \{\Delta (t)\}}{\mathrm{mean}\{N_{\mathrm{err}}(t)\}},t\in [0,t_{\max }].$$

$$N_{\mathrm{err}}(t)$$

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

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

$${\hat{a}}_k^{†}=N^{-1/2}\sum _n{e}^{-ikn}{\hat{a}}_n^{†},{\hat{a}}_n^{†}=N^{-1/2}\sum _k{e}^{ikn}{\hat{a}}_k^{†}.$$

$$c(t)$$

$$\mathrm{MSD}\text{leavevmode}(t)$$

$$c_k(t)$$

$$k-\mathrm{MSD}\text{leavevmode}(t)$$

$$H_{el}=\sum _n\left[-J+\alpha \left(u_{n+1}-u_n\right)\right]\left({\hat{a}}_n^{†}{\hat{a}}_{n+1}+H.c.\right),$$

$$H_{ph}=\frac{K}{2}\sum _n{\left(u_{n+1}-u_n\right)}^2+\frac{M}{2}\sum _n{\dot{u}}_n^2,$$

$$u_n=\sqrt{\frac{1}{2M\omega }}\left(b_n^{†}+b_n\right),$$

$$\gamma =\sqrt{\frac{1}{2M\omega }}\alpha$$

$${\hat{H}}_{ex}=-J\sum _n\left({\hat{a}}_n^{†}{\hat{a}}_{n+1}+H.c.\right),$$

$${\hat{H}}_{ph}=\sum _n{\omega }_0{\hat{b}}_n^{†}\hat{b}{}_n,$$

Full text

I Introduction

Conducting polymers (CPs) are a special class of organic materials with electronic and ionic conductivity, advanced processability and extraordinary wettability xu_2005 ; das_2012 . In $1977$, Heeger et al. reported oxidized iodine-doped polyacetylene as a forerunner of CPs shir_1977 ; heller_2015 . Various experimental strategies have been developed to produce CPs by techniques such as monomer oxidation using chemical oxidative polymerization in solution darm_2014 , electrochemical polymerization on conductive substrates darmanin_2011 , and vapor-phase polymerization im_2008 . Electrical properties of CPs can be tuned by oxidation and reduction, giving rise to rapid growth of applications. Based on their good charge transport property and high quantum efficiency of the luminescence, important utilizations of CPs are found in the large scale organic light-emitting diodes burr_1990 and electronic devices such as field-effect transistors torsi_2000 . CPs have also been used as an electrode material for supercapacitors fusa_2001 ; reddy_2014 . As a logical alternative to conventional inorganic electrode materials, a composite architecture of various CPs have been developed as a cathode for ultrafast rechargeable batteries kim_2014 . In comparison with non-conducting polymers, there are many advantages of CPs with regard to their electronic properties. CPs have also been used for other purposes ates_2012 . For example, because of easy processibility in microsturing processes schultze_2005 , CPs have been considered for a wide range of biomedical and bioengineering applications: artificial muscles otero_1998 , controlled drug release abidian_2006 , and neural recording abidian_2009 ; ravi_2010 . Surface wettability based on CPs can switch between superhydrophilicity and superoleophobicity by surface morphology control at nanoscale darm_2014 , implying usage of CPs in intelligent orthopedic and dental implants liao_2013 .

In the aforementioned applications, the efficiency of charge carrier transport and exciton transport significantly impacts the overall device performance ates_2012 . The carrier or exciton transport in CPs is well described by the Su-Schrieffer-Heeger (SSH) model in which the $\pi$ electrons are treated in a tight-binding approximation and the $\sigma$ electrons are assumed to move adiabatically with the nuclei heeger_1988 . Su et al. convincingly demonstrated that solitons play a critical role in the carrier transport doping mechanism su_1979 . Troisi et al. applied the SSH model to investigate charge carrier dynamics in crystalline organic semiconductors by solving the time-dependent Schrödinger equation for the charge wave function and using the Ehrenfest theorem for classical accelerations of nuclear positions troisi_2006 ; ehrenfest_1927 . Improvements on this semi-classical method has been made to study charge transport in organic materials in recent years hultell_2006 ; troisi_2009 ; si_2015 . Temperature dependent charge carrier mobility has also been considered troisi_2006 ; si_2015 . It is believed that for short times (comparable to the phonon period) the evolution of the system is dominated by semi-classical dynamics. The traditional Ehrenfest dynamics did not well treat the decoherence effect, which is incorporated by an instantaneous decoherence correction (IDC) approach in the framework of semi-classical method si_2015 ; dong_2016 .

Even though the semi-classical dynamics in the SSH model can capture certain features of charge transport, enormous challenges still remain to accurately describe fully quantum dynamical correlations between the electronic and vibrational subsystems guanqi_2013 . In realistic polymer chains, charge transport processes occur on the nano scale and the carriers interact with the environment including the dominant phonon degrees of freedom szymanski_2005 . The SSH model includes off-diagonal exciton-phonon coupling as a nontrivial dependence of the exciton transfer integral on lattice coordinates su_1979 ; zh_97 . Due to inherent difficulties, the off-diagonal coupling is often inadequately treated in theoretical studies. Early treatments of the off-diagonal coupling include the Munn-Silbey theory mu_85 . Recently, the Davydov D2 Ansatz zhao_12 and the multiple Davydov trial states zhou_15 have been developed to study polaron dynamics in the presence of the off-diagonal coupling. However, much awaits to be studied on the rich polaron dynamics with off-diagonal coupling with regard to exciton momentum redistribution and energy relaxations dorfner_2015 .

In this work, In order to offer an accurate description of polaron dynamics including off-diagonal coupling, the Dirac-Frenkel time-dependent variational approach with the multiple Davydov trial states will be employed. We also aim to examine the accuracy of the Ehrenfest dynamics in the SSH model. We first demonstrate that the semi-classical method and the variational method using the single D2 Ansatz are equivalent. Then we check the validity of the semi-classical method (the variational method with the single D2 Ansatz) by examining its deviations from the exact quantum dynamics. The underlying physics is revealed in the real and reciprocal space representations, including the exciton transport, the exciton momentum redistribution and the exciton energy dissipation. At the closing of the paper, we show that the fully quantum mechanical method using our multiple Davydov trial states is also applicable at the finite temperatures.

The reminder of the paper is structured as follows. In Sec. II, we present the model Hamiltonian and the variational wave function, the multi-$\rm D_{2}$ Ansatz used for describing the exciton transport. In Sec. III.1, the accuracy of the variational method using the multi-$\rm D_{2}$ Ansatz is examined by the ansatz deviation, which quantifies how faithfully the trial state follows the Schrödinger equation, and it is shown that large enhancement over that of the semi-classical method have been achieved. Numerical results of polaron dynamics by the variational method using the multi-$\rm D_{2}$ Ansatz are discussed in Sec. III.2. Impacts of the transfer integral and the off-diagonal coupling on the exciton transport are studied in Sec. III.3. Effects of temperature on polaron dynamics is investigated in Sec. III.4. Conclusions are drawn in Sec. IV.

II Methodology

II.1 Model

In presence of only off-diagonal coupling, the Hamiltonian of the one-dimensional Holstein molecular crystal model takes the form

[TABLE]

where $\hat{H}_{\rm ex},\hat{H}_{\rm ph}$ and $\hat{H}_{\rm ex-ph}^{\rm o.d.}$ denote the exciton Hamiltonian, the bath (phonon) Hamiltonian, and the off-diagonal exciton-phonon coupling Hamiltonian, respectively. In the site representation,

[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. In this work, only the anti-symmetric exciton-phonon coupling is considered in Eq. (2). In the phonon momentum space, we can rewrite $\hat{H}_{\rm ph}$ and $\hat{H}_{\rm ex-ph}^{\rm o.d.}$ as,

[TABLE]

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

[TABLE]

The parameters $J$ and $\phi$ represent the transfer integral and the off-diagonal coupling strength, respectively. A linear phonon dispersion is assumed,

[TABLE]

where $\omega_{0}$ denotes a central phonon frequency, $W$ is a constant between [math] and $1$, the bandwidth of the phonon frequency is $2W\omega_{0}$, and $q=2\pi l/N$ represents the momentum index with $l=-\frac{N}{2}+1,\ldots,\frac{N}{2}$. In the rest of the paper, $\omega_{0}$ is set to unity as the energy unit, and a dispersionless optical phonon band with $W=0$ is used.

II.2 Multiple Davydov trial states

In this work, we employ the Dirac-Frenkel variational principle to obtain quantum dynamics. We use the multiple Davydov trial states with multiplicity $M$, which are essentially $M$ copies of the corresponding single Davydov Ansatz. The multi-${\rm D}_{2}$ Ansatz has less variational variables than the multi-${\rm D}_{1}$ Ansatz when a same $M$ is used, but performs better in illuminating the polaron dynamics with the off-diagonal coupling zhou_15 . The multi-${\rm D}_{2}$ state with the multiplicity $M$, can be written as

[TABLE]

where $\psi_{in}$ and $\lambda_{iq}$ are the exciton amplitudes and the phonon displacements, respectively, $n$ is the site index of the molecular ring, and $i$ labels $i$-th $\rm D_{2}$ state in the coherent superposition. If $M=1$, the multi-${\rm D}_{2}$ Ansatz reduces to the original single Davydov $\rm D_{2}$ trial state. Equations of motion for 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, the Lagrangian $L$ is given as

[TABLE]

where the first term yields

[TABLE]

and the second term is

[TABLE]

Detailed derivations of the equations of motion for the variational parameters are given in Appendix C, together with discussions on initial conditions and numerical details.

To quantify the accuracy of the variational dynamics based on the multiple Davydov trial states, we introduce a deviation vector $\vec{\delta}(t)$ defined as

[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={-i}\hat{H}|\Psi(t)\rangle$ and the Dirac-Frenkel variational dynamics $\vec{\gamma}(t)=\partial|{\rm D}^{M}_{2}\rangle/\partial t$ in Eq. (II.2), respectively. The deviation vector $\vec{\delta}(t)$ can be calculated as

[TABLE]

Thus, the accuracy of the trial state is 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,\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]

With the wave function $|{\rm D}^{M}_{2}(t)\rangle$ obtained, the total energy $E_{\rm total}=E_{\rm ex}+E_{\rm ph}+E_{\rm ex-ph}$ is calculated, where $E_{\rm ex}=\langle{\rm D}^{M}_{2}(t)|\hat{H}_{\rm ex}|{\rm D}^{M}_{2}(t)\rangle,\leavevmode\nobreak\ E_{\rm ph}=\langle{\rm D}^{M}_{2}(t)|\hat{H}_{\rm ph}|{\rm D}^{M}_{2}(t)\rangle$, and $E_{\rm ex-ph}=\langle{\rm D}^{M}_{2}(t)|\hat{H}_{\rm ex-ph}^{\rm o.d.}|{\rm D}^{M}_{2}(t)\rangle$ (see Eq. (32)). Additionally, the exciton probability $P_{\rm ex}(t,n)$ and the exciton momentum probability $P_{\rm ex}(t,k)$ are also calculated

[TABLE]

where $\hat{a}_{k}^{{\dagger}}$ ($\hat{a}_{k}$) is the creation (annihilation) operator of the exciton with the exciton momentum $k$,

[TABLE]

We then calculate the mean square displacement ${\rm MSD}$ $(t)$ of the exciton probability $P_{\rm ex}(t,n)$ as a function of time $t$,

[TABLE]

where $c(t)$ describes the centroid motion of the exciton probability. In the exciton momentum representation, the counterpart, $k-{\rm MSD}$ $(t)$ denotes the degree of deviation of the state at the time $t$ from the initial state, as shown in the following,

[TABLE]

where $c_{k}(t)$ illustrates the centroid motion of the exciton momentum probability.

III Results and discussions

III.1 The multi-D2 Davydov Ansatz

In this subsection, dynamics of Hamiltonian (1) is described fully quantum mechanically using the multi-$\rm D_{2}$ Ansatz with sufficiently large multiplicity $M$, yielding numerically accurate quantum dynamics at zero temperature zhou_15 .

We first test the accuracy of our multi-D2 Ansatz with parameters extracted from Refs. [22; 26] (this parameter set was extensively used to study realistic models of pentacene molecules). As shown in Fig. 1, the relative deviation $\sigma$ goes to zero as the multiplicity $M$ approaches infinity. A log-log plot of ($\sigma,1/M$) (inset) indicates a power-law relationship with an exponent of $\mu=0.29(1)$, further inferring a numerically exact solution in the limit of $M\to\infty$. The largest relative deviation $\sigma$ is found for the single D2 Ansatz. As presented in Appendix A, the SSH Hamiltonian is equivalent to the Holstein Hamiltonian with off-diagonal coupling. The equivalence between the semiclassical method and the variational method using the single $\rm D_{2}$ Ansatz is shown in Appendix B. Therefore, this implies that the accuracy of the semi-classical Ehrenfest dynamics can be quantified by the relative deviation of the single D2 Ansatz. The variational method with sufficiently large $M$ fully takes into account the quantum effects, yielding a much more accurate result than that with the single D2 Ansatz, which is equivalent to the semi-classical method. For example, $\sigma$ of the D${}_{2}^{M=16}$ Ansatz in Fig. 1 is smaller than $0.1$, thus the multiplicity of $M=16$ is employed to explore polaron dynamics in following subsections, unless otherwise specified.

In order to further compare the performance of our variational method using the multi-D2 Ansatz and that of the semi-classical method, the exciton movement is studied by calculating the mean square displacement ${\rm MSD}$ $(t)$. As shown in Fig. 2, the amplitudes of ${\rm MSD}$ $(t)$ from the fully quantum variational method using the $\rm D_{2}^{M=8}$, the $\rm D_{2}^{M=16}$, and the $\rm D_{2}^{M=22}$ Ansatz are smaller than that from the semi-classical Ehrenfest method (equivalent to the single $\rm D_{2}$ Ansatz), and ${\rm MSD}$ $(t)$ shows apparent convergence as the multiplicity $M$ exceeds $16$. This result is in agreement with that by the IDC approach, which the carrier is found to be less mobile in comparison with that of original Ehrenfest method yao_2012 ; dong_2016 . In this case, the transfer integral is much larger than the exciton-phonon coupling and makes more contribution to the movement of the wave front in the carrier propagation. Consequently, the exciton-phonon coupling leads to localization of the wave front. The Ehrenfest method treats the phonons semi-classically and underestimate the confinement effect of the exciton-phonon coupling on the wave function. Therefore the reduction of the mobility is attributed to the quantum mechanical description of the phonons and the electron-phonon coupling. We note, however, the change of ${\rm MSD}$ $(t)$ depends on parameter regimes. In some other cases (e.g., $J=0$ and $\phi=1.0$), phonon assisted transport dominates the exciton movement as discussed in Ref. [33].

III.2 Polaron dynamics in exciton momentum representation

In this subsection, we explore impacts of off-diagonal coupling on the exciton movement in the exciton momentum representation by using the multi-D2 Ansatz. Without the exciton-phonon coupling, the Hamiltonian of the bare exciton can be described by the first term of Hamiltonian (1), $\hat{H}_{\rm ex}$, and the energy band is $E(k)=-2J\cos(k)$. The exciton energy and the exciton momentum are constants of motion. However, in the presence of the exciton-phonon coupling, the exciton momentum may move away from initial values and the exciton energy would be dissipated.

The left and the right columns of Figs. 3(a)-(d) present the time evolution of the exciton momentum probability $P_{\rm ex}(t,k)$, starting from initial conditions $P_{\rm ex}(0,k)=1$ at $k=0$ and $\pi$, respectively. $P_{\rm ex}(t,k)$ redistributes toward a quasi stationary state, where no net energy transfer takes place between the exciton and the phonons, as shown in Figs. 3(e) and (f). Notwithstanding the difference of initial excitonic conditions, $P_{\rm ex}(t,k)$ in the case of Figs. 3(a) and (b) still relaxes to the same stationary regime, where final $P_{\rm ex}(t,k)$ is centered around $k=0$. As for $J=-0.6$, $P_{\rm ex}(t,k)$ is centered around $k=\pm\pi$, as shown in Figs. 3(c) and (d). Moreover, the exciton momentum pattern in Fig. 3(a) is shifted by $\pi$ comparing to that in Fig. 3(d) because the Brillouin zone of the former is $\pi$ shifted from that of the latter, and the shift also occurs for $P_{\rm ex}(t,k)$ in Figs. 3(b) and (c).

The energy relaxation process is known to be accompanied with a redistribution of the exciton momentum probability dorfner_2015 . With regard to the Holstein model with diagonal coupling, the exciton kinetic energy is transferred into the phonons, ending up with a constant value of $E_{\rm ex}+E_{\rm ex-ph}$ janez_2012 . However, energy relaxation in the Holstein model with off-diagonal coupling is still not well understood. In order to clarify this issue, we consider time evolution of the exciton and the phonon energy. Figs. 3(e) and (f) show the time evolution of energies in the case of $J=0.6$ and $\phi=1$ for $P_{\rm ex}(0,k=0)=1$ and $P_{\rm ex}(0,k=\pi)=1$, respectively. Under this parameter set, the initial exciton energy of $P_{\rm ex}(0,k=0)=1$ is the lowest, that of $P_{\rm ex}(0,k=\pi)=1$ is the highest, and those of other initial conditions fall in between. After the transfer integral is changed to $J=-0.6$, due to a phase shift of the Brillouin zone in the exciton momentum space, the initial exciton energy of $P_{\rm ex}(0,k=0)=1$ becomes the maximal while that of $P_{\rm ex}(0,k=\pi)=1$ turns into the minimal for all initial excitonic conditions. As a result, identical energy relaxation processes occur despite that transfer integrals have opposite signs. Thus, only the case of $J=0.6$ and $\phi=1$ is displayed for simplicity. At $t=0$, the phonons are in their vacuum states. Later, the incident exciton wave fronts generate phonons via the exciton-phonon coupling. As a consequence, the exciton energy is transferred to the phonon degrees of freedom. After a fast relaxation process, both the energies of the exciton and the phonons reach steady values. $E_{\rm ex}+E_{\rm ex-ph}$ in the steady state settles around $-2|J|$, which corresponds to the energy minimum of the exciton in the absence of the exciton-phonon coupling. In order to identify the energy contribution of each exciton momentum, we also investigate $E_{\rm ex}+E_{\rm ex-ph}$ in the exciton momentum representation. As plotted in Figs. 3(g) and (h), the initial $E_{\rm ex}+E_{\rm ex-ph}$ is $-1.2$ and $1.2$, respectively. After relaxation, the momentum of $k=0$ becomes the main contributor of $E_{\rm ex}+E_{\rm ex-ph}$ for both cases, and also determines the locations of the quasi stationary regime after the exciton momentum redistribution.

Fig. 4 presents the time evolution of the exciton momentum probability in the absence of the transfer integral. We set the initial excitonic conditions $P_{\rm ex}(0,k)=1$ of $k=0$ and $\pi$ in the left ((a),(c) and (e)) and the right ((b),(d) and (f)) column of Fig. 4, respectively. Akin to the cases of $J=0.6$ and $\phi=1$ in Fig. 3, by comparing $P_{\rm ex}(t,k)$ with two types of initial conditions, it is found that the exciton momentum probabilities redistribute and become centered around the same regimes, as shown in Figs. 4(a) and (b). Even in the absence of the transfer integral, the exciton can still be transported by the off-diagonal coupling. Figs. 4(c) and (d) plot the time evolution of the phonon energy and the exciton-phonon interaction energy. As also shown in Figs. 4(c) and (d), for $0<t\leqslant t_{0}$, the amplitudes of both $E_{\rm ph}$ and $E_{\rm ex-ph}$ reach their peaks and fluctuate until the exciton and the phonons cease to exchange energy at $t=10t_{0}$. The energy relaxation process only involves $E_{\rm ex-ph}$ because $E_{\rm ex}$ is always zero. As presented in Fig. 4(e), $E_{\rm ex-ph}$ of each exciton momentum undergoes three stages during the energy relaxation process. During $0<t\leqslant t_{0}$, they all show strong oscillations with largest amplitudes. At the intermediate stage of $t_{0}<t\leqslant 10t_{0}$, the energies of $k=\pm\pi/2$ compete with that of $k=0$. For $t>10t_{0}$, the contribution of the energy of $k=0$ to $E_{\rm ex-ph}$ reduces to almost zero, leaving the energy of $k=\pm\pi/2$ to be the dominant energy contributor. As for the case of the initial condition $P_{\rm ex}(0,k=\pi)=1$ as shown in Fig. 4(f), the competition at the second stage of $t_{0}<t\leqslant 10t_{0}$ occurs between the energies of $k=\pm\pi$ and $\pm\pi/2$ instead, and the energy of $k=\pm\pi/2$ also turns out to be the prominent contributor to $E_{\rm ex-ph}$. Consequently, the exciton momentum probability finally becomes centered around $k=\pm\pi/2$ as shown in Figs. 4(a) and (b).

III.3 Effect of transfer integral and off-diagonal coupling on exciton transport

In this subsection, we investigate the influence of the transfer integral and the off-diagonal coupling on the exciton transport of Hamiltonian (1).

By tuning the transfer integral, contributions of the transfer integral and the off-diagonal coupling on the exciton movement are examined, as shown in Fig. 5. It can be shown in the site representation that the off-diagonal exciton-phonon coupling plays a crucial role in polaron transport tamura_2012_85 ; tamura_2012_86 . As shown in Fig. 5(b), the off-diagonal coupling is the only agent for exciton movement in the absence of the transfer integral, also known as phonon-assisted transport zhao_1994 . When both the off-diagonal coupling and the direct, phonon free exchange transfer are present, because of the competition between them, the extion transport may be inhibited, as shown in Fig. 5(a). The self-trapping phenomenon is expected due to the competition between the off-diagonal coupling and the transfer integral when the energy bands is flattened at the Brollouin zone center zh_97 . In this work, the Toyozawa Ansatz is adopted to study the ground state energy bands of the Holstein model using the variational method. As presented in Fig. 5(i), the lowest energy band of $J=0.5$ and $\phi=1.0$ meets the self-trapping condition, and we can thereby take this case as an example to study the self-trapped exciton from the perspective of dynamics. In agreement with our expectation, $P_{\rm ex}(t,n)$ turns out to be localized in Fig. 5(a). By directly flipping the sign of the transfer integral to $J=-0.5$, the exciton wave fronts are found to move considerably, as shown in Fig. 5(c). Via ${\rm MSD}\leavevmode\nobreak\ (t)$ as defined in Eq. (17), the expansion of the exciton wave packets is further investigated for $J=-0.6$, $-0.5$, [math], $0.5$ and $0.6$. As plotted in Fig. 5(g), the amplitude of ${\rm MSD}\leavevmode\nobreak\ (t)$ for $J=0$ and $\phi=1$ is smaller than those of other cases with non-zero transfer integrals, except the self-trapped case of $J=0.5$ and $\phi=1.0$.

In the crystal momentum representation, the underlying physics of the ground states can be elucidated, where the crystal momentum is denoted as $K$ (see Eq. (38)). The Toyozawa Ansatz is a time independent translationally invariant trial state, viewed as a superposition of the replicas of the D2 Ansatz displaced to every lattice site, weighed by a phase factor of the total momentum zh_97 . We analyzed the energy bands of the ground states obtained from the Toyozawa Ansatz (see Appendix D). In the the off-diagonal coupling only case ($J=0$), the minima of the band are located at $K=\pm\pi/2$. The addition of positive (negative) transfer integrals moves the minima towards the center (boundary). In particular, as mentioned above, the case of $J=0.5$ flattens the band at the center of the Brillouin zone, leading to the largest effective mass of all studied cases, in accord with the self-trapping of $P_{\rm ex}(t,n)$ in Fig. 5(a).

The effect of the transfer integrals on the exciton movement in the presence of the off-diagonal coupling is further examined in the exciton momentum representation in Figs. 5(d)-(f) and (h). The exciton is created in the profile of $\left(2+\cos k\right)/2N$ in the momentum space as we excite two nearest neighboring sites initially (see Appendix C). In the subsequent relaxation process, $P_{\rm ex}(t,k)$ redistributes and becomes localized in a quasi stationary region, and the mean square displacement of the exciton momentum $k-{\rm MSD}\leavevmode\nobreak\ (t)$ approaches a plateau, as shown in Fig. 5(h). After the relaxation process, the final $P_{\rm ex}(t,k)$ is found to be determined by a combination of the transfer integral and the off-diagonal coupling strength. For the off-diagonal coupling only case, $P_{\rm ex}(t,k)$ progressively becomes localized around $k=\pm\pi/2$ (Fig. 5(e)). In the case of $J=0.5$ and $0.6$, $P_{\rm ex}(t,k)$ aggregates toward $k=0$, as seen in Fig. 5(d). Similarly, $P_{\rm ex}(t,k)$ of both $J=-0.5$ and $-0.6$ correspond to $\pm\pi$ in Fig. 5(f). In addition, $k-{\rm MSD}\leavevmode\nobreak\ (t)$ for the extreme cases of $P_{\rm ex}(t,k)=\delta_{k,0}$ and $\delta_{k,\pm\pi}$ are [math] and $2\pi^{2}$, respectively. As shown in Fig. 5(h), $k-{\rm MSD}\leavevmode\nobreak\ (t)$ of $J=0.6$ is closer to the analytical value of [math] than that of $J=0.5$, indicating that $P_{\rm ex}(t,k)$ of $J=0.6$ is more localized around the zone center than that of $J=0.5$. Likewise, $k-{\rm MSD}\leavevmode\nobreak\ (t)$ of $J=-0.6$ is nearer to the limited value of $2\pi^{2}$ than that of $J=-0.5$, illustrating that $P_{\rm ex}(t,k)$ of the former case is more localized around $k=\pm\pi$.

In the site representation, the off-diagonal coupling is known to play the role of assisting the transport of the exciton tamura_2012_86 ; zhao_1994 . In Fig. 6(a), in the absence of transfer integral ($J=0$), we study the dependence of ${\rm MSD}\leavevmode\nobreak\ (t)$ on the off-diagonal coupling strength. It is found that the exciton propagation is facilitated by the off-diagonal coupling, as shown by the site-space ${\rm MSD}\leavevmode\nobreak\ (t)$ in Fig. 6(a). However, the off-diagonal coupling can be simultaneously an agent for exciton localization. The localization effect of $\phi$ gradually increases with the coupling strength $\phi$ if $\phi$ is greater than a critical value $\phi_{c}$ zhao_12 . As shown in Fig. 6(a), the amplitude of ${\rm MSD}\leavevmode\nobreak\ (t)$ decreases with the off-diagonal coupling strength $\phi$ for $\phi>\phi_{c}=1.0$.

In the exciton momentum representation, $P_{\rm ex}(t,k)$ all ends up around $k=\pm\pi/2$ for a variety of off-diagonal coupling strengths, and the corresponding $k-{\rm MSD}\leavevmode\nobreak\ (t)$ approaches the same narrow regime around $0.25\pi^{2}$, which is the theoretical value of $k-{\rm MSD}\leavevmode\nobreak\ (t)$ for $P_{\rm ex}(t,k)=(\delta_{k,\pi/2}+\delta_{k,-\pi/2})/2$, as shown in Fig. 6(b). However, the relaxation time diverges due to the variance of the off-diagonal coupling. The time for the exciton momentum to reach the stationary regime is inversely related to the off-diagonal coupling strength, because the first stage of the time evolution (t$<0.5$t0) is accompanied by the fast exciton movement in the case of large off-diagonal coupling as presented in Fig. 6. In addition, the energy bands of various off-diagonal coupling strengths imply that the band width and effective mass are largest for $\phi_{c}^{\rm static}=1.4$ and get smaller as $\phi$ moves away from $\phi_{c}^{\rm static}$ zhao_12 . The localization feature is found both in static and dynamic calculations although the value of $\phi_{c}$ differs slightly. The off-diagonal exciton-phonon coupling leads to exciton energy dissipation and redistribution of exciton momentum in three typical scenarios corresponding to completely distinguishable band structures (this conclusion is independent of the system size), which may be formed due to a variety of compositions and geometrical structures of the organic materials, defects, doping mechanisms and deformations of CPs bredas_1984 ; kuklja_2007 .

III.4 Temperature effects

In this subsection, we extend the work to study the effect of finite temperatures on polaron dynamics. The conductivity of polymers has been measured by many workers as a function of temperature heeger_1988 ; ishiguro_1992 ; chen_1994 . The temperature effects have been in contention from a theoretical point of view. For example, Cruzeiro et al. claims that Davydov soliton is stable at $T=310$ K cruz_1988 . Later, a quantum Monte Carlo treatment has shown that the Davydov soliton is unstable above $7$ K wang_1989 . In this work, several approaches are used to study the temperature effects: a variational method with importance sampling (see Appendix E.2), the hierarchical equations of motion (HEOM) method Tanimura1 ; ch_15 , and the averaged Hamiltonian method (see Appendix E.1). The variational method with importance sampling developed by Wang et al. simulates thermal fluctuation of phonon modes by sampling the initial phonon displacements based on the Bose distribution, and thus it can deal with Holstein polaron dynamics at both low and high temperatures wanglu_un . The HEOM method is numerically exact and is capable to treat any finite temperature, serving as a benchmark here. However, the HEOM method is also numerically expensive and thus impractical when the system size is large. The variational approach can treat large systems once a proper trial wave function is adopted. In order to compare to previous attempts in the literature, the averaged Hamiltonian method has also been used, and we found that this approach is not even suitable for the spin-boson model (i.e..,  N=$2$) as shown in Appendix E.

Fig. 7 shows polaron dynamics calculated by the multi-D2 Ansatz with importance sampling and the HEOM method for two temperatures. The calculations are carried out for $J=0.8$ and $\phi=0.3$ in a ring of $N=8$ sites. $P_{\rm ex}(t,n)$ outputs obtained from the D${}_{2}^{M=16}$ Ansatz and the HEOM method at $T=0$ are shown and compared in Figs. 7(a) and (b), respectively. As revealed in Fig. 7(c), $\Delta P_{\rm ex}(t,n)$, i.e., the difference between the two methods, are two orders of magnitude smaller than the value of $P_{\rm ex}(t,n)$, indicating that variational method can be numerically exact at low temperatures with sufficient multiplicity $M$ of the multi-$\rm D_{2}$ Ansatz. The phonon displacement $\lambda_{i,q}(t=0)$ is set to zero at $T=0$, while importance sampling is used at $T=2/k_{B}$ ($\beta=0.5$) to simulate the finite temperature effects with the result displayed in Fig. 7(d). Similarly, as shown in Fig. 7(f), differences between results from the two methods are two orders of magnitude smaller than the value of $P_{\rm ex}(t,n)$ in Figs. 7(d) end (e), inferring that the variational method with importance sampling provides numerically exact results at high temperatures with sufficiently large $M$.

Next, we investigate the influence of thermal fluctuations on exciton transport. For both low and high temperatures, the exciton wave fronts depart from the site of exciton creation and propagate in opposing directions until they meet at the opposite side of the ring. During the time evolution, distinct features observed at zero temperature (Figs. 7(a) and (b)) are now significantly smeared due to the thermal fluctuations (Figs. 7(d) and (e)). As shown in Figs. 7(d) and (e), during $0.2t_{0}<t<t_{0}$ the exciton probability is more centered around the site of creation than those of Figs. 7(a) and (b). For $t\geqslant t_{0}$, the bright spots shown in Figs. 7(a) and (b) are significantly quenched in Figs. 7(d) and (e). These results indicate that the exciton transport is weakened when the temperature is increased, in line with Refs. [ 22; 49].

IV Conclusion

In this work, we have studied the dynamics of the Holstein molecular crystal model with off-diagonal coupling 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 state from the soliton literature. Traditionally used to simulate such dynamics is the semi-classical Ehrenfest method, which has been shown to be equivalent to our time-dependent variational method with the single D2 Ansatz. Calculation of the relative deviation, which quantifies the Anstaz accuracy, demonstrates that the variational method with the multi-D2 Ansatz presents much more accurate results than the semi-classical Ehrenfest dynamics. With a sufficiently large multiplicity, our variational method using the multi-D2 Ansatz can offer numerically exact solutions. We further compare $\rm MSD$ $(t)$ obtained from the semi-classical method and our variational method with that from the multi-D2 Ansatz, and find that the mobility is overestimated by the semi-classical method. These results indicate that the description beyond the semi-classical method is essential to quantitatively capture the dynamics of the SSH model.

Secondly, we have explored the underlying physics from the accurate dynamics data for the Holstein model with the off-diagonal coupling. The energy and the momentum of the bare exciton are constants of motion. However, in the presence of the exciton-phonon coupling, the exciton momentum probability is found to redistribute and become centered in stationary regions. We reveal that the momentum redistribution is only determined by the combination of the transfer integral and the off-diagonal strength, and is independent of the initial excitonic conditions used. In addition, in order to study the competition between the transfer integral and the off-diagonal coupling, we investigate the exciton transport within the exciton site and the exciton momentum representation, and the crystal momentum representation. The results show that the combination of the transfer integral and the off-diagonal coupling do not necessarily play a role in enhancing the exciton transport. Moreover, the off-diagonal coupling is demonstrated to be the simultaneous agent of transport and localization in dynamical calculations.

Lastly, the temperature effects are studied using the variational method with importance sampling by employing the multi-D2 Ansatz. In both the low and high temperature regimes, the time evolution of the exciton probability calculated from the variational method with importance sampling agrees well with that from the numerically exact HEOM method, and can be obtained much more efficiently. The results at the finite temperatures show that fast delocalization of the exciton wave is quenched due to the thermal fluctuations, indicating the weakening of the exciton transport by increasing the temperature.

Acknowledgments

We thank Yuta Fujihashi and Zheng Fulu for helpful discussion. 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 Fully quantum description of the semiclassical Hamiltonian

The semi-classical Hamiltonian is composed of the electronic and the phonon part $H=H_{el}+H_{ph}$, the electronic part is

[TABLE]

where $J$ is the transfer integral, $\alpha$ the electron-phonon coupling constant, $u_{n}$ the displacement of phonon on the $n$-th site, and $\hat{a}{}_{n}^{\dagger}\left(\hat{a}_{n}\right)$ the creation (annihilation) operators of electron. The phonon part is

[TABLE]

in which, $K$ denotes the force constant originating from the $\sigma$ bond between carbon atoms, $K=M\omega_{0}^{2}$ and $M$ the total mass of a CH-unit for trans-polyacetylene. The combination of the two parts above is identical to Su-Schrieffer-Heeger (SSH) model used for conductive polymers heeger_1988 .

Using the quantum mechanical creation and annihilation operators to describe the displacement of the phonon bath,

[TABLE]

we get $\hat{H}=\hat{H}_{ex}+\hat{H}_{ph}+\hat{H}_{int}$, with the electronic part

[TABLE]

the phonon part

[TABLE]

and the electron-phonon interaction part

[TABLE]

Fourier transforming the phonon operators into momentum space,

[TABLE]

we get

[TABLE]

just being the off-diagonal Holstein polaron model.

Appendix B Comparison between the variational method using the single $D_{2}$ Ansatz and the semi-classical method

In this part, it is shown that the dynamics obtained from the semi-classical method and the variational method using only the single $D_{2}$ Ansatz are equivalent for the spin-boson model ($N=2$).

B.1 The variational method

The Hamiltonian of the spin-boson model is

[TABLE]

where $\epsilon$ and $V$ are the spin bias and the tunneling constant, respectively. $\lambda$ ($\phi$) is the diagonal(off-diagonal) coupling strength. $\sigma_{x}$ and $\sigma_{z}$ are THE Pauli matrices, $\hat{a}^{\dagger}(\hat{a})$ is the boson creation (annihilation) operator for the phonon of frequency $\omega_{0}$.

Using the variational principle and the $D_{2}$ Ansatz, $\left|D_{2}\right\rangle\\ =A(t)\left|+\right\rangle e^{\left[f(t)\hat{a}^{\dagger}-f^{*}(t)\hat{a}\right]}\left|0\right\rangle+B(t)\left|-\right\rangle e^{\left[f(t)\hat{a}^{\dagger}-f^{*}(t)\hat{a}\right]}\left|0\right\rangle$ the equations of motion can be obtained,

[TABLE]

B.2 The semi-classical method

The semi-classical Hamiltonian can also be written as

[TABLE]

the electronic state is described by the wave function $\left|\psi\right\rangle=A\left(t\right)\left|+\right\rangle+B\left(t\right)\left|-\right\rangle$, $m$ is the effective mass of the phonon. The equations of motion from the semi-classical formalism are

[TABLE]

B.3 Comparison

We now compare the equations of motion from the variational method and the semi-classical method. From Eq. (A$7$) of Ref. [50], we get $f=\frac{1}{\sqrt{2m\omega}}\left(m\omega x+ip\right)$ and $f+f^{*}=\sqrt{{2m\omega}}x$. After we put this into the last equation of Eq. (28), we get

[TABLE]

The real part of Eq. (31) agrees with the fourth equation of Eq. (30), and imaginary part of Eq. (31) is equal to the third equation of Eq. (30), proving the equivalence of the semi-classical method and variational method using only the single $D_{2}$ Ansatz.

In conclusion, the expectation value of position $x$ and momentum $p$ obtianed from the semi-classical method agree with the ones from the variational method using THE single $D_{2}$ Ansatz.

Appendix C The multi-${\rm D}_{2}$ trial states

The energies of the system are given by the following equations,

[TABLE]

where the Debye-Waller factor is formulated as

[TABLE]

In addition, the energies can be converted to the exciton momentum representation by using

[TABLE]

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

[TABLE]

and

[TABLE]

It should be noted that the main results of this work are calculated from the above equations of motion. The equations of motion are solved numerically by means of the fourth-order Runge-Kutta method. The exciton initially sits on two nearest-neighboring sites, i.e  $\psi_{n}=(\delta_{n,N/2}+\delta_{n,N/2+1})/\sqrt{2}$. At $T=0$, the phonon state is at the vacuum state, i.e  $\lambda_{iq}=0$. In order to avoid singularity, the uniformly distributed noise $[-10^{-5},10^{-5}]$ is added to the initial variational parameters $\psi_{in}$ and $\lambda_{iq}$ at $t=0$. More than one hundred samples are averaged to get rid of the influence of the noise and reach convergent results in the simulations.

Appendix D The Toyozawa Ansatz

Our interest here includes the polaron ground-state energy band, computed as

[TABLE]

where$|\Psi(K)\rangle$ is an appropriately normalized, delocalized trial state, and $\hat{H}$ is the system Hamiltonian. The joint crystal momentum is indicated by the $K$. It should be noted that the crystal momentum operator commutes with the system Hamiltonian, and energy eigenstates are also eigenfunctions of the crystal momentum. Therefore, variations for distinct $K$ are independent. The set of $E(K)$ constitutes a variational estimate (an upper bound) for the polaron energy band. The relaxation iteration technique, viewed as an efficient method for identifying energy minima of a complex variational system, is adopted to obtain numerical solutions to a set of self-consistency equations derived from the variational principle. To achieve efficient and stable iterations toward the variational ground state, one may take advantage of the continuity of the ground state with respect to small changes in system parameters over most of the phase diagram and may initialize the iteration using a reliable ground state already determined at some nearby points in parameter space. Starting from those limits where exact solutions can be obtained analytically and executing a sequence of variations along well-chosen paths through the parameter space using solutions from one step to initialize the next, the whole parameter space can be explored.

As the D2 Ansatz is a localized state from the soliton theory. It can be delocalized into the Toyozawa Ansatz, which is the Bloch state with the designated crystal momentum, via a projection operator $\hat{P}_{K}$.

[TABLE]

where

[TABLE]

After the delocalization onto the usual $\rm D_{2}$ Ansatz, the Toyozawa Ansatz is given by

[TABLE]

[TABLE]

where $\psi_{n_{1}-n}^{K^{\prime}}$ is the exciton amplitude and $\lambda_{n_{2}-n}^{K^{\prime}}$ is the phonon displacement.

Appendix E Alternative approaches to temperature effects

We aim to investigate the effect of the temperature by comparing following approaches: the averaged Hamiltonian (see Appendix E.1), the variational method with importance sampling (see Appendix E.2), and the numerically exact HEOM method. The spin-boson model, i.e,  a Holstein model with N$=2$, is taken as the simplest example.

The variational method with importance sampling is simulated by initially employing random number generators to investigate the temperature effects using the multiple Davydov trial wave states. The influence of the temperature on the dynamical behavior is also studied using the method of averaged Hamiltonian with the Davydov D1 Ansatz developed here (see Eq.  (E.1)) and the HEOM method. Population difference $P_{z}(t)$ (see Eq. (E.1)) obtained from the $\rm D_{2}^{M=2}$ Ansatz, the $\rm D_{2}^{M=4}$ Ansatz, and the other two methods are plotted in Fig. 8 using parameters $V=-0.05,\lambda=0.5,\beta=0.5$. Unfortunataley, the more complex D1 Ansatz does not show an improvement over the multi-D2 Ansatz and the HEOM method, at finite $T$ the distinct damping out of the oscillations is not observed, in contrast to the case of the variational method with importance sampling and the HEOM method. Moreover, with more D2 states used, the results for the variational method with importance sampling come closer to those of the HEOM method.

E.1 Averaged Hamiltonian

According to the papers by Cruzeiro et al. cruz_1988 and by Förner forner_1991 , which are based on earlier work by Davydov and coworkers, temperature effects can be taken into account approximately, by using a generalized Davydov-Ansatz

[TABLE]

with the normalized excited states

[TABLE]

In this way, one can view the ’thermally averaged state’ as a linear combination of all states with a fixed phonon distribution in the lattice, where the weight factors of the individual states follow Bose-Einstein statistics. Then we can get a thermally averaged Hamiltonian

[TABLE]

with

[TABLE]

where $\beta=1/kT$ is proportional to the inverse temperature and

[TABLE]

in which $\hat{H}$ is given in Eq. (27). Thus,

[TABLE]

where average phonon number is $\bar{\nu}={1}/{(e^{\beta\omega}-1)}$ and

[TABLE]

From the Dirac-Frenkel variational principle, we get the equations of motion

[TABLE]

In the spin-boson model, physical variables of interest are

[TABLE]

Here $P_{z}(t)$ describes the population difference. With the above trial state, we obtain

[TABLE]

E.2 Variational method with importance sampling

The variational method with importance sampling is used to obtain the dynamics of the Holstein model with the off-diagonal coupling, where initial phonon displacements are chosen according to the Bose distribution. Using only two sites for simplicity, the Holstein model with the off-diagonal coupling can be reduced to a spin-boson Hamiltonian (Eq. (27)). We solve the dynamics by variational method using the multi-D2 Ansatz

[TABLE]

The temperature effects are included by considering the initial displacements based on the Bose distribution hillery_1984 . The initial bath can be expressed as

[TABLE]

where $\left|\alpha\right\rangle\equiv e^{\alpha\hat{a}-\alpha^{*}\hat{a}}\left|0\right\rangle$ and the distribution $P(\alpha)$ is

[TABLE]

it is shown to be a well behaved Gaussian function and has no singularity. Numerically, let $2\sigma^{2}=1/(e^{\beta\omega}-1)$ and $\alpha=x+iy$,

[TABLE]

Then, we can generate the configuration $\alpha$ for the bath according to $P(\alpha)$ by Monte Carlo method. The initial displacements in the trial states is determined by setting $f_{i}\left(t=0\right)=\alpha+\epsilon_{0}$, where a small noise $\epsilon_{0}\in\left[-10^{-2},10^{-2}\right]$ is added to increase the numerical stability. According to the equations of motion obtained from the Dirac-Frenkel variational principle, the dynamics of the system can be obtained. The final result is averaged over enough realizations (more than $50000$) to ensure the convergence of relevant physical quantities. In the same way, initial displacements are also chosen according to the temperature in the fully quantum description of the SSH model.

Following are the corresponded equations of motion

[TABLE]

[TABLE]

and

[TABLE]

where $S_{ki}=e^{-\frac{1}{2}(f_{k}^{*}f_{k}+f_{i}^{*}f_{i})+f_{k}^{*}f_{i}}$.

—