Semiclassical treatments of electron transfer rate from weak to strong electronic coupling regime

Abstract

Electron transfer (ET) rate is a fundamental parameter to characterize ET processes in physical, chemical, material and biologic sciences. It is affected by a number of quantum phenomena, such as nuclear tunneling, curve crossing, quantum interference, and the coupling to the environment. It is thus a challenge to accurately evaluate the ET rate since one has to incorporate both quantum effects and dissipation. In this review article, we present several semiclassical theories proposed in our group to cover the regime from weak to strong electronic coupling. Their applications to some concrete systems are also shown.

Keywords

electron transfer (ET) / solvent dynamic effect / electronic coupling / theoretical study

Introduction

In electron transfer (ET) process, two types of ETs, nonadiabatic and adiabatic reactions, are often named depending on the weak and strong electronic (diabatic) coupling strengths between the donor and acceptor states. In the early treatments of the ET rate, the environmental fluctuation is assumed to be faster than the ET dynamics and the thermal equilibrium distribution is maintained in the donor state. Under this assumption, the nonadiabatic and adiabatic ET is relatively well understood. In the weak electronic coupling regime, for instance, the reactions can be well treated by the Fermi Golden rule [1,2] or the well-known Marcus formula [3,4] at high temperatures. In the strong electronic coupling limit, on the other hand, the transition state theory (TST) [4,5] is applicable in the Marcus normal region. To cover the weak and strong coupling limits, Marcus and Hush [4,6] have used the thermally averaged Landau-Zener (LZ) transmission probability [7–9] to perform the Marcus-Hush formula. The nuclear tunneling effect is additionally incorporated independently from the electronic transition.

Although€ those€ approaches€ have€ been€ successfully applied to explain ET in amount systems, their limitations are also obvious, especially in the intermediate coupling regime. For instance, Marcus-Hush formula becomes invalid to explain ET at low temperatures because the electronic and nuclear tunneling are coupled together and the simple multiplication form cannot represent the effects correctly [10,11]. Also, the tunneling corrections frequently used come from a parabolic approximation [12,13] which is divergent at low temperatures. The transition should be treated properly as a nonadiabatic transition as a whole. Besides, the LZ formula does not work at all at energies close to and lower than the potential crossing. Furthermore, if the solvent relaxation is very slow compared with the ET process, multiple crossings of the transition region become possible even in the weak coupling regime. Ultimately, the reactions can become independent of the electronic coupling, i.e., they become solvent controlled adiabatic reactions, although they are still nonadiabatic in the absence of solvent dynamics.

Based on the original works independently done by Zusman [14] and Burshtein and coworkers [15], a large number of theoretical approaches have been proposed to treat the competition between solvent relaxation and electronic transition (see for instance, Ref. [16]. and references therein). Yet, the true adiabatic process is different from the nonadiabatic one with slow solvent polarization modes. Several works have clarified the differences [17–19].

When the intramolecular vibrational modes are incorporated in ET one has to deal with multi-dimensional diffusive equations because the motions of these modes may be also damped by the surroundings. Under the assumption that the fluctuations of the intramolecular vibrational modes are much faster than that of the solvent mode, the multi-dimensional equations could be reduced to a one-dimensional one with a sink function, as described by the Sumi-Marcus (SM) theory [20]. In the original SM theory [20] the sink function is dealt with in the classical limit of the local Golden rule. To incorporate the quantum effect, Bixon and Jortner (BJ) [21] proposed an alternative approach to treat the high-frequency vibrational modes. The ideas of SM and BJ have been subsequently combined in a hybrid model by Barbara and coworkers [22]. The hybrid model is a minimal model which envisages ET to occur on a three-dimensional surface with respect to a solvent polarization coordinate, a low frequency classical vibrational and a high frequency quantum mechanical vibrational ones.

Note that the alternative studies [23–38] based on the multi-dimensional dynamical approaches have also been carried out on considering the electronic population relaxation and its transient non-equilibrium properties. These approaches are useful for the dynamical investigation, but when the intramolecular modes are explicitly incorporated their numerical convergence for the rate calculation is still a dilemma [39].

In the present work, we start from ET in which the donor population has a thermal equilibrium distribution to investigate the ET from weak to strong electronic couplings [31,40–44], and then incorporate the solvent effect [32,45–49].

When solvent dynamics is very fast, ET can be thought as a conventional nonadiabatic chemical reaction. Thus, TST idea can be employed to treat ET. Therefore, we utilize the generalized NA-TST [40] and the Zhu-Nakamura (ZN) nonadiabatic transition probability to formulate a formula similar to the Marcus formula. The purpose is to propose a simple yet accurate formula directly applicable to explaining experimental data, e.g., the localized to delocalized ET in mixed-valence molecules [50,51]. The formulation can also be used to define the ‘sink’ function [22,52,53] along the fast vibrational modes in the treatment of solvent controlled ET reactions. The difference of the present formula from the Marcus one is the pre-factor which is defined by the thermally averaged ZN formula. The generalized NA-TST is formulated based on the Miller’s reactive flux-flux correlation function approach. The ZN theory [10,11,54–56], on the other hand, is practically free from the drawbacks of the LZ theory mentioned above. It covers all the energy range in the two-state curve crossing problem. Numerical tests [40] have also shown that it is very essential for accurate evaluation of the thermal rate constant to take into account the multi-dimensional topography of the seam surface and treat the nonadiabatic electronic transition and nuclear tunneling effects properly. The presently proposed formula of ET rate is applicable to multidimensional systems and the Monte-Carlo approach can be utilized.

As solvent dynamics begins to control ET, the thermal equilibrium distribution of the donor population becomes doubtful. The case is similar to adiabatic chemical reactions in solvent. The Kramers theory has been developed to study the thermal rate of escape from a metastable state coupled to a dissipative thermal environment. Pollak, Grabert and Häggi [57] formulate a unified theory–Kramers turnover theory, which covers the whole range of friction strength and is applicable to an arbitrary memory friction. Later, quantum tunneling was incorporated by Rips and Pollak [58] with the use of the parabolic approximation. Our strategy for ET is to extend the quantum Kramers theory [58–60] because the Kramers theory can be considered as a special case of studying nonadiabatic chemical reactions where the nonadiabatic coupling is so strong that the reaction is controlled by the lower adiabatic potential energy surface (PES). However, for ET the contribution of the upper PES cannot be neglected in generality. Therefore, one has to consider both the electronic and nuclear tunneling effects. In the present work, those effects are incorporated using the formulas proposed by Zhu and Nakamura [56].

If the solvent dynamic motion is much slower than intramolecular vibrational motions, the better approach is so-called the Sumi-Marcus theory. As already mentioned, the original Sumi-Marcus theory is limited to the weak electronic coupling. We extend the Sumi-Marcus theory by introducing two approaches to investigate electron transfer reactions from weak to strong electronic coupling regime. One of these approaches is the quantum R-matrix theory which is useful in dealing with the intramolecular vibrational motions in the whole electronic coupling domain. The other is the split operator approach that is employed to solve the reaction-diffusion equation.

In the following, we present the details of different approaches. Section 2 describes the semiclassical approximation of thermal rate constants for electron transfer in fast solvent relaxation limit. Section 3 discusses the solvent controlled ET theories. Section 4 gives several applications. The concluding remarks are given in Section 5.

Semiclassical approximation of thermal rate constants for electron transfer in fast solvent relaxation limit

Before showing our approaches, we discuss the fundamental idea for ET process, as well as simple approaches to estimate the ET rate. ET reactions can be considered as a special case of the nonadiabatic chemical reactions. Therefore, we consider general two-state reactions. The potential for these reactions is written in the diabatic representation

(1)V=|\mathrm{1}\rangle V_{\mathrm{11}}\langle \mathrm{1}|+|\mathrm{2}\rangle V_{\mathrm{22}}\langle \mathrm{2}|+|\mathrm{1}\rangle V_{\mathrm{12}}\langle \mathrm{2}|+|\mathrm{2}\rangle V_{\mathrm{21}}\langle \mathrm{1}|,

where V_ii represent the diabatic potentials of the i-th state, and V_ij are electronic couplings between the i-th and j-th states. One may alternatively express Eq. (1) in the adiabatic representation by

(2)V=|{\mathrm{1}}_a\rangle E_{\mathrm{11}}\langle {\mathrm{1}}_a|+|{\mathrm{2}}_a\rangle E_{\mathrm{22}}\langle {\mathrm{2}}_a|.

Here, adiabatic potentials E_ii are related with the diabatic potentials by

(3)E_{\mathrm{1},\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[(V_{\mathrm{11}}+V_{\mathrm{22}})\pm \sqrt{{(V_{\mathrm{11}}-V_{\mathrm{22}})}^{\mathrm{2}}+\mathrm{4}{V}_{\mathrm{12}}^{\mathrm{2}}}\right].

In this case, the electronic coupling is transformed to the kinetic operator. Figure 1 displays the schematic potentials for one-dimensional nonadiabatic chemical reactions.

The easy implementation for the thermal reaction rate is starting from micro-canonical reaction rate. Assuming that the nonadiabatic transition probability P(E) from the left to the right in Fig. 1 is known at a given energy E, the thermal rate k can be obtained by an average with the Boltzmann thermal distribution as a weight function

(4)k={(\mathrm{2}\mathrm{\pi }\mathit{Ћ}Z)}^{-\mathrm{1}}{\displaystyle {\int }_{\mathrm{0}}^{\infty }{\mathrm{e}}^{-E/(k_{B}T)}P(E)},

where Z is the partition function of the reactant, T is the temperature. Therefore, the key point for the calculation of the rate becomes to evaluate P(E). The investigation of the analytical expression for P(E) has a long history. Landau and Zener in 1932 [7–9] have deduced it which now becomes the well-known Landau-Zener (LZ) formula. The defect of LZ formula was overcome by Zhu and Nakamura [10,11,54–56] only recently. We have also proposed a rigorous R-matrix approach numerically based on quantum scattering theory [44].

Equation (4) can be straightforwardly used to ET, where the diabatic potentials are harmonic oscillators. When electronic coupling is weak, the rate is alternatively calculated by Fermi Golden rule. For example, the rate is given by [61,62]

(5)k={\displaystyle \frac{\mathrm{2}\mathrm{\pi }{\left|H_{\mathrm{12}}\right|}^{\mathrm{2}}}{\mathit{Ћ}\mathrm{2}\omega }}{\left({\displaystyle \frac{\overline{v}+\mathrm{1}}{\overline{v}}}\right)}^{p/\mathrm{2}}\exp \left[-S\left(\mathrm{2}\overline{v}+\mathrm{1}\right)\right]I_p\left[\mathrm{2}S{\left(\overline{v}\left(\overline{v}+\mathrm{1}\right)\right)}^{\mathrm{1}/\mathrm{2}}\right],

where ω is the vibrational frequency, \overline{v}={[\exp (\mathit{Ћ}\omega /T)-\mathrm{1}]}^{-\mathrm{1}} is the thermal population, S is the Huang-Rays factor, and I_p is the Bessel function with p=\mathrm{\Delta }G/(\mathit{Ћ}\omega ). At high temperature approximation, Eq. (5) is just well-known Marcus formula

(6)k_{\mathrm{M}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{u}\mathrm{s}}={\displaystyle \frac{H_{AB}^{\mathrm{2}}}{\mathit{Ћ}}}\sqrt{{\displaystyle \frac{\mathrm{\pi }\beta }{\lambda }}}{\mathrm{e}}^{-\frac{\beta {(\lambda +\mathrm{\Delta }G)}^{\mathrm{2}}}{\mathrm{4}\lambda }},

where \lambda \left(={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}S\omega \right) is reorganization energy, and β = 1/(k_BT).

With the above knowledge, let us consider ET in complex systems. In these systems, Eq. (4) cannot be straightforwardly used because one has to consider multi-dimensional problems and does not know the reaction direction. Despite that the perturbation theory is applicable for multi-dimensional systems, it is essentially limited to the weak electronic coupling. Therefore, our strategy is to start from a generous quantum rate expression. It is known that the rigorous chemical reaction rate constants can be evaluated by the flux-side correlation function formula [63]

(7)k=\mathrm{1}/Z_r\underset{t\to \infty }{\lim }\mathrm{t}\mathrm{r}\left[{\mathrm{e}}^{-\beta \hat{H}}\hat{F}{\mathrm{e}}^{i\hat{H}}{}^{t/\mathit{Ћ}}\hat{h}{\mathrm{e}}^{-i\hat{H}}{}^{t/\mathit{Ћ}}\right],

where Z_r is the partition function of the reactant, \hat{H} is the Hamiltonian of the system, \hat{h} is Heaviside function, \hat{F} is the flux operator defined as \hat{F}={\displaystyle \frac{i}{\mathit{Ћ}}}\left[\hat{H},\hat{h}\right]. For two-state ET, the reactions involve two PESs. To enable the use of Eq. (7), one may transform the diabatic PESs of ET to adiabatic ones. Consider that the upper adiabatic PES has no open reactive pathway (a typical case for ET). The reaction occurs only on the lower adiabatic PES, but the probability is very much affected by the nonadiabatic coupling. Using this simple judgment, one may obtain a classical expression of Eq. (7) [40],

(8)k=\mathrm{1}/Z_r\underset{t\to \infty }{\lim }\mathrm{1}/h^N\int \mathrm{d}\mathbf{P}\mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta H(\mathbf{P},\mathbf{Q})}{\mathbf{P}}^{\mathbf{T}}\nabla (S(\mathbf{Q})-{\xi }_{\mathrm{0}})\delta (S(\mathbf{Q})-{\xi }_{\mathrm{0}})h[S(\mathbf{Q}(t))-{\xi }_{\mathrm{0}}]

with replacements of the time dependent Heaviside function by its classical analog and the quantum trace by the phase space integration. Here, S(\mathbf{Q})-{\xi }_{\mathrm{0}}=\mathrm{0} determines the generalized seam surface, and N is the degrees of freedom.

Using standard classical molecular dynamics approaches for infrequent events (see, e.g., Ref. [64].), Eq. (8) can be cast into

(9)k(t)=k_{\mathrm{T}\mathrm{S}\mathrm{T}}\kappa (t),

where k_TST represents the rate of the nonadiabatic transition state theory and it has the following form

(10)k_{\mathrm{T}\mathrm{S}\mathrm{T}}\equiv k(\mathrm{0})=\frac{\mathrm{1}}{Z_r}\frac{\mathrm{1}}{h^N}\int \mathrm{d}\mathbf{P}\mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta H(\mathbf{P},\mathbf{Q})}{\mathbf{P}}^{\mathbf{T}}\nabla S(\mathbf{Q})\delta (S(\mathbf{Q})-{\xi }_{\mathrm{0}})h[S(\mathbf{Q}(\mathrm{0}))-{\xi }_{\mathrm{0}}],

and κ(t) is the dynamic transmission coefficient

(11)\kappa (t)={\displaystyle \frac{{\displaystyle \int \mathrm{d}\mathbf{P}\mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta H(\mathbf{P},\mathbf{Q})}{\mathbf{P}}^{\mathbf{T}}\nabla (S(\mathbf{Q})-{\xi }_{\mathrm{0}})\delta (S(\mathbf{Q})-{\xi }_{\mathrm{0}})h[S(\mathbf{Q}(t))-{\xi }_{\mathrm{0}}]}}{{\displaystyle \int \mathrm{d}\mathbf{P}\mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta H(\mathbf{P},\mathbf{Q})}{\mathbf{P}}^{\mathbf{T}}\nabla (S(\mathbf{Q})-{\xi }_{\mathrm{0}})\delta (S(\mathbf{Q})-{\xi }_{\mathrm{0}})h[S(\mathbf{Q}(\mathrm{0}))-{\xi }_{\mathrm{0}}]}}}.

For most realistic systems, the re-crossing effect of ET reaction can be neglected, i.e., κ(t) ≡ 1. In such a case, k_TST itself should predict reasonable ET rate. However, when the transferred electron in the acceptor state has not enough time to relax its energy to environment, it can be back to the donor state. In this energy diffusion limit, the dynamic transmission coefficient κ(t) may significantly reduce the ET rate.

Nonadiabatic transition state theory and Marcus formula

Let us consider k_TST first. It is known that h[S(Q(0)) - ξ₀] in Eq. (10) represents the reactive component starting from the transition state to product region without considering the recrossing dynamic effect. When the nuclear tunneling effect is incorporated, Eq. (10) becomes [40]

(12)k=Z^{-\mathrm{1}}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\beta }}}{\displaystyle \int \mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta V_{\mathrm{1}}(\mathbf{Q})}P(\beta ,\mathbf{Q})\left|\nabla S(\mathbf{Q})\right|\delta (S(\mathbf{Q})-{\xi }_{\mathrm{0}})},

where Q represents the collective mass-scaled nuclear coordinates of N degrees of freedom, Z={\displaystyle \int \mathrm{d}\mathbf{Q}\exp (-\beta V_{\mathrm{1}}(\mathbf{Q}))(\beta =\mathrm{1}/(k_{B}T))} is the classical partition function of the diabatic potential energy surface V₁(Q) for the reactant, and the effective transition probability P(β, Q) is evaluated with the use of the ZN formula P_ZN(E) [10,56] by

(13)P(\beta ,\mathbf{Q})=\beta {\displaystyle {\int }_{\mathrm{0}}^{\infty }\mathrm{d}E{\mathrm{e}}^{-\beta (E-V_{\mathrm{1}}(\mathbf{Q}))}{P}_{\mathrm{Z}\mathrm{N}}(E,\mathbf{Q})}.

By introducing the free energy profile F(ξ)

(14){\mathrm{e}}^{-\beta F(\xi )}={\displaystyle \int \mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta V_{\mathrm{1}}(\mathbf{Q})}\left|\nabla (S(\mathbf{Q}))\right|\delta (S(\mathbf{Q})-\xi )},

Equation (12) can be cast into

(15)k=Z^{-\mathrm{1}}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\beta }}}{\displaystyle \int \mathrm{d}\xi \delta (\xi -{\xi }_{\mathrm{0}}){\mathrm{e}}^{-\beta F(\xi )}\overline{P}(\beta ,\xi )},

where the average transition probability \overline{P}(\beta ,\xi ) is defined by

(16)\overline{P}(\beta ,\xi )={\displaystyle \frac{\beta {\displaystyle \int \mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta V_{\mathrm{1}}(\mathbf{Q})}\left|\nabla S(\mathbf{Q})\right|\delta (\xi -S(\mathbf{Q}))P(\beta ,\mathbf{Q})}}{{\displaystyle \int \mathrm{d}\mathbf{Q}{\mathrm{e}}^{-\beta V_{\mathrm{1}}(\mathbf{Q})}\left|\nabla (S(\mathbf{Q}))\right|\delta (\xi -S(\mathbf{Q}))}}}.

It is worth mentioning that the Hamiltonian in Eq. (15), in principle, can include intramolecular high frequency modes and solvent low frequency modes together. The crossing seam surface ξ₀ = S(Q) of the potentials of donor and acceptor is taken as the nonadiabatic transition state. It can be shown that Eqs. (15) and (16) give essentially the same expression obtained in Refs. [65–67], (For instance, Eq. (3.18) of Ref. [66]), when the one dimensional reaction coordinate is assumed and the “nuclear tunneling” effect is neglected. In our treatment the multi-dimensionality is taken into account and the instantaneous normal mode analysis is made at the transition point to determine the normal direction to the seam, i.e., the reaction coordinate. In this way, the reaction coordinate has the maximum mean free path and thus the theory is supposed to be applicable in a wide range of friction.

To find the relationship between Eq. (15) and the Marcus theory, we consider the linear response limit. In this case, the free energy F(ξ) becomes a parabolic function of ξ and Eq. (15) can be recast into

(17)k=\kappa k_{\mathrm{M}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{u}\mathrm{s}},

where the prefactor κ is given by

(18)\kappa ={\displaystyle \frac{\mathit{Ћ}\omega }{\mathrm{2}\mathrm{\pi }{H}_{AB}^{\mathrm{2}}}}\sqrt{{\displaystyle \frac{\lambda }{\mathrm{\pi }\beta }}}\overline{P}(\beta ,{\xi }_{\mathrm{0}}).

In Eq. (17), the effects of nonadiabatic transition including the “nuclear tunneling” are properly taken into account by κ and naturally the main task is to evaluate the thermally averaged transition probability \overline{P}(\beta ,\xi ), which has to be evaluated with the Monte-Carlo technique for multi-dimensional systems.

It is easily shown that the Marcus-Hush formula in adiabatic and nonadiabatic limits can be recovered from Eq. (17) - Eq. (18) within the high temperature approximation.

To demonstrate the accuracy of Eq. (15), Figure 2 shows the ET rate vs. electronic coupling for a simple model system, and it confirms that the present approach works quite well from weak to strong electronic coupling regime.

Figure 3 shows the ET rate vs. electronic coupling for a multi-dimensional system in Marcus inverted regime [43]. It is shown that the rate initially increases quadratically in the weak coupling regime, as predicted by the perturbation theory. As the coupling increases, the perturbation theory always overestimates the rate. This is different from the Marcus normal case where the perturbation theory does not necessarily provide a larger rate than the real one [41]. The rate predicted by the present approach reaches maximum at a certain value of coupling strength and decreases with further increase of the coupling. This behavior of the adiabatic suppression is the property peculiar to the inverted case, since the large electronic coupling makes adiabatic potentials separate far from each other and the nonadiabatic transition probability decreases as a result. Compared with the prediction from the LZ formula, the ZN formulas give much larger rate especially at low temperatures. With the temperature decreasing, the coupling strength at the maximal rate predicted by the LZ formula shifts toward smaller value, while it keeps nearly unchanged in the case of the ZN formulas. This interesting behavior can be understood from the energy dependence of the nonadiabatic transition probability. In the case of the ZN formulas, we find that the probability reaches maximum at a certain intermediate coupling strength at energies in the region around the crossing point.

Dynamic transmission coefficient

In this section we shall present the detailed description for the evaluation of the transmission coefficient. In the expression of k_TST, we have incorporated the contribution from the nuclear tunneling effect where kinetic energies are lower than the barrier of the lower adiabatic PES. One may safely neglect this contribution to the transmission coefficient because of the small tunneling recrossing probability. In this case, Eq. (11) is straightforwardly used for the numerical simulation by using the surface hopping approach. In the Tully’s surface hopping technique [68] it is assumed that the hop from one diabatic PES to another occurs when trajectories pass through the seam surface and the hopping probability is given by the single crossing probability. For infrequent events, the transition probability may be given by the total crossing probability from reactant to product regime. Thus, one may implement transmission coefficient with standard reactive flux methods [69–74] for adiabatic chemical reactions. But, when the trajectories meet the seam surface, one has to judge whether they cross the seam surface or are reflected in terms of the total transition probability. The method [31] is summarized as follows: After initializing the coordinates and momentum from the Boltzmann distribution on the seam surface, the trajectories are propagated without judgment of the hop, (since the hopping probability has been involved in the denominate in Eq. (11) we only consider the trajectories which have already crossed the seam surface.) When the trajectories meet the seam surface during the propagation, the hopping probability is evaluated by the ZN formula [10,11,54–56] with the momentum component P_n along the normal direction of the seam surface. When we use the ant-eater method [68], the hop occurs as long as the probability is greater than a random number generated by the computer. When the hop happens, we propagate the trajectories continuously without making any change for the momenta. Otherwise, we change the direction of the momentum P_n and propagate them as though they still stay on the same PES and are reflected. To set new momenta after changing the sign of momentum P_n, we use a relation

(19){\mathbf{P}}_{\mathrm{o}\mathrm{l}\mathrm{d}}={\mathbf{P}}_{||}+P_n\mathbf{n}

with {\mathbf{P}}_{||}={\mathbf{P}}_{\mathrm{o}\mathrm{l}\mathrm{d}}-P_n\mathbf{n}. Change the sign of the second term in Eq. (19), one gets the new momenta

(20){\mathbf{P}}_{\mathrm{n}\mathrm{e}\mathrm{w}}={\mathbf{P}}_{\mathrm{o}\mathrm{l}\mathrm{d}}-\mathrm{2}{P}_n\mathbf{n}.

During the propagation, once the trajectories meet the seam surface we use the same hopping procedure as in the starting step. Finally, the transmission coefficient is obtained from the plateau value of Eq. (11).

To numerically check the present hopping approach, we take a spin-boson model as an example. Figure 4 displays the dynamic transmission coefficient vs. friction coefficient of environment. As the coherent hopping is incorporated [31], the present results are in agreement with those from the rigorous path integral calculation [25] in the whole friction regime.

Electron transfer rate theory incorporated solvent dynamic effect

Like adiabatic chemical reactions in solvent where the TST idea becomes invalid, the approaches based on the fast relaxation solvent approximation are also invalid for ET because solvent dynamics plays a significant role and it even controls ET. In such a case, we have proposed two kinds of approaches, quantum Kramer-like theory and extended Sumi-Marcus approach, to investigate ET. The former is suitable for low to intermediate solvent friction limits while the latter is useful for the case of very strong friction.

Quantum Kramer-like theory

We start from the reaction coordinate model of ET. Its Hamiltonian in the diabatic representation has a conventional form. (The mass scaled coordinate and \mathit{Ћ}=\mathrm{1} are used.)

(21)H=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }E{\hat{\sigma }}_z+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{H}_{AB}{\sigma }_x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{p}_q^{\mathrm{2}}\mathbf{I}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\omega }_{\mathrm{0}}^{\mathrm{2}}{(q\mathbf{I}+\kappa q_{\mathrm{0}}{\sigma }_z)}^{\mathrm{2}}+H_{\mathrm{0}}\mathbf{I}.

Here σ_x and σ_z are the Pauli matrices. I is a unit matrix. The potentials along the reaction coordinate q are characterized by two shifted harmonic oscillators with frequency ω₀, shift 2κq₀, energy bias \mathrm{\Delta }E and electronic coupling element H_AB/2. The environment Hamiltonian H₀ is modeled with a collection of harmonic oscillators

(22)H_{\mathrm{0}}={\displaystyle \sum _j{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}}\left[p_j^{\mathrm{2}}+{\omega }_j^{\mathrm{2}}{\left(x_j+{\displaystyle \frac{c_{j}q}{{\omega }_j^{\mathrm{2}}}}\right)}^{\mathrm{2}}\right],

where x_j, p_j, and ω_j are, respectively, the coordinate, momentum, and frequency of the i-th harmonic oscillator representing the environment, and c_j is the coupling strength of the reaction coordinate to the i-th oscillator. The complete information needed in Eq. (22) is encapsulated in the single spectral density J(ω), defined by the expression

(23)J(\omega )={\displaystyle \frac{\mathrm{\pi }}{\mathrm{2}}}{\displaystyle \sum _j{\displaystyle \frac{c_j^{\mathrm{2}}}{{\omega }_j}}[\delta (\omega -{\omega }_j)-\delta (\omega +{\omega }_j)]}.

Thus, the model is completely defined by the parameters ω₀, κq₀, \mathrm{\Delta }E, H_AB and function J(ω).

The Hamiltonian, Eq. (21), may be also written in the adiabatic representation. One first diagonalizes the potential terms and this leads to two adiabatic PESs. In the adiabatic ET limit, the ET process is solely controlled by the lower PES. The corresponding adiabatic Hamiltonian is given by

(24)H_{\mathrm{a}\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{p}_q^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\omega }_{\mathrm{0}}^{\mathrm{2}}[q^{\mathrm{2}}+{(\kappa q_{\mathrm{0}})}^{\mathrm{2}}]-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{[{(E_{r}q/\kappa q_{\mathrm{0}}-\mathrm{\Delta }E)}^{\mathrm{2}}+H_{AB}^{\mathrm{2}}]}^{\mathrm{1}/\mathrm{2}}+{\displaystyle \sum _j{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[p_j^{\mathrm{2}}+w_j^{\mathrm{2}}{\left(x_j+{\displaystyle \frac{c_{j}q}{w_j^{\mathrm{2}}}}\right)}^{\mathrm{2}}\right]}

where E_r is the reorganization energy. In this case, the ET can be viewed as the overbarrier motion of a particle exerted by the friction and random forces, and the Kramers theory is straightforwardly used to evaluate the ET rate [17].

For the nonadiabatic ET reactions the upper adiabatic PES cannot be neglected. However, we have found in recent work [31,40,41] that such reactions themselves may be described still by the adiabatic Hamiltonian Eq. (24), but the nonadiabatic tunneling effects should be properly incorporated. The physical basis behind is the surface hopping proposed by Tully [68]. One, therefore, may extend the quantum Kramers theory to the ET process. To do so, we make the normal mode analysis in the vicinity of the barrier of the potential in Eq. (24). The standard procedure reads [75]

(25)H_{\mathrm{a}\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\dot{\rho }}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\lambda }_{\mathrm{0}}^{\mathrm{2}}{\rho }^{\mathrm{2}}+{\displaystyle \sum _j{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\dot{y}}_j^{\mathrm{2}}+{\lambda }_j^{\mathrm{2}}{y}_j^{\mathrm{2}}\right)+V(\rho ,y_i)},

where ρ and y_i are mass weighted coordinates of the unstable normal mode and the stable modes, respectively, and V (ρ, y_i) represents the nonlinear interaction between the modes outside the barrier region. As pointed out by Grabert [76], the crucial quantity needed for an estimate of the escape rate from the donor well is the steady energy distribution n(є) in the unstable normal mode ρ. In terms of n(є), the forward ET rate can be written as

(26)k_f={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }\beta }}{\displaystyle {\int }_{-\infty }^{\infty }\mathrm{d}\mathit{є}{P}_{\mathrm{Z}\mathrm{N}}}(\mathit{є})n(\mathit{є}),

since all particles reaching the crossing barrier hop with nonadiabatic transition probability P_ZN(є). In Eq. (26), β is the inverse of temperature. To distinguish it from the Kramers theory, we name Eq. (26) as the quantum Kramers-like formula. In the present work, we use the ZN formulas [11,56] to evaluate P_ZN(є).

What remains to be done is to determine the steady energy distribution function n(є). Following Mel’nilov and Meshkov [77], we introduce two distributions f₁(є) and f₂(є) of the particles, moving toward the crossing point from donor and acceptor wells, respectively. They satisfy the set of Mel’nilov’s equations

(27)f_1\left(\mathit{є}\right)={\displaystyle {\int }_{-\infty }^{+\infty }d{\mathit{є}}^{\prime }{P}_1\left(\mathit{є}|{\mathit{є}}^{\prime }\right)\left[\begin{array}{c}\hfill \left(1-P_{\mathrm{Z}\mathrm{N}}\left({\mathit{є}}^{\prime }\right)\right)f_1\left({\mathit{є}}^{\prime }\right)\hfill \\ \hfill +P_{\mathrm{Z}\mathrm{N}}\left({\mathit{є}}^{\prime }\right)f_2\left({\mathit{є}}^{\prime }\right)\hfill \end{array}\right]},

(28)f_2\left(\mathit{є}\right)={\displaystyle {\int }_{-\infty }^{+\infty }\mathrm{d}{\mathit{є}}^{\prime }{P}_2\left(\mathit{є}|{\mathit{є}}^{\prime }\right)\left[\begin{array}{c}\hfill \left(1-P_{\mathrm{Z}\mathrm{N}}\left({\mathit{є}}^{\prime }\right)\right)f_2\left({\mathit{є}}^{\prime }\right)\hfill \\ \hfill +P_{\mathrm{Z}\mathrm{N}}\left({\mathit{є}}^{\prime }\right)f_1\left({\mathit{є}}^{\prime }\right)\hfill \end{array}\right]},

where P_i(\mathit{є}|{\mathit{є}}^{\prime })(i=\mathrm{1},\mathrm{2}) is the condition probability that the particle leaves the crossing barrier region with energy є' and returns after a round trip with energy є. n(є) is related to f₁(є) and f₂(є) through the relation

(29)n(\mathit{є})=f_{\mathrm{1}}(\mathit{є})-f_{\mathrm{2}}(\mathit{є}).

For a symmetric double well, inserting Eq. (27), Eq. (28) into Eq. (29) leads to the fundamental integral equation

(30)n(\mathit{є})={\displaystyle {\int }_{-\infty }^{+\infty }\mathrm{d}{\mathit{є}}^{\prime }}P(\mathit{є}|{\mathit{є}}^{\prime })(\mathrm{1}-\mathrm{2}{P}_{\mathrm{Z}\mathrm{N}}({\mathit{є}}^{\prime }))n({\mathit{є}}^{\prime }).

Equation (30) has a boundary condition

(31)n(\mathit{є})={\displaystyle \frac{C}{\mathrm{2}\mathrm{\pi }\beta }}{\mathrm{e}}^{-\mathit{є}},\mathit{є}\to -\infty

where C is the normalized constant and it has a form

(32)C={\displaystyle \frac{\mathrm{2}{\omega }_{\mathrm{0}}}{{\omega }_b}}\sin \left({\displaystyle \frac{\beta {\lambda }_{\mathrm{0}}}{\mathrm{2}}}\right)\mathrm{\Xi }{\mathrm{e}}^{-\beta V_b}.

Here, ω_b is the barrier frequency; V_b is the barrier height. The ratio of quantum partition functions \mathrm{\Xi } is given by

(33)\mathrm{\Xi }\equiv {\displaystyle \prod _{n=\mathrm{1}}^{\infty }{\displaystyle \frac{{\omega }_{\mathrm{0}}^{\mathrm{2}}+\hat{\omega }{}_n^{\mathrm{2}}+{\hat{\omega }}_n\hat{\gamma }\left(\hat{\omega }n\right)}{-{\omega }_b^{\mathrm{2}}+\hat{\omega }{}_n^{\mathrm{2}}+{\hat{\omega }}_n\hat{\gamma }\left(\hat{\omega }n\right)}}},

where {\hat{\omega }}_n=\mathrm{2}\mathrm{\pi }n/\beta is the Matsubara frequencies [78], and \hat{\gamma }(s) is the Laplace transform of the time-dependent friction

(34)\hat{\gamma }(s)={\displaystyle {\int }_{\mathrm{0}}^{\infty }\mathrm{d}t{\mathrm{e}}^{-st}}\gamma (t).

In principle, combining Eq. (26) and Eq. (30) with its boundary condition Eq. (31) can predict the ET rate from the weak to strong friction regime as well as from the nonadiabatic to adiabatic limit once the P(\mathit{є}|{\mathit{є}}^{\prime }) is known.

The value of P(\mathit{є}|{\mathit{є}}^{\prime }) dominates the biased distribution of n(є) from the thermal one. It is explicitly needed only when the unstable mode and stable modes are weakly coupled in the well region, because otherwise the mode energy is thermalized and n(є) satisfies the Boltzmann distribution. When the stable modes are thermally distributed the condition probability in the classical limit satisfies a Gaussian distribution whose width is decided by the average energy loss of the unstable-mode motion [59]. In the present work, we focus on the strong friction region, which implies that the mode energy is thermalized. (The energy loss parameter is much larger than 1.)

To demonstrate the possible applications of the current approach for the whole electronic coupling regimes, we plot the ET rate as a function of the coupling strength in Fig. 5 [32]. We set the parameters \mathrm{\Delta }=\mathrm{1}, ω_c = 0.1 and T = 1 to keep the barrier relatively high in the lower adiabatic PES. The current approach predicts correctly the rate in the nonadiabatic limit (H_AB<0.006) and the rate does change quadratically with H_AB. It recovers the Kramers theory in the adiabatic limit (H_AB>0.03) and becomes consistent with the Zusman theory in between. As expected, Zusman theory cannot describe the ET rate in the adiabatic limit. It predicts a much larger rate than the Kramers theory and tends to the TST.

Extended Sumi-Marcus theory [46,49]

Since the fluctuations of the high frequency modes are assumed to be faster than the solvent one the thermal equilibrium is maintained for those fast modes during the course of reactions. The nonadiabatic transition state theories, therefore, are applicable to obtaining the sink function. Sumi and Marcus [20] have treated the sink function with the semiclassical local Golden rule for the nonadiabatic ET, and the conventional transition state theory for the adiabatic reactions. For many systems, however, the electronic coupling strength may locate the regimes between the nonadiabatic and adiabatic transitions, making the perturbation theory and adiabatic transition state theory invalid. It is still a challenge to get the sink function in the crossover regime from the nonadiabatic to the adiabatic limit.

For the ET affected by the motions of the high-frequency vibrational modes and the solvent polarization fluctuation, the potential energy surfaces of the donor (i = 1) and the acceptor (i = 2) can be modeled as

(35)V_{\mathrm{1}}(q,x)=V_{\mathrm{1}}(q)+V_{\mathrm{1}}(x),

(36)V_{\mathrm{2}}(q,x)=V_{\mathrm{2}}(q)+V_{\mathrm{2}}(x)+\mathrm{\Delta }{G}_{\mathrm{0}}.

Here q represents the collective coordinate of the intramolecular vibrational modes; x is the solvent coordinate that is defined either by the mean square fluctuation of the excess solvent polarization or by the energy gap of the solvation energy upon moving an electron from the donor to the acceptor; \mathrm{\Delta }{G}_{\mathrm{0}} is the exothermicity of the reaction. Generally, one has to consider two-dimensional diffusion equations to incorporate the relaxation dynamics both from slow x and fast q coordinates. In the case where the fluctuations of the intramolecular vibrational modes are quite fast compared with the solvent relaxation, the adiabatic elimination procedure [33,35,37] can lead the two-dimensional diffusive equations to one-dimensional ones. Denoting P₁(x, t) and P₂(x, t) as the population distributions of the donor and the acceptor states at a given x and time t, respectively, one obtains the one-dimensional coupled diffusion-reaction Eqs. [20].

(37)\partial P_{\mathrm{1}}(x,t)/\partial t=[L_{\mathrm{1}}-k_{\mathrm{1}}(x)]P_{\mathrm{1}}(x,t)+k_{\mathrm{2}}(x)P_{\mathrm{2}}(x,t)

(38)\partial P_{\mathrm{2}}(x,t)/\partial t=[L_{\mathrm{2}}-k_{\mathrm{2}}(x)]P_{\mathrm{2}}(x,t)+k_{\mathrm{1}}(x)P_{\mathrm{1}}(x,t)\mathrm{.}

Here k_i(x) is the sink function that invokes the contribution from the fast modes. The generalized Smoluchowski operator L_i is given by

(39)L_i=D(t)\{{\partial }^{\mathrm{2}}/\partial x^{\mathrm{2}}+\beta \partial /\partial x[d{V}_i(x)/dx]\}

with a time-dependent diffusion coefficient D(t). D(t) is related to the solvent coordinate correlation function \mathrm{\Delta }(t) by

(40)D(t)=-{\displaystyle \frac{\mathrm{1}}{\beta }}{\displaystyle \frac{\mathrm{d}\ln \mathrm{\Delta }(t)}{\mathrm{d}t}},

where β = 1/(k_bT), and T is the temperature. When the backward ET from the acceptor to donor is ignored, i.e., for irreversible reactions, Eq. (37) and Eq. (38) become a single reaction-diffusion equation

(41){\displaystyle \frac{\partial P\left(x,t\right)}{\partial t}}=D\frac{{\partial }^{2}P\left(x,t\right)}{\partial x^2}+\frac{D}{k_{B}T}\frac{\partial }{\partial x}\left[P\left(x,t\right){\displaystyle \frac{d{V}_1}{dx}}\right]-k_1\left(x\right)P\left(x,t\right)\mathrm{.}

Note that the sink function represents the local reaction rate along the intramolecular vibrational coordinate q at a given solvent coordinate x. It can be obtained from Eq. (4) together with ZN formula or R-matrix approach for wide electronic coupling strengths.

The remaining problem is how to effectively solve the reaction-diffusion equations. It is noted that the coupled equations Eq. (37) and Eq. (38) can be transformed into Hermitian forms. To do so, one may use the substitution

(42)p_i(x,t)=P_i(x,t)/g_i(x),

where g_i(x) is the square root of equilibrium solution of the coupled equations in the absence of the sink functions

(43)g_i(x,t)={[P_i^{\mathrm{e}\mathrm{q}}(x)]}^{\mathrm{1}/\mathrm{2}}=\exp [-\beta V_i(x)/\mathrm{2}]/{\displaystyle \int dx\exp [-\beta V_i(x)/\mathrm{2}]}\mathrm{.}

Insert Eq. (42) into Eq. (37) and Eq. (38), one gets

(44){\displaystyle \frac{\partial }{\partial t}}\left(\begin{array}{c}\hfill p_{\mathrm{1}}(x,t)\hfill \\ \hfill p_{\mathrm{2}}(x,t)\hfill \end{array}\right)=\left[\begin{array}{cc}\hfill H_{\mathrm{1}}(t)\hfill & \hfill k_{\mathrm{2}}^{\prime }(x)\hfill \\ \hfill k_{\mathrm{1}}^{\prime }(x)\hfill & \hfill H_{\mathrm{2}}(t)\hfill \end{array}\right]\left(\begin{array}{c}\hfill p_{\mathrm{1}}(x,t)\hfill \\ \hfill p_{\mathrm{2}}(x,t)\hfill \end{array}\right).

Here

(45)k_{\mathrm{1}}^{\prime }(x)=k_{\mathrm{1}}(x)g_{\mathrm{1}}(x)/g_{\mathrm{2}}(x),

(46)k_{\mathrm{2}}^{\prime }(x)=k_{\mathrm{2}}(x)g_{\mathrm{2}}(x)/g_{\mathrm{1}}(x)\mathrm{.}

H_i is the Hermitian form of the Smoluchowski diffusion operator

(47)H_i(t)=T_i(t)+V_i(t),

where T_i is the time-dependent kinetic operator

(48)T_i(t)=D(t)\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}},

and V_i(t) is time-dependent potential

(49)V_i(t)={\displaystyle \frac{\beta D(t)}{\mathrm{2}}}\left[{\displaystyle \frac{\beta }{\mathrm{2}}}{\left({\displaystyle \frac{\mathrm{d}{V}_i(x)}{\mathrm{d}x}}\right)}^{\mathrm{2}}-{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{V}_i(x)}{\mathrm{d}{x}^{\mathrm{2}}}}\right]-k_i(x)\mathrm{.}

It is easy to see that Eq. (44) becomes standard coupled Schrödinger equations by replacing t with it in the derivative, i.e., replacing the real time by an imaginary one. Therefore, the numerical approaches solving the Schrödinger equation can be straightforwardly used. In this work, we choose the imaginary-time split operator approach (ITSOA) which is so named to get distinguished from the real time one in solving Schrödinger equations.

To apply the ITSOA, we write the formal solution of Eq. (37) as

(50)\mathbf{p}(x,t)={\exp }_{+}\left({\displaystyle {\int }_{\mathrm{0}}^t\mathrm{d}\tau \mathbf{H}}(\tau )\right)\mathbf{p}(x,\mathrm{0}),

where p(x, t) represents the array of p₁ and p₂. Under the Trotter approximation for a short time dt propagation, the population at t + dt can be evaluated by

(51)\mathbf{p}\left(x,t+\mathrm{d}t\right)=\exp (\mathbf{T}(t)\mathrm{d}t/\mathrm{2})\exp (\mathbf{V}(t)\mathrm{d}t)\exp (\mathbf{T}(t)\mathrm{d}t/\mathrm{2})\mathbf{p}(x,t),

where T(t) is the kinetic operator

(52)\mathbf{T}(t)=\left[\begin{array}{cc}\hfill T_{\mathrm{1}}(t)\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill T_{\mathrm{2}}(t)\hfill \end{array}\right],

and V represents the potentials

(53)\mathbf{V}(t)=\left[\begin{array}{cc}\hfill V_{\mathrm{1}}(t)\hfill & \hfill k_{\mathrm{2}}^{\prime }\hfill \\ \hfill k_{\mathrm{1}}^{\prime }\hfill & \hfill V_{\mathrm{2}}(t)\hfill \end{array}\right].

Thus, standard procedures in the real time split operator can be employed to calculate Eq. (51).

As a numerical example, we investigate the ET in the Marcus inverted regime. The potentials in the fast (q) and slow (x) coordinates are modeled by two shifted harmonic oscillators, although the approaches can be straightforwardly applied to anharmonic systems. The parameters, ω = 500€cm^-1, λ_q = 24ω, λ_x = 2ω, and \mathrm{\Delta }{G}_{\mathrm{0}}=-\mathrm{40}\omega are used to make the reaction occur in the inverted regime. Here, λ_q and λ_x represent reorganization energies in the fast and slow coordinates, respectively.

The ET rate is unambiguously defined only from a single exponential decay of the population at the donor state. However, the reaction-diffusion equation Eq. (41) does not always predict population decay exponentially. Alternatively, two kinds of average survival times are introduced as

(54){\tau }_a={\displaystyle {\int }_{\mathrm{0}}^{\infty }\mathrm{d}tS(t)}

(55){\tau }_b={\displaystyle {\int }_{\mathrm{0}}^{\infty }\mathrm{d}ttS(t)/{\tau }_a}\mathrm{.}

Here S(t)={\displaystyle \int P(x,t)\mathrm{d}x} represents total population left in the reactive regime at time t. The quantity τ_a monitors the short-time dynamics of the population and is equivalent to the mean first passage time. τ_b detects information of a time region later than τ_a. Thus the correspondence rates are defined as 1/τ_a and 1/τ_b, respectively. Both the rates get equal only when the population decays single-exponentially.

To investigate the single exponential decay, one can define a time-dependent rate

(56)k_f(t)=-{\displaystyle \frac{\mathrm{d}\ln S(t)}{\mathrm{d}t}}.

If the rate constant is meaningful, k(t) will reach a plateau after a long time. This plateau value corresponds to the true rate constant.

Applications

The first example is the ET in 2,7-Dinitronaphthalene anion radical in different solvents. Nelsen et al. [79] have measured the ET rate and noted that the solvent dynamic effect is not important. Thus, they tried to use the Marcus formula and the Bixon-Jortner theory [1] in the perturbation limit to explain experimental results. However, both of them did not work because the electronic coupling was not weak enough. Since solvent dynamics is fast, here, we can use the rate expression (Eq. (15)) with Zhu-Nakamura transition probability. Indeed, the predicted rates are in excellent agreement with the experiment. The results are shown in Fig. 6. The same reorganization energy and electronic coupling estimated by Nelsen’s measurement are used in the calculations.

Nelsen et al. have also measured ET in the 2,2-dimethyl-4,4-dinitrobiphenyl (4) in several solvents [80]. The optical data reveal that this system has a very large electronic coupling (3000–4000 cm^-1). In such a case, solvent dynamics may play a significant role. Therefore, the quantum Kramers-like theory is used to estimate the rate. We use the temperature dependence of solvent relaxation time τ

(57)\ln \tau =\ln b+H/RT,

where b and H are two parameters whose values have been recommended by Grampp and Jaenicke [81]. Figure 7 shows the calculated and experimental rates, respectively. It is found that they agree with each other very well and even the slop tendency of the experimental rates with temperature can be correctly predicted.

It is more interesting to predict purely theoretically the rates without using any parameters from optical measurements. Following this clue, one may combine quantum chemistry calculation used to predict the parameters with the current approaches of obtaining the rates. We have taken [22/hex/22]⁺ as an example [48]. This system has different isomerized geometries, which challenge the experimental measurements to distinguish the individual contributions to the ET rates from the different isomers. Several approaches are proposed to calculate the parameters of controlling ET: reorganization and electronic coupling. We adopt the two methods to calculate the reorganization energy. One is the four-point technique introduced by Nelsen [82]. It can directly and easily predict the total reorganization energy. The other is based on the summarization of the reorganization energy of each mode, and it enables one to identify which mode plays a more important role in ET. While as to the calculation of the electronic coupling, several different approaches have been proposed. For instance, the energy difference of two molecular orbitals (MOs) based on Koopmans’ theorem (KT) [83–85], direct calculation of the coupling of two diabatic states (two-state model, (TM)) which represent left- and right-localized states [86,87], and the energy gap between two adiabatic states with Spin-Flip (SF) strategy [88,89]. All of these approaches have been used in our calculations for a purpose of comparison. The obtained the ET rate of oo+ [22/hex/22]⁺ conforms to the experiment. The results also reveal that the ET rate is highly dependent on the isomer geometries. For instance, the rate in ii+ [22/hex/22]⁺ is about 1000 times larger than that in oo+ [22/hex/22]⁺. Comparing with different ab initial approaches, we have confirmed that the DFT with 50-50 functional is a reasonable approach for the current system of predicting both the reorganization energy and electronic coupling. Surely, the accuracy of those parameters also depends on further strategies for calculation. In terms of the reorganization energies, both four point approach and summarizing the reorganization energy of each mode work well, but one has to be cautious about applying four point method to isomers. In the calculation of the electronic coupling, the two state model is still a good choice compared with Spin-Flip and Koopmann theorem for the current systems with strong electronic coupling.

Concluding remarks

We have presented several approaches to calculate the ET rates for weak to strong electronic coupling regime. In a fast solvent relaxation limit, the approach based on the nonadiabatic transition state theory can be adopted. It is related with Marcus formula by a prefactor and can be referred to as a modified Marcus formula. When solvent dynamics begins to play a role, the quantum Kramers-like theory is applicable. For the case in which the intramolecular vibrational motions is much faster than solvent motion, the extended Sumi-Marcus theory is a better choice. Several applications suggest that these proposed approaches, together with quantum chemistry calculations for the parameters controlled ET, can be used to predict the ET rates correctly theoretically.

References

[1] Bixon, M.; Jortner, J., Adv. Chem. Phys. 1999, 106, 35–202

[2] Levich, V. G., Adv. Electrochem. Electrochem. Eng1965, 4, 249

[3] Marcus, R. A., J. Chem. Phys. 1956, 24, 966

[4] Marcus, R. A.; Sutin, N., Biochim. Biophys. Acta1985, 811, 265–322

[5] Marcus, R. A., Rev. Mod. Phys. 1993, 65, 599–610

[6] Hush, N. S., Coord. Chem. Rev. 1985, 64, 135–157

[7] Landau, L. D., Phys. Z. Sowjetunion1932, 2, 46–51

[8] Zener, C., Proc. R. Soc. Ser. A1932, 137, 696–702

[9] Stueckelberg, E. C. G., Helv. Phys. Acta1932, 5, 369

[10] Nakamura, H., Nonadiabatic Transition: Concepts, Basic Theories and Applications; World Scientific Pub. Co. Inc., 2002

[11] Nakamura, H., J. Theor. Comput. Chem. 2005, 4, 127–137

[12] Hush, N. S., Electrochim. Acta1968, 13, 1005–1023

[13] Holstein, T., Philos. Mag. 1978, 37, 49–62

[14] Zusman, L. D., Chem. Phys. 1980, 49, 295–304

[15] Burshtein, A. I.; Kofman, A. G., Chem. Phys. 1979, 40, 289–300

[16] Barzykin, A. V.; Frantsuzov, P. A.; Seki, K.; Tachiya, M., Adv. Chem. Phys. 2002, 123, 511–616

[17] Calef, D. F.; Wolynes, P. G., J. Phys. Chem. 1983, 87, 3387–3400

[18] Hynes, J. T., J. Chem. Phys. 1986, 90, 3701–3706

[19] Rips, I.; Jortner, J., J. Chem. Phys. 1988, 88, 818

[20] Sumi, H.; Marcus, R. A., J. Chem. Phys. 1986, 84, 4894

[21] Jortner, J.; Bixon, M., J. Chem. Phys. 1988, 88, 167

[22] Walker, G. C.; Akesson, E.; Johnson, A. E.; Levinger, N. E.; Barbara, P. F., J. Phys. Chem. 1992, 96, 3728–3736

[23] Sparpaglione, M.; Mukamel, S., J. Chem. Phys. 1988, 88, 3263

[24] Jean, J. M.; Friesner, R. A.; Fleming, G. R., J. Chem. Phys. 1992, 96, 5827

[25] Topaler, M.; Makri, N., J. Phys. Chem. 1996, 100, 4430–4436

[26] Evans, D. G.; Nitzan, A.; Ratner, M. A., J. Chem. Phys. 1998, 108, 6387

[27] Wynne, K.; Hochstrasser, R. M., Adv. Chem. Phys. 1999, 107, 263–309

[28] Wang, H.; Song, X.; Chandler, D.; Miller, W. H., J. Chem. Phys. 1999, 110, 4828

[29] Thoss, M.; Wang, H.; Miller, W. H., J. Chem. Phys. 2001, 115, 2991

[30] Mühlbacher, L.; Egger, R., J. Chem. Phys. 2003, 118, 179

[31] Zhao, Y.; Li, X.; Zheng, Z.; Liang, W. Z., J. Chem. Phys. 2006, 124, 114508

[32] Zhao, Y.; Liang, W. Z., Phys. Rev. A2006, 74, 032706

[33] Berezhkovskii, A. M.; Zitserman, V. Y., Chem. Phys. Lett. 1989, 158, 369–374

[34] Coalson, R. D.; Evans, D. G.; Nitzan, A., J. Chem. Phys. 1994, 101, 436

[35] Spirina, O. B.; Cukier, R. I., J. Chem. Phys. 1996, 104, 538

[36] Cho, M. H.; Silbey, R. J., J. Chem. Phys. 1997, 106, 2654

[37] Bicout, D. J.; Szabo, A., J. Chem. Phys. 1998, 109, 2325

[38] Golosov, A. A.; Reichman, D. R., J. Chem. Phys. 2001, 115, 9848

[39] Golosov, A. A.; Frisner, R. A.; Pechukas, P., J. Chem. Phys. 2000, 112, 2095

[40] Zhao, Y.; Mil’nikov, G.; Nakamura, H., J. Chem. Phys. 2004, 121, 8854–8860

[41] Zhao, Y.; Liang, W. Z.; Nakamura, H., J. Phys. Chem. A2006, 110, 8204–8212

[42] Zhao, Y.; Nakamura, H., J. Theor. Comput. Chem. 2006, 5, 299–306

[43] Zhao, Y.; Han, M. M.; Liang, W. Z.; Nakamura, H., J. Phys. Chem. A2007, 111, 2047–2053

[44] Zhao, Y.; Mil’nikov, G., Chem. Phys. Lett. 2005, 413, 362–366

[45] Zhu, W. J.; Han, M. M.; Zhao, Y., Chinese J. Chem. Phys. 2007, 20, 217

[46] Zhu, W. J.; Zhao, Y., J. Chem. Phys. 2007, 126, 184105

[47] Zhao, Y., J. Theor. Comput. Chem. 2008, 7, 869–877

[48] Zhang, W. W.; Zhu, W. J.; Liang, W. Z.; Zhao, Y.; Nelsen, S. F., J. Phys. Chem. B2008, 112, 11079–11086

[49] Zhu, W. J.; Zhao, Y., J. Chem. Phys. 2008, 129, 184111

[50] Nelsen, S. F.; Ismagilov, R. F. Trieber, D. A. II., Science1997, 278, 846–849

[51] Demadis, K. D.; Hartshorn, C. M.; Meyer, T. J., Chem. Rev. 2001, 101, 2655–2686

[52] Ovchinnikova, M. Y., Theor. Exp. Chem. 1982, 17, 507–509

[53] Sumi, H.; Marcus, R. A., J. Chem. Phys. 1986, 84, 4272–4276

[54] Zhu, C.; Nakamura, H., J. Chem. Phys. 1994, 101, 10630

[55] Zhu, C.; Nakamura, H., J. Chem. Phys. 1995, 102, 7448

[56] Zhu, C.; Nakamura, H., Adv. Chem. Phys. 2001, 117, 127–233

[57] Pollak, E.; Grabert, H.; Hänggi, P., J. Chem. Phys. 1989, 91, 4073

[58] Rips, I. I; Pollak, E., Phys. Rev. A1990, 41, 5366–5382

[59] Hanggi, P.; Talkner, P.; Borkovec, M., Rev. Mod. Phys. 1990, 62, 251–341

[60] Griff, U.; Grabert, H.; Hänggi, P.; Riseborough, P. S., Phys. Rev. B1989, 40, 7295–7297

[61] Jortner, J., J. Chem. Phys. 1976, 64, 4860

[62] Bixon, M.; Jortner, J., J. Phys. Chem. 1986, 90, 3795–3800

[63] Miller, W. H.; Schwartz, S. D.; Tromp, J. W., J. Chem. Phys. 1983, 79, 4889

[64] Fenkel, D.; Smit, B., Understanding Molecular Simulation (Computational Science Series), Vol.1; Academy Press: USA, 2002.

[65] Rips, I.; Pollak, E., J. Chem. Phys. 1995, 103, 7912

[66] Rips, I., J. Chem. Phys. 1996, 104, 9795

[67] Rips, I., J. Chem. Phys. 2004, 121, 5356–5371

[68] Tully, J. C.; Preston, R. K., J. Chem. Phys. 1971, 55, 562–572

[69] Keck, J. C., J. Chem. Phys. 1960, 32, 1035

[70] Anderson, J. B., J. Chem. Phys. 1973, 58, 4684

[71] Chandler, D., J. Chem. Phys. 1978, 68, 2959

[72] Bergsma, J. P.; Reimers, J. R.; Wilson, K. R.; Hynes, J. T., J. Chem. Phys. 1986, 85, 5625

[73] Straub, J. E.; Borkovec, M.; Berne, B. J., J. Chem. Phys. 1988, 89, 4833

[74] Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R., J. Chem. Phys. 1990, 93, 7137

[75] Pollak, E., J. Chem. Phys. 1986, 85, 865

[76] Grabert, H., Phys. Rev. Lett. 1988, 61, 1683–1686

[77] Mel’nikov, V. I.; Meshkov, S. V., J. Chem. Phys. 1986, 85, 1018

[78] Wolynes, P. G., Phys. Rev. Lett. 1981, 47, 968–971

[79] Nelsen, S. F.; Weaver, M. N.; Konradsson, A. E.; Telo, J. P.; Clark, T., J. Am. Chem. Soc. 2004, 126, 15431–15438

[80] Telo, J. P.; Nelsen, S. F.; Zhao, Y., J. Phys. Chem. A2009, 113, 7730–7736

[81] Grampp, G.; Jaenicke, W., Ber. Bunsenges. Phys. Chem1991, 95, 904

[82] Nelsen, S. F.; Blackstock, S. C.; Kim, Y., J. Am. Chem. Soc. 1987, 109, 677–682

[83] Koopmans, T., Physica1934, 1, 104–113

[84] Paddon-Row, M. N.; Wong, S. S., Chem. Phys. Lett. 1990, 167, 432–438

[85] Jordan, K. D.; Paddon-row, M. N., Chem. Rev. 1992, 92, 395–410

[86] Farazdel, A.; Dupuis, M.; Clementi, E.; Aviram, A., J. Am. Chem. Soc. 1990, 112, 4206–4214

[87] Rodriguez-Monge, L.; Larsson, S., Int. J. Quantum Chem. 1997, 61, 847–857

[88] You, Z. Q.; Shao, Y.; Hsu, C. P., Chem. Phys. Lett. 2004, 390, 116–123

[89] Yang, C. H.; Hsu, C. P., J. Chem. Phys. 2006, 124, 244507

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg