We present a modified approach for simulating electronically nonadiabatic dynamics based on the Nakajima-Zwanzig generalized quantum master equation (GQME). The modified approach utilizes the fact that the Nakajima-Zwanzig formalism does not require casting the overall Hamiltonian in system-bath form, which is arguably neither natural nor convenient in the case of the Hamiltonian that governs nonadiabatic dynamics. Within the modified approach, the effect of the nuclear degrees of freedom on the time evolution of the electronic reduced density operator is fully captured by a memory kernel super-operator. A methodology for calculating the memory kernel from projection-free inputs is developed. Simulating the electronic dynamics via the modified approach, with a memory kernel obtained using exact or approximate methods, can be more cost effective and/or lead to more accurate results than direct application of those methods. The modified approach is compared to previously proposed GQME-based approaches, and its robustness and accuracy are demonstrated on a benchmark spin-boson model with a memory kernel which is calculated within the Ehrenfest method.

## I. INTRODUCTION

Quantum dynamical effects play a central role in a variety of important processes that take place in molecular condensed phase systems, including vibrational and electronic relaxation, charge transfer, optical response, and photovoltaics.^{1–3} As a result, the simulation of quantum dynamics in such systems remains one of the most important challenges facing computational chemistry. However, the exponential scaling of the computational cost with system dimensionality makes the numerically exact simulation of quantum dynamics in complex molecular systems non-feasible, with the important exception of a subclass of Hamiltonians whose form makes such an exact simulation possible.^{4–10}

Nonadiabatic dynamics corresponds to an important class of inherently quantum-mechanical dynamical processes that plays a central role in a variety of important chemical processes, ranging from redox reactions to photovoltaics.^{11–22} A variety of mixed quantum-classical and semiclassical approximate methods for simulating nonadiabatic dynamics have been proposed over the years, including Ehrenfest,^{23} surface hopping,^{24–35} mixed quantum-classical Liouville (MQCL),^{36–42} and the symmetrical quasiclassical (SQC) method.^{43–48} What those methods have in common is that the dynamics of the nuclear degrees of freedom (DOF) is described in terms of classical-like trajectories. However, the classical-like nature of the dynamics of the nuclear DOF within those approximate methods often causes their accuracy to decrease and/or their computational cost to increase with increasing simulation time.

Another approach for simulating nonadiabatic dynamics in complex condensed phase systems is based on the fact that the dynamics of the electronic populations and coherences can be expressed in terms of quantities like correlation functions and memory kernels.^{1–3,49–69} An advantage of describing the dynamics in terms of such quantities is that they are much simpler than the overall system density operator. Furthermore, for many complex condensed-phase systems of interest, such quantities are often found to be characterized by finite memory times. This implies that one should be able to follow a multi-scale approach, where the finite-lifetime input quantities necessary for simulating the dynamics of experimentally relevant observables on longer time scales are obtained via exact or mixed quantum-classical/semiclassical approximations.

Methods based on Fermi’s golden rule (FGR) and the Redfield equation are examples of the aforementioned approach.^{1,65–78} In those methods, the impact of the nuclear DOF on the electronic DOF is given in terms of typically short-lived two-time correlation functions. However, those methods are based on treating the coupling between electronic states or between the electronic and nuclear DOF as a small perturbation. By contrast, methods based on the generalized quantum master equation (GQME)^{50,52,64,79–84} can describe coupling of any strength, making the memory kernel of the GQME arguably the most general and rigorous form of such a short-time input quantity.

In this paper, we focus on the case where the memory kernel captures the impact of the nuclear DOF on the electronic DOF, for a system whose overall Hamiltonian has the following general form (in what follows, vector quantities are boldfaced, e.g., **R**, and a hat over a symbol, e.g., *Â*, indicates an operator quantity):

Here, $\u0124j=P^2/2+VjR^$ is the nuclear Hamiltonian when the system is in the electronic state | *j*⟩, $R^=R^1,\u2026,R^N$ and $P^=P^1,\u2026,P^N$ are the mass-weighted position and momentum operators of the *N* nuclear DOF, and $V^jk=VjkR^$ couple electronic states to each other.

We restrict ourselves to the case where the electronic states, {| *j*⟩}, are independent of $R^$ (e.g., the so-called crude adiabatic basis^{85}). The index *j* in Eq. (1) runs over the *M* electronic states. For example, a two-state donor-acceptor system would correspond to *M* = 2. Many processes involving nonadiabatic dynamics in condensed phase systems are described in terms of a Hamiltonian of the form of Eq. (1). Furthermore, describing nonadiabatic dynamics in terms of a crude adiabatic representation, rather than in terms of Born-Oppenheimer representation, is not an approximation since the dynamics is independent of the representation employed.^{85}

Importantly, ${R^,P^}$ are meant to correspond to the Cartesian positions and momenta of the individual atoms in a complex molecular system that would typically consist of a large number (>10^{2}) of atoms. Thus, $VjR^$ would typically correspond to the electronic-state-specific potential energy surface (PES) that describes the interaction between the atoms in electronic state | *j*⟩. Similarly, the electronic coupling terms, $VjkR^$, can be $R^$-dependent (assuming that they are $R^$-independent corresponds to the Condon approximation). Our working hypothesis is that $VjR^$ and $VjkR^$ can be obtained from on-the-fly electronic structure calculations and/or semi-empirical force fields (e.g., see Ref. 69).

In most cases of practical importance, the initial state of the overall (nuclear + electronic) system would be given by a density operator of the following form:

Here, $\rho ^n(0)=Tre{\rho ^(0)}$ is the reduced density operator that describes the initial state of the nuclear DOF, where Tr_{e}{·} stands for partially tracing over the electronic Hilbert space. Similarly, the reduced density operator that describes the initial state of the electronic DOF is obtained by partially tracing over the nuclear Hilbert space,

The state of the overall system at a later time *t* would then be given by

Here, *Ĥ* is the overall Hamiltonian, Eq. (1), and $L(\u22c5)=[\u0124,\u22c5]$ is the corresponding Liouvillian (calligraphic font, e.g., $L$, indicates that the quantity is a super-operator). The nuclear and electronic states at time *t* are described by the corresponding reduced density operators,

Importantly, knowledge of $\sigma ^(t)$ would allow for the evaluation of both the electronic populations, {*σ*_{jj}(*t*)}, and coherences, {*σ*_{jk}(*t*)| *j* ≠ *k*}.

It should be pointed out that extending the methodology to non-single-product overall system initial states is rather straightforward. For example, given that $\rho ^(0)=p\rho ^na(0)\u2297\sigma ^a(0)+(1\u2212p)\rho ^nb(0)\u2297\sigma ^b(0)$, with 0 < *p* < 1, the reduced density operator at time *t* would be given by

Thus, each of the two reduced density operators, $\sigma ^a(t)$ and $\sigma ^b(t)$, can be propagated independently of the other starting at its own single-product overall system initial state, and the overall reduced density operator, $\sigma ^(t)$, can be obtained from the linear combination. The same procedure can be applied to other forms of the initial state. To this end, it should be noted that any initial state of the overall system, $\rho ^(0)$, can be expressed as $\u2211j,k\rho ^jk(0)|jk|$, with $\rho ^jk(0)=\u27e8j|\rho ^jk(0)|k\u27e9$. Thus, the computational cost scales linearly with respect to the number of product terms with independent $\rho ^jk(0)$, which is rather favorable.

The *modified GQME-based (M-GQME)* scheme presented herein builds upon the *Nakajima-Zwanzig GQME (NZ-GQME)*,^{3,50,52,53,55–60,63,64} but it is distinctly different from previously proposed NZ-GQME-based schemes^{64,79,86} (see Appendix B), as will be shown below. The M-GQME scheme avoids the commonly employed assumption that the overall Hamiltonian is cast in a system-bath form,

Here, *Ĥ*_{S} is the system Hamiltonian, *Ĥ*_{B} is the bath Hamiltonian, and *Ĥ*_{BS} is the coupling between the system and bath. In the context of the system under consideration here, the system would stand for the electronic DOF, the bath would stand for the nuclear DOF, and the system-bath coupling, *Ĥ*_{BS}, would stand for the coupling between the nuclear and electronic DOF.

While casting the overall Hamiltonian in the system-bath form of Eq. (7) has proven to be extremely useful in many other contexts, it is neither natural nor convenient when dealing with an overall Hamiltonian of the form of Eq. (1). This is because the first term in Eq. (1), *∑*_{j}*Ĥ*_{j}| *j*⟩⟨ *j*|, associates a *different* nuclear Hamiltonian, *Ĥ*_{j}, with each electronic state, | *j*⟩, thereby making it impossible to come up with a uniquely defined bath Hamiltonian, *Ĥ*_{B}. It should be noted that while it is in principle possible to cast the Hamiltonian in Eq. (1) in the form of Eq. (7), the fact that there is no one unique way of accomplishing this can complicate the implementation of a GQME-based approach which is based on the system-bath form in a number of ways. First, when using an approximate method to evaluate the memory kernel, different choices of *Ĥ*_{B}, *Ĥ*_{S}, and *Ĥ*_{BS} may lead to different results without a clear criterion for choosing between them. Second, the nuclear DOF are often assumed to start out at equilibrium with respect to *Ĥ*_{B} such that $\rho ^n(0)=\rho ^Beq=exp[\u2212\beta \u0124B]/ZB$, where *Z*_{B} is the canonical partition function to *Ĥ*_{B}, which means that the definition of *Ĥ*_{B} needs to change whenever the nuclear initial conditions, i.e., $\rho ^n(0)$, do. Third, the system-bath coupling, *Ĥ*_{BS}, is often defined so that $TrB{\u0124BS\rho ^Beq}=0$, which implies that the definition of *Ĥ*_{BS} will also be dependent on the choice of *Ĥ*_{B}. Fourth, the projection operator used to derive the GQME is often defined as $P(\u22c5)=\rho ^Beq\u2297TrB{\u22c5}$ and would also need to be modified when the definition of *Ĥ*_{B} changes. Fifth, the second term in Eq. (1), $\u2211j,k\u2260jV^jk|jk|$, is purely electronic in the Condon approximation, $V^jk\u2192Vjk$, and therefore is part of the system Hamiltonian. However, this term becomes a system-bath coupling term in the non-Condon case, thereby making it difficult to create a unified computational framework that can treat both Condon and non-Condon cases.

In this paper, we pursue a GQME-based approach that is free of the requirement to cast the overall system Hamiltonian in system-bath form. The remainder of this paper is organized as follows. The M-GQME approach for the nonadiabatic dynamics of the electronic DOF, which does not require an overall Hamiltonian of system-bath form, is derived in Sec. II. The methodology for calculating the memory kernel within the M-GQME approach is developed in Sec. III. A protocol for using the Ehrenfest method for calculating the projection-free quantities needed for calculating the memory kernel is outlined in Sec. IV. Demonstrative applications based on using the Ehrenfest method for calculating the memory kernel within the M-GQME approach on benchmark spin-boson systems are presented and analyzed in Sec. V. The main conclusions are outlined and discussed in Sec. VI. Technical aspects are discussed in Appendixes A–D.

## II. THE MODIFIED GENERALIZED QUANTUM MASTER EQUATION (M-GQME)

We start out with a system whose overall Hamiltonian is given by Eq. (1), rather than by Eq. (7). As is well known, the dynamics of the projected state, $P\rho ^(t)$, for any projection operator $P$ that satisfies idempotence ($P2=P$) is given by the NZ-GQME^{1,2,87,88}

Here, $Q=1\u2212P$ is the complimentary projection operator to $P$ (i.e., $Q$ projects-in what $P$ projects-out).

Next, we explicitly define the projection operator as follows:

Here, $\rho ^nref$ is a reference nuclear density operator of one’s choice (as long as $Trn\rho ^nref=1$, which is required for $P2=P$). Substituting the projection operator in Eq. (9) into Eq. (8) and tracing over the nuclear Hilbert space then leads to the following GQME for the reduced electronic density operator $\sigma ^(t)$,

Here,

accounts for the Hamiltonian and Markovian dynamics generated by the Hamiltonian $\u27e8\u0124\u27e9nref=Trn\u0124\rho ^nref$, while the remaining two terms on the R.H.S. account for the non-Hamiltonian and non-Markovian dynamics generated by the coupling between the electronic and nuclear DOF. More specifically, the memory kernel, $K(\tau )$, which accounts for the effect of the nuclear DOF within the time interval (0, *t*) on the electronic DOF at time *t*, is given by

and the inhomogeneous term, *Î*(*t*), which accounts for the effect of the initial state of the nuclear DOF on the electronic DOF at time *t*, is given by

It should be noted that both $K(\tau )$ and *Î*(*t*) would typically vanish at *τ*, *t* > *τ*_{c}, where *τ*_{c} is the characteristic finite correlation or memory time of the electronic DOF.

As is well known, there is no one unique choice of $\rho ^nref$ in Eq. (9) and the specific choice is dictated by the questions one is asking and convenience.^{79,89,90} As a result, different choices of $\rho ^nref$ would lead to different versions of the GQME. In this sense, the equation of motion that governs the dynamics of the electronic DOF is not unique, although the different equations of motion must all reproduce the same electronic dynamics (as long as the quantum-mechanically exact memory kernel and inhomogeneous term can be obtained). In practice, it is convenient to choose $\rho ^nref$ in a manner that will simplify the resulting GQME. The assignment $\rho ^nref=\rho ^n(0)$ (the initial state of the nuclear DOF) is such a convenient choice and leads to the following definition of the projection operator:

This choice is convenient because it leads to the elimination of the inhomogeneous term from the GQME.^{79} However, it should be noted that this choice also implies that the memory kernel will be explicitly dependent on the initial state of the nuclear DOF. In other words, changing the state of the nuclear DOF at the initial time *t* = 0 [$\rho ^n(0)$] would alter the equation of motion [see Eq. (15)] that dictates the dynamics of the electronic DOF at later times *t* > 0. It should also be noted that the specific form of $\rho ^n(0)$ is dictated by how the system is prepared, which is ultimately dependent on the experimental setup. Importantly, $\rho ^n(0)$ need not be of the form $\rho ^Beq=ZB\u22121e\u2212\beta \u0124B$. It should be noted that this is also not required within the Zhang-Ka-Geva (ZKG-NZ) scheme (see Appendix B).^{79}

Substituting the projection operator in Eq. (14) into Eq. (10) leads to the following GQME for the electronic reduced density operator, $\sigma ^(t)$:

Here, $\u27e8L\u27e9n0$ (the overall Liouvillian averaged over the initial state of the nuclear DOF, resulting in a super-operator in the electronic Liouville-subspace) and $K(\tau )$ (the memory kernel super-operator) are given by

and

respectively.

Importantly, evaluation of the Liouvillian and memory kernel terms in Eqs. (16) and (17), respectively, does not require casting the Hamiltonian in system-bath form or that the initial state of the nuclear DOF corresponds to thermal equilibrium with respect to the bath Hamiltonian. In what follows, we refer to Eq. (15) as the *modified GQME (M-GQME)*, in order to distinguish it from versions of the GQME which are based on casting the overall Hamiltonian in system-bath form and assuming that the initial state of the nuclear DOF corresponds to equilibrium with respect to the bath Hamiltonian (see Appendix B).

It should be noted that the expression for the memory kernel, Eq. (17), can be further simplified by introducing the Condon approximation, $V^jk\u2192Vjk$ (see Appendix A for proof),

where $Lzero(\u22c5)=\u2211j=1M\u0124j|jj|,\u22c5=\u0124zero,\u22c5$.

## III. EVALUATION OF THE M-GQME MEMORY KERNEL

Simulating the dynamics of the electronic DOF based on Eq. (15) requires knowledge of $\u27e8L\u27e9n0$ and $K(\tau )$. Obtaining $\u27e8L\u27e9n0$ requires the evaluation of the time-independent averages over the nuclear DOF at the initial time, $\u27e8\u0124j\u27e9n0$ and $\u27e8V^jk\u27e9n0$, which are relatively straightforward to calculate either fully quantum-mechanically, semiclassically, or fully classically.^{91}

Assuming that $\u27e8L\u27e9n0$ can be obtained, the memory kernel of the M-GQME, Eq. (17), is the key quantity needed for simulating the dynamics of the electronic DOF based on Eq. (15). However, unlike $\u27e8L\u27e9n0$, the evaluation of $K(\tau )$ is made challenging by the fact that it is time dependent. Furthermore, the time-dependence of $K(\tau )$ is nontrivial because it is dictated by the projection-dependent propagator, $e\u2212iQL\tau /\u210f$, rather than by the projection-independent propagator, $e\u2212iL\tau /\u210f$. One strategy for overcoming the latter difficulty is by using a scheme for evaluating $K(\tau )$ from projection-free inputs, i.e., inputs that involve $e\u2212iL\tau /\u210f$ rather than $e\u2212iQL\tau /\u210f$. Combining exact or approximate methods for evaluating those projection-free inputs can then lead to a methodology that can be applied to complex molecular systems.

A scheme for evaluating $K(\tau )$ from projection-free inputs can be developed by using the following general operator identity:^{64,79}

Substituting $A=L$ and $B=QL$ into Eq. (19), we obtain

Substituting Eq. (20) into Eq. (17) then leads to the following Volterra equation of the second-kind for $K(\tau )$:

Here,

Thus, given the projection-free quantity $F(\tau )$, Eq. (21) can be solved numerically for the projection-dependent $K(\tau )$ (see Appendix D). Hence, the problem of calculating $K(\tau )$ translates into that of calculating $F(\tau )$ [see Eq. (22)].

It should be noted that $F(\tau )=iU\u0307(\tau )$, where $U(\tau )$ is the time evolution operator of the electronic reduced density operator,

Thus, Eq. (21) can be rewritten in the following form:

This implies that the memory kernel of the M-GQME, $K(\tau )$, can be obtained directly from the time evolution operator of the reduced dynamics, $U(\tau )$. The choice between $F(\tau )$ and $U(\tau )$ as the projection-free input depends on the trade-off between evaluating the commutator explicitly [$F(\tau )$] and performing higher order time derivatives numerically [$U(\tau )$]. In what follows, we choose to proceed with $F(\tau )$.

In practice, both $K(\tau )$ and $F(\tau )$ are represented by *M*^{2} × *M*^{2} matrices in terms of the electronic basis {| *j*⟩|*j* = 1, …, *M*} (e.g., 4 × 4 matrices in the case of a system with two electronic states). The corresponding matrix elements of $F(\tau )$ are given by

Thus, $Fjk,uv(\tau )$ can be given in terms of expressions of the following form:

where the nuclear operator $\Gamma \u2009R^$ is either

$VjR^\u2212VkR^$, with

*a*=*j*and*b*=*k*,$VjaR^$, with

*a*≠*j*and*b*=*k*, or$VbkR^$, with

*a*=*j*and*b*≠*k*,

and terms with *a* ≠ *j* and *b* ≠ *k* do not occur.

The number of quantities of the form of Eq. (26) that needs to be calculated scales rather favorably with the dimensionality of the electronic Hilbert space (∼*M*^{4}). It should be noted that those quantities only need to be calculated once for a given initial state and that they can be calculated independently in a trivially parallelized manner.

## IV. CALCULATION OF $Fjk,uv(\tau )$ VIA THE EHRENFEST METHOD

The methodology outlined in Sec. III is general and can be used for calculating the memory kernel of the M-GQME via any exact or approximate method of one’s choice.^{79,82,90,92,93} In this section, we demonstrate this by outlining a protocol for calculating $Fjk,uv(\tau )$ via the Ehrenfest method. We start out by noting that $Fjk,uv(\tau )$ can be given in terms of expressions of the form of Eq. (26). The expression in Eq. (26) can be interpreted as the expectation value of $\Gamma \u2009R^|ba|$ at time *t* when the initial state is given by $\rho ^n(0)|uv|$ and the dynamics is dictated by the overall Hamiltonian, *Ĥ* [see Eq. (1)]. However, while $\rho ^n(0)$ is by definition a legitimate nuclear density operator, associating |*u*⟩⟨*v*| with an initial electronic density operator, $\sigma ^(0)$, is not possible when *u* ≠ *v*, since in this case |*u*⟩⟨*v*| is not Hermitian, has zero trace, and does not satisfy the Schwarz inequality [$|\sigma jk|\u2264\sigma jj\sigma kk1/2$].

The fact that |*u*⟩⟨*v*| is not a legitimate density operator when *u* ≠ *v* can become an obstacle when one attempts to evaluate Eq. (26) via semiclassical or mixed quantum-classical methods. In this section, we demonstrate this point in the context of the Ehrenfest method. Within this method, the nuclear DOF are treated classically, the electronic DOF are treated quantum-mechanically, and the impact of the electronic DOF on the nuclear DOF is treated in a mean-field manner. In practice, initial positions and momenta of the nuclear DOF are sampled based on the Wigner transform of $\rho ^n(0)$, *ρ*_{n,W}(**R**, **P**; *t* = 0), or its classical limit, and expectation values are obtained by averaging over the corresponding ensemble of classical trajectories. The effect of the nuclear DOF on the electronic DOF is accounted for by the fact that each classical trajectory of the nuclear DOF, **R**(*t*), corresponds to a different realization of the time-dependent Hamiltonian that governs the dynamics of the electronic density operator,

The effect of the electronic DOF on the nuclear DOF is accounted for by propagating the nuclear DOF on the mean-field PES,

Attempting to use the Ehrenfest method when the initial electronic density operator is non-Hermitian, e.g., $\sigma ^(0)=|uv|$ when *u* ≠ *v*, results in a complex mean-field PES, which in turn leads to nonphysical complex classical positions and momenta of the nuclear DOF. This problem can be bypassed by switching to a basis of the electronic Liouville space, consisting of operators that satisfy the conditions for a density operator. The choice of basis is not unique, but a relatively unbiased choice that satisfies hermiticity, trace one, and Schwarz inequality corresponds to

The results reported in this paper were obtained based on this choice. It should be noted that Montoya-Castillo and Reichman^{82} proposed an alternative approach for resolving the above mentioned discrepancies which is based on the identity |*u*⟩⟨*v*| + |*v*⟩⟨*u*| = *∑*_{k}*λ*_{k}|*λ*_{k}⟩⟨*λ*_{k}| and separately simulating Ehrenfest dynamics for each |*λ*_{k}⟩⟨*λ*_{k}| (here, {|*λ*_{k}⟩} are the eigenfunctions of |*u*⟩⟨*v*| + |*v*⟩⟨*u*| and {*λ*_{k}} are the corresponding eigenvalues).

In practice, one starts with $X^uv$ and *Ŷ*_{uv} instead of |*u*⟩⟨*v*| and |*v*⟩⟨*u*| as initial electronic states, to obtain the Ehrenfest approximations of

The corresponding results for |*u*⟩⟨*v*| and |*v*⟩⟨*u*| as the initial electronic states can then be expressed as linear combinations of the results in Eq. (30). For example,

## V. ILLUSTRATIVE APPLICATIONS

In this section, we demonstrate the applicability and robustness of the M-GQME by obtaining the memory kernel via Ehrenfest-based projection-free inputs and applying the M-GQME to a spin-boson model system with two electronic states [donor (*D*) and acceptor (*A*)], harmonic electronic PESs which are shifted in equilibrium energy and geometry, and an electronic coupling coefficient which is assumed to be constant (Condon approximation). The overall Hamiltonian is given by

where

Here, Γ is a positive constant, 2*ϵ* is the shift in equilibrium energy between the donor and acceptor states (*ϵ* = 0 and *ϵ* ≠ 0 correspond to the unbiased and biased cases, respectively), and $2cj/\omega j2$ is the corresponding shift in equilibrium geometry along the *j*th mode coordinate. Since this system satisfies the Condon approximation, we use the projection-free inputs $F1(\tau )$ and $F2(\tau )$ to obtain the memory kernel, as shown in Appendix A.

We also compare the results obtained based on the M-GQME scheme to those obtained from GQME-based schemes that start out with the overall Hamiltonian in a system-bath form. We used the Shi-Geva^{64} (SG-NZ) and Zhang-Ka-Geva^{79} (ZKG-NZ) system-bath-based schemes (see Appendix B), which have the form of Eq. (7) with the system, bath, and system-bath terms given by

The values of {*ω*_{j}} and {*c*_{j}} (*j* = 1, …, *N*) are obtained based on an Ohmic spectral density with an exponential cutoff

Here, *ξ* is the Kondo parameter (a measure of the shift in equilibrium geometry) and *ω*_{c} is the cutoff frequency. The discrete set of nuclear mode frequencies, {*ω*_{j}}, and coupling coefficients, {*c*_{j}}, were obtained from the continuous spectral density, Eq. (35), following the procedure described in Appendix C. The initial state of the nuclear DOF was chosen as

where *Ĥ*_{B} is as in Eq. (34) and the initial nuclear position and momenta are sampled based on the Wigner transform of Eq. (36),

It should be noted that this particular choice is dictated by our desire to compare the M-GQME scheme with the SG-NZ and ZKG-NZ schemes, which require that the initial nuclear state corresponds to thermal equilibrium with respect to the bath Hamiltonian, by definition, and to eliminate the inhomogeneous term, respectively. At the same time, it is also important to emphasize that the M-GQME is designed to accommodate arbitrary initial nuclear states of one’s choice.

Calculations were carried out for five different sets of parameter values (see Table I) averaged over 10^{5} trajectories. The memory kernel and projection-free input elements were found to have the following properties:

$Kjj,uv(\tau )=F1,jj,uv(\tau )=F2,jj,uv(\tau )=0\u2009$ and

$Kjk,uv(\tau )=Kkj,vu*(\tau ),F1,jk,uv(\tau )=F1,kj,vu*(\tau ),F2,jk,uv(\tau )=\u2212F2,kj,vu*(\tau ),$.

The nonvanishing matrix elements of the memory kernel super-operators for models *#*1–5, which were calculated using the Ehrenfest method, are shown in Figs. 1–5, respectively. The population difference between donor and acceptor states, which corresponds to the expectation value of $\sigma ^z=|D\u27e8D|\u2212|A\u27e9A|$, for models *#*1–5 is shown in Figs. 6–10, respectively. Exact results were adopted from Ref. 93 for models *#*1–4 and from Ref. 94 for model *#*5. The nonvanishing matrix elements of the projection-free input super-operators $F1(\tau )$ and $F2(\tau )$ for models *#*2 and *#*3 are given in Appendix A.

. | Model parameters . | Numerical parameters . | ||||||
---|---|---|---|---|---|---|---|---|

Model #
. | ϵ
. | Γ . | β
. | ξ
. | ω_{c}
. | ω_{max}
. | N
. | Δt
. |

1 | 1.0 | 1.0 | 5.0 | 0.1 | 1.0 | 5 | 400 | 0.02 |

2 | 1.0 | 1.0 | 5.0 | 0.1 | 2.0 | 10 | 400 | 0.02 |

3 | 1.0 | 1.0 | 5.0 | 0.1 | 7.5 | 36 | 400 | 0.02 |

4 | 1.0 | 1.0 | 5.0 | 0.4 | 2.0 | 10 | 400 | 0.02 |

5 | 0.0 | 0. $33$ | 3.0 | 0.1 | 1.0 | 5 | 200 | 0.02 |

. | Model parameters . | Numerical parameters . | ||||||
---|---|---|---|---|---|---|---|---|

Model #
. | ϵ
. | Γ . | β
. | ξ
. | ω_{c}
. | ω_{max}
. | N
. | Δt
. |

1 | 1.0 | 1.0 | 5.0 | 0.1 | 1.0 | 5 | 400 | 0.02 |

2 | 1.0 | 1.0 | 5.0 | 0.1 | 2.0 | 10 | 400 | 0.02 |

3 | 1.0 | 1.0 | 5.0 | 0.1 | 7.5 | 36 | 400 | 0.02 |

4 | 1.0 | 1.0 | 5.0 | 0.4 | 2.0 | 10 | 400 | 0.02 |

5 | 0.0 | 0. $33$ | 3.0 | 0.1 | 1.0 | 5 | 200 | 0.02 |

One observation that can be gleamed from Figs. 1–10 is that the M-GQME and ZKG-NZ schemes produce memory kernels that are better behaved at long time than those produced by the SG-NZ scheme. More specifically, with the exception of model #1, the memory kernels obtained via the SG-NZ scheme are observed to oscillate asymptotically, rather than vanish.

The instabilities of the Ehrenfest-based SG-NZ memory kernels have been reported in previous studies.^{82,93} In one previous study,^{93} they were dealt with by truncating the memory kernel at short times. This indeed reproduces the population dynamics reported in Ref. 93, which also happens to be in excellent agreement with the exact result (see Fig. 11). For example, in the case of model #4, this meant truncating the memory kernel at *t* = 1.5Γ^{−1}.^{93} However, truncating the memory kernel at *t* = 1.5Γ^{−1} also causes the M-GQME and ZKG-NZ to disagree with the exact result. A closer inspection of Fig. 4 reveals that the memory kernel is actually longer lived and that truncating it at *t* = 10.0Γ^{−1} would be more reasonable.^{95} Indeed, truncating the memory kernel at *t* = 10.0Γ^{−1} rather than at *t* = 1.5Γ^{−1} significantly improves the agreement between the population dynamics produced by M-GQME and ZKG-NZ and the exact result (see Fig. 9). At the same time, it also causes the population dynamics produced by the SG-NZ to oscillate asymptotically around the exact result, which is consistent with a similar observation made in Ref. 93. In another previous study,^{82} the memory time was determined by a “plateau of stability” found in the *σ*_{z}(*t*) dynamics with respect to the memory time. However, as shown in Fig. 12 and acknowledged in Ref. 82, this plateau can be short-lived or nonexistent. Additionally, determination of the plateau of stability without knowledge of the exact results can be challenging. In comparison, *σ*_{z}(*t*) dynamics within the M-GQME and ZKG-NZ schemes converges smoothly, as seen in Fig. 12, which makes finding a plateau of stability unnecessary. The M-GQME and ZKG-NZ convergence are obtained with a memory time equal to the lifetime of the memory kernel, e.g., *t*_{mem} = 10.0Γ^{−1} in the case of model #4.

Another observation is that the population dynamics produced by M-GQME, ZKG-NZ, and SG-NZ with memory kernels obtained via the Ehrenfest method are in much better agreement with the exact result than the population dynamics obtained via direct application of the Ehrenfest method. At first sight, this is somewhat surprising, given that the memory time, *t* = 10.0Γ^{−1}, is comparable to the population relaxation time scale. However, it should be noted that within the GQME, the effect of the density operator at time *t* − *τ* on its dynamics at time *t* decreases with increasing *τ* due to the finite lifetime of the memory kernel. Thus, as the Ehrenfest method becomes less accurate with increasing time, its contribution to the dynamics through $K(\tau )$ diminishes as well. As a result, using the Ehrenfest method to calculate the memory kernel leads to significantly more accurate results than using the Ehrenfest method to calculate the population dynamics directly. It should be noted that the authors of Ref. 82 also argued that the improvement of the GQME over direct Ehrenfest could be due to the additional information about the electronic-nuclear coupling gained through the sampling of nuclear operators within the inputs for the memory kernel.

Yet another interesting observation is the loss of accuracy and stability of the memory kernels with increasing cutoff frequency, *ω*_{c} (see Figs. 1–3 and 6–8). This can be traced back to the treatment of the nuclear DOF as classical within the Ehrenfest method. More specifically, increasing *ω*_{c} corresponds to increasing the frequency of the nuclear modes and thereby making the assumption that they can be treated classically less valid. Along with the increasing instability, another trend seen in Figs. 1–3 is that with the increasing cutoff frequency, the scale of the memory kernels also increases.

Finally, it is interesting to contrast the biased case (*ϵ* ≠ 0, Figs. 1–4) with the unbiased case (*ϵ* = 0, Fig. 5). While direct application of the Ehrenfest method appears to produce rather accurate results in the unbiased case, it is observed to lead to significant deviations in the biased case. This can be traced back to the Ehrenfest method’s failure to capture detailed balance. Interestingly, restricting the use of the Ehrenfest method to calculating the memory kernel and simulating the electronic dynamics through the GQME gives rise to significantly more accurate results in the biased case. It should be noted that given the quantum-mechanically exact memory kernel, the GQME is guaranteed to satisfy detailed balance since it corresponds to the exact equation of motion of the electronic DOF. The fact that it still does rather well even when the memory kernel is calculated via the Ehrenfest method should be considered as yet another advantage of the GQME approach.

## VI. CONCLUDING REMARKS

Although the system-bath paradigm has been a central theme in the study of quantum open systems, there are cases where it is not desirable to cast the overall Hamiltonian in system-bath form, Eq. (7). A prime example is nonadiabatic dynamics, where it is neither natural nor convenient to cast the Hamiltonian in terms of a system Hamiltonian, which only depends on the electronic DOF, a bath Hamiltonian, which only depends on the nuclear DOF, and a system-bath interaction term, which couples them. This is because the overall Hamiltonian underlying nonadiabatic dynamics associates a *different* nuclear Hamiltonian with each electronic state, thereby making the definition of a single bath Hamiltonian non-unique and essentially arbitrary. The lack of a unique system-bath Hamiltonian can be particularly problematic when approximate methods are used to evaluate the memory kernel, as would often be the case when dealing with realistic molecular models, since different choices of bath Hamiltonian may lead to different results without a clear criterion for choosing between them.

In this paper, we utilized the fact that the GQME, which represents the exact equation of motion of the electronic DOF during nonadiabatic dynamics, does not in fact need to be based on casting the overall Hamiltonian in system-bath form.^{82,90,96} We refer to this form of the GQME as the M-GQME. We also presented a practical scheme for calculating the memory kernel of the M-GQME, either exactly or approximately, that does not rely on the system-bath form. In doing so, we end up with a natural and convenient GQME-based approach for simulating the dynamics of the electronic DOF during nonadiabatic dynamics.

It should be noted that unlike other methods for simulating nonadiabatic dynamics, such as Ehrenfest, surface hopping, MQCL, and SQC, the approach based on the M-GQME is focused on the dynamics of the electronic DOF. The dynamics of the nuclear DOF is only captured to the extent that it impacts the electronic DOF. The memory kernel represents the minimum input of the nuclear DOF that is required in order to account for their effect on the dynamics of the electronic DOF. In this respect, the GQME can be thought of as going beyond approaches based on FGR, where the coupling between the electronic and nuclear DOF is assumed weak and the impact of the nuclear DOF on the electronic DOF is captured by the two-time autocorrelation function of the coupling between nuclear and electronic DOF.^{1,65–71} However, unlike FGR-based approaches, the GQME does not require assuming weak coupling between electronic states and describes the electronic dynamics in terms of the full electronic density matrix, rather than in terms of the electronic populations, which correspond to its diagonal elements.

On the one hand, the loss of more detailed information on the dynamics of the nuclear DOF may be viewed as a disadvantage of the GQME-based approach to nonadiabatic dynamics. On the other hand, focusing on the memory kernel rather than on a complete description of the nuclear DOF offers several important advantages. First, it is often the case that the only interesting aspect of the nuclear dynamics is its impact on the electronic dynamics. Thus, the compactness of the memory kernel provides an elegant way for focusing on this aspect without needing to figure out whether or not a given detail of the nuclear dynamics impacts the electronic DOF. Second, the compactness of the memory kernel and its finite memory time also imply that calculating it via a given method, either exact or approximate, can be more cost-effective and/or lead to more accurate results than a direct application of the same method. Third, it should be remembered that most useful approximate methods describe nuclear dynamics in terms of an ensemble of classical-like trajectories and are constructed in such a way that only the ensemble average, rather than individual trajectories, can be related to physically meaningful measurable quantities like electronic populations and coherences. The fact that the memory kernel is defined in terms of a trace over the nuclear DOF implies that it incorporates this ensemble-averaging automatically and is therefore directly related to the only relevant measurable quantities.

In summary, the M-GQME provides a rigorous and general approach for simulating electronically nonadiabatic dynamics. Within this approach, the memory kernel super-operator is the key quantity which dictates how the electronic dynamics is impacted by the nuclear DOF, regardless of the strength or type of coupling. What makes this approach particularly appealing is the fact that calculating the memory kernel via exact or approximate methods can be more cost-effective and/or accurate than direct application of those methods. Given the non-uniqueness associated with the choice of basis in Eq. (29), which appears to be inherent to the Ehrenfest method, it would also be highly desirable to explore calculating the memory kernel via approximate methods other than the Ehrenfest method and apply the GQME-based approach to systems other than the spin-boson model. Work on such extensions is underway and will be reported in future publications.

## ACKNOWLEDGMENTS

A.S. is grateful for support from an ICAM fellowship, awarded by the Kent State University and University of Michigan ICAM branches. B.D.D. is grateful for support from the NSF via Grant No. CHE-1362504. E.G. is grateful for support from the NSF via Grant Nos. CHE-1464477, CHE-1800325, and CHE-1663636. B.D.D. and E.G. are grateful for support from the DOE via Award No. DE-SC0016501. The computational resources and services were provided by Advanced Research Computing at the University of Michigan, Ann Arbor.

### APPENDIX A: THE MEMORY KERNEL OF THE M-GQME IN THE CONDON APPROXIMATION

In this appendix, we show that within the Condon approximation, $V^jk\u2192Vjk$, the memory kernel in Eq. (17) can be simplified into that in Eq. (18). To this end, let $L=Lzero+Lint$, where

[see Eq. (1)] and note that $Lint$ becomes a purely electronic super-operator in the Condon approximation. As a result,

This implies that one can replace $L$ by $Lzero$ on the left side of the exponent in Eq. (17).

$L$ can also be replaced by $Lzero$ on the right side of the exponent in Eq. (17) since

Therefore, in the Condon limit,

From here, substituting Eq. (20) into Eq. (A4) leads to the following Volterra equation of the second-kind for $K(\tau )$ in the Condon limit:

Here,

In contrast to the memory kernel, $F1(\tau )$ and $F2(\tau )$ are not required to have finite lifetimes. Examples of $F1(\tau )$ and $F2(\tau )$ for models #2 and #3 in Table I are given in Figs. 13–16.

### APPENDIX B: DIFFERENT SCHEMES FOR EVALUATING THE MEMORY KERNEL OF THE GQME WHEN THE OVERALL HAMILTONIAN IS IN A SYSTEM-BATH FORM

In this appendix, we outline previously proposed schemes for evaluating the memory kernel of the GQME for a system with an overall Hamiltonian in a system-bath form [see Eq. (7)]. The reader is referred to Refs. 64, 79, and 86 for further details. These schemes have also been studied previously by Markland *et al.*,^{93,96} Montoya-Castillo and Reichman,^{82,90} and Rabani *et al.*^{81,95,97–99}

In those schemes, one often starts with the overall Hamiltonian in the system-bath form of Eq. (7) and an initial state which is given by Eq. (2),

It is also often assumed, without loss of generality, that *Ĥ*_{BS} is defined such that

Using a projection operator of the form

the quantum-mechanically exact dynamics of the system reduced density operator can be shown to be governed by a GQME of the following form:

Here, $\u2212iLS\sigma ^(t)/\u210f=\u2212i[\u0124S,\sigma ^(t)]/\u210f$ and $\u2212\u222b0td\tau K(\tau )\sigma ^(t\u2212\tau )$ correspond to the bath-free and bath-induced contributions to the system’s reduced dynamics, respectively.

The bath-induced component is dictated by the memory kernel super-operator, $K(\tau )$, which, under the above mentioned conditions, can be written in a variety of different, yet equivalent, forms

Here, $LBS(\u22c5)=[\u0124BS,\u22c5]$ and $Q=1\u2212P$. Appendixes B 1 and B 2 outline two of the previously proposed schemes for calculating $K(\tau )$, which we will be comparing to in this paper. The differences between the schemes can be generally traced back to which of the forms of the memory kernel, Eq. (B5), is chosen as the starting point.

#### 1. The Shi-Geva (SG-NZ) scheme

The original scheme for calculating the memory kernel of the GQME^{64} was based on the following expression for the memory kernel [see Eq. (B5)]:

Substituting $A=L$ and $B=L\u2212LBSP$ into the operator identity in Eq. (19), one can then obtain

Here,

and

are auxiliary kernels that are needed in order to calculate the memory kernel.

It should be noted that unlike $K1(\tau )$, which is projection-free, $K2(\tau )$ is projection-dependent. As such, calculating $K2(\tau )$ involves a similar challenge to that of calculating $K(\tau )$. However, $K2(\tau )$ can be evaluated from the following Volterra equation, obtained by substituting (B7) into Eq. (B10):

Here,

#### 2. The Zhang-Ka-Geva (ZKG-NZ) scheme

An alternative scheme for calculating the memory kernel of the GQME^{79} was based on writing the memory kernel in the following form [see Eq. (B5)]:

Here, Φ(*τ*) is a projection-free input,

### APPENDIX C: BATH DISCRETIZATION PROCEDURE

The discrete set of *N* nuclear mode frequencies, {*ω*_{1}, …, *ω*_{N}}, and coupling coefficients, {*c*_{1}, …, *c*_{N}}, for the Ohmic spectral density with exponential cutoff, Eq. (35), were obtained following the procedure described in Ref. 100,

Here,

where *ω*_{N} = *ω*_{max} is the frequency of the highest frequency mode.

The value of *ω*_{max} is determined using the following procedure. First we define the parameter *α*,

The parameter *α* controls the discretized spectral density. A value close to one yields a spectral density that covers high frequencies but at the cost of an overall coarse-grained frequency distribution. This could be compensated by an increased number of modes at the expense of increased computational costs. In practice, the actual value of *α* is determined in a manner that balances accuracy with cost. For the calculations reported in this paper, the value of *α* was set to 0.95.

Given the value of *α*, *ω*_{max} can be determined from Eq. (C3),

where *W*(*x*) is the Lambert *W* function, *W*(*xe*^{x}) = *x*. *W*(*x*) was calculated numerically using the python library function *scipy*.*special*.*lambertw*(*x*, *i*, *tol*) with *x* = (*α* − 1)/*e*, *α* = 0.95, *i* = −1, and *tol* = 10^{−6}, with the result rounded to the next whole integer for *ω*_{max}.

### APPENDIX D: NUMERICAL SOLUTION OF EQ. (21)

In this appendix, we outline the iterative algorithm used for solving Eq. (21) numerically. It should be noted that Eq. (21) is a Volterra equation of the second kind and as such has the following general form:

Given that *h*(*t*, *s*) and *g*(*t*) are known, this equation is solved for *f*(*t*). Comparing Eq. (D1) with Eq. (21) shows that in our case, this is an operator equation with *t*_{0} → 0, *t* → *τ*, *s* → *τ*′, $f(t)\u2192K(\tau )$, $h(t,s)\u2192iF(\tau \u2212\tau \u2032)$, and $g(t)=iF\u0307(\tau )\u22121\u210fF(\tau )\u27e8L\u27e9n0$. The iterative algorithm starts out with *f*(*t*) = *g*(*t*) as the initial guess. Substituting this initial guess into Eq. (D1) generates another estimator of *f*(*t*), which is then substituted back into Eq. (D1). This procedure is repeated until convergence, where the estimators obtained in two subsequent steps are indistinguishable within a prescribed tolerance.

In practice, *f*(*t*) is a matrix and time is discretized, *t*_{n} = *n*Δ*t* with *n* = 0, 1, 2, …, *N*. Let $fjki(n\Delta t)$ be the value of ( *j*, *k*)-th matrix element of *f* after the *i*th iteration,

The time integral in Eq. (D2) is calculated via the extended trapezoidal rule. The results reported in this paper were obtained using the following criterion for convergence: $fjki(n\Delta t)\u2212fjki\u22121(n\Delta t)\u2009\u226410\u221210$ (for all time steps *n* and matrix elements *jk*). For the applications reported in this paper, the typical number of iterations necessary for obtaining converged results was 4-16.