Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems

TL;DR

This paper develops quantum Fokker-Planck and Smoluchowski equations with low-temperature corrections for simulating coupled electronic and vibrational dynamics in open quantum systems, extending to multi-state cases.

Contribution

It introduces rigorous quantum equations with non-Markovian low-temperature corrections and extends them to multi-state systems, enabling efficient non-adiabatic dynamics simulations.

Findings

Numerical validation for single-state Brownian systems.

Application to multi-state double-well systems comparing with other methods.

Quantum low-temperature corrections are crucial for accurate transient spectra.

Abstract

Simulating electron-nucleus coupled dynamics poses a non-trivial challenge and an important problem in the investigation of ultrafast processes involving coupled electronic and vibrational dynamics. Because irreversibility of the system dynamics results from thermal activation and dissipation caused by the environment, in dynamical studies, it is necessary to include heat bath degrees of freedom in the total system. When the system dynamics involves high-energy electronic transitions, the environment is regarded to be in a low-temperature regime and we must treat it quantum mechanically. In this paper, we present rigorous and versatile approaches for investigating the dynamics of open systems with coupled electronic and vibrational degrees of freedom within a fully quantum mechanical framework. These approaches are based on a quantum Fokker-Planck equation and a quantum Smoluchowski equation employing a heat bath with an Ohmic spectral density, with non-Markovian low-temperature correction terms, and extensions of these equations to the case of multi-state systems. The accuracy of these equations was numerically examined for a single-state Brownian system, while their applicability was examined for multi-state double-well systems by comparing their results with those of the fewest-switch surface hopping and Ehrenfest methods with a classical Markovian Langevin force. Comparison of the transient absorption spectra obtained using these methods clearly reveals the importance of the quantum low-temperature correction terms. These equations allow us to treat non-adiabatic dynamics in an efficient way, while maintaining numerical accuracy. The C++ source codes that we developed, which allow for the treatment of the phase and coordinate space dynamics with any single-state or multi-state potential forms, are provided as Supporting Information.

Figures

Tables

Table 1. Table 1: PSD[N-1/N] coefficients used for K = 1 𝐾 1 K=1 – 5 5 5 .

$K$	$k$	$\beta \hslash {\nu }_k$	${\eta }_k$
$1$	$1$	$7.745967$	$2.5$
2	$1$	$6.305939$	$1.032824$
	$2$	$19.499618$	$5.967176$
3	$1$	$6.2832903$	$1.000227$
	$2$	$12.9582867$	$1.300914$
	$3$	$36.1192894$	$11.198859$
4	$1$	$6.283185$	$1.000000$
	$2$	$12.579950$	$1.015314$
	$3$	$20.562598$	$1.905605$
	$4$	$57.787940$	$18.079081$
5	$1$	$6.283185$	$1.000000$
	$2$	$12.566542$	$1.000262$
	$3$	$19.004690$	$1.113033$
	$4$	$29.579276$	$2.800147$
	$5$	$84.536926$	$26.586558$

Table 2. Table 2: The parameter values for the numerical tests.

Symbol	Value		Symbol	Value		Symbol	Value
${\omega }_0$	$400$	${\mathrm{cm}}^{-1}$	${\omega }_e$	${\omega }_0$		${\omega }_f$	$1.5{\omega }_0$
$\Delta E$	$2,000$	${\mathrm{cm}}^{-1}$	$E_{\mathrm{gap}}^{e-g}$	$1,000$	${\mathrm{cm}}^{-1}$	$E_{\mathrm{gap}}^{f-e}$	$5,000$	${\mathrm{cm}}^{-1}$
$L_0$	$10$		$q^{†}$	$1$		${\sigma }^{†}$	$1$

Equations

$${\hat{H}}_{\mathrm{tot}}(p,q;\vec{p},\vec{x})$$

$$\hat{H}(p,q)\equiv \frac{\hslash {\omega }_0}{2}{\hat{p}}^2+\sum _{j,k}|j⟩U_{jk}^{\mathrm{d}}(\hat{q})⟨k|.$$

$$\hslash \bar{\omega }\simeq \sqrt{\hslash {\omega }_0{\frac{{\partial }^2}{\partial q^2}{U}_{j_0{j}_0}(q)|}_{q=q_0}},$$

$${\hat{H}}_{\mathrm{B}}(\vec{p},\vec{x};q)$$

$$R(t)=\frac{2}{\pi }{\int }_0^{\infty }d\omega \frac{\mathcal{J}(\omega )}{\omega }\cos \omega t$$

$$C(t)=\frac{2}{\pi }{\int }_0^{\infty }d\omega \mathcal{J}(\omega )\left(n(\omega )+\frac{1}{2}\right)\cos \omega t,$$

$$n(\omega )+\frac{1}{2}\simeq \frac{1}{\beta \hslash }\frac{1}{\omega }+\sum _k^K\frac{2{\eta }_k}{\beta \hslash }\frac{\omega }{{\omega }^2+{\nu }_k^2}$$

$$\mathcal{J}(\omega )$$

$$R(t)=\frac{\zeta }{{\omega }_0}\cdot 2\delta (t)$$

$$\begin{array}{rl}C(t)\simeq C_K(t)& =\frac{\zeta }{{\omega }_0}\left(\frac{1}{\beta \hslash }+\sum _k^K\frac{2{\eta }_k}{\beta \hslash }\right)\cdot 2\delta (t)\\ & -\sum _k^K\frac{\zeta }{{\omega }_0}\frac{2{\eta }_k}{\beta \hslash }\cdot {\nu }_k{e}^{-{\nu }_k|t|}.\end{array}$$

$${\rho }_{jk}^{\mathrm{d}}(z,z^{\prime },t)$$

$${\mathbf{W}}^{\mathrm{d}}(p,q,t)$$

$$\begin{array}{rl}\frac{\partial }{\partial t}{\mathbf{W}}_{\vec{n}}^{\mathrm{d}}(p,q,t)& =-\left({\mathcal{L}}_{\mathrm{qm}}^{\mathrm{d}}(p,q)+\sum _k^K{n}_k{\nu }_k+{\hat{\Xi }}_K^{\mathrm{d}}(p,q)\right){\mathbf{W}}_{\vec{n}}^{\mathrm{d}}(p,q,t)\\ & -\sum _k^K{\hat{\Phi }}^{\mathrm{d}}(p,q){\mathbf{W}}_{\vec{n}+{\vec{e}}_k}^{\mathrm{d}}(p,q,t)-\sum _k^K{n}_k{\nu }_k{\hat{\Theta }}_k^{\mathrm{d}}(p,q){\mathbf{W}}_{\vec{n}-{\vec{e}}_k}^{\mathrm{d}}(p,q,t),\end{array}$$

$${\mathcal{L}}_{\mathrm{qm}}^{\mathrm{d}}(p,q)$$

$$\begin{array}{rl}\mathcal{K}(p,q)\mathbf{W}(p,q)& \equiv {\omega }_{0}p\frac{\partial }{\partial q}\mathbf{W}(p,q)\end{array}$$

$$\displaystyle\begin{aligned} \mathcal{U}_{\mathrm{qm}}^{\mathrm{d}}(p,q)\bm{W}(p,q)&\equiv\frac{i}{\hbar}\Bigl{(}\bm{U}^{\mathrm{d}}(q){\star}\bm{W}(p,q)-\bm{W}(p,q){\star}\bm{U}^{\mathrm{d}}(q)\Bigr{)}\end{aligned}$$

$$\star$$

$$\text{underaccent}\to {\partial }_{x}f(x)=f(x)\text{underaccent}\leftarrow {\partial }_x$$

$${\hat{\Phi }}^{\mathrm{d}}(p,q)$$

$${\hat{\Theta }}_k^{\mathrm{d}}(p,q)$$

$$\displaystyle\begin{split}\hat{\Xi}_{K}^{\mathrm{d}}(p,q)&\equiv-\frac{\zeta}{\omega_{0}}\frac{\partial}{\partial p}\biggl{(}\omega_{0}p+\frac{1}{\beta\hbar}\frac{\partial}{\partial p}\biggr{)}+\sum_{k}^{K}\hat{\Phi}^{\mathrm{d}}(p,q)\hat{\Theta}_{k}^{\mathrm{d}}(p,q).\end{split}$$

$${\mathbf{f}}^{\mathrm{d}}(q,t)$$

$$\begin{array}{rl}\frac{\partial }{\partial t}{\mathbf{f}}_{\vec{n}}^{\mathrm{d}}(q,t)& =-\left[{\mathcal{E}}^{\mathrm{d}}(q)+\sum _k^K{n}_k{\nu }_k+\frac{{\omega }_0}{\zeta }\left({\mathcal{F}}^{\mathrm{d}}(q)+{\hat{\Xi }}_K^{\mathrm{od},\mathrm{d}}(q)\right)\right]{\mathbf{f}}_{\vec{n}}^{\mathrm{d}}(q,t)\\ & -\sum _k^K{\hat{\Phi }}^{\mathrm{od},\mathrm{d}}(q){\mathbf{f}}_{\vec{n}+{\vec{e}}_k}^{\mathrm{d}}(q,t)-\frac{{\omega }_0}{\zeta }\sum _k^K{n}_k{\nu }_k{\hat{\Theta }}_k^{\mathrm{od},\mathrm{d}}(q){\mathbf{f}}_{\vec{n}-{\vec{e}}_k}^{\mathrm{d}}(q,t),\end{array}$$

$${\mathcal{E}}^{\mathrm{d}}(q)\mathbf{f}(q,t)$$

$${\mathcal{F}}^{\mathrm{d}}(q)\mathbf{f}(q,t)$$

$${\hat{\Phi }}^{\mathrm{od},\mathrm{d}}(q)$$

$${\hat{\Theta }}_k^{\mathrm{od},\mathrm{d}}(q)$$

$${\hat{\Xi }}_K^{\mathrm{od},\mathrm{d}}(q)$$

$$U(q)$$

$$f(q)$$

$$⟨q^2{⟩}_{\beta ,\zeta }$$

$$=\frac{{\omega }_0}{\beta \hslash }\sum _{k=-\infty }^{\infty }\frac{1}{{\omega }_0^2+\left|{\tilde{\nu }}_k\right|\zeta +{{\tilde{\nu }}_k}^2}$$

$${\tilde{C}}_q(\omega )$$

$${\Delta }_{\vec{n}}$$

$$C_q(t)$$

Full text

I INTRODUCTION

Understanding non-adiabatic dynamics in electronic and bio-nanomaterials is fundamentally important in the study of many types of phenomena, ranging from photo-isomerization to spintronics. Recent advances in experimental technologies have made it possible to observe such non-adiabatic processes that take place on very short timescales Polli et al. (2010); Iwamura et al. (2011); Lewis and Fleming (2016); Dean et al. (2016); Miyata et al. (2017). Theoretical input regarding the complex profiles of potential energy surfaces (PESs) and the non-adiabatic coupling among PESs is important for analyzing such ultrafast transport processes, in particular, those in materials involving biomolecular aggregates and crystalline solar cells Prokhorenko et al. (2016); Kato and Ishizaki (2018). For such systems, the surrounding molecules act as a heat bath and also play an essential role in determining the nature of the transport processes, because they either promote or suppress the wavepacket motion of the system through thermal activation and relaxation.

In the study of systems of the type considered here, while nuclear motion is often treated using a semiclassical approach, which is applicable to the case of heavy nuclei, non-adiabatic transition processes must be described using a purely quantum mechanical approach, because transitions between discretized electronic states are, in their essence, quantum dynamics. For this reason, the effect of the environment should be treated using an open quantum model, even if the dynamics of the nuclei are semiclassical, because otherwise the quantum nature of the electronic transitions is not properly accounted for. Indeed, ignoring the quantum effects of the environment, in particular in the low temperature regime, in which quantum effects become very important, leads to unphysical behavior. For example, in the case that we employ a classical description of the environment in the low temperature regime, while electronic transitions and the motion of wavepackets are described by quantum mechanics, the positivity of the probability distributions of the electronic states cannot be maintained. This is a fundamental complication, known as the “positivity problem,” which imposes a well-known limitation on the applicability of the quantum master equation without the rotating wave approximation Dümcke and Spohn (1979); Pechukas (1994). The positivity problem arises because the classical treatment of the environment leads to the violation of the quantum fluctuation-dissipation (QFD) theorem Ford et al. (1988); Frantsuzov (1997); Dammak et al. (2009); Tanimura (2006).

The excited state dynamics of systems exhibiting ultrafast coupled electronic and vibrational dynamical processes have been investigated with models that explicitly take into account nuclear degrees of freedom and electronic states through approaches employing equations of motion for wave functions, density matrices, phase space distributions Thoss et al. (2000); Kühl and Domcke (2002); Tanimura and Mukamel (1994a); Tanimura and Maruyama (1997); Maruyama and Tanimura (1998); Ikeda and Tanimura (2017, 2018); Kapral and Ciccotti (1999); Xie et al. (2014), and Gaussian quantum wavepackets Ben-Nun and Martínez (1998); Ben-Nun et al. (2000); Makhov et al. (2014), and approaches utilizing mixed quantum-classical trajectories Tully (1990); Coker and Xiao (1995); Subotnik et al. (2016). However, many of these approaches were developed for isolated systems and were verified within the system which have a few degrees of freedom. Moreover, varieties of assumptions (in particular, assumptions regarding the quantum dynamical treatment of the couplings between the electronic states and the nuclear coordinates) were introduced in such approaches and the assumptions severely limit their range of applicability. Contrastingly, the multi-state quantum hierarchical Fokker-Planck equation (MS-QHFPE) approach, which is an extension of the quantum Fokker-Planck equation for the Wigner distribution function Caldeira and Leggett (1983); Tanimura and Wolynes (1991, 1992); Tanimura (2015) to multi-state systems Tanimura and Mukamel (1994a); Tanimura and Maruyama (1997); Maruyama and Tanimura (1998); Ikeda and Tanimura (2017, 2018) and is a variant of the hierarchical equations of motion (HEOM) theories Tanimura and Kubo (1989); Ishizaki and Tanimura (2005); Tanimura (2006), can treat any types of diabatic coupling and PES profiles with non-Markovian system-bath interactions described by a Drude spectral density. However, although the MS-QHFPE approach allows us to compute the dynamics described by a multi-state system-bath Hamiltonian numerically rigorously, integrating the equations of motion is very computationally intensive, in particular, for a system described by multi-dimensional PESs. Hence, presently, calculations carried out for two-dimensional systems are limited to the high-temperature Markovian case described by the MS-QFPE Ikeda and Tanimura (2018).

While it has been found that non-Markovian effects arising from non-Ohmic environments are important in the description of exciton/electron transfer phenomena Ishizaki and Fleming (2009); Fujihashi et al. (2015); Prokhorenko et al. (2016); Kato and Ishizaki (2018), the Ohmic heat bath model for nuclear dynamics has been (implicitly) employed in many investigations for models described by PESs that further coupled to a heat-bath, including a model that causes a Brownian/Drude spectral density after reducing the nuclear degrees of freedom Leggett (1984); Garg et al. (1985); Tanimura and Mukamel (1994b); Tanimura (2012). This results from the fact that the non-Markovian effects on the nuclear motion so far studied are regarded to be insignificant in such systems, in particular, when the damping on the nuclear motion is quite strong. For this reason, although there have been several investigations employing Drude environments for the nuclear dynamics carried out on systems including nonlinear vibrational responses Ishizaki and Tanimura (2006); Sakurai and Tanimura (2011); Ikeda et al. (2015), a ratchet system Kato and Tanimura (2013), and a resonant tunneling diode system Sakurai and Tanimura (2013, 2014), in this paper, we derive equations of motion for single-state and multi-state systems employing the Ohmic environment: low-temperature quantum Fokker-Planck equations (LT-QFPE) and low-temperature quantum Smoluchowski equations (LT-QSE) and their extensions to multi-state (MS) systems, MS-LT-QFPE and MS-LT-QSE. As seen in the theory of quantum Brownian motion, within a quantum mechanical description, an Ohmic bath exhibits peculiar behavior in momentum space Grabert et al. (1988); Pechukas et al. (2000); Ankerhold et al. (2001). We show that this difficulty can be avoided by properly treating the low temperature correction terms in the LT-QFPE and MS-LT-QFPE. In the case of a heat bath with an Ohmic spectral density, the LT-QFPE and LT-QSE are sufficiently accurate, while also being sufficiently simple in comparison to the QHFPE. These features make the LT-QFPE and LT-QSE suited for describing slowly decaying systems and systems rendered in multidimensional phase spaces. Also, it is noteworthy that many of the existing formalisms, including those of the quantum Fokker-Planck equation Caldeira and Leggett (1983); Tanimura and Wolynes (1991, 1992); Tanimura (2015) and Zusman equation Zusman (1980); Shi et al. (2009); Shi and Geva (2009), can be derived from the (MS-)LT-QFPE and (MS-)LT-QSE under certain conditions.

The organization of this paper is as follows. In Sec. II, we introduce a model with multiple electronic states described by the PESs that are coupled to a harmonic heat bath with an Ohmic spectral density. Then, we present the MS-LT-QFPE and MS-LT-QSE and their single PES forms, LT-QFPE and LT-QSE. In Sec. III, we present numerical results for single-state Brownian and multi-state double-well systems to illustrate the validity and applicability of these approaches. Section IV is devoted to concluding remarks. The C++ source codes that we developed are provided as Supporting Information.

II HAMILTONIAN AND FORMALISM

II.1 Model

Because the LT-QFPE and LT-QSE are the simpler, single-potential forms of the MS-LT-QFPE and MS-LT-QSE, we start with a multi-potential system. We consider a molecular system with multiple electronic states $\{|j\rangle\}$ coupled to the nuclear coordinates. For simplicity, we represent the nuclear coordinates by a single dimensionless coordinate, $q$. Here and hereafter, we employ a dimensionless coordinate and a dimensionless momentum defined in terms of the actual coordinate and momentum, $\bar{q}$ and $\bar{p}$, as $q\equiv\bar{q}\sqrt{m\omega_{0}/\hbar}$ and $p\equiv\bar{p}/\sqrt{m\hbar\omega_{0}}$, where $\omega_{0}$ is the characteristic vibrational frequency of the system and $m$ is the effective mass. The reaction coordinate is also bilinearly coupled to the harmonic bath coordinates, $\vec{x}\equiv(\dots,x_{\xi},\dots)$. The Hamiltonian of the total system is expressed as Caldeira and Leggett (1983)

[TABLE]

where the system Hamiltonian, $\hat{H}(p,q)$, is defined as

[TABLE]

Here, the nuclear and electronic operators are denoted by hats, and the direct products with the unit operator in the kinetic and bath terms ($\otimes\hat{1}$) are omitted. The diagonal element $U_{jj}^{\mathrm{d}}(q)$ is the diabatic PES of $|j\rangle$, and the off-diagonal element $U_{jk}^{\mathrm{d}}(q)$ $(j{\neq}k)$ represents the diabatic coupling between $|j\rangle$ and $|k\rangle$. The vibrational frequency $\bar{\omega}$ at a local minimum of the potential $q_{0}$ is determined by the curvature of the PESs as

[TABLE]

where $j_{0}$ is a primary state of the vibrational dynamics. Therefore, the frequency $\mathstrut\omega_{0}$ is chosen to be $\hbar\omega_{0}\simeq\partial^{2}U_{j_{0}j_{0}}(q)/\partial q^{2}|_{q=q_{0}}$ in order to have $\bar{\omega}\simeq\omega_{0}$. The bath Hamiltonian $\hat{H}_{\mathrm{B}}(\vec{p},\vec{x};q)$ is defined as

[TABLE]

where $\omega_{\xi}$, $p_{\xi}$, and $g_{\xi}$ are the vibrational frequency, conjugate momentum, and system-bath coupling constant of the $\xi$th dimensionless bath mode, $x_{\xi}$.

[TABLE]

We choose the coefficients $\nu_{k}$ and $\eta_{k}$ so as to realize the relation

[TABLE]

for finite $K$, where the first term on the right-hand side is the classical contribution from the temperature, and the remaining terms are the quantum low-temperature (QLT) corrections. The Matsubara decomposition scheme (MSD) can be applied straightforwardly to the above. In this scheme, we set $\eta_{k}=1$ and $\nu_{k}=\tilde{\nu}_{k}$, where $\tilde{\nu}_{k}\equiv 2\pi k/\beta\hbar$ is the $k$th Matsubara frequency Tanimura (1990); Ishizaki and Tanimura (2005). In this paper, we employ the Padé spectral decomposition $[N{-}1/N]$ (PSD$[N{-}1/N]$) scheme to enhance the computational efficiency while maintaining the accuracy Hu et al. (2010, 2011); Ding et al. (2012).

In order to reduce the computational times for the computations of the non-adiabatic dynamics with any forms of the PESs, here we employ an Ohmic spectral density, expressed as

[TABLE]

where $\zeta$ is the system-bath coupling strength. Then we have

[TABLE]

In the case of a harmonic PES with frequency $\omega_{0}$, the conditions $\zeta<2\omega_{0}$, $\zeta=2\omega_{0}$, and $\zeta>2\omega_{0}$ correspond to the underdamped, critically damped, and overdamped cases, respectively.

In Fig. 1, we plot $C_{K}(t)$ for various values of the cutoff $K$, at the temperature $\beta\hbar\omega_{0}=7.47$. As this figure and Eq. (8b) indicate, the fluctuation term is always non-Markovian due to the quantum nature of the noise, and it can be regarded as Markovian only in the high-temperature limit, $\beta\hbar\omega_{\mathrm{0}}\ll 1$, in which the heat bath exhibits classical behavior Tanimura (2006); Weiss (2011). This is an important conclusion obtained from the QFD theorem, namely, that the negative non-Markovian terms appear in the case that we do not use a time-coarse-grained, Markovian description. Note that, in the case that we employ an Ohmic spectral density without cutoff functions (e.g. Lorentzian cutoff and exponential cutoff), some of the physical observables, including the mean square of the momentum, $\langle p^{2}\rangle$, diverge due to the divergence of the first and second terms in Eq. (8b) under the infinite summation of the Matsubara frequencies, if there is no finite cutoff function. This divergence is often referred to as the ultra-violet divergence Grabert et al. (1988); Pechukas et al. (2000); Ankerhold et al. (2001). In practice, we can ignore QLT correction terms whose frequencies are sufficiently greater than the characteristic frequency of the system, because the random force generated by such terms is averaged out over a sufficiently short timescale that its influence on the dynamics of interest is negligible. In this way, we are able to calculate non-diverging physical observables, e.g.  the mean square of the coordinate, $\langle q^{2}\rangle$, by simply ignoring the contribution from the high-frequency QLT correction terms by implementing the cutoff $K$, while diverging physical observables still tend to diverge.

II.2 Multi-State Low-Temperature Quantum Fokker-Planck Equations

The state of the total system is represented by the density operator, $\hat{\rho}_{\mathrm{tot}}(z,z^{\prime},\vec{x},\vec{x}\,^{\prime})\equiv\langle\vec{x}|\langle z|\hat{\rho}_{\mathrm{tot}}(t)|z^{\prime}\rangle|\vec{x}\,^{\prime}\rangle$, where $|z\rangle$ and $|\vec{x}\rangle$ are the eigenstate of the system and bath coordinate operators, respectively. We consider the reduced density matrix in the diabatic representation of the system subspace, defined as

[TABLE]

where $\hat{\rho}(z,z^{\prime},t)\equiv\mathrm{Tr}_{\mathrm{B}}\{\hat{\rho}_{\mathrm{tot}}(z,z^{\prime},\vec{x},\vec{x}^{\prime},t)\}$ is the reduced density operator and $\mathrm{Tr}_{\mathrm{B}}\{\dots\}\equiv\int\!d\vec{x}\,\int\!d\vec{x}\,^{\prime}\delta(\vec{x}-\vec{x}\,^{\prime})\{\dots\}$ represents the trace operation in the bath subspace. The diagonal and off-diagonal elements, $\rho^{\mathrm{d}}_{jj}$ and $\rho^{\mathrm{d}}_{jk}$ ($j{\neq}k$), represent the population of $|j\rangle$ and the coherence between $|j\rangle$ and $|k\rangle$, respectively. Hereafter, we employ the matrix forms of the reduced density matrix and the diabatic PESs: $\{\bm{\rho}^{\mathrm{d}}(z,z^{\prime})\}_{jk}\equiv\rho^{\mathrm{d}}_{jk}(z,z^{\prime})$ and $\{\bm{U}^{\mathrm{d}}(z)\}_{jk}\equiv U^{\mathrm{d}}_{jk}(z)$.

We now introduce the Wigner distribution function, which is the quantum analogy of the classical distribution function in phase space. For a multi-state system, the multi-state Wigner distribution function (MS-WDF) is defined as Tanimura and Mukamel (1994a); Tanimura and Maruyama (1997); Maruyama and Tanimura (1998); Ikeda and Tanimura (2017, 2018)

[TABLE]

where $q\equiv(z+z^{\prime})/2$ and $r\equiv z-z^{\prime}$. Both $p$ and $q$ are now c-numbers in this phase space representation.

The reduced dynamics of $\bm{\rho}^{\mathrm{d}}(z,z^{\prime},t)$ and $\bm{W}^{\mathrm{d}}(p,q,t)$ are expressed in the path integral framework using the Feynman-Vernon influence functional Feynman and Vernon (1963). Their time evolutions can be described by a set of time differential equations in the HEOM form (see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems in the Supporting Information). In the present case, these equations are the following:

[TABLE]

where $\vec{n}\equiv(\dots,n_{k},\dots)$ is a $K$-dimensional multi-index whose components are all non-negative integers and $\vec{e}_{k}\equiv(0,\dots,1,0,\dots)$ is the $k$th unit vector. The multi-index $\vec{n}$ represents the index of the hierarchy, and physically, the first hierarchical element, $\bm{W}_{\vec{0}}^{\mathrm{d}}(p,q,t)$, corresponds to the MS-WDF, $\bm{W}^{\mathrm{d}}(p,q,t)$. The rest of the hierarchical elements serve only to facilitate treatment of the non-Markovian system-bath interaction that arises from the QLT effects.

[TABLE]

Here, we have introduced the star operator, $\star$, which represents the Moyal product, defined as Moyal (1949); İmre et al. (1967)

[TABLE]

The differentiation operators from the left and right appearing here are defined as

[TABLE]

[TABLE]

The derivation of Eq. (11) is presented in Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information. Because Eq.  (11) is a generalization of the multi-state quantum Fokker-Planck equation (MS-QFPE) Tanimura and Mukamel (1994a); Tanimura and Maruyama (1997); Maruyama and Tanimura (1998); Ikeda and Tanimura (2017) valid in the low-temperature regime, we refer to these equations as the multi-state low-temperature quantum Fokker-Planck equations (MS-LT-QFPE).

For a single-state system, the matrices $\bm{W}^{\mathrm{d}}(p,q,t)$ and $\bm{U}^{\mathrm{d}}(q)$ reduce to scalar functions, $W(p,q,t)$ and $U(q)$. In this case, we refer to Eq. (11) as the low-temperature quantum Fokker-Planck equations (LT-QFPE). These equations can be understood as an extension of the quantum Fokker-Planck equation (QFPE) Caldeira and Leggett (1983); Tanimura and Wolynes (1991, 1992); Tanimura (2015).

The conventional (multi-state) quantum hierarchical Fokker-Planck equations ((MS-)QHFPE) with a Drude spectral density, $\mathcal{J}^{\mathrm{D}}(\omega)\propto\omega\gamma_{\mathrm{D}}^{2}/(\omega^{2}+\gamma_{\mathrm{D}}^{2})$, where $\gamma_{\mathrm{D}}$ is the cutoff frequency, are capable of treating systems subject to non-Markovian noise, and are not limited to the case of an Ohmic spectral density Sakurai and Tanimura (2011); Kato and Tanimura (2013); Sakurai and Tanimura (2014); Tanimura (2015). However, the (MS-)QHFPE require a ($K+1$)-dimensional multi-index $\vec{n}^{\prime}\equiv(n_{0},\vec{n})$ (i.e. the additional index $n_{0}$) to describe non-Markovian dynamics caused by a finite value of $\gamma_{\mathrm{D}}$, and therefore computationally more expensive than the (MS-)LT-QFPE. Moreover, the (MS-)QHFPE become unstable in the Ohmic limit (i.e. $\gamma_{\mathrm{D}}\gg 1$) at low temperatures due to the fast decaying terms with $\gamma_{\mathrm{D}}$, while the (MS-)LT-QFPE is sufficiently accurate and also being sufficiently simple in comparison to the (MS-)QHFPE. Thus, although applicability of these equations is limited to the Ohmic case, the computational cost to solve the (MS-)LT-QFPE is suppressed than that to solve the (MS-)QHFPE. These features make the (MS-)LT-QFPE and (MS-)LT-QSE suited for describing slowly decaying systems and systems rendered in multidimentional phase spaces. Note that in the case that the diabatic PESs of the system are harmonic, the MS-LT-QFPE yields the same results as the HEOM for a reduced electronic system with a Brownian spectral density Tanimura and Mukamel (1994b); Tanaka and Tanimura (2009). In Appendices A and B, we present a stochastic Liouville description of the (MS-)LT-QFPE and Langevin description of the LT-QFPE, respectively.

II.3 Multi-State Low-Temperature Quantum Smoluchowski Equations

In this section, we present the asymptotic form of Eq. (11) in the Smoluchowski limit, i.e. in the case $\zeta\gg\omega_{0}~{}\text{and}~{}\omega_{e}$, where $\omega_{e}$ is the characteristic frequency of the electronic transition dynamics. We introduce the following probability distribution in coordinate space:

[TABLE]

In the Smoluchowski limit, the equations of motion for $\bm{f}^{\mathrm{d}}(q)$ are

[TABLE]

[TABLE]

[TABLE]

The derivation of Eq. (17) is given in Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information. In the case of a single-state system, the matrices $\bm{f}^{\mathrm{d}}(q,t)$, $\bm{U}^{\mathrm{d}}(q)$ and $\bm{F}^{\mathrm{d}}(q)$ reduce to scalar functions, $f(q,t)$, $U(q)$ and $F(q)$, respectively. The relationship between Eq. (11) and Eq. (17) is similar to the relationship between the Fokker-Planck (Kramers) equation and the Smoluchowski equation Davies (1954); Risken (1989). For this reason, we refer to Eq. (17) as the (multi-state) low-temperature quantum Smoluchowski equations ((MS-)LT-QSE), while we refer to this as the (multi-state) Smoluchowski equation ((MS-)SE) in the high-temperature limit.

A quantum mechanical extension of the Smoluchowski equation valid in the low-temperature regime is known as the quantum Smoluchowski equation (QSE), which treats QLT effects in the framework of the Markovian approximation Ankerhold (2007); Pechukas et al. (2000); Ankerhold et al. (2001); Ankerhold and Grabert (2008); Maier and Ankerhold (2010a, b). However, because QLT corrections are in principle non-Markovian as shown in Eq. (8b), when we lower the bath temperature or we study a system with high energy, the QSE becomes inaccurate. Contrastingly, the (MS-)LT-QSE can describe non-Markovian terms that is necessary to satisfy the QFD theorem, the (MS-)LT-QSE is applicable to a wider range of physical conditions than the QSE, as shown in Appendix C. For electron transfer problems with harmonic PESs, an extension of the Smoluchowski equation to multi-state systems has been carried out as the Zusman equation (ZE) Zusman (1980); Roy and Bagchi (1994); Zusman and Beratan (1996); Shi et al. (2009). However, the original ZE theory does not treat quantum dynamical effects from electronic states properly, as shown in Appendix D. The MS-SE can be regarded as a generalization of the ZE for arbitrary PESs with describing the quantum dynamical effects from the electric states accurately. The MS-LT-QSE can be regarded as a generalization of the MS-SE with the QLT correction terms. Several extensions of the ZE theory valid in the low-temperature regime have been carried out as the following: The generalized ZE is constructed as an extension of the quantum Smoluchowski equation to multi-state systems Ankerhold and Lehle (2004). In the modified ZE theory Shi and Geva (2009) and stochastic ZE theory Lili et al. (2013), the effects of non-Markovian QLT terms are incorporated using an integral-differential equation similar to that used in second-order perturbation theories and using a stochastic differential equation, respectively. Contrastingly, the MS-LT-QSE is a non-Markovian, non-perturbative, and deterministic approach in the framework of the HEOM formalism. Note that when the diabatic PESs of the system are harmonic, the MS-LT-QSE gives the same result as the HEOM for a reduced electronic system with an overdamped Brownian/Drude spectral density Tanimura and Mukamel (1994b); Ishizaki and Fleming (2009). In Appendices A and B, we also present a stochastic Liouville description of the (MS-)LT-QSE and Langevin description of the LT-QSE, respectively.

III NUMERICAL RESULTS

In principle, with the (MS-)LT-QFPE, we are able to calculate various physical quantities with any desired accuracy by adjusting the number of low temperature corrections terms, while the (MS-)LT-QSE is sufficient for computing physical quantities under overdamped conditions. Here, we first examine the validity of Eqs. (11) and  (17) by presenting the results obtained from numerical integrations of these equations in the case of a Brownian oscillator, for which exact solutions are known. Then, we demonstrate the applicability of these equations by using them to compute the population dynamics and transient absorption spectrum for a multi-state double-well system.

III.1 Single-State Case: Brownian Oscillator

Here, we consider the case of a single PES described by the harmonic potential

[TABLE]

The validity of the reduced equation of motion for a non-Markovian system can be examined by comparing the results obtained from a set of numerical tests (non-Markovian tests) to the analytically derived solution in the case of the Brownian oscillator Tanimura (2015).

First, we study the equilibrium distribution function in the case of the above harmonic potential. For this PES, we have the following analytical expression for the equilibrium distribution Grabert et al. (1988); Weiss (2011):

[TABLE]

Here,

[TABLE]

is the mean square of the coordinate $q$, and

[TABLE]

is the symmetrize-correlation function of the coordinate $q$.

To obtain the thermal equilibrium state numerically, we integrated Eqs. (11) and (17) from a temporal initial state to a time sufficiently long that all of the hierarchical elements reached the steady state. For all of our computations, we fixed the oscillator frequency to $\omega_{0}=400~{}\mathrm{cm}^{-1}$. We consider the underdamped ($\zeta=0.1\,\omega_{0}$), critically-damped ($\zeta=2\,\omega_{0}$), and overdamped cases ($\zeta=10\,\omega_{0}$) at the temperature $\beta\hbar\omega_{0}=7.47$ ($T=77~{}K$), which is in the low-temperature regime.

Because both Eqs. (11) and (17) consist of sets of infinitely many differential equations, we need to truncate $\vec{n}$ to carry out numerical calculations. Here, we adopted the truncation scheme proposed in Refs. Härtle et al., 2013, 2015 with modifications: The hierarchy is truncated in accordance with the condition that $\vec{n}$ satisfies the relation $\Delta_{\vec{n}}\omega_{0}/\Gamma_{\vec{n}}>\delta_{\mathrm{tol}}$, where $\delta_{\mathrm{tol}}$ is the tolerance of the truncation, with $\Gamma_{\vec{n}}\equiv\sum_{k}^{K}n_{k}\nu_{k}$ and

[TABLE]

For details, see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems in the Supporting Information.

Numerical calculations were carried out to integrate Eqs. (11) and (17) using the fourth-order low-storage Runge-Kutta (LSRK4) method Yan (2017). The time step for the LSRK4 method was chosen between $\delta t=0.1\times 10^{-2}~{}\mathrm{fs}$ and $\delta t=0.5\times 10^{-2}~{}\mathrm{fs}$. Uniform meshes were employed to discretize the Wigner function and probability distribution function, and the mesh sizes were set to $N_{q}=64$ and $N_{p}=64$ in the $q$ and $p$ directions, respectively. The mesh ranges of the Wigner function and probability distribution function in the $q$ direction were chosen between $-4\leq q\leq+4$ and $-12\leq q\leq+12$. The mesh range of the Wigner function in the $p$ direction was chosen between $-4\leq p\leq+4$ and $-15\leq p\leq+15$. The finite difference calculations for $q$ and $p$ derivatives in Eqs. (11) and (17) were performed using the central difference method with tenth-order accuracy. For the kinetic term of the Liouvillian in Eq. (12a), the upwind difference method with ninth-order accuracy was employed for the $q$ derivative Frensley (1987); Kato and Tanimura (2013); Sakurai and Tanimura (2013, 2014). The Moyal products in Eq. (12c) were evaluated using the discretized convolution representation described in Refs. Frensley, 1987, 1990; Kim, 2007 with modifications for multi-state systems (for details of the modifications, see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information). The number of QLT correction terms was chosen from $K=3$ and $4$ for the low-temperature case ($\beta\hbar\omega_{0}=7.47$ $(77~{}\mathrm{K})$), and $K=2$ was employed for the high-temperature case ($\beta\hbar\omega_{0}=1.92$ $(300~{}\mathrm{K})$). The tolerance of the truncation was chosen between $\delta_{\mathrm{tol}}=10^{-4}$ and $\delta_{\mathrm{tol}}=10^{-6}$. In the case of $\delta_{\mathrm{tol}}=10^{-4}$, the number of total hierarchical elements were $34$ and $14$ for the low-temperature and high-temperature cases, respectively. The C++ source codes, which allow for the treatment of the phase and coordinate space dynamics with any single-state or multi-state potential forms, are provided as Supporting Information. The actual numerical integrations for the present calculations were carried out using C++/CUDA codes with cuBLAS and cuFFT libraries to enhance the computational speed with graphics processing unit (GPU).

In Fig. 2, we compare the analytic equilibrium distribution with that obtained from the LT-QFPE under several conditions for the system-bath coupling at $\beta\hbar\omega_{0}=7.47$. The results obtained from the LT-QFPE are overlapped to analytically exact solutions. In the numerical calculation, the larger number $K$ we used, the more accurate results we had.

In the case of a weak interaction, both the numerical and analytical forms are closer to that for the isolated harmonic oscillator, which is given by Eq. (21) with $\langle q^{2}\rangle_{\beta,0}=(1/2)\coth(\beta\hbar\omega_{0}/2)$. This expression can also be derived from the Boltzmann summation of the eigenstates.

As mentioned above, the mean-square momentum, $\langle p^{2}\rangle$, diverges in the present Ohmic case Grabert et al. (1988); Weiss (2011). This is because high frequency quantum noise destroys the quantum coherence between the “bra” and “ket” wave functions, which results in the condition $\rho(z,z^{\prime})=0$ for $z\neq z^{\prime}$. As a result, the Wigner distribution function in momentum space, which is the Fourier transform of the quantum coherence $r\equiv z-z^{\prime}$, is flat in this case (see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems in the Supporting Information). However, even in such situations, we can use the Wigner function, because the dynamics of the system are controlled by the low-frequency Matsubara terms or the low-frequency QLT correction terms.

In Fig. 3, we plot stable solutions for the Wigner distribution function calculated with the LT-QFPE using several values for the number of QLT correction terms. It is seen that the width of Wigner distributions ($\propto\langle p^{2}\rangle$) increases as the number of QLT correction terms increases, while the $q$ probability distributions converge to on the analytically derived solution. The Gaussian-like profile in the $p$ direction arises from QLT correction terms for finite $K$. Although we observe the larger flat distribution for larger $K$, we can still use the Wigner function by ignoring these contribution for the calculation of the non-diverging physical variables. This is the reason that we can calculate the physical variable using Wigner distribution, while $\langle p^{2}\rangle$ diverges with $K\rightarrow\infty$.

We next study the symmetrized correlation function of the system coordinate, defined as

[TABLE]

where $\mathcal{G}_{\mathrm{tot}}(t)$ is Green’s function for the total system and $W^{\mathrm{eq}}(p,q,\vec{x},\vec{x})$ is the stationary solution of $\mathcal{G}(t)$. In the HEOM formalism, the time evolution of the total system, described by $\rho_{\mathrm{tot}}(z,z^{\prime},\vec{x},\vec{x}\,^{\prime},t)$, is replaced by that of the hierarchical elements, described by $\rho_{\mathrm{H}}(z,z^{\prime},t)\equiv\{\rho_{\vec{n}}(z,z^{\prime},t)\,|\,\vec{n}\in\mathbb{N}^{K}\}$. It has been found that these yield the same reduced dynamics in the system subspace Tanimura (2014, 2015). After the Wigner transformation, Eq. (25) is evaluated as

[TABLE]

where $\mathcal{G}_{\mathrm{H}}(t)$ is Green’s function evaluated from Eq. (11) or (17), and $W_{\mathrm{H}}^{\mathrm{eq}}(p,q)$ is the stationary solution of $\mathcal{G}_{\mathrm{H}}(t)$. We define the Fourier transform of Eq. (25) as

[TABLE]

In the quantum case, Eq. (27) can be analytically evaluated as Eq. (23). Under the condition $\zeta\gg\omega_{0}$, this function asymptotically approaches Garg et al. (1985); Tanimura and Mukamel (1993)

[TABLE]

where $\tilde{\gamma}\equiv\omega_{0}^{2}/\zeta$. The classical high-temperature limit of the above result is obtained by replacing $\coth(\beta\hbar\omega/2)$ with $2/\beta\hbar\omega$.

In Fig. 4, we depict the symmetrized-correlation function calculated under several damping conditions in (i) low temperature ($\beta\hbar\omega_{0}=7.47$) and (ii) high temperature ($\beta\hbar\omega_{0}=1.92)$ cases. As seen there, the numerical results obtained from Eqs. (11) and (17) are close to the exact analytical solutions. Generally, in the HEOM formalism, we are able to obtain as accurate results as we need by employing larger hierarchical space (i.e. by increasing $K$ and by decreasing $\delta_{\mathrm{tol}}$). It should be noted that, although the LT-QSE well predicts the quantum dynamics under strong friction, the stable solution of the LT-QSE depend upon the number of QLT correction terms. As shown in Appendix C, because the analytical solution for the overdamped case, Eq. (28), also diverges under the infinite summation of the Matsubara frequencies, careful verification is important when we use the LT-QSE to calculate observable strongly depend on the equilibrium distribution.

As shown in this section, the LT-QFPE and LT-QSE can describe accurate dynamics with a properly truncated hierarchical space. This finding is important for numerical calculations. As shown in Appendix B, the LT-QFPE and LT-QSE are equivalent to the Langevin expressions for a single harmonic potential case. As analytically calculated symmetrized correlation functions from the Langevin equations indicate, we can obtain Eqs. (23) and (28) from the present approach, under $\delta_{\mathrm{tol}}\rightarrow 0$ and $K\rightarrow\infty$. This result demonstrates the reliability of the LT-QFPE and LT-QSE theories.

III.2 Multi-State Case: Double-Well PESs with Gaussian Adiabatic Coupling

Next, we present our numerical results for multi-state systems. For convenience, we describe our system using adiabatic electronic states, $|\Phi_{a}(z)\rangle$. In the following, $a$, $b$ and $c$ refer to adiabatic electronic states, and $j$, $k$ and $l$ refer to diabatic electronic states.

III.2.1 Adiabatic and diabatic bases

The $a$th adiabatic electronic state is an eigenfunction of the time-independent Schrödinger equation, and thus we have

[TABLE]

where $\hat{U}(z)\equiv\sum_{j,k}|j\rangle U^{\mathrm{d}}_{jk}(z)\langle k|$, and $U^{\mathrm{a}}_{a}(z)$ is the $a$th adiabatic Born-Oppenheimer (BO) PES. The diabatic and adiabatic states are related through the transformation matrix given by

[TABLE]

We introduce the unitary matrix $\bm{Z}(z)$ defined as $\{\bm{Z}(z)\}_{ja}=Z_{ja}(z)$, which satisfies the relation $\bm{Z}(z)^{\dagger}\bm{Z}(z)=\bm{Z}(z)\bm{Z}(z)^{\dagger}=\bm{1}$. Then, Eq. (29) can be expressed in diagonal matrix form as

[TABLE]

where $\{\bm{U}_{\mathrm{a}}(z)\}_{ab}\equiv\delta_{ab}U_{\mathrm{a}}^{a}(z)$, and $\delta_{ab}$ is the Kronecker delta.

[TABLE]

and therefore we have

[TABLE]

where $\mathop{\exp}_{\rightarrow}$ is the ordered exponential in coordinate space. Thus, the transformation matrix, Eq.(30), can be constructed from $\bm{d}(z)$, and hence the diabatic PESs can be obtained from the adiabatic PESs using the inverse of the transformation in Eq. (31).

[TABLE]

Next, we introduce the reduced density matrix in the adiabatic representation, defined as

[TABLE]

where the diagonal element $\rho_{aa}^{\mathrm{a}}(z,z,t)$ and the off-diagonal element $\rho_{ab}^{\mathrm{a}}(z,z^{\prime},t)$ $(a{\neq}b)$ represent the population of $|\Phi_{a}(z)\rangle$ and the coherence between $|\Phi_{a}(z)\rangle$ and $|\Phi_{b}(z^{\prime})\rangle$, respectively. The adiabatic representation of the density matrix, $\{\bm{\rho}^{\mathrm{a}}(z,z^{\prime},t)\}_{ab}=\rho_{ab}^{\mathrm{a}}(z,z^{\prime},t)$, can be obtained from $\bm{\rho}^{\mathrm{d}}(z,z^{\prime},t)$ through application of the transformation matrix $\bm{Z}(z)$ as

[TABLE]

This representation is related to the Wigner representation as

[TABLE]

Although we can construct the equations of motion for $\bm{W}^{\mathrm{a}}(p,q,t)$ directly, the numerical integrations are complicated, because the number of terms that include the Moyal product becomes large (see Appendix E). For this reason, we integrate the equations of motion in the diabatic representation. After obtained the numerical results, we convert these to the adiabatic representation.

III.2.2 Tilted double-well model

[TABLE]

We assume that the non-adiabatic coupling has the Gaussian form

[TABLE]

and $d_{gg}(q)=d_{ee}(q)=0$, where $\sigma^{\dagger}$ is the width of the crossing region. The integral of $d_{eg}(q)$ is given by

[TABLE]

where $\mathrm{erfc(x)}\equiv 1-\mathrm{erf}(x)$ is the complementary error function. Because $d_{eg}(q)$ is normalized with respect to $D_{eg}(\infty)=\pi/2$, the adiabatic bases are exchanged with the change of position from $z=-\infty$ to $z=+\infty$.

Although we can construct the MS-LT-QFPE and their variant equations in the adiabatic representation (see Appendix E), in the present study, we performed the numerical calculation using the diabatic representation, because in this case, the equations are simpler and easier to solve. Then, after we obtained the numerical results, we carried out the inverse transformation in order to convert these to the adiabatic representation. We employ the diabatic basis defined as $|0\rangle\equiv|\Phi_{g}(-\infty)\rangle$ and $|1\rangle\equiv|\Phi_{e}(-\infty)\rangle$. Then Eq. (33) is solved as

[TABLE]

and the diabatic PESs and coupling are given by

[TABLE]

Figure 5 presents the adiabatic BO PESs and diabatic PESs for the parameter values listed in Table 2.

We set the initial distribution as

[TABLE]

and as $W_{gg}^{\mathrm{a}}(p,q,0)=W_{eg}^{\mathrm{a}}(p,q,0)=W_{ge}^{\mathrm{a}}(p,q,0)=0$ in Eq. (43), where $q_{\mathrm{i}}=-L_{0}/2$ and $\mathcal{Z}$ is the partition function. This is the Wigner transformation of the Boltzmann distribution for the harmonic oscillator centered at $q=q_{\mathrm{i}}$. In this demonstration, we ignore the initial correlation at $t=0$.

We performed the numerical calculations to integrate equations of motion using the finite difference method with mesh sizes $N_{q}=256$ and $N_{p}=64$ and mesh ranges $-12\leq q\leq+12$ and $-12\leq p\leq+12$. The other calculation conditions were the same as in Sec. III.1. For comparison, we display the results calculated using the fewest switch surface hopping (FSSH) Tully (1990); Hammes-Schiffer and Tully (1994) and Ehrenfest methods Stock and Thoss (2005); May and Kühn (2008) with a classical Markovian Langevin force under the same conditions. In both methods, the adiabatic electronic basis was employed. In the Ehrenfest methods, the state of the system is described using the electronic density matrix (or the electronic wavefunction), $\bm{\rho}^{\mathrm{el}}(t)$, and a trajectory, $\{p(t),q(t)\}$, which is determined following the mean-field force calculated from the electronic PESs (i.e. averaged force from the populations in $\bm{\rho}^{\mathrm{el}}(t)$). In the FSSH methods, the state is also described using $\bm{\rho}^{\mathrm{el}}(t)$, but its trajectory, $\{p(t),q(t)\}$, now follows the force calculated from the active PES, $U_{\lambda(t)}^{\mathrm{a}}(q)$, where $\lambda(t)$ is a randomly changing index (i.e. $g$ or $e$) in time whose hopping rate is calculated using $\bm{\rho}^{\mathrm{el}}(t)$. In each case, the time evolution of the particle trajectories under the Langevin force was calculated using the Vanden-Eijnden-Ciccotti methods Tuckerman (2010). The time evolution of $\bm{\rho}^{\mathrm{el}}(t)$ was calculated by the numerical integration of the Schrödinger equation (NISE) method Torii (2006); Jansen and Knoester (2006), in which the coefficients of the time evolution operator of the electronic density matrix were held constant during each time step evaluation. We used $\delta t=0.1~{}\mathrm{fs}$ as the time step for the integrations using the FSSH and Ehrenfest methods, and we employed $N\simeq 10,000$ trajectories for the average calculations.

III.2.3 Population dynamics

In Fig. (6), we present plots of the excited population in the adiabatic representation defined as,

[TABLE]

In the FSSH and Ehrenfest cases, we defined $u_{e}(t)$ as an ensemble average of $\rho^{\mathrm{el}}_{ee}(t)$ for all trajectories. In the FSSH (active) case, we further introduced the excited population, $\bar{u}_{\lambda(t){=}e}(t)$, calculated as an ensemble average of the population with the active PES index $\lambda(t)=e$ for all trajectories.

As demonstrated in Sec. III.1, the MS-LT-QFPE is numerically accurate, and therefore the calculated results obtained with these equations can be used as a reference to verify other results from different formalisms. We now discuss the numerical results. A wavepacket centered at $q=q_{\mathrm{i}}$ in the excited state as the initial state. The wavepacket then moves in the direction of the crossing region, i.e. near $q=q^{\dagger}$. Then the population of $|\Phi_{e}(q)\rangle$ decreases due to the non-adiabatic transition. In the case that $\zeta$ is sufficiently small, the de-excited wavepacket moves among the double minima of the adiabatic ground state PES while maintaining a large kinetic energy, and it traverses the crossing region repeatedly through non-adiabatic transitions. As a result, oscillatory behavior is observed in Figs. (6)(i-a) and (6)(ii-a). As the coupling constant $\zeta$ increases, the wavepacket motion becomes slower, and the population of $|\Phi_{e}(q)\rangle$ decays more slowly. While classical treatments of the heat bath produce pathological negative populations (red and green curves), in particular in the low temperature case, the presently investigated equations of motion, Eqs. (11) and (17), accurately describe the population dynamics (the blue lines). This demonstrates the importance of the QLT correction terms.

While the results obtained with the FSSH and Ehrenfest methods do not exhibit negative populations because of the assumptions for a decomposition of the distribution function into the trajectories, they do differ significantly from both the quantum (MS-LT-QFPE) and semi-classical (MS-FPE) results, particularly in the low temperature case. Because the FSSH and Ehrenfest methods were originally developed to study isolated systems, which have a few degrees of freedom, these methods may not be proper for calculations of systems in dissipative conditions that are essentially many-body problem when we include the bath degrees of freedom. It should be noted that, in the FSSH results, $u_{e}(t)$ and $\bar{u}_{\lambda(t){=}e}(t)$ disagrees, while both correspond the electronic excited state. The similar disagreement is reported in Ref. Tempelaar et al., 2013.

Note that, as the coupling constant $\zeta$ increases, the results obtained from the quantum Liouvillian (red curves) and semi-classical Liouvillian (green curves) treatments of the PES become similar. This is because the higher-order differential operators in the Moyal products vanish in the overdamped limit, as shown in the construction of the LT-QSE (see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems in Supporting Information). Moreover, when the energy gap between the ground and excited states in the crossing region is small, the negative population of the semi-classical results is suppressed. This is because, the characteristic frequency of the electronic transition dynamics becomes small in such situation, and therefore the high-temperature approximation of the bath works well, as the authors’ previous investigation for conical intersection problem Ikeda and Tanimura (2018).

III.3 Transient Absorption Spectra

The presently investigated formalisms are capable of calculating non-linear response functions. By calculating the transient absorption spectrum, here we demonstrate the importance of the QLT correction terms for nonlinear optical spectra. In order to include an excited state absorption (ESA) process, we add a second excited state, $|\Phi_{f}(q)\rangle$, with the PES

[TABLE]

to the present model. Because we assume $d_{fg}(q)=d_{fe}(q)=0$, a spontaneous transition between $|\Phi_{f}(q)\rangle$ and the subspace $\{|\Phi_{g}(q)\rangle,|\Phi_{e}(q)\rangle\}$ is prohibited. Thus, only optically stimulated transitions between $|\Phi_{g}(q)\rangle$ and $|\Phi_{e}(q)\rangle$ occur. The diabatic PESs of this three-state system are those given in Eq. (42), along with $U_{22}(q)=U_{f}(q)$, $U_{20}(q)=U_{21}(q)=0$.

The transient absorption (TA) spectrum from the initial state Eq. (43) is given by Ikeda and Tanimura (2017)

[TABLE]

where

[TABLE]

is the response function of the TA and $\{\bm{\mu}^{\mathrm{d}}(q)\}_{jk}=\langle j|\hat{\mu}(q)|k\rangle$ is the dipole operator in the matrix representation. Here, $A^{\times}B\equiv AB-BA$ is the commutator, and

[TABLE]

is the trace over the system subspace. In the HEOM formalism, Eq. (47) becomes

[TABLE]

where $\bm{W}_{\mathrm{H}}^{\mathrm{d}}(p,q,0)\equiv\{\bm{W}_{\vec{n}}^{\mathrm{d}}(z,z^{\prime},t)\,|\,\vec{n}\in\mathbb{N}^{K}\}$ is initialized as $\bm{W}_{\vec{n}}^{\mathrm{d}}(p,q,0)=\bm{W}^{\mathrm{d}}(p,q,0)$ for $\vec{n}=\vec{0}$ and $\bm{W}_{\vec{n}}^{\mathrm{d}}(p,q,0)=0$ otherwise.

Hereafter, for $\hat{\mu}(q)$, we assume the form $\hat{\mu}(q)=|\Phi_{e}(q)\rangle\langle\Phi_{g}(q)|+|\Phi_{f}(q)\rangle\langle\Phi_{e}(q)|+\mathrm{c.~{}c}$. While $\hat{\mu}(q)$ may induce non-vertical transitions among adiabatic electronic states, here we consider the vertical transition only for simplicity (i.e. the laser interaction is described by the commutator, $\bm{\mu}^{\mathrm{a}}(p,q)^{\times}\bm{W}^{\mathrm{a}}(p,q)$ among the electronic states). Note that because the distribution, Eq. (43), is set in the excited state at the initial time, ground state bleaching (GSB) is not observed.

Figure 7 displays the calculated results in the overdamped case at $T=77~{}\mathrm{K}$. Here, the parameter values are the same as in Fig. 6(i-c’). For comparison, we display the TA spectra obtained using the FSSH and Ehrenfest methods with a classical Markovian Langevin force under the same conditions. As shown in Fig. 6, non-adiabatic transitions in the FSSH and Ehrenfest calculations are slow. This is reflected by very sharp peaks in the TA spectra. To avoid numerical errors in the computation of the Fourier transform due to these sharp peaks, we employed an exponential filter, $\exp(-\tau/\tau_{\mathrm{d}})$, with decay time constant $\tau_{\mathrm{d}}=200~{}\mathrm{fs}$ in the all calculations. The procedure of calculating non-linear response functions using the FSSH and Ehrenfest methods following that for the two-dimensional electronic spectra given in Ref. Tempelaar et al., 2013 is presented in Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information.

To carry out FSSH calculations, the effects of the non-adiabatic transition dynamics during waiting time was taken into account using the ensemble average of the active PES index, $\bar{u}_{\lambda(t){=}e}(t)$ (i.e. using the auxiliary wavefunction given in Ref. Tempelaar et al., 2013).

First, we discuss the TA spectrum plotted in Fig. 7(i). At $t=0$, ESA and stimulated emission (SE) are observed at frequencies corresponding to the Frank-Condon point, $U_{f}(q_{\mathrm{i}})-U_{e}(q_{\mathrm{i}})=14,000~{}\mathrm{cm}^{-1}$ and $U_{e}(q_{\mathrm{i}})-U_{g}(q_{\mathrm{i}})=10,552~{}\mathrm{cm}^{-1}$, respectively. The ESA signals is labeled “$\mathrm{A}$” and the SE signal is labeled “$\mathrm{B}$”. Because both $U_{f}(q)-U_{e}(q)$ and $U_{e}(q)-U_{g}(q)$ decrease in accordance with the wavepacket motion from $q=q_{\mathrm{i}}$ to $q=z^{\dagger}$, the ESA and SE peaks move toward the frequencies $\omega=U_{f}(q^{\dagger})-U_{e}(q^{\dagger})=E_{\mathrm{gap}}^{f-e}=5,000~{}\mathrm{cm}^{-1}$ and $U_{f}(q^{\dagger})-U_{e}(q^{\dagger})=E_{\mathrm{gap}}^{e-g}=1,000~{}\mathrm{cm}^{-1}$ as functions of time, respectively. The intensity of the ESA signal decreases due to the de-excitation from $|\Phi_{e}(q)\rangle$ to $|\Phi_{g}(q)\rangle$. The intensity of the SE signal also decreases, and becomes zero near $t=300~{}\mathrm{fs}$, because the SE from $|\Phi_{e}(q)\rangle$ and the absorption from $|\Phi_{g}(q)\rangle$ cancel each other. Then, the absorption signal from the crossing region labeled “$\mathrm{C}$” appears and the position of the peak moves toward the resonant frequency of the product ground state labeled “$\mathrm{D}$”, following the wavepacket motion.

In Fig. 7(ii), the positive peak labeled “$\mathrm{E}^{\prime}$” appears at $\omega=E_{\mathrm{gap}}^{f-e}$ after $t=400$ fs. Its appearance is due to the violation of the positivity of the excited state population of $|\Phi_{e}(q^{\dagger})\rangle$; the $\mathrm{E}^{\prime}$ peak in ESA signal arises from a negative population of $|\Phi_{e}\rangle$, which is observed as emission. The negative population of $|\Phi_{e}(q^{\dagger})\rangle$ also appears as an absorptive contribution in the SE, and the intensity of the negative peak labeled “$\mathrm{C}^{\prime}$” is enhanced. When this negative population propagates to the $|\Phi_{g}(q^{\dagger})\rangle$ state, the absorption from this state is observed as the emission spectrum labeled “$\mathrm{F}^{\prime}$”. These results indicate that the classical treatment does not have the capability to predict optical signals, in particular in the low-temperature regime.

In Figs. 7(iii) and 7(iv), the calculated results from the phenomenological FFSH and Ehrenfest approaches are presented. These results do not agree with the MS-LT-QSE results not only quantitatively but also qualitatively, in particular in the Ehrenfest case. This is due to the poor estimation of non-adiabatic transition rate in these two approaches, as illustrated in Fig. 6. Although the MS-SE approach exhibits unphysical emission or absorption peaks, the MS-SE result is closer to the accurate MS-LT-QSE results than the FSSH and Ehrenfest results. This indicates that the MS-SE approach is physically more consistent than the FSSH and Ehrenfest approaches, even the positivity problem occurs. Although the FSSH and Ehrenfest approaches are simple and easy to implement molecular dynamics simulations, the verification of the calculated results, in particular for calculation of nonlinear response function, should be made before comparing the experimental results. Note that, when we calculate the effects of the non-adiabatic transition dynamic using the ensemble average of the electronic density matrix (i.e. using the primary wavefunction given in Ref. Tempelaar et al., 2013), the FSSH method produces a similar spectrum to the Ehrenfest result.

IV CONCLUSION

In this paper, we investigated (MS-)LT-QFPE and (MS-)LT-QSE that include QLT correction terms to satisfy the QFD theorem. The (MS-)LT-QFPE and (MS-)LT-QSE were rigorously derived from path integral formalism. We found that the (MS-)LT-QFPE can be used to obtain correct numerical descriptions of dynamics of a system coupled to an Ohmic bath when the QLT correction terms are treated properly, even in the strong coupling, low-temperature regime. In the overdamped case, we can further reduce the momentum degrees of freedom from these equations, thereby obtaining the (MS-)LT-QSE. Although the applicability of these equations is limited to the Ohmic case, they are significantly less computationally intensive than (MS-)QHFPE approaches in particular for the case of (MS-)LT-QSE. Moreover, because structures of the PESs play essential roles in nonadiabatic transition phenomena, and because difference between the Markovian and non-Markovian noise cases is minor, the present formalism is sufficient for studying nonadiabatic transition phenomena. Applications of this approach to the study of a molecular motor system will be presented in forthcoming papers. The MS-LT-QFPE and MS-LT-QSE are also helpful for identifying purely quantum effects, because they allow us to compare the quantum results with the classical results obtained in the classical limit of the MS-LT-QFPE and MS-LT-QSE.

As shown in Appendix C, although the LT-QSE accurately predicts the quantum dynamics in the case of strong friction, while we must truncate the number of QLT correction terms properly with estimating the timescales of each term given by, $\nu_{k}$, in comparison with the timescale of the system dynamics. A correction of the LT-QSE using conventional QSE theories may suppress this ambiguity. This is left for future investigations.

Because the (MS-)LT-QFPE is derived using a technique similar to that used in the conventional HEOM approach, an extension of the present formalism to the imaginary time formalism in calculations of the partition function should be straightforward Tanimura (2014, 2015).

Appendix A Stochastic Liouville Description of (MS-)LT-QFPE and (MS-)LT-QSE

Because each contribution from the QLT terms is a Gaussian process, we can construct equations of motion in terms of continuous stochastic variables Tanimura and Kubo (1989); Tanimura (2006). We introduce a set of stochastic variables $\vec{\Omega}\equiv(\dots,\Omega_{k},\dots)$, where $\Omega_{k}$ is the auxiliary stochastic variable for the description of the $k$th QLT correction term. Then, the Wigner distribution function is expressed as

[TABLE]

Here, the functions $\phi_{\vec{n}}(\vec{\Omega})$ and $\phi_{\vec{n}}^{(-1)}(\vec{\Omega})$ are defined as

[TABLE]

where $\psi_{\vec{n}}(\vec{\Omega})$ is the Hermite function, $\psi_{\vec{n}}(\vec{\Omega})\equiv\prod_{k}^{K}\psi_{n_{k}}(\Omega_{k})$, with

[TABLE]

and the $n$th Hermite polynomial, $H_{n}(z)\equiv(-1)^{n}e^{z^{2}}(\partial^{n}/\partial z^{n})e^{-z^{2}}$. The coefficients $\alpha_{k}\neq 0$ and $b_{k}\neq 0$ are real numbers. Then, the MS-WDF, $\bm{W}^{\mathrm{d}}(p,q,t)=\bm{W}_{\vec{0}}^{\mathrm{d}}(p,q,t)$, can be expressed as

[TABLE]

because of the orthogonality of the Hermite functions.

While the coefficients $\alpha_{k}$ and $b_{k}$ in Eqs. (51a) and (51b) can be chosen in an arbitrary manner, we found that the choice $\alpha_{k}=\sqrt{\mathstrut 4\eta_{k}\zeta/\beta\hbar\omega_{0}\nu_{k}}$ and $b_{k}=\sqrt{\mathstrut 2}/\alpha_{k}\nu_{k}$ makes the equation of motion for $\bm{W}^{\mathrm{d}}(p,q,\vec{\Omega},t)$ simple. Thus we obtain the stochastic Liouville description of the (MS-)LT-QFPE as

[TABLE]

where

[TABLE]

For details of Eq. (54), see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information. Note that, while Eq. (54) is similar to the Fokker-Planck equation for classical non-Markovian dynamics via Markovian-type Fokker-Planck equations with “virtual variables” Ferrario and Grigolini (1979); Marchesoni and Grigolini (1983), our variables $\vec{\Omega}=(\dots,\Omega_{k},\dots)$ are introduced to describe the QLT correction terms from the Bose-Einstein distribution function. Equation (50b) is similar to the discretized representation of a phase-space distribution of the classical Kramers equation (the Brinkman hierarchy) Risken (1989); Ikeda and Tanimura (2017). Therefore, the (MS-)LT-QFPE can be regarded as the Brinkman hierarchy representation of Eq. (54).

Similarly, the (MS-)LT-QSE, Eq. (17), is equivalent to the equation of motion,

[TABLE]

where

[TABLE]

Further details of Eq. (56) are given in Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information.

Appendix B Langevin Descriptions of LT-QFPE and LT-QSE

In the case that the system has only a single-state (i.e. in the case of that $\bm{W}^{\mathrm{d}}(p,q,t)$ and $\bm{U}^{\mathrm{d}}(q)$ reduce to scaler functions, $W(p,q,t)$ and $U(q)$) and the quantum nature of the nuclear dynamics is weak (i.e. anharmonicity of the potential is weak or friction is strong), higher-order terms of the Moyal product, Eq. (13), can be omitted and Eq. (12c) is approximated as

[TABLE]

Here, we have introduced the force $F(q)\equiv-(1/\hbar)\partial U(q)/\partial q$. In a harmonic potential case, the above expression is exact. Then, Eqs. (11) and (54) can be decomposed into a set of Langevin equations,

[TABLE]

Here, $\Omega_{k}(t)$ is an auxiliary stochastic variable for the $k$th index of $\vec{n}$, and $\tilde{R}_{q}(t)$ and $\tilde{R}^{(\nu)}_{k}(t)$ are Gaussian-white forces that satisfy the relations,

[TABLE]

and $\langle\tilde{R}_{q}(t)\tilde{R}_{k}^{(\nu)}(t^{\prime})\rangle=\langle\tilde{R}_{k}^{(\nu)}(t)\tilde{R}_{k}^{(\nu)}(t^{\prime})\rangle=0$ ($k\neq k^{\prime}$). For details of Eqs. (59a)–(59c), see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information. By introducing a quantum random force as

[TABLE]

the set of Langevin equations, Eqs. (59a), (59b), and (59c), can be rewritten as

[TABLE]

The random force $\tilde{R}_{\mathrm{qm}}(t)$ satisfies the QFD theorem as

[TABLE]

The set of equations (62a) and (62b) are the quantum Langevin equation Ford et al. (1988); Weiss (2011) for the c-number variables, $p(t)$ and $q(t)$ (i.e. the quasi-classical Langevin equation Senitzky (1960); Schmid (1982)). In the cases that anharmonicity of the system is strong, we cannot employ the above equations. In such cases, we have to evaluate the quantum Langevin equation described by the operators, $\hat{p}(t)$ and $\hat{q}(t)$, or to employ Wigner description (i.e. Eqs. (11) and (54)) with the quantum Liouvillian (12a). This indicates that, the LT-QFPE, Eqs. (11) and (54), can be regarded as the Fokker-Planck equations equivalent to the quantum Langevin equation. Note that, while Eqs. (59a), (59b), and (59c) are similar to the generalized Langevin equation for classical non-Markovian dynamics via Markovian-type Lanvegin equations with virtual variables Ferrario and Grigolini (1979); Marchesoni and Grigolini (1983); Martens (2002); Kawai (2015), our stochastic variables $\vec{\Omega}=(\dots,\Omega_{k},\dots)$ are introduced to describe the QLT correction terms from the Bose-Einstein distribution function.

It should be noted that, for the LT-QFPE with a single PES, Eq. (17) can also be decomposed into a set of Langevin equations,

[TABLE]

By introducing $\tilde{R}_{\mathrm{qm}}(t)$, the set of Langevin equations, Eqs. (64a) and (64b), can be rewritten as

[TABLE]

This is the “overdamped” quantum Langevin equation for the c-number variable, $q(t)$ (i.e. the inertia term $\propto\ddot{q}(t)$ is omitted). For details of Eqs. (64a) and (64b), see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multistate Systems of the Supporting Information.

Appendix C Comparison of QSE and LT-QSE

In this appendix, we compare our LT-QSE theory with conventional QSE theories.

The QSE is proposed in Refs. Pechukas et al., 2000; Ankerhold et al., 2001; Ankerhold and Grabert, 2008. It is given by

[TABLE]

in terms of our dimensionless coordinate $q$. Here, we have

[TABLE]

In the case of the harmonic oscillator (20), Eq. (66) can be written

[TABLE]

Better quantum corrections for the QSE are presented in Refs. Maier and Ankerhold, 2010a, b. For a harmonic potential, this becomes

[TABLE]

where $\Omega\equiv\zeta/2-\sqrt{\mathstrut(\zeta/2)^{2}-\omega_{0}^{2}}$, and $\langle q^{2}\rangle_{\beta,\zeta}$ is given in Eq. (22b). Hereafter, we refer to Eqs. (68) and (69) as the QSE and QSE’, respectively.

In Fig. 8, we display the steady-state solutions of the QSE, QSE’, and LT-QSE for several values of the damping strength at low temperature, $\beta\hbar\omega_{0}=7.47$. The other parameter values are the same as in Sec. III.1. In the overdamped case, all calculated results are qualitatively similar to the analytical result. This is because the QSE and QSE’ are constructed so as to reproduce the steady-state solution. The QSE’, in particular, utilizes $\langle q^{2}\rangle_{\beta,\zeta}$, for this reason it reproduces the analytical result even in the critical-damping case.

Note that if we assume the overdamped approximation Eq. (28) for the analytical equilibrium distribution Eq. (21) instead of Eq. (23), the integral in Eq. (22b) diverges under the infinite summation of the Matsubara frequencies as

[TABLE]

The calculated results using the LT-QSE theory become close to the above analytical result by employing larger hierarchical space (i.e. by increasing $K$). Thus, the steady-state solution predicted by the LT-QSE theory deviates as $K$ increases. This divergence is similar to the ultra-violet divergence of $\langle p^{2}\rangle$ . This indicates that, in order to calculate physical quantities on the basis of LT-QSE, we must truncate the number of QLT correction terms properly by estimating the timescales of them, $\nu_{k}$, in comparison with the timescale of the system dynamics. Note that a Drude spectral density model for reduced electronic states has also this problem, because Eq. (28) is equivalent to the symmetrized correlation function of the collective noise coordinate in the Drude spectral density model.

In Fig. 9, we display the symmetrized correlation functions calculated with the QSE, QSE’ and LT-QSE under the same conditions as in the case of Fig. 8. Because the QSE and QSE’ attempt to account for all quantum corrections with Markovian terms ignoring the non-Markovian nature of the quantum noise, they cannot reproduce the dynamics at low temperature. Contrastingly, the calculated results with the LT-QSE theory become close to the exact solution in the overdamped case, Eq. (23), by increasing $K$, and both numerical results and Eq. (23) approach the solution Eq. (28) for large $\zeta$.

Appendix D Zusman Equation

In this appendix, we show the relation between the present equation and Zusman equation. For simplicity, we assume a two-level system described by the diabatic harmonic PESs as

[TABLE]

with matrices

[TABLE]

Here, $V$ is the diabatic coupling constant between $|0\rangle$ and $|1\rangle$, $d$ is the dimensionless displacement between the minima of two PESs, and $E_{0}\equiv\Delta E+\lambda$ is the summation of the driving force $\Delta E$ and the reorganization energy $\lambda\equiv\hbar\omega_{0}d^{2}/2$. In this case, the MS-SE can be rewritten as

[TABLE]

Here, we have introduced commutation and anti-commutation hyper operators, $\bm{O}^{\times/\circ}\bm{f}\equiv\bm{O}\bm{f}\mp\bm{f}\bm{O}$, and

[TABLE]

Note that Eq. (73) is equivalent to the stochastic Liouville representation of the HEOM for a Drude spectral density in the high-temperature limit, given in Refs. Tanimura and Kubo, 1989; Shi et al., 2009. When we neglect the electronic-nuclear interaction from the force term (i.e. $\bm{F}^{\mathrm{d}}(q)\rightarrow-\omega_{0}q$ and the anti-commutation term ($\propto\bm{B}^{\circ}$) is omitted from Eq. (73)), we obtain the ZE in a matrix representation as

[TABLE]

The equations of motion for diabatic matrix elements are expressed as

[TABLE]

The above equations are equivalent to the original ZE in Ref. Zusman, 1980 except shift and scaling factors of $q$. As this derivation indicates, the MS-SE can be regarded as a generalization of the ZE for arbitrary PESs with including the quantum dynamical effects arising from the electric states. Note that, the ZE for harmonic PESs including this dynamical effect is given in Refs. Garg et al., 1985; Shi et al., 2009. This is also sometimes referred to as the ZE.

Appendix E Adiabatic Representations of MS-LT-QFPE and MS-LT-QSE

In this appendix, we present the adiabatic representations of MS-LT-QFPE (11) and MS-LT-QSE (17).

We can rewrite Eqs. (36) and (37) in terms of the Moyal star product, Eq. (13), as

[TABLE]

Inserting this into Eq. (11), we obtain the MS-QFP-LT in the adiabatic representation as

[TABLE]

[TABLE]

Here, we have introduced the commutation and anti-commutation hyper operators with the Moyal product, $\bm{O}^{\times/\circ_{(\star)}}\bm{f}\equiv\bm{O}\star\bm{f}\mp\bm{f}\star\bm{O}$. The non-Markovian noise terms are $\hat{\Phi}^{\mathrm{a}}(p,q)\equiv\hat{\Phi}^{\mathrm{d}}(p,q)$ and $\hat{\Theta}_{k}^{\mathrm{a}}(p,q)\equiv\hat{\Theta}_{k}^{\mathrm{d}}(p,q)$, and the Markovian noise term, Eq. (15c), becomes

[TABLE]

In the Smoluchowski limit, we introduce the adiabatic representation of $\bm{f}^{\mathrm{d}}(q)$:

[TABLE]

Here, we have omitted the higher-order contributions from the Moyal product in Eq. (77), because such contributions from the quantum coherence, $z-z^{\prime}$, have been removed in the Smoluchowski limit. Then the MS-LT-QSE (17) becomes

[TABLE]

where the operators appearing in Eq. (17) are transformed as

[TABLE]

[TABLE]

Here, we have introduced the force acting on the adiabatic states,

[TABLE]

and its derivative,

[TABLE]

Acknowledgements.

This work was supported by JSPS KAKENHI Grant Number JP26248005.

  

Supporting Information Available

ct8b01195_si_001.zip

The C++ source codes that we developed, which allow for the treatment of the phase and coordinate space dynamics with any single-state or multi-state potential forms.

ct8b01195_si_002.pdf

The derivations of the (MS-)LT-QFPE and (MS-)LT-QSE, the derivation of the stochastic Liouville and Langevin descriptions of (MS-)LT-QFPE and (MS-)LT-QSE, the truncation scheme of hierarchy, note on the ultra-violent divergence, and the details of the numerical calculations.

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b01195.