Activated electron transfer at zero reorganization energy induced by a fluctuating donor–acceptor coupling

 View Online  Export Citation RESEARCH ARTICLE | JUNE 25 2025 Activated electron transfer at zero reorganization energy induced by a fluctuating donor–acceptor coupling Special Collection: Festschrift for Abraham Nitzan Stephanie Valianti ; Spiros S. Skourtis  J. Chem. Phys. 162, 244116 (2025) https://doi.org/10.1063/5.0272917 Articles You May Be Interested In Machine learning to accelerate screening for Marcus reorganization energies J. Chem. Phys. (August 2021) A molecular Debye-Hückel approach to the reorganization energy of electron transfer reactions in an electric cell J. Chem. Phys. (October 2014) Reorganization energy of electron transfer processes in ionic fluids: A molecular Debye-Hückel approach J. Chem. Phys. (March 2013) 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp Activated electron transfer at zero reorganization energy induced by a fluctuating donor–acceptor coupling Cite as: J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 Submitted: 27 March 2025 •Accepted: 29 May 2025 • Published Online: 25 June 2025 Stephanie Valianti and Spiros S. Skourtisa) AFFILIATIONS Department of Physics, University of Cyprus, 1678 Nicosia, Cyprus Note: This paper is part of the JCP Special Topic, Festschrift for Abraham Nitzan. a)Author to whom correspondence should be addressed: [email protected] ABSTRACT Westudyanalyticallyandnumericallyaratemodelfordonor-to-acceptorchargetransferinthecaseofafluctuatingdonor–acceptorcoupling and of zero reorganization energy. The model applies to situations where the donor and acceptor reorganization energies are very low such that there is no polaron formation and the donor-to-acceptor transition arises purely from dynamic disorder in the coupling. We consider both quantum and classical limits for the reaction coordinate that modulates the coupling and describe the adiabaticity parameter for chargetransfer induced by coupling fluctuations, analogous to the Landau–Zener parameter used in the case of transfer induced by energy-level fluctuations. Our purpose is to explore the magnitudes of the charge-transfer rate for realistic parameter values, the transition from nonadiabatic to adiabatic transfer, and the behavior of the rate as a function of the donor–acceptor energy gap and temperature. We find that the coupling-fluctuation mechanism can lead to fast rates. Furthermore, the energy-gap and temperature dependencies can differ from Marcus theory describing polaron transfer. Therefore, these dependencies may be used to identify experimentally whether charge transfer is induced by coupling fluctuations. ©2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercialNoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/). https://doi.org/10.1063/5.0272917 I. INTRODUCTION Molecular electron-transfer reactions are ubiquitous in biological, chemical, and molecular-junction systems.1–4 A theoretical framework that is most often used to describe molecular electron transfer is Marcus theory and its extensions. The theory applies to thephysical situation wheretransfer involves localizedpolarons that exist even at zero donor(D)–acceptor(A) free-energy gap, due to an inner-sphere reorganization energy that is sufficiently large to localize the electronic wave function. The earlier form of Marcus theory consideredchargetransferinducedbyelectronicenergy-levelfluctuations.4It was later recognized that in addition to D(A) energy-level fluctuations, D–A coupling fluctuations can also enhance the transfer rate. To treat these cases, non-adiabatic Marcus theory was modified by introducing a coordinate-dependent D–A coupling. The fluctuations of the coupling coordinates widen the region of D–A resonance, introducing two types of thermal activation arising from energy-level and coupling fluctuations. Such effects are often labeled non-Condon in the context of non-adiabatic electron transfer. Early analytical work (e.g., Refs. 5–10) and computational work that includes coupling fluctuations in rate simulations (e.g., Refs. 11–16) have led to extensive research activity on this subject over the years, as described in several reviews.17–21 In most cases, D(A) energy-level fluctuations to resonance are the dominant mechanism inducing transfer, even if couplingfluctuations are present. This happens especially when the D–A free-energy gap ΔGis the largest energy parameter, such that the thermalactivationtotheD–Aresonanceconformationistheslowest (rate-limiting) process. In cases with zero free-energy gap, the activation energy to D–A resonance is determined by the inner-sphere reorganization energy λand the D–A coupling V(Eact =∣V∣−λ/4). For polymeric systems with multiple and identical D(A) units, Marcus theory is not applicable for describing the hopping steps in the limit of ∣V∣≃λ/4. The reason is that the charge-carrier wave J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-1 © Author(s) 2025 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 1. Types of D–A charge-transfer systems considered in this work are comprised of identical D and A moieties that involve rigid delocalized aromatic structures, giving rise to negligible inner-sphere reorganization energies. D-to-A charge transfer at zero energy gap is induced by a relative-coordinate motion Q that transiently enhances the D-to-A coupling for a limited range of Qvalues that are accessed by thermal activation. Two examples of the reaction coordinate: (a) Qis a torsional coordinate and (b) Qis a slipping coordinate. functions may delocalize over multiple monomers within the polymer. Other theoretical approaches have been developed for such situations, including flickering-resonance theory for molecular electron transfer22,23 and transient delocalization theories for electron transfer in crystals of organic semiconductors.24,25 For a D–A dimer system with fluctuating D–A coupling, the limit of low reorganization energy has been studied in the context of non-adiabatic (Fermi’s golden rule) rate theory that includes non-Condon effects. It was found that the behavior of the rate as a function of energy gap in this regime can be very different from non-adiabatic Marcus theory.26 We revisit this type of system and consider the case of D-to-A electron transfer at zero inner-sphere reorganization energy, where the transfer is induced by coupling fluctuations. Figure 1 describes the generic characteristics of the types of systems relevant to our work. They are comprised of identical D and A moieties that are large and rigid with delocalized electronic wave functions, thus having negligible inner-sphere reorganization energies. Examples of such moieties relevant to both artificial and biological charge-transfer systems include the longer acenes27 and planar metal-organic systems such as metal-porphyrin fused rings28 and nickel-dithiolene fused rings.29 The reaction coordinate that modulates the coupling (Q)may be a relative-motion coordinate such as a slippage displacement that enhances/reduces pi-stacking30,31 or a torsional angle, as in oligophenylene wire systems32–34 and edge-fused nickel-porphyrin nanoribbons.35 Qmay have quantum or semiclassical/classical dynamics, and we consider all of these cases. For the quantum case, we use a density-matrix approach36–38 and explore regimes that go beyond the non-adiabatic limit described by Fermi’s golden rule. We do not assume an approximate linear dependence of the D–A coupling on Qthat is often used in the literature, as we obtain exact analytical expressions for the vibrationally modulated coupling within the nonlinear model of Fig. 2. In addition, we describe the energy-gap dependence of the coupling-fluctuation rate for quantum dynamics of Qin different regimes of the quantum dynamics and show FIG. 2. Top: Harmonic oscillator potential for reaction coordinate Qthat modulates the D–A coupling. Eact is the activation energy to reach the value of Qthat maximizes the D–A coupling. hω(nact +1/2)is the harmonic oscillator eigenenergy closer to Eact.hω(nRT +1/2)is the harmonic oscillator eigenenergy closer to kBT300 K =25.0 meV(nRT ≡n300 K). Bottom: D–A coupling as a function of Q, showing Qmax and δof Eq. (6). that the energy-gap dependence deviates from that of non-adiabatic quantum Marcus theory. For classical/semiclassical dynamics of Q, we develop a transition-state-theory (TST) approach that enables the analysis of the transition from the non-adiabatic to the adiabatic regimes, analogous to the analysis often used for classical Marcus theory. Our approach leads to the derivation of the adiabaticity parameter for the coupling-fluctuation model and to the prediction of a different temperature dependence of the rate in the non-adiabatic limit, as compared to Marcus theory. Although the theory developed relates to charge transfer in a D–A dimer, it is also applicable to linear polymeric systems in solution or in a molecular matrix and in molecular junctions. The requirement is that transport can be described by a Markov chain comprised of forward and backward rates between nearest-neighbor D–A monomers. II. FLUCTUATING-COUPLING MODEL WITH ZERO REORGANIZATION ENERGY The Hamiltonian of the model is given by ˆ H=ˆ Hel ˆ Ivib +ˆ Iel ˆ Hvib +ˆ Hel−vib,(1) where ˆ Hel =E0∣D⟩⟨D∣+(E0+ΔG)∣A⟩⟨A∣,(2) ˆ Hel−vib =V(Q)(∣D⟩⟨A∣+h.c.),(3) J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-2 © Author(s) 2025 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp and ˆ Hvib =1 2μˆ P2+1 2μω2ˆ Q2 =∞ ∑ n=0εn∣n⟩⟨n∣,εn=hω(n+1/2).(4) E0(E0+ΔG)is the time-independent D(A)-state energy and V(Q) is a D–A fluctuating coupling. Qis a reaction coordinate that modulates the coupling. We assume that Qis a harmonic oscillator coordinate described by the vibrational Hamiltonian ˆ Hvib, where ω is the vibrational frequency and μis an effective mass. This model describes a situation where the D and A energy fluctuations are not significant on the timescale of charge transfer as compared to the coupling fluctuations, and there is not significant reorganization energythatcouldinducepolaronlocalizationatΔG=0.Wewill first assume that the dynamics of Qare quantum and develop a theory for the charge-transfer rate. Then, we will consider an approximate transition-state theory (TST) applicable to the case of classical and semiclassical underdamped dynamics of Q. Suppose that V(Q)peaks at a value of Q=Qmax,V(Qmax) =Vmax =V0, and it is non-zero in a region ΔQ:Qmax −δ<Q <Qmax +δ. This is a very generic behavior. Semiclassically, it leads to thermally activated charge transfer if the rms position of the harmonic oscillator ˜ σQ(T,ω)2=σ2 0coth(hω/2kBT),σ2 0=h 2μω,(5) is not in the region ΔQat lower temperatures (Tlow)and enters the region at higher temperatures (Thigh). That is, Qmax −δ >˜ σQ(Tlow,ω)and Qmax −δ<˜ σQ(Thigh,ω). However, this intuitive image can be inaccurate in the deep quantum regime (hω≥kBT). In the following, we first consider a purely quantum theory for the transfer rate and then discuss classical and semiclassical transitionstate theory approximations. In our analysis, we use a Gaussian model for the D–A coupling, i.e., V(Q)=V0e−(Q−Qmax)2/2δ2,(6) (see Fig. 2). The Gaussian model is a fair representation of the nonzero coupling region described earlier that is characterized by a maximum coupling V0and a width δ. A. Quantum theory for Q The Hamiltonian is expressed in the basis of vibronic states ∣i,n⟩=∣i⟩el∣n⟩vib, where ndenotes the nth eigenstate of ˆ Hvib. In this representation, ⟨i,n∣ˆ H∣j,m⟩=Eel iδi,jδn,m+(1−δi,j)VDn,Am,(7) where Eel i=E0,(E0+ΔG)for i=D,(A), VDn,Am =∫dQ ψ∗ n(Q)V(Q)ψm(Q),(8) where ψn(Q)is the wave function for a harmonic oscillator eigenstate, ψn(Q)=(1 2πσ2 0)1/4(1 2nn!)1/2e−Q2/4σ2 0[Hn(1 √2σ0Q)].(9) σ2 0=h/(2μω)is the ground-state rms fluctuation of Qand Hnis a Hermite polynomial of order n. Evaluating this matrix element for the Gaussian V(Q), we obtain VDn,Am =VD0,A0×min(m,n) ∑ k=0fn,m,k(σ0,δ,Qmax),(10) where VD0,A0=∫dQ V(Q)∣ψn=0(Q)∣2 =V0[δ Σ]e−Q2 max/2Σ2,Σ2=δ2+σ2 0,(11) is the vibronic coupling for the ground vibrational state, and fn,m,k(σ0,δ,Qmax)=√m!√n! (m−k)!(n−k)!k![σ0 √2Σ]m+n−2k ×[Hm+n−2k(Qmax √2Σ)].(12) At fixed δ(fixed Σ), as a function of Qmax, the coupling VDn,Am is enhanced with respect to VD0,A0when Qmax >√2Σ since Hm+n−2k(Qmax/√2Σ)∝(Qmax/√2Σ)m+n−2kin the abovementioned equation. Increasing δreduces this enhancement, since (σ0/√2Σ)m+n−2k≪1. Examples of the behavior of VDn,Am as a function of these parameters are given in Figs. 8 and 10. To simulate the D-to-A transition probability, we consider the reduced density matrix of the system, ˜ ρ(t), represented in the electronic-vibrational basis (matrix elements ρi,n;j,m). If we include Noscillator states in our computations, ˜ ρ(t)is a 2N×2Nmatrix, which is written in Liouville space as a (2N)2×1 vector ρ(t). The time evolution of ρ(t)is governed by the equation d dt ρ(t)=−i h˜ Lρ(t),(13) where ˜ Lis the (2N)2×(2N)2Liouvillian matrix, consisting of a coherent and a dissipative component, ˜ L=˜ Lcoh +˜ Ldiss.(14) The matrix elements of the coherent Hamiltonian are given by Lcoh in,jm;kn′,lm′=Hin,kn′δjm,lm′−Hjm,lm′δin,kn′,(15) i.e., they describe the dynamics arising from [ˆ H,ˆ ρ]. The dissipative part of the Liouvillian arises from tracing out the bath (environmental) degrees of freedom of the system–bath Hamiltonian when describing the time evolution of the system variables. The effect of the environmental (solvent environment) degrees of freedom J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-3 © Author(s) 2025 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp is to modulate the system vibronic energy levels (pure dephasing rates) and to induce vibrational relaxation and activation rates (population exchange) between the system vibronic states. The latter rates describe transitions that involve the exchange of energy between the system and the environment. In summary, the matrix elements of the dissipative Liouvillian contain population transfer rates, Γin←im (n≠mand i=Dor A), and population relaxation rates Γtot in =∑n≠mΓin←im. In particular, the relevant non-zero matrix elements of ˜ Ldiss are Ldiss in,in;im,im =ihΓin←im (n≠m),Ldiss in,in;in,in =−ihΓtot in ,(16) with Γβ←α/Γα←β=exp(−(Eβ−Eα)/kBT). The rates satisfy detailed balance such that the steady-state probabilities of the vibronic states obey Boltzmann statistics. The non-zero matrix elements describing the dissipative time evolution of coherences are Ldiss in,jn;in,jn =−ih(Γtot in /2+Γtot jn /2+γ(1−δi,j)),(17) where γdenotes a pure electronic-dephasing rate between the electronic D and A states. The electronic dephasing rates are the fastest ratesintheLiouvillian,andinoursimulations,wechoosetypicalvalues of hγ≤7×10−3eV (timescales of hundreds of fsec or slower). If we associate the energy-level fluctuations with an outer-sphere reorganization energy of the environment, i.e., hγ≃√λoutkBT, then the chosen values of the pure dephasing rates correspond to a very low outer-sphere reorganization energy at room temperature, λout ≤10−3eV. Equation (13) is solved using an initial density matrix vector ρ(t=0)describing a system in the D electronic state where the vibrational states of Qare equilibrated, i.e., the non-zero-matrix elements of ρ(0)are ρDn,Dn(0)=e−εn/kBT/Zvib,(18) where εn=hω(n+1/2)and Zvib =∑ne−εn/kBT. The D-to-A transition probability as a function of time is given by PD→A(t)=∑ nρAn,An(t).(19) The estimation of the thermally averaged D-to-A rate is discussed in Sec.III.InTableI,wesummarizethedifferentlimitsofthequantum model that are explored in this work. B. Different limits of PD→A(t) We now consider regimes for which we can derive quasianalyticalexpressionsforPD→A(t).WhenΔG=0(elastictransport), TABLE I. Coupling-fluctuation rate model for quantum dynamics of the reaction coordinate Q. Regimes of the model under study. Strongly coupled: Weakly coupled: hω>hγ>hΓtot n2VDn,An >hγ/2∀n2VDn,An <hγ/2∀n [Figs. 4(b) and 4(c)] [Fig. 4(a)] Strongly coupled: Weakly coupled: hγ>hω>hΓtot n2VDn,An >hγ/2∀n2VDn,An <hγ/2∀n [Figs. 5(b) and 5(c)] [Fig. 5(a)] each ∣Dn⟩is resonant to ∣An⟩and coupled via the matrix element VDn,An. If, in addition, the off-resonant matrix elements VDn,Am of each∣Dn⟩to∣An ≠m⟩are weak, ∣VDn,Am∣<(n−m)hωforn≠m,the exactPD→A(t)canbeapproximatedasasumofindependentchargetransport channels Pchan Dn→An(t)that involve only the diagonal VDn,An elements. In particular, for the case of slow vibrational relaxation with respect to electronic dephasing, hΓtot n≪hγ, we obtain PD→A(t)≃Pind−chan D→A(t)=∞ ∑ n=0 1 Zvib e−εn/kBTPchan Dn→An(t).(20) The analytical form of Pchan Dn→An(t)depends on the relative magnitudes of γ/2 and VDn,An. If 2VDn,An <hγ/2, Pchan Dn→An(t)=Pchan(I) Dn→An(t), where Pchan(I) Dn→An(t)≃1 2(1−e−γt/2[hγ/2 2˜ Ansinh(˜ An ht)+cosh(˜ An ht)]), (21) with ˜ An=√(hγ/2)2−(2VDn,An)2. If 2VDn,An >hγ/2, Pchan Dn→An(t) =Pchan(II) Dn→An (t), where Pchan(II) Dn→An (t)≃1 2(1−e−γt/2[hγ/2 2Ansin(An ht)+cos(An ht)]),(22) with An=√(2VDn,An)2−(hγ/2)2. For the regimes explored in Sec. III, these approximate forms can be quite accurate (see the supplementary material, Sec. I). The case 2VDn,An <hγ/2 gives slower charge transfer as compared to the case 2VDn,An >hγ/2. It is possible that at lower temperatures (low n), VDn,An <hγ/2, and at higher temperatures the thermally accessible higher −noscillator eigenstatesgiveVDn,An >hγ/2, leading toalargeenhancement of the rate. For the inelastic case of ΔG≠0, channels involving VDn,Am (n≠m)become important as they conserve energy, and the abovementioned approximation fails (see Fig. 3). For example, when ΔG=−j(hω), where jis an integer, inelastic channels involving VDn,An+jbecome important. Finally, the quantum non-adiabatic regime arises when hω>hγ,hΓtot n, and in addition each ∣Dn⟩is FIG. 3. Schematic diagram of the parameters of the quantum coupling-fluctuation model for the charge-transfer rate, in the basis of vibronic states ∣Dn⟩and ∣An⟩. Thediagram does notshow thepure electronic dephasing rate γ.Three cases are considered, ΔG=0, ΔG<0, and ΔG>0. J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-4 © Author(s) 2025 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp weakly and off-resonantly coupled to all final ∣Am⟩. In this situation, the time evolution of each Dn →Am channel is described by (2VDn,Am/(εn−εm−ΔG))2sin2((εn−εm−ΔG)t/2h). This behavior leads to the standard Fermi’s golden-rule rate dependence of V2 Dn,Amδ(εn−εm−ΔG)39 that is identified with the quantum non-adiabatic limit. The inelastic case is analyzed in Sec. III. C. Transition state theory for classical and semiclassical dynamics of Q We first use classical transition-state theory (TST) to derive the thermally averaged D-to-A rate. Our derivation follows the same logic as the derivation of classical Marcus theory from TST and the Landau–Zener (LZ) probability. This approach clarifies the meanings of the classical-dynamics adiabatic and non-adiabatic regimes for the coupling-fluctuation model. Our aim is to compare classical Marcus theory for pure polaron transport, i.e., the case of ΔG=0 and λ≠0, to the analogous theory for coupling-fluctuation induced transportintheabsenceofpolaronformation,i.e.,thecaseofΔG=0 and λ=0, where λis the inner-sphere reorganization energy. We then derive the classical TST from the quantum model. The results are summarized in Table II. Finally, we discuss a semiclassical TST for the rate. V(Q)isgivenbyEq. (6), but itisassumedthatQobeysclassical dynamics, governed by the Hamiltonian Hvib(Q,˙ Q)=μ 2˙ Q2+μω2 2Q2.(23) To apply TST as used to derive Marcus theory, we assume underdamped oscillations of Qso that the non-zero coupling region is localized, in the sense that δ≪Qmax (see Fig. 2). Therefore, V(Q)is significant only in a small region around Qmax ±δ, which is considered the transition-state (TS) region, i.e., Qmax =QTS. The activation energy for reaching the TS region is Eact =μω2Q2 max/2. Classical trajectories Q(t)with Etot =Hvib(Q,˙ Q)<Eact do not reach QTS =Qmax and do not contribute to charge transfer since V(Q(t))≃0. Trajectories with Etot =Hvib(Q,˙ Q)>Eact cross QTS twice per oscillator period,inducingchargetransfersinceV(Q(t))≠0intheTSregion. The transfer probability from D to A at every crossing of QTS can be approximatedbyadoptingaconstantvelocityapproximationforthe trajectory near QTS, i.e., Q(t)−QTS =v×(t−tTS),(24) TABLE II. Top: The regime of validity of classical (ballistic) TST for the couplingfluctuation model in terms of the parameters of the quantum model (see also Table I). Bottom: Adiabatic and non-adiabatic limits of classical TST in terms of the adiabaticity parameter of the coupling-fluctuation model. σvis the rms velocity of the reaction coordinate Q.V0and δare the maximum coupling and the width of the coupling region, respectively. kBT≫hω>hΓtot n/2, 2VDn,An >hγ/2, hΓtot n/2 Adiabatic limit: Non-adiabatic limit: 1<(V0  h√πδ σv)2(V0  h√πδ σv)2<10−1 where tTS is a time for which Q=QTS, and v=[˙ Q(t)]t=tTS is the velocity at Q=QTS,v=hstep(E−Eact)×√(2/μ)(Etot −Eact). Rewriting V(Q)as a function of time, we have V(Q(t))=V0e−(t−tTS)2/2σ2 t,(25) where σt=δ/v. The D-to-A transition probability for every complete crossing of QTS is PD→A(v)=sin2(1 h∫∞ −∞ dt′V(t′)) =sin2(√2πV0 hσt)=sin2(√2πV0 h(δ v)).(26) Using standard classical TST,38,40 the thermally averaged D-to-A rate is ⟨kDA⟩cl TST =2∫∞ 0ρcl(QTS,v)v PD→A(v)dv,(27) where ρcl(Q,v)=e−Hvib(Q,v)/kBT/Zcl vib =ρcl(Q)ρcl(v), with ρcl(Q)=1 √2πσQe−Q2/2σ2 Q,ρcl(v)=1 √2πσve−v2/2σ2 v,(28) and σQ=¿ Á Á ÀkBT μω2,σv=√kBT μ.(29) σQis the high-temperature limit (hω/2≪kBT)of ˜ σQin Eq. (5). Therefore, the thermally averaged D-to-A ET rate becomes ⟨kDA⟩cl TST =1 πσQσve−Q2 TS/2σ2 Q∫∞ 0e−v2/2σ2 vv PD→A(v)dv ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ I(v) ,(30) where PD→A(v)=sin2(b/v),b=√2πV0 hδ.(31) The exponential factor is the Arrhenius factor for activation to the transition state, e−Q2 TS/2σ2 Q=e−Eact/kBT. 1. Non-adiabatic and adiabatic contributions to the rate Analogous to the LZ argument, the non-adiabatic limit is reached when the transition state QTS is crossed very rapidly (at high velocity vsuch that b/v≪1). The non-adiabatic limit is when PD→A(v)can be approximated well by its lowest-order Taylor expansion with respect to b/v, i.e., PD→A(v)≃(b/v)2. To make the distinctionbetweenadiabaticandnon-adiabaticlimitsmoreprecise, we rewrite the rate in Eq. (30) in terms of a dimensionless velocity parameter, x=v/b, ⟨kDA⟩cl TST =1 πσQσve−Q2 TS/2σ2 Qb2I(c) =μω kBTe−Eact/kBT2V2 0 h2δ2I(c),(32) J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-5 © Author(s) 2025 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp where I(c)=∫∞ 0e−c2x2x PD→A(x)dx =∫∞ 0e−c2x2xsin2(1/x)dx.(33) The parameter ccan be written as c=b √2σv=τcoup τRabi where τRabi ≡πh V0,τcoup =δcoup σv,(34) where δcoup =π3/2δand τRabi is the Rabi time for a complete D-to-A transitionat D–A resonance. σvisthe thermal rms velocityof Q, and δis a measure of the Qregion that gives non-zero coupling (see also Table II). Therefore, τcoup can be interpreted as the time spent in the non-zero coupling region. The non-adiabatic contribution to I(c)is defined here as the part of the integral (the xregion) for which the lowest-order Taylor expansion of sin2(1/x), i.e., (1/x)2, gives 90% or more of sin2(1/x). This holds for xc<x<∞where xc=0.8. Therefore, we define the non-adiabatic and the remaining contributions as follows: I(c)=∫xc 0e−c2x2xsin2(1/x)dx ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ ˜ I(c)+∫∞ xce−c2x2xsin2(1/x)dx ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Inad(c) . (35) The rate is considered fully non-adiabatic when I(c)≃Inad(c). For c2<10−1, I(c)≈Inad(c)=∫∞ 0.8 e−c2x21 xdx =−1 2Ei(−(0.8c)2)+h.c., (36) where Ei(z)is the exponential integral. We can expand the exponential integral in the region as Ei(−z)=γ′+ln(z)+∞ ∑ k=1(−z)k kk!,(37) where γ′=0.57721 is the Euler–Mascheroni constant.41 SubstitutingEqs.(36)and(37)intotherateexpressionofEq.(32)andkeeping terms to first order in z, we obtain ⟨kDA⟩nad(cl) TST =μω kBTe−Eact/kBTV2 0 h2δ2 ×(−γ′−ln(0.64μb2 2kBT)+0.64μb2 2kBT).(38) Therefore,thenon-adiabaticratehasatemperaturedependencethat is very different from that of classical Marcus theory (see below). We also find numerically that for 1 <c2,Inad(c)makes a negligible contribution to the sum I(c)=˜ I(c)+Inad(c), and ˜ I(c)=∫0.8 0e−c2x2xsin2(1 x)dx ≈1 4c2.(39) Substituting into the rate equation [Eq. (32)], we obtain the adiabatic-limit rate, ⟨kDA⟩ad(cl) TST =ω 2πe−Eact/kBT.(40) Therefore, the purely adiabatic and non-adiabatic regimes appear for 1 <c2and c2<10−1, respectively (see Table II). In the region 10−1<c2<1, the rate is in an intermediate regime. 2. Comparison to the Landau–Zener model and classical Marcus theory In contrast to Eqs. (2) and (3), the Hamiltonian for deriving the classical non-adiabatic and adiabatic Marcus theory considers fluctuating D and A energies and a constant coupling, i.e., ˆ Hel +ˆ Hel−vib =ED(Q)∣D⟩⟨D∣+EA(Q)∣A⟩⟨A∣ +V0(∣A⟩⟨D∣+h.c.).(41) Qis a classical harmonic oscillator coordinate whose equilibrium position changes upon the change of electronic state, i.e., QA=QD+δQand ED(Q)=μω2(Q−QD)2/2, EA(Q)=μω2 (Q−QA)2/2+ΔG. For the transition state, Q=QTS is the value of Qat which there is D–A resonance, ED(QTS)=EA(QTS)=Eres. The D-to-A transition probability is obtained from the Landau–Zener (LZ) model by setting ED(A)(Q)=Eres +[dED(A)(Q)/dQ]QTS (Q−QTS)in ˆ Hel(Q)and assuming constant velocity motion, Q(t)=v×(t−tTS), where v=[˙ Q]Q=QTS (transforming the electronic Hamiltonian to a time-dependent Hamiltonian, ˆ Hel(Q(t)). The D-to-A transition probability is2,38,40 PLZ D→A(v)=1−exp(−˜ b/v),˜ b=2πV2 0 h√2μω2λ,(42) and the thermally averaged D-to-A rate according to TST is computed by substituting PLZ D→A(v)for PD→A(v)into Eq. (27). Following the same steps as before, we transform the integral over the dimensionless velocity x=v/˜ b. The result [analogous to Eq. (32)] is ⟨kLZ DA⟩cl TST =1 πσQσve−Q2 TS/2σ2 Q˜ b2ILZ(˜ c),(43) where ILZ(˜ c)=∫∞ 0e−˜ c2x2x PLZ D→A(x)dx =∫∞ 0e−˜ c2x2x(1−exp(−1/x))dx,(44) and ˜ c=˜ b √2σv=τLZ τRabi where τRabi =πh V0,τLZ =δLZ σv.(45) δLZ =π2V0/√μω2λis the length of the LZ region and τLZ is the time spent in this region [the analogous quantities for the fluctuating-coupling model are δcoup and τcoup in Eq. (34)]. In this model, the non-adiabatic contribution to ILZ is defined as coming from the region of xfor which J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-6 © Author(s) 2025 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 4. D-to-A transition probability as a function of time for ΔG=0 at high (T=300 K)and low temperatures (T=35 K)for (a) V0=5.0 meV, (b) V0=10.0 meV, and (c) V0=0.1 eV. For (a) 2VDn,An <hγ/2 for all thermally occupied n. For (b) and (c), 2VDn,An >hγ/2 for all thermally occupied n. Parameter values: N=20, hω=18.0 meV, γ=(0.1 ps)−1,Γn−1←n=(1.0 ps)−1, Qmax =7.0○,σ0=2.9○, and δ=2.175○. FIG. 5. D-to-A transition probability for ΔG=0 as a function of time at high (T=300 K)and low temperatures (T=35 K)for (a) V0=5.0 meV, (b) V0=10.0 meV, and (c) V0=0.1 eV. For (a) and (b), 2VDn,An <hγ/2 for all thermallyoccupiedn.For(c),2VDn,An >hγ/2forn≥3.Parametervalues:N=20, hω=6.0 meV, γ=(0.1 ps)−1,Γn−1←n=(1.0 ps)−1,Qmax =8.5○,σ0=2.0○, and δ=0.5○. J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-7 © Author(s) 2025 25 June 2025 10:49:30 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp PLZ D→A(x)≃1/x. Following the same steps as before, we can write ILZ(˜ c)=˜ ILZ(˜ c)+Inad LZ (˜ c). Using the classical expressions for σQand σv[see Eq. (29)], we find that for ˜ c≪1 (nonadiabatic limit), ILZ(˜ c)≃Inad LZ (˜ c)=√π/2˜ c. Substituting into Eq. (43) gives the classical non-adiabatic Marcus-rate expression ⟨kLZ DA⟩nad(cl) TST =(2π/h)(V2 0/√4πλkBT)×exp(−Eact/kBT). For ˜ c>1 (adiabatic limit), ILZ(˜ c)≃1/2˜ c2and substitution into Eq. (43), leads to the adiabatic rate expression ⟨kLZ DA⟩ad(cl) TST =(ω/2π) ×exp(−Eact/kBT). To summarize, we compare the important parameters and parameter relations of classical Marcus theory to those of classical TST for the coupling-fluctuation model (both for ΔG=0). In the case of Marcus’s theory, these are ∣V0∣<λ/4, Eact =λ/4−∣V0∣, ˜ c=τLZ/τRabi,τLZ =δLZ/σv, δLZ =π2V0/√μω2λ, (46) where σv=√kBT/μ,τRabi =πh/V0, and ˜ cis the adiabaticity parameter. The relation ∣V0∣<λ/4 ensures that a localized D or A polaron image is valid. For the coupling-fluctuation model, the corresponding parameters are Qmax −δ<σQ,Eact =μω2Q2 max/2, c=τcoup/τRabi,τcoup =δcoup/σv, δcoup =π3/2δ, (47) where cis the adiabaticity parameter and δcoup is the effective width of the non-zero coupling region. The condition Qmax −δ<σQ(T) ensures that at temperature T, the localized D and A state image is valid, in the sense that the D–A coupling is small. Observe that there is no equivalent to λin the case of coupling fluctuations. In both theories the adiabatic limit (˜ c, or c>1) is (ω/2π)e−Eact/kBT. However, the non-adiabatic limit (˜ c, or c<10−1) differs. In the case of Marcus theory, the pre-exponential factor is proportional to T−1/2, whereas for the fluctuating-coupling model, the pre-exponential factor is proportional to T−1×(const1+const2ln(T−1)+const3T−1). 3. Derivation of the classical transition-state rate from the quantum model In order to understand the regime of validity of the classical TST results for the fluctuating-coupling model, it is necessary to derive them from quantum theory. To this end, we rewrite Eq. (26) as PD→A(v)=sin2(1 h∫∞ −∞ dt′V(t′))=PD→A(E), where PD→A(E)=sin2(1 h∫∞ −∞ dQ 1 v(E,Q)V(Q)), v(E,Q)=dQ(t)/dt =√(2/μ)(Etot −Eact). (48) Equation (27) is also rewritten as ⟨kDA⟩cl TST =∫∞ Eact dE e−E/kBT Zcl vib kDA(E), where kDA(E)=2PD→A(E) Tper =ω πPD→A(E).(49) The energy-resolved rate kDA(E)is twice the D-to-A transition probability for a single QTS crossing, divided by the oscillator period (the factor of 2 arises because there are two crossings per period for underdamped harmonic motion). Suppose that the quantum thermally averaged probability PD→A(t)is well approximated by Pind−chan D→A(t)in Eq. (20). If, in addition, 2VDn,An >hγ/2, Pind−chan D→A(t)is the thermal average of FIG. 6. Equilibrium occupation probability for the nth vibrational state as a function of nforhω=6.0 meV and (a) T=35 K and (b) T=300 K. J. Chem. Phys. 162, 244116 (2025); doi: 10.1063/5.0272917 162, 244116-8 © Author(s) 2025 25 June 2025 10:49:30