A comparative study of different methods for calculating electronic transition rates

Nonadiabatic dynamics plays a central role in a variety of important chemical processes that range from redox reactions to photovoltaics.^1–12 The inherently quantum nature of the underlying molecular dynamics makes simulating such processes in complex molecular systems challenging.^2,13–21 One approach for addressing this challenge is based on employing a mixed quantum-classical strategy, where the electronic degrees of freedom (DOF) are treated quantum-mechanically, while the nuclear DOF are treated in a classical-like manner. A variety of such mixed quantum-classical methods has been employed for this purpose over the years. The Ehrenfest method,²² which treats the feedback between electronic and nuclear DOF in a mean-field manner, is arguably the simplest among them. The fewest switches surface-hopping (FSSH) approach,^23,24 in its many different implementations,^25–34 represents another family of widely used mixed quantum-classical methods that attempt to go beyond the mean-field approximation. A more rigorous approach is based on the mixed quantum-classical Liouville (MQCL) equation,^35–39 which can be solved, e.g., via the trajectory-based Sequential Short-Time Propagation (SSTP)⁴⁰ or Trotter-Based Surface-Hopping (TBSH)⁴¹ algorithms.

An alternative approach for estimating the rate of nonadiabatic transitions is provided by Fermi’s golden rule (FGR), which is based on treating the electronic coupling within the framework of second-order time-dependent perturbation theory. The equilibrium FGR (E-FGR) allows one to calculate electronic transition rate constants when the system starts out at equilibrium on the donor state and electronic transitions occur on a time scale that is slower than the lifetime of the electronic coupling autocorrelation function.^10,42–45 The nonequilibrium FGR (NE-FGR) extends the applicability of FGR to nonequilibrium initial states and to situations where the above-mentioned separation of time scales is not necessarily valid.⁴⁶ Recently proposed implementations of E-FGR and NE-FGR within the framework of the linearized semiclassical (LSC) approximation have given rise to trajectory-based algorithms for calculating E-FGR and NE-FGR electronic transition rates in complex molecular systems.^47–51 

In this paper, we present a comprehensive comparison of the four methods, NE-FGR, MQCL, Ehrenfest, and FSSH, for estimating electronic transition rates. To this end, we compare the donor state population relaxation dynamics for the Garg-Onuchic-Ambegaokar (GOA) benchmark charge-transfer model,⁵² as obtained via these four methods, over a broad range of temperatures and electronic coupling strengths, with different nonequilibrium initial states, in the normal and inverted regimes. In Sec. II, we briefly review the NE-FGR, MQCL, Ehrenfest, and FSSH methods. In Sec. III, we describe the model. Simulation techniques are outlined in Sec. IV. Results are presented and discussed in Sec. V. Conclusions are provided in Sec. VI.

We consider a two-state system with the overall Hamiltonian $Ĥ= H ^ 0 + H ^ I$⁠, where the zero Hamiltonian, $H ^ 0$⁠, and the perturbation, $H ^ I$⁠, are given by

Here, $|D$ and $|A$ denote the donor and acceptor electronic states (assumed independent of the nuclear DOF), respectively, $H ^ D ∕ A$ are the corresponding nuclear Hamiltonians,

and $Γ ^ D A$ is the electronic coupling between donor and acceptor states (assumed constant from now on, within the Condon approximation, $Γ ^ D A ≡ Γ D A = Γ A D *$⁠). Here, $R ^ = R ^ 1 , … , R ^ N$ and $P ^ = P ^ 1 , … , P ^ N$ are the mass-weighted nuclear coordinates and momenta, and $V D ∕ A R ^$ are the potential energy surfaces (PESs) of the donor (D) and acceptor (A) states.

The system is assumed to start out at the donor electronic state with the nuclear DOF described by the density operator $ρ ^ D ( 0 )$⁠,

The donor state population at time $t>0$ is then given by

where $T r [ ⋯ ] = T r N T r e [ ⋯ ]$ denotes the trace over both nuclear (N) and electronic (e) Hilbert spaces. The different methods described below aim at employing a mixed quantum-classical approach to evaluate P_D(t) in complex molecular systems.

The NE-FGR approach is based on the assumption that the electronic coupling term, $H ^ I$⁠, constitutes a small perturbation in comparison to the zero Hamiltonian, $H ^ 0$⁠. Within the framework of second-order perturbation theory, the time evolution operator in the interaction picture is given by⁵³ 

where $exp +$ stands for positively time-ordered exponential, and the perturbation operator in the interaction picture is given by

Substituting Eq. (6) and $exp ( − i Ĥ t ∕ ℏ ) =exp ( − i H ^ 0 t ∕ ℏ ) U ^ I ( t )$ into Eq. (5), we obtain

Next, inserting Eq. (7) into Eq. (6) and truncating the expansion after the second-order term lead to the following approximate expression for P_D(t):⁴⁶ 

where the time-dependent rate coefficient $k ( t ′ )$ is given by

Here,

The NE-FGR, Eq. (10), is based on treating the nuclear DOF fully quantum-mechanically. As such, its use is limited to the very few model systems where exact fully quantum-mechanical real-time dynamics is accessible. Thus, from the point of view of computational feasibility and tractability, it is highly desirable to formulate an approximate NE-FGR which would be based on treating the nuclear DOF in a classical-like manner. However, doing so in a self-consistent manner is nontrivial since the $τ$-dependence of $C ( t ′ , τ )$⁠, Eq. (11), lacks a well-defined classical limit.⁴⁸ The LSC approximation^47,48 provides a rigorous approach for obtaining a uniquely defined classical-like expression for $C ( t ′ , τ )$⁠,^49,51

Here,

The evaluation of $C L S C ( t ′ , τ )$ calls for sampling initial nuclear coordinates and momenta, (R₀, P₀), based on the Wigner phase-space density of $ρ ^ D ( 0 )$⁠, $ρ ^ D ( 0 ) W ( R 0 , P 0 )$⁠, followed by propagating classically forward in time on the donor PES, V_D(R), from time 0 to time $t ′$⁠, and then backward in time on the average PES, V_av(R) = [V_D(R) + V_A(R)]/2, from time $t ′$ to time $t ′ −τ$⁠. $∫ t ′ t ′ − τ dtU ( R t a v )$ is evaluated during the $t ′ → t ′ −τ$ segment of the trajectory, where U(R_t) = V_D(R_t) − V_A(R_t) is the instantaneous potential energy gap between the donor and acceptor states. $C L S C ( t ′ , τ )$ is then obtained by averaging $exp − i ℏ ∫ t ′ t ′ − τ d t U ( R t a v )$ over a large number of trajectories until a converged average is obtained.

The derivation of the MQCL equation^37,54 starts out by performing a partial Wigner transform of the overall (electronic+nuclear) density operator, $ρ ^ ( t )$⁠, with respect to the nuclear DOF,

It should be noted that $ρ ^ W ( R , P , t )$ is an operator in the electronic DOF and a classical-like phase-space distribution in the nuclear DOF. The MQCL equation is then obtained by expanding the exact equation of motion for $ρ ^ W ( R , P , t )$ to first order in m_e/M_N, where m_e and M_N are the electronic and nuclear masses, respectively,

Here, the square brackets $⋅ , ⋅$ represent the commutator, the curly braces ${ ⋅ , ⋅ }$ denote the classical Poisson bracket, $H ^ W$ is the partial Wigner transform of the overall Hamiltonian [Eq. (1) + Eq. (2)], and $iL$ is the quantum-classical Liouville operator

The MQCL equation is given in an operator form in Eq. (15) and is therefore representation-independent. In practice, it is convenient to rewrite it in a matrix form that depends on the choice of a basis for the Hilbert space of the electronic DOF. To this end, the adiabatic basis is a convenient choice. The adiabatic basis vectors, ${ | α ; R }$⁠, are defined as the eigenstates of the electronic Hamiltonian,

Within the adiabatic representation, $−iL$ is given by an $N e 2 × N e 2$ matrix (N_e is the dimension of the electronic DOF Hilbert space). The corresponding matrix element is given by

Here, $δ ( α , β )$ is the Kronecker delta, $ω α β ( R ) = ( E α ( R ) − E β ( R ) ) ∕ℏ≡Δ E α β ∕ℏ$⁠, $i L α α ′$ is a classical-like Liouville operator

where

and

where $S α β =Δ E α β d α β R ̇ ⋅ d α β − 1$ and $d α β =α;R|∂∕∂R|β;R$⁠. Importantly, $J α α ′ , β β ′$ accounts for nonadiabatic transitions.

Several trajectory-based numerical algorithms for solving the MQCL equation in the adiabatic representation have been proposed. The TBSH algorithm is arguably the most accurate and numerically stable to date.⁴¹ Within TBSH, the expectation value of an overall system operator, $Â$⁠, at time $t N =NΔt$⁠, where $Δt$ is the time step, is given by

Here, the prefactor $N e 2$ arises from the uniform sampling over all initial quantum states, $( α 0 β 0 )$⁠, $κ$ denotes the trajectory index and results are averaged over N_traj trajectories. To simplify the notation, the quantum state of trajectory $κ$ at the jth time step $( α β ) j κ$ is replaced by a single index $s j κ = α j κ N e + β j κ$ and at j = 0 by $s ̃ 0 κ = β 0 κ N e + α 0 κ$⁠, respectively. The quantity $W s j ( t j , t j + Δ t ∕ 2 ) =exp i ∫ t t + Δ t ∕ 2 d τ ω s j ( R s j ( τ ) )$ stands for the phase factor and $R s j ( τ )$ are the coordinates of the classical DOF at time $τ$⁠. Finally, Q₁ is the matrix of transition probabilities and $M$ is the matrix containing the momentum shifts for every pair of quantum states. Expressions for Q₁ and $M$ can be found in Ref. 41.

Evaluation of Eq. (22) involves Monte Carlo sampling of the initial quantum state of the system, $s ̃ 0 = β 0 N e + α 0$⁠, and classical DOF, (R₀, P₀), adiabatic propagation for a half-time step $e i L s ̃ 0 Δ t ∕ 2 ( R , P ) → ( R ′ , P ′ )$⁠, and a hopping attempt to another quantum state $s 1 = α 1 N e + β 1$ with probability $| Q 1 s ̃ 0 s 1 Δ t |∕| ∑ s Q 1 s ̃ 0 s Δ t |$⁠. If the kinetic energy along the direction of coupling vector $d α β$ is sufficient, a momentum jump is performed, $( R ′ , P ′ ) → ( R ′ , P ′ + Δ P α β )$⁠, with

where $d ¯ α β = d α β ∕| d α β |$⁠. Following this, the classical DOF are propagated for the second half-time step, $e i L s 1 Δ t ∕ 2 ( R ′ , P ′ + Δ P α β ) → ( R ″ , P ″ )$⁠. The procedure is then repeated for $N−1$ time steps.

The Ehrenfest and FSSH methods are two other commonly used mixed quantum-classical methods.⁵⁵ For both methods, one starts out by assuming that the nuclear DOF can be treated classically, while the electronic DOF are treated quantum-mechanically. The nuclear DOFs follow a yet to be specified classical trajectory, {R(t), P(t)} (see below). The dynamics of the electronic DOF is governed by an explicitly time-dependent quantum-mechanical Hamiltonian, which coincides with $h ^ W [ R ( t ) ]$ [see Eq. (17)]. The time-dependence of the electronic Hamiltonian originates from the coupling between the nuclear and electronic DOF, through the time-dependence of the classical nuclear coordinates, R(t).

Assuming that the initial electronic state is pure (e.g., $|D$ for the case considered in this paper) and therefore can be described as a ket, $|ψ ( 0 )$⁠, the electronic state at time t is also given by a ket, $|ψ ( t )$⁠, which is obtained by solving the following electronic time-dependent Schrödinger equation:

Within the Ehrenfest and FSSH methods, initial conditions of the nuclear DOF, {R(0), P(0)}, are sampled from an initial phase-space probability density. Different choices of initial conditions result in different realizations of the nuclear trajectory, R(t), and therefore give rise to an ensemble of different independent nuclear trajectories, ${ R κ ( t ) }$⁠. This in turn will give rise to an ensemble of time-dependent electronic Hamiltonians, ${ h ^ W [ R κ ( t ) ] }$⁠, and thereby an ensemble of different electronic states at time t, ${ | ψ κ ( t ) }$⁠. The expectation value of any observable quantity, $Â$⁠, is then obtained by ensemble averaging over independent trajectories

Here, $A ^ W [ R , P ]$ is the partial Wigner transform of $Â$ with respect to the nuclear DOF (an electronic operator which depends parametrically on the nuclear phase-space variables). The independent trajectory approximation is made for both the Ehrenfest and the FSSH methods, throughout the manuscript. Alternative approaches that include cross talk between trajectories have been proposed but are beyond the scope of this paper (e.g., Ref. 56). Thus, all quantities depending on $R κ ( t )$ and $P κ ( t )$⁠, respectively, are evaluated along single trajectories. For readability, we drop the trajectory index $κ$ in the following.

In practice, the electronic ket, $|ψ ( t )$⁠, is represented in terms of a basis of one’s choice, ${ | ϕ a [ R ( t ) ] }$⁠,

The basis states may or may not be R-dependent. Furthermore, the R-dependence does not need to be the same as that of the adiabatic basis states (i.e., eigenstates of $h ^ W [ R ]$⁠).⁵⁷ Within the ${ | ϕ a [ R ( t ) ] }$ basis representation, the electronic time-dependent Schrödinger equation, Eq. (24), is given by

It should be noted that when the basis is adiabatic, the first term on the R.H.S. of Eq. (26) becomes diagonal (⁠$ϕ b | h ^ W | ϕ a =δ ( a , b ) ϕ a | h ^ W | ϕ a$⁠), while the second term becomes purely off-diagonal (⁠$ϕ b | ϕ ̇ a = [ 1 − δ ( a , b ) ] ϕ b | ϕ ̇ a$⁠) and accounts for nonadiabatic transitions. In contrast, when the basis is R-independent, the second term on the R.H.S. of Eq. (26) vanishes and the first term gives rise to nonvanishing diagonal and off-diagonal contributions.²³ 

The fundamental difference between the Ehrenfest and the FSSH methods is the way the classical nuclear trajectories, {R(t), P(t)}, are generated. Within the Ehrenfest method, {R(t), P(t)} are obtained via the mean-field approximation, i.e., by assuming that the classical-like force is given by $F M F ( t ) =− ψ ( t ) ∂ h ^ W ∂ R ψ ( t )$⁠. Note that for each trajectory, $|ψ ( t )$ is the full electronic state vector so that each trajectory’s F^MF(t) becomes independent of the chosen electronic basis. Thus, the dynamics of the nuclear DOF are governed by the following classical-like equations of motion:

Within FSSH, the nuclear DOF are subject to a classical-like force which is dictated by one of the basis functions, the so-called active state a, $F a S H ( t ) =− ϕ a [ R ( t ) ] ∂ h ^ W ∂ R ϕ a [ R ( t ) ]$ so that the dynamics of the nuclear DOF are governed by the following classical-like equations of motion:

Which basis function is used at a given time is determined by a stochastic algorithm that allows for hopping between electronic states as a function of time. The stochastic algorithm is set up so that the fraction of trajectories that employ the state $| ϕ a$ at time t, n_a(t)/N_traj, are equal to the electronic population of this state, |c_a(t)|², as obtained from solving the electronic time-dependent Schrödinger equation.^23,58 This is achieved by assuming that the transition probability at time t from electronic state $| ϕ a$ to electronic state $| ϕ b$ is given by

Here, $Δt$ is the time step used for propagating the nuclear DOF,

and $Θ ( x )$ is the Heaviside function $Θ ( x < 0 ) = 0$ and $Θ ( x ≥ 0 ) = 1$⁠, which guarantees positivity of the transition probability. It should be noted that the first term in Eq. (30) vanishes when the adiabatic basis is used, whereas the second term vanishes when using an R-independent basis.

In practice, a random number $ζ∈ [ 0,1 ]$ is generated from a uniform distribution at each point in time. A hop changing the current active electronic state, $| ϕ a$⁠, to a different state, $| ϕ b$⁠, only takes place if the following condition is satisfied:

Within FSSH, a hopping event between electronic states is typically accompanied by a change in the nuclear momenta (so-called momentum adjustment). This appears to be based on the somewhat ad-hoc expectation that energy needs to be conserved along individual nuclear trajectories. However, it should be noted that a surface hopping scheme that avoids this constraint was recently proposed in Ref. 56. However, if energy conservation is imposed, the corresponding shift in the kinetic energy needs to compensate for the change in the potential energy. When the adiabatic representation is used, it was argued that the momentum adjustment to each of the nuclear DOF needs to be proportional to the projection of the nonadiabatic coupling vector, $d a b ≡ ϕ a ( R ) | ∂ ∂ R | ϕ b ( R )$⁠, onto this DOF⁵⁸ 

where $γ$ is dictated by the above-mentioned energy conservation.

When a representation in terms of an R-independent basis is used, it was proposed that the momentum adjustment is introduced via a uniform velocity rescaling^59,60

Finally, we note that if energy is to be conserved along individual trajectories, a hop to an energetically higher state can only occur if sufficient classical kinetic energy is available. Otherwise, the transition needs to be aborted. In the case of such frustrated hops, some recent studies proposed a partial momentum reversal.^61,62 However, we found that such modifications have a negligible effect on the results presented below. Hence, in the results reported below, the momentum remains unchanged during frustrated hops, in accordance with earlier studies.⁶³ 

In this paper, both basis schemes are examined, i.e., active surfaces are either given by R-independent states {D, A} (FSSH-D) or by adiabatic states {g, e} (FSSH-A), respectively. The quantities of interest in this paper are populations of R-independent states. Since the calculation of R-independent state populations from adiabatic active states is ambiguous,^64–66 we report in the following only population measures based on wavefunction coefficients as detailed in Sec. IV. A brief discussion about alternative population measures can be found in the  Appendix.

The main goal of this paper is to provide a comprehensive comparison between NE-FGR, MQCL, Ehrenfest, and FSSH by evaluating the corresponding predictions of the donor state population relaxation dynamics within the GOA benchmark charge-transfer model.⁵² The GOA model has been used extensively in recent years for benchmarking numerous approximate methods for describing quantum dynamics.^{1,10,46,48–50,54,61,65,67–85} However, a similar benchmarking has not been previously reported to the best of our knowledge.

The overall Hamiltonian of the GOA model is given by⁵² 

Here, $σ ^ x =|DA|+|AD|$ and $σ ^ z =|DD|−|AA|$⁠, where $|D$ and $|A$ are the (R-independent) donor and acceptor electronic states; $Γ D A = Γ A D$ is the electronic coupling coefficient (constant within the Condon approximation); $−ℏ ω D A$ is the donor-to-acceptor reaction energy; y, P_y, and $Ω$ are the primary mode mass-weighted coordinate, momentum, and frequency, respectively; ${ x i , p i , ω i }$ are the mass-weighted coordinates, momenta, and frequencies of the N − 1 secondary bath modes, respectively; {c_i} are the coupling coefficients between the secondary modes and the primary mode; $2 y 0$ is the shift in equilibrium geometry of the primary mode between the donor and acceptor states such that the reorganization energy is given by

The spectral density of the secondary modes is given by

Here, $η$ is the friction coefficient, which determines the coupling strength between the primary mode and secondary modes, and $ω c$ is the cutoff frequency.

The initial state is obtained by shifting the donor thermal equilibrium state by s along the primary mode coordinate

where Z_D is the donor-state partition function and k_BT is the temperature with the Boltzmann constant k_B.

The fact that the primary and secondary modes are harmonic and bilinearly coupled implies an equivalent form of the GOA Hamiltonian in terms of the corresponding normal modes. The latter are obtained by diagonalizing the Hessian matrix,

such that $T T DT= d i a g ( ω ̃ 1 2 , … , ω ̃ N 2 )$⁠, with the normal mode frequencies $ω ̃ j$ and T being the unitary diagonalizing transformation matrix. It should be noted that the normal modes are obtained via numerical diagonalization of the Hessian matrix rather than by using the asymptotic effective spectral density, which is only known in closed form in the limit $Ω<< ω c$⁠.⁸⁶ 

The donor and acceptor nuclear Hamiltonians in terms of normal modes are given by

Here, ${ R j , P j , ω ̃ j }$ are the normal mode coordinates, momenta, and frequencies, respectively, and ${ R j e q }$ are the donor-to-acceptor shifts in equilibrium geometry along the normal mode coordinates,

The same set of normal mode coordinates and momenta is used for all methods in this study.

The initial state in terms of normal modes is given by

where out-of-equilibrium shifts S_j are given by

The fully quantum-mechanical correlation function required in the NE-FGR approach, Eq. (11), can be derived for the GOA model in closed form⁴⁹ 

It should be noted that the LSC approximation for the NE-FGR, Eq. (12), has been recently shown to coincide with the fully quantum-mechanical result, Eq. (43), for the GOA model.⁴⁹ 

MQCL, Ehrenfest, and FSSH calculations on the GOA model were performed in the normal mode representation so that both the primary and the secondary nuclear modes are treated as classical-like. It should be noted that the MQCL equation is known to reproduce the exact fully quantum-mechanical dynamics when the quantum system is bilinearly coupled to a harmonic bath, such as the GOA model considered here.⁸⁷ Although MQCL is formally exact for this system, the TBSH algorithm invokes the momentum jump approximation.^37,54 However, provided, a sufficiently large number of trajectories are used, this approximation has been observed to have almost negligible effect on the density matrix evolution.^38,54,87

Dimensionless variables were used throughout such that energy and time are given in units of $ℏ ω c$ and $ω c − 1$⁠, respectively. The primary mode was chosen to have a frequency of $Ω∕ ω c =0.5$ and an equilibrium displacement of $2 y 0 =2$⁠, resulting in a reorganization energy of $E r ∕ ω c =0.5$⁠. Two values of the donor-to-acceptor reaction energy were considered: $ω D A ∕ ω c ∈ { 0,2 }$⁠, corresponding to the normal (⁠$ℏ ω D A < E r$⁠) and inverted (⁠$ℏ ω D A > E r$⁠) regimes, respectively (see Fig. 1). Three different values of electronic coupling were examined, $Γ D A ∕ℏ ω c ∈ { 0.05 , 0.1 , 0.2 }$⁠, ranging from weak to strong coupling, where we classified the strength retrospectively by its ability to cause significant deviations between the FGR approach and the MQCL reference method (see below). Three different temperatures, $k B T∕ℏ ω c ∈ { 0.2 , 0.5 , 1 }$ (low to high), and four different nonequilibrium initial states, $s ω c ∕ ℏ ∈ { − 1,1,3,5 }$ [see Eq. (37)], were considered. For an analysis of NE-FGR under initially equilibrated conditions, we refer the reader to Ref. 49. The calculations reported below were performed with 200 secondary bath modes sampled from an Ohmic spectral density with $η∕ℏ=1.0$ and a maximum frequency of $15 ω c$⁠. For the sake of clarity, we set $ℏ=1$ and $ω c =1$ in the following.

The initial partially Wignerized density operator $ρ ^ W ( R , P , t = 0 )$ is given by

$ρ b , W e q ( R , P )$ is given as a product of N Wigner-transformed harmonic oscillator density operators, Eq. (41), where each DOF is shifted by S_j given by Eq. (42).

In MQCL simulations, the time-dependent population of the donor state was calculated using the TBSH algorithm, Eq. (22), with $A W α α ′ ( R ) =α;R|DD| α ′ ;R$⁠. All TBSH results presented in this work were obtained by averaging over 5 $×1 0 8$ trajectories. Error bars were generated by performing five independent stochastic TBSH simulations and calculating the standard deviations of these groups from the combined average. Additionally, we verified that converged results can be obtained with a time step of $Δt=0.01$⁠. It should be noted that techniques for improving convergence, such as transition filtering or observable cutting,^88–92 were not employed. In the following, the term MQCL refers synonymously to this implementation of the TBSH algorithm.