The biological functions of photoenzymes are often triggered by photoinduced electron transfer (ET) reactions. An ultrafast backward ET (BET) reaction follows the initial photoinduced forward ET (FET), which dissipates the energy of absorbed photons and terminates the biological function in vain. Based upon our previous works, we reasoned that the dynamics of the BET is coupled with that of the FET and other local motions. In this work, the dynamics of the FET and BET is modeled as the master equation of the reduced density operator of a three-state system coupled with a classical harmonic reservoir. The coupling of the FET and BET is reflected in the time-evolution of the charge-transfer state’s population, which is generated by a source, the reaction flux for the FET, and annihilated by a sink, the reaction flux for the BET. Surprisingly, numerical simulations show that when the BET is in the Marcus normal region, the BET can be accelerated by nonequilibrium local motions and becomes faster than what is predicted by the Marcus theory. The experimental confirmation of this novel dynamics would provide qualitative evidence for nonequilibrium effects on ultrafast ET dynamics. Additionally, the effects of quantum vibrational modes on the dynamics are discussed. This work can help understand the dynamical interactions between the chain of ultrafast reactions and the complex local environmental motions, revealing the physical nature underlying biological functions.

## I. INTRODUCTION

The electron transfer (ET) reaction is one of the fundamental reactions in chemistry and a major building block for biological functions.^{1–3} Within this class of ET reactions, the photoinduced forward ET (FET) reaction, the transfer of an electron triggered by the electron donor or acceptor absorbing a photon, has been shown to play a critical role in the repair of DNA photo-damage by photolyase^{4–8} and other biological processes.^{9,10} These ET reactions are known to be *nonadiabatic* because the coupling between the electron donor and the acceptor is weak enough such that the states participating in the reaction can be approximated as diabatic states.^{11} For ET reactions in condensed matter, the reaction is intrinsically coupled with motions of surrounding molecules. Depending on the relative time scales of these processes, the physical picture of an ET reaction can be vastly different.^{12,13} When the local motions are much faster than the ET reaction such that the reactant and product states can be *instantly* equilibrated, the reaction dynamics is finely described by the celebrated Marcus theory.^{14} On the other end of the spectrum, known as the *solvent-controlled regime*, where the reaction rate is large, the determining factor of the overall ET dynamics becomes the local motions that bring the reactants to the reaction region.^{15–17}

However, in many biological systems, ultrafast photoinduced ET reactions, with reaction time scales in the femtosecond (fs) to picosecond (ps) range, are coupled to a complex local environment in which the motions are characterized by multiple time scales.^{18–21} In these systems, there is no separation of time scales between the ET reaction and local motions. By modeling the environmental motions as a diffusive process along a reaction coordinate while the ET reaction rate being dependent on the reaction coordinate,^{22,23} it is shown that the coupling between ET reactions and local motions can generate a variety of reaction dynamics, whose behaviors are not anticipated by the Marcus theory or solvent-controlled models.^{24} It is further elucidated that because of the slow components of local motions, the product state of the ET reaction is not fully equilibrated during the reaction. As a result, the driving force of the reaction Δ*G* and the reorganization energy contributed by the environment can be significantly modified.^{25,26}

The photoinduced forward ET (FET) reaction is often accompanied by a *backward* electron transfer (BET) from the charge-transfer (CT) state to the ground state of the donor and acceptor.^{27–29} The BET is often a futile step that causes dissipation of energy of the absorbed photon without triggering downstream reactions in the biological machinery.^{4,7} Since the population of the CT state is prepared by the FET and is not equilibrated, it is reasonable to infer that unless the BET happens much slower than other dynamical processes in the system, the dynamics of the BET is coupled to the FET along with local environmental relaxations. To establish a physical picture for the BET, in this work, a model with a system of three states is proposed, which include the ground state of the electron donor and acceptor (|*g*⟩), the photoexcited state (|*e*⟩), and the charge-transfer state (|*ct*⟩). The environment, as well as the low-frequency vibrational modes inside the donor and acceptor, is modeled as a reservoir made of classical harmonic oscillators. The system is assumed to be linearly coupled to the reservoir. Similar approaches have been applied to model ET reactions by different authors using path integrals^{30–32} and master equations.^{33–35} The equations of motion for the reduced density operator of the system are developed. Within the context of ultrafast short-range ET reactions, the physical picture underlying this model is comprehensively discussed.

## II. DEVELOPMENT OF THEORY

### A. The quantum master equation

#### 1. The model Hamiltonian

We assume that the ET system consists of three electronic states, |*g*⟩, |*e*⟩, and |*ct*⟩. It is linearly coupled with a classical harmonic reservoir. The Hamiltonian of the combined system $H0$ reads

$HS$ is the system’s Hamiltonian,

where the ground state energy is shifted to 0 and *U*_{e} and *U*_{ct} are the relative energies of the states |*e*⟩ and |*ct*⟩, respectively. The couplings between the diabatic states, *V*_{f} and *V*_{b}, are real values ($Vf,Vb\u2208R$) and are assumed to be small compared with other energy scales such as |*U*_{e}| and |*U*_{ct}| such that the dynamics can be solved perturbatively with respect to the couplings.^{31} We assume that at *t* = 0, the wavepacket of the system is prepared at the excited state. The process of photo-excitation is assumed to happen instantly and the coupling between |*g*⟩ and |*e*⟩ is assumed to be 0.

The Hamiltonian of the unperturbed reservoir is given by

The interaction between the system and the reservoir is given by^{31}

where *c*_{α} are the couplings between the states and each harmonic mode in the reservoir. It has been shown that the influences of the reservoir on the system’s dynamics can be quantified by the spectral density of the reservoir, *J*(*ω*), defined as^{11}

In the classical limit, the spectral density *J*(*ω*) is related to the correlation function of the unperturbed reservoir’s polarization, $C(t)$, through

We define the real part and imaginary part of $C(t)$ as

in which *C*′(*t*) is proportional to the derivative of *C*(*t*),

In many complex systems, there exist motions with multiple time scales. The correlation function *C*(*t*) of these systems is often modeled as a linear combination of exponential functions,^{36,37}

which gives the spectral density *J*(*ω*) as a linear combination of Lorentzian functions,

We assume that *k*_{j} are ordered such that *k*_{0} > *k*_{1} > ⋯ > *k*_{N}.

#### 2. The reaction coordinates

Assume that we are able to separate the harmonic oscillators of the reservoir into (*N* + 1) groups such that each group of harmonic oscillators has a spectral density of the Lorentzian form. That is, we define

and

We then introduce the reaction coordinates *X*_{j}, *j* = 0, 1, …, *N*,

where the index *α* runs over all the harmonic oscillators in the *j*th group and *j* ∈ {0, 1, …, *N*}. As a side note, in the later discussion, we will assume that the correlation function *C*_{0}(*t*), corresponding to the reaction coordinate *X*_{0}, characterizes the classical intramolecular vibrations inside the donor and acceptor pair, which decays much faster than other *C*_{i}(*t*), *i* = 1, 2, …, *N*, characterizing other motions of the environment.

Then, the Hamiltonian $H0$ [Eq. (1)] can be canonically transformed into the Hamiltonian $H$ with reaction coordinates,^{30}

$HX$ is the *new system*’s Hamiltonian which includes the reaction coordinates,

where $HS$ is defined in Eq. (2) and *P*_{j} is the conjugate momentum for the reaction coordinate *X*_{j}. $HjI$ is the coupling term between the electronic states and the reaction coordinate *X*_{j}, given by

$HR$ is the new reservoir after the canonical transformation,

*V*_{X} gives the coupling between the harmonic modes of the new reservoir with the reaction coordinates,

After the canonical transformation, the ET system, characterized by the Hamiltonian $HS$, is only coupled to the (*N* + 1) reaction coordinates, while each reaction coordinate is coupled to an isolated (transformed) harmonic bath.

#### 3. The master equation

If we assume that the coupling between $HX$ and $HR$ is weak, we can perform a perturbative expansion with respect to the coupling coefficients in *V*_{X} and average over the reservoir’s degrees of freedom. We then obtain the quantum master equation for the reduced density operator *ρ*_{X}(*t*), which describes the evolution of the system $HX$,^{38}

where $CjX(\tau )$ is the correlation function of the part of harmonic oscillators in the new reservoir that couples to *X*_{j} and *U*_{X}(*t*) is the system’s evolution operator, given by

### B. The equations of motion (EOMs) for the populations

To reach the working model of ultrafast ET reactions, we will perform a series of approximations on the master equation [Eq. (19)]. We only outline the results after the approximations are performed. For mathematical details, please refer to the supplementary material.

#### 1. Quantum-classical approximation

Since the reaction coordinates *X*_{j} are linear combinations of classical harmonic oscillators, they are classical variables. The operators $HX$ and *ρ*_{X}(*t*) are both dependent on the classical degrees of freedom, *P*_{j} and *X*_{j}. We, therefore, can make the quantum-classical approximation on the quantum master equation [Eq. (19)].^{39,40}

#### 2. The Smoluchowski–Kramers approximation

If the original correlation function *C*_{j}(*t*) is an exponential function, the transformed correlation function $CjX(t)$ [the real part of $CjX(t)$] becomes a delta function.^{30} In other words, the reaction coordinates *X*_{j}, *j* = 0, 1, …, *N*, are overdamped such that we can perform the Smoluchowski–Kramers approximation, which gives the Fokker–Planck-type equation (shown after the next assumption).^{41}

#### 3. Assume *Q*_{e} = 0

In practice, the excited state is closely aligned with the ground state, which means

As a first-order approximation, we assume *Q*_{e} = 0. Then, the Fokker–Planck equation of the density operator $\rho X(x\u2192,t)$ reads

where [*A*,*B*]_{+} = *AB* + *BA* is the anti-commutator for the operators *A* and *B* and $HS$ is given by Eq. (2). To simplify the notations, we defined unitless reaction coordinates *x*_{j} as

such that the minimum of the states |*g*⟩ and |*e*⟩ along *x*_{j} is at *x*_{j} = 0, while that of |*ct*⟩ is at *x*_{j} = 1. The vector $x\u2192$ is short for the (*N* + 1) reaction coordinates, {*x*_{0}, *x*_{1}, …, *x*_{N}}. We also defined the Liouville operator along *x*_{j},

We defined the reorganization energy *λ* as

*λ*_{j} can be understood as the reorganization energy for the ET reaction along the reaction coordinate *X*_{j}. The diffusion coefficient along the coordinate *x*_{j}, *D*_{j}, is given by

The operator *K* is given by

#### 4. Fast motion along *x*_{0}

It turns out that under some general conditions, the off-diagonal terms of the density operator $\rho X(x\u2192,t)$ (coherences) have much shorter lifetimes compared to the time scale of the system’s evolution. Hence, they can be solved as expressions of the diagonal terms (see Sec. I C in the supplementary material). Then, since we are interested in the coupling of the system’s evolution with environmental relaxations when there is no separation of time scales, we assume that there exists at least one reaction coordinate that has much faster relaxation compared to the evolution of populations. We let it be *x*_{0}. With this assumption, we finally obtain a closed form of equations of motion for the populations of the three states, which are the diagonal terms of $\rho X(X\u2192,t)$ [Eq. (S37) in the supplementary material],

where the superscript *N* means that the equations are solved for the *N* coordinates, *x*_{1}, …, *x*_{N}. Since $\rho g(x\u2192,t)$ can be given by the following identity:

the EOM of $\rho g(x\u2192,t)$ does not offer new information. The dynamics of the photoinduced ET system is governed by these equations.

The reaction kernel $Kf(x\u2192(N))$ determines the distribution of the forward reaction rate along the *N* dimensions and is given by

where Δ*U*_{f} = *U*_{ct} − *U*_{e}. Similarly, we can compute the kernel of the *reverse* forward reaction, $Kfr(x\u2192(N))$, which determines the flow of population from |*ct*⟩ back to |*e*⟩.^{42} $Kfr(x\u2192(N))$ has the following expression:

Following the same procedure, we get the backward reaction kernel $Kb(x\u2192(N))$ and the *reverse* backward reaction kernel $Kbr(x\u2192(N))$,

where Δ*U*_{b} = −*U*_{ct}. $Kb(x\u2192(N))$ determines the distribution of the backward reaction rate, |*ct*⟩ → |*g*⟩, along the *N* reaction coordinates, while $Kbr(x\u2192(N))$ determines the reverse flow from |*g*⟩ back to |*ct*⟩.

### C. Reduction to the Marcus theory

It is easy to see that when *N* = 1, except for the reverse reaction kernel, the EOM for *ρ*_{e}(*x*, *t*) in Eq. (28) is similar to the ET model developed by Sumi and Marcus because their model was developed by assuming the Debye relaxation for the solvent coordinate.^{23} In the case where the rates of the FET and BET are both much smaller than those of the local motions, $\rho e(N)(x\u2192(N),t)$ and $\rho c(N)(x\u2192(N),t)$ are instantly equilibrated along all dimensions. As a result, the reaction kernels of the FET and BET are reduced to constant rates, given by the Marcus formula,

The Δ*G*_{f} is the forward reaction free energy and is related to Δ*U*_{f} by Δ*G*_{f} = Δ*U*_{f} − *λ*. Similarly, Δ*G*_{b} is the backward reaction free energy, given by Δ*G*_{b} = Δ*U*_{b} + *λ*. Similarly, the reverse reaction kernels are reduced to

### D. Dynamical suppression of reverse reactions

If we compare the reverse reaction kernel $Kar(x\u2192(N))$ with $Ka(x\u2192(N))$, *a* = *f*, *b*, $Kfr(x\u2192)$ is shifted negatively (toward the state |*e*⟩’s equilibrium position) along all dimensions compared to $Kf(x\u2192)$, while $Kbr(x\u2192)$ is shifted positively (toward the state |*ct*⟩’s equilibrium position) compared to $Kb(x\u2192)$. Since the width of each reaction kernel is proportional to $2kBT\lambda 0$, we find that if the following condition is satisfied:

the reverse forward kernel $Kfr(x\u2192)$ [Eq. (31)] is dynamically suppressed. Similarly, when

$Kbr(x\u2192)$ [Eq. (32)] is dynamically suppressed.

Under these conditions, we can further simplify Eq. (28) and get

where

and

However, these conditions of dynamical suppression do not always hold. In fact, reverse reactions have been observed in biological ET systems when the magnitude of the reaction free energy, |Δ*G*_{f}| or |Δ*G*_{b}|, is small.^{43}

In Sec. III, we discuss some qualitative results of this model using the ultrafast short-range ET reaction in proteins as an example.

## III. DISCUSSIONS

In some biological systems, the process of vibrational relaxations has been observed in experiments and has a time scale from sub-ps to a few ps.^{28,29} If the ET rate is smaller than the rate of intramolecular vibrational relaxations inside the donor and acceptor pair, *x*_{0} can represent the vibrational relaxations. In this case, the reorganization energy along *x*_{0}, *λ*_{0}, becomes the inner reorganization energy contributed by intramolecular vibrations of the donor and acceptor. In experiments, the reservoir correlation function *C*(*t*) is usually fitted by a sum of several exponential functions, which is also the assumed expression of *C*(*t*) [Eq. (11)] in this work.^{36,37}

### A. Solvation-enhanced backward electron transfer

An unusual phenomenon is observed when we use the model [Eq. (37)] to perform simulations of the coupled dynamics of the FET and BET. We find that when BET is in the Marcus normal region (|Δ*G*_{b}| < *λ*), the rate of the back electron transfer can become faster than the rate predicted by the Marcus formula.

We illustrate this phenomenon by the following example. We assume that the reservoir is described by two reaction coordinates (*N* = 1), *x*_{0} and *x*_{1}, in which *x*_{0} describes the fast intramolecular vibrational motions. The correlation function of the environment, *C*_{1}(*t*) [Eq. (9)], is an exponential function,

where *τ*_{D} is the solvation time constant for the single reaction coordinate describing the motions of the environment, *x* ≡ *x*_{1}. The parameters of the ET model, Δ*G*_{f}, Δ*G*_{b}, *V*_{f}, *V*_{b}, and *λ*_{0} and *λ*_{1}, are chosen such that both *K*_{f}(*x*) and *K*_{b}(*x*) [Eq. (39)] are in the Marcus normal region (Fig. 1). Given the values of these parameters (see the caption of Fig. 1), it is easy to check that both $Kfr(x)$ and $Kbr(x)$ are dynamically suppressed. The dynamics of the FET and BET can be represented by the time-evolution of the populations of the excited state and the CT state, *ρ*_{e}(*t*) and *ρ*_{c}(*t*), which are generally given by

where $\rho e(x\u2192,t)$ and $\rho c(x\u2192,t)$ are the solutions of Eq. (37).

Next, we simulate the time-evolution of *ρ*_{e}(*t*) and *ρ*_{c}(*t*) with different values of *τ*_{D}, while fixing other parameters (Fig. 2). The dynamics of the FET is monotonically slowed down as *τ*_{D} increases, with the rate predicted by the Marcus formula being the upper limit of the forward ET rate [Fig. 2(a)].^{24} However, the dynamics of the BET is more involved. When the solvation is much faster than the FET and BET, it is expected that the dynamics of the BET is close to the Marcus picture [red line vs black line in Fig. 2(b)]. This is because when the solvation is fast, the distribution of the CT state, *ρ*_{c}(*x*, *t*), always stays at equilibrium [Fig. 3(a)]. Interestingly, when the solvation gets slower and becomes comparable to the ET reaction rates, the BET is accelerated [orange line and green line vs red line in Fig. 2(b)]. This is because the maximum of the backward kernel *K*_{b}(*x*) is in between the free energy minima of the excited and CT states (orange line in Fig. 1). *ρ*_{c}(*x*, *t*), once generated by the FET, has to *gradually cross* the maximal region of *K*_{b}(*x*) before reaching its equilibrium centered at *x* = 1 [Figs. 3(b) and 3(c)], as compared to Fig. 3(a) in which *ρ*_{c}(*x*, *t*) *instantly crosses* the maximal region of *K*_{b}(*x*) and equilibrates. However, when the solvation becomes much slower than the ET reaction rates, the distribution *ρ*_{c}(*x*, *t*) would be depleted by the BET without diffusing significantly along the solvent coordinate [Fig. 3(d)]. This results in slower dynamics of the BET [purple line in Fig. 2(b)]. Moreover, during the evolution, *ρ*_{c}(*x*, *t*) can deviate significantly from the Gaussian distribution, which is a strong indicator of nonexponential dynamics of the BET.

Based on the analysis, we find that the enhanced rate of the BET is a result of the nonequilibrium motion along the solvent coordinate, *x*. Hence, we name this phenomenon the *solvation-enhanced backward electron transfer*. It is easy to see that the solvation-enhanced BET can exist regardless of whether the FET is in the Marcus normal or inverted region. In general, this type of reaction dynamics can be observed in an environment in which local motions with larger and smaller rates, compared to the reaction rates, coexist.

### B. The effects of quantum vibrational modes

In several ET systems examined experimentally, the forward ET reaction is triggered by the donor or acceptor absorbing light that is in the UV or visible range. In these systems, the FET is usually in the normal region defined by the Marcus theory in which |Δ*G*_{f}| < *λ*, while the BET is in the inverted region.^{4,43} It is well known thatET reactions in the inverted region can be significantly assisted by high-frequency vibrational modes.^{44–46} Indeed, vibrationally excited states have been observed in some photoinduced ET systems.^{28,29} Hence, in this section, quantum vibrational modes are included into the photoinduced ET model.

#### 1. Single effective vibrational mode

We assume that there exists an *effective* harmonic vibrational mode inside the donor and acceptor, which participates in the ET reactions. The system’s Hamiltonian $HS$ in Eq. (2) is, therefore, modified,

where $Hva$, *a* = *g*, *c*, is the free vibrational Hamiltonian for each state and |*e*⟩ has the same vibrational mode as |*g*⟩ does. $Hvg$ is given by

where *b*^{†} and *b* are the creation and annihilation operators for the vibrational mode. The zero-point energy for the vibrational mode is neglected since it is the same for all the states. $Hvc$ has the following expression:

where *b*_{c}^{†} is related to *b* through

where $D$ is the displacement operator,

#### 2. Intramolecular vibrational redistribution

Since there exist other intramolecular vibrational modes, vibrationally excited states can quickly decay due to anharmonic coupling with other vibrational modes. If the dominant channel of the decay is due to cubic anharmonic coupling and through intramolecular vibrational redistributions (IVR), the decay rate for a particle in the *n*th vibrational level, $knv$, is approximately given by

where $k1v$ is the decay rate of |1⟩ → |0⟩.^{47,48} We assume that this decay rate applies to the vibrational states of all electronic states.

#### 3. The equations of motion

Similarly, by using the approach of quantum master equations and performing the quantum-classical approximation, we obtain the equations of motion for the diagonal density matrix elements,

where $Le(x\u2192)$ and $Lc(x\u2192)$ are defined as

and $Lg(x\u2192)=Le(x\u2192)$. $\rho a,n(x\u2192,t)$, *a* = *g*, *e*, *ct*, means the population of the *n*th vibrational level of the electronic state |*a*⟩. The reaction kernels $Kfn,m(x\u2192)$ and $Kbm.n(x\u2192)$, which are dependent on the vibrational quantum numbers, are given by

It is relatively straightforward to write out the expressions of vibrationally modified reverse reaction kernels, which are not provided here. The coupling factor ⟨*n*, *g*|*m*, *c*⟩ is given by

where

$\lambda v$ is the reorganization energy contributed by the vibrational mode $\omega v$ and is included in the inner reorganization energy, *λ*_{0}. The |0, *c*⟩ → |0, *g*⟩ Franck–Condon factor $\u27e80,g|0,c\u27e9=exp(\u221212S)$ has been included into the coupling constants *V*_{f} and *V*_{b}.

#### 4. Qualitative analysis of the effects of vibrational modes

Before closing the section, we want to discuss qualitatively the effects of high-frequency vibrational modes on ET dynamics. The quantum vibrational effects on ET reactions have been discussed extensively in the literature.^{13,45,46,49} However, the existence of nonequilibrium dynamics introduces more complexity. We assume that the FET is in the Marcus normal region, whose dynamics is generally not affected by quantum vibrational modes, while the BET is in the Marcus inverted region, on which the vibrational effects cannot be neglected. We assume that the nonequilibrium solvent motions can be projected to a single coordinate (*N* = 1). The forward and backward reaction kernels, *K*_{f}(*x*) and *K*_{b}(*x*), are shown in Fig. 4 (Solid lines), where *K*_{b}(*x*) (in orange) is the sum of all the backward reaction kernels, $Kb0,n$. *K*_{b}(*x*) represents the cumulative distribution of the backward ET reaction rate. The dashed lines in Fig. 4 represent the backward reaction kernels with different vibrational quantum numbers, $Kb0,n(x)$ [Eq. (50)], multiplied by the factor

where only the vibrational “hot” states of the ground electronic state (|*g*⟩) are considered. Since the BET is in the Marcus inverted region, the center of $Kb0,0(x)$ (corresponding to the 0 → 0 transition for the vibrational levels) is on the right of *x* = 1. The center of $Kb0,n(x)$ moves toward *x* = 0 as *n* increases until the quantity Δ*G*_{b} − *λ* + 2*λ*_{0} + *nℏω*_{v} in Eq. (50) becomes positive. Since the relaxation along *x* happens on the similar time scale as the reaction does, the population of the CT state, *ρ*_{c}(*x*, *t*), is generated close to the equilibrium position of the excited state *x* = 0 and relaxes toward the equilibrium of the CT state *x* = 1, similar to the dynamics displayed in Fig. 3. Clearly, depending on the relative rate of the reaction and the motion along *x*, we expect to see a variety of BET dynamics that display nonexponential behaviors with different average reaction rates.

Finally, we want to add that although they have not been analyzed in detail, the effects of relaxation dynamics characterized by multiple time scales (*N* > 1) are not simple generalizations of the case of *N* = 1. We think that much of the complexity in ultrafast ET dynamics is a result of these effects. This is a subject we plan to study in the future work.

## IV. CONCLUSIONS

In order to understand the photoinduced ET reactions, a model of a three-state system coupled with a harmonic reservoir is developed. The spectral density of the reservoir is assumed to be a sum of multiple Lorentzian functions such that reaction coordinates with overdamped motions can be introduced. The effects of quantum vibrational modes are discussed. The connections with previous theoretical models are established. There are a few free parameters in this model, including the forward and backward reaction free energies, the reorganization energies, and the correlation function of environmental relaxations, many of which are available experimentally.^{24} We thus believe that this model can serve as an analytical tool to understand experimental data of photoinduced ET reactions in complex environments. Additionally, we predicted that when the BET is in the Marcus normal region and some relaxation rates of the environment are comparable to the ET reaction rates, the reaction rate of the BET can be enhanced by the solvation process and becomes faster than the rate given by the Marcus formula. The existence of this novel reaction dynamics, if confirmed experimentally, would demonstrate the effects of nonequilibrium local motions on the (ultrafast) ET dynamics qualitatively.

In this work, we show that the backward ET reaction is coupled with the forward ET reaction as well as the nonequilibrium vibrational and environmental relaxations. Biological functions often involve a series of ultrafast reactions, whose dynamics are correlated with each other as well as motions in the complex biological environments.^{4,5} An understanding of these intricate interplays is key to unveiling the fundamental mechanisms of biological functions.

There are a few limitations in this work. First, the dynamics in the sub-ps time scale is not well described. Within this time scale, the ET dynamics can be correlated with the dynamics of photo-excitation and the dephasing of vibrational excited states.^{48} If the ET time scale is comparable to the time scale of nuclear relaxations inside the donor and acceptor, the reverse reactions should also be considered.^{43,50} Second, in general, the spectral density can be more complicated than a linear combination of Lorentzian functions, which, as a result, will lead to a master equation with memory kernels.^{51} Furthermore, it has been discussed that the assumption of a harmonic reservoir may not be justified for relaxations on the ps time scale.^{52–54}

## SUPPLEMENTARY MATERIAL

See the supplementary material for the detailed derivation of Eq. (28).

## ACKNOWLEDGMENTS

This work was supported in part by the National Science Foundation (Grant No. CHE1412033) and the National Institutes of Health (Grant No. GM118332). We also acknowledge the support of the National Natural Science Foundation of China for the collaborative effort with a visit of Y.L. and a sabbatical stay of D.Z. in Shanghai Jiao Tong University to finally finish this work.