Semiclassical Theory of Multistage Nonequilibrium Electron Transfer in Macromolecular Compounds in Polar Media with Several Relaxation Timescales

Many specific features of ultrafast electron transfer (ET) reactions in macromolecular compounds can be attributed to nonequilibrium configurations of intramolecular vibrational degrees of freedom and the environment. In photoinduced ET, nonequilibrium nuclear configurations are often produced at the stage of optical excitation, but they can also be the result of electron tunneling itself, i.e., fast redistribution of charges within the macromolecule. A consistent theoretical description of ultrafast ET requires an explicit consideration of the nuclear subsystem, including its evolution between electron jumps. In this paper, the effect of the multi-timescale nuclear reorganization on ET transitions in macromolecular compounds is studied, and a general theory of ultrafast ET in non-Debye polar environments with a multi-component relaxation function is developed. Particular attention is paid to designing the multidimensional space of nonequilibrium nuclear configurations, as well as constructing the diabatic free energy surfaces for the ET states. The reorganization energies of individual ET transitions, the equilibrium energies of ET states, and the relaxation properties of the environment are used as input data for the theory. The effect of the system-environment interaction on the ET kinetics is discussed, and mechanisms for enhancing the efficiency of charge separation in macromolecular compounds are analyzed.


Introduction
Intramolecular electron transfer (ET) is an essential component of many biological processes including photosynthesis, cellular respiration, activation of sensory proteins, and others [1][2][3][4][5][6]. In ET reactions, quantum tunneling of an electron between the redox sites within a macromolecule is largely controlled by the environment that creates the necessary preconditions for an electron jump. Both fast fluctuations of solvent polarization and relatively slow and large-scale reorganization of the environment around the redox sites are important in these processes [7][8][9][10].
The role of the environment is even more pronounced in ultrafast photoinduced reactions, where the evolution of the Franck-Condon state is strongly coupled to the medium [11][12][13][14][15]. These reactions often exhibit unusual behavior that does not correspond to the predictions of classical ET theories. For example, photoinduced ET in donor-acceptor compounds may proceed faster than the characteristic timescale of solvent relaxation [16], which indicates the inapplicability of the Zusman model of diffusion along the ET coordinate [17]. Ultrafast reactions are also sensitive to the properties of the pumping pulse [18,19], and often manifest oscillatory kinetics due to vibrational coherence effects [20][21][22].
One of the promising areas of research in the field of photovoltaics and solar energy conversion is the design of macromolecular compounds aimed at efficient photochemical separation of charges [23][24][25]. These compounds commonly involve several redox sites organized in a way that ensures the stabilization of separated charges. This stabilization is achieved by fast and efficient ET from a photosensitizer and primary electron donor to the distant electron-accepting centers, thus preventing loss of energy due to charge recombination. Fast ET transitions in these compounds lead to the formation of nonequilibrium nuclear configurations and strong coupling of the subsequent ET steps to relaxation processes in the environment. To describe these phenomena, in this paper we develop a semiclassical model of multistage intramolecular electron transfer reactions controlled by the medium with a multi-component dielectric relaxation function. Particular attention is paid to the design of the configuration space of the medium, as well as methods for constructing the diabatic free energy surfaces corresponding to the ET states. The study is mainly focused on the role of the polar environment in these processes and does not take into account the ET-active high-frequency (hω k B T) intramolecular vibrational modes. However, an appropriate extension of the theory is straightforward and can be easily introduced to the model. In what follows we assume that the macromolecule is «rigid» in the course of ultrafast ET, so that we can omit reorganization of slow large-amplitude intramolecular nuclear coordinates and consider only fast dielectric response both from the outer environment and internal low-frequency (classical) degrees of freedom.
Many theoretical approaches describe electron transfer in condensed phase in terms of the energy-gap coordinate-the quantity evaluated by mapping the system's nuclear degrees of freedom onto the energy gap ∆E between the ET reactant and product states [26][27][28]. Using the ∆E quantity as the reaction coordinate allows one to reduce the dimensionality of the nuclear configuration space, and ultimately admit the one-dimensional picture of ET with the two intersecting electronic surfaces for the reactant and product states in the diabatic presentation. This approach was suggested by Marcus to describe intermolecular electron transfer in polar liquids [29,30], and later was applied to intramolecular processes, including those initiated by optical excitation of reactants [31][32][33][34][35].
When applied to multistage ET reactions involving several redox centers within a macromolecule, the Marcus approach suggests using several reaction coordinates ∆E k , each controlling ET between the two centers. This extension is straightforward and enables one to consider the multistage reaction as a set of elementary ET steps occurring sequentially and/or in parallel. The key problem of this approach is the non-orthogonality of the ∆E k coordinates [36][37][38]. This non-orthogonality implies some important features of the reaction, for example, the mutual influence of consequent ET steps on each other [39]. However, the multi-∆E k approach leads to apparent technical difficulties associated with the use of a non-orthogonal basis. At present, the method has been implemented only for the three-center model systems of the type A 1 -D-A 2 or D-A 1 -A 2 , where D is the electron donor, and A 1 /A 2 are the acceptors [38,40].
An alternative approach was proposed by Tang and Norris [41], who introduced two independent solvent coordinates (x, y) to describe ET in a model system involving the initial (reactant) |1 , intermediate |2 , and final (product) |3 states. This three-component molecular system reproduces the structure of the photosynthetic reaction center, where primary photochemical charge separation proceeds as ultrafast ET from a bacteriochlorophyll special pair dimer P * (donor) through a bridging accessory bacteriochlorophyll monomer B L (primary acceptor) to a distant bacteriopheophytin H L (secondary acceptor). The use of the independent solvent coordinates in this model implies certain assumptions about the medium, namely, its linear response to charge redistribution between the molecular centers. Now, this model is a common framework for the description of ET in molecular triads [42][43][44][45].
Recently a more general theory of solvent-controlled ET in the Debye polar liquids was developed [46], applicable to multistage reactions with an arbitrary number of elementary ET steps. Following Tang and Norris, the configuration space of multistage ET was constructed with the use of independent solvent coordinates Q k . The number K of the Q k coordinates was shown to be related to the number N of the active redox sites according to the equation K = N − 1. It was also shown that the ET energy gaps ∆E k in this model can be calculated by simple linear transformations over the Q k coordinates [46].
In the present paper, the theory developed in ref. [46] is extended to the non-Debye environments, i.e., media characterized by several dielectric relaxation timescales τ i . The aim of this extension is to broaden the scope of the previous theory, and make it applicable to «complex» environments, as one may expect in molecular electronics devices, biomolecules, and other macromolecular compounds [34,47,48]. As a method, we employ «splitting» of the polarization coordinates Q k into the relaxation components, each corresponding to the individual τ i mode. This approach is similar to that used by Zusman [49], but it is applied here to independent polarization coordinates rather than the energy-gap coordinate ∆E k . In the next section, we formulate the general idea of the method and introduce the corresponding mathematical framework. After that, we consider an illustrative example by applying the general theory to a simple three-center model in a medium with the two-component relaxation function. Finally, the proposed general approach is verified by analyzing the two special cases in which the results of the well-known Najbar/Tachiya and Zusman models of ET are reproduced.

General Formulation of the Theory
We consider a macromolecule containing N redox centers and denote |ϕ n (n = 1, · · · , N) the state of the system with the electron localized on the nth redox center. We admit here a general case when ET transitions are possible between any |ϕ n and |ϕ n states. We introduce the free energy of medium reorganization λ (nn ) that quantifies the response of the environment to the shift of the electronic density |ϕ n |ϕ n . Since λ (nn) = 0 and λ (nn ) = λ (n n) , the λ (nn ) values form a square symmetric matrix of the size N with zero diagonal elementŝ λ (12) . . . λ (1N) λ (12) 0 The number of independent parameters here is N(N − 1)/2. Complex dynamics of medium polarization in view of multiple redox centers and ET transitions can be considered by extending the method developed in ref [46]. First, we introduce the set of independent polarization coordinates Q k , where k ranges from 1 to K = N − 1. In the linear-response approximation, the diabatic free energy surfaces (FESs) of the |ϕ n electronic states in the Q k coordinates are written as [46] HereQ (n) k are the coordinates of the FES minimum (unknown model parameters so far), andǦ (n) is the equilibrium free energy of the nth electronic state. The mutual arrangement of the G (n) surfaces in the Q space is fully determined by theλ matrix. To show this, one can use the definition of λ (nn ) and Equation (2) to obtain It is easy to see from this equation that D (nn ) is a vector in the Q space connecting the G (n) and G (n ) FES minima [46]. The distance between theQ (n ) andQ (n) points is thus equal to √ λ (nn ) , that is determined only by the free energy of the medium reorganization for the |ϕ n → |ϕ n ET transition.
Dynamic properties of the medium are commonly introduced to the model using the energy-gap autocorrelation function X(t) = ∆E(0)∆E(t) / ∆E(0)∆E(0) . This function is an experimentally measured quantity and is often approximated by a sum of several exponentials [50,51] The X i (t) terms in Equation (5) are often referred to as the relaxation components of the medium, with x i and τ i being the weight and the relaxation timescale of the ith component. The x i quantities satisfy the normalization condition ∑ x i = 1. It should be noted here that polar solvents generally reveal 2-3 relaxation components with the τ i values differing by several times. In mixtures and structured environments, the range of τ i s can be even larger. An early-time X(t) dependence in liquids is often determined using a combination of ultrafast spectroscopic techniques and computational methods (see, e.g., [52]).
Due to the linearity of the medium response, the X(t) splitting into the relaxation components X(t) = (X 1 (t), X 2 (t), . . . , X R (t)) in Equation (5) is expected to be valid, not only for the energy-gap coordinate ∆E (the Zusman method [49]), but for any polarization coordinate as well. Keeping this in mind, we replace the single Q k value with an Rdimensional vector as follows where q ki is the ith component of Q k corresponding to the ith relaxation mode. The coordinate splitting (6) extends the Q space and produces a D-dimensional configuration space, where D = (N − 1)R.
By construction, this multidimensional space (denoted in what follows as q) is a direct product of the «polarization» and «relaxation» subspaces, Q and R, with the coordinates q ki = (Q k , R i ), where R i (i = 1, . . . , R) are the coordinates of R. Figure 1 illustrates this basic idea. We will show later that the splitting method allows one to get rather simple model equations in spite of the expanded dimensionality of the configuration space. The coordinate splitting (6) is an essentially new element of the theory compared to that presented in ref. [46]. The extended configuration space enables one to describe not only the nonequilibrium polarization of the environment, but also the multi-component relaxation of this polarization to a new equilibrium state. Since the q ki coordinate is directly related to the i-th mode of the environment, the system's relaxation along q ki proceeds with the timescale τ i . This approach allows one to overcome the limitation of the previous theory and develop a more general framework applicable to non-Debye solvents, solvent mixtures, and other environments with complex relaxation functions. Now we explore the properties of the q space. In the q ki coordinates, Equation (2) becomes Using (3) and (8), one can find the following expression for the reorganization free energy This expression is similar to Equation (4), but is written for the D-dimensional vector d (nn ) in the extended configuration space. Calculating the projection of the d (nn ) vector Here λ (nn ) i is the part of λ (nn ) corresponding to the i-th relaxation mode, (10) can be easily verified by direct summation over the Q k coordinates.
Equations (9) and (10)  ki . It follows from Equation (9), thatq (n) ki can be found using the distances between the FESs minima in the q space, d (nn ) = √ λ (nn ) . This expression allows us to formulate a general scheme for theq (n) ki evaluation. As an illustration, we consider a simple three-center model of ET in a medium with two relaxation components. The choice of this model is due to the following reasons: (1) theq (n) ki values in this case can be found analytically, (2) the model to be verified for compatibility with earlier results.

Three-Center Molecular System in a Two-Component Environment
We consider a molecular triad of the type DA 1 A 2 , where D is the electron donor and photosensitizer, and A 1 /A 2 are the electron acceptors. Optical excitation of this compound leads to a population of the locally excited state DA 1 A 2 → D * A 1 A 2 , and triggers a series of electronic transitions, including charge separation to the D + A − 1 A 2 and D + A 1 A − 2 states, charge shift between the A 1 and A 2 units, charge recombination back to the DA 1 A 2 state, and radiationless deactivation of the excited state. The general scheme of the reaction is shown in Figure 2 and involves the following transitions It should be noted here that many molecular compounds have the same structure as the redox centers, and can be considered prototype systems. A common and well-known example is the photosynthetic reaction center (RC), where primary charge separation proceeds as a two-step ET involving a BChl special pair dimer (D), an accessory BChl monomer (A 1 ) and a BPh molecule (A 2 ). The spatial arrangement of the D, A 1 , A 2 units in RCs do not facilitate direct electron transfer between D and A 2 , therefore, the CS 2 and CR 2 transitions in these compounds are effectively suppressed. This fact, however, does not change the arrangement of the electronic FESs in the q space, since positions of the FESs minima are determined solely by the medium reorganization. The molecular triad is assumed to interact with a polar environment with the twocomponent relaxation function X(t). The relaxation parameters of the medium (x 1 , τ 1 , x 2 , τ 2 ) are considered to be known. Our goal here is to construct a general model of ET in this system taking into account both multiple redox centers and the multi-component relaxation of the medium. It is important to note that this is the simplest model, which has not yet been studied within semiclassical theories. On the one hand, the two-stage ET models in the Debye solvents are known [36][37][38][39]44], and, on the other hand, there are general models of single-step ET in liquids with several relaxation timescales [34,49]. The theory presented here combines both these approaches.
According to Equation (7), the configuration space for a three-center (N = 3) molecular system in a two-mode (R = 2) medium should minimally involve (N − 1)R = 4 independent coordinates. It is convenient to represent the 4-dimensional vector q as a 2 × 2 matrix q = q 11 q 12 q 21 q 22 .
As before, the first index here relates to the polarization coordinate Q k , and the second one-to the relaxation component R i . We assume that optical excitation does not produce a significant redistribution of the electronic density, so medium reorganization at the stage of excitation can be neglected. According to Equation (9), the positions of the ground-state and the excited-state FESs minima, in this case, coincide,q (ex) ki =q (gr) ki . To simplify the description we omit the ground state from consideration and restrict our model to electronic transitions with the charge-transfer character only. The following notation will be used for the ET states The equilibrium free energiesǦ (1) ,Ǧ (2) ,Ǧ (3) , and the reorganization free energies λ (12) , λ (13) and λ (23) are considered as known parameters.
Equations (12)- (14) can be easily verified by direct evaluation of the d The G (n) (q) FESs, constructed using Equations (8) and (12)- (15), completely determine the energy barriers for all ET steps. As an example we consider the |ϕ 1 → |ϕ 2 transition, i.e., transfer of an electron from D * to A 1 . The ET is possible only if the reactant and product energies are equal, G (1) (q) = G (2) (q). The solution of this equation is a hyperplane in the 4-dimensional space q 11 λ (12) x 1 + q 12 λ (12) x 2 = λ (12) that describes the intersection of the two parabolic FESs. The section itself is a threedimensional paraboloid, and the minimum value on this surface is the saddle point between the reactant and product FESs (we denote it as q ). The height H 12 of the energy barrier between the equilibrium states then can be calculated using the q value. This parameter is very important for the analysis of thermal (quasi-equilibrium) reactions and included many expressions for the ET rate constants [26,28,33]. At the same time, the reaction pathway in the case of nonequilibrium ET may not be limited to the saddle point, especially in ultrafast processes, where the reaction flux deviates from the saddle point significantly [53]. The effective energy barrier in ultrafast ET may therefore differ from H 12 . These peculiarities of nonequilibrium ET can be taken into account by considering the motion of the wave packets on the corresponding FESs along with ET transitions. We introduce the time-dependent density function ρ n (q, t) for the nth ET state of the system. According to the Zusman method [17], the q ki dynamics can be considered as diffusion, and chemical transformations as reversible quantum transitions localized at the intersection regions of the reactant and product FESs. Applying this method, one obtains the following equation for ρ 1 (q, t) The last two terms on the right-hand side describe the ET transitions |ϕ 1 ↔ |ϕ 2 and |ϕ 1 ↔ |ϕ 3 in the nonadiabatic limit. V 12 and V 13 are the coupling energies between the corresponding electronic states,L Similar kinetic equations can be written for the ρ 2 (q, t) and ρ 3 (q, t) densities, they have the same structure as Equation (17). The resulting set of equations accomplishes the formulation of the model. This set can also be supplemented with the initial conditions, which usually reflect the conditions of the initial state formation. In the case of activated ET from the |ϕ 1 state, the initial distribution on the G (1) FES is often assumed to be thermodynamically equilibrium In ultrafast photoreactions, on the contrary, the initial distribution in the excited state is generally nonequilibrium and depends both on the energetic parameters of the system and spectral characteristics of the pumping pulse (see, e.g., [54]).
The proposed theory can be used as a framework for numerical simulations and analysis of experimental data on ultrafast ET in macromolecules. The theory can be linked to experiments, for example, by calculating the populations of the electronic states and comparing them with the observed ET kinetics. In the present model, the |ϕ n state population is evaluated as the integral of the density function over the entire configuration space Another important quantity of ultrafast ET is the quantum yield of charge separation Y CS , which is commonly identified with the total population of charge-transfer states at the time T when relaxation processes are over [55,56] It is important to note here, that the ρ n (q, t) functions can also be used for simulations of luminescent properties of the system in the course of ET. Recently, specific computational methods have been developed for simulations of the reaction kinetics and transient absorption/emission spectra of nonequilibrium macromolecules [46,57,58], but within simpler models that do not take into account complex dynamics of medium relaxation. Similar algorithms can also be developed for the generalized theory presented here, but these results will be reported elsewhere.

Verification of the Theory
The proposed theory can be verified by comparing it with earlier models, especially those that can be considered as special cases with respect to the present general approach. Hereafter, the following well-known frameworks are used for this verification: (1) the Najbar/Tachiya theory of two-stage ET in the Debye polar solvent [37], (2) the Zusman theory of ET in solvents with a two-component relaxation function [49]. Both these theories are special cases: the Najbar/Tachiya theory was developed for the three-center molecular compounds (N = 3), but only in a single-component solvent (R = 1); the Zusman theory is applicable to the multi-component environments (any R), but involves ET between only two redox centers (N = 2) and thus does not provide appropriate description of multistage reactions. The difference between these two models can also be clarified by comparing their configuration spaces. Both spaces are two-dimensional, but their nature is completely different. In the Najbar/Tachiya model, the free energy surfaces are constructed as functions of the polarization coordinates, while the Zusman model employs the relaxation coordinates. The two models thus operate in different subspaces as shown in Figure 1.
Since the theory presented in this paper is a generalization of the two approaches, we expect the new theory to reproduce the results of the previous ones in domains of their applicability according to the correspondence principle.
To demonstrate the correspondence to the Najbar/Tachiya model, we set τ 1 = τ 2 and eliminate the relaxation subspace R from consideration by projecting the composite space q onto the subspace Q This R-space folding gives the following coordinates of the FES minimǎ The G (n) FESs in the (Q 1 , Q 2 ) coordinates take the form of Equations (2). This result reproduces the Najbar/Tachya model [37], where two-stage ET is considered in terms of two-dimensional parabolic surfaces. It should also be noted that two-dimensional presentation of FESs in the form (2), (28) has been repeatedly used earlier (see, for example, [34,45]).
Consider the 3-center molecular system in more detail, assuming the electronic densities on D, A 1 , A 2 to be described by the effective radii r 1 , r 2 , r 3 , and the distances r 12 , r 13 , r 23 between the centers. This simple model is shown in Figure 2a, and allows us to estimate the system's energetics using the Marcus formula for the reorganization energy in a continuous dielectric medium Here e is the charge of an electron, c p = −1 o − −1 s is the Pekar factor, r i and r j are the effective radii of the donating and accepting sites, and r ij is the center-to-center distance between them. We adopt here the following model parameters that provide relatively high yields of charge separation [44] r 1 = r 2 = 4, r 3 = 8, r 12 = 10.4, r 23 = 12, r 13 = 19.4 (in angstroms).
Using these model parameters one can estimate the reorganization energies from Equation (24) λ (12) = 1.17 eV, λ (13) = 1.03 eV, λ (23) = 0.79 eV, θ = 50 • . (26) and calculate the diabatic FESs minima from Equations (12) To illustrate the R-space «folding» (Equation (22)), we calculate the projections of thě q (n) points into the Q subspace. They arě Figure 3 shows parabolic FESs calculated using Equation (28). Panel b) also shows the displacement vectors D (nn ) and the angle θ between the directions of the |ϕ 1 → |ϕ 2 and |ϕ 1 → |ϕ 3 transitions. The standard one-dimensional ET picture of two intersecting parabolic curves can be obtained by cutting the surfaces in Figure 3a by vertical planes passing through the minima of the reactant and product FESs. Now we check the present theory for its correspondence to the Zusman model of ET in solvents with two relaxation timescales [49]. To do this, we eliminate the polarization subspace Q by projecting the q space onto the relaxation subspace R. This projection is done as follows From Equations (12)-(14) one getš In the Zusman model, the FESs displacements in relaxation coordinate R i are proportional to the corresponding weight factors x i [49]. The same result is valid for our model, where theŘ (n) i coordinates are located along a straight line given by the equation (31)). Such an arrangement provides the required ratio of the relaxation components for all ET transitions in the system. Equation (17) for the system evolution is transformed into the Smoluchowski equation of diffusion in a two-dimensional parabolic potential. Diffusion coefficients along the R 1 and R 2 coordinates, however, are different: D 1 = k B T/τ 1 and D 2 = k B T/τ 2 . These equations reproduce the results of ref [49] as well.
Applying the Q-folding transformation to our model system (Equation (27)) we find the FESs minima in the R subspacě The G (n) curves on the (R 1 , R 2 ) plane are pictured in Figure 4.

Conclusions
Many studies of photoinduced ET in macromolecules are related to the prospects for their use as components of molecular electronics devices, in particular, solar cells, organic light-emitting diodes, optical sensors, switches, and others. Ultrafast photochemical reactions on the femtosecond timescale are often accompanied by the nonequilibrium states of the reactants and the environment. Multistage photoinduced ET in non-Debye environments, such as polymer mixtures, protein matrices, nanoaggregates, etc., is of particular interest. To describe the nonequilibrium effects in these reactions, we propose a method based on the splitting of the independent polarization coordinates into the relaxation components. The method can be considered as a generalization of the two wellknown approaches, one of which is used to describe multistage processes, and the other, to account for the multi-component relaxation of the environment. In this study, the properties of the extended space are studied, and the relationship between the reorganization free energy, the medium relaxation parameters, and the arrangement of the diabatic FESs in the extended configuration space is established. The proposed general approach is applied to a three-center molecular system in a medium with a two-component relaxation function. The algorithm for the diabatic FESs construction is described in detail, and a set of kinetic equations for the electronic density functions is specified. The general approach is verified by demonstrating its correspondence to the well-known Najbar/Tachiya and Zusman models.