The role of different reorganization energies within the Zusman theory of electron transfer

J. Chem. Phys. 118, 291 (2003); https://doi.org/10.1063/1.1525799 118, 291 © 2003 American Institute of Physics. The role of different reorganization energies within the Zusman theory of electron transfer Cite as: J. Chem. Phys. 118, 291 (2003); https://doi.org/10.1063/1.1525799 Submitted: 05 August 2002 . Accepted: 07 October 2002 . Published Online: 16 December 2002 Jesús Casado-Pascual, Manuel Morillo, Igor Goychuk, and Peter Hänggi ARTICLES YOU MAY BE INTERESTED IN Electron transfer dynamics: Zusman equation versus exact theory The Journal of Chemical Physics 130, 164518 (2009); https://doi.org/10.1063/1.3125003 Effect of friction on electron transfer in biomolecules The Journal of Chemical Physics 83, 4491 (1985); https://doi.org/10.1063/1.449017 Dynamic solvent effects on outer-sphere electron transfer The Journal of Chemical Physics 87, 2090 (1987); https://doi.org/10.1063/1.453184 The role of different reorganization energies within the Zusman theory of electron transfer Jesu ´s Casado-Pascual and Manuel Morillo Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, Sevilla 41080, Spain Igor Goychuk and Peter Ha ¨nggi Institut fu ¨r Physik, Universita ¨t Augsburg, Universita ¨tsstraße 1, D-86135, Augsburg, Germany 共Received 5 August 2002; accepted 7 October 2002兲 We consider the kinetics of electron transfer reactions in condensed media with different reorganization energies for the forward and backward processes. The starting point of our analysis is an extension of the well-known Zusman equations to the case of parabolic diabatic curves with different curvatures. Ageneralized master equation for the populations as well as formal expressions for their long-time limit is derived. We discuss the conditions under which the time evolution of the populations of reactants and products can be described at all times by a single exponential law. In the limit of very small tunnel splitting, a novel rate formula for the nonadiabatic transitions is obtained. It generalizes previous results derived within the contact approximation. For larger values of the tunnel splitting, we make use of the consecutive step approximation leading to a rate formula that bridges between the nonadiabatic and solvent-controlled adiabatic regimes. Finally, the analytical predictions for the long-time populations and for the rate constant are tested against precise numerical solutions of the starting set of partial differential equations. © 2003 American Institute of Physics. 关DOI: 10.1063/1.1525799兴 I. INTRODUCTION Electron transfer reactions are of prime importance in many physicochemical and biological processes.1At a very fundamental level, an electron transfer step is essentially an electron tunneling event in the presence of a medium 共solvent兲. A non-negligible tunneling probability requires resonance in the energy of localized electronic states. The solvent thermal fluctuations provide the necessary energy for the resonance condition. Thus, the kinetics of electron transfer reaction requires an adequate description of the medium thermal fluctuations that mediate the electronic charge redistribution in an electron transfer.2 In the classical theory of Marcus,3Hush,4and Levich and Dogonadze5the solvent fluctuations are described by an equilibrium probability law. Thus, the knowledge of the free energy as a function of an appropriate reaction coordinate is all that is needed to evaluate the rate constant. About 20 years ago, Zusman6and Alexandrov7introduced into the theory the idea that nonequilibrium effects associated with the relaxation of solvent fluctuations could also affect the rate. In a phenomenological way, Zusman proposed a set of four partial differential equations that incorporated the relaxation of the nonequilibrium probability laws of the reaction coordinate in each diabatic state and the tunneling transitions between them. In the original model, the diabatic states are parabolic functions of the reaction coordinate, with equal curvatures for the reactant and product states. The curvature is related to the reorganization energy in such a way that equal curvatures implies that the values of the reorganization energy for the forward and backward reactions are identical. Later on, Garg et al.8presented a derivation of the Zusman equations starting from a Hamiltonian model. They use techniques of functional integrals to carry out the elimination of the bath degrees of freedom from the total density operator. Previous analytical and numerical studies indicate that the reorganization energies for the direct and inverse reactions might indeed have different values.9–14 It is, therefore, interesting to extend the classical Zusman–Alexandrov formulation to the case of parabolas with different curvatures. A few years ago, Tang15 discussed such an extension. From the very beginning in his analysis, Tang made use of the socalled contact approximation. Namely, he assumed that tunneling between diabatic curves takes place strictly at their crossing points, thus neglecting any delocalization effects. Even for diabatic surfaces with equal curvatures, the fact that tunneling transitions are somewhat delocalized around the crossing point is important, especially in the inverted regime.16,17 Delocalization leads to modifications of the rate expression with respect to the typical Marcus–Arrhenius structure. In previous work,18 we have also found strong indications that the contact approximation is not adequate to describe strongly biased electron transfer processes. The starting point of this paper is a set of four partial differential equations similar to those used by Tang.15 The structure of such equations is the same as in the original Zusman model, but with parabolic potentials with different curvatures. The equations describe the dynamics of diagonal and off-diagonal matrix elements of the reduced density operator. We accept here the validity of Zusman model. Frantsuzov19 has correctly pointed out the limitations of Zusman equations to describe electronic transfer processes in strongly polar solvents. Indeed, if the solvent polarity is too strong, the conditions under which Zusman equations are derived from a more microscopic point of view20–23 might be JOURNAL OF CHEMICAL PHYSICS VOLUME 118, NUMBER 1 1 JANUARY 2003 2910021-9606/2003/118(1)/291/13/$20.00 © 2003 American Institute of Physics violated. This is an interesting point that we propose to address in the future. In the present work, we will concentrate on the derivation of the usual macroscopic kinetic description of electron transfer reactions from a Zusman-like model and on the determination of suitable expressions for the rate constant and the long-time populations of reactants and products. In our derivation, we will not assume the contact approximation. The scheme of the paper is as follows: In Sec. II we set up the model and the notation. In Sec. III we derive formal expressions for the populations on the diabatic surfaces in Laplace space. This is achieved by using Green functions and projection operator techniques on the Zusman equations. An alternative method of solution of the Zusman equations has been put forward recently by Cao and Jung.24 It is based on the spectral properties of the evolution operator for the density matrix in the Zusman approximation. As far as we know, this alternative has only been applied to diabatic parabolas of equal curvatures. The long-time limit is explicitly obtained. It is also shown that the populations satisfy integrodifferential equations with a complicated kernel. Under suitable conditions which are discussed in Sec. IV, we prove that the populations satisfy a single exponential relaxation law for all relevant time scales, characterized by the rate constant and the long-time values of the populations. In Sec. V, we present approximate analytical expressions for those two kinetic parameters. In Sec. VI we present a detailed comparison of the analytical predictions and precise numerical solution of the Zusman equations. Finally, we conclude with comments about the main findings in this work. Some of the calculations are very involved and they are presented in the Appendices. II. THE ZUSMAN EQUATIONS FOR DIABATIC POTENTIALS WITH DIFFERENT CURVATURES The basic elements to describe electron transfer processes are two diabatic electronic energy curves Vj(x), j ⫽1, 2, and a generalized one-dimensional reaction coordinate xwith effective mass m. The electronic states before and after the charge transfer will be denoted as donor, 兩1典, and acceptor, 兩2典, respectively. The reaction coordinate represents a combination of the selected nuclear modes coupled directly to the electronic transfer system.1The reaction coordinate is also coupled to the rest of nuclear modes. This coupling introduces friction in the dynamics of the reaction coordinate with a phenomenological friction coefficient ␩ . This is a well established starting point in the microscopic treatment of intramolecular electron transfer.1,8,14 It should be noted that our Vj(x) are potential energy curves. They should not be confused with the parabolic free energy curves appearing in alternative descriptions of electron transfer, as those relying on computer simulations.10 In the overdamped limit, Zusman equations provide an appropriate description for the time evolution of the matrix elements ␳ jk(x,t) ª 具 j,x 兩 ␳ ˆ(t) 兩 x,k 典 of the reduced density operator in the electron and reaction coordinate Hilbert space. Zusman equations and their validity conditions have been repeatedly derived and discussed in the literature.8,19–23 These equations read ⳵ ⳵ t ␳ 11共x,t兲⫽L ˆ1 ␳ 11共x,t兲⫹i⌬ 2ប关 ␳ 12共x,t兲⫺ ␳ 21共x,t兲兴,共1兲 ⳵ ⳵ t ␳ 22共x,t兲⫽L ˆ2 ␳ 22共x,t兲⫺i⌬ 2ប关 ␳ 12共x,t兲⫺ ␳ 21共x,t兲兴,共2兲 ⳵ ⳵ t ␳ 12共x,t兲⫽ 再 L ˆ⫺i ប关V1共x兲⫺V2共x兲兴 冎 ␳ 12共x,t兲 ⫹i⌬ 2ប关 ␳ 11共x,t兲⫺ ␳ 22共x,t兲兴,共3兲 ⳵ ⳵ t ␳ 21共x,t兲⫽ 再 L ˆ⫹i ប关V1共x兲⫺V2共x兲兴 冎 ␳ 21共x,t兲 ⫺i⌬ 2ប关 ␳ 11共x,t兲⫺ ␳ 22共x,t兲兴.共4兲 Here L ˆjare the Smoluchowski operators describing diffusion on each diabatic potential: L ˆj⫽D ⳵ ⳵ x 冋 ⳵ ⳵ x⫹Vj ⬘共x兲 kBT 册 .共5兲 The macroscopic diffusion constant Dis connected with the friction coefficient ␩ , which is assumed to be identical in both diabatic states, and the temperature Tby the Einstein relation D⫽kBT/ ␩ . The operator L ˆ⫽(L ˆ1⫹L ˆ2)/2 describes diffusion on the average potential 关V1(x)⫹V2(x)兴/2. Finally, ⌬denotes the electronic coupling matrix element, and it characterizes the degree of overlap of the donor and acceptor wave functions. Here, we will take ⌬to be independent of the nuclear coordinates 共Condon approximation兲. In this work, we will assume parabolic diabatic curves of the form Vj共x兲⫽m ␻ j 2 2共x⫺x0 ␦ j,2兲2⫺ ⑀ 0 ␦ j,2 ,共6兲 where x0and ⑀ 0are the horizontal and vertical shifts, respectively, between the minima of the parabolas 共cf. Fig. 1兲. The frequencies ␻ jcharacterize their curvatures, and they are FIG. 1. Parabolic diabatic surfaces with different curvatures as a function of the reaction coordinate x. Notice that the number of crossing points depends on the value of the energy bias ⑀ 0. 292 J. Chem. Phys., Vol. 118, No. 1, 1 January 2003 Casado-Pascual et al. related to the reorganization energies ␭jby the expression ␭j⫽m ␻ j 2x0 2/2. Therefore, the fact that the curvatures are different implies that these reorganization energies for the forward and backward reactions are also different. To avoid confusion, we want to emphasize that the curvatures of our diabatic potential energies, ␻ j, can be different. If one describes electron transfer processes in terms of free energy profiles, then, as pointed out by Tachiya,25 the two free energy curves are not independent, and they cannot be strictly parabolic when the curvatures at their minima are different. Obviously, our potential energies, not being free energies, are not tied up by such a restriction. The difference in curvatures also yields a difference in the phenomenological relaxation times corresponding to each diabatic curve, ␶ j⫽ ␩ m ␻ j 2⫽kBT m ␻ j 2D.共7兲 Notice that these relaxation times are related to the reorganization energies by ␶ 1/ ␶ 2⫽␭2/␭1. For later convenience, we will also introduce the relaxation time of the overdamped oscillator on the averaged potential 关V1(x)⫹V2(x)兴/2, ␶ ⫽2kBT mD共 ␻ 1 2⫹ ␻ 2 2兲.共8兲 From the above expressions, it follows that 2/ ␶ ⫽1/ ␶ 1 ⫹1/ ␶ 2. Electron tunneling is most effective near the crossing points of the diabatic curves, which are given by xj *⫽x0 ␭2⫺␭1 冋 ␭2⫹共⫺1兲j 冑 冉 1⫺ ⑀ 0 ⑀ c 冊 ␭1␭2 册 ,共9兲 where ⑀ c⫽␭1␭2/(␭1⫺␭2). Depending on the relative values of ⑀ 0and ⑀ cthere can be two, one, or no crossing points. III. FORMAL SOLUTION OF THE ZUSMAN EQUATIONS In this section we obtain some exact, though formal, analytical results for the time evolution of the populations Pj(t). They are obtained by integrating the corresponding probability densities over configuration space, i.e., Pj(t) ⫽ 兰 ⫺⬁ ⬁dx ␳ jj(x,t). The first step is to reduce the four Zusman equations to just two integral equations for the diagonal elements ␳ jj(x,t). This is achieved by first replacing ␳ 12(x,t) and ␳ 21(x,t) in Eqs. 共1兲and 共2兲by the result of formally solving the two off-diagonal equations 共3兲and 共4兲. This yields ⳵ ⳵ t ␳ jj共x,t兲⫽共⫺1兲j⌬ ប 冕 ⫺⬁ ⬁dx1Im关God共x,t 兩 x1兲 ␳ 12共x1,0兲兴 ⫹共⫺1兲j⌬2 2ប2 冕 0 tdt1 冕 ⫺⬁ ⬁dx1 ⫻Re关God共x,t⫺t1 兩 x1兲兴 ⫻关 ␳ 11共x1,t1兲⫺ ␳ 22共x1,t1兲兴⫹L ˆj ␳ jj共x,t兲,共10兲 where God(x,t 兩 x⬘), the off-diagonal Green function, is the solution of the partial differential equation ⳵ ⳵ tGod共x,t 兩 x⬘兲⫽ 再 L ˆ⫺i ប关V1共x兲⫺V2共x兲兴 冎 God共x,t 兩 x⬘兲, 共11兲 with initial condition God共x,0 兩 x⬘兲⫽ ␦ 共x⫺x⬘兲,共12兲 and boundary conditions lim x→⫾⬁ God共x,t 兩 x⬘兲⫽0. 共13兲 An evaluation of this Green function for the harmonic potentials with different curvatures can be found in Appendix A 1. From now on, we shall assume that ␳ 12(x,0)⫽0, so that the first term on the right-hand side of expression 共10兲will not be present. This initial condition describes the usual situation in which the electronic coherences between donor and acceptor states are initially neglected. Next, formal solution of Eq. 共10兲in terms of the diagonal Green functions Gd (j)(x,t 兩 x⬘) leads to ␳ jj共x,t兲⫽共⫺1兲j 冕 0 tdt1 冕 ⫺⬁ ⬁dx1aj共x,t⫺t1 兩 x1兲 ⫻关 ␳ 11共x1,t1兲⫺ ␳ 22共x1,t1兲兴 ⫹ 冕 ⫺⬁ ⬁dx1Gd (j)共x,t 兩 x1兲 ␳ jj共x1,0兲.共14兲 The propagators aj(x,t 兩 x⬘) in Eq. 共14兲are given in terms of the Green functions by aj共x,t 兩 x⬘兲⫽⌬2 2ប2 冕 0 tdt⬘ 冕 ⫺⬁ ⬁dx⬙Gd (j)共x,t ⫺t⬘ 兩 x⬙兲Re关God共x⬙,t⬘ 兩 x⬘兲兴.共15兲 The diagonal Green functions describe the diffusive motion on the diabatic curves Vj(x). They are the solution of the partial differential equations ⳵ ⳵ tGd (j)共x,t 兩 x⬘兲⫽L ˆjGd (j)共x,t 兩 x⬘兲,共16兲 with the same type of initial and boundary conditions as God(x,t 兩 x⬘). Explicit expressions for Gd (j)(x,t 兩 x⬘) can be found in Appendix A 2. Typically, a reaction starts from a situation where the solvent is thermally equilibrated. Thus, we will restrict our study to initial conditions for the diagonal terms of the form ␳ jj(x,0)⫽gj(x)Pj(0), where gj共x兲⫽exp关⫺Vj共x兲/kBT兴 冕 ⫺⬁ ⬁dx⬘exp关⫺Vj共x⬘兲/kBT兴 共17兲 are the equilibrium distributions on each of the diabatic curves, and Pj(0) are the initial conditions for the populations. Obviously, with these initial conditions, the second term on the right-hand side of expression 共14兲reduces to gj(x)Pj(0). 293J. Chem. Phys., Vol. 118, No. 1, 1 January 2003 Electron transfer with different reorganization energies For later convenience, we will rewrite Eq. 共14兲in matrix notation by introducing the one-column vector %(x,t) with components %j(x,t)⫽ ␳ jj(x,t), and the 2⫻2 matrices U, A(x,t 兩 x⬘), and g(x) with matrix elements Ujk⫽(⫺1)j⫹k, Ajk(x,t 兩 x⬘)⫽aj(x,t 兩 x⬘) ␦ j,k, and gjk(x)⫽gj(x) ␦ j,k, respectively. The solution of Eq. 共14兲is best sought by using the Laplace transform f ˜ (s)⫽ 兰 0 ⬁dt f(t)exp(⫺st). Thus, one finds % ˜ 共x,s兲⫽s⫺1g共x兲P共0兲⫺ 冕 ⫺⬁ ⬁dx1A ˜ 共x,s 兩 x1兲U% ˜ 共x1,s兲.共18兲 We will use standard projection operator techniques to get an expression for the Laplace transform of the populations P ˜ (s). We define the projection operators ⌸and Qas ⌸F共x兲⫽F储共x兲⫽g共x兲 冕 ⫺⬁ ⬁dx⬘F共x⬘兲,共19兲 QF共x兲⫽F共x兲⫽F共x兲⫺F 储 共x兲,共20兲 F(x) being an arbitrary one-column vector or 2⫻2 matrix depending on x. By acting with ⌸on Eq. 共18兲and integrating over x, one obtains after some simplifications P共0兲⫽ 冋 sI⫹ 冕 ⫺⬁ ⬁dx1K ˜ 共x1,s兲Ug共x1兲 册 P ˜ 共s兲 ⫹ 冕 ⫺⬁ ⬁dx1K ˜ 共x1,s兲U% ˜ 共x1,s兲,共21兲 where Iis the 2⫻2 unit matrix and K ˜ 共x,s兲⫽⌬2 2ប2 冕 ⫺⬁ ⬁dx⬘Re关G ˜ od共x⬘,s 兩 x兲兴.共22兲 Using Eq. 共A16兲in Eq. 共22兲, we obtain K ˜ 共x,s兲⫽⌬2 2ប2 冕 0 ⬁dt Re 兵 exp关⫺st⫺c共t,x兲兴 其 ,共23兲 with c(t,x) given by Eq. 共A11兲. The action of Qon Eq. 共18兲 leads to % ˜ 共x,s兲⫽⫺ 冕 ⫺⬁ ⬁dx1A ˜ 共x,s 兩 x1兲U关g共x1兲P ˜ 共s兲⫹% ˜ 共x1,s兲兴, 共24兲 where A ˜ 共x,s 兩 x⬘兲⫽⌬2 2ប2g共x兲 冕 ⫺⬁ ⬁dx⬙J ˜ 共x,s 兩 x⬙兲 ⫻Re关G ˜ od共x⬙,s 兩 x⬘兲兴,共25兲 with J ˜ (x,s 兩 x⬘) being the 2⫻2 diagonal matrix with matrix elements J ˜ jk共x,s 兩 x⬘兲⫽ ␦ j,k 兵 关gj共x兲兴⫺1G ˜ d (j)共x,s 兩 x⬘兲⫺s⫺1 其 .共26兲 A formal solution for % ˜ (x,s) is obtained by solving iteratively the integral equation 共24兲. Substitution of the result in Eq. 共21兲leads to P共0兲⫽关sI⫹Uk ˜ 共s兲兴P ˜ 共s兲,共27兲 where k ˜ (s) is the diagonal matrix obtained by summing the series k ˜ 共s兲⫽兺 n⫽1 ⬁ k ˜ (n)共s兲.共28兲 The terms of the series k ˜ (n)(s) are given by k ˜ (1)共s兲⫽ 冕 ⫺⬁ ⬁dx1K ˜ 共x1,s兲g共x1兲共29兲 and k ˜ (n)共s兲⫽共⫺1兲n⫺1 冕 ⫺⬁ ⬁dx1••• 冕 ⫺⬁ ⬁dxnK ˜ 共x1,s兲 ⫻ 兿 j⫽2 n Tr关A ˜ 共xj⫺1,s 兩 xj兲兴g共xn兲共30兲 for n⭓2. According to Eq. 共27兲, the Laplace transform of the populations can be expressed in terms of the diagonal matrix k ˜ (s)as P ˜ j共s兲⫽Pj共0兲⫹s⫺1关 ␦ j,1k ˜ 22共s兲⫹ ␦ j,2k ˜ 11共s兲兴 s⫹k ˜ 11共s兲⫹k ˜ 22共s兲.共31兲 After carrying out the inverse Laplace transform of the above expression, one obtains the time evolution of the populations, Pj(t). It follows from expression 共31兲that the longtime limit of the populations is given by Pj共⬁兲ªlim t→⫹⬁ Pj共t兲⫽lim s→0⫹ sP ˜ j共s兲 ⫽ ␦ j,1k ˜ 22共0兲⫹ ␦ j,2k ˜ 11共0兲 k ˜ 11共0兲⫹k ˜ 22共0兲.共32兲 Notice that these values are independent of the initial conditions Pj(0). Finally, rearranging Eq. 共27兲and carrying out the inverse Laplace transform, we find that the populations Pj(t) satisfy the set of generalized master equations d dtP1共t兲⫽⫺ 冕 0 tdt⬘关k11共t⫺t⬘兲P1共t⬘兲⫺k22共t⫺t⬘兲P2共t⬘兲兴, 共33兲 d dtP2共t兲⫽⫺ 冕 0 tdt⬘关k22共t⫺t⬘兲P2共t⬘兲⫺k11共t⫺t⬘兲P1共t⬘兲兴, where kjj(t) is the inverse Laplace transform of k ˜ jj(s). From the above set of equations, it follows immediately the conservation of probability, i.e., P1(t)⫹P2(t)⫽1. Thus, the set 共33兲reduces to a single integrodifferential equation for, say, the population P1(t). This equation can be conveniently written as d dtP1共t兲⫽⫺ 冕 0 tdt⬘关k11共t⬘兲⫹k22共t⬘兲兴P1共t⫺t⬘兲 ⫹ 冕 0 tdt⬘k22共t⬘兲.共34兲 Then, after solution of Eq. 共34兲, the evolution of P2(t) follows immediately. Up to now, the formal results that we have obtained are exact. We have only assumed the convergence of the series 共28兲and a specific family of initial conditions for the densi294 J. Chem. Phys., Vol. 118, No. 1, 1 January 2003 Casado-Pascual et al. ties ␳ jk(x,t). The result 共33兲is important. The populations on the donor and acceptor diabatic curves satisfy first order integrodifferential equations with a convolution structure. The convolution kernels, kjj(t), are rather complicated. The next important task is to analyze under which conditions can the solutions of Eq. 共33兲be properly approximated by a single exponential time evolution. IV. VALIDITY CONDITIONS FOR THE RATE REGIME In kinetics, the time evolution of the populations on the donor, P1(t), and acceptor, P2(t), diabatic curves is usually given by Pj共t兲⫽Pj共⬁兲⫹关Pj共0兲⫺Pj共⬁兲兴e⫺⌫t,共35兲 where ⌫is the total rate constant. In this section we analyze the conditions under which this rate description of the time evolution of the populations can be obtained from the Zusman equations. The starting point is the integrodifferential equation 共34兲. In this equation one can distinguish two different clear-cut time scales. The first one, ␶ a, is associated with the relaxation time of the kernels kjj(t). As all the time dependence of kjj(t) is through the diagonal and offdiagonal Green functions, ␶ adepends mainly on the relaxation times of these Green functions. The second one is given by ␶ b⫽关k ˜ 11(0)⫹k ˜ 22(0)兴⫺1and, as we will see below, it is associated with the relaxation time of the populations. In general, ␶ bdepends on the relaxation times of the Green functions and also on the characteristic tunneling time scale ប/⌬. In the next section, we will obtain approximate expressions for ␶ b ⫺1. In order to study the validity conditions for the rate regime, it is convenient to express Eq. 共34兲in dimensionless form as d dY P1共Y; ␨ 兲⫽⫺ 冕 0 Y/ ␨ dy⬘关 ␬ 11共y⬘兲⫹ ␬ 22共y⬘兲兴P1共Y⫺ ␨ y⬘; ␨ 兲 ⫹ 冕 0 Y/ ␨ dy⬘ ␬ 22共y⬘兲.共36兲 Here we have introduced the dimensionless quantities Y⫽t/ ␶ b,y⬘⫽t⬘/ ␶ a, ␨ ⫽ ␶ a/ ␶ b, and ␬ jj(y⬘) ⫽ ␶ a ␶ bkjj( ␶ ay⬘), and we have indicated explicitly the dependence of P1(Y; ␨ ) on the parameter ␨ . Notice that the dimensionless integration variable y⬘has been chosen so that the contributions of the integrand for values of y⬘much larger than unity can be safely neglected. Let us assume that we are in a regime in which ␶ a Ⰶ ␶ b. To find the leading-order approximation to the solution of Eq. 共36兲as ␨ →0⫹, i.e., P1(Y;0)⫽lim ␨ →0⫹P1(Y; ␨ ), we take the limit ␨ →0⫹in Eq. 共36兲while keeping Y⫽0 fixed. The result is d dY P1共Y;0兲⫽⫺ 冕 0 ⬁dy⬘关 ␬ 11共y⬘兲⫹ ␬ 22共y⬘兲兴P1共Y;0兲 ⫹ 冕 0 ⬁dy⬘ ␬ 22共y⬘兲⫽⫺P1共Y;0兲⫹P1共⬁兲, 共37兲 where we have taken into account that the values of the integrand for y⬘much larger than unity are negligible. The solution of Eq. 共37兲is P1(Y;0)⫽Ce⫺Y⫹P1(⬁), Cbeing an unknown constant of integration. This constant is determined by making use of the initial condition, i.e., limY→0⫹P1(Y;0)⫽C⫹P1(⬁)⫽P1(0). Therefore, after going back to the original variable t, we find that the leadingorder approximation to the solution of Eq. 共36兲as ␨ →0⫹is given by Eq. 共35兲with j⫽1 and ⌫⫽ ␶ b ⫺1. From the conservation of probability, it follows immediately that, in this limit, the population P2(t) is also of the form 共35兲. In conclusion, we have proved that, when ␶ aⰆ ␶ b, the values of the populations obtained from the Zusman equations can be properly approximated by the rate expression 共35兲with total rate ⌫⫽k ˜ 11共0兲⫹k ˜ 22共0兲.共38兲 Notice that, when ␶ aⰆ ␶ b, Eq. 共35兲describes properly the relaxation of the populations for all the relevant time scales associated with Pj(t), even for short times. When the condition ␶ aⰆ ␶ bis violated, then the discrepancies between the predictions of Eq. 共35兲and those of the Zusman equations might be relevant at all times, as we will see later, when we compare with numerical solution of Zusman equation 共see Sec. VI兲. V. ANALYTICAL EXPRESSIONS FOR THE LONG-TIME POPULATIONS AND THE TOTAL RATE CONSTANT From now on, we will assume that a rate description of the time evolution of the populations is appropriate. In that case, the parameters Pj(⬁) and ⌫can be expressed in terms of the diagonal matrix k ˜ (0), according to Eqs. 共32兲and 共38兲. The evaluation of this matrix entails the summation of the series 共28兲. In order to do so, one has to resort to approximations. The nature of the approximations is dictated by the relative values of the tunneling frequency and the characteristic frequencies of the solvent dynamics. A. The nonadiabatic limit If the characteristic tunneling frequency ⌬/បis very small relative to the relaxation frequencies associated with the Green functions, then tunneling becomes the limiting step mechanism of the rate process. In this nonadiabatic regime, the matrix elements of k ˜ (0) can be approximated as k ˜ jj共0兲⬇kNA (j)ª⌬2lim ⌬→0 k ˜ jj共0兲 ⌬2⫽k ˜ jj (1)共0兲.共39兲 After inserting the expressions 共17兲and 共23兲, with s⫽0, into Eq. 共29兲and integrating over x1, we find kNA (j)⫽ ⌬2⌳j 1/2 ប2 冕 0 ⬁dt Re 兵 Nj共t兲exp关Rj共t兲兴 其 ,共40兲 where 295J. Chem. Phys., Vol. 118, No. 1, 1 January 2003 Electron transfer with different reorganization energies Nj共t兲⫽ ␣ 1/2exp关共1⫺ ␣ 兲t/共2 ␶ 兲兴 兵 共 ␣ ⫹1兲关2⌳j⫹共 ␣ ⫺1兲/2兴⫹共 ␣ ⫺1兲关2⌳j⫺共 ␣ ⫹1兲/2兴exp关⫺2 ␣ t/ ␶ 兴 其 1/2 ,共41兲 Rj共t兲⫽i 冋 ␭1␭2共⌳1⫺⌳2⫹i ␹ 兲 ប共␭1⫹␭2兲 ␣ 2⫺ ⑀ 0 ប 册 t⫹4 ␶ ␭1 2␭2 2 ប␭j共␭1⫹␭2兲2 ␣ 3 ⫻ 再 关i共2⌳j⫹1兲共⫺1兲j⫺1⫹ ␹ 兴 ␣ sinh共 ␣ t/ ␶ 兲 4⌳j ␣ cosh共 ␣ t/ ␶ 兲⫹关4⌳j⫹ ␣ 2⫺1兴sinh共 ␣ t/ ␶ 兲 ⫹2⌳j关2i共⫺1兲j⫺1⫹ ␹ 兴关cosh共 ␣ t/ ␶ 兲⫺1兴 4⌳j ␣ cosh共 ␣ t/ ␶ 兲⫹关4⌳j⫹ ␣ 2⫺1兴sinh共 ␣ t/ ␶ 兲 冎 ,共42兲 and the dimensionless parameters ⌳j, ␹ , and ␣ are defined in Eqs. 共A12兲–共A14兲. The remaining time integral in Eq. 共40兲 can be calculated by a numerical quadrature. The expressions for the long-time populations and the rate constant in the nonadiabatic limit, Pj (NA)(⬁) and ⌫NA , are obtained by replacing k ˜ jj(0) with kNA (j)in Eqs. 共32兲and 共38兲. To the best of our knowledge, Eqs. 共40兲–共42兲have never been derived previously in the literature. These expressions for the nonadiabatic parameters are one of the main results in this paper. They constitute a generalization of the typical golden rule rate expressions to the case of diabatic parabolas with different curvatures. Notice that we are not assuming that the tunneling transitions are exactly localized at the crossing points of the parabolas as it is assumed in the so-called contact approximation.15,18 Actually, our expressions reflect the delocalization induced by the solvent dynamics. This delocalization might be important. As we have previously shown,18 detailed comparison of the results obtained with the rate formulas, with and without the contact approximation, and the results of the numerical solution of the Zusman equations indicates that the rate expressions with the contact approximation become invalid for large values of the energy bias, ⑀ 0. For the case of equal curvatures, ␭1⫽␭2⫽␭, Eq. 共40兲 reduces to the known result21 kNA (j)⫽⌬2 2ប2 冕 0 ⬁dt exp 冋 2␭kBT ␶ 2 ប2 冉 1⫺e⫺t/ ␶ ⫺t ␶ 冊 册 ⫻cos 冋 ␭ ␶ ប共1⫺e⫺t/ ␶ 兲⫹共⫺1兲j ⑀ 0 បt 册 .共43兲 In this case, the time integral can be evaluated explicitly in terms of the complete 关⌫(x)兴and incomplete 关⌫(x,y)兴 Gamma functions26 as kNA (j)⫽⌬2 ␶ 2ប2Re 兵 ebb⫺aj关⌫共aj兲⫺⌫共aj,b兲兴 其 ,共44兲 where we have introduced the dimensionless parameters aj ⫽ ␶ 关2␭kBT ␶ /ប⫺(⫺1)ji ⑀ 0兴/បand b⫽ ␶ ␭关2kBT ␶ /ប⫹i兴/ប. By the relation ␥ (a,x)ª⌫(a)⫺⌫(a,x)⫽a⫺1xae⫺xM(1,1 ⫹a,x),27 where M(a,b,c) is the Kummer’s function, our Eq. 共44兲is equivalent to Eq. 共3.13兲in Ref. 19. It should be noticed that, as Frantsuzov19 has pointed out, Eq. 共43兲can lead to nonphysical predictions, such as violation of the detailed balance and negative values for the rates in strongly polar solvents with large reorganization energies. As we discuss in Sec. VI, in the case of different reorganization energies and ⑀ 0⬎ ⑀ c, we have also observed deviations between the values of the long-time populations predicted by Eqs. 共32兲and 共40兲–共42兲and those expected from statistical thermodynamical considerations in the semiclassical limit. These anomalous results, which are intrinsic to the Zusman equations, can be used as a numerical criterion in order to test the validity of the Zusman description of ET reactions. The expressions for the nonadiabatic rate constants 共40兲–共42兲关or Eq. 共43兲in the case of equal curvatures兴simplify considerably if one makes use of the contact approximation. Within this approximation, one assumes that the electronic transitions take place precisely at the crossing points, so that, the function K ˜ (x,0) in Eq. 共29兲can be approximated by K ˜ 共x,0兲⯝ ␲ ⌬2 2ប ␦ 关V1共x兲⫺V2共x兲兴.共45兲 Then, the nonadiabatic rate constant, for ⑀ 0⬍ ⑀ c, can be expressed as15,18 kNA (1)⫽⌬2 4ប 冑 ␲ kBT␭2共1⫺ ⑀ 0/ ⑀ c兲 再 exp 冋 ⫺共␭2⫺ ⑀ 0兲2 4␭⫹共 ⑀ 0兲kBT 册 ⫹exp 冋 ⫺共␭2⫺ ⑀ 0兲2 4␭⫺共 ⑀ 0兲kBT 册 冎 ,共46兲 kNA (2)⫽ 冑 ␭2 ␭1exp 冉 ⫺ ⑀ 0 kBT 冊 kNA (1) ,共47兲 where we have defined the auxiliary, bias-dependent quantities ␭⫾共 ⑀ 0兲⫽关␭2⫾ 冑 共1⫺ ⑀ 0/ ⑀ c兲␭1␭2兴2 4␭1.共48兲 This contact approximation plays an essential role in Tang’s analysis of the Zusman equations.15 In the case of equal curvatures, the nonadiabatic rate constants in Eqs. 共46兲and 共47兲 reduce to the celebrated Marcus–Levich–Dogonadze rate3,5 and, therefore, they can be considered as its natural generalization to the case of different curvatures. B. The consecutive step approximation In order to go beyond the nonadiabatic limit, we need to evaluate the series in Eq. 共28兲. We will do this by extending the consecutive step approximation20,21 to the case of differ296 J. Chem. Phys., Vol. 118, No. 1, 1 January 2003 Casado-Pascual et al. ent reorganization energies for the forward and backward reactions. In this approximation, the terms of the series are simplified after disentangling the dynamical effects associated with diffusion from those relying on tunneling. We will assume that the function J ˜ (x,0 兩 x⬙) in Eq. 共26兲varies in x⬙ with a characteristic scale much larger than the width of the interval around x⬘where G ˜ od(x⬙,0 兩 x⬘) differs appreciably from zero. Then, according to Eqs. 共22兲and 共25兲, one can approximate A ˜ 共x,0 兩 x⬘兲⬇g共x兲J ˜ 共x,0 兩 x⬘兲K ˜ 共x⬘,0兲.共49兲 With this simplified expression for A ˜ (x,0 兩 x⬘), the exact expression in Eq. 共30兲for the terms k ˜ (n)(0) in the series expansion can be approximated by k ˜ (n)共0兲⬇共⫺1兲n⫺1 冕 ⫺⬁ ⬁dx1... 冕 ⫺⬁ ⬁dxnK ˜ 共xn,0兲g共xn兲 ⫻ 兿 j⫽2 n Tr关K ˜ 共xj⫺1,0兲g共xj⫺1兲J ˜ 共xj⫺1,0 兩 xj兲兴共50兲 for n⭓2. Henceforth, we will assume that there are two crossing points, i.e., ⑀ 0⬍ ⑀ c. Then, as it can be checked by numerical integration of Eq. 共23兲with s⫽0, the function K ˜ (x,0) shows peaks of similar heights and widths centered at the crossing points, at least when they are well separated. Assuming that the characteristic scale of variation of the functions J ˜ (x,0 兩 x⬘) and g(x) are also much larger than the widths of those peaks, we finally obtain that, for n⭓2, the matrices k ˜ (n)(0) can be well approximated by k ˜ (n)共0兲⬇ ␯ (n)k ˜ (1)共0兲,共51兲 where ␯ (n)⫽共⫺1兲n⫺1兺 j1⫽1 2 ••• 兺 jn⫽1 2 rj1•••rjn ⫻ 兿 l⫽2 n Tr关k ˜ (1)共0兲J ˜ 共xjl⫺1 *,0 兩 xjl *兲兴.共52兲 In the above expression, xj *represents the coordinate of the jth crossing point given by Eq. 共9兲, and the coefficients rj⫽g1共xj *兲 g1共x1 *兲⫹g1共x2 *兲 ⫽g2共xj *兲 g2共x1 *兲⫹g2共x2 *兲共53兲 denote the equilibrium weights of the two crossing points contributions. Then, according to Eq. 共28兲, we find that k ˜ 共0兲⬇ ␯ k ˜ (1)共0兲,共54兲 where ␯ is the result of summing up the series ␯ ⫽兺 n⫽1 ⬁ ␯ (n),共55兲 with ␯ (1)⫽1. In Appendix B we carry out the summation of this series explicitly 关cf. Eq. 共B7兲兴. The traces appearing in Eq. 共B7兲共the coefficients ⌶l,m) can be expressed in terms of the nonadiabatic rate constants, kNA (j),as Tr关k ˜ (1)共0兲J ˜ 共xl *,0 兩 xm *兲兴⫽兺 n⫽1 2kNA (n) kDlm (n).共56兲 Here, according to Eqs. 共26兲and 共A18兲, we have defined the coefficients kDlm (j)as 1 kDlm (j)ªlim s→0⫹ J ˜ jj共xl *,s 兩 xm *兲⫽ ␶ j 冕 0 ⬁dz 再 exp 冋 ␭j 2kBT 冉 共yl⫹ym⫺2 ␦ j,2兲2 ez⫹1⫺共yl⫺ym兲2 ez⫺1 冊 册 共1⫺e⫺2z兲1/2 ⫺1 冎 ,共57兲 where we have expressed the time integral in dimensionless units, and we have introduced the dimensionless coordinates of the crossing points yl⫽xl */x0. The coefficients kDlm (j)arise from the diffusional dynamics along the diabatic surfaces. The diagonal terms kDll (j)can be expressed through the generalized hypergeometric functions,26 2F2(a,b;c,d;z), as21,18 1 kDll (j)⫽ ␶ j 冋 ln2⫹ 2Eal (j) kBT2F2 冉 1,1; 3 2,2; Eal (j) kBT 冊 册 ,共58兲 where Eal (j)are the activation energies measured from the bottom of the diabatic potential Vj(x) to the crossing point xl *, i.e., Eal (j)⫽Vj(xl *)⫹ ⑀ 0 ␦ j,2 . Finally, taking into account Eqs. 共54兲,共56兲, and 共B7兲,we can conclude that, in the consecutive step approximation 共CSA兲, the matrix elements of k ˜ (0) can be approximated as k ˜ jj共0兲⬇kCSA (j) ª 1⫹r1r2兺 n⫽1 2 兺 l⫽1 2 兺 m⫽1 2 共⫺1兲l⫹mkNA (n) kDlm (n) 兿 l⫽1 2 冋 1⫹rl兺 n⫽1 2kNA (n) kDll (n) 册 ⫺r1r2 冋 兺 n⫽1 2kNA (n) kD12 (n) 册 2kNA (j). 共59兲 The expressions for the long-time populations and the rate constant in the consecutive step approximation, Pj (CSA)(⬁) and ⌫CSA , are obtained by replacing k ˜ jj(0) with kCSA (j)in Eqs. 共32兲and 共38兲. Notice that the equilibrium populations in the consecutive step approximation coincide with those obtained 297J. Chem. Phys., Vol. 118, No. 1, 1 January 2003 Electron transfer with different reorganization energies within the nonadiabatic limit, namely, Pj (CSA)(⬁) ⫽Pj (NA)(⬁). This can be easily seen after substitution of Eq. 共59兲into Eq. 共32兲. As we have previously analyzed,18 if the crossing point x2 *is much higher in energy than x1 *, then one can replace in Eq. 共59兲r1⇒1 and r2⇒0. In this case, expression 共59兲simplifies considerably to kCSA (j)⬇kNA (j) 1⫹kNA (1)/kD11 (1) ⫹kNA (2)/kD11 (2) .共60兲 The above formula has the same structure as the one used in the literature for equal curvatures.20,21 Equation 共59兲and its simplified version, Eq. 共60兲, describe in a unified way the different rate regimes, ranging from nonadiabatic to solvent controlled adiabatic reactions, depending upon the relative values of the system parameters characterizing tunneling and diffusion. The derivation of Eq. 共59兲is one of the main results of this paper, and as far as we know, it has never been obtained before. A few years ago, Tang15 arrived to an expression somewhat similar to Eq. 共59兲. A detailed comparison of our work and that of Tang reveals, nonetheless, some important differences. First, Tang neglects the off-diagonal diffusion terms, kD12 (j), present in our Eq. 共59兲. Second, the nonadiabatic rate constants appearing in Tang’s expression are the ones obtained within the contact approximation 关cf. Eqs. 共46兲–共48兲兴. VI. COMPARISON WITH NUMERICAL RESULTS AND DISCUSSION In this section we shall compare our analytical results with those provided by numerical integration of the Zusman equations. The latter has been carried out using the standard numerical algorithm group routine D03PCF on a LINUX PC with an Intel 800 MHz processor. In the numerical procedure, artificial absorbing boundary conditions have been properly superimposed far away from the reaction region, in order to model the natural boundary conditions, ␳ ij(x,t) →0atx→⫾⬁. Such a modeling did not affect the quality of the numerics, which was controlled by the numerical conservation of the total probability P1(t)⫹P2(t)⫽1 on the whole time scale. Namely, the deviation of the total probability from unity did not exceed 3⫻10⫺7for the mesh of 1500 space points and the single time step accuracy parameter of 10⫺7. We have adjusted both the number of mesh points and the time accuracy in order to achieve convergence of the results within the width of the plotted curves. Depending upon the values of the parameters, the calculation of a relaxation curve involving 100 time points took from about several seconds to about half an hour. The long-time population, P1(⬁), and the rate constant, ⌫, have been extracted from the numerical P1(t) making use of a nonlinear, singleexponential fitting procedure in GNUPLOT. The following set of parameter values is kept fixed in the calculations: ␭1⫽800 cm⫺1,␭2⫽200 cm⫺1,T⫽300 K. The other parameters, ⑀ 0,⌬, and ␶ , have been varied. Strongly different values for the reorganization energies ␭1and ␭2 have been chosen on purpose, in order to demonstrate the quality of our analytical results. In realistic situations, the difference between the reorganization energies may not be so dramatically large. For example, it was found in Ref. 10 that the reaction of primary charge separation in the bacterial photosynthetic center immersed in a nonpolar lipid membrane occurs with ␭1⬇1.45 kcal/mol⬇507.22 cm⫺1,␭2 ⬇1.55 kcal/mol⬇542.20 cm⫺1. In such a case, we expect our approximate results to work even better for a similar set of the remaining parameters. In Fig. 2 we show two typical numerical evolutions of P1(t) and their corresponding single-exponential fitting curves for a fixed value ⌬⫽20 cm⫺1and two different values of the relaxation time 共a兲 ␶ ⫽0.2 ps and 共b兲 ␶ ⫽2 ps. Figure 2共a兲demonstrates that the evolution is single exponential to a very good degree. The increase of ␶ by one order of magnitude 关cf. Fig. 2共b兲兴 introduces visible deviations from the strictly exponential behavior. In the following, we restrict our analysis to the case ⌬⭐10 cm⫺1and ␶ ⭐2.5 ps in order to ensure the strictly exponential character of the evolution. In Figs. 3, 4, and 5 we depict the numerical and analytical results for a fixed value ␶ ⫽1 ps and three different values of the tunneling matrix element, ⌬⫽1, 5, and 10 cm⫺1, respectively. In the three figures we find an excellent agreement between the numerics and our analytical theory. In particular, for weak tunneling, ⌬⫽1cm ⫺1, the transfer is nonadiabatic and the numerical transfer rate ⌫is perfectly reproduced by the nonadiabatic rate expression, Eqs. 共38兲 and 共40兲–共42兲, in the whole range of the electronic energy bias ⑀ 0关cf. Fig. 3共a兲兴. When ⌬increases, the nonadiabatic rate expression starts to fail 关cf. Figs. 4共a兲and 5共a兲兴, especially in the vicinity of the decoupling point ⑀ 0⫽ ⑀ cof the two diabatic energy surfaces ( ⑀ c⬇266 cm⫺1for the present parameters兲. Here, the adiabatic corrections due to the sluggish dynamics of the reaction coordinate become increasingly important as the nonadiabatic tunneling gets drastically accelerated. However, the numerical results are still pretty well reproduced by the consecutive step rate given in Eqs. 共38兲,共59兲,共40兲–共42兲,共53兲, and 共57兲. This agreement holds only in the range ⑀ 0⬍ ⑀ c, since for ⑀ 0⭓ ⑀ cthe consecutive FIG. 2. Comparison between the numerical results for the evolution of the donor population, P1(t), 共solid lines兲and their single-exponential fitting curves 共dashed lines兲for two different values of the relaxation time ␶ . The parameter values are ␭1⫽800 cm⫺1,␭2⫽200 cm⫺1,⌬⫽20 cm⫺1,T ⫽300 K, 共a兲 ␶ ⫽0.2 ps and 共b兲 ␶ ⫽2ps. 298 J. Chem. Phys., Vol. 118, No. 1, 1 January 2003 Casado-Pascual et al.