Finite-temperature time-dependent variation with multiple Davydov states

# Finite-temperature time-dependent variation with multiple Davydov states

## Abstract

The Dirac-Frenkel time-dependent variational approach with Davydov Ansätze is a sophisticated, yet efficient technique to obtain an accurate solution to many-body Schrödinger equations for energy and charge transfer dynamics 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 order 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ätze 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.

## 1Introduction

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 complexes[1]. 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 approach[3] and path integral methods such as iterative quasiadiabatic propagator path integral (QUAPI) technique[8] 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 strengths[10] although it becomes extremely difficult to predict quantum dissipative dynamics at very low temperatures.[11] Recent progress in ultrafast time resolved spectroscopy has stimulated methodological developments, and a large number of efficient, approximate reduced density matrix approaches are available.[13]

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.[2] Alternatively, the wave-function approach such as the multiconfigurational time-dependent Hartree (MCTDH) approach[21] and its multilayer variant (ML-MCTDH)[23] 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,[24] and for including finite temperature effects with the help of statistical sampling of the bath initial conditions.[25]

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, [27] 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.[30] 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.[31] 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.[34] Increasing the multiplicity of these Ansätze leads to numerically exact dynamics in the Holstein model and the sub-Ohmic spin boson model.[36] 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.[25] 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.

## 2Model

### 2.1Hamiltonian

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.[42] 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.[1] 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 )

where

Here, () are Pauli operators defined as and with and representing two electronic localized states. and is the energy bias and the coupling constant between two electronic states, respectively. () is creation (annihilation) operator for the bosonic bath mode of frequency , and is the strength of the coupling between the system and the th mode.

Owing to Wick’s theorem,[44] all the fluctuations and the dissipation process of can be specified by a two-body correlation function:[46]

where denotes averaging over with being the inverse temperature. The bath spectral density is defined in terms of the coupling strength as[2]

Several forms of are employed in the literature based on model assumptions. In this study, we adopt the following spectral density form

where represents the strength of the coupling, and provides a phenomenological frequency cutoff. Those spectral densities with , , and are referred to as sub-Ohmic, Ohmic, and super-Ohmic, respectively.[42]

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,[47] 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.[23] Following Ref. , we introduce a density of frequencies defined on , where is the upper bound of the frequency, and discretize the continuum of frequencies as

where is the number of discrete bath modes and . The parameter for each is then given by . The precise functional form of does not affect the final outcome if sufficient bath modes are included. For efficient numerical calculations, however, the form of can be chosen as

where

Here, is the factor guaranteeing . For the Ohmic bath (), is expressed as . For the sub-Ohmic bath (), can only be given implicitly through Eq. (Equation 2).

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. (Equation 1) since the goal of our study is to demonstrate the validity of our finite temperature approach.

### 2.2Time dependent variational approach with multiple Davydov trial states

In general, the time evolution of the wave function for the total Hamiltonian, , is described with the schrödinger equation,

In this work, we employ two trial wave functions, namely, the multiple and Ansatz, to solve the above schrödinger equation. Both are known as the multiple Davydov Ansätze. An extension to the single Ansatz, the multiple 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.[35] The time-dependent version of the multi- Ansatz can be written as

where is the vacuum state of the boson bath. The variational parameters and are the amplitudes in states and , respectively. and represent bath-mode displacements with and denoting the th coherent state and the th bath mode, respectively. is the multiplicity, which denotes the number of single states in Eq. (Equation 5). For , the multi- Ansatz reduces to the standard Davydov Ansatz. Similarly, the multi- Ansatz is given by[39]

where are the displacements of th bath mode on any electronic states. Thus, the multi- Ansatz is a simplified version of the multi- Ansatz, since the bath displacements of the multiple () trial state is site-dependent (site-independent). For , the multi- Ansatz reduces to the single Ansatz.

The time-dependent variational parameters, , , and are determined by adopting the Lagrangian formalism of the Dirac-Frenkel time-dependent variational principle. The Lagrangian associated with the trial state is given by

where the operator denotes . The Dirac-Frenkel time-dependent variational principle yields the equations of motion for the variational parameters,

where are the variational parameters, i.e., , , and , and is the complex conjugate of . Similarly, time evolution for the multi- 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. .

The expectation value of an observable can be expressed as , where is the total density operator, and is a time independent operator. Using the multi- or Ansatz, the expectation value of the observable of interest at zero temperature can be calculated as

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

### 2.3Observables 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.[24] The initial density matrix for the entire system is assumed to have a factorized form, i.e. , where . The extension to more general initial conditions with superposition of and is straightforward, which is important for modeling nonlinear spectroscopy.

The expectation value of an observable at finite temperatures can be expressed as

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

where denotes a direct product of coherent states , , , , for the discrete bath modes, and is expressed as . Each runs over all of the feasible coherent states. The element of area on the complex plane of denotes , which and are the real and imaginary part of , respectively. The equilibrium density matrix of the bath at a finite temperature is a diagonal matrix, and can be expressed as[48]

where represents the diagonal elements of the density matrix in the coherent state representation and can be expressed as[50]

As shown in Eq. (Equation 11), is a positive defined function of and can be seen as a probability density. Substituting Eq. (Equation 10) into Eq. (Equation 9), the observables at finite temperature can be obtained by the average according to the probability density as

For the second equality, we have use the multi- or Ansatz, , and denotes a trial state with initial bath displacements of at . For the case of the multi- Ansatz, initial condition parameters are , , for and for all and . Likewise, initial parameters of the multi- Ansatz are , , for and for all and .

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

where is the sampling number. The configuration for the bath is numerically generated according to by importance sampling, where is the Boltzmann distribution used as the weighting function in the importance sampling procedure. Letting and , in Eq. (Equation 11) can be partitioned into two independent Gaussian distribution as

where can be taken as the variance of the Gaussian distribution. To avoid singularity, the initial displacements in the trial states is determined by setting , where noise satisfying the uniform distribution is added to the variational parameters of the initial states. From the definition of , a lower temperature or the higher frequency gives a smaller . 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, for all . In this case, the observable expression of Eq. (Equation 12) or Eq. (Equation 13) reduces to Eq. (Equation 8).

## 3Results and discussion

In this section, we compare results from our variational approach with those of the QUAPI method[8] 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, , calculated as

All units in the numerical calculations are scaled by the cutoff frequency . In this section, for simplicity, we choose zero bias as a large energy bias requires a large multiplicity.[51] The electronic coupling constant and the coupling strength of spectral density are set as and , respectively. The number of the discrete bath modes is fixed at . 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 . 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.[8] The short-time density matrix propagator describes time evolution over a time slice . This splitting is by construction exact in the limit . With finite time slice in practical calculations, the numerical error can be eliminated by choosing sufficiently small 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 , thus a memory time window, , should be defined to handle the bath correlation truncation beyond a certain time span. Neglected beyond , all correlations are included exactly within a finite memory time of . To reach convergence, the memory time window should be enlarged by increasing 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 in Figs. Figure 1 and Figure 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 and , respectively. The convergence failure of the QUAPI in deep sub-Ohmic regime is discussed in Subsection IIIB.

### 3.1The Ohmic Regime

Figures Figure 1(a) and Figure 1(b) display time evolution of calculated by the multi- Ansätze, the multi- Ansätze and the QUAPI approach in the Ohmic case () at two temperatures ( and ), respectively. There is perfect agreement between the Ansatz and the QUAPI approach at the low temperature , as shown in the right panel of Figure 1(a). On the other hand, the agreement between the Ansatz and the QUAPI approach is unsatisfactory. As discussed in Ref. , 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- Ansatz provides accurate population dynamics by , and the agreement between the multi- Ansatz and the QUAPI approach is good even at . This clearly demonstrates that the multi- Ansätze is more flexible than the multi- Ansätze.

In the left panel of Figure 1(b), dynamics obtained by both the Ansatz and the QUAPI approach agree with each other at high temperature () despite small difference between the Ansatz and QUAPI populations. The result by the Ansatz also appears to be similar to that by the case. The small difference between the multi- Ansätze and the QUAPI approach may be due to insufficient number of bath modes. From Eq. (Equation 14) and the definition of the variance of the Gaussian distribution , the value of becomes large at high temperature. The increase of leads to large values of and even for high frequency modes which can be ignored at low temperatures. Therefore, the number of bath modes required as well as the sampling number become large due to large value of at high temperatures. The accuracy of the Ansatz with at high temperature is improved significantly unlike in the low temperature regime, and are similar to that of the 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.

### 3.2The Sub-Ohmic Regime

In this subsection, we focus on the sub-Ohmic regime. Figures Figure 2(a) and Figure 2(b) present results of at and , respectively. These results demonstrate that the variational approach with the multi- Ansätze provides excellent dynamics simulation over a range of temperatures and bath spectral exponents . Calculated dynamics in Figure 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, , which represents the dissipated environment energy after electronic excitation in accordance with vertical Franck-Condon transition in the electron transfer theory.[14] From the definition of spectral density of Eq. (Equation 1), the value of in the case of is large compared with the Ohmic bath case when other parameters except the value of are fixed. A large corresponds to larger fluctuations according to the fluctuation-dissipation relation,[14] 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- and multi- Ansätze in the case of is not dissimilar to the Ohmic case.

Finally, we discuss potential applicability of the multiple Davydov Ansätze in the deep sub-Ohmic regime (), 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. . To obtain correct long time dynamics, one would need to increase the memory window size , making QUAPI convergence more difficult to achieve. It has been demonstrated that the multi- 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 at zero temperature, and is a reliable, efficient method for describing quantitatively accurate dynamics of the deep sub-Ohmic regime.[36] From Figs. Figure 1 and Figure 2, the validity of the multiple Davydov Ansätze at finite temperatures also seems to be independent of the exponent , 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 smaller than that in Figure 2. As a helpful reference, we also compare our results with those from the thermal field dynamics (TFD),[52] similar to the approach of Borrelli and Gelin.[53] Figure 3 presents result of at . The dynamics by Ansatz is in perfect agreement with those of QUAPI and TFD, as shown in left panel of Figure 3. As demonstrated in the right panel of Figure 3, it is found that the multi- Ansatz slightly overestimates the amplitude of population oscillation at , and one needs to increase the multiplicity further to for numerical convergence. Fortunately, no large increase in the multiplicity of the multi- Ansatz is needed for adequate convergence with decreasing . To explore how the required multiplicity is affected by thermal fluctuations, we have considered a case of relatively small system- bath coupling strength in all numerical calculations of this study. From some previous studies of the multiple Davydov Ansätze,[36] 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- and the multi- 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 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.

## 4Concluding 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.[54] 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.[54] 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. . 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.[6] 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 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 approach[53], 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.[52]

### References

1. A. Nitzan, Chemical dynamics in condensed phases: relaxation, transfer and reactions in condensed molecular systems (Oxford University Press, New York, 2006).
2. V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems (Wiley-VCH, Weinheim, 2004).
3. Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989).
4. A. Ishizaki and Y. Tanimura, J. Phys. Soc. Jpn. 74, 3131 (2005).
5. A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234111 (2009).
6. Y. Fujihashi, G. R. Fleming, and A. Ishizaki, J. Chem. Phys. 142, 212403 (2015).
7. L. Chen, Y. Zhao, and Y. Tanimura, J. Phys. Chem. Lett. 6, 3110 (2015).
8. N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4600 (1995).
9. N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4611 (1995).
10. P. Nalbach, A. Ishizaki, G. R. Fleming, and M. Thorwart, New J. Phys. 13, 063040 (2011).
11. Z. Tang, X. Ouyang, Z. Gong, H. Wang, and J. Wu, J. Chem. Phys. 143, 224112 (2015).
12. H. T. Chen, G. Cohen, and D. R. Reichman, e-print arXiv:1610.09402.
13. Y.-C. Cheng and G. R. Fleming, Annu. Rev. Phys. Chem. 60, 241 (2009).
14. A. Ishizaki, T. R. Calhoun, G. S. Schlau-Cohen, and G. R. Fleming, Phys. Chem. Chem. Phys. 12, 7319 (2010).
15. S. F. Huelga and M. B. Plenio, Contemp. Phys. 54, 181 (2013).
16. A. Chenu and G. D. Scholes, Annu. Rev. Phys. Chem. 66, 69 (2015).
17. L. Chen, P. Shenai, F. Zheng, A. Somoza, and Y. Zhao, Molecules 20, 15224 (2015).
18. M. K. Lee, P. Huo, and D. F. Coker, Annu. Rev. Phys. Chem. 67, 639 (2016).
19. P. Kukura, D. W. McCamant, and R. A. Mathies, Annu. Rev. Phys. Chem. 58, 461 (2007).
20. T. A. Oliver, N. H. Lewis, and G. R. Fleming, Proc. Natl. Acad. Sci. USA 111, 10061 (2014).
21. H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990).
22. F. Matzkies and U. Manthe, J. Chem. Phys. 110, 88 (1999).
23. H. Wang, M. Thoss, and W. H. Miller, J. Chem. Phys. 115, 2979 (2001).
24. H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003).
25. H. Wang and M. Thoss, J. Chem. Phys. 124, 034114 (2006).
26. U. Manthe, J. Chem. Phys. 128, 164116 (2008).
27. M. J. Škrinjar, D. V. Kapor, and S. D. Stojanović, Phys. Rev. A 38, 6402 (1988).
28. W. Förner, Phys. Rev. B 50, 6291 (1996).
29. J. Sun, B. Luo, and Y. Zhao, Phys. Rev. B 82, 014305 (2010).
30. A. S. Davydov, Phys. Scr. 20, 387 (1979).
31. T. D. Huynh, K. W. Sun, M. F. Gelin, and Y. Zhao, J. Chem. Phys. 139, 104103 (2013).
32. V. Chorošajev, A. Gelzinis, L. Valkunas, and D. Abramavicius, J. Chem. Phys. 140, 244108 (2014).
33. L. Chen, M. F. Gelin, W. Domcke, Y. Zhao, J. Chem. Phys. 142 164106 (2015).
34. N. Zhou, L. Chen, Y. Zhao, D. Mozyrsky, and V. Chernyak, Phys. Rev. B 90, 155135 (2014).
35. N. Zhou, L. Chen, D. Xu, V. Chernyak, and Y. Zhao, Phys. Rev. B 91, 195129 (2015).
36. L. Wang, L. Chen, N. Zhou, and Y. Zhao, J. Chem. Phys. 144, 024101 (2016).
37. N. Zhou, L. Chen, Z. Huang, K. Sun, Y. Tanimura, and Y. Zhao, J. Phys. Chem. A 120, 1562 (2016).
38. T. Deng, Y. Yan, L. Chen, and Y. Zhao, 144, 144102 (2016).
39. N. Zhou, Z. Huang, J. Zhu, V. Chernyak, and Y. Zhao, J. Chem. Phys. 143, 014113 (2015).
40. Z. Huang, L. Wang, C. Wu, L. Chen, F. Grossmann, and Y. Zhao, Phys. Chem. Chem. Phys. 19, 1655 (2017).
41. Y. Fujihashi, L. Chen, A. Ishizaki, J. Wang, and Y. Zhao, J. Chem. Phys. 146, 044101(2017).
42. A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
43. U. Weiss, Quantum dissipative systems (World scientific, Singapore, 1999).
44. G. C. Wick, Phys. Rev. 80, 268 (1950).
45. A. L. Fetter and J. D. Walecka, Quantum theory of many-particle systems (McGraw-Hill, New York, 1971).
46. R. Kubo, M. Toda, and N. Hashitsume, Statical physics II: nonequilibrium statical mechanics (Springer, Berlin, Heidelberg, 1985).
47. R. Bulla, N. H. Tong, and M. Vojta, Phys. Rev. Lett. 91, 170601 (2003).
48. R. J. Glauber, Phys. Rev. Lett. 131, 2766 (1963).
49. E. C. G. Sudarshan, Phys. Rev. Lett. 10, 277 (1963).
50. M. O. S. M. Hillery, R. F. O’connell, M. O. Scully, and E. P. Wigner, Phys. Rep. 106, 121 (1984).
51. R. Borrelli and A. Peluso, J. Chem. Phys. 144, 114102 (2016).
52. L. Chen and Y. Zhao, “Finite-temperature dynamics of a Holstein polaron: thermofield dynamics approach” (unpublished).
53. R. Borrelli and M. F. Gelin, J. Chem. Phys. 145, 224101 (2016).
54. J. Schulze and O. Kühn, J. Phys. Chem. B 119, 6211 (2015).
55. J. Schulze, M. F. Shibl, M. J. Al-Marri, and O. Kühn, J. Chem. Phys. 144, 185101 (2016).
13718