text:
# Electronic coherence in mixed-valence systems: Spectral analysis

## I INTRODUCTION

Quantum coherence in the dynamics of condensed phase systems has become a subject of recent experimental and theoretical studies. A central issue is the observability of electronic coherence in electron transfer systems given the fast dephasing time in many-body quantum systems. Experimentally, with the advance in ultrafast laser technology, oscillations in electronic dynamics have been observed in photo-synthetic reaction centers and other electron transfer systems and are believed to arise from vibrational and/or electronic coherence. Accurate measurements on photo-induced electron transfer in mixed-valence compounds have demonstrated oscillations in electronic populations on the femtosecond time-scale. Theoretically, detailed path-integral simulations suggest that such oscillations take place in electron transfer systems with large electronic coupling constants and are sensitive to the initial preparation of the bath modes associated with the transfer processes. Lucke et al. extended the non-interacting blip approximation to incorporate the non-equilibrium initial preparation and carried out extensive path-integral quantum dynamics simulations for electron transfer reactions. According to their findings, large-amplitude oscillations are most likely to be observed in symmetric mixed-valence systems that are nearly adiabatic and with initial configurations that are centered in the Landau-Zener crossing region. Using the transfer matrix technique, Evans, Nitzan, and Ratner calculated short-time evolution for the photo-induced electron transfer reaction in $`(\mathrm{NH}_3)_5\mathrm{Fe}^{\mathrm{II}}(\mathrm{CN})\mathrm{Ru}^{\mathrm{III}}(\mathrm{CN})_5`$. Their results show fast oscillations in the electronic population on the short time-scale(20 fs) followed by a slower population relaxation on the long time-scale(100 fs). They pointed out that these fast oscillations arise as the wave-function oscillates coherently between the donor and acceptor states. The calculated long-time decay rate is considerably smaller than the prediction by the golden-rule formulae, confirming the inadequacy of non-adiabatic rate theory in studying mixed-valence systems.

In fact, a simple classical argument helps understand the nature of the observed oscillations. As a function of the ratio between $`\lambda `$ (the bath reorganization energy) and $`V`$ (the electronic coupling constant), there is a thermodynamic transition from the localized electronic state in a double-well potential to the delocalized electronic state in a single well potential. (i) In the localized regime ($`\lambda V`$), the large reorganization energy destroys electronic coherence; hence, electron transfer is an incoherent rate process, which can be described by the non-interacting blip approximation or golden-rule rate in the non-adiabatic limit and by transition state theory in the adiabatic limit. (ii) In the delocalized regime ($`\lambda V`$), the electronic wave function extends to both the donor and acceptor states and electronic coherence persists over several oscillations. For mixed-valence compounds, the electronic coupling constant is estimated to be in the range of $`10^3\text{cm}^1`$, which is in the same order as the reorganization energy. Therefore, the observed oscillations and relaxation in mixed-valence systems are the consequence of a highly non-equilibrium coherence transfer process.

Due to the delocalization nature of electronic states, an adiabatic picture is more useful than the diabatic representation for analyzing the short-time dynamics in strongly-coupled systems. In this picture, electronic coherence arises from Rabi oscillations between two adiabatic surfaces and decays because of electronic dephasing. Further, initial preparation and wave-packet dynamics can modulate Rabi oscillations and the overall electronic dynamics. Thus, the adiabatic representation provides a simple picture for mixed-valence systems as well as a simple analytical method to model fast electron dynamics initiated by laser pulses.

As a general approach to describe condensed phase dynamics, we recently proposed a spectral analysis method, which is based on eigen-structures of dissipative systems instead of dynamic trajectories. An important application of the approach is to analyze a set of two-state diffusion equations, which was first used by Zusman to treat solvent effects on electron transfer in the non-adiabatic limit. The analysis allows us to characterize multiple time-scales in electron transfer processes including vibrational relaxation, electronic coherence, activated curve crossing or barrier crossing. With this unified approach, the observed rate behavior, bi-exponential and multi-exponential decay, and population oscillations are different components of the same kinetic spectrum. Thus, several existing theoretical models, developed for limited cases of electron transfer, can be analyzed, tested, and extended. In particular, rate constants extracted from the analysis bridge smoothly between the adiabatic and non-adiabatic limits, and the kinetic spectrum in the large coupling regime reveals the nature of the localization-delocalization transition as the consequence of two competing mechanisms.

In this paper, the spectral analysis approach developed in Ref. is employed to study the electron transfer dynamics in mixed-valence systems. We invoke the non-adiabatic diffusion equation proposed by Zusman to describe the electron transfer process in the over-damped solvent regime. As discussed earlier, electron transfer in mixed-valence systems takes place in a different kinetic regime from the thermal activated regime described by Marcus theory. Thus, the time-scale separation is not satisfied, and multi-exponential decay and oscillations are intrinsic nature of electron transfer kinetics. As a result, the kinetic spectra exhibit bifurcation, coalescence, and other complicated patterns. Careful examination of these patterns reveals the underlying mechanisms in mixed-valence systems.

The rest of the paper is organized as follows: The spectral structure of the non-adiabatic diffusion equation is formulated in Sec. II. Numerical examples of the spectral structure of strongly mixed electron transfer systems are presented and discussed in Sec. III and concluding remarks are given in Sec. IV.

## II THEORY

There have been extensive studies of the solvent effect on electron transfer dynamics in literature with various approaches. One of the most extensively studied models for quantum dissipation is the spin-boson Hamiltonian,

$`H_{SB}={\displaystyle \frac{ϵ}{2}}\sigma _z+V\sigma _x+{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{p_\alpha ^2}{2m_\alpha }}+{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{1}{2}}m_\alpha \omega _\alpha ^2\left(x_\alpha \sigma _z{\displaystyle \frac{c_\alpha }{m_\alpha \omega _\alpha ^2}}\right)^2,`$ (1)
where $`ϵ`$ is the energy bias between the two electronic states, $`V`$ is the electronic coupling constant, $`\sigma _z`$ and $`\sigma _x`$ are the usual Pauli matrices, and $`\{x_\alpha `$,$`p_\alpha \}`$ represents the bath degree of freedom with mass $`m_\alpha `$, frequency $`\omega _\alpha `$, and the coupling constant $`c_\alpha `$. In this model effects of the bath modes on the dynamics of the system can be described via the spectral density defined by,

$`J(\omega )={\displaystyle \frac{\pi }{2}}{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{c_\alpha ^2}{m_\alpha \omega _\alpha }}\delta (\omega \omega _\alpha ).`$ (2)
Equivalently, the spin-boson Hamiltonian in Eq. (1) can be separated into the electronic two-level part $`H_{TLS}`$ and the nuclear bath part $`H_B`$,

$`H_{SB}=H_{TLS}+H_B.`$ (3)
The two-level part of the Hamiltonian can be explicitly written as

$`H_{TLS}(E)=U_1(E)|11|+U_2(E)|22|+V(|12|+|21|),`$ (4)
where the diabatic energy surfaces $`U_1(E)`$ and $`U_2(E)`$ are functions of the stochastic variable $`E`$, which represents the polarization energy for a given solvent configuration. The transformation from the spin-boson Hamiltonian to the two-level system Hamiltonian has been shown in the literature by the identity,

$`E(\{x_\alpha \})={\displaystyle \underset{\alpha }{}}c_\alpha x_\alpha .`$ (5)
It is worthwhile to mention that the polarization energy $`E`$ was recognized as the reaction coordinate by Marcus in formulating non-adiabatic electron transfer theory. Since the electron transfer process involves the collective motion of a large number of solvent degrees of freedom and the two-level system is linearly coupled to the harmonic bath modes in the spin-boson Hamiltonian in Eq. (1), the functional form for the free energy surface is harmonic, thus giving

$`U_1(E)`$ $`=`$ $`{\displaystyle \frac{(E+\lambda )^2}{4\lambda }},`$ (6)
$`U_2(E)`$ $`=`$ $`{\displaystyle \frac{(E\lambda )^2}{4\lambda }}+ϵ,`$ (7)
where $`\lambda `$ is the reorganization energy, which is related to the parameters in Eq. (1),

$`\lambda ={\displaystyle \underset{\alpha }{}}{\displaystyle \frac{c_\alpha ^2}{2m_\alpha \omega _\alpha ^2}}={\displaystyle \frac{1}{\pi }}{\displaystyle 𝑑\omega \frac{J(\omega )}{\omega }}.`$ (8)
Considering the fact that electron transfer processes are usually probed at room temperature in polar solvents, we can treat the bath degrees of freedom in $`H_B`$ classically. Then, the spin-boson Hamiltonian in Eq. (3) can be used to derive a two-level classical equation of motion,

$`i{\displaystyle \frac{}{t}}\rho (t)=\rho (t)=(_B+_{TLS})\rho (t),`$ (9)
where $`i_B=\{H_B,\}`$ is the Poisson operator for the classical bath and $`_{TLS}=[H_{TLS},]/\mathrm{}`$ is the Liouville operator for the two level system. Explicitly, we express Eq. (9) in terms of the density matrix elements,

$`\dot{\rho }_1`$ $`=`$ $`_1\rho _1+iV(\rho _{12}\rho _{21}),`$ (11)
$`\dot{\rho }_2`$ $`=`$ $`_2\rho _2iV(\rho _{12}\rho _{21}),`$ (12)
$`\dot{\rho }_{12}`$ $`=`$ $`_{12}\rho _{12}i\omega _{12}\rho _{12}+iV(\rho _1\rho _2),`$ (13)
$`\dot{\rho }_{21}`$ $`=`$ $`_{21}\rho _{21}+i\omega _{12}\rho _{21}iV(\rho _1\rho _2),`$ (14)
where the Planck constant $`\mathrm{}`$ is set to unity for simplicity, $`\rho _i`$ is the diagonal matrix element for electronic population, and $`\rho _{ij}`$ is the off-diagonal matrix element for electronic coherence. Here, $``$ describes the relaxation process of classical bath, with $`_i`$ defined on the free energy surface for the $`i`$th electronic state, and with $`_{12}`$ and $`_{21}`$ defined on the averaged free energy surface. This set of semi-classical two-state equations has been previously derived in different context by several authors.  It should be mentioned that the mapping from the spin-boson Hamiltonian into the Zusman model requires the Lorentzian form of the spectral density,

$`J(\omega )=2\lambda {\displaystyle \frac{\omega \omega _c}{\omega ^2+\omega _c^2}}.`$ (15)
Furthermore, we note that many chemically and biologically important electron transfer processes take place in the over-damped solvent environment. Therefore, to describe the density matrix evolution in the electron transfer kinetics in the mixed-valence system, we invoke the non-adiabatic diffusion equation proposed by Zusman. Then, the bath relaxation operators in Eq. (9) are one-dimensional Fokker-Planck operators $`_{ij}`$,

$`_i=D_E{\displaystyle \frac{}{E}}\left({\displaystyle \frac{}{E}}+\beta {\displaystyle \frac{U_i(E)}{E}}\right),`$ (16)
$`_{12}=_{21}={\displaystyle \frac{_{11}+_{22}}{2}}=D_E{\displaystyle \frac{}{E}}\left({\displaystyle \frac{}{E}}+\beta {\displaystyle \frac{\overline{U}(E)}{E}}\right).`$ (17)
where $`\beta =1/k_BT`$, $`\overline{U}`$ and $`\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}},`$ (18)
$`\omega _{12}(E)=U_1(E)U_2(E).`$ (19)
The energy diffusion constant $`D_E`$ is defined as

$`D_E=\mathrm{\Omega }\mathrm{\Delta }_E^2,`$ (20)
where $`\mathrm{\Delta }_E^2`$ is the mean square fluctuation of the solvent polarization energy

$$\mathrm{\Delta }_E^2=E^2=2\lambda k_BT,$$

and $`\tau _D=1/\mathrm{\Omega }`$ is the the characteristic timescale of the Debye solvent. The correlation function of the solvent polarization energy is given by

$`C(t)=E(t)E(0)=\mathrm{\Delta }_E^2\mathrm{exp}(\mathrm{\Omega }t).`$ (21)
Note that since the nuclear dynamics is modeled by the Fokker-Planck operator, the possibility of the vibrational coherence is excluded in this model of electron transfer dynamics. It is worthwhile to mention that one can obtain the non-adiabatic diffusion equation starting from the spin-boson Hamiltonian, by first deriving the evolution equation for the quantum dissipative dynamics, and then taking the semi-classical limit using the Wigner distribution functions, and finally assuming the over-damped diffusion limit.

We investigate the spectral structure of the non-adiabatic diffusion operator by calculating the eigenvalues $`\{Z_\nu \}`$ and the corresponding eigen-functions $`\{|\psi _\nu \}`$. Hereafter we use Greek indices to denote the eigenstates and Latin indices to denote the basis states of the non-adiabatic diffusion operator. Because the non-adiabatic Liouville operator is non-Hermitian, the eigenvalues are generally given by complex values, and the right and left eigen-functions corresponding to the same eigenvalue are not simply the Hermitian conjugate to each other. For a given eigen-value $`Z_\nu `$, the right and left eigen-functions of the non-adiabatic diffusion operator are obtained from

$`|\psi _\nu ^R`$ $`=`$ $`Z_\nu |\psi _\nu ^R,`$ (22)
$`\psi _\nu ^L|`$ $`=`$ $`Z_\nu \psi _\nu ^L|.`$ (23)
The method of eigenfunction solution is well known for the diffusion process on the harmonic potential energy surface. For a single quadratic potential $`U(x)=\frac{1}{2}m\omega ^2x^2`$, the one-dimensional Fokker-Planck operator $`_{FP}=D(\frac{^2}{x^2}+\beta \frac{}{x}U^{})`$ can be transformed into the quantum mechanical Hamiltonian in imaginary time,

$`H_s=e^{\beta U(x)/2}_{FP}e^{\beta U(x)/2}={\displaystyle \frac{1}{2\mu }}{\displaystyle \frac{^2}{x^2}}+V_s(x),`$ (24)
where $`\mu ^1=2D`$, and the quadratic potential is

$`V_s(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}},`$ (25)
with $`\gamma =Dm\omega ^2/k_BT`$. Since the transformed potential in Eq. (25) is just the same form as for a simple harmonic oscillator with zero point energy compensation, the eigenvalues and the eigen-functions for the original Fokker-Planck operator can be constructed immediately from the eigen-solutions of the harmonic oscillator Hamiltonian. Unlike the diffusion problem on the single potential energy surface, there have been limited studies on the non-adiabatic diffusion problem involving more than one potential energy surface. In this aspect, Cukier and co-workers have calculated the electron transfer rate by calculating the lowest eigenvalue of the non-adiabatic diffusion equation; however, their calculation was limited to the weak-coupling regime where the Zusman rate is applicable.

An important issue in solving the non-adiabatic diffusion equation for electron transfer is the choice of the basis functions since three different free energy surfaces are involved in Eq. (9): two diabatic surfaces for the population density matrix elements and one averaged surface for the coherence density matrix element. In this paper, the eigen-functions of $`_{12}`$ are used as our basis set to represent the non-adiabatic diffusion equation. In principle, one could have chosen the eigen-functions of $`_1`$ or $`_2`$ as basis functions, however, in that case one has to evaluate appropriate Franck-Condon factors when calculating the coupling matrix elements even with the Condon approximation. The Fokker-Planck operator $`_{12}`$ is defined on the averaged harmonic potential centered at $`E=0`$, and its eigen-solutions are

$`_{12}|\varphi _n^R`$ $`=`$ $`n\mathrm{\Omega }|\varphi _n^R,`$ (26)
$`\varphi _n^L|_{12}`$ $`=`$ $`n\mathrm{\Omega }\varphi _n^L|,`$ (27)
where the right and left eigen functions are

$`\varphi _n^R(E)={\displaystyle \frac{1}{(2^nn!)^{\frac{1}{2}}(2\pi \mathrm{\Delta }_E^2)^{\frac{1}{4}}}}\mathrm{exp}\left({\displaystyle \frac{E^2}{2\mathrm{\Delta }_E^2}}\right)H_n\left({\displaystyle \frac{E}{\sqrt{2}\mathrm{\Delta }_E}}\right),`$ (28)
and

$`\varphi _n^L(E)={\displaystyle \frac{1}{(2^nn!)^{\frac{1}{2}}(2\pi \mathrm{\Delta }_E^2)^{\frac{1}{4}}}}H_n\left({\displaystyle \frac{E}{\sqrt{2}\mathrm{\Delta }_E}}\right),`$ (29)
where $`H_n`$ is the $`n`$th order Hermite polynomial. As shown below, this choice of the basis set is convenient for our purpose.

To be consistent with the $`_{12}`$ basis set, we separate the real and imaginary parts of the coherence density matrix, namely, $`u=`$Re$`\rho _{12}`$ and $`v=`$Im$`\rho _{12}`$, and rewrite Eq. (9) as

$`\dot{\rho }_1`$ $`=`$ $`(_{12}+\delta )\rho _12Vv,`$ (31)
$`\dot{\rho }_2`$ $`=`$ $`(_{12}\delta )\rho _2+2Vv,`$ (32)
$`\dot{u}`$ $`=`$ $`_{12}u+\omega _{12}v,`$ (33)
$`\dot{v}`$ $`=`$ $`_{12}v\omega _{12}u+V(\rho _1\rho _2),`$ (34)
where we have defined $`\delta `$ as

$`\delta ={\displaystyle \frac{_{11}_{22}}{2}}.`$ (35)
Then, all the relevant operators in Eqs. (31)-(34) can be evaluated in terms of the right and left eigen-functions of $`_{12}`$, giving

$`\varphi _n^L|_{12}|\varphi _m^R`$ $`=`$ $`n\mathrm{\Omega }\delta _{nm},`$ (36)
$`\varphi _n^L|\delta |\varphi _m^R`$ $`=`$ $`\mathrm{\Omega }\sqrt{{\displaystyle \frac{\lambda }{2k_BT}}}\sqrt{m+1}\delta _{n,m+1},`$ (37)
$`\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},`$ (38)
$`\varphi _n^L|V|\varphi _m^R`$ $`=`$ $`V\delta _{nm},`$ (39)
where we assume the Condon approximation, i.e., the electronic coupling matrix element is independent of the solvent degrees of freedom. With the basis set, we can expand the density matrix elements as

$`\rho _1(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}a_n(t)\varphi _n^R(E),`$ (41)
$`\rho _2(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}b_n(t)\varphi _n^R(E),`$ (42)
$`u(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}c_n(t)\varphi _n^R(E),`$ (43)
$`v(E,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}d_n(t)\varphi _n^R(E).`$ (44)
Substituting Eqs. (41)-(44) into the eigenvalue equation Eq. (22), we have the following coupled linear equations

$`Z_\nu a_n`$ $`=`$ $`n\mathrm{\Omega }a_n\mathrm{\Omega }\sqrt{{\displaystyle \frac{\lambda }{2k_BT}}}\sqrt{n}a_{n1}2Vd_n,`$ (46)
$`Z_\nu b_n`$ $`=`$ $`n\mathrm{\Omega }b_n+\mathrm{\Omega }\sqrt{{\displaystyle \frac{\lambda }{2k_BT}}}\sqrt{n}b_{n1}+2Vd_n,`$ (47)
$`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,`$ (48)
$`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),`$ (49)
which is an explicit basis set representation for the two-state diffusion operator in Eq. (9). The linear equations for the left eigen-solution as defined by Eq. (23) can be written by the transpose of Eqs. (46)-(49). Diagonalizing the $`4N\times 4N`$ matrix ($`N`$ = number of basis functions) defined in Eqs. (46)-(49), we obtain the eigenvalues $`Z_\nu `$ and the corresponding eigenvectors of the non-adiabatic diffusion operator,

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

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

$`{\displaystyle \underset{n=0}{}}L_{\nu n}R_{n\nu ^{^{}}}=\delta _{\nu \nu ^{^{}}},`$ (52)
for the orthogonality and

$`{\displaystyle \underset{\nu }{}}R_{n\nu }L_{\nu m}=\delta _{nm},`$ (53)
for the completeness. Using these properties, we can construct the real time propagator for the operator $``$ as

$`G(t)={\displaystyle \underset{\nu }{}}|\psi _\nu ^R\psi _\nu ^L|e^{Z_\nu t},`$ (54)
and express the time evolution of the density matrix by projecting a given initial distribution onto the eigenstates, giving

$`|\rho (t)=G(t)|\rho (0)={\displaystyle \underset{\nu }{}}|\psi _\nu ^R\psi _\nu ^L|\rho (0)e^{Z_\nu t}.`$ (55)
Hence, the eigen-solution to the two-state non-adiabatic diffusion equation leads to a complete description of electron transfer dynamics.

## III RESULTS AND DISCUSSIONS

In the section, we present the spectral structure of the non-adiabatic diffusion operator by diagonalizing its matrix representation in Eqs. (46)-(49). In principle, we need infinite number of basis functions to diagonalize the non-adiabatic diffusion operator, however, in practice, we have to truncate our basis set at some finite number. In all the calculations below, we have used $`N=50200`$ to diagonalize the $`4N\times 4N`$ matrix and the effect of finite number basis on the spectral structure has been carefully examined.

### A Spectral Structure

#### 1 Mixed-valence systems

In the mixed-valence compounds, the electronic coupling constant has the same order of magnitude as the reorganization energy and the electron transfer dynamics is usually probed experimentally at room temperature in polar solvents. To study this process, Evans, Nitzan, and Ratner carried out real time path-integral simulations for the photo-induced electron transfer reaction in $`(\mathrm{NH}_3)_5\mathrm{Fe}^{\mathrm{II}}(\mathrm{CN})\mathrm{Ru}^{\mathrm{III}}(\mathrm{CN})_5`$. Based on their model, we chose the parameters for the calculation shown in Fig. 1 as $`\beta \mathrm{\Omega }=0.6716`$, $`\beta \lambda =18.225`$, $`\beta V=11.99`$, and $`\beta ϵ=18.705`$. As mentioned in the introduction the mapping between the spin-boson Hamiltonian and the semi-classical Zusman equation is not rigorously defined. For example, for the non-adiabatic diffusion equation, the solvation energy correlation function takes an exponential form with the rate $`\mathrm{\Omega }`$, whereas, for the spin-boson model Hamiltonian, it depends on the functional form of the spectral density. It can be shown that the Ohmic spectral density with an exponential cut-off $`\omega _c`$

$`J(\omega )=\eta \omega \mathrm{exp}(\omega /\omega _c),`$ (56)
used in the calculation of Evans et al., leads to an energy correlation function with a Lorentzian form at high temperature,

$`C_{SB}(t){\displaystyle \frac{2\eta \omega _ck_BT}{\pi }}{\displaystyle \frac{1}{1+(\omega _ct)^2}}.`$ (57)
Then, the relaxation rate $`\mathrm{\Omega }`$ used in our calculation is taken as the inverse of the mean survival time of $`C_{SB}(t)`$, which is $`\mathrm{\Omega }=2\omega _c/\pi `$.

In Fig. 1 the spectral structure of the non-adiabatic operator is shown in complex space. We have used $`N=200(4N=800)`$ basis functions to calculate the eigenvalues. To remove the effect of finite basis set from the resulting spectral structure, we only show the first 400 eigenvalues in the complex plane. Since the non-adiabatic diffusion operator is non-Hermitian, the resulting spectrum shows complex conjugate paired eigenvalues as well as real eigenvalues, giving rise to the tree structure with three major branches (which we will call the eigen-tree). In Fig .1, we separate the real and imaginary parts of eigenvalue by

$`Z_\nu =k_\nu i\omega _\nu .`$ (58)
Obviously, the real part, $`k_\nu `$, is always negative as all non-equilibrium physical quantities decay to zero at time infinity, and it scales linearly with the index $`\nu `$ since the relaxation rate corresponding to the $`n`$th basis state $`\varphi _n`$ is proportional to $`n`$. In general, the relative magnitudes of real and imaginary parts of eigenvalues determine the time-evolution of the density matrix: the real eigenvalues correspond to the simple exponential decay components and the complex conjugate paired eigenvalues correspond to the damped oscillation components.

To classify the eigenvalues quantitatively according to their dynamic behavior, we introduce the dimensionless quantity $`\theta _\nu `$

$`\theta _\nu {\displaystyle \frac{2\pi k_\nu }{|\omega _\nu |}},`$ (59)
where $`k_\nu `$ is the decay rate and $`2\pi /\omega _\nu `$ is the oscillation period. The time-evolution of the density matrix component associated with the eigenvalue $`Z_\nu `$ is an exponential decay if $`\theta _\nu =\mathrm{}`$, an under-damped oscillation if $`\theta _\nu >1`$, and a damped oscillation if $`\theta _\nu 1`$. The relative amplitude of the each component depends on the overlap matrix element between the initial density matrix and the eigenstate. As an approximate criterion for the classification of the eigenvalues, the slope corresponding to $`\theta _\nu =1`$ is shown in the eigen-tree diagram in Fig. 1. There are a few eigenstates around and below the $`\theta _\nu =1`$ line, with a typical rate of $`\beta k_\nu 5`$. For the parameters used in the calculation, $`\beta `$ corresponds to $``$170 fs in real time, and, therefore, these eigenstates exhibits damped oscillations with a period and a decay time in the femtosecond regime. In their real-time path integral simulations, Evans et al. showed that the population in the acceptor state oscillates with a few femtosecond period and these oscillation decays in within 20 femtoseconds. Thus, qualitative features of the electron transfer dynamics can be predicted and understood from a careful examination of the spectral structure. Since the spectral analysis presented here is based on the semi-classical diffusion equation while the path-integral study is based on the quantum mechanical spin-boson Hamiltonian, the comparison between the two approaches is expected to be qualitative. In the following subsection, further analysis reveals the nature of these oscillations.

#### 2 Dependence on the coupling constant $`V`$

To examine the underlying spectral structure in more details, eigenvalues of the non-adiabatic diffusion operator are plotted as functions of the electronic coupling constant in Fig. 2. All the parameters except for the electronic coupling constant are the same as used in Fig. 1.

In Fig. 2(a), the real parts of the first 20 eigenvalues are shown as functions of the electronic coupling constant. Note that eigenvalues corresponding to complex conjugate pairs have the same real part, thus they coalesce in the real eigenvalue diagram. When the coupling constant is very small $`(\beta V1)`$, the real part of the first non-zero eigenvalue is very well separated from the eigenvalues of excited states, so the dynamics of electron transfer can be considered as a incoherent rate process with a well-defined rate constant, $`k_1`$. When the coupling constant is larger $`(\beta V1)`$, the first excited state becomes close to the second excited state, and they start to merge into a complex conjugate pair. If the coupling constant increases further, eigen-values show a bifurcation behavior at $`\beta V10`$. Therefore, in this regime, the electron transfer kinetics show multiple time-scale relaxation as well as coherent oscillation. The complicated behavior of coalescence and bifurcation in the real eigenvalue appears more frequently at higher states.

Another interesting feature of the real eigenvalue diagram is that a set of real eigenvalues decreases consistently as the coupling constant increases from zero. It turns out that these eigenstates take on large imaginary parts, which are responsible for the onset of the imaginary branches of the eigen-tree. In Fig. 2(b), the imaginary parts of the lowest 30 eigenvalues are plotted as functions of the coupling constant. Interestingly, the imaginary part of the eigenvalue increases approximately linearly with the coupling constant at large coupling regime. In fact, the dependence on the coupling constant is similar to that of the Rabi frequency for the two-level system,

$`\mathrm{\Omega }_R=\sqrt{ϵ^2+4V^2},`$ (60)
which is shown in Fig. 2(b). As pointed out in a recent paper, electronic coherence in mixed-valence systems arises from Rabi oscillations between two adiabatic surfaces and decays because of dephasing.

To demonstrate the correlation of the real and imaginary parts of the eigenvalues as functions of the coupling constant, we present a three dimensional plot of the spectral structure in Fig. 2(c). For clarity, only the positive branches of the imaginary eigenvalues are shown. If we compare Fig. 2(c) with Fig. 2(a), the very rapidly decaying states shown in Fig. 2(a) take on large imaginary parts corresponding to the Rabi oscillations as the coupling constant increases, and these states are responsible for the onset of the imaginary branches in the eigen-tree for the mixed-valence system shown in Fig. 1.

### B Density Matrix Propagation

To check the validity of the spectral analysis as a density matrix propagation scheme, we calculated the time-evolution of the density matrix by applying the propagator defined by Eq. (54) to the initial density matrix for various energy biases. Although it may seem straightforward to use the spectral method as a propagation scheme, the case for a non-Hermitian operator is not trivial and has not been explored. The main reason is that though the left and right eigen-functions of a non-Hermitian operator can be shown to form a bi-orthogonal set for the non-degenerate eigenvalue case, numerically these eigen-functions may not be stable enough to be used as a complete orthonormal basis for the density matrix propagation, especially in the nearly degenerate eigenvalue case. We can understand the situation as follows: When the two nearly degenerate eigenvalues $`Z_1`$ and $`Z_2`$ are obtained from a non-Hermitian operator, the orthogonality implies that $`L_2|`$ and $`|R_1`$ are orthogonal to each other as well as $`L_1|`$ and $`|R_2`$. When two eigenvalues become very close to each other, unlike the Hermitian operator case, $`L_1|`$ and $`L_2|`$ almost coincide and so do $`|R_1`$ and $`|R_2`$, so that $`L_1|`$ and $`|R_1`$ become almost orthogonal to each other. To still satisfy the normalization condition $`L_n|R_n=1`$ in this case, the eigenfunction should be scaled up, thus making the spectral structure very sensitive to the numerical error involved in the calculation of eigenfunctions. For an interesting discussion on this point, one may refer to the work by Nelson and co-workers. Due to this numerical instability, the use of the spectral method as a density matrix propagation scheme is not without limitation.

Figure 3a shows the spectral structure and the time-evolution of the density matrix propagation for the case of $`\beta \mathrm{\Omega }=1`$, $`\beta \lambda =15`$, $`\beta V=12`$, and $`\beta ϵ=5`$. Generally, when the energy bias is small ($`\beta ϵ5`$), the left and right eigenfunctions can form a complete orthonormal basis set, so the spectral method is stable and can be used as a numerical propagation method for the density matrix. With a large energy bias, however, the calculated eigenfunctions may not form a complete orthonormal basis. To model for the photo-induced back electron transfer experiment in the mixed-valence compounds the initial density matrix is chosen as a thermal equilibrium distribution of the donor state(i.e. 1-state) pumped to the acceptor state(i.e. 2-state),

$`\rho _i(E,0)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }\mathrm{\Delta }_E}}\mathrm{exp}\left({\displaystyle \frac{(E+\lambda )^2}{2\mathrm{\Delta }_E^2}}\right)\delta _{i2},`$ (62)
$`\rho _{12}(E,0)`$ $`=`$ $`\rho _{21}(E,0)=0.`$ (63)
It would be straightforward to calculate the spatial distribution of the density matrix in time $`\rho (E,t)`$ by applying the propagator in Eq. (54) to the above initial density matrix; however, to demonstrate the overall temporal behavior only the time evolution of the total population in the acceptor state is calculated,

$`P_2(t)={\displaystyle 𝑑E\rho _2(E,t)}.`$ (64)
In order to check the validity of the spectral method as a propagation scheme in this case, we also calculated the time evolution of the density matrix by directly solving the $`4N`$ differential equations for the expansion coefficients of the density matrix using the Bulirsh-Stoer algorithm, and the comparison in Fig. 3(a) shows a perfect agreement. If only the transient behavior is concerned with, the direct propagation method would be preferred over the spectral method, however, the spectral propagation has the advantage when calculating the long time behavior once the complete spectrum is known. Overall, the computational costs for two method are comparable to each other. As expected from the spectral structure shown in the previous section the population in the acceptor state shows an underdamped coherent oscillation behavior at initial times followed by a damped oscillation behavior at later times.

Further, we have also studied the density matrix propagation for different energy biases to examine the electronic dephasing effect. As seen from Fig. 4(a), the increase in energy bias destroys the electronic coherence dramatically. Another interesting observation is the phase shift in the population dynamics as the energy bias is varied, and it is because the Rabi oscillation frequency increases with energy bias. We can confirm the temporal behavior of the density matrix propagation by examining the spectral structure shown in Fig. 4(b). The period of the initial coherence is estimated to be $`\tau _{osc}0.25\beta `$ from Fig. 4(a). In comparison, the Rabi frequency for the corresponding adiabatic two-level system is given by $`\mathrm{\Omega }_R=\sqrt{ϵ^2+4V^2}25\beta ^1`$, which can also be obtained from the onset of imaginary branches in the eigen-tree shown in Fig. 4(b), and the estimation is consistent with the oscillation period observed in the dynamics since $`\tau _{osc}2\pi /\mathrm{\Omega }_R`$. The real eigenvalues of the lowest excited states in the the imaginary branches are estimated to be $`\beta k12`$, and they agree with the decay time of the oscillation amplitude in Fig. 4(a), confirming the validity of the spectral method as a density matrix propagation scheme. Even though it has been well known in the literature that the damping of population is enhanced with increased energy asymmetry, we have also confirmed this through the spectral analysis method.

As an example of the eigenfunction responsible for the coherent oscillation behavior observed in Fig. 4 (b), we show the left and right eigenfunctions corresponding to a complex eigenvalue $`\beta Z=2.6228\pm i26.394`$ for a symmetric case $`(\beta ϵ=0)`$ and $`\beta Z=2.8057\pm i26.466`$ for an asymmetric case$`(\beta ϵ=5)`$ in Figs. 5 and 6. The eigenfunctions corresponding to a complex conjugate pair of eigenvalues are also complex conjugate to each other; therefore, the frequency spectrum of the density matrix evolution is proportional to the norm of wavefunction. We note that the left eigenfunction is more extended than the right eigenfunction. Although the population distribution in the donor and acceptor states corresponding to coherent oscillation is inverted with respect to the Boltzmann distribution, it does not contribute to the steady-state population distribution due to the transient nature.

## IV CONCLUDING REMARKS

In this paper we have applied the spectral analysis method to the non-adiabatic two-state diffusion equation, that describes electron transfer dynamics in Debye solvents. In particular, we have examined electronic coherence in mixed-valence compounds, and demonstrated that underdamped Rabi oscillations are observed in an overdamped solvent environment. Detailed study of the spectral structure of the non-adiabatic operator for various energy biases and coupling constants allows us to determine the underlying mechanisms of electron transfer kinetics. Eigenvalues form three branches in the eigen-diagram: a single branch of real eigenvalues and two symmetric branches of complex conjugate eigenvalues. In strongly coupled systems, all three branches have a similar order of magnitude, indicating that both multiple-exponential decay and coherent oscillations can be observed experimentally.

We have investigated the dependence of the spectral structure on the coupling constant. In the very weak coupling regime, the lowest excited state is well separated from higher states, which makes the electron transfer dynamics a well-defined rate process. In the strong coupling regime, however, the eigenvalue diagram shows coalescence/bifurcation behavior in the complex plane. We have used the spectral method to calculate the time-evolution of the density matrix, and indeed, observed electronic coherence in the temporal behavior of population in the acceptor state for non-equilibrium initial distributions. We also found a good agreement between results of the spectral propagation method and of the numerical propagation method for small energy bias cases. Due to non-Hermitianity of the non-adiabatic operator, the spectral propagation method was not numerically stable for large energy bias cases.

For an isolated quantum system, the eigen-solution to the Schrödinger equation completely determines its dynamics. In a similar fashion, the eigen-solution to the non-adiabatic diffusion operator completely characterizes the dynamics of a dissipative system and thus provides a powerful tool to analyze dissipative dynamics. It is well known that quantum dynamics comes from the underlying spectra, especially in gas-phase chemical systems; however, the spectral aspect of condensed phase dissipative systems has not been well recognized yet and deserves further investigation. Though the analysis presented here is restricted to semi-classical dissipative systems, it may also be applied to quantum dissipative dynamics. In principle, we can derive the evolution equation for quantum dissipative systems either from first principles or through numerical reduction, and then pose the quantum dissipative equation of motion as an eigen-value problem. Along this line, the dissipative dynamics of the spin-boson Hamiltonian, which has been studied mostly as a dynamic problem, can also be explored as a spectral problem in the future.

## ACKNOWLEDGMENTS

The authors would like to thank NSF for financial support. One of us (YJ) would like to thank the Korean Foundations for Advanced Studies for financial support.