Finite-temperature time-dependent variation with multiple Davydov states

TL;DR

This paper extends the Dirac-Frenkel variational approach with multiple Davydov states to finite-temperature spin-boson dynamics, demonstrating improved accuracy over single Davydov Ansatz and comparing favorably with exact methods across various temperatures.

Contribution

The authors develop a finite-temperature variational method using multiple Davydov states for the spin-boson model, enhancing accuracy especially at low temperatures.

Findings

Multiple Davydov Ansatz improves dynamics accuracy at low temperatures.

The variational approach agrees well with QUAPI at high temperatures.

The method is computationally efficient across different spectral densities.

Abstract

The Dirac-Frenkel time-dependent variational approach with Davydov Ans\"atze is a sophisticated, yet efficient technique to obtain an acuurate solution to many-body Schr\"odinger equations for energy and charge transfer dy- namics in molecular aggregates and light-harvesting complexes. We extend this variational approach to finite temperatures dynamics of the spin-boson model by adopting a Monte Carlo importance sampling method. In or- der to demonstrate the applicability of this approach, we compare real-time quantum dynamics of the spin-boson model calculated with that from numerically exact iterative quasiadiabatic propagator path integral (QUAPI) technique. The comparison shows that our variational approach with the single Davydov Ans\"atze is in excellent agreement with the QUAPI method at high temperatures, while the two differ at low temperatures. Accuracy in dynamics calculations employing a multitude of Davydov trial states is found to improve substantially over the single Davydov Ansatz, especially at low temperatures. At a moderate computational cost, our variational approach with the multiple Davydov Ansatz is shown to provide accurate spin-boson dynamics over a wide range of temperatures and bath spectral densities.

Equations

$$\hat{H}={\hat{H}}_{\mathrm{S}}+{\hat{H}}_{\mathrm{SB}}+{\hat{H}}_{\mathrm{B}},$$

$${\hat{H}}_{\mathrm{S}}=\frac{ϵ}{2}{\sigma }_z-\frac{\Delta }{2}{\sigma }_x,$$

$${\hat{H}}_{\mathrm{SB}}=\frac{{\sigma }_z}{2}\sum _l{\lambda }_l(b_l^{†}+b_l),$$

$${\hat{H}}_{\mathrm{B}}=\sum _l{\omega }_l{b}_l^{†}{b}_l.$$

$$C(t)=$$

$$=$$

$$J(\omega )=\sum _l{\lambda }_l^2\delta (\omega -{\omega }_l).$$

$$J(\omega )=2\alpha {\omega }_c^{1-s}{\omega }^s{e}^{-\omega /{\omega }_c},$$

$${\int }_0^{{\omega }_l}d\omega \Xi (\omega )=l,l=1,2,\dots ,N_b,$$

$$\Xi (\omega )=\frac{1}{{\Gamma }_{{\omega }_{\max }}}\frac{J(\omega )}{\omega },$$

$${\Gamma }_{{\omega }_{\max }}=\frac{1}{N_b}{\int }_0^{{\omega }_{\max }}d\omega \frac{J(\omega )}{\omega }.$$

$$i\frac{\partial }{\partial t}|\psi (t)⟩=\hat{H}|\psi (t)⟩.$$

$$|{\mathrm{D}}_1^M(t)⟩$$

$$|{\mathrm{D}}_2^M(t)⟩$$

$$L=\frac{i}{2}⟨{\mathrm{D}}_1^M(t)|\frac{\overleftrightarrow{\partial }}{\partial t}|{\mathrm{D}}_1^M(t)⟩-⟨{\mathrm{D}}_1^M(t)|\hat{H}|{\mathrm{D}}_1^M(t)⟩,$$

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{u}}_n^{\ast }}\right)-\frac{\partial L}{\partial u_n^{\ast }}=0,$$

$$⟨\hat{O}(t)⟩=⟨{\mathrm{D}}_{1,2}^M(t)|\hat{O}|{\mathrm{D}}_{1,2}^M(t)⟩.$$

$$⟨\hat{O}(t)⟩=\mathrm{Tr}\{\hat{O}{e}^{-i\hat{H}t}{\hat{\rho }}_{\mathrm{B}}^{\mathrm{eq}}|1⟩⟨1|e^{i\hat{H}t}\}.$$

$$\displaystyle\langle\hat{O}(t)\rangle=\pi^{-N_{b}}\int d^{2}\mbox{\boldmath$\alpha$}\langle\mbox{\boldmath$\alpha$}|\langle 1|\hat{\rho}_{\rm B}^{\rm eq}e^{i\hat{H}t}\hat{O}e^{-i\hat{H}t}|1\rangle|\mbox{\boldmath$\alpha$}\rangle,$$

$$\displaystyle\hat{\rho}_{\rm B}^{\rm eq}=\int d^{2}\mbox{\boldmath$\alpha$}\,p(\mbox{\boldmath$\alpha$};\,\beta)|\mbox{\boldmath$\alpha$}\rangle\langle\mbox{\boldmath$\alpha$}|,$$

$$p(\mbox{\boldmath$\alpha$};\,\beta)=\prod_{l}^{N_{b}}\left[\frac{e^{\beta\omega_{l}}-1}{\pi}\exp\left(-|\alpha_{l}|^{2}(e^{\beta\omega_{l}}-1)\right)\right].$$

$$⟨\hat{O}(t)⟩=$$

$$=$$

$$\displaystyle\langle\hat{O}(t)\rangle=\frac{1}{N_{s}}\sum_{i}^{N_{s}}\langle\mathrm{D}_{1,2}^{M}(t;\,\mbox{\boldmath$\alpha$}_{i})|\hat{O}|\mathrm{D}_{1,2}^{M}(t;\,\mbox{\boldmath$\alpha$}_{i})\rangle,$$

$$\displaystyle p(\mbox{\boldmath$\alpha$};\,\beta)=\prod_{l}^{N_{b}}\frac{1}{\sqrt{2\pi}\sigma_{l}}e^{-\frac{x_{l}^{2}}{2\sigma_{l}^{2}}}\frac{1}{\sqrt{2\pi}\sigma_{l}}e^{-\frac{p_{l}^{2}}{2\sigma_{l}^{2}}},$$

$$\displaystyle P_{z}(t)=\langle\sigma_{z}(t)\rangle=\frac{1}{N_{s}}\sum_{i}^{N_{s}}\langle\mathrm{D}_{1,2}^{M}(t;\,\mbox{\boldmath$\alpha$}_{i})|\sigma_{z}|\mathrm{D}_{1,2}^{M}(t;\,\mbox{\boldmath$\alpha$}_{i})\rangle.$$

Full text

I Introduction

Responsible for electronic dephasing and energy relaxation, interplay between electronic and nuclear degrees of freedom (DOFs) is a fundamental aspect of dynamical processes in condensed phases such as molecular aggregates and light-harvesting complexesnitzan2006chemical ; may2011charge . An accurate description of quantum dissipative dynamics in the condensed phases remains a challenging problem. One of main schemes for treating quantum dissipation is the reduced density matrix approach, with a focus on truncated system dynamics in the presence of a macroscopic thermal bath. The second-order cumulant time-nonlocal quantum master equation approachtanimura1989time ; ishizaki2005quantum ; ishizaki2009unified ; fujihashi2015impact ; chen2015dynamics and path integral methods such as iterative quasiadiabatic propagator path integral (QUAPI) techniquemakri1995tensor1 ; makri1995tensor2 ; nalbach2011iterative are examples of numerically exact methods for propagating the reduced density matrix. These nonperturbative and non-Markovian approaches allow for exploration of a full range of system-bath coupling and electronic coupling strengthsnalbach2011iterative ; ishizaki2009unified although it becomes extremely difficult to predict quantum dissipative dynamics at very low temperatures.tang2015extended ; chen2016inchworm Recent progress in ultrafast time resolved spectroscopy has stimulated methodological developments, and a large number of efficient, approximate reduced density matrix approaches are available.cheng2009dynamics ; ishizaki2010quantum ; huelga2013vibrations ; chenu2015coherence ; chen2015optimal ; lee2016semiclassical

In the aforementioned reduced density matrix approaches, coupled quantum dynamics of electronic and bosonic degrees of freedom (DOFs) is obtained explicitly only for electronic systems, whereas the bath DOFs are traced out in the reduced density matrix.may2011charge Alternatively, the wave-function approach such as the multiconfigurational time-dependent Hartree (MCTDH) approachmeyer1990multi ; matzkies1999accurate and its multilayer variant (ML-MCTDH)wang2001systematic ; wang2003multilayer ; manthe2008multilayer ; wang2006quantum can describe time propagation of the wave function of all DOFs explicitly, with equations of motion determined by the Dirac-Frenkel time-dependent variational principle. These approaches are shown to be powerful tools for obtaining numerically exact quantum dynamics at very low temperatures,wang2003multilayer and for including finite temperature effects with the help of statistical sampling of the bath initial conditions.wang2006quantum

One of the established approximate methods for describing time evolution of many-body wave functions is to employ the Dirac-Frenkel time-dependent variational approach with the Davydov Ansätze, vskrinjar1988classical ; forner1996accuracy ; sun2010dynamics which consist of sums of direct product states of localized electronic state and coherent states of the bath modes as proposed by Davydov in 1970s to describe a soliton motion in molecular chains.davydov1979solitons The time-dependent variational approach with these Ansätze has been widely used for describing excitation dynamics and nonlinear optical spectra of Holstein polaron and molecular aggregates.huynh2013polaron ; chorovsajev2014dynamics ; chen2015theory Recently it has been shown that superpositions of the Davydov Ansätze, also known as the multiple Davydov Ansätze, provide significant improvements in the accuracy of the time-dependent variation.zhou2014ground ; zhou2015symmetry ; zhou2015polaron ; wang2016variational ; huang2016polaron ; zhou2016fast ; deng2016dynamics ; fujihashi2016effect Increasing the multiplicity of these Ansätze leads to numerically exact dynamics in the Holstein model and the sub-Ohmic spin boson model.wang2016variational ; huang2016polaron ; zhou2016fast ; deng2016dynamics However, the usual variational approach with the Davydov Ansätze has been restricted to zero temperature case because of the wave function formalism.

In this study, we remove the restriction of zero temperature of the variational approach with the Davydov Ansätze in order to include the effect of finite temperature. For this purpose, we adopt a Monte Carlo like method, inspired by the ML-MCTDH approach.wang2006quantum The reliability of the extended variational approach with the Davydov Ansatz at finite temperature is identified through comparison with benchmark results obtained from the QUAPI approach. The rest of the paper is organized as follows. In Sec. II, we introduce the spin-boson Hamiltonian, Davydov Ansätze, and Monte Carlo importance sampling method. In Sec. III, we present a comparison of population dynamics between our theoretical approach and the numerically accurate QUAPI results at finite temperatures. Finally, Sec. IV is devoted to concluding remarks.

II Model

II.1 Hamiltonian

In this study, we consider a two-level system coupled to a single dissipative bath composed of harmonic oscillators as described by the spin-boson model.leggett1987dynamics The spin-boson model is a paradigm for studying a variety of physical and chemical phenomena such as electron transfer and exciton dynamics in condensed phase systems.nitzan2006chemical ; weiss2012quantum Thus it is suitable for a benchmarking system for investigating the validity of the newly developed quantum dynamical approach. The Hamiltonian of the spin-boson model is written as (we set $\hbar=1$)

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Here, $\sigma_{i}$ ($i=x,z$) are Pauli operators defined as $\sigma_{x}=|1\rangle\langle 2|+|2\rangle\langle 1|$ and $\sigma_{z}=|1\rangle\langle 1|-|2\rangle\langle 2|$ with $|1\rangle$ and $|2\rangle$ representing two electronic localized states. $\epsilon$ and $\Delta$ is the energy bias and the coupling constant between two electronic states, respectively. $b_{l}^{\dagger}$ ($b_{l}$) is creation (annihilation) operator for the bosonic bath mode of frequency $\omega_{l}$, and $\lambda_{l}$ is the strength of the coupling between the system and the $l$th mode.

Owing to Wick’s theorem,wick1950evaluation ; fetter2003quantum all the fluctuations and the dissipation process of $\hat{u}(t)\equiv e^{i\hat{H}_{\rm B}t}\sum_{l}\lambda_{l}(b_{l}^{\dagger}+b_{l})e^{-i\hat{H}_{\rm B}t}$ can be specified by a two-body correlation function:kubo1985statistical

[TABLE]

where $\langle\dots\rangle_{\rm B}$ denotes averaging over $\hat{\rho}_{\rm B}^{\rm eq}=e^{-\beta\hat{H}_{\rm B}}/{\rm Tr_{B}}\{e^{-\beta\hat{H}_{\rm B}}\}$ with $\beta=k_{\rm B}T$ being the inverse temperature. The bath spectral density is defined in terms of the coupling strength $\lambda_{l}$ asmay2011charge

[TABLE]

Several forms of $J(\omega)$ are employed in the literature based on model assumptions. In this study, we adopt the following spectral density form

[TABLE]

where $\alpha$ represents the strength of the coupling, and $\omega_{c}$ provides a phenomenological frequency cutoff. Those spectral densities with $s<1$, $s=1$, and $s>1$ are referred to as sub-Ohmic, Ohmic, and super-Ohmic, respectively.leggett1987dynamics

To obtain numerical solutions to the equations of motion for the variational parameters, the continuum spectral density needs to be discretized. The conventional logarithmic discretization method has been employed to investigate the sub-Ohmic spin-boson model at zero temperature,bulla2003numerical ; zhou2014ground ; zhou2015symmetry ; wang2016variational but may not be appropriate for finite temperatures. The logarithmic discretization method samples more in the low frequency domain, and can easily characterize the low frequency bath modes at zero temperature. On the other hand, the high frequency bath modes should also be of importance at finite temperatures due to their initial excitations according to the Bose statistics. For this reason, we follow the spectral density discretization procedure in the ML-MCTDH approach.wang2001systematic ; wang2003multilayer Following Ref. wang2003multilayer, , we introduce a density of frequencies $\Xi(\omega)$ defined on $[0,\omega_{\rm max}]$, where $\omega_{\rm max}$ is the upper bound of the frequency, and discretize the continuum of frequencies as

[TABLE]

where $N_{b}$ is the number of discrete bath modes and $\omega_{N_{b}}=\omega_{\max}$. The parameter $\lambda_{l}$ for each $\omega_{l}$ is then given by $\lambda_{l}=\sqrt{J(\omega_{l})/\Xi(\omega_{l})}$. The precise functional form of $\Xi(\omega)$ does not affect the final outcome if sufficient bath modes are included. For efficient numerical calculations, however, the form of $\Xi(\omega)$ can be chosen as

[TABLE]

where

[TABLE]

Here, $\Gamma_{\omega_{\max}}$ is the factor guaranteeing $\int_{0}^{\omega_{\max}}d\omega\ \Xi(\omega)=N_{b}$. For the Ohmic bath ($s=1$), $\omega_{l}$ is expressed as $\omega_{l}=-\omega_{c}\ln\left(1-l\Gamma_{\omega_{\max}}/\omega_{c}\right)$. For the sub-Ohmic bath ($s<1$), $\omega_{l}$ can only be given implicitly through Eq. (2.8).

It is noted that an arbitrary spectral density can be employed in our variational approach. However, handling the spectral density with a complex structure may require a large number of discretized bath modes or a sophisticated discretization procedure. Thus, for simplicity we restrict our numerical calculations to the functional form of Eq. (2.7) since the goal of our study is to demonstrate the validity of our finite temperature approach.

II.2 Time dependent variational approach with multiple Davydov trial states

In general, the time evolution of the wave function $|\psi(t)\rangle$ for the total Hamiltonian, $\hat{H}$, is described with the schrödinger equation,

[TABLE]

In this work, we employ two trial wave functions, namely, the multiple $\mathrm{D}_{1}$ and $\mathrm{D}_{2}$ Ansatz, to solve the above schrödinger equation. Both are known as the multiple Davydov Ansätze. An extension to the single $\mathrm{D}_{1}$ Ansatz, the multiple $\mathrm{D}_{1}$ Ansatz has been employed to investigate both static and dynamical properties of many-body quantum systems such as the spin-boson model and the Holstein molecular crystal model.zhou2015symmetry ; wang2016variational ; zhou2016fast The time-dependent version of the multi-$\mathrm{D}_{1}$ Ansatz can be written as

[TABLE]

where $\ket{0}_{\rm B}$ is the vacuum state of the boson bath. The variational parameters $A_{n}(t)$ and $B_{n}(t)$ are the amplitudes in states $\ket{1}$ and $\ket{2}$, respectively. $f_{nl}(t)$ and $g_{nl}(t)$ represent bath-mode displacements with $n$ and $l$ denoting the $n$th coherent state and the $l$th bath mode, respectively. $M$ is the multiplicity, which denotes the number of single $\mathrm{D}_{1}$ states in Eq. (LABEL:D1fun). For $M=1$, the multi-$\mathrm{D}_{1}$ Ansatz reduces to the standard Davydov $\mathrm{D}_{1}$ Ansatz. Similarly, the multi-$\mathrm{D}_{2}$ Ansatz is given byzhou2015polaron ; huang2016polaron

[TABLE]

where $f_{nl}(t)$ are the displacements of $l$th bath mode on any electronic states. Thus, the multi-$\mathrm{D}_{2}$ Ansatz is a simplified version of the multi-$\mathrm{D}_{1}$ Ansatz, since the bath displacements of the multiple $\mathrm{D}_{1}$ ($\mathrm{D}_{2}$) trial state is site-dependent (site-independent). For $M=1$, the multi-$\mathrm{D}_{2}$ Ansatz reduces to the single $\mathrm{D}_{2}$ Ansatz.

The time-dependent variational parameters, $A_{n}(t)$, $B_{n}(t)$, $f_{nl}(t)$ and $g_{nl}(t)$ are determined by adopting the Lagrangian formalism of the Dirac-Frenkel time-dependent variational principle. The Lagrangian associated with the trial state $|{\rm D}_{1}^{M}(t)\rangle$ is given by

[TABLE]

where the operator $\overset{\leftrightarrow}{\partial}/\partial t$ denotes $\overset{\rightarrow}{\partial}/\partial t-\overset{\leftarrow}{\partial}/\partial t$. The Dirac-Frenkel time-dependent variational principle yields the equations of motion for the variational parameters,

[TABLE]

where $u_{n}$ are the variational parameters, i.e., $A_{n}(t)$, $B_{n}(t)$, $f_{nl}(t)$ and $g_{nl}(t)$, and $u_{n}^{\ast}$ is the complex conjugate of $u_{n}$. Similarly, time evolution for the multi-$\mathrm{D}_{2}$ Ansatz, can be derived using the Dirac-Frenkel time-dependent variational principle. Detailed derivations of the equations of motion for the variational parameters of multiple Davydov Ansätze can be found in Ref. wang2016variational, .

The expectation value of an observable can be expressed as $\langle\hat{O}(t)\rangle=\mathrm{Tr}\{\hat{O}\hat{\rho}_{\rm tot}(t)\}$, where $\hat{\rho}_{\rm tot}(t)=e^{-i\hat{H}t}\hat{\rho}_{\rm tot}(0)e^{i\hat{H}t}$ is the total density operator, and $\hat{O}$ is a time independent operator. Using the multi-$\mathrm{D}_{1}$ or $\mathrm{D}_{2}$ Ansatz, the expectation value of the observable of interest at zero temperature can be calculated as

[TABLE]

In the next subsection, this expression is extended to the finite temperatures.

II.3 Observables at finite temperature

The conventional time-dependent variation described above is only applicable at zero temperature. In this subsection, we extend this variational approach to finite temperatures by adopting a Monte Carlo like method, similar to what is used in the ML-MCTDH approach.wang2003multilayer ; wang2006quantum The initial density matrix for the entire system is assumed to have a factorized form, i.e. $\hat{\rho}_{\rm tot}(0)=\hat{\rho}(0)\hat{\rho}_{\rm B}^{\rm eq}$, where $\hat{\rho}(0)=|1\rangle\langle 1|$. The extension to more general initial conditions with superposition of $|2\rangle$ and $|1\rangle$ is straightforward, which is important for modeling nonlinear spectroscopy.

The expectation value of an observable $\hat{O}(t)$ at finite temperatures can be expressed as

[TABLE]

In principle, the observables at $t$ can be calculated in any representations. We employ the coherent state representation to calculate the observable as

[TABLE]

where |\mbox{\boldmath\alpha}\rangle denotes a direct product of coherent states $(\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$, $\cdots$, $\alpha_{N_{b}})$ for the $N_{b}$ discrete bath modes, and is expressed as |\mbox{\boldmath\alpha}\rangle=\exp(\sum_{l}\alpha_{l}b_{l}^{\dagger}-\mathrm{{H.c.}})\ket{0}_{\rm B}. Each $\alpha_{i}$ runs over all of the feasible coherent states. The element of area $d^{2}\alpha_{i}$ on the complex plane of $\alpha_{i}$ denotes $d\mathrm{Re}(\alpha_{i})\cdot d\mathrm{Im}(\alpha_{i})$, which $\mathrm{Re}(\alpha_{i})$ and $\mathrm{Im}(\alpha_{i})$ are the real and imaginary part of $\alpha_{i}$, respectively. The equilibrium density matrix of the bath at a finite temperature is a diagonal matrix, and can be expressed asglauber1963coherent ; sudarshan1963equivalence

[TABLE]

where p(\mbox{\boldmath\alpha};\,\beta) represents the diagonal elements of the density matrix in the coherent state representation and can be expressed ashillery1984distribution

[TABLE]

As shown in Eq. (2.20), p(\mbox{\boldmath\alpha};\,\beta) is a positive defined function of $\alpha$ and can be seen as a probability density. Substituting Eq. (2.19) into Eq. (2.18), the observables $\braket{O(t)}$ at finite temperature can be obtained by the average according to the probability density p(\mbox{\boldmath\alpha};\,\beta) as

[TABLE]

For the second equality, we have use the multi-$\mathrm{D}_{1}$ or $\mathrm{D}_{2}$ Ansatz, \ket{\mathrm{D}_{1,2}^{M}(t;\,\mbox{\boldmath\alpha})}=e^{-i\hat{H}t}\ket{1}\ket{\mbox{\boldmath\alpha}}, and \ket{\mathrm{D}_{1,2}^{M}(0;\,\mbox{\boldmath\alpha})} denotes a trial state with initial bath displacements of $\alpha$ at $t=0$. For the case of the multi-$\mathrm{D}_{1}$ Ansatz, initial condition parameters are $A_{1}(0)=1$, $B_{1}(0)=0$, $A_{n}(0)=B_{n}(0)=0$ for $n\neq 1$ and $f_{nl}(0)=g_{nl}(0)=\alpha_{l}$ for all $n$ and $l$. Likewise, initial parameters of the multi-$\mathrm{D}_{2}$ Ansatz are $A_{1}(0)=1$, $B_{1}(0)=0$, $A_{n}(0)=B_{n}(0)=0$ for $n\neq 1$ and $f_{nl}(0)=\alpha_{l}$ for all $n$ and $l$.

The expectation value of the observable at finite temperature can numerically be calculated by the technique of Monte Carlo importance sampling as

[TABLE]

where $N_{s}$ is the sampling number. The configuration \mbox{\boldmath\alpha}_{i} for the bath is numerically generated according to p(\mbox{\boldmath\alpha};\,\beta) by importance sampling, where p(\mbox{\boldmath\alpha};\,\beta) is the Boltzmann distribution used as the weighting function in the importance sampling procedure. Letting $2\sigma_{l}^{2}=1/(e^{\beta\omega_{l}}-1)$ and $\alpha_{l}=x_{l}+ip_{l}$, p(\mbox{\boldmath\alpha};\,\beta) in Eq. (2.20) can be partitioned into two independent Gaussian distribution as

[TABLE]

where $\sigma_{l}$ can be taken as the variance of the Gaussian distribution. To avoid singularity, the initial displacements in the trial states is determined by setting $f_{nl}(0)=g_{nl}(0)=\alpha_{l}+\epsilon_{0}$, where noise $\epsilon_{0}$ satisfying the uniform distribution $[-10^{-2},10^{-2}]$ is added to the variational parameters of the initial states. From the definition of $\sigma_{l}$, a lower temperature or the higher frequency $\omega_{l}$ gives a smaller $\sigma_{l}$. The zero temperature case corresponds to every bath mode being in the ground state initially, and it is equivalent to a coherent state with displacement parameter, $\alpha_{l}=0$ for all $l$. In this case, the observable expression of Eq. (2.21) or Eq. (2.22) reduces to Eq. (2.16).

III Results and discussion

In this section, we compare results from our variational approach with those of the QUAPI methodmakri1995tensor1 ; makri1995tensor2 ; nalbach2011iterative in order to demonstrate the applicability of our approach to unraveling many-body dynamics. In the spin-boson model, a principal observable of interest is the population difference between the electronic states, $P_{z}(t)$, calculated as

[TABLE]

All units in the numerical calculations are scaled by the cutoff frequency $\omega_{c}$. In this section, for simplicity, we choose zero bias as a large energy bias requires a large multiplicity.borrelli2016quantum The electronic coupling constant and the coupling strength of spectral density are set as $\Delta=0.1\omega_{c}$ and $\alpha=0.05$, respectively. The number of the discrete bath modes is fixed at $N_{b}=250$. In all numerical calculations of this study, statical averages are taken over a maximum of 400 realizations, which is the number sufficient for convergence at $k_{\rm B}T=0.2\omega_{c}$. For lower temperatures the observables converge with less realizations.

Here we briefly mention the numerical convergence of the QUAPI procedure. The basic idea here is to consider a Trotter splitting of the short-time propagator for the total Hamiltonian into one part depending on the system Hamiltonian and another involving the bath and the system-bath coupling term.makri1995tensor1 ; makri1995tensor2 ; nalbach2011iterative The short-time density matrix propagator describes time evolution over a time slice $\Delta t$. This splitting is by construction exact in the limit $\Delta t\rightarrow 0$. With finite time slice $\Delta t$ in practical calculations, the numerical error can be eliminated by choosing sufficiently small $\Delta t$ to achieve convergence. On the other hand, the bath DOFs generate bath correlations in the Feynman-Vernon influence functional. For any finite temperature, these correlations decay with a dephasing time constant of $\omega_{c}^{-1}$, thus a memory time window, $\tau_{\rm mem}=k_{\rm max}\Delta t$, should be defined to handle the bath correlation truncation beyond a certain time span. Neglected beyond $\tau_{\rm mem}$, all correlations are included exactly within a finite memory time of $\tau_{\rm mem}$. To reach convergence, the memory time window should be enlarged by increasing $k_{\rm max}$ such that all memory effects are taken into account up to a desired accuracy. The QUAPI procedure fails to converge when the memory time is too long. The accuracy of the truncation of $k_{\rm max}$ in Figs. 1 and 2 were checked to make sure that the numerical results are converged. In this work, the time slice and the number of memory parameter are set to be $\Delta t=0.125$ and $k_{\rm max}=7$, respectively. The convergence failure of the QUAPI in deep sub-Ohmic regime is discussed in Subsection IIIB.

III.1 The Ohmic Regime

Figures 1(a) and 1(b) display time evolution of $P_{z}(t)$ calculated by the multi-$\rm D_{1}$ Ansätze, the multi-$\rm D_{2}$ Ansätze and the QUAPI approach in the Ohmic case ($s=1$) at two temperatures ($k_{\rm B}T/\omega_{c}=0.01$ and $k_{\rm B}T/\omega_{c}=0.2$), respectively. There is perfect agreement between the ${\rm D}_{2}^{M=3}$ Ansatz and the QUAPI approach at the low temperature $k_{\rm B}T/\omega_{c}=0.01$, as shown in the right panel of Fig. 1(a). On the other hand, the agreement between the ${\rm D}_{2}^{M=1}$ Ansatz and the QUAPI approach is unsatisfactory. As discussed in Ref. wang2016variational, , this is because time evolution of the phonon wave functions show the plane wave like behavior in the weak coupling regime and is difficult to be described by superposition of coherent states. Thus, more phonon coherent states are needed to capture the accurate dynamics. The multi-${\rm D}_{1}$ Ansatz provides accurate population dynamics by $M=2$, and the agreement between the multi-${\rm D}_{1}$ Ansatz and the QUAPI approach is good even at $M=1$. This clearly demonstrates that the multi-${\rm D}_{1}$ Ansätze is more flexible than the multi-${\rm D}_{2}$ Ansätze.

In the left panel of Fig. 1(b), dynamics obtained by both the ${\rm D}_{1}^{M=2}$ Ansatz and the QUAPI approach agree with each other at high temperature ($k_{\rm B}T/\omega_{c}=0.2$) despite small difference between the ${\rm D}_{1}^{M=2}$ Ansatz and QUAPI populations. The result by the ${\rm D}_{1}^{M=1}$ Ansatz also appears to be similar to that by the $M=2$ case. The small difference between the multi-${\rm D}_{1}$ Ansätze and the QUAPI approach may be due to insufficient number of bath modes. From Eq. (2.23) and the definition of the variance of the Gaussian distribution $\sigma_{l}$, the value of $\sigma_{l}$ becomes large at high temperature. The increase of $\sigma_{l}$ leads to large values of $f_{nl}(0)$ and $g_{nl}(0)$ even for high frequency modes which can be ignored at low temperatures. Therefore, the number $N_{b}$ of bath modes required as well as the sampling number become large due to large value of $\sigma_{l}$ at high temperatures. The accuracy of the ${\rm D}_{2}$ Ansatz with $M=1$ at high temperature is improved significantly unlike in the low temperature regime, and are similar to that of the ${\rm D}_{1}^{M=2}$ Ansatz. These results indicate that the superposition of the Davydov trial states is more fragile against thermal fluctuations at higher temperatures, and thus the trial states with small multiplicity are sufficient for description of dynamics in this regime. The increased computational cost due to additional bath modes and extended sampling is offset by the reduced Ansatz multiplicity, and thus our variational approach with importance sampling remains efficient even at high temperatures.

III.2 The Sub-Ohmic Regime

In this subsection, we focus on the sub-Ohmic regime. Figures 2(a) and 2(b) present results of $s=0.8$ at $k_{\rm B}T/\omega_{c}=0.01$ and $k_{\rm B}T/\omega_{c}=0.2$, respectively. These results demonstrate that the variational approach with the multi-${\rm D}_{1}$ Ansätze provides excellent dynamics simulation over a range of temperatures and bath spectral exponents $s$. Calculated dynamics in Fig. 2 shows the fast coherence decay compared to the Ohmic bath case, although the sub-Ohmic regime corresponds to the slow bath regime. This fast coherence decay can be explained as follows. A direct measure of the coupling strength between the system and the bath is the bath reorganization energy, $E_{r}=\int_{0}^{\infty}d\omega J(\omega)/\omega$, which represents the dissipated environment energy after electronic excitation in accordance with vertical Franck-Condon transition in the electron transfer theory.ishizaki2010quantum From the definition of spectral density of Eq. (2.7), the value of $E_{r}$ in the case of $s=0.8$ is large compared with the Ohmic bath case when other parameters except the value of $s$ are fixed. A large $E_{r}$ corresponds to larger fluctuations according to the fluctuation-dissipation relation,ishizaki2010quantum ; kubo1985statistical which eradicate electronic coherence similar to what occurs in the Ohmic case. This physical explanation is consistent with the numerical fact that the required multiplicity in both multi-${\rm D}_{1}$ and multi-${\rm D}_{2}$ Ansätze in the case of $s=0.8$ is not dissimilar to the Ohmic case.

Finally, we discuss potential applicability of the multiple Davydov Ansätze in the deep sub-Ohmic regime ($s<0.5$), which corresponds to situations with the ultra slow bath relaxation. In the slow bath regime, it is known that some notable discrepancies exist between exact results and those from the QUAPI procedure at long times, as discussed in Ref. chen2016inchworm, . To obtain correct long time dynamics, one would need to increase the memory window size $k_{\rm max}$, making QUAPI convergence more difficult to achieve. It has been demonstrated that the multi-${\rm D}_{1}$ Ansatz is consistent with results from the second-order cumulant time-nonlocal quantum master equation and the real-time path integral Monte Carlo approaches for $s=0.25$ at zero temperature, and is a reliable, efficient method for describing quantitatively accurate dynamics of the deep sub-Ohmic regime.wang2016variational From Figs. 1 and 2, the validity of the multiple Davydov Ansätze at finite temperatures also seems to be independent of the exponent $s$, but the require multiplicity may depend on it. In order to clarify the relation between the spectral exponent and the required multiplicity, we explore a case with an exponent $s$ smaller than that in Fig. 2. As a helpful reference, we also compare our results with those from the thermal field dynamics (TFD),chen2017 similar to the approach of Borrelli and Gelin.borrelli2016 Figure 3 presents result of $s=0.6$ at $k_{\rm B}T/\omega_{c}=0.01$. The dynamics by ${\rm D}_{1}^{M=2}$ Ansatz is in perfect agreement with those of QUAPI and TFD, as shown in left panel of Fig. 3. As demonstrated in the right panel of Fig. 3, it is found that the multi-${\rm D}_{2}$ Ansatz slightly overestimates the amplitude of population oscillation at $M=3$, and one needs to increase the multiplicity further to $M=4$ for numerical convergence. Fortunately, no large increase in the multiplicity of the multi-${\rm D}_{2}$ Ansatz is needed for adequate convergence with decreasing $s$. To explore how the required multiplicity is affected by thermal fluctuations, we have considered a case of relatively small system- bath coupling strength $\alpha$ in all numerical calculations of this study. From some previous studies of the multiple Davydov Ansätze,wang2016variational ; zhou2016fast it has been found that the required multiplicity is reduced by increasing the system-bath coupling strength because the large bath-induced fluctuations eradicate the superposition of single Davydov states. Especially in the mod- erate to strong coupling regimes, our results suggest that both the multi-${\rm D}_{1}$ and the multi-${\rm D}_{2}$ Ansätze hold an advantage over the QUAPI method in the sub-Ohmic regime, where our approach is not burdened with significantly increased computational cost as the spectral density exponent $s$ is reduced. Therefore, the multiple Davdov Ansätze including temperature effects is expected to be a reliable, efficient tool for exploring quantum dynamics in the sub-Ohmic regime at both low and high temperatures.

IV Concluding remarks

In this study, we have extended the Dirac-Frenkel time-dependent variational approach with the Davydov Ansätze to finite temperatures by adopting the Monte Carlo importance sampling method. To demonstrate the applicability of this approach, we have compared real-time quantum dynamics of the spin-boson model calculated by our approach with that from the numerically exact QUAPI approach. It is shown that our variational approach with the multiple Davydov Ansätze is accurate over a range of temperatures and bath spectral densities. Variational dynamics with the single Davydov Ansatz shows excellent agreement with the QUAPI results at high temperatures, while the difference between both approaches becomes significant at low temperatures. Accuracy in dynamical calculations employing a multiple Davydov trial states is improved significantly over the single Davydov Ansatz, especially in the low temperature regime. The reduction in the Ansatz multiplicity due to thermal fluctuations cancels out the computational cost increase due to an increased number of bath modes and extended sampling at high temperatures, and thus our variational technique with importance sampling remains efficient even at an elevated temperature. Our results in the sub-Ohmic regime demonstrate great advantage of the variational approach with the multiple Davydov Ansätze because conventional perturbative approaches fail to describe strongly non-Markovian dynamics due to the long-time tail of the time correlation function of the sub-Ohmic bath. Our variational approach including temperate effects is expected to open up new avenues of probing quantum dynamics in the sub-Ohmic regime.

The novel advantage of the wave function propagation methods is to give access to the dynamics of all bath DOFs explicitly.schulze2015explicit ; schulze2016multi ; sun2010dynamics ; fujihashi2016effect ; zhou2015polaron Kühn and coworker have investigated impacts of the quantum mechanically mixed electronic and vibrational states on electronic energy transfer dynamics in the Fenna-Matthews-Olson pigment-protein complex by tracking time evolution of bath DOFs based on the ML-MCTDH approach.schulze2015explicit ; schulze2016multi Their calculations clearly showed the importance of vibrational motion on the local electronic ground states in the quantum mixing of electronic and vibronic excitations, which is consistent with the argument in Ref. fujihashi2015impact, . However, the zero temperature assumption may lead to unreliable predictions for their role at physiological temperatures because the mixed electronic and vibrational states are fragile against thermal fluctuations.fujihashi2015impact Our finite-temperature time-dependent variational approach with the multiple Davydov states is demonstrated to remain efficient even at an elevated temperature, and thus can explore effects of thermal fluctuations on the mixed electronic-vibrational states by tracking dynamics of vibrational DOFs.

The time-dependent variational approach with importance sampling requires averaging over a large number of realizations at high temperatures, which may increase the computational cost in comparison with zero-temperature cases despite that the multiplicity $M$ required for convergence decreases with the increasing temperature. Developing efficient techniques of importance sampling holds the key to improved statistics. By employing the thermo field dynamics approachborrelli2016 , the variational approach with the Davydov Ansätze can be applied to finite temperature scenarios while avoiding the sampling procedure. A comprehensive study along this direction is currently in progress.chen2017

Acknowledgements.

One of us (L.W.) would like to thank Ke-Wei Sun for helpful discussion. This work was supported by the Singapore National Research Foundation through the Competitive Research Programme (CRP) under Project No. NRF-CRP5-2009-04 and the Singapore Ministry of Education of Education through the Academic Research Fund (Tier 2) under Project No. MOE2014-T2-1-099.