format: text

text:
# Spectral analysis of electron transfer kinetics II

## I Introduction

Since the seminal work of Marcus on the nonadiabatic electron transfer reactions,  a great amount of effort has been made in the studies of the electron transfer reaction with a variety of tools such as time-resolved spectroscopies,  computer simulation methods, and analytical theories.  It is not surprising that there have been great interests in the studies of the electron transfer reaction, considering that it is involved in many important chemical and biological systems (for the most recent reviews, see Ref. ). For example, in the photo-synthetic reaction center, electron transfer process creates the initial charge separation which will eventually lead to the production of the adenosine triphosphate.  Also, recent studies of the molecular electronics will depend crucially on the complete understanding and controlling of the electron transfer in the chemical systems. Another kind of electron transfer reactions which are currently subject to extensive studies is proton-coupled electron transfer reactions, and many theoretical  and experimental  studies have been performed on this subject.

One of important and ubiquitous aspects of electron transfer kinetics is the dynamic solvent effect on the electron transfer rate.  As both experimental and theoretical investigations have been carried out on electron transfer reactions in solutions, many diverse phenomena including predictions of the original Marcus theory have been revealed, depending on physical parameters involved in the electron transfer kinetics.

Electronic coupling constant, $`V`$, given by the interaction matrix element between the electron donor and acceptor wavefunctions is one of the most important physical parameters, and depending on its magnitude compared with other parameters, various kinetic regimes are exhibited in the electron transfer process. When the electronic coupling constant is the smallest parameter of the electron transfer process, the electron transfer rate is well described by perturbation theory, which predicts the golden-rule rate,

$`kk_{\mathrm{GR}}={\displaystyle \frac{2\pi V^2}{\mathrm{}}}\rho _c,`$ (1)
where $`\rho _c`$ is the equilibrium population of the reactant state in the crossing regime, and this is a well known result of the Marcus electron transfer theory. When the electronic coupling constant large enough, the overall reaction process is not determined by the Marcus rate, Eq. (1), but by the solvent diffusion rate describing the polarization dynamics of the solvent molecules,

$`kk_\mathrm{D}\mathrm{\Omega }\lambda \rho _c,`$ (2)
where $`\mathrm{\Omega }`$ is the solvent relaxation rate and $`\lambda `$ is the classical reorganization energy, and this case is called the solvent-controlled limit.

About two decades ago, Zusman investigated the crossover between the Marcus and solvent-controlled regimes in the studies of the electron transfer reaction in Debye solvents. He used a mixed quantum-classical approach where the stochastic operator is introduced to describe bath relaxation processes occurring on two diabatic surfaces and treated the electron transfer process as a nonadiabatic transition between them. Zusman solved the nonadiabatic diffusion equation in the weak coupling limit $`(Vk_BT)`$, and obtained the expression for the overall electron transfer rate,

$`k^1=k_{\mathrm{GR}}^1+k_\mathrm{D}^1,`$ (3)
which shows a transition from the Marcus to the solvent-controlled limits in the nonadiabatic regime.

As the electronic coupling constant is increased further to be comparable to or larger than the thermal energy, it is expected that the electron transfer process involves an adiabatic barrier crossing. Finally, when the electronic coupling constant is even larger so that it has the same order of magnitude as the solvent reorganization energy, $`\lambda `$, which is indeed the case for the mixed valence compounds, the electronic states are delocalized on the lower adiabatic surface. Due to the delocalization nature of electronic states, an adiabatic picture is more useful than the diabatic one for analyzing the short-time dynamics in strongly coupled systems. In this picture, electronic coherence arises from Rabi oscillation between two adiabatic surfaces.

Although there have been several studies to bridge between the Marcus regime and the solvent-controlled regime using various approaches,  few studies have discussed diverse kinetic regimes in a unified way, and often different approaches are taken in different regimes. It is desirable to investigate the effects of solvent dynamics on the electron transfer process in a unified approach for various parameter regimes.

As a general approach to describing condensed phase dynamics, we recently proposed a spectral analysis method. Instead of focusing on dynamical trajectories of the reduced density matrix for dissipative systems, this methodology investigates the spectral structure of the evolution operator for dissipative systems, and it has been applied to the electron transfer process in mixed valence compounds to investigate the possibility of electronic coherence in those systems.

In this paper we present a thorough analysis of the electron transfer kinetics in Debye solvents based on the spectral analysis method. Electron transfer rate constant extracted from the spectral analysis is compared with other previous results both in the nonadiabatic and adiabatic regimes, and transition from the incoherent to coherent regimes is demonstrated by the spectral analysis method. When the solvent relaxation rate is very fast, it is found that the solvent dynamics has a significant effect on the Marcus curve. The spectral analysis method is also utilized as a density matrix propagation scheme. Preliminary results of the spectral analysis method focusing on symmetric reaction cases have been reported.

The rest of the paper is organized as follows. The spectral analysis method of the nonadiabatic diffusion equation is formulated in Sec. II. Two important limiting cases are discussed in Sec. III. In Sec. IV comprehensive analysis of the spectral structure of the nonadiabatic diffusion equation is performed for a broad range of parameters, and diverse kinetic behaviors in electron transfer reactions are identified and characterized. We conclude in Sec. V by summarizing results obtained in this work.

## II Theory

We consider two electronic states, $`|1`$ and $`|2`$, which represent the electron donor and acceptor sites of the electron transfer system, respectively, and they are coupled to each other via the electronic coupling matrix element, $`V`$. Moreover, each electronic state is coupled to bath degrees of freedom.

There have been extensive studies of the solvent effect on electron transfer dynamics in literature with various approaches.  Considering that electron transfer processes are usually probed at room temperature in polar solvents, one can treat the bath degrees of freedom classically. Zusman and Yakobson-Burshtein  proposed a mixed quantum-classical evolution equation of the reduced density matrix, $`\rho (E,t)`$, independently, to investigate the solvent effect on electron transfer,

$`{\displaystyle \frac{}{t}}\rho (E,t)=\rho (E,t)=(_B+i𝒱)\rho (E,t).`$ (4)
Here, $`E`$ is the solvent polarization energy which plays a role of the reaction coordinate as first noticed by Marcus. $``$ and $`𝒱`$ represent operators for the solvent relaxation dynamics and for the electronic transition between two states, respectively. Explicitly, Eq. (4) is written in terms of the density matrix elements,

$`\dot{\rho }_{11}`$ $`=`$ $`_{11}\rho _{11}+iV(\rho _{12}\rho _{21}),`$ (5)
$`\dot{\rho }_{22}`$ $`=`$ $`_{22}\rho _{22}iV(\rho _{12}\rho _{21}),`$ (6)
$`\dot{\rho }_{12}`$ $`=`$ $`_{12}\rho _{12}i\omega _{12}\rho _{12}+i{\displaystyle \frac{V}{\mathrm{}}}(\rho _{11}\rho _{22}),`$ (7)
$`\dot{\rho }_{21}`$ $`=`$ $`_{21}\rho _{21}+i\omega _{12}\rho _{21}i{\displaystyle \frac{V}{\mathrm{}}}(\rho _{11}\rho _{22}),`$ (8)
Here, the diagonal and off-diagonal matrix elements of the reduced density matrix $`\rho (E,t)`$ represent populations of the electronic states and coherences between them, respectively, and $`_{ij}`$’s describe the relaxation process of classical bath over the free energy surfaces, with $`_{ii}`$ defined on the free energy surface for the $`i`$th electronic state, and $`_{12}`$ and $`_{21}`$ defined on the averaged free energy surface. The functional form for the free energy surface in the electron transfer system is usually harmonic,

$`U_1(E)`$ $`=`$ $`{\displaystyle \frac{(E+\lambda )^2}{4\lambda }},`$ (9)
$`U_2(E)`$ $`=`$ $`{\displaystyle \frac{(E\lambda )^2}{4\lambda }}+ϵ,`$ (10)
where $`\lambda `$ is the reorganization energy, and we assume $`ϵ<0`$ without loss of generality. It is convenient to define $`\overline{U}`$ and $`\mathrm{}\omega _{12}`$ are the average and the difference of the two free energy surfaces, respectively,

$`\overline{U}(E)={\displaystyle \frac{U_1(E)+U_2(E)}{2}}={\displaystyle \frac{E^2+\lambda ^2}{4\lambda }}+{\displaystyle \frac{ϵ}{2}},`$ (11)
$`\mathrm{}\omega _{12}(E)=U_1(E)U_2(E)=Eϵ.`$ (12)
We set $`\mathrm{}=1`$ for simplicity henceforth. This set of a mixed quantum-classical two-state equation has been previously derived by several authors starting from the spin-boson Hamiltonian.

We note that many chemically and biologically important electron transfer processes take place in an overdamped solvent environment. Then, the bath relaxation operators in Eqs. (5)-(8) are modeled by one-dimensional Fokker-Planck operators $`_{ij}`$,

$`_{ii}=D_E{\displaystyle \frac{}{E}}\left({\displaystyle \frac{}{E}}+\beta {\displaystyle \frac{U_i(E)}{E}}\right),`$ (13)
$`_{12}=_{21}={\displaystyle \frac{_{11}+_{22}}{2}}=D_E{\displaystyle \frac{}{E}}\left({\displaystyle \frac{}{E}}+\beta {\displaystyle \frac{\overline{U}(E)}{E}}\right),`$ (14)
where $`\beta =1/(k_BT)`$.

The Fokker-Planck equation models the relaxation process of the density matrix element as a diffusion process in the energy space and various parameters are identified; energy diffusion constant, $`D_\mathrm{E}=\mathrm{\Omega }\mathrm{\Delta }^2`$, fluctuation of the solvent polarization energy, $`\mathrm{\Delta }^2=E^2=2\lambda k_\mathrm{B}T`$, and characteristic timescale of a Debye solvent $`\tau _\mathrm{D}=1/\mathrm{\Omega }`$. The correlation function of the solvent polarization energy is given by a single exponential form in a Debye solvent,

$`C(t)=E(t)E(0)=\mathrm{\Delta }^2\mathrm{exp}(\mathrm{\Omega }t).`$ (15)
Note that since the nuclear dynamics is modeled by the Fokker-Planck operator, the possibility of the vibrational coherence is not considered in this model of electron transfer dynamics.

We investigate the spectral structure of the nonadiabatic diffusion operator, $``$, by calculating eigenvalues, $`\{Z_\nu \}`$, and corresponding right and left eigenfunctions, $`\{|\psi _\nu ^R\}`$ and $`\{\psi _\nu ^L|\}`$,

$`|\psi _\nu ^R`$ $`=`$ $`Z_\nu |\psi _\nu ^R,`$ (16)
$`\psi _\nu ^L|`$ $`=`$ $`Z_\nu \psi _\nu ^L|.`$ (17)
Because the nonadiabatic Liouville operator is non-Hermitian, the eigenvalues are generally given by complex values, and the right and left eigenfunctions corresponding to the same eigenvalue are not simply the Hermitian conjugate to each other.

The method of eigenfunction solution is well known for the diffusion process on the harmonic potential energy surface as discussed in App. A.  Unlike the diffusion problem on the single potential energy surface, however, there have been limited studies on the nonadiabatic diffusion problem involving more than a single potential energy surface. Cukier and co-workers have calculated the electron transfer rate by calculating the lowest eigenvalue of the nonadiabatic diffusion equation in the weak-coupling regime.

In this paper, the eigenfunctions of $`_{12}`$ are used as our basis set to represent the nonadiabatic diffusion equation. In principle, one could have chosen the eigenfunctions of $`_{11}`$ or $`_{22}`$ as basis functions; however, in that case one has to evaluate appropriate Franck-Condon factors when calculating the coupling matrix elements. The Fokker-Planck operator $`_{12}`$ is defined on the averaged single harmonic potential centered at $`E=0`$, and its eigenvalue solutions are obtained following a similar procedure in App. A,

$`_{12}|\varphi _n^R`$ $`=`$ $`n\mathrm{\Omega }|\varphi _n^R,`$ (18)
$`\varphi _n^L|_{12}`$ $`=`$ $`n\mathrm{\Omega }\varphi _n^L|,`$ (19)
where $`n=0,1,2,\mathrm{}`$, and the $`n`$th right and left eigenfunctions for $`_{12}`$ are given by

$`\varphi _n^R(E)`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}\left(\frac{E^2}{2\mathrm{\Delta }^2}\right)}{(2^nn!)^{\frac{1}{2}}(2\pi \mathrm{\Delta }^2)^{\frac{1}{4}}}}H_n\left({\displaystyle \frac{E}{\sqrt{2}\mathrm{\Delta }}}\right),`$ (20)
$`\varphi _n^L(E)`$ $`=`$ $`{\displaystyle \frac{1}{(2^nn!)^{\frac{1}{2}}(2\pi \mathrm{\Delta }^2)^{\frac{1}{4}}}}H_n\left({\displaystyle \frac{E}{\sqrt{2}\mathrm{\Delta }}}\right),`$ (21)
with $`H_n`$ being the $`n`$th order Hermite polynomial.

We separate real and imaginary parts of the coherence density matrix, namely, $`u=`$Re$`\rho _{12}`$ and $`v=`$Im$`\rho _{12}`$, and rewrite Eqs. (5)-(8) as

$`\dot{\rho }_{11}`$ $`=`$ $`(_{12}+\delta )\rho _{11}2Vv,`$ (22)
$`\dot{\rho }_{22}`$ $`=`$ $`(_{12}\delta )\rho _{22}+2Vv,`$ (23)
$`\dot{u}`$ $`=`$ $`_{12}u+\omega _{12}v,`$ (24)
$`\dot{v}`$ $`=`$ $`_{12}v\omega _{12}u+V(\rho _{11}\rho _{22}),`$ (25)
where $`\delta =(_{11}_{22})/2`$. Then, all the relevant operators in Eqs. (22)-(25) can be evaluated in terms of the basis functions,

$`\varphi _n^L|_{12}|\varphi _m^R`$ $`=`$ $`n\mathrm{\Omega }\delta _{nm},`$ (26)
$`\varphi _n^L|\delta |\varphi _m^R`$ $`=`$ $`\mathrm{\Omega }\sqrt{{\displaystyle \frac{\lambda }{2k_BT}}}\sqrt{m+1}\delta _{n,m+1},`$ (27)
$`\varphi _n^L|\omega _{12}|\varphi _m^R`$ $`=`$ $`\sqrt{2\lambda k_BT}(\sqrt{m}\delta _{n,m1}+\sqrt{m+1}\delta _{n,m+1})ϵ\delta _{nm},`$ (28)
$`\varphi _n^L|V|\varphi _m^R`$ $`=`$ $`V\delta _{nm}.`$ (29)
With this basis set, we can expand the density matrix elements as

$`\rho _{11}(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}a_n(t)\varphi _n^R(E),`$ (30)
$`\rho _{22}(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}b_n(t)\varphi _n^R(E),`$ (31)
$`u(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}c_n(t)\varphi _n^R(E),`$ (32)
$`v(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}d_n(t)\varphi _n^R(E).`$ (33)
Substituting Eqs. (30)-(33) into the right eigenvalue equation Eq. (16), we have the following coupled eigenvalue equations,

$`Z_\nu a_n`$ $`=`$ $`n\mathrm{\Omega }a_n\mathrm{\Omega }\sqrt{{\displaystyle \frac{\lambda }{2k_BT}}}\sqrt{n}a_{n1}2Vd_n,`$ (34)
$`Z_\nu b_n`$ $`=`$ $`n\mathrm{\Omega }b_n+\mathrm{\Omega }\sqrt{{\displaystyle \frac{\lambda }{2k_BT}}}\sqrt{n}b_{n1}+2Vd_n,`$ (35)
$`Z_\nu c_n`$ $`=`$ $`n\mathrm{\Omega }c_n+\sqrt{2\lambda k_BT}(\sqrt{n+1}d_{n+1}+\sqrt{n}d_{n1})ϵd_n,`$ (36)
$`Z_\nu d_n`$ $`=`$ $`n\mathrm{\Omega }d_n\sqrt{2\lambda k_BT}(\sqrt{n+1}c_{n+1}+\sqrt{n}c_{n1})+ϵc_n+V(a_nb_n),`$ (37)
which is an explicit basis set representation for the two-state diffusion operator in Eqs. (5)-(8). Eigenvalue equations for the left eigenvector in Eq. (17) can be written by making the transpose of Eqs. (34)-(37). Diagonalizing the $`4N\times 4N`$ matrix ($`N`$ = number of basis functions) defined in Eqs. (34)-(37), we obtain the eigenvalues $`Z_\nu `$ and the corresponding eigenvectors of the nonadiabatic diffusion operator,

$`|\psi _\nu ^R`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}R_{n\nu }|\varphi _n^R,`$ (38)
$`\psi _\nu ^L|`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}L_{\nu n}\varphi _n^L|,`$ (39)
where $`R_{n\nu }`$ and $`L_{\nu n}`$ are elements of the transformation matrices.

In general, due to the non-Hermitian nature of the nonadiabatic diffusion operator, the left and right eigenfunctions do not form an orthogonal set by themselves. However, when the eigenvalues are all nondegenerate, the left and right eigenfunctions form an orthogonal and complete set in a dual Hilbert space. Explicitly, we have

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}L_{\nu n}R_{n\nu ^{^{}}}`$ $`=`$ $`\delta _{\nu \nu ^{^{}}},`$ (40)
$`{\displaystyle \underset{\nu =0}{\overset{\mathrm{}}{}}}R_{n\nu }L_{\nu m}`$ $`=`$ $`\delta _{nm},`$ (41)
for the orthogonality and the completeness relation, respectively. Using these properties, we can construct the real time propagator for the operator $``$ as

$`G(t)={\displaystyle \underset{\nu =0}{\overset{\mathrm{}}{}}}|\psi _\nu ^R\psi _\nu ^L|e^{Z_\nu t},`$ (42)
and express the time evolution of the density matrix as an eigenfunction expansion,

$`|\rho (t)=G(t)|\rho (0)={\displaystyle \underset{\nu =0}{\overset{\mathrm{}}{}}}|\psi _\nu ^R\psi _\nu ^L|\rho (0)e^{Z_\nu t}.`$ (43)
Note that the right and left eigenfunctions play asymmetric roles in the construction of the propagator for the non-Hermitian operator.

## III Limiting Cases

In order to compare the eigenvalue solution developed in this work with previous theoretical predictions available in different kinetic regimes, we briefly discuss two limiting cases which have been studied extensively in the literature. Instead of giving detailed derivations, we will briefly mention the solutions in these limiting cases, relegating details to the App. B.

### A Nonadiabatic Regime : Weak Coupling Case

When the electronic coupling matrix element $`V`$ is very small, we can reduce the full nonadiabatic diffusion equation, Eqs. (5)-(8), into the population evolution equation, and detailed derivations can be found in Refs.  (see also App. B). Kinetic equations for the population elements of $`\rho `$ can be approximately written as

$`{\displaystyle \frac{\rho _{11}}{t}}`$ $`=`$ $`_{11}\rho _{11}K(\rho _{11}\rho _{22}),`$ (44)
$`{\displaystyle \frac{\rho _{22}}{t}}`$ $`=`$ $`_{22}\rho _{22}K(\rho _{22}\rho _{11}),`$ (45)
where $`K(E)`$ is the rate kernel given in Eq. (B7). A dynamical quantity usually measured in the electron transfer kinetics experiment is the total population, $`P_i(t)`$, in each electronic state rather than the polarization energy dependent population, $`\rho _{ii}(E,t)`$,

$`P_i(t)={\displaystyle _{\mathrm{}}^{\mathrm{}}}dE\rho _i(E,t).`$ (46)
Using the projection operator method and making a time-scale separation approximation, it can be shown that the kinetic process between two electronic states is described by the time-independent rate constant instead of the rate kernel in the weak coupling limit, (see also App. B). Then the kinetic equation for $`P_i`$ is given by

$`\left(\begin{array}{c}\dot{P}_1(t)\\ \dot{P}_2(t)\end{array}\right)=\left(\begin{array}{cc}k^1& k^2\\ k^1& k^2\end{array}\right)\left(\begin{array}{c}P_1(t)\\ P_2(t)\end{array}\right).`$ (53)
Here $`k^1`$ and $`k^2`$ are forward and backward rate constants and they are related to the nonadiabatic transition rates $`k_{\mathrm{NA}}^i`$ and solvent diffusion rates $`k_\mathrm{D}^i`$ by

$`k^i`$ $`=`$ $`{\displaystyle \frac{k_{\mathrm{NA}}^i}{1+k_{\mathrm{NA}}^1/k_\mathrm{D}^1+k_{\mathrm{NA}}^2/k_\mathrm{D}^2}},`$ (54)
$`k_{\mathrm{NA}}`$ $`=`$ $`k_{\mathrm{NA}}^1+k_{\mathrm{NA}}^2,`$ (55)
$`k_{\mathrm{ND}}`$ $`=`$ $`k^1+k^2={\displaystyle \frac{k_{\mathrm{NA}}}{1+k_{\mathrm{NA}}^1/k_\mathrm{D}^1+k_{\mathrm{NA}}^2/k_\mathrm{D}^2}},`$ (56)
Recently, Eq. (56) has also been obtained for the symmetric reaction case in the weak coupling limit by calculating the first excited eigenvalue explicitly via the Goldstone theorem.

Eq. (56) has a form of an overall relaxation rate for consecutive reactions, and it involves two different types of rate processes; nonadiabatic transition and solvent diffusion rate. Nonadiabatic transition rate constants $`k_{\mathrm{NA}}^1`$ and $`k_{\mathrm{NA2}}^2`$ are forward and backward quantum transition rates between the two electronic states and are calculated in terms of the coherent Green’s functions, $`G_{12}(E,t|E_0)`$ and $`G_{21}(E,t|E_0)`$(see also App. B),

$`k_{\mathrm{NA}}^i=2\pi V^2\mathrm{Re}{\displaystyle _0^{\mathrm{}}}dte^{i(ϵ+\lambda )tg(t)},`$ (57)
where the bath correlation function $`g(t)`$ is given by

$`g(t)=\left[\left({\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Omega }}}\right)^2i\pm i{\displaystyle \frac{\lambda }{\mathrm{\Omega }}}\right]\left[\mathrm{exp}(\mathrm{\Omega }t)+\mathrm{\Omega }t1\right].`$ (58)
with $`+()`$ sign for $`i=1(2)`$. When the bath dynamics is slow such that $`\mathrm{\Omega }\mathrm{\Delta }`$, Eq. (57) reduces to the standard Marcus result(see App. B),

$`k_{\mathrm{NA}}^i`$ $``$ $`k_{\mathrm{GR}}^i,=2\pi V^2\rho _i^{\mathrm{eq}}(ϵ),`$ (59)
$`k_{\mathrm{NA}}`$ $``$ $`k_{\mathrm{GR}}^1+k_{\mathrm{GR}}^2.`$ (60)
where $`\rho _i^{\mathrm{eq}}(ϵ)=e^{(\lambda \pm ϵ)^2/2\mathrm{\Delta }^2}/\sqrt{2\pi \mathrm{\Delta }^2}`$ is the equilibrium population of the $`i`$th state at the crossing point $`E=ϵ`$. When $`\mathrm{\Omega }\mathrm{\Delta }`$, the nonadiabatic transition rate shows a deviation from the Marcus rate, which will be discussed later in details.

Solvent diffusion rate constants, $`k_\mathrm{D}^1`$ and $`k_\mathrm{D}^2`$, describe solvent relaxation processes that equilibrate the nonequilibrium wavepacket created at the crossing region, $`E_c=ϵ`$, during the electron transfer process. They can be written in terms of the population Green’s functions $`G_{11}(E,t|E_0)`$ and $`G_{22}(E,t|E_0)`$,

$`(k_\mathrm{D}^i)^1`$ $`=`$ $`{\displaystyle \frac{1}{\rho _{ii}^{\mathrm{eq}}(ϵ)}}{\displaystyle _0^{\mathrm{}}}dt[G_{ii}(ϵ,t|ϵ)G_{ii}(ϵ,\mathrm{}|ϵ)].`$ (61)
Noting that $`G_{ii}(ϵ,\mathrm{}|ϵ)=\rho _i^{\mathrm{eq}}(ϵ)`$ is given by the equilibrium population distribution of the $`i`$th state at $`E=ϵ`$, the solvent diffusion rates are identified as the inverse of the mean survival time of the relative nonequilibrium population created at the crossing point, $`E_c=ϵ`$.

The boundary between the Marcus regime and the solvent-controlled regime in the nonadiabatic limit is estimated when comparing the Marcus rate in Eq. (59), and the solvent-diffusion rate in Eq. (B19), which leads to

$`2\pi V^2|ϵ\pm \lambda |\mathrm{\Omega }.`$ (62)
The nonadiabatic regime is established when the timescale of the off-diagonal density matrix is much faster than that of the diagonal ones, yielding the following condition,

$`VD_\mathrm{E}^{1/3}=(2\lambda k_BT\mathrm{\Omega })^{1/3}.`$ (63)
in addition to the weak coupling condition, $`Vk_BT`$, for the validity of the nonadiabatic regime.

### B Adiabatic Regime : Strong Coupling Case

When the electronic coupling constant is comparable to or larger than the thermal energy $`\beta V1`$, it is more adequate to describe electron transfer in an adiabatic representation. Given a specific bath configuration characterized by the value of $`E`$, the electronic part of the Hamiltonian in the diabatic representation is given by

$`H_\mathrm{d}(E)=\left(\begin{array}{cc}U_1(E)& V\\ V& U_2(E)\end{array}\right).`$ (66)
After the diagonalization of $`H_\mathrm{d}(E)`$ for a given bath configuration $`E`$, we can obtain two adiabatic surfaces,

$`U_\pm (E)`$ $`=`$ $`\overline{U}(E)\pm {\displaystyle \frac{1}{2}}\left[\omega _{12}(E)^2+4V^2\right]^{1/2}.`$ (67)
The separation between two adiabatic surfaces is much larger than $`k_BT`$, $`U_+U_{}Vk_BT`$. When the energy bias $`ϵ`$ and the electronic coupling $`V`$ are not so large compared with the reorganization energy, $`\lambda `$, the lower adiabatic surface, $`U_{}`$, has a well-defined double well structure. In this case, we can consider the electron transfer process as a diffusional barrier crossing process occurring on the lower adiabatic surface $`U_{\mathrm{AD}}=U_{}`$ described by the following diffusion equation in the energy space,

$`{\displaystyle \frac{\rho }{t}}=D_\mathrm{E}\left({\displaystyle \frac{^2}{E^2}}+\beta {\displaystyle \frac{}{E}}U_{\mathrm{AD}}^{}\right)\rho =_{\mathrm{AD}}\rho .`$ (68)
To calculate the reaction rate for the adiabatic barrier crossing described by Eq. (68), one may take several different routes. For example, one can calculate the barrier crossing rate as the first nonzero eigenvalue of the Fokker-Planck operator, $`_{\mathrm{AD}}`$, corresponding to the lower adiabatic surface,

$`_{\mathrm{AD}}\psi =Z\psi k_{\mathrm{AD}}^{\mathrm{ev}}=Z_1`$ (69)
Another popular approach is a population-over-flux method, which obtains the barrier crossing rate by calculating the steady-state flux at the crossing region, $`E_c=ϵ`$, divided by equilibrium population in the reactant region,

$`k_{\mathrm{AD}}^{\mathrm{fp}}=D_\mathrm{E}\left[{\displaystyle _{\mathrm{}}^{E_c}}dE_1{\displaystyle _{\mathrm{}}^{E_1}}dE_2e^{\beta \left[U_{}(E_1)U_{}(E_2)\right]}\right]^1`$ (70)
In the symmetric reaction case ($`ϵ=0`$), the adiabatic rate given in Eq. (70) can be evaluated approximately in the strong coupling limit ($`\beta V1`$) by expanding the integrand at $`E=0`$ and performing a Gaussian integration approximately to yield

$`k_{\mathrm{AD}}^{\mathrm{fp}}{\displaystyle \frac{\mathrm{\Omega }}{\pi }}\sqrt{{\displaystyle \frac{\lambda }{2V}}}e^{\beta (\frac{\lambda }{4}V)},`$ (71)
which is the Kramer’s adiabatic reaction rate in the strong damping regime.

## IV Results of Spectral Analysis

### A General Feature of Spectra

Figure 1 shows the eigenvalue spectra in three different cases of the energy bias for various electronic couplings. Since the nonadiabatic diffusion operator is not Hermitian, the eigenvalues are complex in general,

$`\psi _\nu =Z_\nu \psi _\nu =(Z_\nu ^{}+iZ_\nu ^{\prime \prime })\psi _\nu .`$ (72)
In order for the density matrix relax to the stationary state at long times, the real part of the eigenvalue should satisfy $`Z_\nu ^{}0`$. Eigenvalues of the nonadiabatic diffusion operator are classified into three different cases with their associated kinetic behaviors in the density matrix evolution;

(1) $`Z_0^{}=Z_0^{\prime \prime }=0`$ : equilibrium,

(2) $`Z_\nu ^{}>0,Z_\nu ^{\prime \prime }=0`$ : exponential decay,

(3) $`Z_\nu ^{}>0,Z_\nu ^{\prime \prime }0`$ : damped oscillation.
The eigenstate with zero eigenvalue \[case (1)\] corresponds to the equilibrium solution of the density matrix while those with real \[case (2)\] and complex eigenvalues \[case (3)\] correspond to exponential decays and to damped oscillations in the evolution of the density matrix, respectively. When the first nonzero real eigenvalue $`Z_1^{}`$ is well separated from the other higher eigenvalues, $`Z_1^{}Z_{\nu 2}^{}`$, the dynamics of the density matrix will be described by a single exponential process, and the relaxation rate is well defined as $`k=Z_1^{}`$. However, if this is not the case, the dynamics of the density matrix will generally involve multiple timescales and coherent oscillations.

We used 800 basis functions to calculate the eigenvalues, and to remove the finite basis set effect at higher eigenvalues, we present only lower 400 eigenvalues in Fig. 1. Three different cases of the energy bias chosen in Fig. 1(a)-(c) correspond to different regimes in nonadiabatic electron transfer theory, i.e., normal ($`|ϵ|<\lambda `$ in Fig. 1(a)), activationless ($`|ϵ|=\lambda `$ in Fig. 1(b)), and inverted regimes ($`|ϵ|>\lambda `$ in Fig. 1(c)), respectively. As a general feature, all the spectra show tree-like structures with three major characteristic branches; one branch for the real eigenvalues in the middle, real axis, and two other branches for the complex conjugate eigenvalues. Therefore, the spectral structure indicates that multiple exponential decays (real eigenvalues) and damped oscillations (complex eigenvalues) are inherent features in electron transfer processes in the overdamped solvent. In order to infer the dynamical behavior of the density matrix, we only need to focus on the lower part of the eigenvalue diagrams in Fig. 1 corresponding to the eigenstates with $`Z^{}/\mathrm{\Omega }1`$ since those with $`Z^{}/\mathrm{\Omega }1`$ will decay out very quickly.

(i) symmetric reaction case($`ϵ=0`$) : In Fig. 1(a), one can notice that the real parts of the eigenvalues, $`Z_\nu ^{}`$’s, located at the branching regime in the eigenvalue tree decrease as the coupling constant increases up to $`\beta V1`$, which means that the oscillatory components will persist during a more extended period $`t1/Z_\nu ^{}`$ as the coupling constant increases. When the coupling constant increases further, $`\beta V1`$, the branching point of the complex branches start to separates along the imaginary axis and the separation between the two complex branches along the imaginary axis is given by $`2V`$, and this corresponds to the Rabi oscillation frequency for the two adiabatic states in the strong coupling regime

(ii) asymmetric reaction cases($`ϵ0`$) : In Figs. 1(b) ($`\beta |ϵ|=10`$) and (c) ($`\beta |ϵ|=20`$), we notice that spectral structures show more branched behaviors in the asymmetric reaction cases than in the symmetric reaction case in Fig. 1(a). When comparing the zero electronic coupling cases, $`\beta V=0`$, in Figs. 1(a), (b), and (c), we notice that in the asymmetric reaction cases two complex branches separate out along the imaginary axis by the amount of the energy bias, $`ϵ`$ unlike in the symmetric reaction case. This corresponds to the natural oscillation frequency of the off-diagonal density matrix elements.

In Fig. 2 we have followed the evolution of the lowest 30 eigenvalues in both the real and imaginary axes as a function of the coupling constant for different energy bias cases. In all the three cases real parts of the eigenvalues start as pairs of doubly degenerate states separated from each other by the solvent relaxation rate, $`\beta \mathrm{\Omega }=1`$, when $`V=0`$. As the coupling constant increases, the degeneracies of the real eigenvalues are first lifted at small coupling constants, and bifurcation and coalescence of eigenvalues occur at large coupling constants. The coalescence corresponds to the case where two closely separated real eigenvalues become complex conjugate with the same $`Z^{}`$, and this indicates the transition from the incoherent, rate process to coherent oscillatory behavior in the dynamics of the density matrix, which will be further discussed in Sec. IV C. In the imaginary part diagram, we find an interesting behavior that imaginary parts of the eigenvalues asymptotically follow the Rabi oscillation frequency,

$`Z^{\prime \prime }\omega _{\mathrm{Rabi}}=\sqrt{4V^2+ϵ^2},`$ (73)
in the strong coupling limit, and this demonstrates the signature of adiabatic Rabi oscillation picture of the electron transfer in the strong electronic coupling case.

### B Crossover from Nonadiabatic to Adiabatic Regimes

In Fig. 3 we compared the first nonzero eigenvalue of the nonadiabatic diffusion operator with several predictions from existing theories in order to characterize various kinetic regimes in the nonadiabatic diffusion equation. The filled circles are calculation results of the first nonzero eigenvalue obtained from the spectral analysis method. We compare our eigenvalue solution with previous theoretical results, each of which is applicable for different limits;

(i) weak coupling limit : The eigenvalue solution is compared with the purely nonadiabatic transition rate without the solvent diffusion effects, Eq. (55) (solid line) and with the nonadiabatic-diffusion rate with the solvent diffusion effects included, Eq. (56) (long dashed line). In any case, we compare the eigenvalue solution with the sum of the forward and backward rates. since the eigenvalue solution gives the overall relaxation rate of the density matrix.

(ii) strong coupling limit : We compare the eigenvalue solution with the adiabatic barrier crossing rates as discussed in Sec. III B. We calculated two different adiabatic reaction rates, the lowest nonzero eigenvalue given in Eq. (69) (dot-dashed line) and the flux-over-population rate given in Eq. (70) (short dashed line), respectively.

In the slow bath relaxation case ($`\beta \mathrm{\Omega }=0.05`$) given in Fig. 3(a), when the the coupling constant is very small ($`\beta V0.1`$) (Marcus regime), the eigenvalue solution agrees well both with the nonadiabatic transition rate and with the nonadiabatic-diffusion rate. The nonadiabatic transition is the rate limiting step in this case. As the coupling constant becomes large and so does the nonadiabatic transition rate, the overall rate is now affected by the solvent relaxation rate, and the eigenvalue solution follows the nonadiabatic-diffusion rate, Eq. (56), and this is the solvent-controlled regime.

As the coupling constant becomes much larger than the thermal energy, the eigenvalue solution shows a saturation behavior at $`\beta V1`$ first, and starts to increase rapidly as the electronic coupling becomes stronger, and this indicates the adiabatic barrier crossing in the strong coupling regime. In this case, $`\beta V1`$, the eigenvalue shows a qualitative agreement with the adiabatic reaction rates, Eq. (69) and Eq. (70).

As we increase the solvent relaxation rate $`(\beta \mathrm{\Omega }=1)`$ in Fig. 3(b), the distinction between the nonadiabatic and the adiabatic regime in the eigenvalue solution becomes less obvious than in the slow bath case, Fig. 3(a). In the fast bath relaxation case, as the nonadiabatic transition rate becomes large to be in the solvent-controlled regime, the adiabatic barrier crossing is already appreciable, so the transition between the solvent-controlled and adiabatic barrier crossing regimes is not so clear. This is related to a dynamical modulation of the nonadiabatic transition rate in the fast bath case, and will be discussed in details in Sec. IV D.

The agreement between the eigenvalue solution and the two adiabatic rate calculations is qualitative in the strong coupling case. Although it is reasonable to treat the electron transfer process as an adiabatic barrier crossing reaction on the lower adiabatic surface when the coupling constant is large $`(\beta V1)`$, it should be mentioned that a rigorous mathematical proof has not been given to show that the electron transfer rate from the nonadiabatic diffusion equation really corresponds to the adiabatic reaction rate in the strong coupling limit. Therefore, it may not be surprising that the agreement between the eigenvalue solution and the adiabatic solutions is only qualitative. Related to this problem, it should be mentioned that the original nonadiabatic diffusion equation is based on diabatic representations of electronic states, therefore it may not give quantitatively the same rate constant as that from the adiabatic diffusion equation in the strong coupling case.

The transition from the Marcus to solvent-controlled regimes has been studied theoretically by using different theoretical methods such as the projection operator techniques and the path-integral methods. These theoretical approaches do not cover the adiabatic barrier crossing regime. The crossover behavior from the nonadiabatic to adiabatic regimes has been observed in quantum mechanical approaches to electron transfer theories based on the instanton solution and the diagrammatic technique which do not take into account the solvent dynamics effects. The eigenvalue solution presented here clearly demonstrates transitions between three different regimes.

### C Transition from Incoherent to Coherent Regimes

For the symmetric electron transfer reaction case$`(ϵ=0)`$ the lower adiabatic surface, $`U_{}`$, has a double well structure if $`V<\lambda /2`$. When the electronic coupling constant becomes even larger than $`\lambda /2`$, the lower adiabatic surface becomes a single well without any barrier. We consider this situation as a thermodynamic transition from the localized to delocalized electronic states when viewed from a perspective of the diabatic states. Even though the free energy surfaces of the donor and acceptor states support the large energy barrier between them in the diabatic picture, a large coupling constant in the same order of magnitude as the reorganization energy makes the distinction between the donor and acceptor states inappropriate. Instead, two new single well potential surfaces are obtained, whose eigenfunctions are neither donor nor acceptor wavefunctions, and rather linear combinations of them. This situation occurs in mixed valence compounds and other strongly coupled systems where $`V\lambda k_BT`$.  The critical value of $`V=V_c`$ at which the double well structure disappears in the lower adiabatic surface, $`U_{}`$ can be obtained from Eq. (67) when the reorganization energy is $`\lambda `$ and the energy bias is $`|ϵ|`$,

$`(2V_c)^{2/3}+|ϵ|^{2/3}=\lambda ^{2/3},`$ (74)
which indicates that the localization-delocalization transition occurs at a smaller value of $`V_c<\lambda /2`$ in the asymmetric reaction case. Recently, Cukier and co-workers, and Hänggi and co-workers have studied such strong coupling regimes based on the semi-classical approaches.

To investigate the feature of the localization-delocalization transition in the spectral structure, we show the lowest 5 eigenvalues as a function of the coupling constant in the strong coupling regime when $`ϵ=0`$ in Fig. 4. When the electronic coupling is about the half of the reorganization energy as predicted by Eq. (74) $`\beta V5`$, the first nonzero eigenvalue (solid line) is no longer well separated from higher complex eigenvalues, it is no longer valid to describe the dynamics of the density matrix by a incoherent rate process. Also, there is a rapid drop of the first excited eigenvalue around $`\beta V6`$ due to the crossing of the first excited state eigenvalue and higher complex eigenvalues (note that $`\beta V_c=5`$ in this case). The crossing behavior of the first and higher excited state eigenvalues is a signature of the transition between the incoherent to coherent regimes in electron transfer processes.

A crossing of the two lowest nonzero eigenvalues can be used as an approximate criterion for the transition from the incoherent to coherent regimes. In Fig. (5), we plotted the value of the electronic coupling constant, $`V_c`$, where a crossing of the two lowest nonzero eigenvalues occur for a given energy bias, $`|ϵ|`$, and the prediction of Eq. (74) for the case $`\beta \mathrm{\Omega }=0.01`$ and $`\beta \lambda =10`$. They show the same trend that the coherent regime appears earlier in the asymmetric reaction case than in the symmetric reaction case.

### D Dynamic Bath Effect

In addition to the fact that the overall electron transfer rate appears a combined result of the nonadiabatic transition and solvent diffusion processes in Eq. (56), the nonadiabatic transition rate itself, Eq. (57), is modulated by the dynamic nature of the solvent relaxation process that is characterized by the bath correlation function in Eq. (58). We will call this a dynamic bath effect. This will be evident if the bath relaxation rate is comparable to the amount of the polarization energy fluctuation,

$`\mathrm{\Omega }\mathrm{\Delta }=\sqrt{2\lambda k_BT}.`$ (75)
To investigate the dynamic bath effect, we compare the eigenvalue solution and two different reaction rates, Marcus rate in Eq. (60) (static bath limit) and the nonadiabatic transition rate, Eq. (55) (dynamic bath effect included) in the weak coupling limit as we vary the energy bias in Fig. 6. All the calculation results show a characteristic behavior of the Marcus curve; they show a maximum value when $`|ϵ|\lambda `$. In the slow bath case, Fig. 6(a), all the three curves agree with each other since the dynamic bath effect is negligible in this case. In the fast bath case, Fig. 6(b), we notice that the eigenvalue and the nonadiabatic transition rate agree with each other, and show a deviation from the Marcus rate. In particular, in the fast bath case, the maximum position of the eigenvalue as well as of the nonadiabatic transition rate shift toward a smaller value of the energy bias compared with the Marcus rate.

In the static bath limit, the nonadiabatic transition rate is reduced to the Marcus rate, exhibiting the maximum position at the activationless case, $`|ϵ|=\lambda `$. When the bath relaxation is fast, the solvent reorganization energy is effectively reduced, and the maximum electron transfer rate appears at a smaller value of the bare solvent reorganization energy, exhibiting the peak position in the normal regime, $`|ϵ|<\lambda `$, and this results in a significant deviation of the nonadiabatic transition rate from the Marcus rate. The amount of the peak shift as a function of the solvent relaxation rate can also be estimated from Eq. (57), and it is shown in Fig. 7 that the peak shift increases linearly with $`\mathrm{\Omega }`$.

The peak shift observed in the eigenvalue solution has a close analogy to the Stokes shift observed in the condensed phase spectroscopy. It is worthwhile to mention that $`g(t)`$ obtained in Eq. (58) has exactly the same form as the line broadening function obtained in the hight temperature limit of the overdamped Brownian oscillator model in the condensed phase spectroscopy, and the nonadiabatic transition rates given in Eq. (57) have the same functional forms as the line shapes of the chromophore interacting with the Debye solvent. As discussed in the App. B, the Stokes shift between the absorption and the emission profiles disappears in the fast bath case and it has the same origin as the peak shift observed here.

Implication of the dynamic bath effect on the experimental observables may be significant. The reorganization energy of the solvent can be estimated from the maximum of the electron transfer rates based on a series of the measurements of electron transfer rates for different systems of which the energy bias and the electronic coupling constants are known. However, the dynamic bath effect can make the value of the reorganization energy measured in the experiment smaller than its true value. Therefore, care should be taken to disentangle the dynamical bath effect from the measured value of the solvent reorganization energy.

### E Density Matrix Propagation

The spectral method can be used as a propagation method of the density matrix as shown in Eq. (43). We demonstrate the time evolution of the donor and acceptor populations, $`P_1(t)`$ and $`P_2(t)`$ based on the spectral method. Parameters are chosen as $`\beta \lambda =10`$, $`\beta \mathrm{\Omega }=1`$, and $`\beta ϵ=0`$. The electronic coupling constant will be varied. The initial condition is chosen such that at $`t=0`$ only the donor state is populated with the equilibrium distribution in the diabatic state,

$`\rho _{11}(E,0)`$ $`=`$ $`\rho _{11}^{\mathrm{eq}}(E),`$ (76)
$`\rho _{22}(E,0)`$ $`=`$ $`\rho _{12}(E,0)=\rho _{21}(E,0)=0.`$ (77)
This initial condition corresponds to that of the photo-induced back electron transfer reaction.

In Fig. 8 we show the time evolution of the donor (solid line) and the acceptor (dashed line) state populations with various electronic coupling constants. At the smallest coupling constant, $`\beta V=0.5`$, the donor and the acceptor populations relax exponentially to the equilibrium at long times. The electronic coupling constant is small such that the time evolution is governed by the lowest nonzero eigenvalue of the nonadiabatic diffusion equation. As the electronic coupling constant is increased, the temporal behavior of the density matrix undergoes the incoherent-coherent transition as discussed based on the spectral analysis in Sec. IV D. When the coherent oscillation is observed in the temporal evolution of the density matrix($`\beta V=2`$ and $`\beta V=4`$), the oscillation frequency matches the Rabi oscillation frequency given in Eq. (73). The time-dependent study of the nonadiabatic diffusion equation based on the spectral method was also performed for back electron transfer dynamics occurring in the mixed-valence system.

### F Nonphysical Behavior in Nonadiabatic Diffusion Equation

We found that the nonadiabatic diffusion equation can yield nonphysical behavior in some parameter regimes. We checked the case of large reorganization energy, $`\beta \lambda 1`$, and found that in this case the smallest eigenvalue can be negative instead of zero, which yields an exponential growth of the density matrix in the long time limit. Fig. 9 demonstrates how two lowest real eigenvalues change as the reorganization energy is varied in the weak coupling case. When $`\beta \lambda 12`$, the lowest nonzero eigenvalue becomes negative instead of positive, thus violating the positive definiteness of the density matrix in the nonadiabatic diffusion equation.

There are several possible reasons for this problematic behavior. One of them is an intrinsic limitation in consistent treatments of the electronic and nuclear degrees of freedom in the mixed quantum-classical approach. Similar problem in the mixed quantum-classical approach has been pointed out in other contexts. The other possible reason for the positivity violation may be the Markovian approximation used in the calculation of the frictional kernel in the overdamped regime for the nonadiabatic dynamics. Although the Markovian approximation for the friction kernel has been successfully applied to the adiabatic reaction rate theory where the single potential energy surface is involved, there have been limited studies on the validity of this approximation in the nonadiabatic reaction cases. Along this line, similar nonphysical behaviors in the nonadiabatic diffusion equation yielding negative electron transfer rates have been found recently. The precise origin of this artifact is still not clear, and more studies need to be done to clarify this issue in a decisive way.

## V Conclusion

In this paper we studied the electron transfer reaction kinetics in Debye solvents described by the two-state nonadiabatic diffusion equation. By applying the spectral analysis method developed recently, we investigated various kinetic regimes of the electron transfer in Debye solvents in a unified way. The results obtained in this work are summarized;

(i) The general spectral structure of the nonadiabatic diffusion equation has been revealed. It exhibits a tree-like structure with three major branches; a single branch of real eigenvalues corresponding to multi-exponential decays and two symmetric branches of complex conjugate eigenvalues corresponding to damped oscillations.

(ii) Several kinetic regimes in the electron transfer reactions in solutions have been identified. The first nonzero eigenvalue solution bridges between the Marcus, the solvent-controlled, and the adiabatic crossing regime in the incoherent kinetics regime. It agrees well with nonadiabatic reaction rates from previous theoretical results in the weak coupling regime. In the adiabatic crossing regime the eigenvalue solution agrees with adiabatic reaction rates qualitatively.

(iii) When the electronic coupling constant is about half the reorganization energy in the symmetric reaction case (less than half in the asymmetric reaction case), the eigenvalue diagram shows coalescence/bifurcation behavior, where the two lowest nonzero eigenvalues cross each other and become complex conjugate. This indicates the damped oscillation behavior in the overdamped solvent, and correlates with the localization-delocalization transition in the lower adiabatic surface.

(iv) Due to the finite timescale of the solvation, the eigenvalue solution is modulated in the weak coupling regime such that the maximum position of the electron transfer rate appears at a smaller value of the bare solvent reorganization energy in the Marcus curve, which shows a substantial deviation from the standard Marcus theory. The consequence of the dynamic bath effect may be significant in the experimental determination of the solvent reorganization energy.

(v) The spectral analysis method was used as a propagation scheme to study the time-dependent behavior of the density matrix, and it has confirmed previous predictions made based on the investigation of the spectral structure.

(vi) Nonphysical behavior of the nonadiabatic diffusion equation, violating the positivity, has been identified in the large reorganization energy case by the spectral analysis method, and the precise origin of this problem should deserve further investigations.

In the isolated quantum system the eigenvalue solution method offers a powerful way to analyze the dynamics. In comparison, the spectral structure of the dissipative quantum dynamics has not been fully explored yet, and most studies on the quantum dissipative system employ trajectory methods. In this regard, our studies on the spectral structure of the electron transfer kinetics can be a good motivation to alternative studies of various dissipative systems. The nonadiabatic diffusion equation studied in this work is based on the mixed quantum-classical approach. However, the spectral analysis method itself can be applied to other quantum dissipative dynamical equations in general such as the Redfield equation, quantum Fokker-Planck equation, and quantum master equation, and the studies based upon the spectral analysis will bring more direct insights on the dissipative dynamics.

## Acknowledgments

YJ thanks Dr. Seogjoo Jang and Prof. Robert J. Silbey for useful discussions and Prof. Robert J. Silbey for his support. This work was supported by grants from Research Corporation and National Science Foundation.

## A Spectral analysis for Fokker-Planck operator on a single harmonic potential

For a single harmonic potential, $`U(x)=\frac{1}{2}m\omega ^2x^2`$, the Fokker-Planck operator is given by

$`_{\mathrm{FP}}=D\left({\displaystyle \frac{^2}{x^2}}+\beta {\displaystyle \frac{}{x}}U^{}\right),`$ (A1)
and it is possible to transfrom $`_{\mathrm{FP}}`$ into a quantum mechanical Hamiltonian corresponding to a fictitious harmonic oscillator in imaginary time,

$`H=e^{\beta U(x)/2}_{\mathrm{FP}}e^{\beta U(x)/2}={\displaystyle \frac{1}{2\mu }}{\displaystyle \frac{^2}{x^2}}+V(x),`$ (A2)
where $`\mu =(2D)^1`$ and the potential for the fictitious harmonic oscillator is given by

$`V(x)=D\left[{\displaystyle \frac{1}{4}}(\beta U^{^{}}(x))^2{\displaystyle \frac{1}{2}}\beta U^{^{\prime \prime }}(x)\right]={\displaystyle \frac{1}{2}}\mu \gamma ^2x^2{\displaystyle \frac{\gamma }{2}},`$ (A3)
with $`\gamma =Dm\omega ^2/k_BT`$. Since the transformed potential in Eq. (A3) is in the same form as that of a harmonic oscillator model with an off-set of the zero point energy, eigenvalues and eigenfunctions for the original Fokker-Planck operator can be obtained immediately from the eigenvalue solutions of the harmonic oscillator Hamiltonian,

$`_{\mathrm{FP}}|\psi _n^R=n\gamma |\psi _n^R,`$ (A4)
$`\psi _n^L|_{\mathrm{FP}}=n\gamma \psi _n^L|,`$ (A5)
where $`\psi _n^L(x)=H_n(x/\sqrt{2}x_0)/((2\pi )^{1/2}2^nn!x_0)^{1/2}`$ and $`\psi _n^R(x)=e^{x^2/2x_0^2}\psi _n^L(x)`$ with $`x_0=\sqrt{k_BT/m\omega ^2}`$ and $`H_n`$ being the $`n`$th order Hermite polynomial.

## B Rate Equations in Nonadiabatic Regime

We briefly review the rate expressions in the nonadiabatic limit ($`\beta V1`$).  First, we obtain a formal solution for $`\rho _{12}(E,t)`$ in the time domain from Eq. (7),

$`\rho _{12}(E,t)=iV{\displaystyle _0^t}d\tau {\displaystyle _{\mathrm{}}^{\mathrm{}}}dE_0G_{12}(E,\tau |E_0)\left[\rho _{11}(E_0,t\tau )\rho _{22}(E_0,t\tau )\right],`$ (B1)
where $`G_{12}(E,t|E_0)`$ is the Green’s function for the operator, $`_{12}i\omega _{12}`$,

$`G_{12}(E,t|E_0)=E|e^{(_{12}i\omega _{12})t}|E_0.`$ (B2)
Substituting $`\rho _{12}`$ and $`\rho _{21}=\rho _{12}^{}`$ from Eq. (B1) into the equations for $`\rho _{11}`$ and $`\rho _{22}`$ in Eqs. (5) and (6), we can obtain a closed set of integro-differential equations for $`\rho _{11}`$ and $`\rho _{22}`$,

$`{\displaystyle \frac{}{t}}\rho _{ii}(E,t)`$ $`=`$ $`_{ii}\rho _{ii}(E,t)`$ (B3)
$``$ $`2V^2\mathrm{Re}{\displaystyle _{\mathrm{}}^{\mathrm{}}}dE_0{\displaystyle _0^t}d\tau G_{12}(E,\tau |E_0)[\rho _{ii}(E_0,t\tau )\rho _{jj}(E_0,t\tau )].`$ (B4)
In the weak coupling regime, $`\beta V1`$, the electron transfer process is characterized as nonadiabatic, and the transition between two electronic states takes place only when the bath polarization energy is close to the free energy gap, and $`\rho _{11}`$ and $`\rho _{22}`$ vary slowly compared with $`\rho _{12}`$ and $`\rho _{21}`$. We can assume that the coherent Green’s function $`G_{12}`$ varies much faster than the population difference $`\rho _{11}\rho _{22}`$ in this case and that dominant contribution to the integral in Eq. (B4) comes from the short time region. Using these arguments we can approximately deconvolute the integrals and extend the upper limit of the integral to the infinity,

$`2V^2\mathrm{Re}{\displaystyle _{\mathrm{}}^{\mathrm{}}}dE_0{\displaystyle _0^t}d\tau G_{12}(E,\tau |E_0)[\rho _{11}(E_0,t\tau )\rho _{22}(E_0,t\tau )]`$ (B5)
$``$ $`K(E)[\rho _{11}(E,t)\rho _{22}(E,t)],`$ (B6)
where $`K(E)`$ is the rate kernel and is given in terms of the coherent Green function $`G(E,t|E_0)`$ defined in Eq. (B2),

$`K(E)=2V^2\mathrm{Re}{\displaystyle _0^{\mathrm{}}}d\tau {\displaystyle _{\mathrm{}}^{\mathrm{}}}dE_0G_{12}(E,\tau |E_0).`$ (B7)
The approximations made here are well established in the normal regime of the electron transfer, but they are questionable in the inverted or activationless regime, and improvements on these regimes have been made.

The nonadiabatic rates are given in terms of the coherent Green’s function,

$`k_{\mathrm{NA}}^i=2V^2\mathrm{Re}{\displaystyle _0^{\mathrm{}}}dt{\displaystyle _{\mathrm{}}^{\mathrm{}}}dE{\displaystyle _{\mathrm{}}^{\mathrm{}}}dE_0G_{ij}(E,t|E_0)\rho _{ii}^{\mathrm{eq}}(E_0).`$ (B8)
The coherent Green’s function, $`G_{12}(E,t|E_0)`$, can be evaluated first by transforming the Fokker-Planck operator into the harmonic oscillator Hamiltonian as done in Eq. (A2),

$`H_{12}(E)`$ $`=`$ $`e^{\beta \overline{U}(E)/2}(_{12}(E)i\omega _{12}(E))e^{\beta \overline{U}(E)/2}`$ (B9)
$`=`$ $`{\displaystyle \frac{1}{2m}}{\displaystyle \frac{^2}{\stackrel{~}{E}^2}}+{\displaystyle \frac{1}{2}}m\mathrm{\Omega }^2\stackrel{~}{E}^2i\stackrel{~}{ϵ},`$ (B10)
where $`m=1/(2D_E)=1/(2\mathrm{\Omega }\mathrm{\Delta }^2)`$. Two complex energy variables are defined by $`\stackrel{~}{E}=E+i\left(\frac{2\mathrm{\Delta }^2}{\mathrm{\Omega }}\right)`$ and $`\stackrel{~}{ϵ}=ϵ+i\left(\frac{\mathrm{\Delta }^2}{\mathrm{\Omega }}\frac{\mathrm{\Omega }}{2}\right)`$. Since $`H_{12}`$ is a simple harmonic oscillator Hamiltonian with respect to the complex energy variable $`\stackrel{~}{E}`$, we can calculate its propagator, $`\stackrel{~}{E}|e^{H_{12}t}|\stackrel{~}{E}_0`$. Then, we can express $`G_{12}(E,t|E_0)`$ in Eq. (B2) explicitly,

$`G_{12}(E,t|E_0)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{4\pi \mathrm{\Delta }^2\mathrm{sinh}(\mathrm{\Omega }t)}}}\mathrm{exp}[{\displaystyle \frac{1}{4\mathrm{\Delta }^2\mathrm{sinh}(\mathrm{\Omega }t)}}(e^{\mathrm{\Omega }t}E^2e^{\mathrm{\Omega }t}E_0^22EE_0)`$ (B12)
$`+({\displaystyle \frac{2\mathrm{\Delta }^2}{\mathrm{\Omega }^2}}i{\displaystyle \frac{(E+E_0)}{\mathrm{\Omega }}})\mathrm{tanh}\left({\displaystyle \frac{\mathrm{\Omega }t}{2}}\right)+i\stackrel{~}{ϵ}t].`$
Straightforward Gaussian integrations over the space coordinates in Eq. (B8) result in Eq. (57).

Depending on the timescale of the bath, $`\mathrm{\Omega }`$, and the amount of the polarization energy fluctuation, $`\mathrm{\Delta }=\sqrt{2\lambda k_BT}`$, the nonadiabatic reaction rates in Eq. (57) have two different limits.  In the static bath limit $`(\mathrm{\Omega }\mathrm{\Delta })`$, we take the short time limit $`\mathrm{\Omega }t1`$ of $`g(t)`$ in Eq. (58), $`g(t)\frac{1}{2}\mathrm{\Delta }^2t^2`$, then the nonadiabatic transition rate yields the Marcus rate given in Eq. (59) after a Gaussian integration. The maxima of $`k_{\mathrm{NA}}^1`$ and $`k_{\mathrm{NA}}^2`$ are separated by $`2\lambda `$ in this case and this is analogous to the Stokes shift observed in condensed phase spectroscopy. When the dynamics of bath degrees of freedom is very fast $`(\mathrm{\Omega }\mathrm{\Delta })`$, we can take the long time limit $`\mathrm{\Omega }t1`$ of $`g(t)`$, $`g(t)(\mathrm{\Delta }^2/\mathrm{\Omega }i\lambda )t`$, which leads to,

$`k_{\mathrm{NA}}^1k_{\mathrm{NA}}^2{\displaystyle \frac{2V^2\gamma }{ϵ^2+\gamma ^2}},`$ (B13)
where $`\gamma =\mathrm{\Delta }^2/\mathrm{\Omega }`$. The Stokes shift between $`k_{\mathrm{NA}}^1`$ and $`k_{\mathrm{NA}}^2`$ has now disappeared as a result of fast fluctuation of the solvent polarization energy.

The solvent diffusion rate is given in terms of the population Green’s function $`G_{ii}(E,t|E_0)`$ at the crossing point $`E=E_0=ϵ`$ in Eq. (61). Since $`_{ii}`$ is given by the Fokker-Planck operator with a displaced harmonic potential for each electronic state, the population Green’s functions $`G_{ii}`$ are easily obtained

$`G_{ii}(E,t|E_0)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{\Delta }^2(1e^{2\mathrm{\Omega }t})}}}\mathrm{exp}\left[{\displaystyle \frac{(E\pm \lambda (E_0\pm \lambda )e^{\mathrm{\Omega }t})^2}{2\mathrm{\Delta }^2(1e^{2\mathrm{\Omega }t})}}\right],`$ (B14)
with $`+()`$ sign for $`i=1(2)`$, and the solvent diffusion rates can be calculated.

In the cases of high and low barrier limits, the solvent diffusion rate given in an integral form in Eq. (61) can be approximately evaluated,

$`k_\mathrm{D}^i\{\begin{array}{cc}\frac{\mathrm{\Omega }}{\mathrm{ln}2}\left[1\left(\frac{2}{\mathrm{ln}2}\right)\beta E_\mathrm{b}^i+\left(\left(\frac{2}{\mathrm{ln}2}\right)^2\frac{2}{3\mathrm{ln}2}\right)(\beta E_\mathrm{b}^i)^2+\mathrm{}\right]\hfill & \beta E_\mathrm{b}^i1\hfill \\ \mathrm{\Omega }\sqrt{\frac{\beta E_\mathrm{b}^i}{\pi }}\mathrm{exp}(\beta E_\mathrm{b}^i)\hfill & \beta E_\mathrm{b}^i1\hfill \end{array}`$ (B17)
where $`E_\mathrm{b}^i`$ is the forward ($`i=1`$) or backward ($`i=2`$) reduced barrier height, $`E_\mathrm{b}^i=(ϵ\pm \lambda )^2/(4\lambda )`$.

The Zusman’s solution given in Ref.  amounts to two approximations made in the nonadiabatic solution. First, it is assumed that electron transfer occurs only when the solvation polarization energy $`E`$ matches with the energy bias $`ϵ`$, which yields,

$`K(E)2\pi V^2\delta (Eϵ).`$ (B18)
In this case, Eq. (45) is further reduced to the localized transition model, and the nonadiabatic transition rates $`k_{\mathrm{NA}}^i`$ reduces to the original Marcus expressions given in Eq. (59). The second approximation made by Zusman is that the solvent diffusion rates are approximated by their high-barrier limits given in Eq. (B17) for $`\beta E_\mathrm{b}^i1`$, which can be written as

$`k_\mathrm{D}^i\mathrm{\Omega }|\lambda \pm ϵ|\rho _{ii}^{\mathrm{eq}}(ϵ).`$ (B19)