Nonequilibrium dynamics of photoinduced forward and backward electron transfer reactions

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.¹¹ 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.¹⁴ 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.²⁴ 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.

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∈R$⁠) 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.³¹ 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³¹ 

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¹¹ 

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₀ > k₁ > ⋯ > k_N.

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 jth group and j ∈ {0, 1, …, N}. As a side note, in the later discussion, we will assume that the correlation function C₀(t), corresponding to the reaction coordinate X₀, 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,³⁰ 

$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.

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$⁠,³⁸ 

where $CjX(τ)$ 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

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.

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

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.³⁰ 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).⁴¹ 

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 $ρX(x→,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→$ is short for the (N + 1) reaction coordinates, {x₀, x₁, …, 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

It turns out that under some general conditions, the off-diagonal terms of the density operator $ρX(x→,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₀. 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 $ρX(X→,t)$ [Eq. (S37) in the supplementary material],

where the superscript N means that the equations are solved for the N coordinates, x₁, …, x_N. Since $ρg(x→,t)$ can be given by the following identity:

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

The reaction kernel $Kf(x→(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→(N))$⁠, which determines the flow of population from |ct⟩ back to |e⟩.⁴² $Kfr(x→(N))$ has the following expression:

Following the same procedure, we get the backward reaction kernel $Kb(x→(N))$ and the reverse backward reaction kernel $Kbr(x→(N))$⁠,

where ΔU_b = −U_ct. $Kb(x→(N))$ determines the distribution of the backward reaction rate, |ct⟩ → |g⟩, along the N reaction coordinates, while $Kbr(x→(N))$ determines the reverse flow from |g⟩ back to |ct⟩.

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.²³ In the case where the rates of the FET and BET are both much smaller than those of the local motions, $ρe(N)(x→(N),t)$ and $ρc(N)(x→(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

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

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

$Kbr(x→)$ [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.⁴³ 

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

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₀ can represent the vibrational relaxations. In this case, the reorganization energy along x₀, λ₀, 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

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₀ and x₁, in which x₀ describes the fast intramolecular vibrational motions. The correlation function of the environment, C₁(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₁. The parameters of the ET model, ΔG_f, ΔG_b, V_f, V_b, and λ₀ and λ₁, 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 $ρe(x→,t)$ and $ρc(x→,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)].²⁴ 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.

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.

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,

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 nth 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.

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→)$ and $Lc(x→)$ are defined as

and $Lg(x→)=Le(x→)$⁠. $ρa,n(x→,t)$⁠, a = g, e, ct, means the population of the nth vibrational level of the electronic state |a⟩. The reaction kernels $Kfn,m(x→)$ and $Kbm.n(x→)$⁠, 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

$λv$ is the reorganization energy contributed by the vibrational mode $ωv$ and is included in the inner reorganization energy, λ₀. The |0, c⟩ → |0, g⟩ Franck–Condon factor $⟨0,g|0,c⟩=exp(−12S)$ has been included into the coupling constants V_f and V_b.

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