Relaxation in charge-transfer systems with very large tunnel splitting: A semiclassical stochastic approach

J. Chem. Phys. 113, 11176 (2000); https://doi.org/10.1063/1.1326907 113, 11176 © 2000 American Institute of Physics. Relaxation in charge-transfer systems with very large tunnel splitting: A semiclassical stochastic approach Cite as: J. Chem. Phys. 113, 11176 (2000); https://doi.org/10.1063/1.1326907 Submitted: 15 August 2000 . Accepted: 27 September 2000 . Published Online: 13 December 2000 J. Casado-Pascual, C. Denk, M. Morillo, and R. I. Cukier ARTICLES YOU MAY BE INTERESTED IN Nonadiabatic semiclassical dynamics in the mixed quantum-classical initial value representation The Journal of Chemical Physics 148, 102326 (2018); https://doi.org/10.1063/1.5005557 Quantum mechanical transition state theory and a new semiclassical model for reaction rate constants The Journal of Chemical Physics 61, 1823 (1974); https://doi.org/10.1063/1.1682181 Is the direct observation of electronic coherence in electron transfer reactions possible? The Journal of Chemical Physics 107, 8397 (1997); https://doi.org/10.1063/1.475040 Relaxation in charge-transfer systems with very large tunnel splitting: A semiclassical stochastic approach J. Casado-Pascual and C. Denk Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, Sevilla 41080, Spain M. Morillo Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, Sevilla 41080, Spain, and Department of Chemistry and Center for Fundamental Materials Research, Michigan State University, East Lansing, Michigan 48824-1322 R. I. Cukier Department of Chemistry and Center for Fundamental Materials Research, Michigan State University, East Lansing, Michigan 48824-1322 共Received 15 August 2000; accepted 27 September 2000兲 Electron transfer in strongly coupled systems, appropriate to mixed-valence compounds, is studied to explore the competition between electronic coherence and dissipation. A set of stochastic equations is derived for a spin-boson Hamiltonian with large tunneling coupling matrix element 共adiabatic regime兲and strong system-bath-coupling. The bath dynamics is treated classically while the quantum character of the system is maintained. The bath dynamics is affected by the system dynamics, the effect being included by a mean-field description, valid for the adiabatic regime. Numerical solutions of the stochastic equations are presented and compared with exact quantum mechanical results. The numerical implementation of the method is straightforward and the long-time behavior of the system can be accessed. Analytic equilibrium solutions for the adiabatic regime are obtained, and we find good agreement between the long-time solution of the stochastic equations and these equilibrium solutions. We examine the dependence of the electronic population on the initial preparation of the bath and find that the proportion between oscillation 共coherence兲and decay 共dissipation兲is quite sensitive to this initial condition. © 2000 American Institute of Physics. 关S0021-9606共00兲50448-6兴 I. INTRODUCTION The phenomena of decoherence and relaxation in dissipative quantum systems are a general problem in condensedphase chemistry and physics. Their study has stimulated the development of several methods of analysis that rely on the regime to be studied.1–18 The quantum nature of the interrogated system is always assumed, while the bath that couples to the quantum system may be classical or quantum, depending on the values of its characteristic parameters, the systembath-coupling strength and the temperature, T. In many cases, the potential energy relief of the quantum system supports two minima well separated by a barrier. Then, if the temperature is low enough, a reduced description of the quantum system in terms of its first two localized states is adequate. The strength of the coupling between the two states of the system, បV, relative to the cutoff energy of the system-bath interaction ប ␻ cand to kBTare key to the kind of approximation scheme to be used. In this work, we continue the development of a semiclassical method19–21 that we have found to be useful for studying the competition between electronic coherence and dephasing in optical spectroscopy, where a strong external field is applied to a quantum system with just a few degrees of freedom. The quantum system is coupled bi-linearly to a bath of harmonic oscillators, providing a relaxation mechanism. When restricted to a two-state system, the Hamiltonian of the isolated system under the influence of an external field may be written in pseudospin form, HS⫽⫺ ⌬G0 2 ␴ z⫹បV共t兲 ␴ x,共1兲 where ⌬G0is the free energy difference 共asymmetry兲and V(t) is proportional to the external field that induces transitions between the two states. In the case of strong field spectroscopy,19,20 the parameter V(t) cannot be considered small, and a treatment based on perturbation theory in the strength of the external field is, therefore, not feasible. When V(t) is large, there are multiple crossings between the states and quantum coherence effects are expected to be important. Then, a semiclassical approach can only be accurate under special circumstances. In a previous study, restricted to a ‘‘bath’’ consisting of one oscillator, we were able to show by numerical calculation and Landau–Zener analysis that semiclassical and quantum mechanical methods will agree in an adiabatic limit that is guaranteed when V(t) is large enough.21 The approximation consists of replacing expectation values of the products of system-oscillator operators with products of their expectation values in the exact equations of motion for the spin operators. In this fashion, the differences in the force exerted on the two states is effectively averaged by the fast transitions induced by the large value of V(t). We found that the more classical the oscillaJOURNAL OF CHEMICAL PHYSICS VOLUME 113, NUMBER 24 22 DECEMBER 2000 111760021-9606/2000/113(24)/11176/11/$17.00 © 2000 American Institute of Physics tor, the more adiabatic the regime. Furthermore, V(t) does not have to be extremely large 共in terms of the value of the corresponding Landau–Zener parameter兲in order for this semiclassical procedure to work quite well. A similar system Hamiltonian has been applied by many authors to study electron transfer 共ET兲reactions, with V(t) ⫽Vbeing the tunneling coupling element.3,2,4,6 ET reactions in mixed-valence compounds22,23 are characterized by large couplings between the electronic states 共large V) and large system-bath-couplings, rendering perturbation techniques inadequate. The case of large VET has stimulated several theoretical developments.24–29 The spectral density describing the system-bath-coupling is characterized by a cutoff frequency, ␻ c. For typical polar liquids and experiments at T⫽300 K, kBT⬎ប ␻ c.3Thus, a classical treatment of the bath seems adequate for these situations, although the quantum character of the system has to be maintained. Therefore, the one-oscillator semiclassical method just discussed should provide a good starting point for this investigation, as it is best suited to large V(t) and to oscillators that can be treated classically from the thermal (kBT⬎ប ␻ c) and quantum (V(t)⬎ ␻ c) points of view. Our calculations will be relevant in the analysis of photoinitiated ET reactions30 involving three electronic states. Before exciting with a laser pulse, the ground-state solute ( 兩 g 典 ) is stable and the solvent is thermally equilibrated to the current electronic charge distribution. A laser pulse puts the solute in an electronically excited state, 兩 0 典 , with the solvent distribution unchanged. This state will be our initial state for a subsequent electron transfer process to a third electronic state, 兩 1 典 . The tunneling matrix element Vand the exothermicity ⌬G0characterize the ET between the 兩 0 典 and 兩 1 典 states, 共cf. Fig. 1兲. It is important to notice that the solvent is initially in a nonequilibrium situation with respect to the charge distribution of the initial state pertinent to the ET process. There are two possibilities to consider: 共a兲V⬍ ␻ c 共fast bath兲. Here, the bath relaxation toward equilibrium takes place on a time scale that is short compared to the time scale of the electron transfer. Then, the usual expression for the ET rate obtained from a Golden Rule 共GR兲approach under the assumption of an initially equilibrated solvent 关for a large enough reorganization energy (Er)] will suffice to describe the process, at least during a relevant time range. On the other hand, 共b兲V⬎ ␻ c共slow bath兲. The appropriate bath initial condition is the nonequilibrium one. There can be several different stages in the time evolution of the electronic population. For short times, oscillations of the electronic populations, arising from the coupling between the 兩 0 典 and 兩 1 典 states can exist. These oscillations are probably too fast to be observed experimentally.30 The amplitude of the oscillations will decay due to the solvent fluctuations. As we shall see, the proportion between oscillation and decay in the overall time behavior will depend quite sensitively on the initial condition of the bath. The influence of the initial condition on population evolution has been pointed out by other investigators.31,24 When VⰆ ␻ c共the fast bath limit兲, approximate perturbative schemes have been developed to deal with the dynamics of the reduced density operator for the system.7,8,15,9,10,32 Even though Vis small, if the system-bath-coupling strength is not too strong and for low enough temperature, coherent oscillations can be found. If Ergets large enough, then a rate regime is obtained where the population decays exponentially with rate constants given by the Marcus–Levich electron transfer 共Golden Rule兲expression. Various numerical methods have also been developed to obtain the dynamics of the matrix elements of the reduced system density operator. The Tensor multiplication 共TensMult兲method of Makri and Makarov33,34 is based on a real time path-integral formulation. The practical applicability of the TensMult scheme relies on the fast decay of the memory kernel of the influence functional and this is favored by fast baths 共large ␻ c). The TensMult method describes correctly the long-time behavior of dissipative quantum systems; the density matrix has proven to reproduce the correct equilibrium values. In this paper, we mainly focus on slow baths, V⬎ ␻ c. In this case, it is difficult to obtain reliable results from the TensMult scheme. Stock35,36 recently presented a semiclassical method appropriate to the slow bath limit. In his method, the overall density operator is assumed to factorize at all times into a product of system and bath reduced density operators. The system is described by a finite number of oscillators whose frequencies are distributed in the same form for both diabatic states. Consequently, the scheme is valid up to times t Ⰶ1/⌬ ␻ , where ⌬ ␻ is the frequency separation between the bath oscillators. In order to access the long-time behavior, a large number of bath oscillators has to be considered, which increases the computational effort. Pechukas et al.18 have combined the Makri–Makarov and Stock procedures to give promising results, when the bath spectral density can be partitioned into separate contributions from slow and fast bath modes. They have also implemented a memory equation algorithm for the fast bath part that introduces little loss of accuracy and is numerically efficient.37 In this paper we focus on the regime of classical bath dynamics (kBT⬎ប ␻ c) and mainly on the slow baths, V FIG. 1. Schematic solvent free energy surfaces for a photoinitiated electron transfer reaction. The lower surface 兩 g 典 corresponds to a neutral solute state where the solvent equilibrium is independent of the solute’s presence, ␣ ⫽0. The solvent equilibrium points for the excited reactant, 兩 0 典 ,共product, 兩 1 典 ) states are displaced, relative to this neutral state, as parametrized by the values ␣ ⫽1(⫺1), cf. Eq. 共5兲. Photochemical excitation 兩 g 典 → 兩 0 典 initiates the reaction 兩 0 典 ↔ 兩 1 典 . If the excited state reaction is fast compared with solvent equilibration times, then an appropriate solvent initial condition is that of the ground equilibrium state, ␣ ⫽0. 11177J. Chem. Phys., Vol. 113, No. 24, 22 December 2000 Relaxation in charge-transfer systems ⬎ ␻ c. This second condition precludes the use of perturbation theory in the strength of V. Furthermore, we will assume that the system-bath-coupling energy, as characterized by Er, is not necessarily small either. Our aim is to construct an approximate stochastic description based on the idea that due to the large value of V, transitions among the diabatic states are so fast on the time scale of the bath dynamics, that the bath itself can only follow the system in an average sense, as discussed above for the one-oscillator ‘‘bath.’’21 This procedure provides a feedback of the electronic population into the dynamics of the bath that is absolutely essential for obtaining reasonable behavior for all times. In particular, the long-time average behavior of the system density matrix agrees with the predictions of an analytic expression for the equilibrium density matrix that we derive in this work. This latter expression is based on tracing out the bath degrees of freedom from the overall canonical density operator in the hightemperature limit. The use of a stochastic description permits the introduction of an infinite number of oscillators in a compact fashion and, as we shall see, leads to a very efficient numerical scheme. The plan of the rest of the paper is as follows. In Sec. II, we use the spin-boson Hamiltonian as a starting point to derive a set of stochastic equations of motion for the matrix elements of the system density operator and the stochastic process representing the bath reaction coordinate. In Sec. III, we obtain an expression for the equilibrium density matrix of the system that is a generalization to larger values of Vof the equilibrium density matrix adequate for small Vthat is given in terms of the standard reaction free energy ⌬G0and temperature T.4In Sec. IV, we discuss the solution of the stochastic equations for a wide range of parameter values. Our concluding remarks are presented in Sec. V. II. THE HAMILTONIAN AND EFFECTIVE EQUATIONS OF MOTION A model Hamiltonian frequently used to characterize an electron transfer between donor and acceptor centers in condensed media is1,2,4,3 HT⫽兺 j 冋 pj 2 2mj ⫹1 2mj ␻ j 2 冉 qj⫺ ␥ j mj ␻ j 2 ␴ z 冊 2 册 ⫺⌬G0 2 ␴ z⫹បV ␴ x.共2兲 The two electronic states are denoted by the kets 兩 0 典 and 兩 1 典 . The solvent is represented by a set of independent harmonic oscillators with shifted centers of oscillations. These shifts depend substantially upon the electronic state of the solute and reflect the differing interaction energies between the two charge distributions of the solute with the solvent. The mj’s and ␻ j’s are, respectively, the oscillator masses and frequencies, and the ␥ j’s are the solute-solvent coupling constants. The ␴ i(i⫽x,y,z) are the Pauli spin operators. Once the solute is reduced to only two active states, it can be described by a spin variable, and we shall use solute and spin language interchangeably. The quantity ⌬G0is the standard free-energy difference between reactants and products. The electronic coupling element Vis responsible for the ET between reactants and products. The equation of motion for the density operator of the total system ␳ T(t)is iប ␳ ˙T共t兲⫽关HT, ␳ T共t兲兴.共3兲 We assume that the initial density operator can be represented as ␳ T共0兲⫽ ␳ S共0兲丢 ␳ B共0兲,共4兲 where ␳ S(0) describes the initial state of the system, and ␳ B(0), theinitial density operator for the bath, is of the form ␳ B共0兲⫽1 ZB共 ␤ 兲exp 冋 ⫺ ␤ 兺 j 冉 pj 2 2mj ⫹1 2mj ␻ j 2 ⫻ 冉 qj⫺ ␣␥ j mj ␻ j 2 冊 2 冊 册 .共5兲 Here, ZB( ␤ ) is the normalization function. This distribution describes an equilibrium ensemble of independent harmonic oscillators with center of oscillations shifted by a quantity parameterized by ␣ . This shift reflects the fact that the solvent equilibrium depends upon the solute charge distribution. Adjusting ␣ will permit us to fix the initial bath preparation. The value ␣ ⫽1 corresponds to the standard initial condition of nonadiabatic ET reactions, where the bath is equilibrated to the charge distribution of the reactant state, while the value ␣ ⫽0 would describe a bath at equilibrium with a neutral reactant state. Finally, we shall consider that the system is initially in a pure state of the form 兩 ␾ 典具 ␾ 兩 . We will be interested in situations where the temperature is large, compared to the characteristic energy of the bath. Under these conditions, a classical description for the bath dynamics, where the qj(t) are treated as ordinary functions of time, is expected to be a good approximation. On the other hand, the system evolution must be described quantum mechanically. When employing this approximation, a difficulty arises: The coupling term ⫺ ␴ z兺j ␥ jqjin the Hamiltonian, Eq. 共2兲, is an operator in the system Hilbert space that influences the bath dynamics. We will introduce an approximation for this term, in order to ensure that the bath variables will maintain their classical character. Following these considerations, our first step is to introduce an effective Hamiltonian for the bath dynamics of the form HB eff共q,p,t兲⫽兺 j 冋 pj 2 2mj ⫹1 2mj ␻ j 2qj 2 册 ⫺f共t兲兺 j ␥ jqj. 共6兲 The idea is to approximately incorporate the feedback of the spin dynamics into the bath evolution equations with the term ⫺f(t)兺j ␥ jqj, where f(t) is a function of time that contains information about ␴ z(t). The explicit form of f(t) will be specified below. Similar approximations have been considered by other authors to analyze the influence of classical variables on the tunneling of quantum objects.16,35 This approximation amounts to neglecting quantum correlations 11178 J. Chem. Phys., Vol. 113, No. 24, 22 December 2000 Casado-Pascual et al. among the spin operators in the Heisenberg picture at different times. According to Eq. 共6兲, the Hamiltonian equations for the classical variables q(t),p(t) are q ˙j共t兲⫽ ⳵ HB eff共q,p,t兲 ⳵ pj ⫽pj共t兲 mj,共7兲 p ˙j共t兲⫽⫺ ⳵ HB eff共q,p,t兲 ⳵ qj ⫽⫺mj ␻ j 2qj共t兲⫹ ␥ jf共t兲.共8兲 This system of differential equations can be easily solved, and the result for qj(t)is qj共t兲⫽qj共0兲cos共 ␻ jt兲⫹pj共0兲 mj ␻ jsin共 ␻ jt兲 ⫹ ␥ j mj ␻ j 冕 0 tdsf共s兲sin共 ␻ j共t⫺s兲兲.共9兲 The influence of the dynamical evolution of the bath into the dynamics of the spin subsystem can then be described in terms of an effective Hamiltonian for the spin variables, given by HS eff共t兲⫽⫺ 冉 ⌬G0 2⫹兺 j ␥ jqj共t兲 冊 ␴ z⫹បV ␴ x ⫽⫺ ប 2 冉 ⌬G0 ប⫹ ␩ ␣ 共t兲⫹2 ␣ ប兺 j ␥ j 2 mj ␻ j 2 ⫻cos共 ␻ jt兲⫹ ␺ 共t兲 冊 ␴ z⫹បV ␴ x,共10兲 where ␩ ␣ 共t兲⫽2 ប兺 j ␥ j 冋 冉 qj共0兲⫺ ␣␥ j mj ␻ j 2 冊 cos共 ␻ jt兲 ⫹pj共0兲 mj ␻ jsin共 ␻ jt兲 册 ,共11兲 and ␺ 共t兲⫽2 ប 冕 0 tds 兺 j ␥ j 2 mj ␻ jsin共 ␻ j共t⫺s兲兲f共s兲.共12兲 If we introduce the spectral density J( ␻ ) J共 ␻ 兲⫽ ␲ 2兺 j ␥ j 2 mj ␻ j ␦ 共 ␻ ⫺ ␻ j兲,共13兲 we will have 兺 j ␥ j 2 mj ␻ j 2cos共 ␻ jt兲⫽2 ␲ 冕 ⫺⬁ ⬁d ␻ J共 ␻ 兲 ␻ cos共 ␻ t兲.共14兲 From now on, we shall assume that the coupling constants and the solvent frequencies are distributed with the Debye spectral density38 J共 ␻ 兲⫽Er 2 ␻ / ␻ c 1⫹共 ␻ / ␻ c兲2 ␪ 共 ␻ 兲,共15兲 and hence 兺 j ␥ j 2 mj ␻ j 2cos共 ␻ jt兲⫽Er 2e⫺ ␻ ct.共16兲 With this assumption, one obtains HS eff共t兲⫽⫺ ប 2 冉 ⌬G0 ប⫹ ␩ ␣ 共t兲⫹ ␣ Er បe⫺ ␻ ct⫹ ␺ 共t兲 冊 ␴ z ⫹បV ␴ x,共17兲 and ␺ 共t兲⫽Er ␻ c ប 冕 0 tdse⫺ ␻ c(t⫺s)f共s兲,共18兲 which is the solution of the differential equation ␺ ˙共t兲⫽⫺ ␻ c ␺ 共t兲⫹Er ␻ c បf共t兲,共19兲 with the initial condition ␺ (0)⫽0. The system density operator ␳ S(t) develops in the Schro ¨dinger picture according to the effective dynamical equation iប ␳ ˙S共t兲⫽关HS eff共t兲, ␳ S共t兲兴,共20兲 with the initial condition 兩 ␾ 典具 ␾ 兩 . This differential equation depends on the initial values 兵 q(0),p(0) 其 through the stochastic process ␩ ␣ (t). According to Eqs. 共5兲and 共11兲, ␩ ␣ (t) is the stationary, Markovian and Gaussian process known as the Ornstein–Uhlenbeck 共OU兲process,39 and hence it is fully specified by its average 具 ␩ ␣ 共t兲 典 B⫽0, 共21兲 and its second moment 具 ␩ ␣ 共t兲 ␩ ␣ 共t⬘兲 典 B⫽⌬2e⫺ ␻ c 兩 t⫺t⬘ 兩 ,共22兲 where 具典 Bmeans the average over the bath degrees of freedom, and ⌬⫽ 冑 2Er/(ប2 ␤ ). This process can be generated by the solution of the stochastic differential equation ␩ ˙ ␣ 共t兲⫽⫺ ␻ c ␩ ␣ 共t兲⫹w共t兲,共23兲 where w(t) is a zero-average white noise with correlation function 具 w(t)w(s) 典 ⫽2⌬2 ␻ c ␦ (t⫺s). The initial value ␩ ␣ ⫽ ␩ ␣ (0) is distributed according to the Gaussian distribution P共 ␩ ␣ 兲⫽1 冑 2 ␲ ⌬2e⫺共 ␩ ␣ 2/2⌬2兲.共24兲 It is important to note that, for a given realization of the OU process ␩ ␣ (t), the solution of Eq. 共20兲can be formally written as ␳ S共t兲⫽U共t兲 兩 ␾ 典具 ␾ 兩 U†共t兲,共25兲 where U(t) is the unitary evolution operator that satisfies the equation iបU ˙共t兲⫽HS eff共t兲U共t兲.共26兲 The initial condition is U(0)⫽I, with Ithe identity operator in the system Hilbert space. Hence, Eq. 共20兲maintains the normalization condition TrS( ␳ S(t))⫽1, where TrSis the 11179J. Chem. Phys., Vol. 113, No. 24, 22 December 2000 Relaxation in charge-transfer systems trace in the system Hilbert space. As ␳ S(t) is an operator in the system Hilbert space, it can be expanded in terms of the identity and Pauli operators ␳ S共t兲⫽1 2关I⫹x共t兲 ␴ x⫹y共t兲 ␴ y⫹z共t兲 ␴ z兴.共27兲 The above system of equations can be closed with a choice for the function f(t). Motivated by the large Vadiabatic regime, where the coupling is large, compared with the bath fluctuation speed, we choose21 f共t兲⫽TrS共 ␳ S共t兲 ␴ z兲⫽ 具 0 兩 ␳ S共t兲 兩 0 典 ⫺ 具 1 兩 ␳ S共t兲 兩 1 典 ⫽z共t兲. 共28兲 Using Eqs. 共27兲and 共28兲, in conjunction with Eqs. 共20兲and 共19兲,weget x ˙共t兲⫽ 冋 ⌬G0 ប⫹ ␩ ␣ 共t兲⫹ ␣ Er បe⫺ ␻ ct⫹ ␺ 共t兲 册 y共t兲,共29兲 y ˙共t兲⫽⫺ 冋 ⌬G0 ប⫹ ␩ ␣ 共t兲⫹ ␣ Er បe⫺ ␻ ct⫹ ␺ 共t兲 册 x共t兲 ⫺2Vz共t兲,共30兲 z ˙共t兲⫽2Vy共t兲,共31兲 ␺ ˙共t兲⫽⫺ ␻ c ␺ 共t兲⫹Er ␻ c បz共t兲.共32兲 These differential equations, together with Eq. 共23兲, constitute the set of stochastic differential equations that yield the evolution of ␳ S(t). They must be solved with the initial conditions x(0)⫽2Re 关 具 0 兩 ␾ 典具 1 兩 ␾ 典 *兴,y(0)⫽⫺2Im 关 具 0 兩 ␾ 典 ⫻ 具 1 兩 ␾ 典 *兴,z(0)⫽2 兩 具 0 兩 ␾ 典 兩 2⫺1, ␺ (0)⫽0, and Eq. 共24兲. Obviously, the actual system density operator is the one that is obtained after averaging over the bath variables, 具 ␳ S(t) 典 , or, equivalently, after averaging over the realizations of the OU process ␩ ␣ (t). If the initial system density operator does not correspond to a pure state but to a mixture of pure states of the form ␳ S(0)⫽c1 兩 ␾ 1 典具 ␾ 1 兩 ⫹c2 兩 ␾ 2 典 ⫻ 具 ␾ 2 兩 , with c1⫹c2⫽1, the method we have just described must be applied to each pure state independently, obtaining 具 ␳ S (j)(t) 典 (j⫽1,2), and the final result will be 具 ␳ S共t兲 典 ⫽c1 具 ␳ S (1)共t兲 典 ⫹c2 具 ␳ S (2)共t兲 典 .共33兲 Finally, for the purposes of numerical analysis and presentation of the results, it is convenient to define the dimensionless parameters t ˜ ⫽ ␻ ct,E ˜ r⫽Er/(ប ␻ c), T ˜ ⫽( ␤ ប ␻ c)⫺1,⌬G ˜ 0⫽⌬G0/(ប ␻ c), V ˜ ⫽V/ ␻ c, and the dimensionless functions ␩ ˜ ␣ (t ˜ )⫽ ␩ ␣ (t ˜ / ␻ c)/⌬and ␺ ˜ (t ˜ ) ⫽ប ␺ (t ˜ / ␻ c)/Er. Using these definitions, Eqs. 共29兲–共32兲and Eq. 共23兲can be expressed as x ˙共t ˜ 兲⫽关⌬G ˜ 0⫹共2E ˜ rT ˜ 兲1/2 ␩ ˜ ␣ 共t ˜ 兲⫹ ␣ E ˜ re⫺t ˜ ⫹E ˜ r ␺ ˜ 共t ˜ 兲兴y共t ˜ 兲,共34兲 y ˙共t ˜ 兲⫽⫺关⌬G ˜ 0⫹共2E ˜ rT ˜ 兲1/2 ␩ ˜ ␣ 共t ˜ 兲⫹ ␣ E ˜ re⫺t ˜ ⫹E ˜ r ␺ ˜ 共t ˜ 兲兴x共t ˜ 兲⫺2V ˜ z共t ˜ 兲,共35兲 z ˙共t ˜ 兲⫽2V ˜ y共t ˜ 兲,共36兲 ␺ ˜ ˙共t ˜ 兲⫽⫺ ␺ ˜ 共t ˜ 兲⫹z共t ˜ 兲,共37兲 ␩ 8 ␣ 共t ˜ 兲⫽⫺ ␩ ˜ ␣ 共t ˜ 兲⫹w ˜ 共t ˜ 兲,共38兲 where w ˜ (t ˜ ) is a zero-average white noise with correlation function 具 w ˜ (t ˜ )w ˜ (t ˜ ⬘) 典 ⫽2 ␦ (t ˜ ⫺t ˜ ⬘), and the initial value ␩ ˜ ␣ (0) is distributed according to Eq. 共24兲with ⌬⫽1. The interpretation of the above equations of motion is as follows: Eq. 共38兲describes the autonomous part of the bath dynamics that gives rise to a level broadening in the effective Hamiltonian Eq. 共17兲. The initial displacement provides the exponentially decaying terms proportional to ␣ Erin the effective Hamiltonian and the corresponding equations of motion. The dynamics of the isolated system is determined by the parameters V ˜ and ⌬G ˜ 0. Equation 共37兲describes the part of the bath dynamics affected by the system dynamics in terms of the population difference z(t ˜ ) and enters the system equations as a feedback term. The presence of the feedback ␺ ˜ (t ˜ ) in the equations of motion is essential for the system to reach an equilibrium that depends on ⌬G ˜ 0. When the feedback is ignored 关 ␺ ˜ (t ˜ ) ⫽0兴, we find, by numerically solving Eqs. 共34兲–共38兲, that 具 z(t ˜ ) 典 →0 for all values of ⌬G ˜ 0. The inclusion of ␺ ˜ (t ˜ ) allows the population difference to relax to a value that depends on ⌬G ˜ 0, only for ⌬G ˜ 0⫽0 do we obtain limt ˜ →⬁ 具 z(t ˜ ) 典 →0. An analytic analysis of the behavior of the equations of motion seems difficult, as the noise term ␩ ˜ ␣ (t ˜ ) is multiplicative. We have not succeeded in finding closed equations for the mean values. III. THE EQUILIBRIUM DENSITY OPERATOR In this section we shall obtain an approximate expression for the systems’s equilibrium density operator that should be correct for a classical bath. Let us stress that, in general, especially when the bath is not classical, there is no rigorous way to disentangle the system’s density operator from that of the bath. Let us assume that the total equilibrium density operator has the canonical form ␳ T eq共 ␤ 兲⫽1 ZT共 ␤ 兲e⫺ ␤ HT,共39兲 where HTis the total Hamiltonian given by Eq. 共2兲and ZT( ␤ ) is the total partition function. Our goal now is to trace out the bath variables by considering them as classical quantities. Hence, in the definition of the reduced density operator of the system at equilibrium, we must replace the trace in the bath Hilbert space, TrB, by an integral over phase space ␳ S eq共 ␤ 兲⫽1 Z共 ␤ 兲TrB关e⫺ ␤ HT兴→1 Z共 ␤ 兲 冕冕 dqdpe⫺ ␤ HT. 共40兲 Integrating over the pvariables, one finds ␳ S eq共 ␤ 兲⫽C共 ␤ 兲 冕冕 dqe⫺ ␤ /2兺jmj ␻ j 2qj 2 ⫻e共 ␤ /2兲((⌬G0⫹2兺j ␥ jqj) ␴ z⫺2បV ␴ x),共41兲 11180 J. Chem. Phys., Vol. 113, No. 24, 22 December 2000 Casado-Pascual et al. where C( ␤ ) is a normalization factor. If we define ␩ ⫽(2/ប)兺j ␥ jqj, then the integration over the qvariables can be easily carried out by inserting an appropriate delta function to obtain ␳ S eq共 ␤ 兲⫽C共 ␤ 兲 冕冕 dqe⫺共 ␤ /2兲兺jmj ␻ j 2qj 2 冕 ⫺⬁ ⬁d ␩ ⫻ ␦ 共 ␩ ⫺共2/ប兲兺j ␥ jqj兲e共 ␤ /2兲(A( ␩ ) ␴ z⫺2បV ␴ x) ⫽C⬘共 ␤ 兲 冕 ⫺⬁ ⬁d ␩ P共 ␩ 兲e共 ␤ /2兲(A( ␩ ) ␴ z⫺2បV ␴ x),共42兲 where C⬘( ␤ ) is another normalization constant, A( ␩ ) ⫽⌬G0⫹ប ␩ , and P( ␩ ) is given in Eq. 共24兲. Using the identity40 e共 ␤ /2兲(A( ␩ ) ␴ z⫺2បV ␴ x)⫽cosh共 ␤ E共 ␩ 兲兲I ⫹sinh共 ␤ E共 ␩ 兲兲 E共 ␩ 兲 冋 A共 ␩ 兲 2 ␴ z⫺បV ␴ x 册 , 共43兲 where E( ␩ )⫽关(A( ␩ )/2)2⫹(Vប)2兴1/2, the normalization constant C⬘( ␤ ) can be expressed in the form C⬘共 ␤ 兲⫽ 冉 2 冕 ⫺⬁ ⬁d ␩ P共 ␩ 兲cosh共 ␤ E共 ␩ 兲兲 冊 ⫺1.共44兲 Hence, if we write ␳ S eq( ␤ )⫽1 2关I⫹Xeq( ␤ ) ␴ x⫹Yeq( ␤ ) ␴ y ⫹Zeq( ␤ ) ␴ z兴, we obtain Xeq共 ␤ 兲⫽⫺បV 兰 ⫺⬁ ⬁d ␩ P共 ␩ 兲sinh共 ␤ E共 ␩ 兲兲/E共 ␩ 兲 兰 ⫺⬁ ⬁d ␩ P共 ␩ 兲cosh共 ␤ E共 ␩ 兲兲 ,共45兲 Yeq共 ␤ 兲⫽0, 共46兲 Zeq共 ␤ 兲⫽ 兰 ⫺⬁ ⬁d ␩ P共 ␩ 兲A共 ␩ 兲sinh共 ␤ E共 ␩ 兲兲/E共 ␩ 兲 2 兰 ⫺⬁ ⬁d ␩ P共 ␩ 兲cosh共 ␤ E共 ␩ 兲兲 ,共47兲 or, in dimensionless form, Xeq共T ˜ 兲⫽⫺V ˜ 兰 ⫺⬁ ⬁d ␩ ˜ P共 ␩ ˜ 兲sinh共E ˜ 共 ␩ ˜ 兲/T ˜ 兲/E ˜ 共 ␩ ˜ 兲 兰 ⫺⬁ ⬁d ␩ ˜ P共 ␩ ˜ 兲cosh共E ˜ 共 ␩ ˜ 兲/T ˜ 兲 ,共48兲 Yeq共T ˜ 兲⫽0, 共49兲 Zeq共T ˜ 兲⫽ 兰 ⫺⬁ ⬁d ␩ ˜ P共 ␩ ˜ 兲A ˜ 共 ␩ ˜ 兲sinh共E ˜ 共 ␩ ˜ 兲/T ˜ 兲/E ˜ 共 ␩ ˜ 兲 2 兰 ⫺⬁ ⬁d ␩ ˜ P共 ␩ ˜ 兲cosh共E ˜ 共 ␩ ˜ 兲/T ˜ 兲,共50兲 where ␩ ˜ ⫽ ␩ /⌬,A ˜ ( ␩ ˜ )⫽⌬G ˜ 0⫹(2E ˜ rT ˜ )1/2 ␩ ˜ ,E ˜ ( ␩ ˜ ) ⫽关(A ˜ ( ␩ ˜ )/2)2⫹(V ˜ )2兴1/2, and P( ␩ ˜ ) is the Gaussian distribution defined in Eq. 共24兲with ⌬⫽1. These expressions can be evaluated numerically quite readily. In the following section, we will compare the equilibrium results with those obtained from the long-time solution of the dynamical equations. IV. RESULTS AND DISCUSSION In this section, we will present numerical results based on our equations of motion that illustrate the competition between electronic coherence and relaxation induced by the system-bath-coupling. Comparisons with other available methods are made, when feasible. Throughout this section we will use the dimensionless form of the equations, Eqs. 共34兲–共38兲, i.e., all frequencies are scaled with respect to the bath frequency ␻ c.If ␻ c⫽50.0 cm⫺1, then a V ˜ in the range 0.01–100.0 corresponds to V⫽0.5–5000.0 cm⫺1, which covers the range of typical electronic couplings. The values of reorganization energies, E ˜ r, used below will be characteristic of nonpolar and polar solvents. For convenience, we will drop the tilde in what follows. Our set of stochastic equations has been integrated using a stochastic Runge–Kutta algorithm of second order in the noise and the deterministic parts.41 For each initial value of ␩ ␣ , drawn from the Gaussian distribution Eq. 共24兲, the equations are solved for a realization of the white noise w(t). Mean values have been obtained with Ntr⫽10000 trajectories. The integration time step has to be matched to the system parameters, especially for large values of Vand ⌬G0. The system is initially prepared with x(0)⫽y(0)⫽0, and z(0)⫽1, corresponding to unit population of state 兩 0 典 关cf. Fig. 共1兲兴. In panel 共a兲of Fig. 2 we display the behavior of the population difference 具 z(t) 典 for ⌬G0⫽0.0, V⫽60.0, Er ⫽80.0, T⫽4.0, corresponding to the strictly adiabatic limit VⰇ ␻ c⫽1.0, and for a nonequilibrium initial condition, ␣ ⫽0.0. This figure shows that the decay is slow and highly oscillatory. The decay should be to zero since ⌬G0⫽0共see below兲. The semiclassical method of Stock35 is equivalent to our method in the limit of an infinite number of bath oscillators. One would, therefore, expect perfect agreement between the two methods for time intervals tⰆ1/⌬ ␻ , where ⌬ ␻ is the separation in frequency between the finite number of oscillators used in Stock’s method. We have integrated Stock’s equations with N⫽400 oscillators with frequencies uniFIG. 2. 共a兲 具 z(t) 典 for V ˜ ⫽60, E ˜ r⫽80.0, T ˜ ⫽4.0, ⌬G ˜ 0⫽0.0 and ␣ ⫽0.0. 共All quantities are given in dimensionless units relative to the bath characteristic frequency, ␻ c.) For this large value of V ˜ , there are fast oscillations whose amplitude slowly decays to zero, the equal population state. 共b兲Comparison of our semiclassical stochastic method 共solid line兲and the semiclassical method of Stock 共plus signs兲. 11181J. Chem. Phys., Vol. 113, No. 24, 22 December 2000 Relaxation in charge-transfer systems formly distributed in the interval 关 ␻ min , ␻ max兴⫽关10⫺3,10兴for a Debye spectral density. The initial conditions for the bath oscillators are obtained from the classical action-angle variables qj共0兲⫽ 冑 2nj⫹1sin共 ␾ j兲, 共51兲 pj共0兲⫽ 冑 2nj⫹1cos共 ␾ j兲, with random phases ␾ jand quantum numbers njdistributed according to the Boltzmann distribution. For the integration of the system of equations in Stock’s model, a fixed time step Runge–Kutta method of second order has been used, integrating Ntr⫽2500 trajectories in order to obtain converged results. In panel 共b兲of Fig. 2, we compare our results with those obtained using the scheme of Stock, which is also designed to be accurate in the adiabatic limit. As expected, the agreement between the two methods is very good in the observed time interval. The difference for larger times is due to the limited number of oscillators in the method of Stock; the agreement may be further improved by increasing this number. The CPU time needed for 1000 trajectories on a Pentium II 共400 MHz兲was 1 min for our method and 290 min for the method of Stock with the abovementioned parameters. This difference stems from the fact that, in Stock’s method, the equations of motion for all Nbath oscillators have to be integrated. This results in a system of 2N⫹4 ordinary first order differential equations, rendering the numerical solution of this problem time consuming for a large number of bath oscillators. On the other hand, the stochastic scheme requires the numerical solution of a set of only five stochastic first order differential equations. Although the necessary CPU time for Stock’s method could be reduced by employing faster integration algorithms such as the velocity Verlet algorithm,42 the stochastic method will always be less numerically demanding. When Eris quite large, damping of 具 z(t) 典 will be in evidence, as shown in Fig. 3. Even though Vis large, the initial decay and oscillation goes over to exponential relaxation. Note that, for the equilibrium initial condition ( ␣ ⫽1.0) used here, the solvent configurations are centered around the reactant well minimum, such that their weight around the crossing point of the electronic surfaces is small, thus de-emphasizing the V-dependent oscillatory behavior of 具 z(t) 典 . The ratio of the slopes of the decays for the two values of Vis about 1.5, indicating that they are not independent of the strength of the electronic coupling. That is, if one were using a Marcus adiabatic ET rate expression,5consonant with a large value of V, the rate ratio would be one, since the reorganization energy is fixed here. In contrast with the equilibrium initial condition used to construct Fig. 3, Fig. 4 shows that if we use the nonequilibrium initial condition ( ␣ ⫽0.0) the oscillations are much more prominent and the decay behavior, while still evident, starts from a much smaller amplitude than for ␣ ⫽1.0. The much larger weight of solvent configurations around the crossing point of the surfaces for this initial condition produces this behavior. We also display in Fig. 4 具 z(t) 典 for the larger Ervalue of 40.0. The decay is similar to the Er⫽20.0 case, with a larger initial amplitude reflecting the increased solvent-solute interaction. The feature of more pronounced oscillations for a nonequilibrium initial condition agrees with the findings of other investigators.24,25 Now we consider the predictions made by our stochastic equations for the long-time average behavior, by comparing their results with the analytic equilibrium values obtained by FIG. 3. 具 z(t) 典 for V ˜ ⫽4.0, E ˜ r⫽20.0, T ˜ ⫽2.0, ⌬G ˜ 0⫽0.0 and ␣ ⫽1.0 共dashed line兲. The initial condition has the bath equilibrated to the reactant state. There is an initial fast decay due to the relatively large value of V ˜ followed by an exponential decay, due to the large reorganization energy, E ˜ r. Since initially there is little solvent population at the crossing point of the surfaces 共see text兲, the crossover to relaxation behavior dominates. For a larger electronic coupling, V ˜ ⫽8.0 共solid line兲, there is more initial oscillation, but 具 z(t) 典 still goes over to an exponential decay, with a decay constant about 1.5 times larger than that of the V ˜ ⫽4.0 case. The circles 共boxes兲 are exponential fits to the data for V ˜ ⫽4.0 (V ˜ ⫽8.0). FIG. 4. 具 z(t) 典 for V ˜ ⫽4.0, E ˜ r⫽20.0 共solid line兲and 40.0 共dashed line兲with T ˜ ⫽2.0, G ˜ 0⫽0.0 and ␣ ⫽0.0. For the nonequilibrium initial condition, the initial solvent population at the surface crossing is much larger than for the data in Fig. 3 and the initial decay is more oscillatory. Nevertheless, 具 z(t) 典 does go over to an exponential decay, with the larger E ˜ rcase starting from a larger value. The boxes 共circles兲are exponential fits to the data for E ˜ r ⫽20.0 (E ˜ r⫽40.0). 11182 J. Chem. Phys., Vol. 113, No. 24, 22 December 2000 Casado-Pascual et al. quadrature from Eqs. 共48兲–共50兲. Note that the stochastic equation formulation is not computationally demanding, which permits us to readily obtain the long-time behavior. In Fig. 5, we compare the analytic and numerical results for values of Vthat are sufficiently large to guarantee the accuracy of our theory. In Fig. 6, the comparison is made as a function of ⌬G0. Clearly, the stochastic equations provide a good account of the equilibrium behavior that is predicted by Eqs. 共48兲–共50兲. However, it must be pointed out that there is no a priori reason that demands coincidence between these results, as we discuss further in Sec. V. The key feature of the stochastic equations responsible for obtaining nonzero 具 z(⬁) 典 values, as must be the case when ⌬G0is not zero, is the dynamical feedback from the system to the bath, as embodied in Eqs. 共29兲–共32兲,共28兲, and 共18兲. If inside the integral defining ␺ (t) in Eqs. 共18兲and 共28兲the value of z(s) were set to one, the resulting equations would lead to 具 z(⬁) 典 ⫽0 for all ⌬G0values. This approximation would amount to neglecting the influence of the dynamics of the solute 共spin兲on the solvent 共bath兲fluctuations. Due to this dynamical coupling, the evolution of the bath cannot be specified by an autonomous stochastic process. Our stochastic equations are appropriate when Vis sufficiently large compared with the other characteristic frequencies; nevertheless, it is interesting to investigate if they give reasonable results for smaller values of V. To carry out this task, we will compare our results with numerical results obtained using the tensor multiplication scheme of Makri and Makarov.33,34 The practical implementation of this scheme requires V⭐ ␻ c, in order to obtain convergence when a reasonable number of nonlocal memory terms are taken into account (⌬kmax⬇10). In what follows, we have chosen parameter values where convergence can be achieved with ⌬kmax⫽9 in the TensMult algorithm, so that the results of the TensMult algorithm can be considered exact. In Fig. 7, we study the behavior of 具 x 典 , 具 y 典 , and 具 z 典 for the unbiased case ⌬G0⫽0 with V⫽0.25 and T⫽2.0 for three different values of the reorganization energy Er ⫽0.2, 0.5, and 1.0. These parameter values correspond to relatively high-temperature and weak system-bath-coupling. For all system-bath-couplings we find very good agreement between the stochastic equation and the TensMult schemes for 具 z 典 . A similar conclusion was reached by Golosov et al.18 in comparing the Stock and TensMult schemes. We have also examined 具 y 典 and find good agreement though the deviations are larger here. However, we see that for 具 x 典 , the FIG. 5. The dependence of 具 x(⬁) 典 , 具 y(⬁) 典 ,and 具 z(⬁) 典 on the tunneling coupling element Vfor fixed T ˜ ⫽4.0, E ˜ r⫽80.0, and ⌬G ˜ 0⫽80.0, as obtained from the numerical solution of the stochastic equations, and the dependence of Xeq,Yeq,Zeq obtained from the equilibrium distribution. All the quantities are in dimensionless form. FIG. 6. The dependence of 具 x(⬁) 典 , 具 y(⬁) 典 and 具 z(⬁) 典 ,on⌬G ˜ 0for fixed T ˜ ⫽4.0, E ˜ r⫽80.0, and V ˜ ⫽60.0, as obtained from the numerical solution of the stochastic equations, and the dependence of Xeq,Yeq and Zeq, obtained from the equilibrium distribution. All the quantities are in dimensionless form. FIG. 7. Comparison of the semiclassical stochastic method 共lines兲and the numerically exact tensor multiplication scheme of Makri and Makarov 共symbols兲.共a兲V ˜ ⫽0.5, E ˜ r⫽0.2, T ˜ ⫽2.0, ⌬G ˜ 0⫽0.0, and ␣ ⫽1.0. 共b兲As in 共a兲but Er⫽0.5. 共c兲As in 共a兲but E ˜ r⫽1.0. The agreement between 具 z(t) 典 and 具 y(t) 典 is excellent for all times, while the 具 x(t) 典 behavior is only good at short times. This deviation reflects the breakdown of our adiabatic approximation as V ˜ gets sufficiently small. Nevertheless, the population behavior is well described by our method. As E ˜ rincreases the 具 z(t) 典 behavior changes from an underdamped to an overdamped oscillatory decay. 11183J. Chem. Phys., Vol. 113, No. 24, 22 December 2000 Relaxation in charge-transfer systems