A new method to account for the difference between classical and quantum baths in quantum dissipative dynamics

We investigate the difference between a quantum and classical bath within the Caldeira–Leggett model for dissipative quantum dynamics. It is well known that a Markovian equation of motion for the reduced dynamics can be derived by taking the classical approximation of the harmonic bath correlation function. However, such approximation is only valid at high temperatures, and it is necessary to include the non-Markovian effect of the quantum bath in more general cases. We show that the equation of motion derived for the classical bath can be extended to the exact quantum one, by simply adding a real stochastic process to take account of the difference between the quantum and classical bath correlation functions. Numerical examples in calculating electron and excitation energy transfer dynamics, as well as absorption spectra of molecular aggregates indicate that the proposed method is a valid approach to extend the existing theories to include the quantum effect of the harmonic bath. The possibility of applying a similar idea to account for the difference between zero and high temperature quantum dissipative dynamics is discussed.


Introduction
Quantum dissipative dynamics plays an important role in many physical and chemical phenomena [1][2][3]. Recent problems of interest include coherent excitation energy transfer (EET) processes in photosynthetic systems [4,5], and conjugate polymers [6]. Theoretically, a classical system coupled to a dissipative bath can be described using well established theories such as the generalized Langevin equation [7], or the Fokker-Planck equation [7,8]. However, the problem of dissipative quantum dynamics is more complicated [2,9,10].
In literature, dynamics of a quantum system coupled to a bath is often treated with generalized quantum master equations derived using second order perturbation theory, with or without further Markovian approximations [11][12][13]. However, as the second order perturbation methods are derived under the assumption of weak system-bath coupling, they become invalid in the so called intermediate coupling regime, where the intermolecular couplings and the system-bath couplings are of similar strength [14,15]. For example, it is found that the second order perturbation theories may give incorrect energy transfer rates and spectra [16,17].
Non-perturbative equations of motion for quantum dissipative dynamics also exist. The Caldeira-Leggett equation [18] (also known as the quantum Fokker-Planck equation) was derived for a quantum system coupled to an Ohmic bath, using the high temperature approximation for the bath correlation function. As the Ohmic bath without a cutoff is unphysical and may cause technical difficulties in the quantum regime, Garg et al [19] have derived equations of motion for a quantum systems coupled to a harmonic oscillator which then couples to an Ohmic bath. The system-oscillator-Ohmic bath model is found to be equivalent to a system coupled directly to a bath with the Brownian oscillator spectral density, and again, a Markovian equation of motion can be derived by applying the classical approximation for the bath correlation function [19]. In the over damped limit, the Brownian oscillator spectral density reduces to the Debye spectral density, and the equation of motion is found to be the same as the mixed quantum-classical Zusman equation originally proposed to study electron transfer reactions with an over damped reaction coordinate [20,21].
The classical approximation for the bath correlation function is only valid in the case of high temperature. Its limitation in describing the quantum dissipative dynamics are now well known, which includes violation of the von Neumann condition [22,23], unphysical reaction rates [24][25][26] and inaccurate peaks in spectroscopy calculations [17]. It is further shown that, the criterion for the classical bath approximation being valid is given byhβ 1/ω S , where ω S is determined by the system characteristic frequency and the reorganization energy [24,26]. Hence, equations of motion derived with the classical bath approximation can be rather limited in many applications.
Modification of the quantum Fokker-Planck and Zusman equations have been explored to extend their range of validity beyond the classical bath approximation [10,23,[27][28][29][30][31]. Although different schemes are proposed in recent studies, a complete solution of the existing problems were not yet found [10,23]. An important reason could be that, all these methods have retained the Markovian nature of the original equation, while the quantum bath correlation functions are no longer Gaussian Markovian although their classical counterparts are [24,32]. It is thus important to account for the non-Markovian nature of the quantum bath to recover the correct quantum dissipative dynamics.
The recently derived hierarchical equations of motion (HEOM) method [33,34] have provided a way to take into account the effects of the quantum bath. Recently, we have shown that, the high temperature limit of the HEOM method with the Debye spectral densities leads to the Zusman equation [35]. While in the full version of the HEOM method, the non-Markovian effect of the quantum bath are included by adding the so-called low temperature correction terms to the equations of motion, which originates from the Matsubara terms in an exponential expansion of the quantum bath correlation functions [34,36]. In related works, Shao and co-workers have derived HEOM equations from the stochastic formulation to deal with the dynamics of quantum open systems [37], and devised a mixed random-deterministic method [38,39] to represent the whole or part of the bath correlation function using stochastic processes in combination with the HEOM treatment.
In this paper, inspired by the above developments in the HEOM studies, we investigate a new method to extend the original quantum Fokker-Planck and Zusman equations to account for the effect of the quantum bath. It is shown that, by adding a real Gaussian stochastic process to account for the difference between the quantum and classical bath correlation functions, the resulted stochastic quantum Fokker-Planck and Zusman equations provide a rigorous description of the exact quantum dissipative dynamics. The idea of the combined stochastic-deterministic description of quantum dissipative dynamics has been investigated previously [34,38,39], but it is the first time that existing approximate theories derived with the classical bath approximation are shown to become exact when proper stochastic fluctuations are added. We then show that this new approach can be applied to simulations of electron and EET dynamics, as well as absorption spectra of molecular aggregates. Finally, the possibility of applying a similar idea to account for the difference between high temperature and zero temperature quantum dissipative dynamics is also discussed.
The remaining sections of this paper are arranged as follows. The total Hamiltonian of the Caldeira-Leggett model, derivation of the stochastic quantum Fokker-Planck and Zusman equations are presented in section 2. We also propose an efficient method to generate the stochastic forces in the numerical simulation. In section 3, the new method is tested in applications to quantum dissipative dynamics in model systems, including spin-boson model for electron transfer reactions, absorption spectra of molecular aggregates and the EET process in the FMO complex. Conclusions and discussions are made in section 4.

Model Hamiltonian and existing theories
We start from the Caldeira-Leggett model, which considers a quantum system coupled to a bath of harmonic oscillators [1,2,18]. The total Hamiltonian including the system and bath degrees of freedom can be written aŝ whereĤ S is the system Hamiltonian,Ĥ B is the Hamiltonian of the harmonic bath, andĤ BS couples the system and bath degrees of freedom, which is assumed to be linear in the bath coordinates, In equation (3),q s is a system operator,F = i c ixi is defined as the collective bath coordinate that couples to the system degree of freedom.
The system-bath coupling is often characterized by the spectral density J (ω), which is defined as The bath correlation function C(t) ≡ F (t)F(0) can be calculated using J (ω) as Dynamics of the reduced density matrix for the system degrees of freedom is calculated using the path integral formalism. For a product initial state ρ T (0) = ρ(0) ⊗ e −βĤ B /Tr e −βĤ B , it is found that [2,9,40] In the above equation (6), q ± s (τ ) are the forward and backward paths, S ± [q ± s (τ )] are the associated system actions and where q × s (s) = q + s (s) − q − s (s) and q • s (s) = q + s (s) + q − s (s). α R (t) and α I (t) are the real and imaginary parts of the bath correlation function in equation (5): As stated in the introduction, the current study employs the system-harmonic oscillator-Ohmic bath model in [19], rather than the original system-Ohmic bath model in [18]. This model is found to be equivalent to a harmonic bath with the Brownian oscillator spectral density [19] where It was shown that the equation of motion in high temperature limit can be written as [9,19] where q is the collective bath coordinate defined as q = 1/C j c j x j , and p is the conjugate momentum; ρ w S ( p, q; t) is the partial Wigner transform in the collective bath coordinate. The Fokker-Planck operator L FK in equation (11) is defined as For simplicity, in the following derivations, we consider a harmonic oscillator bath with the Debye spectral density, which corresponds to the over-damped limit of the Brownian oscillator spectral density in equation (10) [19], while extensions to the general Brownian spectral densities is straightforward.
To this end, the Debye spectral density is given by where the cut-off frequency ω c = ω 2 0 /γ . The equation of motion in the high temperature limit is given by [ where Lρ =h −1 [H S , ρ]. The above equation is the same as the mixed quantum-classical Zusman equation [20,21] originally proposed to study electron transfer reactions with an over damped reaction coordinate. Equations (11) and (14) are non-perturbative since the collective bath mode is included explicitly in the equation of motion. They are also Markovian because that, in the high temperate limit, the bath correlation functions are Gaussian Markovian [19].
The HEOM method provides a way to go beyond the high temperature/classical approximation. First, an exponential summation is used to expand the quantum bath correlation function C(t), with ν 0 = ω c and ν k = 2kπ/(hβ); k 1 are the Matsubara frequencies. For the Debye spectral density in equation (13), and The equations of motion are then written as [34] d dt The subscript n denotes the set of indices n ≡ {n k } = {n 0 , n 1 , . . .} and n ± k differs from n only by changing the specified n k to n k ± 1. The ρ 0 with 0 = {0} = {0, . . . , 0, . . .} is the system reduced density operator while the others ρ n are all auxiliary density operators.
In the high temperature limit, by dropping the Matsubara terms and take the classical limit of the bath correlation function C(t), the full HEOM in equation (18) then reduces to It has been shown [26] that the above equation is a basis expansion of the Zusman equation in equation (14), and the relation between the distribution function ρ(q, t) in equation (14) and the auxiliary density operators ρ n in equation (20) is given by where X = ω 0 q/ √ k B T , and φ har n is the normalized harmonic oscillator eigenfunctions for the nth eigenstate.

The stochastic Zusman and quantum Fokker-Planck equations
In the HEOM method, effects of the quantum bath are included by adding the Matsubara terms, which makes the equations of motion more complex. In this subsection, we will derive the stochastic Zusman equation as an alternative way to account for the effect of the quantum bath. The idea is to decompose the influence functional in equation (7) into a stochastic part to account for the difference between the quantum and classical bath correlation functions, and a remaining classical bath correlation part that could be treated using an equation of the HEOM type. This is achieved by using the Hubbard-Stratonovich transformation [41,42]. By following the similar derivations to unravel the real part of the bath correlation function in the full path integral expression in equation (6) [34,43,44], we use a stochastic process (t) linearly coupled to the system operatorq s for unraveling part of the bath correlation function. The evolution of the reduced density operator ρ in equation (6) can then be rewritten as ] is the action caused by the stochastic process (τ ); P[ (τ )] is a Gaussian functional for (τ ) determined by its correlation function; ; t] is a to be determined residual influence functional, such that the path integral calculated via equation (22) is the same as the original result in equation (6). In this way, the influence functional in equation (7) is decomposed into a stochastic part originated from the integration over (t), plus a residual part F . This will guarantee that the stochastic Zusman equation will give the correct result for the quantum bath.
We now proceed to determine (t) and F . If the correlation function of the real Gaussian stochastic process (τ ) is defined as it can be easily shown that the residual influence functional F should take the following form [34,43,44]: such that ρ(t) in equation (22) is equal to that in equation (6). Similar to the derivation of the HEOM equation, the 'residual correlation function' α R (t) − α(t) + iα I (t) in F can be decomposed using a sum of exponential terms, and treated using the HEOM approach. Since (t) enters the equations as a stochastic external force, this will lead to a stochastic HEOM equation. For example, if all the real part of influence functional is unraveled using the stochastic method, α(t) = α R (t), only the imaginary part is left, which results in the following equation [34,38]: where Lq s ρ =h −1 [q s , ρ]. This equation was obtained previously in [34,38], which can be regarded as an analogue to the classical generalized Langevin equation, where (t) resembles the fluctuating force, and the hierarchical equation plays a role similar to the dissipation term.
To obtain the stochastic Zusman equation, however, we use another choice of α(t) as the difference between α R (t) and its classical limit of the correlation function such that the ηk B T e −ω c t term is included in the HEOM treatment. Equation (26) utilizes the additive character of the correlation functions shown in equations (22)- (24), and results in the following equation: Using the equivalence between equations (14) and (20) found previously [26], the corresponding continuous form of stochastic Zusman equation is given by Equation (28) is derived for the Debye spectral density, we note that a similar equation can be derived for the Brownian spectral density For simplicity, the associated HEOM form of the stochastic quantum Fokker-Planck equation is not presented here.
Equations (27)- (29) are the main results of this paper. It means that the quantum effect of the bath can be completely recovered by adding proper stochastic fluctuations in the Zusman and quantum Fokker-Planck equations, where the bath has been treated classically. Through the above derivations we can see that, this is possible due to the fact that the imaginary part of the bath correlation functions are the same in both the quantum and classical cases, such that the full quantum dynamics can be recovered after we have taken account of the difference between the real part of the quantum and classical bath correlation functions.

Generating the stochastic process
Although there are many different methods to generate a colored Gaussian noise (t) according to its time correlation function (see e.g. [39,[43][44][45]) we propose in this subsection, an efficient method to generate the stochastic process (t) whose correlation function is given by equation (26). We first note that generation of the stochastic process depends only on the spectral density of the bath. Using equations (8) and (26), we can calculate α(t) as Since tanh(x) < x when x > 0, the correlation function α(t) is positive definite. This means that we only need a real stochastic process in the stochastic Zusman equation. Next, we write equation (30) as a discrete form using equation (4), where The stochastic process (t) can then be generated by an initial sampling of the positions and momenta of a set of harmonic oscillators using and (t) is then calculated as It can be shown that The simulation of the stochastic Zusman equation in equations (27) or (28) can then be performed in the following way: (i) the stochastic variables x j , p j are first generated using equation (33), which are used to calculate the stochastic process (t) by equation (34); (ii) the reduced and auxiliary density operators are then propagated using equation (27) to obtain a single trajectory; and (iii) repeat steps (i) and (ii) for many times, and average over the trajectories to get converged results.

Results
In this section, using several numerical examples, we show applications of the above theory developed in sections 2.2 and 2.3. For simplicity, we will focus on applications of the stochastic Zusman equation. In previous studies, it has been shown that the original Zusman equation is invalid for large reorganization energy in electron transfer reactions [24,25,35], and gives inaccurate results in calculating absorption line shapes of molecular aggregates [17]. We will show that the stochastic Zusman equation does correct all the problems of the original Zusman equation. In addition, using an example of the EET in the Fenna-Matthews-Olson (FMO) complex, we show that the stochastic Zusman equation can be applied to simulations of quantum dissipative dynamics in the intermediate coupling regime.

The spin-boson model
The spin-boson model [1,2] is widely used as a microscopic model for electron transfer reactions. The system Hamiltonian isĤ S = σ z + Vσ x , and the system operator that couples to the bath isq s =σ z , whereσ x = |1 2| + |2 1| andσ z = |1 1| − |2 2|. The initial density matrix of the total system ρ T (0) is set to be the factorized form as |1 1| ⊗ e −βĤ B /Z B , with Z B = Tr e −βĤ B . Figure 1 shows the population dynamics till the rate dynamics is reached, calculated using the parameters = 0, βη = 10, βV = 0.1,hβω c = 1.0. The result was obtained by averaging over 40 000 trajectories of the reduced density operator using the stochastic Zusman equation. As found out previously [26], when the time scale of the system dynamics is faster than the lowest Matsubara frequency (βη 1 in this case), the non-Markovian effect in the quantum correlation function becomes important. In fact, the original Zusman equation becomes invalid When the coupling between the two states is large, or the reorganization energy is small, coherent dynamics in the spin-boson model will be observed.

Absorption line shape
Another example in the literature that has shown the deficiency of the high temperature approximation in the Zusman type equation is in calculation of the absorption line shapes of molecular aggregates [17]. Optical spectroscopy has been an important tool in investigating structural dynamics of complex systems, especially in revealing the coherent dynamics in the EET processes [4,5]. Hence, we will also test the new method in calculations of the absorption line shape of model molecular aggregates. A Frenkel exciton model for N coupled two-level systems is used to describe the excitations of the molecular aggregate. The general form of the total Hamiltonian is written asĤ tot =Ĥ e +Ĥ ph +Ĥ el-ph .
The first term in the above equation (36) denotes the electronic Hamiltonian. When calculating the absorption spectral, only the ground and first excitonic states are considered where | j represents the state where only the jth molecule is excited; j is the transition energy of the jth molecule, and V i j describes the electronic coupling between the ith and jth molecules. The second termĤ ph in equation (36) is the Hamiltonian of the phonon bath, and is treated as a set of harmonic oscillators. Electronic excitation of each molecule is assumed to be independently coupled to the phonon bath, such that where x m j and p m j are the position and momentum of the jth harmonic oscillator bath mode with frequency ω m j . The electron-phonon coupling termĤ el-ph is written aŝ where the c jk s describe the coupling of the excitation of the jth molecule to the harmonic bath. As stated in [17], the absorption line shape can be obtained by calculating the dipole-dipole autocorrelation where with the total dipole operatorμ defined asμ = N i=1 µ i (|0 i| + |i 0|). As a result, we need to propagate the system reduced density matrix ρ S (t) = Tr B [e −iĤ tμ ρ g (0) e iĤ t ], with the initial electronic state ρ g (0) set to be |0 0| ⊗ e −βĤ B /Z B . We further assume that the bath correlation functions for all the molecules have the same Debye spectral density. The stochastic Zusman equation in the HEOM-like form is given by where ρ 0,0,...,0 corresponds to the system reduced density matrix ρ S (t), j (t) is the stochastic process associated with the jth bath mode. For simplicity, we investigate the case of a J -type excitonic dimer (the electronic couplings are negative). The site energies are 0 = 0, 1 = 2 = and the electronic coupling V 12 = V 21 = −|V |. The bath parameters used arehω c = |V |, η = |V | and β|V | = 3. The dipole autocorrelation function in equation (41) is first calculated using the stochastic Zusman equation, and the absorption line shape is then calculated using equation (40). The exact HEOM and Zusman equation results are also presented for comparison.
The normalized real parts of the dipole-dipole correlation function μ(t)μ(0) g e i t are shown in figure 3, and the calculated spectra are shown in figure 4. It can be seen that classical bath approximation in the Zusman type equation again results in too much coherence by neglecting the quantum effect of the bath, while the stochastic Zusman equation has corrected this problem.

Electron energy transfer dynamics in the Fenna-Matthews-Olson complex
We now consider the EET process in the FMO complex. Ishizaki and Fleming [46] have simulated the EET dynamics of the FMO complex using a high temperature approximate  HEOM augmented with low-temperature correction terms. In this work, the stochastic Zusman equation method will be applied to perform the calculation of EET dynamics. The total Hamiltonian is similar to that in equation (36), except that only the one excitonic part is and the equation of motion is the same as that in equation (42). The site energies j of bacteriochlorophyll (BChl) pigment j, and the electron coupling V i j between pigments i and j are taken from [46]. The reorganization energy and the phonon relaxation time are η j = 70 cm −1 , (ω j c ) −1 = 50 fs and temperature is 77 K [46]. In the numerical calculation, the first pigment BChla 1 is chosen to be the initially excited. The simulation results are shown in figure 5, it can be seen that, the original high temperature approximation actually gives unphysical population on pigment 6 at long times ( figure 5(c)). And the stochastic Zusman equation approach give the correct result by adding the proper stochastic terms.

Conclusions and discussions
In summary, by adding a real stochastic process to account for the difference between the quantum and classical bath correlation functions, we have obtained a new method to extend the quantum Fokker-Planck and the Zusman equations beyond the high temperature approximation. Numerical simulations have shown that the stochastic Zusman equation solves the various problems associated with the classical bath approximation, and can be applied in areas of ET, EET and related spectroscopic studies. As an alternative to the exact HEOM based on the Matsubara and other expansion schemes, applications of the stochastic Zusman and quantum Fokker-Planck equation can be further explored.
The idea of using stochastic approaches to unravel the influence functionals in the Caldeira-Leggett model has been investigated intensively [34,38,39,43,44]. It turns out that the real part of the influence functional can be easily unraveled with a real stochastic processes, while the imaginary part of the bath correlation function becomes harder. Stockburger and co-workers have proposed to use a complex stochastic process to solve the problem, and have derived a stochastic Liouville equation and its more efficient trace-conserving form [43,44]. There are also stochastic sampling methods without using the Hubbart-Stratonovich transformation, which has the advantage to take into account the initial system-bath correlation beyond the product initial state [47].
However, numerical difficulties may still arise in cases of strong system-bath couplings due to the use of a complex stochastic process. Later, Shao and co-workers derived the HEOM using the stochastic method based on stochastic decoupling [37], and it is realized that the imaginary part of the bath correlation function can be treated by a deterministic hierarchical method without resorting to the complex stochastic process, while the real part can be treated stochastically [34,38]. A scheme more closely related to this study is the mixed randomdeterministic method in [39], where a portion of the real part of the bath correlation function can also be included in the HEOM treatment. However, the connections to existing theories of the quantum Fokker-Planck and Zusman equations, and the simple extension by adding proper stochastic process to make them exact in the quantum regime are obtained for the first time in this study.
In this work, derivation of the stochastic Zusman and quantum Fokker-Planck equations has been based on the fact that the imaginary part of the correlation function α I (t) is the same for both the quantum and classical bath, such that the additional stochastic process only needs to account for the difference between the real parts of the quantum and classical bath correlation functions. We further note that this term also does not depend on the temperature, so if one can solve the zero (or low temperature) dynamics (although the low temperature calculations are harder for approaches such as the HEOM method, they may be easier for some wave function based methods, see e.g. [48,49]), similar stochastic dynamics scheme may also be used to extend these methods to higher temperatures, by adding a stochastic process (t) to account for the difference between the high temperature and low temperature bath correlation functions. Similar to the discussions in section 2.3, it can be shown that the correlation function of (t) is also positively definite, and a real stochastic process can be used.
As an illustrative example, we consider a simple system of a one mode spin-boson problem We seek to reproduce the high temperature dynamics from the zero temperature dynamics. As discussed above, a stochastic process (t) is added into the Schrödinger equation to take into account the difference between high temperature and zero temperature bath correlation functions where φ T is the total wave functions. The Gaussian stochastic process (t) is given by the correlation function and is generated in a similar approach in section 2.3. The initial state is set to be φ T (0) = |1 ⊗ |0 B .
The parameters of the one-mode spin-boson model are βhω 0 = 1.0, V = 0.5hω 0 , c = 0.2hω 0 . The calculated population dynamics is shown in figure 6. The zero temperature population dynamics, as well as the high temperature dynamic calculated using the deterministic Schrödinger equation (i.e. by choosing the initial states of the bath based on the Boltzmann distribution e −βĤ B ) are also presented for comparison. It can be seen that including the simple stochastic process described in equation (46) is indeed capable of accounting for the difference between high temperature and zero temperature dynamics.