Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations andTheir Extension to Multi-State Systems

# Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multi-State Systems

Tatsushi Ikeda Department of Chemistry, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan    Yoshitaka Tanimura Department of Chemistry, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan
July 20, 2019
###### Abstract

Simulating ultrafast 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.

## 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 coupled to the nuclear coordinates. For simplicity, we represent the nuclear coordinates by a single dimensionless coordinate, . Here and hereafter, we employ a dimensionless coordinate and a dimensionless momentum defined in terms of the actual coordinate and momentum, and , as and , where is the characteristic vibrational frequency of the system and is the effective mass. The reaction coordinate is also bilinearly coupled to the harmonic bath coordinates, . The Hamiltonian of the total system is expressed as Caldeira and Leggett (1983)

 ^Htot(p,q;→p,→x) ≡^H(p,q)+^HB(→p,→x;q), (1)

where the system Hamiltonian, , is defined as

 ^H(p,q)≡ℏω02^p2+∑j,k|j⟩Udjk(^q)⟨k|. (2)

Here, the nuclear and electronic operators are denoted by hats, and the direct products with the unit operator in the kinetic and bath terms () are omitted. The diagonal element is the diabatic PES of , and the off-diagonal element represents the diabatic coupling between and . The vibrational frequency at a local minimum of the potential is determined by the curvature of the PESs as

 ℏ¯ω≃√ℏ(ω0∂2∂q2Uj0j0(q)∣∣∣q=q0, (3)

where is a primary state of the vibrational dynamics. Therefore, the frequency is chosen to be in order to have . The bath Hamiltonian is defined as

 ^HB(→p,→x;q) ≡∑ξℏωξ2[^p2ξ+(^xξ−gξωξ^q)2], (4)

where , , and are the vibrational frequency, conjugate momentum, and system-bath coupling constant of the th dimensionless bath mode, .

 The bath is characterized by the dissipation and fluctuation that it engenders. These are represented by the relaxation function R(t)=2π∫∞0dωJ(ω)ωcosωt (5a) and the symmetrized correlation function C(t)=2π∫∞0dωJ(ω)(n(ω)+12)cosωt, (5b) where the spectral density is defined as J(ω)≡π∑ξ(g2ξ/2)δ(ω−ωξ), and we have introduced the Bose-Einstein distribution function, n(ω)≡(eβℏω−1)−1, for the inverse temperature divided by the Boltzmann constant, β≡1/kBT.

We choose the coefficients and so as to realize the relation

 n(ω)+12≃1βℏ1ω+K∑k2ηkβℏωω2+ν2k (6)

for finite , 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 and , where is the th Matsubara frequency Tanimura (1990); Ishizaki and Tanimura (2005). In this paper, we employ the Padé spectral decomposition (PSD) 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

 J(ω) =ζω0ω, (7)

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

 R(t)=ζω0⋅2δ(t) (8a) and C(t)≃CK(t)=ζω0(1βℏ+K∑k2ηkβℏ)⋅2δ(t)−K∑kζω02ηkβℏ⋅νke−νk|t|. (8b)

In the case of a harmonic PES with frequency , the conditions , , and correspond to the underdamped, critically damped, and overdamped cases, respectively.

In Fig. 1, we plot for various values of the cutoff , at the temperature . 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, , 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, , 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, , by simply ignoring the contribution from the high-frequency QLT correction terms by implementing the cutoff , 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, , where and 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

 ρdjk(z,z′,t) ≡⟨j|^ρ(z,z′,t)|k⟩, (9)

where is the reduced density operator and represents the trace operation in the bath subspace. The diagonal and off-diagonal elements, and (), represent the population of and the coherence between and , respectively. Hereafter, we employ the matrix forms of the reduced density matrix and the diabatic PESs: and .

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)

 Wd(p,q,t) (10)

where and . Both and are now c-numbers in this phase space representation.

The reduced dynamics of and 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 Multi-State Systems in the Supporting Information). In the present case, these equations are the following:

 ∂∂tWd→n(p,q,t)=−(Ldqm(p,q)+K∑knkνk+^ΞdK(p,q))Wd→n(p,q,t)−K∑k^Φ(p,q)Wd→n+→ek(p,q,t)−K∑knkνk^Θk(p,q)Wd→n−→ek(p,q,t), (11)

where is a -dimensional multi-index whose components are all non-negative integers and is the th unit vector. The multi-index represents the index of the hierarchy, and physically, the first hierarchical element, , corresponds to the MS-WDF, . The rest of the hierarchical elements serve only to facilitate treatment of the non-Markovian system-bath interaction that arises from the QLT effects.

 The quantum Liouvillian for the MS-WDF is given by Ldqm(p,q) ≡K(p,q)+Udqm(p,q), (12a) where K(p,q)W(p,q)≡ω0p∂∂qW(p,q) (12b) and Udqm(p,q)W(p,q)≡iℏ(Ud(q)⋆W(p,q)−W(p,q)⋆Ud(q)) (12c) are the kinetic and potential terms in the diabatic representation, respectively.

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

 ⋆ (13)

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

 \lx@underaccentset→∂xf(x)=f(x)\lx@underaccentset←∂x ≡∂f(x)∂x. (14)
 The operators for the fluctuation and dissipation, ^Φ(p,q), ^Θk(p,q), and ^Ξdk(p,q), appearing in Eq. (11), are defined as ^Φ(p,q) ≡−∂∂p, (15a) ^Θk(p,q) ≡ζω02ηkβℏ∂∂p, (15b) and (15c) The first two operators in the above equations, ^Φ(p,q) and ^Θk(p,q), arises from Eq. (8a) and the first term in Eq. (8b), while the last operator, ^ΞdK(p,q), arises from the second term in Eq. (8b).

The derivation of Eq. (11) is presented in Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multi-State 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 and reduce to scalar functions, and . 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, , where 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 ()-dimensional multi-index (i.e. the additional index ) to describe non-Markovian dynamics caused by a finite value of , and therefore computationally more expensive than the (MS-)LT-QFPE. Moreover, the (MS-)QHFPE become unstable in the Ohmic limit (i.e. ) at low temperatures due to the fast decaying terms with , 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 a 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 , where is the characteristic frequency of the electronic transition dynamics. We introduce the following probability distribution in coordinate space:

 fd(q,t) ≡∫dpWd(p,q,t). (16)

In the Smoluchowski limit, the equations of motion for are

 ∂∂tfd→n(q,t)=−[Ed(q)+K∑knkνk+ω0ζ(Fd(q)+^Ξod,dK(q))]fd→n(q,t)−K∑k^Φod,d(q)fd→n+→ek(q,t)−ω0ζK∑knkνk^Θod,dk(q)fd→n−→ek(q,t), (17)
 where Ed(q)f(q,t) ≡iℏ(Ud(q)f(q,t)−f(q,t)Ud(q)) (18a) corresponds to the Liouville-von Neumann equation for the electronic subspace and Fd(q)f(q,t) ≡∂∂q12(Fd(q)f(q,t)+f(q,t)Fd(q)) (18b) is the drift term that arises from the force Fd(q)≡−(1/ℏ)∂Ud(q)/∂q.
 The operators ^Φod,d(q) =−∂∂q (19a) and ^Θod,dk(q) =2ηkβℏ∂∂q (19b) represent the non-Markovian parts of the noise, while ^Ξod,dK(q) =−1βℏ∂2∂q2+K∑k^Φod,d(q)^Θod,dk(q) (19c) represent the Markovian part of the noise. The superscript “od” means “overdamped”.

The derivation of Eq. (17) is given in Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multi-State Systems of the Supporting Information. In the case of a single-state system, the matrices , and reduce to scalar functions, , and , 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

 U(q) =ℏω02q2. (20)

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):

 f(q) =1√(2π⟨q2⟩β,ζe−q2/2⟨q2⟩β,ζ. (21)

Here,

 ⟨q2⟩β,ζ =2π∫∞0dω~Cq(ω) (22a) =ω0βℏ∞∑k=−∞1ω20+|~νk|ζ+~νk2 (22b)

is the mean square of the coordinate , and

 ~Cq(ω) =12coth(βℏω2)1ω0ζω20ω(ω20−ω2)2+ζ2ω2 (23)

is the symmetrize-correlation function of the coordinate .

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 . We consider the underdamped (), critically-damped (), and overdamped cases () at the temperature (), 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 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 satisfies the relation , where is the tolerance of the truncation, with and

 Δ→n ≡K∏k1nk!(ηkηK)nk. (24)

For details, see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multi-State 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 and . Uniform meshes were employed to discretize the Wigner function and probability distribution function, and the mesh sizes were set to and in the and directions, respectively. The mesh ranges of the Wigner function and probability distribution function in the direction were chosen between and . The mesh range of the Wigner function in the direction was chosen between and . The finite difference calculations for and 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 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 Multi-State Systems of the Supporting Information). The number of QLT correction terms was chosen from and for the low-temperature case ( ), and was employed for the high-temperature case ( ). The tolerance of the truncation was chosen between and . In the case of , the number of total hierarchical elements were and 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 . The results obtained from the LT-QFPE are overlapped to analytically exact solutions. In the numerical calculation, the larger number 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 . This expression can also be derived from the Boltzmann summation of the eigenstates.

As mentioned above, the mean-square momentum, , 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 for . As a result, the Wigner distribution function in momentum space, which is the Fourier transform of the quantum coherence , is flat in this case (see Sec. {NoHyper}Low-Temperature Quantum Fokker-Planck and Smoluchowski Equations and Their Extension to Multi-State 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 () increases as the number of QLT correction terms increases, while the probability distributions converge to on the analytically derived solution. The Gaussian-like profile in the direction arises from QLT correction terms for finite . Although we observe the larger flat distribution for larger , 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 diverges with .

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

 Cq(t) ≡Tr{zGtot(t)z+z′2ρeqtot(z,z′,→x,→x′)}, (25)

where is Green’s function for the total system and is the stationary solution of . In the HEOM formalism, the time evolution of the total system, described by , is replaced by that of the hierarchical elements, described by . 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

 Cq(t) =∫dp∫dqq{GH(t)qWeqH(p,q)}∣∣∣→n=→0, (26)

where is Green’s function evaluated from Eq. (11) or (17), and is the stationary solution of . We define the Fourier transform of Eq. (25) as

 Cq(ω) ≡∫∞0dtCq(t)cosωt. (27)

In the quantum case, Eq. (27) can be analytically evaluated as Eq. (23). Under the condition , this function asymptotically approaches Garg et al. (1985); Tanimura and Mukamel (1993)

 ~Cζ≫ω0q(ω) =12coth(βℏω2)1ω0~γωω2+~γ2, (28)

where . The classical high-temperature limit of the above result is obtained by replacing with .

In Fig. 4, we depict the symmetrized-correlation function calculated under several damping conditions in (i) low temperature () and (ii) high temperature ( 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 and by decreasing ). 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 and . 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, . In the following, , and refer to adiabatic electronic states, and , and refer to diabatic electronic states.

#### iii.2.1 Adiabatic and diabatic bases

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

 ^U(z)|Φa(z)⟩ =Uaa(z)|Φa(z)⟩, (29)

where , and is the th adiabatic Born-Oppenheimer (BO) PES. The diabatic and adiabatic states are related through the transformation matrix given by

 Zja(z)≡⟨j|Φa(z)⟩. (30)

We introduce the unitary matrix defined as , which satisfies the relation . Then, Eq. (29) can be expressed in diagonal matrix form as

 Z(z)†Ud(z)Z(q) =Ua(z), (31)

where , and is the Kronecker delta.

 In adiabatic representations of kinetic equations, non-adiabatic couplings between adiabatic states are characterized by the non-adiabatic coupling matrix, d(z), expressed in terms of the first-order derivative of the coordinate as Stock and Thoss (2005); May and Kühn (2008) {d(z)}ab=dab(z) ≡⟨Φa(z)|∂∂z|Φb(z)⟩. (32a) The non-adiabatic coupling matrix, d, is skew-Hermitian (i.e. d†=−d), and therefore the diagonal elements are zero. This matrix can also be expressed in terms of Z(z) as d(z) =Z(z)†∂∂zZ(z), (32b)

and therefore we have

 Z(z) =Z(−∞)exp→(∫z−∞dz′d(z′)), (33)

where is the ordered exponential in coordinate space. Thus, the transformation matrix, Eq.(30), can be constructed from , and hence the diabatic PESs can be obtained from the adiabatic PESs using the inverse of the transformation in Eq. (31).

 If necessary, we can introduce the non-adiabatic coupling matrix of the second-order, defined as {h(z)}ab=hab(z) ≡⟨Φa(z)|∂2∂z2|Φb(z)⟩, (34a) which can be constructed from d(z) as h(z) =Z(z)†∂2∂z2Z(z)=∂d(z)∂z+d(q)2. (34b)

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

 ρaab(z,z′,t) ≡⟨Ψa(z)|^ρ(z,z′,t)|Ψb(z′)⟩, (35)

where the diagonal element and the off-diagonal element represent the population of and the coherence between and , respectively. The adiabatic representation of the density matrix, , can be obtained from through application of the transformation matrix as

 ρa(z,z′,t) =Z(z)†ρd(z,z′,t)Z(z′). (36)

This representation is related to the Wigner representation as

 Wa(p,q,t) (37)

Although we can construct the equations of motion for 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

 As a schematic model for IC in a photoisomerization process, we adopt the following tilted double-well adiabatic ground BO PES: Uag(q) =ℏω02L20q2(q2−L202)+ΔEL0q. (38a) Here, L0 and ΔE are the displacement between the wells and the difference between their energies, respectively. We use the following harmonic adiabatic excited BO PES: Uae(q) =ℏω2e2ω0(q−q†)2+Uag(q†)+Ee−ggap. (38b) Here, ωe, q† and Egape−g are the vibrational characteristic frequency in the excited state, the position of the crossing region, and the energy gap between the ground and excited BO PES in the crossing region, respectively.

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

 deg(q) =−dge(q)=√(π8σ†2e−(q−q†)2/2σ†2, (39)

where is the width of the crossing region. The integral of is given by

 Deg(q)≡∫q−∞dq′deg(q′)=π4erfc(−q−q†√(2σ†), (40)

where is the complementary error function. Because is normalized with respect to , the adiabatic bases are exchanged with the change of position from to .

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 and . Then Eq. (33) is solved as

 Z(q)=(cosDeg(q)−sinDeg(q)+sinDeg(q)cosDeg(q)), (41)

and the diabatic PESs and coupling are given by

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Ud00(q)=Uae(q)+Uag(q)2−cos(2Deg(q))Uae(q)−Uag(q)2Ud11(q)=Uae(q)+Uag(q)2+cos(2Deg(q))Uae(q)−Uag(q)2Ud10(q)=Ud01(q)=−sin(2Deg(q))Uae(q)−Uag(q)2. (42)

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

We set the initial distribution as

 Waee(p,q,0)=1Ze−tanh(βℏω0/2)[p2+(q−qi)2], (43)

and as in Eq. (43), where and is the partition function. This is the Wigner transformation of the Boltzmann distribution for the harmonic oscillator centered at . In this demonstration, we ignore the initial correlation at .

We performed the numerical calculations to integrate equations of motion using the finite difference method with mesh sizes and and mesh ranges and . 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), , and a trajectory, , which is determined following the mean-field force calculated from the electronic PESs (i.e. averaged force from the populations in ). In the FSSH methods, the state is also described using , but its trajectory,