# Reduced density matrix hybrid approach: An efficient and accurate method for adiabatic and non-adiabatic quantum dynamics

###### Abstract

We present a new approach to calculate real-time quantum dynamics in complex systems. The formalism is based on the partitioning of a system’s environment into “core” and “reservoir” modes, with the former to be treated quantum mechanically and the latter classically. The presented method only requires the calculation of the system’s reduced density matrix averaged over the quantum core degrees of freedom which is then coupled to a classically evolved reservoir to treat the remaining modes. We demonstrate our approach by applying it to the spin-boson problem using the noninteracting blip approximation to treat the system and core, and Ehrenfest dynamics to treat the reservoir. The resulting hybrid methodology is accurate for both fast and slow baths, since it naturally reduces to its composite methods in their respective regimes of validity. In addition, our combined method is shown to yield good results in intermediate regimes where neither approximation alone is accurate and to perform equally well for both strong and weak system-bath coupling. Our approach therefore provides an accurate and efficient methodology for calculating quantum dynamics in complex systems.

## I Introduction

The calculation of real-time quantum dynamics for large molecular systems is a longstanding problem in chemical physics. Of particular interest is the case of a small subsystem embedded in a surrounding thermal bath, which forms the basis for the investigation of condensed phase energy and electron transfer as well as spin and charge transport in nanoscale devicesWeiss (2008); Nitzan (2006). In such situations, one is typically interested in the calculation of the system’s reduced dynamics, averaged over the bath degrees of freedom. Due to the importance of these problems and the absence of a general solution, a variety of methods have been developed which vary as to the regimes in which they are accurate and their applicability to large systems. Many of these approaches rely on either averaging over semi-classical trajectories of bath degrees of freedom or propagating the system’s reduced density matrix (RDM) directly without explicit treatment of the bath, for example by master equations or path integral techniques.

Trajectory based methods include the mean-field EhrenfestTully (1998); Stock (1995), surface-hoppingPechukas (1969); Tully (1998, 1990), semi-classical initial value representationMiller (2001, 2009), linearizationShi and Geva (2003); Poulsen, Nyman, and Rossky (2003); Bonella and Coker (2005); Dunkel, Bonella, and Coker (2008), and classical mappingMeyer and Miller (1979); Stock and Thoss (1997); Thoss and Stock (1999) approaches. These offer the appeal of a transparent link to the classical limit and simple application to complex condensed phase systems and to slow baths. However, these methods typically fail in regimes with strong system-bath coupling or high-frequency baths where quantum effects are known to be important.

Techniques based on direct propagation of the reduced density matrix include the Markovian and non-Markovian Bloch-Redfield equationsBloch (1957); Redfield (1965); Breuer and Petruccione (2002) and the noninteracting blip approximationLeggett et al. (1987); Weiss (2008) which serves as a time-dependent generalization of Förster theoryFörster (1953). To make such approaches manageable, perturbation theory in the system-bath coupling (Bloch-Redfield) and in the non-adiabatic coupling (NIBA, Förster) are typically carried out only to second-order and hence performance degrades when these terms become large. In addition, path integral approaches, such as quantum Monte Carlo Mak and Chandler (1989, 1991); Egger and Mak (1994) and the quasi-adiabatic propagator path integral (QUAPI)Makri (1992) have also been developed to propagate the RDM. While these approaches can be made numerically exact, they are frequently difficult to converge in practice when the system-bath coupling is strong or the bath is slow. Furthermore, many variants scale poorly with the number of system states, drastically limiting the size of treatable systems.

Hence when the bath dynamics are “fast” compared to those of the system (the so-called non-adiabatic regime), both perturbative quantum master equations and numerically exact path integral approaches are typically accurate since the correlation time of the bath is short, rendering non-Markovian memory effects less significant. When the bath dynamics are “slow” compared to those of the system (the adiabatic regime), these methods fail but it becomes reasonable to treat the bath degrees of freedom classically and hence trajectory-based schemes generally yield accurate results. Given this complementarity, it would clearly be desirable to develop a scheme which can accurately treat both regimes by naturally tuning between the above approaches in a consistent manner. In addition, such an approach might effectively treat intermediate regimes where no obvious time-scale separation exists and hence neither RDM nor trajectory-based approaches are accurate.

The “dynamical hybrid” scheme of Wang, Thoss and Miller Wang, Thoss, and Miller (2001); Thoss, Wang, and Miller (2001) offers such an approach in which bath degrees of freedom are partitioned into “core” modes to be treated with quantum mechanics and “reservoir” modes to be treated with classical mechanics. In this scheme, the wavefunction-based multi-configurational time-dependent Hartree (MCTDH)Meyer, Manthe, and Cederbaum (1990); Manthe, Meyer, and Cederbaum (1992); Beck et al. (2000) method is used for the system and core while a mean-field Ehrenfest treatment is applied to the reservoir. When applied to the spin-boson problem, this dynamical hybrid scheme frequently converges to the exact result with only a small percentage of modes included in the core. However, to achieve convergence in the presence of fast baths, up to 80% of the bath modes need to be included in the core. Due to the high cost of treating modes using MCTDH, this severely reduces the efficiency of the dynamical hybrid approach and hence limits the ease with which it can be applied to larger systems.

In this work, we present a methodology which avoids the use of wavefunctions, replacing them with a single reduced density matrix averaged over the quantum core modes. The reservoir modes exert a classical time-dependent driving force on this system RDM, which in turn drives the classical modes via a back-reaction force. This reduced density matrix hybrid formalism (RDM-Hybrid) allows one to use an appropriately chosen partitioning of the bath into core and reservoir modes such that the system RDM can be propagated either by an exact reduced dynamics algorithm with a decreased cost or by an inexpensive approximate master equation with improved accuracy. In other words, by using physical intuition, one can choose the core-reservoir partitioning such that different methods are applied only to the set of bath modes for which they are expected to work.

We demonstrate our approach by applying it to the spin-boson model and treating the system and core using the noninteracting blip approximation (NIBA)Leggett et al. (1987); Weiss (2008) This approach yields a methodology which is numerically cheap while obtaining impressive quantitative accuracy over all regimes investigated as well as offering excellent scaling with the number of system states.

The outline of the paper is as follows. Section II.1 briefly reviews the Ehrenfest methodology. This background is then used in Sec. II.2, to derive the RDM-Hybrid approach. Section III then describes how this scheme can be applied to the spin-boson model. The results of the application are presented in Sec. IV, and Sec. V concludes.

## Ii Theory

### ii.1 Ehrenfest dynamics

While there are many routes to the mean-field Ehrenfest equations of motion and subsequently the RDM-Hybrid method to be presented below, we proceed via the quantum-classical Liouville equation, which provides a particularly simple derivation.

We begin by considering the generic coupled system-bath Hamiltonian,

(1) |

where the free bath Hamiltonian is

(2) |

and the coupling is assumed to be of the form

(3) |

The dynamics of the total density matrix, , is given by the Liouville equation,

(4) |

where henceforth we set . By employing the partial Wigner transform of the density matrix with respect to the environmental degrees of freedom,

(5) |

and taking the classical limit of the corresponding partial Wigner transformed Liouville equation (for more details, see e.g. Ref. Kapral, 2006), one obtains the so-called quantum-classical Liouville equation

(6) |

where denotes a commutator and is the Poisson bracket.

Now, one makes the approximation that the system and bath reduced density matrices (RDMs) decouple at all times,

(7) |

with the system RDM and bath RDM . Inserting this product form into the quantum-classical Liouville equation and noting that the bath reduced density matrix is simply the classical phase space distribution,

(8) |

one obtains the mean-field Ehrenfest equations of motionGrunwald, Kelly, and Kapral (2009) for the system RDM,

(9) |

and the classical bath degrees of freedom,

(10) | ||||

(11) |

The Ehrenfest approximation to an observable of the reduced system is given by

(12) |

In accordance with the classical treatment of the bath, this latter integral is evaluated by molecular dynamics obeying the equations of motion given above. To generate the canonical ensemble average consistent with the initial density matrix

(13) |

one averages over classical trajectories with bath initial conditions sampled from either the Wigner-transformed bath density operator,

(14) |

or more simply from the classical Boltzmann distribution,

(15) |

Although the Wigner distribution gives initial conditions consistent with the exact quantum distribution, this property is not conserved by the classical dynamics.

### ii.2 Reduced density matrix hybrid dynamics

In contrast to the approach taken above, in which the full Hamiltonian is split into a system and bath with coupling between the two, we instead divide the bath degrees of freedom into two sets: the reservoir modes () to be treated classically and the core modes () to be treated quantum mechanically. Hence, we partition the Hamiltonian as

(16) |

where the system-core Hamiltonian is

(17) |

and the system-core and system-reservoir couplings are defined as

(18) | ||||

(19) |

Since nothing in our above discussion of Ehrenfest dynamics utilized the specific properties of the system, we could just as well re-label the combined system and core as the system of the previous section. Likewise, we associate the reservoir with the bath from before. Thus we assume the total (Wigner-transformed) density matrix factorizes at all times into a system-core RDM and reservoir RDM,

(20) |

Taking the classical limit we arrive at the Ehrenfest-like RDM-Hybrid equations of motion for the system-core density matrix,

(21) |

and the classical bath degrees of freedom,

(22) | ||||

(23) |

Although the exact calculation of the system-core density matrix dynamics given by Eq. (21) is extremely demanding, this was essentially the approach taken by Wang, Thoss, and Miller (WTM)Wang, Thoss, and Miller (2001); Thoss, Wang, and Miller (2001) who averaged high-dimensional wavefunction trajectories with initial conditions of the core wavefunction sampled from the exact quantum mechanical core density matrix and classical reservoir degrees of freedom sampled from the quantum-classical Wigner distribution.

The important point we seek to make in this work is that is a pure system operator, i.e. the core and reservoir are not directly coupled, such that the classical reservoir equation of motion above may be written

(24) |

where we have introduced the system RDM, , averaged over the quantum core degrees of freedom. Furthermore, the calculation of any dynamical system variable is given by an expression analogous to the pure Ehrenfest approximation above,

(25) |

which again only requires the system RDM (and not the much higher-dimensional system-core RDM), with classical reservoir trajectories to be sampled from an appropriate distribution.

We have thus arrived at a self-consistent dynamical scheme which treats the core bath modes quantum mechanically yet only requires the propagation of a typical system RDM and classical reservoir degrees of freedom. This RDM-Hybrid formulation presented above constitutes the major theoretical result of this work. We point out that a similar approach appears to have been investigated by Golosov et al.Golosov, Friesner, and Pechukas (2000), though it was restricted to their previously developed memory-equation algorithmGolosov, Friesner, and Pechukas (1999) and only minimally pursued.

When the system RDM averaged over the quantum core modes, , is treated fully quantum mechanically, our methodology yields the exact result when all bath modes are included in the core. However, in practice one expects to achieve numerical convergence before this limit is reached. Of course, the exact quantum dynamical treatment of the system RDM is in general no less numerically challenging than the original problem (a system coupled to a quantum bath). However, the advantage of this approach is that one can tailor the partitioning in Eq. (16) either to alleviate the computational expenses of a numerically exact reduced quantum dynamics method (such as the QUAPI-based iterative tensor propagation scheme of Makri et al.Makri (1992); Makarov and Makri (1994); Makri and Makarov (1995a, b)) or to improve the accuracy of approximated quantum dynamics (such as a perturbative quantum master equation) for the system RDM. In this sense, the above presented formalism outlines a framework for the efficient partitioning and calculation of real-time reduced quantum dynamics where one can employ any exact or approximate dynamical scheme for the system RDM depending on the complexity of the problem. As discussed in the following section, here we adopt an approximate treatment of the system RDM using a perturbative master equation to yield a very cheap hybrid methodology with impressive quantitative accuracy.

## Iii Application to the spin-boson problem

While the above formalism is indeed quite general here we demonstrate its effectiveness at treating the spin-boson Hamiltonian,

(26) |

which describes a two-level system with energy bias and tunneling matrix element , linearly coupled to the coordinates of a bath of harmonic oscillators. The spin-boson Hamiltonian has been used as a model for a variety of physically distinct quantum relaxation and transport processes.

A principal observable of interest is the population variable,

(27) |

which measures the difference in population of sites 1 and 2. Furthermore, we will henceforth work in reduced units with respect to the tunneling matrix element, , in addition setting . To begin our application, we first discuss the separation of modes into the core and the reservoir.

### iii.1 Separation of modes

The harmonic bath in the spin boson problem is fully specified by its spectral density

(28) |

When separating the bath modes into a core and reservoir, we will require that

(29) |

Intending to treat high frequency modes in the quantum mechanical core and low frequency modes in the classical reservoir, we use a switching function, , which switches from 1 to 0 as goes from 0 to . This allows the core and reservoir spectral densities to be defined as,

(30) | ||||

(31) |

An infinite number of switching functions satisfy such a criteria and the best form will depend on the Hamiltonian adopted, the methodology used to treat the core, and the required balance been accuracy and efficiency. The simplest approach would be to use a step function,

(32) |

However in many cases, such a sharp switch will leave a core which is expensive to treat using quantum mechanics. The reason for this difficulty can be seen by considering the bath correlation function, given by

(33) |

whose correlation time determines the range of memory (non-Markovian) effects. It is straightforward to see that using a step function switching like that in Eq. (32) will result in a long-time oscillatory tail in the bath correlation function, , and thus extensive non-Markovian effects which are difficult to treat quantum mechanically. It is therefore advantageous to separate the spectral density using a smooth switching function which can yield a short bath correlation time for the core modes and allow for efficient quantum mechanical treatment.

For our purposes a switching function of the form,

(34) |

was chosen. In Fig. 1, we show the bath correlation function when this function is applied to the Debye spectral density,

(35) |

where is the reorganization energy. As can be seen, this smooth switching does not introduce any long-time tails and allows the bath correlation time for the core modes to be decreased as is increased and more of the slow modes, which give rise to temporally non-local effects, are moved into the reservoir. In addition, removal of slow modes from the core encourages the use of approximate master equations, such as the noninteracting blip approximation discussed in the next section, which complement the Ehrenfest treatment of the slow modes.

### iii.2 Noninteracting blip approximation for the system and core

Although the RDM-Hybrid framework allows for the application of many approaches to treat the system and core modes, in this work we employ the noninteracting blip approximation (NIBA)Leggett et al. (1987); Weiss (2008). To see why such a choice is computationally simple, note that because the fluctuating classical degrees of freedom provide an effective time-dependent bias, the system Hamiltonian no longer commutes with itself at different times, and thus its propagator is a cumbersome time-ordered exponential (Dyson series). However, since NIBA only treats the diagonal part of this Hamiltonian exactly, its propagator is trivially calculated.

Another important reason for our choice is the complementarity between NIBA and Ehrenfest. While Ehrenfest dynamics work best in the adiabatic regime, , NIBA is perturbative in , thus working best in the non-adiabatic regime, . To exemplify this complementarity, in Fig. 2 we show a coherent to incoherent crossover, as a function of in the high-temperature, strong-coupling regime. As can be seen in Fig. 2, the Ehrenfest method yields quantitatively exact population dynamics for an adiabatic choice of parameters, [panels (a) and (b)], whereas NIBA gives the exact result upon crossing over into the non-adiabatic regime, [panels (c) and (d)]. With this in mind we are strongly encouraged to consider a NIBA-Ehrenfest hybridization, the details of which are described next.

### iii.3 Simulation details

The total bath was discretized by the relation,

(36) |

where is a density of frequencies chosen to reproduce the reorganization energy for any number of modes, . For an Ohmic spectral density with cutoff function , we thus require

(37) |

which is clearly satisfied by the choice

(38) |

For the Debye spectral density considered here, the cutoff function is given by .

Having obtained a discrete distribution of frequencies, we then take into account the switching function and recalculate the reservoir couplings as

(39) |

in the process removing all bath modes with from the reservoir. In a complimentary fashion, the core spectral density, , is used in the NIBA equations.

In all results presented, discrete bath modes were found to be sufficient for convergence and averages were performed over reservoir initial conditions. For pure Ehrenfest results, we sample from the Wigner distribution, which for the spin-boson Hamiltonian’s harmonic bath takes the form,

(40) |

As alluded to above, when high frequency baths are present, Wigner sampling can improve the accuracy of Ehrenfest dynamics since the initial conditions are consistent with the exact quantum mechanical distribution. However, since the zero-point energy inserted by Wigner sampling is not conserved by the classical dynamics, long-time predictions can be unreliable (zero-point energy leakage)Stock and Müller (1999); Müller and Stock (1999). Hence, in regimes with strong coupling or high bath frequencies, sampling from either the classical Boltzmann or quantum Wigner distributions can produce inaccurate results. In contrast, for our RDM-Hybrid approach we find that the results are largely insensitive to the choice of reservoir initial conditions, since most quantum mechanical modes are included in the core leading to only low frequencies being present in the reservoir, for which classical and Wigner sampling give identical results. This effect is shown later, in Fig. 4, where pure Ehrenfest dynamics using classical sampling is seen to underestimate the decay of the system population variable while using Wigner sampling overestimates the decay. The RDM-Hybrid approach, however, produces graphically identical results using either sampling scheme.

The time evolution of the coupled classical dynamics and NIBA equations were solved with a second-order Runge-Kutta scheme with a timestep of . For the spin-boson Hamiltonian considered here, the reservoir equations of motion, Eqs. (22)-(24), take the form

(41) | ||||

(42) |

with the system population variable . Taking to be constant over a half-timestep as called for by the Runge-Kutta scheme used here, one has

(43) | ||||

(44) |

For a two-level system coupled to a common bath, the NIBA master equation for the population difference with a time-dependent bias is given simply asGoychuk, Petrov, and May (1995)

(45) |

with the kernels

(46) | ||||

(47) |

The parameter originates from the bath initial condition, with if the system and bath are initially uncoupled and if the bath is initially solvating the donor. The twice-integrated bath correlation function, , is given by

(48) |

The effect of a general time-dependent bias is captured in the function

(49) |

which for the core-reservoir splitting considered here is defined as

(50) |

The true population difference, is finally given by the average of over classical trajectories with initial conditions sampled from the appropriate distribution.

Importantly, the scheme described here yields the correct population dynamics of a two-level system driven by a classical external field, so that by taking the limit , we remove all bath modes from the core, and recover the Ehrenfest solution. In the opposite limit, , we recover the full NIBA treatment.

Because we are not treating the core exactly, it is important to emphasize that is no longer a convergence parameter. Instead, tunes the balance between the two approximate methods, here taken to be Ehrenfest dynamics and NIBA. The optimal a priori choice of this switching frequency is an interesting problem for future work, but here we present a physically motivated prescription for its determination. We again recall that “adiabaticity” is inherently a problem of timescale separation. The adiabatic regime is physically realized when the timescale of the bath is much greater than that of the system, and vice versa for the non-adiabatic regime. A given bath mode could be classified into “core” or “reservoir” by comparing its characteristic timescale, with that of the system, . For the spin-boson Hamiltonian, Eq. (26), the system frequency is given by the Rabi frequency,

(51) |

In practice we use the smooth switching function, Eq. (34), and following the line of argument above, set the switching frequency equal to the Rabi frequency, i.e. . This procedure effectively partitions the bath modes into those that are faster than the system dynamics (comprising the core) and those that are slower (comprising the reservoir). While comparison with existing exact results has shown that this choice, , is not always truly optimal, we will show in the next section that it nevertheless provides a very robust methodology, with quantitative predictive ability.

## Iv Results

In all our results we compare to the numerically exact population dynamics computed by Thoss, Wang, and MillerThoss, Wang, and Miller (2001) for the spin-boson Hamiltonian, Eq. (26), with a Debye spectral density, Eq. (35) and initial condition , i.e. in Eqs. (46)-(47).

We begin in the adiabatic regime, , where Ehrenfest dynamics are known to be relatively accurate. Fig. 3(a) shows that this is indeed the case for high temperature, . Though qualitatively good, the very strong coupling, , degrades the accuracy of the Ehrenfest approach. However, our RDM-Hybrid method yields a long-time population decay in much better agreement with exact results. In Fig. 3(b), we investigate another system with strong coupling, , but with a temperature one order of magnitude lower. The classical mechanical approximation made by using Ehrenfest dynamics intrinsically worsens for lower temperature, where quantum mechanical effects are expected to play a more dominant role. We see that the RDM-Hybrid result is again in better quantitative agreement with the exact long-time limit. Both panels (a) and (b) in Fig. 3 underscore the important point that our RDM-Hybrid approach is not simply an “average” of the two methods, but can in fact generate qualitatively different dynamics which are not situated “in between” those of the individual methods.

The system trapping seen above leading to a quasi-stationary state is nearly impossible to capture with perturbative quantum master equation approaches, as exemplified by the poor NIBA results. For comparison, the WTM MCTDH-Ehrenfest scheme required 25% of the bath-modes to be treated quantum mechanically, supporting our observation of largely classical bath dynamics.

In Fig. 4, we investigate the non-adiabatic regime, , with strong coupling , where NIBA is quantitatively accurate. Here we show the effects of sampling from a classical distribution or a Wigner distribution. The Ehrenfest approximation with Wigner initial conditions yields a decay which is too rapid, while that with classical initial conditions yields a decay which is too slow. Both completely miss the two-step relaxation dynamics at very short times. RDM-Hybrid captures these short time dynamics exactly and shows a robust insensitivity to the initial sampling employed. Although our hybrid method exhibits a slight discrepancy at longer times, it is much more accurate than the Ehrenfest approach.

Having demonstrated that our RDM-Hybrid approach can yield very good agreement with the exact results when only one of its composite methods (Ehrenfest dynamics or NIBA) is successful, we now consider in Fig. 5 the intermediate regime, , where neither method is particularly accurate by itself. Fig. 5(a) depicts the population dynamics at the high temperature for relatively strong coupling, . Here, both NIBA and Ehrenfest dynamics are in very poor numerical agreement with the exact result, though the Ehrenfest dynamics correctly predicts a two-step relaxation, albeit only qualitatively. On the contrary, RDM-Hybrid correctly exhibits a rebound in the population dynamics near (although one that is somewhat overpronounced) and yields dynamics at all later times in excellent agreement with the exact result.

In the lower panel, Fig. 5(b), we consider a weaker coupling and a temperature two orders of magnitude lower. Here, NIBA does quite poorly and Ehrenfest dynamics are mildly better, exhibiting slightly overdamped oscillations and a weak phase-shift at long times. Once again, RDM-Hybrid is nearly exact, despite this severely low temperature.

In both examples above, i.e. Figs. 5(a) and (b), Thoss, Wang, and Miller found it necessary to treat over 50% of the bath modes quantum mechanically with MCTDH. It is remarkable that our computationally inexpensive RDM-Hybrid method employing NIBA for the core modes is able to yield such similarly accurate results for these problems, which serves as a testament to its robustness.

We conclude this section by briefly considering the effect of an energetic bias on the performance of our RDM-Hybrid method. In Figs. 6(a) and (b), we show two sets of population dynamics for the same high-temperature, adiabatic regime, but with an energetic bias in panel (b). As expected, NIBA performs poorly due to the adiabaticity of the conditions under consideration in both panels. As for the Ehrenfest method, while it gives excellent results in the unbiased case [panel (a)], it strongly deviates from the correct dynamics in the presence of a bias [panel (b)]. This deficiency of the Ehrenfest approach and similar quantum-classical methods for biased system is well-known and typically connected to unphysical treatment of zero-point energy in classical treatments. Attempts to correct this behavior have included an adjustable zero-point energy parameter for sampling the initial conditions Müller and Stock (1999) and its determination by comparison with exact short-time moments.Golosov and Reichman (2000)

In spite of this potential concern, we see in Fig. 6 that the RDM-Hybrid approach performs incredibly well and equally so for both unbiased and biased systems. Unfortunately, this success is not universal, as seen in Fig. 7, which considers the effect of a bias on non-adiabatic dynamics. While RDM-Hybrid is very successful for the unbiased case, it performs more poorly in the presence of a bias. This error appears to be attributable to the Ehrenfest dynamics, which RDM-Hybrid tracks quite closely at short times. Nonetheless, the RDM-Hybrid population does appear to approach the correct long-time value, which both Ehrenfest and NIBA incorrectly predict.

## V Conclusions and Future Work

To summarize, we have formulated a quantum-classical methodology for the quantum dynamics of a coupled system-bath Hamiltonian, averaged over the bath degrees of freedom, by employing a core-reservoir partitioning of the bath degrees of freedom. Unlike the technique of Wang, Thoss, and MillerWang, Thoss, and Miller (2001), our RDM-Hybrid approach avoids the use of expensive wavefunction-based quantum mechanics, replacing it by a reduced density matrix calculation. This allows for flexibility to use a variety of inexpensive approximate approaches for different parts of the computation, allowing for combined accuracy and efficiency.

Within the above framework, we then used the noninteracting blip approximation (NIBA) for the reduced dynamics, yielding an efficient, scalable quantum dynamics methodology with excellent accuracy for both adiabatic and non-adiabatic parameter regimes. Specifically, for discrete states, a numerical implementation of NIBA scales only as (see Eq. (45)), which allows for very efficient investigation of bigger systems. The largest discrepancy with exact results was found in the non-adiabatic regime with an energetic bias. To alleviate such deficiencies, we enumerate a number of avenues for future work.

One could imagine pursuing alternative approximate treatments of the quantum core-averaged system RDM. For example, the accuracy of NIBA (employed here) is known to degrade at low temperature, in particular for biased systems. The Redfield equations, on the other hand, have been very successfully applied to such systems, being restricted, however, to weak-coupling and non-adiabatic baths. Thus, an RDM-Hybrid approach employing Redfield and Ehrenfest dynamics may be one way to extend the validity of the Redfield equations into these regimes where it would otherwise fail.

Ideally, one would like an exact treatment of the quantum system RDM, such that the method becomes numerically exact when all modes are included in the core. As alluded to at the end of Sec. II.2, we believe the tensor propagation scheme of Makri et al. Makri (1992); Makarov and Makri (1994); Makri and Makarov (1995a, b) is especially promising. Because the QUAPI algorithm scales exponentially with the memory time required to span the bath correlation function, the partitioning discussed here leading to a reduced bath correlation time would constitute an enormous reduction in computational expense. Work along this direction is currently in progress.

Lastly, for a chosen system-core dynamics method, a more rigorous investigation of the form of the switching function, , and switching frequency, , should be pursued. Ideally these choices should be automated and not left up to physical considerations as was done here. For example, one might try comparing the short-time moments of the population dynamics to those obtained exactly from quantum perturbation theory. The investigation of this point will be a subject of future work.

In spite of these desired improvements, the RDM-Hybrid method employing NIBA has been shown to be very accurate, often correcting the long-time populations of the otherwise accurate Ehrenfest dynamics. Future work will include the application and investigation of both of these methods to models of electronic energy transfer in molecular aggregates, including the multi-site Frenkel exciton Hamiltonian describing the Fenna-Matthews-Olsen complex, a prototypical photosynthetic system.

###### Acknowledgements.

The authors gratefully thank Aaron Kelly for a critical reading of the manuscript and helpful comments. T.C.B. was supported by the Department of Energy Office of Science Graduate Fellowship Program (DOE SCGF), administered by ORISE-ORAU under Contract No. DE-AC05-06OR23100. D.R.R. was supported by the National Science Foundation under Grant No. CHE-0719089.## References

- Weiss (2008) U. Weiss, Quantum Dissipative Systems (World Scientific Publishing, 2008).
- Nitzan (2006) A. Nitzan, Chemical Dynamics in Condensed Phases (Oxford University Press, 2006).
- Tully (1998) J. C. Tully, Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B. J. Berne, G. Ciccotti, and D. F. Coker (World Scientific Publishing Company, 1998).
- Stock (1995) G. Stock, J. Chem. Phys. 103, 1561 (1995).
- Pechukas (1969) P. Pechukas, Phys. Rev. 181, 174 (1969).
- Tully (1990) J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
- Miller (2001) W. H. Miller, J. Phys. Chem. A 105, 2942 (2001).
- Miller (2009) W. H. Miller, J. Phys. Chem. A 113, 1405 (2009).
- Shi and Geva (2003) Q. Shi and E. Geva, J. Chem. Phys. 118, 8173 (2003).
- Poulsen, Nyman, and Rossky (2003) J. A. Poulsen, G. Nyman, and P. J. Rossky, J. Chem. Phys. 119, 12179 (2003).
- Bonella and Coker (2005) S. Bonella and D. F. Coker, J. Chem. Phys. 122, 194102 (2005).
- Dunkel, Bonella, and Coker (2008) E. R. Dunkel, S. Bonella, and D. F. Coker, J. Chem. Phys. 129, 114106 (2008).
- Meyer and Miller (1979) H.-D. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979).
- Stock and Thoss (1997) G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997).
- Thoss and Stock (1999) M. Thoss and G. Stock, Phys. Rev. A 59, 64 (1999).
- Bloch (1957) F. Bloch, Phys. Rev. 105, 1206 (1957).
- Redfield (1965) A. G. Redfield, Adv. Magn. Reson. 1, 1 (1965).
- Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, 2002).
- Leggett et al. (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
- Förster (1953) T. Förster, Discuss. Faraday Soc. 27, 7 (1953).
- Mak and Chandler (1989) C. H. Mak and D. Chandler, Phys. Rev. A 41, 5709 (1989).
- Mak and Chandler (1991) C. H. Mak and D. Chandler, Phys. Rev. A 44, 2352 (1991).
- Egger and Mak (1994) R. Egger and C. H. Mak, Phys. Rev. B 50, 15210 (1994).
- Makri (1992) N. Makri, Chem. Phys. Lett. 193, 435 (1992).
- Wang, Thoss, and Miller (2001) H. Wang, M. Thoss, and W. H. Miller, J. Chem. Phys. 115, 2979 (2001).
- Thoss, Wang, and Miller (2001) M. Thoss, H. Wang, and W. H. Miller, J. Chem. Phys. 115, 2991 (2001).
- Meyer, Manthe, and Cederbaum (1990) H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990).
- Manthe, Meyer, and Cederbaum (1992) U. Manthe, H.-D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 97, 3199 (1992).
- Beck et al. (2000) M. H. Beck, A. Jackle, G. A. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000).
- Kapral (2006) R. Kapral, Annu. Rev. Phys. Chem. 57, 129 (2006).
- Grunwald, Kelly, and Kapral (2009) R. Grunwald, A. Kelly, and R. Kapral, Energy Transfer Dynamics in Biomaterial Systems, edited by I. Burghardt, V. May, D. A. Micha, and E. R. Bittner (Springer-Verlag, 2009).
- Golosov, Friesner, and Pechukas (2000) A. A. Golosov, R. A. Friesner, and P. Pechukas, J. Chem. Phys. 112, 2095 (2000).
- Golosov, Friesner, and Pechukas (1999) A. A. Golosov, R. A. Friesner, and P. Pechukas, J. Chem. Phys. 110, 138 (1999).
- Makarov and Makri (1994) D. E. Makarov and N. Makri, Chem. Phys. Lett. 221, 482 (1994).
- Makri and Makarov (1995a) N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4600 (1995a).
- Makri and Makarov (1995b) N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4611 (1995b).
- Stock and Müller (1999) G. Stock and U. Müller, J. Chem. Phys. 111, 65 (1999).
- Müller and Stock (1999) U. Müller and G. Stock, J. Chem. Phys. 111, 77 (1999).
- Goychuk, Petrov, and May (1995) I. A. Goychuk, E. G. Petrov, and V. May, Phys. Rev. E 52, 2392 (1995).
- Golosov and Reichman (2000) A. A. Golosov and D. R. Reichman, J. Chem. Phys. 114, 1065 (2000).