# The Effect of Duschinskii Rotations on Spin-Dependent Electron Transfer
Dynamics

Suraj S. Chandran, Yanze Wu, Joseph E. Subotnik

###### Abstract

We investigate spin-dependent electron transfer in the presence of a
Duschinskii rotation. In particular, we propagate dynamics for a two-level
model system for which spin-orbit coupling introduces an interstate coupling
of the form $e^{iWx}$, which is both position($x$)-dependent and complex-
valued. We demonstrate that two-level systems coupled to Brownian oscillators
with Duschinskii rotations (and thus entangled normal modes) can produce
marked increases in transient spin polarization relative to two-level systems
coupled to simple shifted harmonic oscillators. These conclusions should have
significant relevance for modeling the effect of nuclear motion on chiral
induced spin selectivity.

##  1 Introduction

Electron transfer is a crucial phenomenon underlying many aspects of
chemistry, highlighting the fundamental energy conversion between electronic
and nuclear degrees of freedom. As such, there has been an enormous push
within the physical chemistry and chemical physics communities over the last
60 years to explore a variety of aspects related to electron transfer,
including barriers to solution-phase thermal electron transfer [1, 2],
transient photo-induced electron transfer dynamics [3], electron-hole pair
generation at metal interfaces[4, 5, 6, 7], and much more.

The basic model for such processes is the shifted harmonic oscillator (SHO)
model[8, 9] whereby we imagine two coupled harmonic surfaces, identical in all
respects except for one: the well minima are shifted in space and found at
different energies. The relevant Hamiltonian can be written as follows[10]:

$\displaystyle\hat{H}$ $\displaystyle=\hat{H}_{nuc}+\hat{H}_{el}+\hat{V}$
(1.1) $\displaystyle\hat{H}_{nuc}$
$\displaystyle=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\frac{1}{2}\hat{\textbf{q}}^{T}\hat{\bm{\Omega}}^{2}\hat{\textbf{q}}$
(1.2) $\displaystyle\hat{H}_{el}$ $\displaystyle=\hat{H}_{D}\lvert
d\rangle\langle d\rvert+\hat{H}_{A}\lvert a\rangle\langle a\rvert$ (1.3)
$\displaystyle\hat{H}_{D}$
$\displaystyle=E_{D}+\bm{c}_{D}^{T}\hat{\textbf{q}}$ (1.4)
$\displaystyle\hat{H}_{A}$
$\displaystyle=E_{A}+\bm{c}_{A}^{T}\hat{\textbf{q}}$ (1.5)
$\displaystyle\hat{V}$ $\displaystyle=V\lvert d\rangle\langle
a\rvert+V^{*}\lvert a\rangle\langle d\rvert$ (1.6)

Here, $\hat{H}_{nuc}$ and $\hat{H}_{el}$ respectively represent nuclear and
electronic Hamiltonians. In $\hat{H}_{nuc}$, $\hat{\textbf{p}}$ and
$\hat{\textbf{q}}$ are normal-mode momentum and position vector operators, and
$\hat{\bm{\Omega}}^{2}$ is the (diagonal) harmonic matrix. In $\hat{H}_{el}$,
$\lvert d\rangle$ and $\lvert a\rangle$ respectively represent donor and
acceptor electronic states. The two states are coupled together by $V$ and
coupled to the nuclear coordinates via the coupling vectors $\bm{c}_{D}$ and
$\bm{c}_{A}$. Taken together, Eqs. 1.1-1.6 embody the SHO model, otherwise
known as the spin-boson model.

The spin-boson model has captivated researchers for decades, owing to the rich
dynamics that emerge from its apparently simple construction[11]. One regime
of particular interest is that where the coupling $V$ between electronic
states is small: the nonadiabatic regime[10, 12]. In a typical treatment of
the nonadiabatic spin-boson model, one takes the initial state to be
equilibrated in the donor well and applies perturbation theory to find the
first-order correction to the acceptor-state population. In other words, the
nuclear density matrix at time 0 is taken to be of the form

$\displaystyle\hat{\rho}(0)$
$\displaystyle=\frac{1}{Z}e^{-\beta(\hat{H}_{nuc}+\hat{H}_{D})}$ (1.7)
$\displaystyle Z$
$\displaystyle=\text{Tr}\left[e^{-\beta(\hat{H}_{nuc}+\hat{H}_{D})}\right]$
(1.8)

In the high-temperature limit, the resulting rate of transfer to the acceptor
well is given by the classic Marcus rate expression:

$k=\frac{2\pi\lvert V\rvert^{2}}{\sqrt{4\pi
E_{r}k_{B}T}}\exp\left(-\frac{(\Delta G+E_{r})^{2}}{4E_{r}k_{B}T}\right)$
(1.9)

Here, $E_{r}$ is the reorganization energy, which defines the energy required
to distort the nuclear geometry from the donor well minimum to the acceptor
well minimum without an electronic transition.

Of course, the nonadiabatic regime is only one limit of the model; the spin-
boson model can be used to investigate many other regimes as well. For
example, one might ask how the dynamics are different if the interstate
coupling is large, a regime which is termed the adiabatic limit[13].
Alternatively, one might interrogate the solvent-controlled regime[14, 15,
16], where strong interactions between the system and bath lead to heavy
damping of the nuclear motions that produce electronic transitions. As yet
another example, the spin-boson model can be extended beyond the basic
formulation above to a model where the interstate coupling is dependent on
nuclear coordinates, breaking the so-called Condon approximation[17, 18].

In recent years, the discovery of the chiral induced spin selectivity (CISS)
effect has begun to motivate investigation into yet another fundamental aspect
of electron transfer[19]: the realm of spin. A host of experiments within
recent decades have shown that electron currents through chiral organic
molecules can be strongly spin polarized[19, 20]. In this context, new
questions arise: how can condensed-phase electron transfer exhibit such
dramatic spin selectivity given the disproportionately small magnitude of
organic spin-orbit couplings? Might such polarization be tied to nuclear
geometry and motion[21, 22, 23] and, if so, how large can the resulting spin
polarization be? And finally, in the spirit of the last few paragraphs, can
the spin-boson model be an effective theoretical setting for probing such
phenomena?

If we indeed seek to use the spin-boson model as a template for addressing the
questions above, we will need to make one crucial modification. Namely, in
order to capture the effects of SOC, we will need to explore regimes where the
Hamiltonian is not real-valued, but rather complex-valued (i.e. where the off-
diagonal coupling elements in Eq. 1.6 have non-zero imaginary components). To
understand why this is so, consider the simple case of a intramolecular
crossing between two doublets. This is a four-state system where a basis can
be written as $\\{\lvert 1,\uparrow\rangle,\lvert 1,\downarrow\rangle,\lvert
2,\uparrow\rangle,\lvert 2,\downarrow\rangle\\}$; in other words, we consider
states formed by a direct product between two electronic molecular orbitals
and two spin states. In this basis, the Hamiltonian can be written as follows
(with $\hat{\textbf{p}}$ and $\hat{\textbf{q}}$ respectively denoting nuclear
momenta and positions):

$\displaystyle\hat{H}$
$\displaystyle=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\hat{\mathcal{V}}(\hat{\textbf{q}})+\begin{bmatrix}E_{1}&0&V&0\\\
0&E_{1}&0&V\\\ V&0&E_{2}&0\\\
0&V&0&E_{2}\end{bmatrix}+\hat{\textbf{L}}\cdot\hat{\textbf{S}}$ (1.10)
$\displaystyle\hat{\textbf{L}}\cdot\hat{\textbf{S}}$
$\displaystyle=\begin{bmatrix}0&0&\langle
1|\hat{\textbf{L}}_{z}|2\rangle&\langle
1|\hat{\textbf{L}}_{x}-i\hat{\textbf{L}}_{y}|2\rangle\\\ 0&0&\langle
1|\hat{\textbf{L}}_{x}+i\hat{\textbf{L}}_{y}|2\rangle&-\langle
1|\hat{\textbf{L}}_{z}|2\rangle\\\ \langle
2|\hat{\textbf{L}}_{z}|1\rangle&\langle
2|\hat{\textbf{L}}_{x}-i\hat{\textbf{L}}_{y}|1\rangle&0&0\\\ \langle
2|\hat{\textbf{L}}_{x}+i\hat{\textbf{L}}_{y}|1\rangle&-\langle
2|\hat{\textbf{L}}_{z}|1\rangle&0&0\end{bmatrix}$ (1.11)

In general, $\hat{\textbf{L}}_{x}$, $\hat{\textbf{L}}_{y}$, and
$\hat{\textbf{L}}_{z}$ are all functions of $\hat{\textbf{q}}$. Notice that
the $\hat{\textbf{L}}$ operators are proportional to momentum, so the matrix
elements $\langle 1|\hat{\textbf{L}}|2\rangle$ are purely imaginary.
Therefore, the two nonzero blocks of $\hat{\textbf{L}}\cdot\hat{\textbf{S}}$
are each anti-Hermitian and there exists a $\hat{\textbf{q}}$-dependent spin
basis
$\\{\tilde{\uparrow}(\hat{\textbf{q}}),\tilde{\downarrow}(\hat{\textbf{q}})\\}$
in which the blocks are diagonal with imaginary eigenvalues. In this basis,
the Hamiltonian reads exactly as in Eq. 1.10, except that the spin-orbit
operator now takes the form

$\hat{\textbf{L}}\cdot\hat{\mathcal{S}}=\begin{bmatrix}0&0&i\alpha(\hat{\textbf{q}})&0\\\
0&0&0&-i\alpha(\hat{\textbf{q}})\\\ -i\alpha(\hat{\textbf{q}})&0&0&0\\\
0&i\alpha(\hat{\textbf{q}})&0&0\end{bmatrix}$ (1.12)

Here, $\alpha$ is a real-valued function of $\hat{\textbf{q}}$. If we now make
the assumption that the spin basis is independent of $\hat{\textbf{q}}$ (i.e.
the Hamiltonian is spin-conserving) and project the Hamiltonian onto the two
spin-states, we obtain two Hamiltonians, one of which governs the evolution of
$\lvert\tilde{\uparrow}\rangle$ states, and another which governs the
evolution of $\lvert\tilde{\downarrow}\rangle$ states:

$\displaystyle\hat{H}(\tilde{\uparrow})$
$\displaystyle=\langle\tilde{\uparrow}|\hat{H}|\tilde{\uparrow}\rangle=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\hat{\mathcal{V}}(\hat{\textbf{q}})+\begin{bmatrix}E_{1}&V+i\alpha(\hat{\textbf{q}})\\\
V-i\alpha(\hat{\textbf{q}})&E_{2}\\\ \end{bmatrix}$ (1.13)
$\displaystyle\hat{H}(\tilde{\downarrow})$
$\displaystyle=\langle\tilde{\downarrow}|\hat{H}|\tilde{\downarrow}\rangle=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\hat{\mathcal{V}}(\hat{\textbf{q}})+\begin{bmatrix}E_{1}&V-i\alpha(\hat{\textbf{q}})\\\
V+i\alpha(\hat{\textbf{q}})&E_{2}\\\
\end{bmatrix}=\hat{H}(\tilde{\uparrow})^{*}$ (1.14)

In other words, for a spin-conserving process occurring in the presence of
spin-orbit coupling, spin-up and spin-down electrons evolve according to
different complex-valued Hamiltonians: if the former evolve according to
$\hat{H}$, the latter evolve according to $\hat{H}^{*}$.

Now what are the implications of propagating a complex-valued Hamiltonian? At
first glance, one might expect equivalent dynamics to arise regardless of
whether one propagates using $\hat{H}$ or $\hat{H}^{*}$. After all, the
eigenvalues of both operators are exactly the same. Indeed, if we consider
equilibrium Fermi Golden Rule rates from donor to acceptor for a general
vibronic system, it is clear that there can be no difference in $\hat{H}$ vs.
$\hat{H}^{*}$ dynamics:

$\displaystyle k_{d\rightarrow a}(\hat{H})$
$\displaystyle=\sum_{\textbf{v},\textbf{v}^{\prime}}P_{\textbf{v}}\lvert\langle
a,\textbf{v}^{\prime}|\hat{V}|d,\textbf{v}\rangle\rvert^{2}\delta(E_{a,\textbf{v}^{\prime}}-E_{d,\textbf{v}})$
(1.15) $\displaystyle k_{d\rightarrow a}(\hat{H}^{*})$
$\displaystyle=\sum_{\textbf{v},\textbf{v}^{\prime}}P_{\textbf{v}}\lvert\langle
a,\textbf{v}^{\prime}|\hat{V}^{*}|d,\textbf{v}\rangle\rvert^{2}\delta(E_{a,\textbf{v}^{\prime}}-E_{d,\textbf{v}})=k_{d\rightarrow
a}(\hat{H})$ (1.16)

Here, $\textbf{v}^{\prime}$ and v respectively index the acceptor and donor
vibrational eigenstates, and $P_{\textbf{v}}=\frac{1}{Z}\exp(-\beta
E_{d,\textbf{v}})$ with $Z$ defined as in 1.8.

However, upon reflection, one should realize that these conclusions hold only
to the lowest order of perturbation theory; furthermore, they rely on the
initial conditions being those of equilibrium[24]. Beyond these limits, a
subtle difference in the properties of $\hat{H}$ and $\hat{H}^{*}$ becomes
critically important; namely the fact that the respective eigenvectors differ
by complex conjugation. From a semiclassical (rather than purely quantum)
perspective, this subtle difference can lead to very different so-called
“Berry forces” for $\hat{H}$ vs $\hat{H}^{*}$ trajectories. As a review for
the reader, in the adiabatic limit, semiclassical trajectories on adiabat $n$
experience a Berry force given by the following equation[25, 26], :

$\textbf{F}_{B}^{(n)}=i\hbar\dot{\textbf{R}}\times(\nabla_{\textbf{R}}\times\textbf{D}_{nn})$
(1.17)

Here, $\nabla_{\textbf{R}}$ is the gradient operator with respect to nuclear
coordinates and $\textbf{D}_{nn}$ is the derivative coupling from adiabatic
surface $n$ to itself. Note that the Berry force vector is orthogonal to the
momentum of the nuclei ($\textbf{F}_{B}\cdot\dot{\textbf{R}}=0$) and can be
understood semiclassically as an effective magnetic field that acts on nuclei
near crossing points. Note also that, according to Eq. 1.17, it follows that
the adiabatic Berry forces for $\hat{H}$ are equal and opposite to those for
$\hat{H}^{*}$. which should appear reasonable insofar as flipping the spin of
the relevant electron(s) effectively reverses all magnetic fields in the
system. The take home point is that, if we run nuclera trajectories for
systems with spin up and spin down electrons, we will find different dynamics;
the magnitude of these differences will necessarily depend on many factors,
including the topography of the potential energy surfaces as well as the size
of the Berry force.

At this point, a few words are appropriate regarding the size and meaning of
the Berry force. In general, one should apply Eq. 1.17 to a trajectory only in
the adiabatic limit; Eq. 1.17 is derivable only when there is a large energy
gap between adiabatic electronic states. That being said, the astute reader
will notice that in the nonadiabatic regime, the derivative couplings (and
therefore Eq. 1.17) happen to be very large near electronic crossing points;
one is then inevitably compelled to ask whether or not such large matrix
elements might have a significant effect on entangled nuclear-electronic-spin
quantum dynamics. Might such Berry forces (arising from Hamiltonians with
small SOC matrix elements) have large effects on dynamics and therefore (at
least partly) mediate the CISS effect? Unfortunately, at present, very little
is known about the meaning of the forces in Eq. 1.17 in the nonadiabatic
regime. To that end, our research group has recently begun constructing novel,
nonadiabatic surface hopping approaches with the eventual goal of treating ab
initio nuclear dynamics in the presence of a set of degenerate electronic spin
states.

Looking forward, until such a semiclassical ab initio algorithm is
available[27, 28], we believe the best approach to evaluating the size and
impact of Berry forces is to develop an analytic model where intuition can be
easily gleaned and relied upon. To that end, in Ref. [29], we studied a spin-
boson model with an Exponential Spin-Orbit Coupling (ESOC). For such a
problem, the Hamiltonian is defined exactly as in Eq. 1.1 with the only change
being that the interstate coupling operator $\hat{V}$ takes the form complex-
valued non-Condon form:

$\hat{V}=Ve^{i\textbf{W}^{T}\hat{\textbf{q}}}\lvert d\rangle\langle
a\rvert+V^{*}e^{-i\textbf{W}^{T}\hat{\textbf{q}}}\lvert a\rangle\langle
d\rvert$ (1.18)

Here, W can be understood as encoding the gradient of the spin-orbit coupling
near a diabatic crossing point. To illustrate this, we expand the coupling
operator in Eq. 1.18 to first order in $\hat{\textbf{q}}$ around a crossing
point $\textbf{q}_{0}$:

$Ve^{i\textbf{W}^{T}\hat{\textbf{q}}}\simeq
Ve^{i\textbf{W}^{T}\textbf{q}_{0}}(1+i\textbf{W}^{T}(\hat{\textbf{q}}-\textbf{q}_{0}))$
(1.19)

Eq. 1.19 can be directly compared with $V+i\alpha(\hat{\textbf{q}})$ from Eq.
1.13, justifying the use of this Hamiltonian.

Using the Hamiltonian defined by Eqs. 1.1 and 1.18, we demonstrated in Ref.
[29] that spin polarization could indeed emerge dynamically as a consequence
of nonequilibrium initial conditions. We note, however, that these conclusions
required broken mirror symmetry; if a system is mirror symmetric, then the
equal and opposite Berry forces produced in the respective $\hat{H}$ and
$\hat{H}^{*}$ systems produce equivalent, albeit mirrored, nuclear dynamics
(see below). These results would appear consistent with the fact that
experimentally, for the CISS effect, one must run current through a chiral
molecule in order to achive spin separation in general (and of course, chiral
molecules, by definition, have no mirror symmetry).

With this background in mind, the goal of the present article is to gain
intuition about the about the extent of the expected transient spin-polarized
electron transfer that will occur in the absence of mirror symmetry when we go
beyond the constraints of the SHO model. In particular, we will study if/how
the presence of intertwined vibrational modes (as modeled with a Duchinskii
rotation) affects spin selectivity. To achieve this goal, we will investigate
the ESOC model in the context of an intuitive model electron-transfer system
where two coupled primary modes are each interacting with a large dissipative
thermal bath [30]. Working within this framework, we will answer two
questions. First, how does the direction of the intramolecular spin-orbit
coupling vector influence the spin-polarization of the final state? Here, the
partitioning between primary and bath modes proposed above is especially
useful, since the primary modes can be thought of as defining a ”molecule of
interest” within our very large system. In this context, the concepts of
intramolecular and intermolecular interactions are well-defined; any coupling
to the primary modes is an intramolecular interaction, whereas any coupling to
the bath modes is an intermolecular interaction. Second, can systems in which
normal mode frequencies change upon electronic transition produce much
stronger spin polarization than simple Shifted Harmonic Oscillator (SHO) model
systems, and if so, when are these effects maximal? Answering these questions
will inevitably give us a great deal of insight into the coupled nature of
nuclear-electronic-spin dynamics.

An outline of this article is as follows. First, in Sec. 2, we present the
relevant theoretical background needed to construct the aforementioned model
system and to model its dynamics using a Duschinskii rotation formalism [31,
32]. Next, in Sec. 3 we present results on how changing the geometry of a two-
state harmonic system via Duschinskii rotations affects spin polarization.
Finally, in Sec. 4 we discuss these results in the context of mirror symmetry
and chirality.

##  2 Theory

###  2.1 The Definition of Duschinskii Rotations for a General Two-State
Quadratic Hamiltonian

Consider a general quadratic Hamiltonian with a potential of the following
form:

$\displaystyle\hat{\mathcal{V}}=\begin{bmatrix}\frac{1}{2}\hat{\textbf{q}}^{T}\hat{\bm{\Omega}}_{g}^{2}\hat{\textbf{q}}+\bm{\lambda}_{g}^{T}\hat{\textbf{q}}+E_{g}&\hat{V}(\hat{\textbf{q}})\\\
\hat{V}(\hat{\textbf{q}})^{\dagger}&\frac{1}{2}\hat{\textbf{q}}^{T}\hat{\bm{\Omega}}_{e}^{2}\hat{\textbf{q}}+\bm{\lambda}_{e}^{T}\hat{\textbf{q}}+E_{e}\end{bmatrix}$
(2.1)

Here, $\hat{\textbf{q}}$ is a vector operator consisting of the position
operators for each mode in the system and $g$, $e$ respectively denote the
ground and excited states. The Hamiltonian defined by Eq. 2.1 is a generalized
spin-boson model where the harmonic surfaces of the ground and excited state
diabats are not required to be identical. The two diabatic harmonic matrices
$\hat{\bm{\Omega}}_{g}^{2}$ and $\hat{\bm{\Omega}}_{e}^{2}$ are symmetric and
can be diagonalized respectively by orthogonal matrices $\textbf{S}_{g}$ and
$\textbf{S}_{e}$:

$\displaystyle\hat{\bm{\Omega}}_{g}$
$\displaystyle=\textbf{S}_{g}\hat{\bm{\Delta}}_{g}\textbf{S}_{g}^{T}$ (2.2)
$\displaystyle\hat{\bm{\Omega}}_{e}$
$\displaystyle=\textbf{S}_{e}\hat{\bm{\Delta}}_{e}\textbf{S}_{e}^{T}$ (2.3)

Therefore, we can write the diabatic Hamiltonians as follows:

$\displaystyle\frac{1}{2}\hat{\textbf{q}}^{T}\hat{\bm{\Omega}}_{g}^{2}\hat{\textbf{q}}+{\bm{\lambda}}_{g}^{T}\hat{\textbf{q}}$
$\displaystyle=\frac{1}{2}(\hat{\textbf{q}}^{T}+{\bm{\lambda}}_{g}^{T}\hat{\bm{\Omega}}_{g}^{-2})\textbf{S}_{g}\hat{\bm{\Delta}}_{g}^{2}\textbf{S}_{g}^{T}(\hat{\textbf{q}}+\hat{\bm{\Omega}}_{g}^{-2}{\bm{\lambda}}_{g})-\frac{1}{2}{\bm{\lambda}}_{g}^{T}\hat{\bm{\Omega}}_{g}^{-2}{\bm{\lambda}}_{g}+E_{g}$
(2.4)
$\displaystyle\frac{1}{2}\hat{\textbf{q}}^{T}\hat{\bm{\Omega}}_{e}^{2}\hat{\textbf{q}}+{\bm{\lambda}}_{e}^{T}\hat{\textbf{q}}$
$\displaystyle=\frac{1}{2}(\hat{\textbf{q}}^{T}+{\bm{\lambda}}_{e}^{T}\hat{\bm{\Omega}}_{e}^{-2})\textbf{S}_{e}\hat{\bm{\Delta}}_{e}^{2}\textbf{S}_{e}^{T}(\hat{\textbf{q}}+\hat{\bm{\Omega}}_{e}^{-2}{\bm{\lambda}}_{e})-\frac{1}{2}{\bm{\lambda}}_{e}^{T}\hat{\bm{\Omega}}_{e}^{-2}{\bm{\lambda}}_{e}+E_{e}$
(2.5)

The matrices $\hat{\bm{\Delta}}$ are diagonal operators, and they encode the
normal mode frequencies for the upper and lower diabats. Now, let us
substitute
$\hat{\textbf{q}}^{\prime}=\textbf{S}_{g}^{T}(\hat{\textbf{q}}+\hat{\bm{\Omega}}_{g}^{-2}{\bm{\lambda}}_{g})$
and
$\hat{\textbf{q}}^{\prime\prime}=\textbf{S}_{e}^{T}(\hat{\textbf{q}}+\hat{\bm{\Omega}}_{e}^{-2}{\bm{\lambda}}_{e})$.
In this new basis, the two blocks become:

$\displaystyle\frac{1}{2}\hat{\textbf{q}}^{\prime
T}\hat{\bm{\Delta}}_{g}^{2}\hat{\textbf{q}}^{\prime}-\frac{1}{2}{\bm{\lambda}}_{g}^{T}\hat{\bm{\Omega}}_{g}^{-2}{\bm{\lambda}}_{g}+E_{g}$
(2.6) $\displaystyle\frac{1}{2}\hat{\textbf{q}}^{\prime\prime
T}\hat{\bm{\Delta}}_{e}^{2}\hat{\textbf{q}}^{\prime\prime}-\frac{1}{2}{\bm{\lambda}}_{e}^{T}\hat{\bm{\Omega}}_{e}^{-2}{\bm{\lambda}}_{e}+E_{e}$
(2.7)

Also note that the position operators $\textbf{q}^{\prime}$ and
$\textbf{q}^{\prime\prime}$ are related by the following transformation:

$\hat{\textbf{q}}^{\prime}=\textbf{S}_{g}^{T}\textbf{S}_{e}\hat{\textbf{q}}^{\prime\prime}+\textbf{S}_{g}^{T}(\hat{\bm{\Omega}}_{g}^{-2}{\bm{\lambda}}_{g}-\hat{\bm{\Omega}}_{e}^{-2}{\bm{\lambda}}_{e})$
(2.8)

The ground state normal mode vector $\hat{\textbf{q}}^{\prime}$ is of the form
$\textbf{S}\hat{\textbf{q}}^{\prime\prime}+\textbf{d}$. The matrix
$\textbf{S}=\textbf{S}_{g}^{T}\textbf{S}_{e}$ is known as a Duschinskii
rotation matrix [31] and describes how the normal modes of one diabat can be
expressed as linear combinations of normal modes of the other diabat. In a
shifted harmonic oscillator system, the Duschinskii matrix would simply be the
identity since both diabats would have identical normal modes. The vector
$\textbf{d}=\textbf{S}_{g}^{T}(\hat{\bm{\Omega}}_{g}^{-2}\bm{\lambda}_{g}-\hat{\bm{\Omega}}_{e}^{-2}\bm{\lambda}_{e})$
in Eq. 2.8 is the shift vector, and describes the displacement between the
well minima of the two diabats. Note that, under the transformation in Eqs.
2.6 and 2.7, $\hat{V}(\hat{\textbf{q}})$ is transformed to
$\hat{V}(\textbf{S}_{e}\hat{\textbf{q}}^{\prime\prime}-\hat{\bm{\Omega}}_{e}^{-2}\bm{\lambda}_{e})$;
in particular, if $\hat{V}=e^{i\textbf{W}^{T}\hat{\textbf{q}}}$, $\hat{V}$ is
mapped to $e^{i\textbf{W}^{T}\textbf{S}_{e}\hat{\textbf{q}}^{\prime\prime}}$
(up to a constant complex phase). Finally, note also that energy difference
between the ground and excited state minima ($\Delta G$) can be written as
$E_{g}-E_{e}+\frac{1}{2}\left(\bm{\lambda}_{e}^{T}\hat{\bm{\Omega}}_{e}^{-2}\bm{\lambda}_{e}-\bm{\lambda}_{g}^{T}\hat{\bm{\Omega}}_{g}^{-2}\bm{\lambda}_{g}\right)$.

We will use the parametrization in Eqs. 2.6 and 2.7 to model the dynamics of
molecular electron transfer systems for which the vibrational modes may or may
not change over the course of the event.

###  2.2 Dynamics of Systems with Duschinskii Rotations

We have applied the Fermi Golden Rule formalism to the Hamitonian above. The
lengthy calculation [31] is provided in the appendix. We focus on a
Hamiltonian of the following form:

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=\begin{bmatrix}\hat{H}_{B}^{(g)}&Ve^{i\textbf{W}^{T}\hat{\textbf{x}}}\\\
V^{*}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}&\hat{H}_{B}^{(e)}\end{bmatrix}$
(2.9) $\displaystyle\hat{H}_{B}^{(g)}$
$\displaystyle=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\frac{1}{2}\left(\textbf{S}\hat{\textbf{x}}+\textbf{d}\right)^{T}\hat{\bm{\Omega}}_{g}^{2}\left(\textbf{S}\hat{\textbf{x}}+\textbf{d}\right)+\Delta
G$ (2.10) $\displaystyle\hat{H}_{B}^{(e)}$
$\displaystyle=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\frac{1}{2}\hat{\textbf{x}}^{T}\hat{\bm{\Omega}}_{e}^{2}\hat{\textbf{x}}$
(2.11)

here, S is the Duschinskii rotation matrix, d is the displacement vector
between the well minima, and the harmonic matrices $\hat{\bm{\Omega}}_{e}$ and
$\hat{\bm{\Omega}}_{g}$ are assumed to be diagonal. Before proceeding further,
it will be helpful to define several matrix functions which will appear below:

$\bm{a}_{e}(t)=\hat{\bm{\Omega}}_{e}[\sin(\hat{\bm{\Omega}}_{e}t)]^{-1}$
(2.12)
$\bm{a}_{g}(t)=\hat{\bm{\Omega}}_{g}[\sin(\hat{\bm{\Omega}}_{g}t)]^{-1}$
(2.13)
$\bm{b}_{e}(t)=\hat{\bm{\Omega}}_{e}[\tan(\hat{\bm{\Omega}}_{e}t)]^{-1}$
(2.14)
$\bm{b}_{g}(t)=\hat{\bm{\Omega}}_{g}[\tan(\hat{\bm{\Omega}}_{g}t)]^{-1}$
(2.15)
$\textbf{A}=\bm{a}_{e}(-\tau-i\beta)+\textbf{S}^{T}\bm{a}_{g}(\tau)\textbf{S}$
(2.16)
$\textbf{B}=\bm{b}_{e}(-\tau-i\beta)+\textbf{S}^{T}\bm{b}_{g}(\tau)\textbf{S}$
(2.17) $\textbf{E}=\bm{b}_{e}(-\tau-i\beta)-\bm{a}_{e}(-\tau-i\beta)$ (2.18)
$\textbf{G}=\bm{b}_{g}(\tau)-\bm{a}_{g}(\tau)$ (2.19)
$\textbf{G}_{\beta}=\bm{b}_{g}(-i\beta)-\bm{a}_{g}(-i\beta)$ (2.20)

####  2.2.1 Equilibrium Dynamics

Let us begin by initializing the system on the excited state with the
equilibrium density matrix

$\displaystyle\hat{\rho}(0)$
$\displaystyle=\frac{1}{Z}e^{-\beta\hat{H}_{B}^{(e)}}\lvert e\rangle\langle
e\rvert$ (2.21) $\displaystyle Z$
$\displaystyle=\text{Tr}_{B}\left[e^{-\beta\hat{H}_{B}^{(e)}}\right]$ (2.22)

The resulting correlation function that defines the rate of electron transfer
to the ground state of the system is given by (see App. 7.1):

$\displaystyle C(\tau)$
$\displaystyle=\langle\hat{V}_{eg}(\tau)\hat{V}_{ge}\rangle$ (2.23)
$\displaystyle=\frac{1}{Z}\text{Tr}_{B}\left[e^{-i\hat{H}_{B}^{(e)}(-\tau-i\beta)}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(g)}\tau}e^{i\textbf{W}^{T}\hat{\textbf{x}}}\right]$
(2.24)
$\displaystyle=\frac{1}{Z}\sqrt{\frac{\det(\bm{a}_{e}(-\tau-i\beta))\det(\bm{a}_{g}(\tau))}{\det(\textbf{B})\det(\textbf{B}-\textbf{AB}^{-1}\textbf{A})}}\exp\left(i\left[\textbf{d}^{T}\textbf{GS}(\textbf{B}-\textbf{A})^{-1}\textbf{ES}^{T}\textbf{d}-\textbf{W}^{T}(\textbf{B}+\textbf{A})^{-1}\textbf{W}\right]\right)$
(2.25)

Notably, the term in the exponential has a bilinear dependence on W, meaning
that the sign of W has no bearing on the equilibrium dynamics. This suggests
that if there were to be a dynamical dependence on the sign of W for
equilibrium dynamics, such an effect should be quite small because it could
reveal itself only at higher orders of perturbation theory. Following standard
perturbation theory, the ground state population at time $t$ is given by:

$P_{g}(t)\simeq 2\lvert
V\rvert^{2}\cdot\text{Re}\left[\int_{0}^{t}dt^{\prime}\,\int_{0}^{t^{\prime}}d\tau\,e^{-i\Delta
G\tau}C(\tau)\right]$ (2.26)

where $\Delta G$ is defined in Eq. 2.10.

####  2.2.2 Nonequilibrium Dynamics

Our primary motivation in developing a nonequilibrium formalism for
Duschinskii systems is to understand the consequences of Duschinskii rotations
with regards to photoexcitation dynamics. To that end, let us next focus on
the case where the system begins in the excited state with the density matrix

$\displaystyle\hat{\rho}_{neq}(0)$
$\displaystyle=\frac{1}{Z}e^{-\beta\hat{H}_{B}^{(g)}}\lvert e\rangle\langle
e\rvert$ (2.27) $\displaystyle Z$
$\displaystyle=\text{Tr}\left[e^{-\beta\hat{H}_{B}^{(g)}}\right]$ (2.28)

Notice here that the thermal distribution is based on the ground state
Hamiltonian; the intuition is that, prior to a fast photoexcitation, we expect
the nuclear modes of the system to be equilibrated in the ground electronic
state. When the fast photoexcitation occurs, the comparatively slow timescale
of the nuclear vibrations implies that the nuclei will not have time to
equilibrate according to the new electronic state; instead, dynamics will
launch from a state where the nuclei are still in the ground state geometry.

In the nonequilibrium case, the correlation function is significantly more
complicated and cannot be reduced to a function of a single time variable but
rather depends on two times:

$\displaystyle C(t^{\prime},t^{\prime\prime})$
$\displaystyle=\text{Tr}_{B}\left[\hat{\rho}_{neq}(0)\hat{V}_{eg}(t^{\prime})\hat{V}_{ge}(t^{\prime\prime})\right]$
(2.29)
$\displaystyle=\text{Tr}_{B}\left[e^{-\beta\hat{H}_{B}^{(g)}}e^{i\hat{H}_{B}^{(e)}t^{\prime}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(g)}(t^{\prime}-t^{\prime\prime})}e^{i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(e)}t^{\prime\prime}}\right]$
(2.30)
$\displaystyle=\frac{1}{Z}\sqrt{\frac{\det\left[\bm{a}_{g}(-i\beta)\bm{a}_{e}(-t^{\prime})\bm{a}_{g}(t^{\prime}-t^{\prime\prime})\bm{a}_{e}(t^{\prime\prime})\right]}{\det\Sigma}}$
$\displaystyle\qquad\times\exp\left(i\textbf{d}^{T}\left[\textbf{G}_{\beta}+\textbf{G}\right]\textbf{d}-\frac{1}{2}i\begin{bmatrix}\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}-\textbf{W}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}+\textbf{W}\end{bmatrix}^{T}\raisebox{-3.0pt}{\huge{$\Sigma$}}^{-1}\begin{bmatrix}\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}-\textbf{W}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}+\textbf{W}\end{bmatrix}\right)$ (2.31)
$\small{\raisebox{-3.0pt}{\huge{$\Sigma$}}=\begin{bmatrix}\textbf{S}^{T}\bm{b}_{g}(-i\beta)\textbf{S}+\bm{b}_{e}(t^{\prime\prime})&-\bm{a}_{g}(-i\beta)&0&-\bm{a}_{e}(t^{\prime\prime})\\\
-\bm{a}_{g}(-i\beta)&\textbf{S}^{T}\bm{b}_{g}(-i\beta)\textbf{S}+\bm{b}_{e}(-t^{\prime})&-\bm{a}_{e}(-t^{\prime})&0\\\
0&-\bm{a}_{e}(-t^{\prime})&\bm{b}_{e}(-t^{\prime})+\textbf{S}^{T}\bm{b}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}&-\textbf{S}^{T}\bm{a}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}\\\
-\bm{a}_{e}(t^{\prime\prime})&0&-\textbf{S}^{T}\bm{a}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}&\bm{b}_{e}(t^{\prime\prime})+\textbf{S}^{T}\bm{b}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}\end{bmatrix}}$
(2.32)

Again, the correlation function can be used to write down an expression for
the approximate ground-state population at time $t$:

$P_{g}(t)\simeq\lvert
V\rvert^{2}\int_{0}^{t}dt^{\prime}\,\int_{0}^{t}dt^{\prime\prime}\,e^{-i\Delta
G(t^{\prime\prime}-t^{\prime})}C(t^{\prime},t^{\prime\prime})$ (2.33)

###  2.3 An Intuitive Model with Two Intertwined Modes and Nuclear Friction

The above formulation is rather abstract, and it is therefore difficult to
properly interpret the obtained equations and their consequences. Therefore,
it will be helpful to define an intuitive test model consisting of two
intertwined nuclear modes (primary modes) each coupled to a large thermal bath
of secondary modes. The exact coupling to the bath of secondary modes will be
defined so as to produce Langevin dynamics for each primary nuclear mode. It
is easy to see how such a model maps onto a real condensed-phase chemical
system; one can imagine that the primary modes represent a molecule of
interest undergoing electron transfer dynamics while interacting with a
solvent environment, which is described by the secondary modes.
Mathematically, we describe this system using the following Hamiltonian:

$\displaystyle\hat{\mathcal{H}}=$
$\displaystyle\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\hat{H}_{B}^{(x)}+\hat{H}_{B}^{(y)}+$
(2.34)
$\displaystyle\begin{bmatrix}\frac{1}{2}\omega_{2}^{2}\left(\hat{x}-d\cos\theta\right)^{2}+\frac{1}{2}\omega_{1}^{2}\left(\hat{y}-d\sin\theta\right)^{2}+E_{g}&Ve^{iW(\hat{x}\cos\eta+\hat{y}\sin\eta)}\\\
V^{*}e^{-iW(\hat{x}\cos\eta+\hat{y}\sin\eta)}&\frac{1}{2}\omega_{1}^{2}(\hat{x}\cos\phi+\hat{y}\sin\phi)^{2}+\frac{1}{2}\omega_{2}^{2}(\hat{y}\cos\phi-\hat{x}\sin\phi)^{2}+E_{e}\end{bmatrix}$
$\displaystyle\hat{H}_{B}^{(x)}=\sum_{\alpha}\frac{1}{2}\omega_{\alpha}^{2}\left(\hat{q}_{\alpha}+\frac{c_{\alpha}\hat{x}}{\omega_{\alpha}^{2}}\right)^{2}\qquad\hat{H}_{B}^{(y)}=\sum_{\beta}\frac{1}{2}\omega_{\beta}^{2}\left(\hat{q}_{\beta}+\frac{c_{\beta}\hat{y}}{\omega_{\beta}^{2}}\right)^{2}$

Here, $\hat{x}$ and $\hat{y}$ are the position operators for the primary
modes, and are respectively coupled to the thermal baths given
by$\hat{H}_{B}^{(x)}$ and $\hat{H}_{B}^{(y)}$ (the exact couplings
$c_{\alpha}$ and $c_{\beta}$ will be defined below). These baths are
responsible for producing random fluctuations as well as an apparent friction
in the motion of the primary modes.

There are three main parameters defining the geometry of the primary-mode
system (i.e. the molecule of interest): $\phi$, $\theta$, and $\eta$. The
excited-state well (the lower right block of Eq. 2.34) is centered on the
origin and has a rotation angle defined by $\phi$ relative to the $x-$axis.
The minimum of the ground-state well (the upper left block of Eq. 2.34) is
shifted away from that of the excited-state minimum along an axis defined by
the angle $\theta$ relative to the $x-$axis. Together, $\phi$ and $\eta$
define the spatial geometry of the primary-mode system; the well geometries
are plotted for a few different values of $\phi$ and $\theta$ in Fig. 1.
Finally, $\eta$ describes the direction of the spin-orbit coupling vector.
These three parameters, together with the magnitudes of the shift and spin-
orbit coupling vectors ($d$ and $W$) are enough to specify the vector
parameters d, S, and W (respectively the shift vector, Duschinskii matrix, and
spin-orbit coupling vector) from Section 2.2.

We have imposed several important and deliberate constraints in the
construction of Eq. 2.34. First, note that the interstate coupling operator
depends purely on $\hat{x}$ and $\hat{y}$; it is wholly independent from the
secondary bath modes. This constraint reflects the physical notion that the
environment of a molecule has no influence on the spin-orbit coupling within
the molecule. Such an assumption, while likely not universally true in real
systems, is a reasonable approximation we have made for the sake of gaining
basic intuition about the problem at hand.

Second, notice that we have deliberately chosen to use the coordinate $\phi$
(which represents rotations of the excited state well) to characterize each
Hamiltonian. In general, note that rotating the excited state does not lead to
drastically different ranges of electron transfer rates; rotating the ground
state, however, would lead to drastically different ranges of electron
transfer rates. Therefore, in order to model dynamics over a range of
different geometries – but with a reasonable range of time scales – for a
fixed $\theta$, we will simply vary $\phi$ below and record the resulting
dynamics. More generally, even for a set of displaced multidimensional
harmonic oscillators with different frequencies, a reasonable energy scale can
be calculated the driving force that renders the electron transfer
barrierless; notice that this driving force depends only on the frequencies of
the acceptor molecule (and not on the frequencies of the donor molecule).
Thus, changing $\phi$ (while keeping $\theta$ fixed) is a reasonable means of
generating a comparable set of Hamiltonians with different geometries but
similar time scales of interest. See Fig. 2.

Now, for the Hamiltonian in Eq. 2.34, we have yet to define the precise form
of the constants $c_{\alpha}$ and $c_{\beta}$. For a Langevin model, we
require each of the two baths to have spectral densities of the form
$J(\omega)=\gamma\omega$, where $\gamma$ is the damping factor and defines the
strength of the apparent friction on the primary modes. If $\rho_{\alpha}$ and
$\rho_{\beta}$ are the densities of states in their respective baths, the
values of the $c_{\alpha}$ and $c_{\beta}$ terms are defined via the spectral
function as follows:

$\displaystyle
J_{\alpha}(\omega_{\alpha})=\frac{\pi}{2}\frac{c_{\alpha}^{2}}{\omega_{\alpha}}\rho_{\alpha}(\omega_{\alpha})=\gamma\omega_{\alpha}\quad\Rightarrow\quad
c_{\alpha}=\sqrt{\frac{2}{\pi}\frac{\gamma\omega_{\alpha}^{2}}{\rho_{\alpha}(\omega_{\alpha})}}$
(2.35) $\displaystyle
J_{\beta}(\omega_{\beta})=\frac{\pi}{2}\frac{c_{\beta}^{2}}{\omega_{\beta}}\rho_{\beta}(\omega_{\beta})=\gamma\omega_{\beta}\quad\Rightarrow\quad
c_{\beta}=\sqrt{\frac{2}{\pi}\frac{\gamma\omega_{\beta}^{2}}{\rho_{\beta}(\omega_{\beta})}}$
(2.36)

If we let $\bm{c}_{\alpha}$ and $\bm{c}_{\beta}$ be the vectors of couplings
$c_{\alpha}$ and $c_{\beta}$, and we let q be the vector containing all the
position operators of the system:

$\hat{\textbf{q}}=\begin{bmatrix}\hat{x}&\hat{y}&\\{\hat{q}_{\alpha}\\}&\\{\hat{q}_{\beta}\\}\end{bmatrix}^{T}$
(2.37)

then we can rewrite Hamiltonian in 2.34 in the form of Eq. 2.1:

$\mathcal{H}=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\begin{bmatrix}\frac{1}{2}\hat{\textbf{q}}^{T}\hat{\bm{\Omega}}^{2}_{g}\hat{\textbf{q}}+\bm{\lambda}_{g}^{T}\hat{\textbf{q}}+\tilde{E}_{g}&Ve^{i\textbf{W}^{T}\hat{\textbf{q}}}\\\
V^{*}e^{-i\textbf{W}^{T}\hat{\textbf{q}}}&\frac{1}{2}\hat{\textbf{q}}^{T}\hat{\bm{\Omega}}^{2}_{e}\hat{\textbf{q}}+\bm{\lambda}_{e}^{T}\hat{\textbf{q}}+E_{e}\end{bmatrix}$
(2.38)
$\hat{\bm{\Omega}}_{g}^{2}=\begin{bmatrix}\omega_{2}^{2}+\sum_{\alpha}\frac{c_{\alpha}^{2}}{\omega_{\alpha}^{2}}&0&\bm{c}_{\alpha}^{T}&\bm{0}\\\
0&\omega^{2}_{1}+\sum_{\beta}\frac{c_{\beta}^{2}}{\omega_{\beta}^{2}}&\bm{0}&\bm{c}^{T}_{\beta}\\\
\bm{c}_{\alpha}&\bm{0}&\bm{\Omega}_{\alpha}^{2}&\bm{0}\\\
\bm{0}&\bm{c}_{\beta}&\bm{0}&\bm{\Omega}_{\beta}^{2}\end{bmatrix}$
$\hat{\bm{\Omega}}_{e}^{2}=\begin{bmatrix}\omega_{1}^{2}\cos^{2}\phi+\omega_{2}^{2}\sin^{2}\phi+\sum_{\alpha}\frac{c_{\alpha}^{2}}{\omega_{\alpha}^{2}}&-\frac{1}{2}\sin(2\phi)(\omega_{2}^{2}-\omega_{1}^{2})&\bm{c}_{\alpha}^{T}&\bm{0}\\\
-\frac{1}{2}\sin(2\phi)(\omega_{2}^{2}-\omega_{1}^{2})&\omega^{2}_{2}\cos^{2}\phi+\omega_{1}^{2}\sin^{2}\phi+\sum_{\beta}\frac{c_{\beta}^{2}}{\omega_{\beta}^{2}}&\bm{0}&\bm{c}^{T}_{\beta}\\\
\bm{c}_{\alpha}&\bm{0}&{\bm{\Omega}}_{\alpha}^{2}&\bm{0}\\\
\bm{0}&\bm{c}_{\beta}&\bm{0}&{\bm{\Omega}}_{\beta}^{2}\end{bmatrix}$
$\bm{\lambda}_{e}=\begin{bmatrix}0&0&\bm{0}&\bm{0}\end{bmatrix}^{T}\qquad\bm{\lambda}_{g}=-d\begin{bmatrix}\omega_{2}^{2}\cos\theta&\omega_{1}^{2}\sin\theta&\bm{0}&\bm{0}\end{bmatrix}^{T}\qquad\textbf{W}=W\begin{bmatrix}\cos\eta&\sin\eta&\bm{0}&\bm{0}\end{bmatrix}^{T}$
$\tilde{E}_{g}=E_{g}+\frac{1}{2}\left(\omega_{1}^{2}\sin^{2}\theta+\omega_{2}^{2}\cos^{2}\theta\right)\cdot
d^{2}$

From here, we simply apply the procedure outlined in Section 2.1 to obtain the
ground- and excited-state normal mode matrices, the Duschinskii matrix S, and
the shift vector d, and then evaluate nonequilibrium Fermi golden rule rates.

In principle, the Hamiltonian that has been defined is infinite dimensional;
it is in this limit that we recover Langevin dynamics. In practice, however,
it is obviously impossible to simulate an infinite-dimensional system. A brief
note on the numerical implementation is provided in Appendix 7.3.

Figure 1: Shown is an example of a system that may be defined by the
Hamiltonian in Eq. 2.34. The parameters $d$, $W$, $\theta$, $\eta$, and $\phi$
are shown in context. $d$ and $\theta$ together define the shift between the
minima of the primary mode wells. Similarly, $W$ and $\eta$ together define
the spin-orbit coupling vector within the primary system. Finally, $\phi$
defines the rotation angle on the excited state well, as measured from the
major axis to the horizontal $x-$axis. Figure 2: Pictured here is a simple
one-dimensional two-state vibronic system. Notice that the
donor$\rightarrow$acceptor activation barrier is more sensitive to the
acceptor state frequency than it is to the donor state frequency. This is
illustrated by the observation that the changes in energy from $\bm{a}$ to
$\bm{c}$ and $\bm{b}$ to $\bm{d}$ are not as significant as those from
$\bm{a}$ to $\bm{b}$ and $\bm{c}$ to $\bm{d}$. By analogy, we expect that the
activation barrier for a much larger multidimensional will be more sensitive
to the geometry of the acceptor state well than that of the donor state well.
Therefore, to ensure a similar timescale for electron transfer across a range
of Hamiltonians, we will use the parameter $\phi$ to rotate the excited state
well rather than the ground state well.

##  3 Results

The model above is characterized by quite a few parameters. A summary of what
the various model parameters mean in the context of the system being studied
is provided in Table 1 for the reader’s convenience.

Parameter | Physical meaning
---|---
$\gamma$ | The damping factor which defines the friction in the system
$V$ | The magnitude of the interstate coupling
$\phi$ | The rotation angle of the excited-state well orientation
$\theta$ | The rotation angle of the primary-mode shift vector
$d$ | The magnitude of the primary-mode shift vector
$\eta$ | The rotation angle of the spin-orbit vector
$W$ | The magnitude of the spin-orbit coupling vector
Table 1: A summary of the parameters that define the system at hand. The first
three terms define the strengths of the various couplings in the system; the
latter four define the physical geometry.

All of the results below use the following parameters (atomic units):

$\omega_{1}=2\cdot 10^{-4}$ $\omega_{2}=4\cdot 10^{-4}$ $\gamma=4\cdot
10^{-4}$

$W=0.05$ $\beta=1000$ $V=1\cdot 10^{-4}$

There are three different values of $\theta$ investigated: $\pi/4$, $\pi/6$,
and 0. For the first case, $d$ is set as 884 au; for the latter two cases, $d$
is set as 769 au. (The value for $d$ has been adjusted for the latter two
cases in order to keep the activation energy approximately constant as
$\theta$ is changed.) Finally, three different values of $\Delta G$ are
probed: 0, -0.01, and -0.029. Note that $\Delta G$ is not the bias between the
primary mode wells, but rather the energy gap between the well minima after
coupling to the bath is accounted for. In other words, $\Delta G$ is as
defined in Eq. 2.10 (and is not the same as $E_{g}$ in Eq. 2.34 or
$\tilde{E}_{g}$ in Eq. 2.38). The rationale behind these choices for $\Delta
G$ will be explained below.

###  3.1 Equilibrium Dynamics: System Characterization

Figure 3: Equilibrium transfer rate plotted against $\Delta G$ for the model
Hamiltonian in Eq. 2.34 with $\theta=\pi/4$. Solid curves are for $W=0$ and
dashed curves are for $W=0.05$ with $\eta=\pi/2$; in practice, the value of
$\eta$ does not lead to to any significant effect here. Note also that
changing the rotation angle $\phi$ does not significantly alter the $\Delta G$
at which the maximal transfer rate is achieved (i.e. the rate-defining
reorganization energy of the system).

We begin with a brief look at the equilibrium behavior of the model as a means
of characterizing the system. This characterization will be performed based on
Marcus-style rate curves, where the equilibrium transfer rate is plotted as a
function of $\Delta G$. Several of these curves for different values of $\phi$
are plotted in Fig. 3. Here, the well geometry used is $\theta=\pi/4$; the
corresponding plots for the other 2 geometries are nearly identical. From
these curves, we can extract information about the standard and geometric
reorganization energies of the system.

First, note that the system has the desired property that changing the
rotation angle $\phi$ does not significantly alter the electron transfer rates
(which are on a similar order of magnitude).

Next, notice that for all rotation angles $\phi$, adding spin-orbit coupling
to the system (making $W$ nonzero) slightly increases the activation barrier,
shifting the Marcus curve to the left with a magnitude of about 0.0012 Hartree
for all $\phi$. By analogy with standard Marcus theory, one might then argue
that the “effective” reorgnization energy has therefore experienced a small
shift, on the order of 3-4 per cent. This behavior is consistent with the
behavior of the simpler shifted harmonic oscillator model seen in Ref. [29]
where one could easily define a total reorganization (without the headache of
having two different reorganization energies from donor-to-acceptor and
acceptor-to-donor).

###  3.2 Photoexcitation Dynamics

Before presenting the main results pertaining to spin polarization, let us
establish a simple picture of the relaxation dynamics in our model. The ground
state population curve can be loosely characterized by three epochs, which can
be observed in the dynamics for any parameter set. First, there is a brief
epoch of no population growth. Then, there is a large spike in the growth
rate. Finally, the rate spike decays leaving a steady, constant growth of the
ground state population. As an example, see Fig. 6. These three epochs can be
understood semiclassically as follows: First, following the photoexcitation,
the nuclei are quite far from the crossing seam between the excited and ground
state diabats. Therefore, as they begin to evolve on the excited state
surface, not much population transfer can take place. However, the nuclei soon
approach the crossing seam. Here, population is readily transferred from the
excited state to the ground state. However, all the while, energy is being
dissipated by the large bath coupled to the primary modes. This dissipation
causes the system to decohere, ultimately resulting in a decay to equilibrium
dynamics, which accounts for the steady and constant growth seen in the third
epoch of the dynamics. For our model, we find that 25,000 a.u. is sufficient
time for the system to fully decohere; the absolute difference in ground state
up- and down-spin populations will remain constant after this time. Because we
are primarily interested in the persistence of spin-polarization, we will
perform all polarization analysis at this time, after decoherence has
completed.

Figure 4: Plotted are the final ground-state spin polarizations ($\chi$ in Eq.
3.1) as functions of $\phi$ and $\eta$ for three different values of $\theta$
and three different values of $\Delta G$. The three columns from left to right
correspond respectively to $\Delta G=0$, $\Delta G=-0.01$, and $\Delta
G=-0.029$. The three rows from top to bottom respectively correspond to
$\theta=\pi/4$, $\theta=\pi/6$, and $\theta=0$ (with respective values of
$d=884$, $d=769$, and $d=769$). Also plotted in solid cyan are the zero-
polarization isolines. Note that the polarization can depend quite sensitively
on the values of $\phi$ and $\eta$ (i.e. the geometry of the system).

In Figs. 4 and 5, we plot the polarization and the ground state populations at
time 25,000 au. If we denote the ground-state populations for spin-up and
spin-down as $P^{+}$ and $P^{-}$, the polarization is defined as:

$\chi=\frac{P^{+}-P^{-}}{P^{+}+P^{-}}$ (3.1)

Similarly, the total ground state population is defined follows:

$P_{g}=\frac{P^{+}+P^{-}}{2}$ (3.2)

Note that we expect a realistic initial state to be completely unpolarized,
i.e. the initial electronic populations should be equally split between spin-
up and spin-down. Therefore within our model, because we ignore spin-flip
processes, one can calculate the total overall population dynamics (with Eq.
3.2) by taking a simple average of the spin-up and spin-down populations.

For both of these data sets (Fig. 4 and 5), dynamics were computed with
different parameter pairs ($\phi$, $\eta$), different well positions ($d$,
$\theta$), and a $\Delta G$ of either 0, -0.01 or -0.029. These three regimes
for $\Delta G$ were chosen based on the Marcus curves in Fig. 3. On the one
hand, -0.01 Hartree is a modest $\Delta G$ for which there is a nonzero donor-
acceptor driving force but still a large activation barrier; the equilibrium
transfer rate is quite sensitive to $\phi$. On the other hand, -0.029 is a
value for $\Delta G$ at which all the rate curves for different values of
$\phi$ nearly intersect, suggesting similarity in activation barrier heights.
Finally, at $\Delta G=0$, the system is exactly symmetric in ground- and
excited-state well heights.

Figure 5: Plotted are the final ground-state populations $P_{g}$ as functions
of $\phi$ and $\eta$ for three different values of $\theta$ and three
different values of $\Delta G$. Each panel uses the same parameter set as the
corresponding panel in Fig. 4. Note that the population can vary greatly with
$\phi$, but shows very little dependence on $\eta$.

Note several general features of the data in Fig. 4. First, note that the data
is $\pi$-periodic in $\phi$ and $\pi$-antiperiodic in $\eta$. The periodicity
in $\phi$ is easily explained by the fact that the well geometry for
$\phi=\phi_{0}$ is exactly the same as the well geometry for
$\phi=\phi_{0}+n\pi$. The $\eta$ antiperiodicity can be understood from the
simple observation that:

$W(\hat{x}\cos\eta+\hat{y}\sin\eta)=-W(\hat{x}\cos(\eta+\pi)+\hat{y}\sin(\eta+\pi))$
(3.3)

Second, by comparison with Fig. 5, note that $\eta$ (which, as a reminder,
defines the direction of the spin-orbit coupling vector) can dramatically
change the spin polarization of the final state without significantly altering
the final population $P_{g}$. This finding indicates that the observed spin
polarization effects do not necessarily correlate with the final state
populations. This clear distinction between polarization and population is
highlighted in Fig. 6, where we plot dynamics for two different pairs of
systems, where the only difference within each pair is the direction of the W
vector.

Finally, in Fig. 7, we investigate the temperature-dependence of the spin-
polarization effects observed above. We find that the well rotation angle
$\phi$ is tightly coupled to the low-temperature behavior of the spin
polarization. As an example, notice in Fig. 7 that for the shifted harmonic
oscillator case $\phi=\pi/2$, the polarization attains a maximum at around 400
K, whereas adding a Duschinskii rotation and switching to the $\phi=0$ case
results in a completely monotonic dependence on temperature.

Figure 6: Plotted are ground state populations (top) and spin polarizations
(bottom) following a fast photoexcitation for two pairs of systems. On the
left are two SHO systems ($\phi=\pi/2$), and on the right are two systems
which have a 90 degree well rotation ($\phi=0$). In each panel, one of the
systems has the spin-orbit coupling vector W directed along the x-axis
($\eta=0$) and the other has has W directed along the y-axis
($\eta=\frac{\pi}{2}$). Note that although the final population averaged over
the up- and down-spin dynamics is approximately the same for the two values of
$\eta$ (i.e. the average of the dashed red and dashed yellow lines is near the
average of the solid red and solid yellow lines), the spin polarization for
$\eta=0$ is very different from the spin polarization for $\eta=\pi/2$
(especially in the case $\phi=0$). Figure 7: Plotted on the left is the
magnitude of the final ground state spin-polarization as a function of $\phi$
and temperature for the configuration $\theta=\pi/4$, $\eta=\pi/2$, $\Delta
G=-0.01$. Note that there are very specific values of $\phi$ which promote
spin polarization at low temperatures. In particular, when $\phi=0$, spin-
polarization increases monotonically as temperature tends to 0. This trend is
in stark contrast to the $\phi=\pi/2$ (shifted harmonic oscillator) case, in
which the spin-polarization attains a maximum around 400 K.

##  4 Discussion

Figure 8: The dynamics in the a panels are the same as those in the
corresponding b and c panels. This is due to the Hamiltonian in Eq. 2.38 being
invariant under reflections and rotations in position space. This allows us to
easily connect chirality to the emergence of spin-dependent transfer effects.
On the one hand, the dynamics of system $\bm{1}\textbf{a}$ will have no spin-
dependence; system $\bm{1}\textbf{c}$ is precisely the $-\textbf{W}$ analog of
system $\bm{1}\textbf{a}$. On the other hand, the dynamics of systems
$\bm{2}\textbf{a}$ and $\bm{3}\textbf{a}$ are not expected to exhibit spin-
independence, since systems $\bm{2}\textbf{c}$ and $\bm{3}\textbf{c}$ are not
superimposable on the $-\textbf{W}$ analogs of $\bm{2}\textbf{a}$ and
$\bm{3}\textbf{a}$. Although system $\bm{2}\textbf{a}$ is spatially achiral,
it possesses asymmetry with respect to W. System $\bm{3}\textbf{a}$ has
spatial chirality, so any spatial reflection will lead to a nonsuperimposable
system.

Overall, the observations above allow us to readily answer both of the two
overarching questions posed in the introduction. First, the direction of the
spin-orbit coupling vector is extremely important in generating spin
polarization in electron transfer, and the optimal direction is tightly
coupled to the geometry of the system. Second, changes in normal mode
frequencies and couplings upon electron transfer can greatly promote spin
polarization over shifted harmonic oscillator systems, and these effects can
be amplified significantly at low temperatures. As can be seen in Fig. 4,
significantly more spin polarization can be observed in systems with a 90
degree well rotation ($\phi=0$) over the SHO case ($\phi=\frac{\pi}{2}$).
Additionally, we see that $\phi=0$ systems more consistently produce sizable
spin polarization than $\phi=\frac{\pi}{2}$ systems. Although we have been
able to produce a satisfying answer to the two overarching questions of this
investigation, we would also ultimately like some intuition as to the general
shape of the polarization landscape. Despite the evident complexity, as well
as the nontrivial relationship the landscape has with $\Delta G$, we will put
forward a hypothesis that at least partially explains the observed results.

Consider the various systems depicted in Fig. 8. The electronic dynamics of
the systems in the a panels are exactly the same as the dynamics in the
corresponding b and c systems. In order to see this equivalence, note that the
b systems are obtained by mirror reflecting the left-hand systems in space,
and that the c panels are subsequently obtained by a simple rotation of the b
systems. The electronic dynamics of the Hamiltonian in 2.38 are invariant to
such transformations, which are equivalent to basis transformations of the
nuclear coordinates. This relationship demonstrates the natural connection
between the spin-dependent effects observed above and system chirality. We can
easily see that system $\bm{1}\textbf{a}$ in Fig. 8 must have exactly the same
dynamics for $\pm\textbf{W}$, since system $\bm{1}\textbf{c}$ is exactly
superimposable onto its $-\textbf{W}$ analog. A chiral system (such as that in
panel $\bm{3}\textbf{a}$) lacks spatial symmetry, and will typically exhibit
dynamical differences for different spin states. In such a system, it is not
possible to map $+\textbf{W}$ dynamics to $-\textbf{W}$ dynamics through
operations which respect the symmetry of the Hamiltonian (i.e. spatial
reflections and rotations).

However, in the systems presently under consideration, spatial symmetry is not
enough to eliminate spin-dependent transfer effects; consider system
$\bm{2}\textbf{a}$ in Fig. 8. Although the potential wells are superimposable
upon one another following a rotation, we find that the W vector of
$\bm{2}\textbf{c}$ does not map to $-\textbf{W}$ in $\bm{2}\textbf{a}$.
Consequently, there is no reason to expect that the dynamics of system
$\bm{2}\textbf{a}$ will be spin-independent. Therefore, for the purposes of
understanding spin-dependent dynamics, we must adopt a much more restrictive
notion of achirality than is standard; for our purposes, an achiral system
must satisy two constraints. First, as usual, $(i)$ there must be a mirror
plane; second, not as usual, $(ii)$ that mirror plane must be orthogonal to W.
Only in such a case, when the reflection maps W to $-\textbf{W}$, will we be
left with a system that is spatially superimposable onto the original system.

Now we return to the systems investigated in Fig. 4. We expect minimum
polarizations to occur at some system configurations where chirality (as per
the above definition) is “minimized”; in other words, the system is as close
to achiral as possible. First, we focus on the case $\theta=\pi/4$. In the
spirit of condition $(i)$, we will first identify values of $\phi$ which
produce well configurations with relatively high spatial symmetry.
$\phi=\frac{\pi}{4}$ and $\phi=-\frac{\pi}{4}$ are two such configurations; in
each of these two systems, the major or minor axis of the excited state
surface is oriented along the axis connecting the well minima (see Fig. 9,
where the relevant configurations are pictured in the lower panels). Next,
following condition $(ii)$, we expect that the minimum chirality configuration
will be achieved with $\eta=\frac{3\pi}{4}$, because this condition places the
mirror plane on the axis connecting the well minima, making the excited state
surface exactly symmetric (again, see Fig. 9). These intuitive arguments are
indeed confirmed by the numerical results in Fig. 4. In Fig. 9, we take the
zero-polarization isolines from Fig. 4 and plot them together for each value
of $\theta$ (i.e. each variation of well position). We see that, even as
$\Delta G$ changes, the $\theta=\pi/4$ zero-polarization level curves pass
nearly exactly through both of the points
$(\phi,\eta)=\left(\frac{\pi}{4},\frac{3\pi}{4}\right),\left(-\frac{\pi}{4},\frac{3\pi}{4}\right)$
(marked in Fig. 9 with black circles). Similarly, the level curves for the
case $\theta=\pi/6$ all nearly pass through the points
$(\phi,\eta)=\left(\frac{\pi}{6},\frac{2\pi}{3}\right),\left(-\frac{\pi}{3},\frac{2\pi}{3}\right)$
and those for the case $\theta=0$ nearly pass through
$(\phi,\eta)=\left(0,\frac{\pi}{2}\right),\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$.

Interestingly, the notion of achirality developed above can also be applied
(roughly) to understand the locations of the “islands” of strong polarization
in Fig. 4. Consider again, for illustrative purposes, the case $\theta=\pi/4$.
If minimum chirality configurations occur at $\phi=\pm\pi/4$, then we would
naturally expect $\phi=0,\pi/2$ to produce maximum chirality configurations
and have the potential for producing large spin polarization. This is again
confirmed numerically in Fig. 4, where we see strong polarization islands
around $\phi=0$ and $\phi=\frac{\pi}{2}$. Similarly, for $\theta=\pi/6$ and
$\theta=0$ respectively, we observe islands of strong polarization around
$\phi=\frac{5\pi}{12},-\frac{\pi}{12}$ and $\phi=\pm\frac{\pi}{4}$. All of
these observations taken together suggest that the strong notion of chirality
developed above may partly account for the polarization landscape in Fig. 4.
However, if chirality were the entire story, we would expect the zero-
polarization level curves to be independent of $\Delta G$, which is clearly
not the case. One must suppose that, for different values of $\Delta G$, the
dynamics become more complicated as one cannot ignore the time scale for
energy dissipation when calculating transient spin polarization.

Figure 9: (Upper) The three zero-polarization isolines for each well
positioning are plotted together. Blue, yellow, and red traces correspond
respectively to $\Delta G$ values of 0, -0.01, and -0.029. Black circles mark
the points $(\phi,\eta)=\left(\pm\frac{\pi}{4},\frac{3\pi}{4}\right)$ on the
left,
$(\phi,\eta)=\left(\frac{\pi}{6},\frac{2\pi}{3}\right),\left(-\frac{\pi}{3},\frac{2\pi}{3}\right)$
in the middle, and
$(\phi,\eta)=\left(0,\frac{\pi}{2}\right),\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$
on the right. Note that, in every panel, all three isolines intersect very
close to these target points which designate expected minimal chirality
configurations. (Lower) Pictured for each value of $\theta$ are the two
minimum chirality configurations designated by the black circles. In each
panel, the ground state (purple) minimum is indicated with an “x” and the
relevant reflection plane (i.e. that plane orthogonal to W which bisects the
excited state well) is indicated with a red dashed line.

As a final remark, it is interesting to note that Duschinskii rotations, in
addition to promoting spin polarization, also appear to produce greater
temporal separation between up- and down-spin dynamics (i.e. greater transient
spin polarization). In other words, as can be seen in Fig. 6, the differences
in populations for the up- and down-spin dynamics in the $\phi=0$ case are
larger than they are for the $\phi=\pi/2$ (SHO) case; in fact, for both the
$\eta=0$ and $\eta=\pi/2$ cases, the spin up and spin down populations largely
overlap for the $\phi=\pi/2$ case. Notably, as $\eta$ varies, Fig. 6 shows
that greater transient polarization at early times in the SHO model correlates
with decreased final spin-polarization, whereas greater transient polarization
at early times in the rotated-well model correlates with increased final spin
polarization. Overall, this result suggests that if one can drive
nonequilibrium electron transfer dynamics with Duschinskii rotations so as to
prevent equilibration, one might recover even larger spin polarizations than
what we have found here. Future work may well investigate the effect of such a
driving force.

##  5 Conclusions

In conclusion, we have used a Duschinskii rotation formalism to investigate
spin-polarization in a Langevin-ESOC model. Our investigation demonstrates
that the final spin polarization can be quite large and highly sensitive to
the exact system geometry. This finding has particular relevance for systems
where changes in electronic state may be accompanied by significant changes in
coupling between vibrational modes. A notion of chirality has been shown to
partially explain how polarization can depend on system geometry, but further
work is needed to understand and predict such polarization in general. While
the model here is very abstract (and we have not yet used ab initio
calculations to extract parameters for a given molecule), our conclusions
highlight that standard Marcus models for electron transfer using only the
donor and acceptor minimum geometries may be very inadequate as far as
reproducing spin polarization of a realistic system. Future extensions of this
work would also benefit by treating multi-electron problems. For instance,
although a great deal of work has been done treating the singlet to triplet
ISC problem, [33, 34], the spin polarization of product states has not yet
been fully explored[21, 35]. In the future, it will also be very helpful to
extend the two-level results presented here to the case of a metal surface
with many electronic states, where the CISS phenomenon is often observed.[36]

##  6 Acknowledgements

This work was supported by the U.S. Air Force Office of Scientific Research
(USAFOSR), AFOSR Grant Nos. FA9550-18-1-0 497 and FA9550-18-1-0 420.

##  7 Appendix

In this appendix we provide details on the derivation of the equilibrium and
nonequilibrium time correlation functions for a harmonic system incorporating
Duschinskii rotations. We will make extensive use of two identities:

$\displaystyle\int_{\mathbb{R}^{n}}dx^{n}\,\exp\left(-\frac{1}{2}\textbf{x}^{\text{T}}\hat{A}\textbf{x}+B^{\text{T}}\textbf{x}\right)=\sqrt{\frac{(2\pi)^{n}}{\text{det}A}}\exp\left(\frac{1}{2}B^{T}\hat{A}^{-1}B\right)$
(7.1)
$\displaystyle\langle\textbf{x}|e^{-i\hat{H}t}|\textbf{y}\rangle=\sqrt{\frac{\det{\bm{a}(t)}}{(2\pi
i)^{N}}}\exp\left(\frac{i}{2}\left[\bm{x}^{T}\bm{b}(t)\bm{x}+\bm{y}^{T}\bm{b}(t)\bm{y}-2\bm{x}^{T}\bm{a}(t)\bm{y}\right]\right)$
(7.2)

The former identity is simply a multivariate Gaussian integral. The latter is
derived from the explicit evaluation of the Green’s function for a harmonic
Hamiltonian $\hat{H}$ composed of $N$ independent harmonic modes. Here,
$\bm{a}(t)=\hat{\bm{\Omega}}[\sin(\hat{\bm{\Omega}}t)]^{-1}$ and
$\bm{b}(t)=\hat{\bm{\Omega}}[\tan(\hat{\bm{\Omega}}t)]^{-1}$, where
$\hat{\bm{\Omega}}$ is the (diagonal) harmonic matrix of $\hat{H}$.

The Hamiltonian of interest is given in Eq. 2.9. We will rewrite it here for
convenience:

$\hat{H}=\left(E_{g}+\hat{H}_{B}^{(g)}\right)\lvert g\rangle\langle
g\rvert+\left(E_{e}+\hat{H}_{B}^{(e)}\right)\lvert e\rangle\langle
e\rvert+Ve^{i\textbf{W}^{T}\hat{\textbf{x}}}\lvert g\rangle\langle
e\rvert+V^{*}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}\lvert e\rangle\langle
g\rvert$ (7.3)
$\hat{H}_{B}^{(e)}=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\frac{1}{2}\hat{\textbf{x}}^{\dagger}\hat{\bm{\Omega}_{e}}^{2}\hat{\textbf{x}}$
$\hat{H}_{B}^{(g)}=\frac{\hat{\textbf{p}}\cdot\hat{\textbf{p}}}{2}+\frac{1}{2}\left(\textbf{S}\hat{\textbf{x}}+\textbf{d}\right)^{\dagger}\hat{\bm{\Omega}_{g}}^{2}\left(\textbf{S}\hat{\textbf{x}}+\textbf{d}\right)$

###  7.1 Derivation of Eq. 2.25

We start the system in the excited state with the density matrix:

$\displaystyle\hat{\rho}(0)$
$\displaystyle=\frac{1}{Z}\exp(-\beta\hat{H}_{B}^{(e)})\lvert e\rangle\langle
e\rvert$ (7.4) $\displaystyle Z$
$\displaystyle=\text{Tr}\left[\exp(-\beta\hat{H}_{B}^{(e)})\right]$ (7.5)

We want to approximate:

$\displaystyle P_{g}(t)$
$\displaystyle=\frac{1}{Z}\text{Tr}_{B}\left[\left\langle
g\middle|e^{-i\hat{H}t}e^{-\beta\hat{H}_{B}^{(e)}}\middle|e\right\rangle\left\langle
e\middle|e^{i\hat{H}t}\middle|g\right\rangle\right]$ (7.6)
$\displaystyle=\frac{1}{Z}\text{Tr}_{B}\left[\left\langle
g\middle|e^{-i\hat{H}t}\middle|e\right\rangle
e^{-\beta\hat{H}_{B}^{(e)}}\left\langle
e\middle|e^{i\hat{H}t}\middle|g\right\rangle\right]$ (7.7)

We make the following definition for convenience:

$\displaystyle\hat{H}_{0}$
$\displaystyle=\left(E_{g}+\hat{H}_{B}^{(g)}\right)\lvert g\rangle\langle
g\rvert+\left(E_{e}+\hat{H}_{B}^{(e)}\right)\lvert e\rangle\langle e\rvert$
(7.8) $\displaystyle\hat{V}$
$\displaystyle=Ve^{i\textbf{W}^{T}\hat{\textbf{x}}}\lvert g\rangle\langle
e\rvert+V^{*}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}\lvert e\rangle\langle
g\rvert$ (7.9)

and expand $e^{-i\hat{H}t}$ in a Dyson series:

$\displaystyle e^{-i\hat{H}t}$ $\displaystyle\simeq
e^{-i\hat{H}_{0}t}\left(\textbf{I}-i\int_{0}^{t}dt^{\prime\prime}\,e^{i\hat{H}_{0}t^{\prime\prime}}\hat{V}e^{-i\hat{H}_{0}t^{\prime\prime}}\right)$
(7.10) $\displaystyle\Rightarrow\left\langle
g\middle|e^{-i\hat{H}t}\middle|e\right\rangle$ $\displaystyle\simeq-
iVe^{-iE_{g}t}\int_{0}^{t}dt^{\prime\prime}\,e^{i(E_{g}-E_{e})t^{\prime\prime}}e^{i\hat{H}_{B}^{(g)}t^{\prime\prime}}e^{i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(e)}t^{\prime\prime}}$
(7.11) $\displaystyle\Rightarrow\left\langle
e\middle|e^{i\hat{H}t}\middle|g\right\rangle$ $\displaystyle\simeq
iVe^{iE_{g}t}\int_{0}^{t}dt^{\prime}\,e^{-i(E_{g}-E_{e})t^{\prime}}e^{-i\hat{H}_{B}^{(g)}t^{\prime}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{i\hat{H}_{B}^{(e)}t^{\prime}}$
(7.12)

$\displaystyle\Rightarrow P_{g}(t)$ $\displaystyle\simeq\frac{\lvert
V\rvert^{2}}{Z}\int_{0}^{t}dt^{\prime}\,\int_{0}^{t}dt^{\prime\prime}\,e^{-i(E_{g}-E_{e})(t^{\prime}-t^{\prime\prime})}\text{Tr}_{B}\left[e^{-(\beta-i(t^{\prime}-t^{\prime\prime}))\hat{H}_{B}^{(e)}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(g)}(t^{\prime}-t^{\prime\prime})}e^{i\textbf{W}^{T}\hat{\textbf{x}}}\right]$
(7.13) $\displaystyle=\frac{2\lvert
V\rvert^{2}}{Z}\text{Re}\left[\int_{0}^{t}dt^{\prime}\,\int_{0}^{t^{\prime}}d\tau\,e^{-i(E_{g}-E_{e})\tau}\text{Tr}_{B}\left[e^{-(\beta-i\tau)\hat{H}_{B}^{(e)}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(g)}\tau}e^{i\textbf{W}^{T}\hat{\textbf{x}}}\right]\right]$
(7.14)

Letting $\tau_{e}=-\tau-i\beta$ and $\tau_{g}=\tau$, we need to compute the
correlation function

$C(\tau)=\text{Tr}_{B}\left[e^{-i\hat{H}_{B}^{(e)}\tau_{e}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(g)}\tau_{g}}e^{i\textbf{W}^{T}\hat{\textbf{x}}}\right]$
(7.15)

We proceed by evaluating the trace using position coordinates:

$\displaystyle C(\tau)$
$\displaystyle=\int_{\mathbb{R}^{n}}dx^{n}\,\left\langle\textbf{x}\middle|e^{-i\tau_{e}\hat{H}_{B}^{(e)}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(g)}\tau_{g}}e^{i\textbf{W}^{T}\hat{\textbf{x}}}\middle|\textbf{x}\right\rangle$
(7.16)
$\displaystyle=\int_{\mathbb{R}^{n}}dx^{n}\,\int_{\mathbb{R}^{n}}dy^{n}\,\int_{\mathbb{R}^{n}}dx^{\prime
n}\,\int_{\mathbb{R}^{n}}dy^{\prime
n}\,\left\langle\textbf{x}\middle|e^{-i\tau_{e}\hat{H}_{B}^{(e)}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}\middle|\textbf{y}\right\rangle\left\langle\textbf{y}\middle|\textbf{y}^{\prime}\right\rangle\left\langle\textbf{y}^{\prime}\middle|e^{-i\hat{H}_{B}^{(g)}\tau_{g}}\middle|\textbf{x}^{\prime}\right\rangle\left\langle\textbf{x}^{\prime}\middle|e^{i\textbf{W}^{T}\hat{\textbf{x}}}\middle|\textbf{x}\right\rangle$
(7.17)

$\textbf{x}^{\prime}$ and $\textbf{y}^{\prime}$ represent ground state
coordinates, while x and y represent excited state position coordinates.
Therefore, we have that
$\langle\textbf{y}|\textbf{y}^{\prime}\rangle=\delta(\bm{y}^{\prime}-(\textbf{S}\bm{y}+\textbf{d}))$.
Using this fact and Eq. 7.2, the integral reduces to

$\displaystyle\sqrt{\frac{\det(\bm{a}_{e})}{(2\pi
i)^{n}}}\sqrt{\frac{\det(\bm{a}_{g})}{(2\pi
i)^{n}}}\int_{\mathbb{R}^{n}}dx^{n}\,\int_{\mathbb{R}^{n}}dy^{n}\,e^{i\left[\frac{1}{2}\left(\bm{x}^{\text{T}}\bm{b}_{e}\bm{x}+\bm{y}^{\text{T}}\bm{b}_{e}\bm{y}\right)-\bm{x}^{\text{T}}\bm{a}_{e}\bm{y}\right]}\times$
(7.18) $\displaystyle
e^{i\left[\frac{1}{2}\left(\left(\textbf{S}\bm{x}+\textbf{d}\right)^{\text{T}}\bm{b}_{g}\left(\textbf{S}\bm{x}+\textbf{d}\right)+\left(\textbf{S}\bm{y}+\textbf{d}\right)^{\text{T}}\bm{b}_{g}\left(\textbf{S}\bm{y}+\textbf{d}\right)\right)-\left(\textbf{S}\bm{x}+\textbf{d}\right)^{\text{T}}\bm{a}_{g}\left(\textbf{S}\bm{y}+\textbf{d}\right)\right]}e^{i\textbf{W}^{T}\bm{x}}e^{-i\textbf{W}^{T}\bm{y}}$

Here, $\bm{a}_{e}$ and $\bm{a}_{g}$ are shorthand for $\bm{a}_{e}(\tau_{e})$
and $\bm{a}_{g}(\tau_{g})$, and these functions are as in Sec. 2.2. Next, we
focus on the integrand. We would like to collect terms and put it in the form
of a multivariate Gaussian. Recall the following operators from Sec. 2.2
(again using the shorthand in Eq. 7.18):

A $\displaystyle=\bm{a}_{e}+\textbf{S}^{T}\bm{a}_{g}\textbf{S}$ (7.19) B
$\displaystyle=\bm{b}_{e}+\textbf{S}^{T}\bm{b}_{g}\textbf{S}$ (7.20) G
$\displaystyle=\bm{b}_{g}-\bm{a}_{g}$ (7.21) E
$\displaystyle=\bm{b}_{e}-\bm{a}_{e}$ (7.22)

We can use these operators to write the integrand as

$\displaystyle\exp\left(i\left[\frac{1}{2}\left(\bm{x}^{T}\textbf{B}\bm{x}+\bm{y}^{T}\textbf{B}\bm{y}\right)-\bm{x}^{T}\textbf{A}\bm{y}+\left(\textbf{d}^{T}\textbf{GS}+\textbf{W}^{T}\right)\bm{x}+\left(\textbf{d}^{T}\textbf{GS}-\textbf{W}^{T}\right)\bm{y}+\textbf{d}^{T}\textbf{G}\textbf{d}\right]\right)$
(7.24)
$\displaystyle=\exp\left(i\left[\frac{1}{2}\begin{bmatrix}\bm{x}^{T}&\bm{y}^{T}\end{bmatrix}\begin{bmatrix}\textbf{B}&-\textbf{A}\\\
-\textbf{A}&\textbf{B}\end{bmatrix}\begin{bmatrix}\bm{x}\\\ \bm{y}\\\
\end{bmatrix}+\begin{bmatrix}\textbf{d}^{T}\textbf{GS}+\textbf{W}^{T}&\textbf{d}^{T}\textbf{GS}-\textbf{W}^{T}\end{bmatrix}\begin{bmatrix}\bm{x}\\\
\bm{y}\\\ \end{bmatrix}+\textbf{d}^{T}\textbf{G}\textbf{d}\right]\right)$
(7.25)

Thus, the double integral above becomes a single integral over
$\mathbb{R}^{2n}$:

$\int_{\mathbb{R}^{2n}}dz^{2n}\,\exp\left(-\frac{1}{2}\left(\bm{z}^{T}\begin{bmatrix}-i\textbf{B}&i\textbf{A}\\\
i\textbf{A}&-i\textbf{B}\end{bmatrix}\bm{z}\right)+i\begin{bmatrix}\textbf{d}^{T}\textbf{GS}+\textbf{W}^{T}&\textbf{d}^{T}\textbf{GS}-\textbf{W}^{T}\end{bmatrix}\bm{z}+i\textbf{d}^{T}\textbf{G}\textbf{d}\right)$
(7.26)

where ${\bm{z}}=\left[\begin{array}[]{c}{\bm{x}}\\\
{\bm{y}}\end{array}\right].$

Using the multivariate Gaussian integral identity, this evaluates to

$\sqrt{\frac{(2\pi)^{2n}}{(-i)^{2n}\det\begin{bmatrix}\textbf{B}&-\textbf{A}\\\
-\textbf{A}&\textbf{B}\end{bmatrix}}}\exp\left(i\textbf{d}^{T}\textbf{G}\textbf{d}-\frac{1}{2}i\begin{bmatrix}\textbf{d}^{T}\textbf{GS}+\textbf{W}^{T}&\textbf{d}^{T}\textbf{GS}-\textbf{W}^{T}\end{bmatrix}\begin{bmatrix}\textbf{B}&-\textbf{A}\\\
-\textbf{A}&\textbf{B}\end{bmatrix}^{-1}\begin{bmatrix}\textbf{S}^{T}\textbf{G}^{T}\textbf{d}+\textbf{W}\\\
\textbf{S}^{T}\textbf{G}^{T}\textbf{d}-\textbf{W}\end{bmatrix}\right)$ (7.27)

First, the determinant in the square root can be expressed as
$\det(\textbf{B})\det(\textbf{B}-\textbf{A}\textbf{B}^{-1}\textbf{A})$.
Second, we can explicitly write out the inverse matrix in the exponential as
follows:

$\begin{bmatrix}\textbf{B}&-\textbf{A}\\\
-\textbf{A}&\textbf{B}\end{bmatrix}^{-1}=\begin{bmatrix}(\textbf{B}-\textbf{A}\textbf{B}^{-1}\textbf{A})^{-1}&(\textbf{B}-\textbf{A}\textbf{B}^{-1}\textbf{A})^{-1}\textbf{A}\textbf{B}^{-1}\\\
(\textbf{B}-\textbf{A}\textbf{B}^{-1}\textbf{A})^{-1}\textbf{A}\textbf{B}^{-1}&(\textbf{B}-\textbf{A}\textbf{B}^{-1}\textbf{A})^{-1}\end{bmatrix}$
(7.28)

Using this identity, we can multiply out the matrix product in the exponential
and obtain the following reduced form of Eq. 7.27:

$\sqrt{\frac{\det(\bm{a}_{e})\det(\bm{a}_{g})}{\det(\textbf{B})\det(\textbf{B}-\textbf{AB}^{-1}\textbf{A})}}\exp\left(i\left[\textbf{d}^{T}\textbf{G}\textbf{d}-\textbf{d}^{T}\textbf{GS}(\textbf{B}-\textbf{A})^{-1}\textbf{S}^{T}\textbf{G}^{T}\textbf{d}-\textbf{W}^{T}(\textbf{B}+\textbf{A})^{-1}\textbf{W}\right]\right)$
(7.29)

Finally, for the sake of further simplification, we note that

$\displaystyle\textbf{E}+\textbf{S}^{T}\textbf{GS}=\textbf{B}-\textbf{A}$
(7.30) $\displaystyle\Rightarrow\quad$
$\displaystyle\textbf{S}^{T}\textbf{G}=(\textbf{B}-\textbf{A})\textbf{S}^{\text{T}}-\textbf{ES}^{T}$
(7.31) $\displaystyle\Rightarrow\quad$
$\displaystyle\textbf{d}^{T}\textbf{GS}(\textbf{B}-\textbf{A})^{-1}\textbf{S}^{T}\textbf{G}^{T}\textbf{d}=\textbf{d}^{T}\textbf{G}\textbf{d}-\textbf{d}^{T}\textbf{GS}(\textbf{B}-\textbf{A})^{-1}\textbf{ES}^{T}\textbf{d}$
(7.32)

Finally, the correlation function evaluates to

$\sqrt{\frac{\det(\bm{a}_{e})\det(\bm{a}_{g})}{\det(\textbf{B})\det(\textbf{B}-\textbf{AB}^{-1}\textbf{A})}}\exp\left(i\left[\textbf{d}^{T}\textbf{GS}(\textbf{B}-\textbf{A})^{-1}\textbf{ES}^{T}\textbf{d}-\textbf{W}^{T}(\textbf{B}+\textbf{A})^{-1}\textbf{W}\right]\right)$
(7.33)

A reader familiar with the standard spin-boson model may wonder why we have
not utilized the polaron transformation[10], which is standard procedure in
the derivation of Marcus rates. In the method outlined above, a polaron
transformation will lead to exactly the same result obtained in Eq. 7.33.
Although the polaron transformation would modify the interstate coupling
operator, it would analogously modify the position basis used in the expansion
of Eq. 7.17.

###  7.2 Derivation of Eq. 2.31

The correlation function we now need to compute is

$C(t^{\prime},t^{\prime\prime})=\text{Tr}\left[e^{-\beta\hat{H}_{B}^{(g)}}e^{i\hat{H}_{B}^{(e)}t^{\prime}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(g)}(t^{\prime}-t^{\prime\prime})}e^{i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(e)}t^{\prime\prime}}\right]$
(7.34)

Here, we must do 8 integrations over space instead of 4, with 4 excited state
coordinate sets x, y, z, w and 4 ground state coordinate sets
$\textbf{x}^{\prime}$, $\textbf{y}^{\prime}$, $\textbf{z}^{\prime}$,
$\textbf{w}^{\prime}$:

$\displaystyle\int_{\mathbb{R}^{n}}dx^{n}\,\int_{\mathbb{R}^{n}}dy^{n}\,\int_{\mathbb{R}^{n}}dz^{n}\,\int_{\mathbb{R}^{n}}dw^{n}\,\int_{\mathbb{R}^{n}}dx^{\prime
n}\,\int_{\mathbb{R}^{n}}dy^{\prime n}\,\int_{\mathbb{R}^{n}}dz^{\prime
n}\,\int_{\mathbb{R}^{n}}dw^{\prime n}\,$ (7.35)
$\displaystyle\left\langle\textbf{x}\middle|\textbf{x}^{\prime}\right\rangle\left\langle\textbf{x}^{\prime}\middle|e^{-\beta\hat{H}_{B}^{(g)}}\middle|\textbf{y}^{\prime}\right\rangle\left\langle\textbf{y}^{\prime}\middle|\textbf{y}\right\rangle\left\langle\textbf{y}\middle|e^{i\hat{H}_{B}^{(e)}t^{\prime}}e^{-i\textbf{W}^{T}\hat{\textbf{x}}}\middle|\textbf{z}\right\rangle\left\langle\textbf{z}\middle|\textbf{z}^{\prime}\right\rangle\left\langle\textbf{z}^{\prime}\middle|e^{-i\hat{H}_{B}^{(g)}(t^{\prime}-t^{\prime\prime})}\middle|\textbf{w}^{\prime}\right\rangle\langle\textbf{w}^{\prime}|\textbf{w}\rangle\left\langle\textbf{w}\middle|e^{i\textbf{W}^{T}\hat{\textbf{x}}}e^{-i\hat{H}_{B}^{(e)}t^{\prime\prime}}\middle|\textbf{x}\right\rangle$
(7.36) $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{n}}dx^{n}\,\int_{\mathbb{R}^{n}}dy^{n}\,\int_{\mathbb{R}^{n}}dz^{n}\,\int_{\mathbb{R}^{n}}dw^{n}$
(7.37)
$\displaystyle\left\langle\textbf{Sx}+\textbf{d}\middle|e^{-\beta\hat{H}_{B}^{(g)}}\middle|\textbf{Sy}+\textbf{d}\right\rangle\left\langle\textbf{y}\middle|e^{i\hat{H}_{B}^{(e)}t^{\prime}}\middle|\textbf{z}\right\rangle\left\langle\textbf{Sz}+\textbf{d}\middle|e^{-i\hat{H}_{B}^{(g)}(t^{\prime}-t^{\prime\prime})}\middle|\textbf{Sw}+\textbf{d}\right\rangle\left\langle\textbf{w}\middle|e^{-i\hat{H}_{B}^{(e)}t^{\prime\prime}}\middle|\textbf{x}\right\rangle
e^{i\textbf{W}^{T}{\textbf{w}}}e^{-i\textbf{W}^{T}{\textbf{z}}}$ (7.38)

Now we look at the individual terms:

$\displaystyle\left\langle\textbf{Sx}+\textbf{d}\middle|e^{-\beta\hat{H}_{B}^{(g)}}\middle|\textbf{Sy}+\textbf{d}\right\rangle$
$\displaystyle=\sqrt{\frac{\det\bm{a}_{g}(-i\beta)}{(2\pi i)^{N}}}\times$
(7.39) $\displaystyle
e^{i\left[\frac{1}{2}\left((\textbf{S}\bm{x}+\textbf{d})^{T}\bm{b}_{g}(-i\beta)(\textbf{S}\bm{x}+\textbf{d})+(\textbf{S}\bm{y}+\textbf{d})^{T}\bm{b}_{g}(-i\beta)(\textbf{S}\bm{y}+\textbf{d})\right)-(\textbf{S}\bm{x}+\textbf{d})^{T}\bm{a}_{g}(-i\beta)(\textbf{S}\bm{y}+\textbf{d})\right]}$
(7.40)
$\displaystyle\left\langle\textbf{y}\middle|e^{i\hat{H}_{B}^{(e)}t^{\prime}}\middle|\textbf{z}\right\rangle$
$\displaystyle=\sqrt{\frac{\det\bm{a}_{e}(-t^{\prime})}{(2\pi
i)^{N}}}e^{i\left[\frac{1}{2}\left(\bm{y}^{T}\bm{b}_{e}(-t^{\prime})\bm{y}+\bm{z}^{T}\bm{b}_{e}(-t^{\prime})\bm{z}\right)-\bm{y}^{T}\bm{a}_{e}(-t^{\prime})\bm{z}\right]}$
(7.41)
$\displaystyle\left\langle\textbf{Sz}+\bm{\gamma}\middle|e^{-i\hat{H}_{B}^{(g)}(t^{\prime}-t^{\prime\prime})}\middle|\textbf{Sw}+\bm{\gamma}\right\rangle$
$\displaystyle=\sqrt{\frac{\det\bm{a}_{g}(t^{\prime}-t^{\prime\prime})}{(2\pi
i)^{N}}}\times$ (7.42) $\displaystyle
e^{i\left[\frac{1}{2}\left((\textbf{S}\bm{z}+\textbf{d})^{T}\bm{b}_{g}(t^{\prime}-t^{\prime\prime})(\textbf{S}\bm{z}+\textbf{d})+(\textbf{S}\bm{w}+\textbf{d})^{T}\bm{b}_{g}(t^{\prime}-t^{\prime\prime})(\textbf{S}\bm{w}+\textbf{d})\right)-(\textbf{S}\bm{z}+\textbf{d})^{T}\bm{a}_{g}(t^{\prime}-t^{\prime\prime})(\textbf{S}\bm{w}+\textbf{d})\right]}$
(7.43)
$\displaystyle\left\langle\textbf{w}\middle|e^{i\hat{H}_{B}^{(e)}t^{\prime}}\middle|\textbf{x}\right\rangle$
$\displaystyle=\sqrt{\frac{\det\bm{a}_{e}(t^{\prime\prime})}{(2\pi
i)^{N}}}e^{i\left[\frac{1}{2}\left(\bm{w}^{T}\bm{b}_{e}(t^{\prime\prime})\bm{w}+\bm{x}^{T}\bm{b}_{e}(t^{\prime\prime})\bm{x}\right)-\bm{w}^{T}\bm{a}_{e}(t^{\prime\prime})\bm{x}\right]}$
(7.44)

Collecting bilinear and linear terms, we can write the integrand (sans
prefactor) as

$\displaystyle\exp\left(i\left[\frac{1}{2}\begin{bmatrix}\bm{x}^{T}&\bm{y}^{T}&\bm{z}^{T}&\bm{w}^{T}\end{bmatrix}\raisebox{-3.0pt}{\huge{$\Sigma$}}\begin{bmatrix}\bm{x}\\\
\bm{y}\\\ \bm{z}\\\
\bm{w}\end{bmatrix}+\begin{bmatrix}\bm{x}^{T}&\bm{y}^{T}&\bm{z}^{T}&\bm{w}^{T}\end{bmatrix}\begin{bmatrix}\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}-\textbf{W}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}+\textbf{W}\end{bmatrix}+\textbf{d}^{T}\left(\textbf{G}_{\beta}+\textbf{G}\right)\textbf{d}\right]\right)$
(7.45)

where

$\displaystyle\raisebox{-3.0pt}{\huge{$\Sigma$}}=\begin{bmatrix}\textbf{S}^{T}\bm{b}_{g}(-i\beta)\textbf{S}+\bm{b}_{e}(t^{\prime\prime})&-\textbf{S}^{T}\bm{a}_{g}(-i\beta)\textbf{S}&0&-\bm{a}_{e}(t^{\prime\prime})\\\
-\textbf{S}^{T}\bm{a}_{g}(-i\beta)\textbf{S}&\textbf{S}^{T}\bm{b}_{g}(-i\beta)\textbf{S}+\bm{b}_{e}(-t^{\prime})&-\bm{a}_{e}(-t^{\prime})&0\\\
0&-\bm{a}_{e}(-t^{\prime})&\bm{b}_{e}(-t^{\prime})+\textbf{S}^{T}\bm{b}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}&-\textbf{S}^{T}\bm{a}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}\\\
-\bm{a}_{e}(t^{\prime\prime})&0&-\textbf{S}^{T}\bm{a}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}&\bm{b}_{e}(t^{\prime\prime})+\textbf{S}^{T}\bm{b}_{g}(t^{\prime}-t^{\prime\prime})\textbf{S}\end{bmatrix}$
(7.46)
$\displaystyle\textbf{G}=\bm{b}_{g}(t^{\prime}-t^{\prime\prime})-\bm{a}_{g}(t^{\prime}-t^{\prime\prime})$
(7.47)
$\displaystyle\textbf{G}_{\beta}=\bm{b}_{g}(-i\beta)-\bm{a}_{g}(-i\beta)$
(7.48)

As before, we switch to a single integral over a $4n$-dimensional space:

$\displaystyle\int_{\mathbb{R}^{4n}}d\alpha^{4n}\,\exp\left(-\frac{1}{2}\bm{\alpha}^{T}\left(-i\raisebox{-3.0pt}{\huge{$\Sigma$}}\right)\bm{\alpha}+i\,\bm{\alpha}^{T}\begin{bmatrix}\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}-\textbf{W}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}+\textbf{W}\end{bmatrix}+i\,\textbf{d}^{T}\left(\textbf{G}_{\beta}+\textbf{G}\right)\textbf{d}\right)$
(7.49)

Finally, the integral (with prefactors) evaluates to

$\displaystyle\sqrt{\frac{\det\left[\bm{a}_{g}(-i\beta)\bm{a}_{e}(-t^{\prime})\bm{a}_{g}(t^{\prime}-t^{\prime\prime})\bm{a}_{e}(t^{\prime\prime})\right]}{\det\Sigma}}\exp\left(i\,\textbf{d}^{T}\left(\textbf{G}_{\beta}+\textbf{G}\right)\textbf{d}-\frac{1}{2}i\begin{bmatrix}\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}-\textbf{W}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}+\textbf{W}\end{bmatrix}^{T}\raisebox{-3.0pt}{\huge{$\Sigma$}}^{-1}\begin{bmatrix}\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}_{\beta}\textbf{d}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}-\textbf{W}\\\
\textbf{S}^{T}\textbf{G}\textbf{d}+\textbf{W}\end{bmatrix}\right)$ (7.50)

###  7.3 A Note on Numerical Implementation

In the Hamiltonian defined in Eq. 2.34, the number of bath coordinates is
theoretically infinite, resulting in a bath with a continuous spectrum.
Furthermore, the frequencies present in that spectrum are unbounded (i.e. the
bath contains modes of arbitrarily large frequency). In order to perform a
Fermi’s Golden Rule calculation, we must truncate and discretize this bath.
The procedure used for our simulations is as follows: First, we choose a
cutoff frequency $\omega_{C}$. Then, we choose $M$ frequencies evenly spaced
in the interval $(0,\omega_{C})$ and calculate the density of modes as
$\omega_{C}/M$. Finally, these parameters can be used to calculate the
harmonic matrices required in Eq. 2.38. Note that one need not use a uniform
patterning of frequency space, and a nonuniform discretization may well
provide faster convergence of computed dynamics as the number of modes is
increased. However, a uniform discretization proved to be adequate for our
needs. In this scheme, good convergence can be achieved with as few as 40 bath
modes.

## References

  * [1] M. Bixon and J. Jortner  “From Isolated Molecules to Biomolecules”  In _Advances in Chemical Physics_ 106, 1999, pp. 35–203
  * [2] P.. Barbara, T.. Meyer and M.. Ratner  “Contemporary Issues in Electron Transfer Research”  In _Journal of Physical Chemistry_ 100, 1996, pp. 13148–13168
  * [3] T. Nelson, S. Fernandez-Alberti, A.. Roitberg and S. Tretiak  “Nonadiabatic Excited-State Molecular Dynamics: Modeling Photophysics in Organic Conjugated Materials”  In _Accounts of Chemical Research_ 47, 2014, pp. 1155–1164
  * [4] E.. d’Agliano, P. Kumar, W. Schaich and H. Suhl  “Brownian motion model of the interactions between chemical species and metallic electrons: Bootstrap derivation and parameter evaluation”  In _Physical Review B_ 11, 1975, pp. 2122–2143  DOI: 10.1103/PhysRevB.11.2122
  * [5] Barton B. Smith and James T. Hynes  “Electronic friction and electron transfer rates at metallic electrodes”  In _Journal of Chemical Physics_ 99.9, 1993, pp. 6517–6530  DOI: 10.1063/1.465843
  * [6] M. Head-Gordon and J.. Tully  “Molecular dynamics with electronic frictions”  In _Journal of Chemical Physics_ 103, 1995, pp. 10137–10145
  * [7] Yuhui Huang, Charles T. Rettner, Daniel J. Auerbach and Alec M. Wodtke  “Vibrational Promotion of Electron Transfer”  In _Science_ 290.5489, 2000, pp. 111–114  DOI: 10.1126/science.290.5489.111
  * [8] A. Garg, J.. Onuchi and V. Ambegaokar  “Effect of friction on electron transfer in biomolecules”  In _Journal of Chemical Physics_ 83, 1985, pp. 4491–4503
  * [9] A. Leggett et al.  “Dynamics of the dissipative two-state system”  In _Reviews of Modern Physics_ 59, 1987, pp. 1–85
  * [10] A. Nitzan  “Chemical Dynamics in Condensed Phases”  USA: Oxford University Press, 2006
  * [11] T.. Berkelbach, T.. Markland and D.. Reichman  “Reduced density matrix hybrid approach: Application to electronic energy transfer”  In _Journal of Chemical Physics_ 136, 2012, pp. 084104
  * [12] J.. Onuchic and D.. Beratan  “A predictive theoretical model for electron tunneling pathways in proteins”  In _Journal of Chemical Physics_ 92, 1990, pp. 722
  * [13] R.. Marcus  “ON THE THEORY OF OXIDATION-REDUCTION REACTIONS INVOLVING ELECTRON TRANSFER. V. COMPARISON AND PROPERTIES OF ELECTROCHEMICAL AND CHEMICAL RATE CONSTANTS”  In _Journal of Physical Chemistry_ 67, 1963, pp. 853–857
  * [14] L.. Zusman  “Outer-sphere electron transfer in polar solvents”  In _Chemical Physics_ 49, 1980, pp. 295–304
  * [15] J.. Hynes  “Outer-sphere electron transfer reactions and frequency-dependent friction”  In _Journal of Physical Chemistry_ 90, 1986, pp. 3701–3706
  * [16] John E. Straub and Bruce J. Berne  “A statistical theory for the effect of nonadiabatic transitions on activated processes”  In _Journal of Chemical Physics_ 87.10, 1987, pp. 6111–6116
  * [17] S. Jang and M.. Newton  “Theory of Torsional non-Condon Electron Transfer: A Generalized Spin-Boson Hamiltonian and Its Nonadiabatic Limit Solution”  In _Journal of Chemical Physics_ 122, 2005, pp. 024501
  * [18] Emile S. Medvedev and Alexei A. Stuchebrukhov  “Inelastic tunneling in long-distance biological electron transfer reactions”  In _Journal of Chemical Physics_ 107.10, 1997, pp. 3821–3831
  * [19] B Göhler et al.  “Spin selectivity in electron transmission through self-assembled monolayers of double-stranded DNA”  In _Science_ 331.6019 American Association for the Advancement of Science, 2011, pp. 894–897
  * [20] Ron Naaman, Yossi Paltiel and David H Waldeck  “Chiral Molecules and the Spin Selectivity Effect”  In _Journal of Physical Chemistry Letters_ 11.9 ACS Publications, 2020, pp. 3660–3666
  * [21] Thomas P Fay  “Chirality-induced spin coherence in electron transfer reactions”  In _The Journal of Physical Chemistry Letters_ 12.5 ACS Publications, 2021, pp. 1407–1412
  * [22] Jonas Fransson  “Vibrational origin of exchange splitting and chiral-induced spin selectivity”  In _Physical Review B_ 102.23 APS, 2020, pp. 235416
  * [23] Lior Opennheim and Karen Michaeli  “Incoherent chiral-induced spin selectivity”, 2021  arXiv:2112.13561
  * [24] Y. Wu, G. Miao and J.. Subotnik  “Chemical Reaction Rates for Systems with Spin-Orbit Coupling and an Odd Number of Electrons: Does Berry’s Phase Lead to Meaningful Spin-Dependent Nuclear Dynamics for a Two State Crossing?”  In _Journal of Physical Chemistry A_ 124, 2020, pp. 7355–7372
  * [25] M.. Berry and J.. Robbins  “Chaotic classical and half-classical adiabatic reactions: geometric magnetism and deterministic friction”  In _Proceedings of the Royal Society of London A_ 442, 1993, pp. 659–672  DOI: https://doi.org/10.1098/rspa.1993.0127
  * [26] Joseph Subotnik et al.  “A demonstration of consistency between the quantum classical Liouville equation and Berry’s phase and curvature for the case of complex Hamiltonians”  In _Journal of Chemical Physics_ 151.7, 2019, pp. 074113  DOI: 10.1063/1.5116210
  * [27] Xuezhi Bian et al.  “Modeling Spin-Dependent Nonadiabatic Dynamics with Electronic Degeneracy: A Phase-Space Surface-Hopping Method” PMID: 35926097  In _The Journal of Physical Chemistry Letters_ 13.32, 2022, pp. 7398–7404  DOI: 10.1021/acs.jpclett.2c01802
  * [28] Roel Tempelaar and David R. Reichman  “Generalization of fewest-switches surface hopping for coherences”  In _Journal of Chemical Physics_ 148.10, 2018, pp. 102309  DOI: 10.1063/1.5000843
  * [29] Suraj S. Chandran et al.  “Electron transfer and spin–orbit coupling: Can nuclear motion lead to spin selective rates?”  In _The Journal of Chemical Physics_ 156.17, 2022, pp. 174113  DOI: 10.1063/5.0086554
  * [30] A.O. Caldeira and A.J. Leggett  “Path integral approach to quantum Brownian motion”  In _Physica A_ 121.3, 1983, pp. 587–616
  * [31] Reuven Ianconescu and Eli Pollak  “Photoinduced Cooling of Polyatomic Molecules in an Electronically Excited State in the Presence of Dushinskii Rotations”  In _The Journal of Physical Chemistry A_ 108.39, 2004, pp. 7778–7784  DOI: 10.1021/jp037739q
  * [32] F. Duschinsky  “Zur Deutung der Elektronenspektren mehratomiger Moleküle”  In _Acta Physicochimica U.R.S.S_ 7.4, 1937, pp. 551–566
  * [33] Aleksandr O. Lykhin et al.  “Nonadiabatic transition state theory: Application to intersystem crossings in the active sites of metal-sulfur proteins”  In _International Journal of Quantum Chemistry_ 116.10, 2016, pp. 750–761  DOI: 10.1002/qua.25124
  * [34] Vsevolod D. Dergachev et al.  “NAST: Nonadiabatic Statistical Theory Package for Predicting Kinetics of Spin-Dependent Processes”  In _Topics in Current Chemistry_ 380.15, 2022
  * [35] Thomas P Fay and David T Limmer  “Origin of Chirality Induced Spin Selectivity in Photoinduced Electron Transfer”  In _Nano letters_ 21.15 ACS Publications, 2021, pp. 6696–6702
  * [36] R Naaman and David H Waldeck  “Chiral-induced spin selectivity effect”  In _Journal of Physical Chemistry Letters_ 3.16 ACS Publications, 2012, pp. 2178–2187