A modified approach for simulating electronically nonadiabatic dynamics via the generalized quantum master equation

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,²³ 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, $Ĥj=P^2/2+VjR^$ is the nuclear Hamiltonian when the system is in the electronic state | j⟩, $R^=R^1,…,R^N$ and $P^=P^1,…,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⁸⁵). 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.⁸⁵ 

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²) 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, $ρ^n(0)=Tre{ρ^(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(⋅)=[Ĥ,⋅]$ 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 $σ^(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 $ρ^(0)=pρ^na(0)⊗σ^a(0)+(1−p)ρ^nb(0)⊗σ^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, $σ^a(t)$ and $σ^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, $σ^(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, $ρ^(0)$⁠, can be expressed as $∑j,kρ^jk(0)|jk|$⁠, with $ρ^jk(0)=⟨j|ρ^jk(0)|k⟩$⁠. Thus, the computational cost scales linearly with respect to the number of product terms with independent $ρ^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 $ρ^n(0)=ρ^Beq=exp[−βĤB]/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., $ρ^n(0)$⁠, do. Third, the system-bath coupling, Ĥ_BS, is often defined so that $TrB{ĤBSρ^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(⋅)=ρ^Beq⊗TrB{⋅}$ and would also need to be modified when the definition of Ĥ_B changes. Fifth, the second term in Eq. (1), $∑j,k≠jV^jk|jk|$⁠, is purely electronic in the Condon approximation, $V^jk→Vjk$⁠, 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.

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ρ^(t)$⁠, for any projection operator $P$ that satisfies idempotence (⁠$P2=P$⁠) is given by the NZ-GQME^1,2,87,88

Here, $Q=1−P$ 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, $ρ^nref$ is a reference nuclear density operator of one’s choice (as long as $Trnρ^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 $σ^(t)$⁠,

Here,

accounts for the Hamiltonian and Markovian dynamics generated by the Hamiltonian $⟨Ĥ⟩nref=TrnĤρ^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(τ)$⁠, 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(τ)$ 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 $ρ^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 $ρ^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 $ρ^nref$ in a manner that will simplify the resulting GQME. The assignment $ρ^nref=ρ^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.⁷⁹ 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 [$ρ^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 $ρ^n(0)$ is dictated by how the system is prepared, which is ultimately dependent on the experimental setup. Importantly, $ρ^n(0)$ need not be of the form $ρ^Beq=ZB−1e−βĤB$⁠. It should be noted that this is also not required within the Zhang-Ka-Geva (ZKG-NZ) scheme (see  Appendix B).⁷⁹ 

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

Here, $⟨L⟩n0$ (the overall Liouvillian averaged over the initial state of the nuclear DOF, resulting in a super-operator in the electronic Liouville-subspace) and $K(τ)$ (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→Vjk$ (see  Appendix A for proof),

where $Lzero(⋅)=∑j=1MĤj|jj|,⋅=Ĥzero,⋅$⁠.

Simulating the dynamics of the electronic DOF based on Eq. (15) requires knowledge of $⟨L⟩n0$ and $K(τ)$⁠. Obtaining $⟨L⟩n0$ requires the evaluation of the time-independent averages over the nuclear DOF at the initial time, $⟨Ĥj⟩n0$ and $⟨V^jk⟩n0$⁠, which are relatively straightforward to calculate either fully quantum-mechanically, semiclassically, or fully classically.⁹¹ 

Assuming that $⟨L⟩n0$ 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 $⟨L⟩n0$⁠, the evaluation of $K(τ)$ is made challenging by the fact that it is time dependent. Furthermore, the time-dependence of $K(τ)$ is nontrivial because it is dictated by the projection-dependent propagator, $e−iQLτ/ℏ$⁠, rather than by the projection-independent propagator, $e−iLτ/ℏ$⁠. One strategy for overcoming the latter difficulty is by using a scheme for evaluating $K(τ)$ from projection-free inputs, i.e., inputs that involve $e−iLτ/ℏ$ rather than $e−iQLτ/ℏ$⁠. 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(τ)$ 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(τ)$⁠:

Here,

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

It should be noted that $F(τ)=iU̇(τ)$⁠, where $U(τ)$ 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(τ)$⁠, can be obtained directly from the time evolution operator of the reduced dynamics, $U(τ)$⁠. The choice between $F(τ)$ and $U(τ)$ as the projection-free input depends on the trade-off between evaluating the commutator explicitly [$F(τ)$] and performing higher order time derivatives numerically [$U(τ)$]. In what follows, we choose to proceed with $F(τ)$⁠.

In practice, both $K(τ)$ and $F(τ)$ are represented by M² × M² 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(τ)$ are given by

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

where the nuclear operator $ΓR^$ is either

• $VjR^−VkR^$⁠, 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⁴). 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.

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(τ)$ via the Ehrenfest method. We start out by noting that $Fjk,uv(τ)$ 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 $ΓR^|ba|$ at time t when the initial state is given by $ρ^n(0)|uv|$ and the dynamics is dictated by the overall Hamiltonian, Ĥ [see Eq. (1)]. However, while $ρ^n(0)$ is by definition a legitimate nuclear density operator, associating |u⟩⟨v| with an initial electronic density operator, $σ^(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 [$|σjk|≤σjjσ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 $ρ^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., $σ^(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⁸² 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,

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/ωj2$ is the corresponding shift in equilibrium geometry along the jth mode coordinate. Since this system satisfies the Condon approximation, we use the projection-free inputs $F1(τ)$ and $F2(τ)$ 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⁶⁴ (SG-NZ) and Zhang-Ka-Geva⁷⁹ (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⁵ trajectories. The memory kernel and projection-free input elements were found to have the following properties:

• $Kjj,uv(τ)=F1,jj,uv(τ)=F2,jj,uv(τ)=0$ and

• $Kjk,uv(τ)=Kkj,vu*(τ),F1,jk,uv(τ)=F1,kj,vu*(τ),F2,jk,uv(τ)=−F2,kj,vu*(τ),$⁠.

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 $σ^z=|D⟨D|−|A⟩A|$⁠, 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(τ)$ and $F2(τ)$ for models #2 and #3 are given in  Appendix A.

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,⁹³ 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Γ⁻¹.⁹³ However, truncating the memory kernel at t = 1.5Γ⁻¹ 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Γ⁻¹ would be more reasonable.⁹⁵ Indeed, truncating the memory kernel at t = 10.0Γ⁻¹ rather than at t = 1.5Γ⁻¹ 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,⁸² 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Γ⁻¹ 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Γ⁻¹, 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(τ)$ 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.

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.

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.

In this appendix, we show that within the Condon approximation, $V^jk→Vjk$⁠, 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