figurec {}$missing\text{\,}\mathrm{missing}$

# Extending non-adiabatic rate theory to strong electronic couplings in the
Marcus inverted regime

Thomas P. Fay<EMAIL_ADDRESS>Aix Marseille University, CNRS, ICR,
13397 Marseille, France

###### Abstract

Electron transfer reactions play an essential role in many chemical and
biological processes. Fermi’s Golden rule, which assumes that the coupling
between electronic states is small, has formed the foundation of electron
transfer rate theory, however in short range electron/energy transfer
reactions this coupling can become very large, and therefore Fermi’s Golden
Rule fails to make even qualitatively accurate rate predictions. In this paper
I present a simple modified Golden Rule theory to describe electron transfer
in the Marcus inverted regime at arbitrarily large electronic coupling
strengths. The theory is based on an optimal global rotation of the diabatic
states, which makes it compatible with existing methods for calculating Golden
Rule rates that can account for nuclear quantum effects with anharmonic
potentials. Furthermore the Optimal Golden Rule (OGR) theory can also combined
with analytic theories for non-adiabatic rates, such as Marcus theory and
Marcus-Levich-Jortner theory, offering clear physical insight into strong
electronic coupling effects in non-adiabatic processes. OGR theory is also
tested on a large set of spin-boson models and an anharmonic model against
exact quantum dynamics calculations, where it performs well, correctly
predicting rate turnover at large coupling strengths. Finally, an example
application to a BODIPY-Anthracene photosensitizer reveals that strong
coupling effects inhibit excited state charge recombination in this system,
reducing the rate of this process by a factor of four. Overall OGR theory
offers a new approach to calculating electron transfer rates at strong
couplings, offering new physical insight into a range of non-adiabatic
processes.

##  I Introduction

Figure 1: An illustration of the diabatic (dashed lines, $V_{\mathrm{A}}$ and
$V_{\mathrm{B}}$) and adiabatic/Born-Oppenheimer (solid lines, $E_{+}$ and
$E_{-}$) potential energy curves for an electron transfer reaction A) in the
Marcus normal regime and B) in the Marcus inverted regime. The reactant and
product wells are also labelled with R and P respectively.

The accurate calculation of electron transfer rates is important in many areas
of chemistry, biology and materials science.[1, 2, 3, 4] Most theories for
calculating these rate constants assume that the coupling between diabatic
electronic states is small, meaning Fermi’s Golden Rule perturbation theory
can be used.[1] However in many real systems, ranging from photo-active
proteins[5] to organic semiconductors[6] to synthetic photocatalysts[7] the
electronic coupling can exceed thermal energy at room temperature,[8]
$\sim\\!\quantity{200}{cm^{-1}}$, sometimes by over a factor of 10.[7, 9] In
these cases strong electronic coupling effects can be expected to become
important, and the naive application of Fermi’s Golden rule can introduce
errors of several orders of magnitude in rate constants, particularly for
reactions in the Marcus inverted regime.[10] Thus the aim of this work is to
address this challenge and develop a simple theory for the calculation of
electron transfer rates at arbitrarily strong electronic coupling in the
inverted regime.

For rates in the Golden Rule (GR) limit, there have been impressive
developments in recent years, with methods emerging that can accurately
predict rate constants including nuclear quantum effects for anharmonic
potential energy surfaces.[11, 12, 13, 14, 15, 16, 17, 18, 19] Many of these
methods are obtained using imaginary time path integral or instanton
approaches, such as Wolynes theory,[20, 21] non-adiabatic instanton
methods,[22, 12] Golden Rule Quantum Transition State Theory (GRQTST),[11] and
the Linear Golden Rule (LGR) method.[14] The success of these methods in the
GR limit motivates the development of theories for non-adiabatic rates beyond
the weak diabatic coupling regime which are amenable to these path integral
and instanton approaches. There have also been significant steps made towards
this goal and methods have been developed which can account for strong
electronic coupling effects, such as higher-order instanton methods[23, 24,
18] and the Non-Adiabatic Quantum Instanton (NAQI) method,[25] as well as the
development of an interpolation formula technique.[10]

Existing methods for calculating strong-coupling effects on electron transfer
rates generally use the fact that in the Marcus normal regime, where
decreasing the free energy change of reaction lowers the activation energy,
the reaction can be described as adiabatic, with the electron transfer
occurring continuously on a single electronic state (this is illustrated in
Fig. 1A).[10] This means strong-coupling effects can be accounted for in the
normal regime by using established techniques for calculating rate constants
on single electronic surfaces,[26] such as ring polymer molecular dynamics
(RPMD),[27] centroid molecular dynamics (CMD),[28] and instanton methods.[29,
30] However the Marcus inverted regime, where decreasing the reaction free
energy increases the activation barrier, still poses a challenge, and the
aforementioned methods cannot be applied in this regime (with the possible
exception of the higher-order instanton approach[23]). This is because in the
inverted regime, even when diabatic coupling is strong, the reactant and
product states still correspond to different adiabatic electronic states, so
in this regime the reaction always involves non-adiabatic transitions between
different electronic states even in the strong electronic coupling limit.[10]
This is illustrated in Fig. 1. Direct semiclassical non-adiabatic dynamics
approaches, such as surface hopping[31, 32, 33, 34, 35, 36] and mean-field-
type methods,[37, 38, 39, 40] could provide access to strong coupling effects
deep in the inverted regime. But these approaches rely on a classical
approximation for the motion of nuclei, and in the inverted regime nuclear
quantum effects can enhance rates by many orders of magnitude.[10, 14, 9]
Furthermore rates of processes in the inverted regime can be very slow, with
lifetimes exceeding {1}μs being typical. These factors pose a significant
limitation on the application of direct semiclassical dynamics methods to
calculate rates in the inverted regime for real molecular systems. Exact
quantum dynamical methods can be used to calculate rates for certain
simplified Hamiltonians, such as the (Multi-Layer) Multi-Configurational Time-
Dependent Hartree method,[41, 42] Hierarchical Equations of Motion,[43, 44]
and real-time path-integral methods,[45, 46] but these methods cannot be
applied in general anharmonic systems.

Given that in real systems electronic couplings can be very large for short-
range electron transfer,[9, 7, 47] there is clearly a need to address strong
diabatic coupling effects in the inverted regime. Currently the only existing
method that can potentially be applied to these problems is Trenins and
Richardson’s fourth-order instanton method (which was originally developed for
the normal regime),[23] although the reformulation of their theory for the
inverted regime has not yet been presented. Ideally any new theory for this
regime should synergise with existing imaginary time path integral and
instanton approaches, to enable the inclusion of nuclear quantum effects,
which are known to be very important in the inverted regime. In this paper, I
outline a theory, based on an optimal rotation of diabatic states, which
enables the calculation of electron transfer rates deep in the inverted regime
with arbitrarily large diabatic coupling, and which is directly compatible
with existing methods for calculating non-adiabatic rates in the Golden Rule
regime. I also show how this approach can be used to obtain analytic modified
Marcus theory and Marcus-Levich-Jortner theory rate expressions including
strong-coupling effects in the inverted regime, and demonstrate how in a
BODIPY-Anthracene photosensitizer strong-coupling effects inhibit charge
recombination.

##  II Theory

###  II.1 Non-adiabatic rates

The aim of this work is to find a rate expression for the transfer between
diabatic electronic states $\mathrm{A}$ and $\mathrm{B}$. The Hamiltonian for
this system is given by

$\displaystyle\hat{H}=\outerproduct{\mathrm{A}}{\mathrm{A}}\hat{H}_{\mathrm{A}}+\outerproduct{\mathrm{B}}{\mathrm{B}}\hat{H}_{\mathrm{B}}+\hat{\Delta}(\outerproduct{\mathrm{A}}{\mathrm{B}}+\outerproduct{\mathrm{B}}{\mathrm{A}})$
(1)

where $\hat{H}_{J}=\hat{T}+V_{J}(\hat{{\bf\it q}})$ is the diabatic state
Hamiltonian for state $J=\mathrm{A},\mathrm{B}$ and
$\hat{\Delta}=\Delta(\hat{{\bf\it q}})$ is the diabatic coupling between these
states. The diabatic potentials $V_{J}({{\bf\it q}})$ and the coupling
${\Delta}({{\bf\it q}})$ are both assumed to be real-valued functions of the
nuclear positions ${\bf\it q}$. We will be interested in the limit where
$\Delta({\bf\it q})$ can be arbitrarily large.

In order to find the rate constant for the transition from $\mathrm{A}$ to
$\mathrm{B}$, we first define the reactant and product side operators, which
measure the populations of the reactants and products, as
$\hat{P}_{\mathrm{r}}=\outerproduct{\mathrm{A}}{\mathrm{A}}$ and
$\hat{P}_{\mathrm{p}}=\outerproduct{\mathrm{B}}{\mathrm{B}}$ respectively.
These operators define a diabatic state dividing surface between the reactant
and product. In the long time limit the expectation values of observables
should obey first-order kinetic equations

$\displaystyle\derivative{t}\expectationvalue{P_{\mathrm{r}}(t)}$
$\displaystyle=-k\expectationvalue{P_{\mathrm{r}}(t)}+k^{\prime}\expectationvalue{P_{\mathrm{p}}(t)}$
(2a) $\displaystyle\derivative{t}\expectationvalue{P_{\mathrm{p}}(t)}$
$\displaystyle=+k\expectationvalue{P_{\mathrm{r}}(t)}-k^{\prime}\expectationvalue{P_{\mathrm{p}}(t)}.$
(2b)

With this the forward rate constant can be found as[10]

$\displaystyle
k=\lim_{t\to\infty}\frac{\derivative{t}\expectationvalue{P_{\mathrm{p}}(t)}}{1-\expectationvalue{P_{\mathrm{p}}(t)}/\expectationvalue{P_{\mathrm{p}}(\infty)}}.$
(3)

$\expectationvalue*{P_{\mathrm{p}}(t)}=\Tr[e^{+i\hat{H}t/\hbar}\hat{P}_{\mathrm{p}}e^{-i\hat{H}t/\hbar}\hat{P}_{\mathrm{r}}^{\mathrm{K}}]/Q_{\mathrm{r}}$
is the Kubo-transformed side-side correlation function, where the Kubo-
transformed operator $\hat{O}^{\mathrm{K}}$ is given by
$\hat{O}^{\mathrm{K}}=(1/\beta)\int_{0}^{\beta}e^{-(\beta-\lambda)\hat{H}}\hat{O}e^{-\lambda\hat{H}}$
and the reactant partition function is
$Q_{\mathrm{r}}={\Tr}[\hat{P}_{\mathrm{r}}e^{-\beta\hat{H}}]$. If the side-
side correlation function derivative plateaus on a time $t_{\mathrm{p}}$ and
there is minimal population transfer up to this time,then the rate constant
can be simplified to[48]

$\displaystyle k\simeq\frac{1}{Q_{\mathrm{r}}}\lim_{t\to
t_{\mathrm{p}}}\derivative{t}\Tr[e^{+i\hat{H}t/\hbar}\hat{P}_{\mathrm{p}}e^{-i\hat{H}t/\hbar}\hat{P}_{\mathrm{r}}^{\mathrm{K}}].$
(4)

With some simple manipulations, the rate constant can alternatively be written
in terms of the flux-flux correlation function as

$\displaystyle
k\simeq\frac{1}{2Q_{\mathrm{r}}}\int_{-t_{\mathrm{p}}}^{t_{\mathrm{p}}}\Tr[e^{+i\hat{H}t/\hbar}\hat{F}e^{-i\hat{H}t/\hbar}\hat{F}^{\mathrm{K}}]\differential{t},$
(5)

where the flux operator
$\hat{F}=\frac{i}{\hbar}[\hat{H},\hat{P}_{\mathrm{r}}]$. The standard Golden
rule rate expression can be obtained by expanding this to lowest order in the
coupling term $\hat{\Delta}$. This can be achieved by partitioning the
Hamiltonian into $\hat{H}_{0}+\hat{V}$ with
$\hat{H}_{0}=\outerproduct{\mathrm{A}}{\mathrm{A}}\hat{H}_{\mathrm{A}}+\outerproduct{\mathrm{B}}{\mathrm{B}}\hat{H}_{\mathrm{B}}$
and
$\hat{V}=\hat{\Delta}(\outerproduct{\mathrm{A}}{\mathrm{B}}+\outerproduct{\mathrm{B}}{\mathrm{A}})$
and noting that
$\hat{F}=(i/\hbar)\hat{\Delta}(\outerproduct{\mathrm{B}}{\mathrm{A}}-\outerproduct{\mathrm{A}}{\mathrm{B}})$
and $e^{z\hat{H}}=e^{z\hat{H}_{0}}+\mathcal{O}(\Delta)$, which yields the
Golden Rule rate expression,[49]

$\displaystyle
k_{\mathrm{GR}}\\!\\!=\\!\\!\frac{1}{\hbar^{2}Q_{\mathrm{A}}}\int_{-\infty}^{\infty}\\!\\!\\!\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{\mathrm{A}}-i\hat{H}_{\mathrm{A}}t/\hbar}\hat{\Delta}e^{+i\hat{H}_{\mathrm{B}}t/\hbar}\hat{\Delta}]\differential{t}\\!,$
(6)

where $\Tr_{\mathrm{n}}[\cdots]$ denotes a trace just over nuclear degrees of
freedom, and $Q_{J}=\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{J}}]$ for
$J=\mathrm{A},\mathrm{B}$. The Golden Rule rate is accurate in the limit where
the diabatic coupling $\Delta$ is small, but can deviate significantly from
the exact rate beyond this limit. The aim of the following is to find an
alternative approximation to the full rate constant expression, Eq. (5), which
accurately captures strong diabatic coupling effects.

###  II.2 A modified Golden Rule rate

For any sensible choice of dividing surface that does not radically change the
character of the reactant and product states, the exact rate constant $k$ is
unchanged. This is because the rate constant is simply the slowest relaxation
rate of the entire system. Thus it is possible to choose an alternative
dividing surface, which defines $\hat{P}_{\mathrm{r}}$ and
$\hat{P}_{\mathrm{p}}=1-\hat{P}_{\mathrm{r}}$, for which the coupling between
states is minimised, and perturbation theory can therefore be applied more
accurately. Here I will define the side operators using rotated electronic
states $\ket{\psi_{\pm}(\theta({\bf\it q}))}$ (although the explicit ${\bf\it
q}$ dependence of these states will be dropped shortly), such that the side
operators are given by,

$\displaystyle\hat{P}_{\mathrm{r}}$
$\displaystyle=\outerproduct{\psi_{+}(\theta({{\bf\it
q}}))}{\psi_{+}(\theta({{\bf\it q}}))}$ (7)
$\displaystyle\hat{P}_{\mathrm{p}}$
$\displaystyle=\outerproduct{\psi_{-}(\theta({\bf\it
q}))}{\psi_{-}(\theta({\bf\it q}))}$ (8)

where the rotated states are defined as

$\displaystyle\ket{\psi_{+}(\theta({\bf\it q}))}$
$\displaystyle=\ket{\mathrm{A}}\cos(\theta({\bf\it
q}))+\ket{\mathrm{B}}\sin(\theta({\bf\it q}))$ (9)
$\displaystyle\ket{\psi_{-}(\theta({\bf\it q}))}$
$\displaystyle=-\ket{\mathrm{A}}\sin(\theta({\bf\it
q}))+\ket{\mathrm{B}}\cos(\theta({\bf\it q})).$ (10)

If the rotation angle $\theta({\bf\it q})$ is chosen to diagonalise the Born-
Oppenheimer Hamiltonian $\hat{H}-\hat{T}$ then these states are the adiabatic
states. In this case the rotation angle is given by

$\displaystyle\tan(2\theta({\bf\it q}))=\frac{2\Delta({\bf\it
q})}{V_{\mathrm{A}}({\bf\it q})-V_{\mathrm{B}}({\bf\it q})}$ (11)

and the eigenvalues of $\hat{H}-\hat{T}$ are given by

$\displaystyle E_{\pm}({\bf\it q})=\frac{V_{\mathrm{A}}({\bf\it
q})+V_{\mathrm{B}}({\bf\it q})}{2}\pm\frac{\sqrt{4\Delta({\bf\it
q})^{2}+(V_{\mathrm{A}}({\bf\it q})-V_{\mathrm{B}}({\bf\it q}))^{2}}}{2}$ (12)

but in this basis the kinetic energy operator $\hat{T}$ acquires off-diagonal
elements which become very large in regions of position space where the
adiabatic energies are close. This suggests that choosing the adiabatic state
dividing surface may not be the optimal choice.

The large off-diagonal elements of the kinetic energy operator in the
adiabatic basis originate in the position dependence of the electronic states.
One simple way to remove these is to choose a global rotation angle,
$\theta({\bf\it q})\to\theta$, to define the rotated diabatic basis states
$\ket{\pm}\equiv\ket{\psi_{\pm}(\theta)}$ which can then be used to define the
dividing surface between reactants and products. The Hamiltonian can be
written in terms of these globally rotated diabatic states as

$\displaystyle\hat{H}$
$\displaystyle=\hat{H}_{+}\outerproduct{+}{+}+\hat{H}_{-}\outerproduct{-}{-}+\hat{\Delta}_{+-}(\outerproduct{+}{-}+\outerproduct{-}{+})$
(13)

where the constituent terms are

$\displaystyle\hat{H}_{+}$
$\displaystyle=\hat{H}_{\mathrm{A}}\cos^{2}\theta+\hat{H}_{\mathrm{B}}\sin^{2}\theta+\hat{\Delta}\sin(2\theta)$
(14) $\displaystyle\hat{H}_{-}$
$\displaystyle=\hat{H}_{\mathrm{A}}\sin^{2}\theta+\hat{H}_{\mathrm{B}}\cos^{2}\theta-\hat{\Delta}\sin(2\theta)$
(15) $\displaystyle\hat{\Delta}_{+-}$
$\displaystyle=\hat{\Delta}\cos(2\theta)+\frac{1}{2}(\hat{H}_{\mathrm{B}}-\hat{H}_{\mathrm{A}})\sin(2\theta).$
(16)

With these rotated diabatic states, a perturbative expression for the rate
constant can be obtained by setting
$\hat{H}_{0}=\hat{H}_{+}\outerproduct{+}{+}+\hat{H}_{-}\outerproduct{-}{-}$
and $\hat{V}=\hat{\Delta}_{+-}(\outerproduct{+}{-}+\outerproduct{-}{+})$ and
expanding the resulting expression to lowest order in $\Delta_{+-}$, which
gives a modified Golden Rule rate expression, $k_{\mathrm{GR,\theta}}$,

$\displaystyle
k_{\mathrm{GR,\theta}}\\!=\\!\frac{1}{\hbar^{2}Q_{+}}\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{+}-i\hat{H}_{+}t/\hbar}\hat{\Delta}_{+-}e^{+i\hat{H}_{-}t/\hbar}\hat{\Delta}_{+-}]\\!\differential{t}\\!,$
(17)

where $Q_{\pm}=\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{\pm}}]$. If $\theta=0$ this
reduces to the standard Golden Rule rate expression, but for $\theta\neq 0$
the coupling term $\hat{\Delta}$ is partially folded into reference
Hamiltonian, which in principle enables strong electronic coupling effects to
be accounted for in the rate constant. The question that remains is what is
the optimal choice of mixing angle $\theta$? This is addressed in the next
section.

###  II.3 The optimal choice of $\theta$

The perturbative expansion of the rate constant in $\Delta_{+-}$ will be most
justified when the magnitude of this term is minimised in some sense. There
are many potential measures of how strong $\Delta_{+-}$ is, but one simple
choice is to minimise the expectation value of the square coupling
${\Delta_{+-}^{2}}$ in either in the local equilibrium distribution of
$\ket{+}$ or $\ket{-}$,
$\hat{\rho}_{\pm}^{\mathrm{eq}}=e^{-\beta\hat{H}_{\pm}}/Q_{\pm}$. Deep in the
inverted regime, the fluctuations in
$\Delta\hat{V}=\hat{V}_{\mathrm{B}}-\hat{V}_{\mathrm{A}}$, denoted
$\delta\Delta\hat{V}=\Delta\hat{V}-\expectationvalue{\Delta\hat{V}}_{\pm}$,
will be small compared to $\expectationvalue{\Delta V}_{\pm}^{2}$, so the
$\expectationvalue{\delta\Delta\hat{V}^{2}}_{\pm}$ terms appear in the above
can be ignored. Likewise if the fluctuations in the coupling are also ignored
(which is not necessary if the Condon approximation is also invoked[10] and
$\hat{\Delta}$ is assumed to be a constant), then
$\expectationvalue{\Delta_{+-}^{2}}_{\pm}$ can be approximated as

$\displaystyle\expectationvalue{\Delta_{+-}^{2}}_{\pm}\approx\left(\expectationvalue{\Delta}_{\pm}\sin(2\theta)+\frac{1}{2}\expectationvalue{\Delta
V}_{\pm}\cos(2\theta)\right)^{2}.$ (18)

Now assuming $e^{-\beta\hat{H}_{+}}\approx e^{-\beta\hat{H}_{\mathrm{A}}}$ and
$e^{-\beta\hat{H}_{-}}\approx e^{-\beta\hat{H}_{\mathrm{B}}}$ we find that
$\expectationvalue{\Delta_{+-}^{2}}_{\pm}$ is given approximately by

$\displaystyle\expectationvalue{\Delta_{+-}^{2}}_{+}$
$\displaystyle\approx\left(\expectationvalue{\Delta}_{\mathrm{A}}\sin(2\theta)+\frac{1}{2}\expectationvalue{\Delta
V}_{\mathrm{A}}\cos(2\theta)\right)^{2}$ (19)
$\displaystyle\expectationvalue{\Delta_{+-}^{2}}_{-}$
$\displaystyle\approx\left(\expectationvalue{\Delta}_{\mathrm{B}}\sin(2\theta)+\frac{1}{2}\expectationvalue{\Delta
V}_{\mathrm{B}}\cos(2\theta)\right)^{2},$ (20)

where $\expectationvalue{\cdots}_{J}=\Tr_{\mathrm{n}}[\cdots
e^{-\beta\hat{H}_{J}}/Q_{J}]$. Finally differentiating these expressions with
respect to $\theta$ yields values of $\theta$ which minimise the coupling in
the A and B equilibrium distributions,

$\displaystyle\tan(2\theta_{J})=\frac{2\expectationvalue{\Delta}_{J}}{\expectationvalue{V_{\mathrm{A}}-V_{\mathrm{B}}}_{J}}$
(21)

for $J=\mathrm{A},\mathrm{B}$. Clearly the coupling cannot be minimised
simultaneously for the $\ket{+}$ and $\ket{-}$ states with a global choice of
$\theta$. This suggests that the optimal value of $\theta$ should lie
somewhere in this range.

Figure 2: An illustration of the rotated diabatic states for a single-mode
spin-boson model with $\epsilon=-\Delta A=5$, $\lambda=2$, $m\omega^{2}=1$ and
$\Delta=1$ with $\theta=\theta_{\mathrm{A}}$.

In order to obtain an optimal value, we note that the Golden Rule rate is
essentially a transition state theory, and the higher order terms in
$\Delta_{+-}$ that are neglected correspond primarily to recrossing effects in
the Marcus inverted regime,[50, 51, 49] and therefore $k_{\mathrm{GR},\theta}$
should be an approximate upper bound on the true rate constant in the inverted
regime. This motivates defining the optimal value of $\theta$ as the value
that minimises the rate in the range
$\theta_{\mathrm{B}}\leq\theta\leq\theta_{\mathrm{A}}$, so the optimal value
of the rate constant is given by

$\displaystyle
k_{\mathrm{GR},\theta^{*}}=\underset{\theta\in[\theta_{\mathrm{B}},\theta_{\mathrm{A}}]}{\mathrm{min}}k_{\mathrm{GR},\theta}.$
(22)

Very deep in the inverted regime, assuming $\mathrm{B}$ is the lower energy
state, we expect the rate constant to be minimised when either
$\theta=\theta_{\mathrm{B}}\approx 0$ in the limit of small diabatic coupling
$\Delta$ or at $\theta=\theta_{\mathrm{A}}$ in the strong diabatic coupling
limit. This motivates choosing the optimal rate constant to be the minimum of
$k_{\mathrm{GR},\theta_{\mathrm{A}}}$ and
$k_{\mathrm{GR},\theta_{\mathrm{B}}}$,

$\displaystyle
k_{\mathrm{GR},\theta^{*}}\approx\min(k_{\mathrm{GR},\theta_{\mathrm{A}}},k_{\mathrm{GR},\theta_{\mathrm{B}}}),$
(23)

if calculations for the full range of $\theta$ are prohibitive.

An example of the rotated diabats with $\theta=\theta_{\mathrm{A}}$ are shown
in Fig. 2, where it can be seen that the rotatation raises the energy of the
crossing point between diabats and thereby decreases the rate constant, so we
see qualitatively how this rotation of the diabats can reduce the rate
constant. Fig. 3A shows the behaviour of $k_{\mathrm{GR},\theta}$ as a
function of $\theta$ for a spin-boson model (see Sec. III.1 for details) deep
in the inverted regime, with $\epsilon=-\Delta A=5\lambda$. We see that the
rate constant is minimised at $\theta=\theta_{\mathrm{A}}$ in the range
$\theta_{\mathrm{B}}\leq\theta\leq\theta_{\mathrm{A}}$, although this is not
unconstrained minimum which lies at $\theta>\theta_{\mathrm{A}}$. This minimum
underestimates the exact HEOM rate constant, and the constrained minimum
provides a better estimate of the true rate constant than the true minimum.
This shows the importance of applying a physically motivated constraint on
optimisation of the global rotation angle. Furthermore we see that in the
range $\theta_{\mathrm{B}}\leq\theta\leq\theta_{\mathrm{A}}$
$k_{\mathrm{GR},\theta}$ is a monotonic function of $\theta$, thus the OGR
rate Eq. (22) and the simplified OGR rate Eq. (23) are equivalent, which
motivates the simplified form Eq. (23).

Figure 3: A) Behaviour of the $k_{\mathrm{GR},\theta}$ as a function of
$\theta$ for a spin-boson model with an underdamped Brownian Oscillator
spectral density with $\beta\lambda=20$, $\beta\epsilon=100$,
$\beta\hbar\Omega=4$, $\beta\hbar\gamma=4$ and $\beta\Delta=20$. The exact
HEOM rate constant and the GR rate constants are also shown, together with the
proposed limits for the optimisation of $k_{\mathrm{GR},\theta}$. B) The
optimal values of $\theta$ for an underdamped spin-boson model close to the
activation-less case $\epsilon=1.25\lambda$ and in the inverted regime
$\epsilon=2\lambda$ as a function of $\beta\Delta$. In this example
$\beta\hbar\Omega=\beta\hbar\gamma=4$ and $\beta\lambda=20$.

To further explore the behaviour of $\theta^{*}$, in Fig. 3B the optimal
values of $\theta$ are plotted as a function of $\beta\Delta$ for a spin-boson
model in the inverted regime (details of the model are given in Sec. III.1).
Deep in the inverted regime, $\epsilon\geq 2\lambda$, the rate is minimised at
$\theta=\theta_{B}$ (which is true for all deep-inverted regime models
explored in this work). Conversely closer to the activation-less case,
$\epsilon=1.25\lambda$, a transition is observed where the optimal value of
$\theta$ changes from close to $\theta_{\mathrm{A}}$ to $\theta_{\mathrm{B}}$
as $\Delta$ increases, with $\theta^{*}$ not lying exactly at one of the
bounds for only a relatively small range of values of $\Delta$. This suggests
that the simplified optimisation, Eq. (23), should accurately approximate the
full version as long as the problem in question lies far from the Marcus
activation-less regime.

This approximate optimisation condition can be viewed through the lens of the
higher-order non-adiabatic rate theory in Ref. 23. Without going into detail,
it has been found that within a semiclassical ($\hbar\to 0$) approximation,
the dominant term in the fourth-order term in $\Delta$ in the perturbative
expansion of the rate constant, denoted $k^{(4)}$, is proportional to
$k_{\mathrm{GR}}^{2}$. This is also true in Zusman theory[52] and its
generalisations[53, 54]. This means that the square ratio of the fourth order
contribution to the GR contribution, $(k^{(4)}/k_{\mathrm{GR}})^{2}$, is
approximately proportional to $k_{\mathrm{GR}}^{2}$. Minimising the magnitude
of higher-order contributions to the rate should make the lowest order rate
constant closer the the true rate. By constraining $\theta$ to lie in the
range $\theta_{\mathrm{B}}\leq\theta\leq\theta_{\mathrm{A}}$ we ensure that
the the approximate minimisation of the fourth-order contribution to the rate
constant is constrained to values of $\theta$ which also minimise the coupling
between diabatic states and which preserve the character of the reactant and
product states.

The modified coupling term $\hat{\Delta}_{+-}$ contains the diabatic energy
gap $\hat{H}_{\mathrm{B}}-\hat{H}_{\mathrm{A}}$, which may complicated
evaluation of the rate constant. We note however that the energy gap operator
can be written as

$\displaystyle\hat{H}_{\mathrm{B}}-\hat{H}_{\mathrm{A}}=(\hat{H}_{-}-\hat{H}_{+})\sec(2\theta)+2\Delta\tan(2\theta),$
(24)

and furthermore any terms proportional to $\hat{H}_{-}-\hat{H}_{+}$ do not
contribute to the rate constant $k_{\mathrm{GR},\theta}$. An intuitation for
this at first surprising result can be gained by considering the classical
limit. In this limit transitions between rotated diabats only occur at points
in configuration space ${\bf\it q}$ where diabatic states, $V_{\pm}({\bf\it
q})$, are equal in energy, $V_{+}({\bf\it q})=V_{-}({\bf\it q})$, so
$V_{-}({\bf\it q})-V_{+}({\bf\it q})$ contributions to the coupling vanish
where diabatic transitions occur and they do not contribute to the rate
constant. The generalisation of this to the quantum case is given in Appendix
A. This means the rate constant can be expressed more simply as

$\displaystyle
k_{\mathrm{GR,\theta}}\\!=\\!\frac{\sec^{2}(2\theta)}{\hbar^{2}Q_{+}}\\!\\!\int_{-\infty}^{\infty}\\!\\!\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{+}-i\hat{H}_{+}t/\hbar}\hat{\Delta}e^{+i\hat{H}_{-}t/\hbar}\hat{\Delta}]\differential{t}\\!.$
(25)

One clear advantage of this modified Golden Rule expression is that it has
exactly the same form as the standard Golden Rule rate expression, which means
that existing path integral approximations to the GR rate, such as Wolynes
theory[20, 21], instanton approaches[12, 23] and the linear Golden Rule
method,[14] can be directly applied with this method. Overall these equations
define the Optimal Golden Rule rate (OGR) method.

###  II.4 Modified Marcus and Marcus-Levich-Jortner theories

Figure 4: A–I) HEOM, GR and OGR rate constants as a function of $\beta\Delta$
for a range of spin-boson models with a Brownian oscillator spectral density.
J–L) HEOM, FG and OGR rate constants as a function of $\beta\Delta$ for spin-
boson models with a Debye spectral density. The model parameters are given in
the individual panels. From A–C the reorganisation energy increases with
$\epsilon=2\lambda$. From D–F, G–I, and J–L the rate moves from the
activation-less regime $\epsilon=\lambda$, to deep in the inverted regime, up
to $\epsilon=5\lambda$ for F and $\epsilon=3\lambda$ for I and L. The exact
results in C are taken from Ref. 10.

It is straightforward to apply the optimal golden rule expression method to
spin-boson models, where the diabatic coupling is treated as a constant
$\hat{\Delta}=\Delta$ the diabatic Hamiltonians are treated as a set of
independent Harmonic oscillators

$\displaystyle\hat{V}_{J}=\sum_{\alpha}\frac{1}{2}m_{\alpha}\omega_{\alpha}^{2}\left(\hat{q}_{\alpha}-\frac{c_{J,\alpha}}{m_{\alpha}\omega_{\alpha}^{2}}\right)^{2}+\epsilon_{J},$
(26)

Here $m_{\alpha}$ and $\omega_{\alpha}$ are the oscillator masses and
frequencies, $\epsilon_{J}$ is the energy of $J=\mathrm{A},\mathrm{B}$ and
$c_{J,\alpha}=\delta_{J,\mathrm{B}}c_{\alpha}$ quantifies how much each
oscillator is displaced by going from diabatic state
$\mathrm{A}\to\mathrm{B}$. This is related to the reorganisation energy
$\lambda$ for the harmonic oscillators by
$\lambda=\sum_{\alpha}{c_{\alpha}^{2}}/{2m_{\alpha}\omega_{\alpha}^{2}}$. The
free energy change for the $\mathrm{A}\to\mathrm{B}$ reaction in this model is
simply $\Delta A=\epsilon_{\mathrm{B}}-\epsilon_{\mathrm{A}}$. Note that this
definition differs slightly from that often used in the literature, since only
state B couples to the bath, but this choice has the advantage of making the
definitions of the reorganisation energy in terms of the spectral density (see
below) and the diabatic energy differences equivalent.

Applying the OGR method requires the rotated diabatic Hamiltonians,
$\hat{H}_{\pm}$. Given that the rotated diabatic Hamiltonians $\hat{H}_{\pm}$,
Eqs. (14) and (15), are linear combinations of the original diabats (with an
additional constant shift) it follows that the new oscillator frequencies
after rotation are unchanged and only the coefficients coupling $c_{\alpha}$
are altered. From this it follows that the rotated harmonic potential energy
surfaces are also harmonic with the same distribution of frequencies, but with
a re-normalised reorganisation energy $\tilde{\lambda}_{\theta}$ and free
energy change $\Delta\tilde{A}_{\theta}$, which after some simple algebra are
found to be

$\displaystyle\tilde{\lambda}_{\theta}$
$\displaystyle=\lambda\cos^{2}(2\theta)$ (27)
$\displaystyle\Delta\tilde{A}_{\theta}$
$\displaystyle=\Delta{A}\cos(2\theta)-2\Delta\sin(2\theta).$ (28)

Deep in the Marcus inverted regime, the optimised value of $\theta$ is
$\theta^{*}=\theta_{\mathrm{A}}=(1/2)\atan(2\Delta/|\lambda+\Delta A|)$. The
Marcus theory rate[55, 56, 57] is obtained from the high-temperature limit of
the exact spin-boson rate constant expression. This means that the modified
Marcus theory rate is simply given by

$\displaystyle
k_{\mathrm{MT},\theta^{*}}=\frac{2\pi\Delta^{2}\sec^{2}(2\theta^{*})}{\hbar\sqrt{{4\pi\tilde{\lambda}_{\theta^{*}}k_{\mathrm{B}}T}}}\exp(-\frac{(\tilde{\lambda}_{\theta^{*}}+\Delta\tilde{A}_{\theta^{*}})^{2}}{4\tilde{\lambda}_{\theta^{*}}k_{\mathrm{B}}T}).$
(29)

So overall the modified Marcus expression is simply the standard Marcus theory
rate expression with modified coupling, $\Delta\to\Delta\sec(2\theta^{*})$,
reorganisation energy, $\lambda\to\tilde{\lambda}_{\theta^{*}}$, and free
energy change, $\Delta A\to\Delta\tilde{A}_{\theta^{*}}$. To lowest order in
the coupling $\Delta$, this can be approximated as

$\displaystyle k_{\mathrm{MT},\theta^{*}}\simeq
k_{\mathrm{MT},\theta=0}\exp(-\frac{2\Delta^{2}}{|\lambda+\Delta
A|k_{\mathrm{B}}T}),$ (30)

where $k_{\mathrm{MT},\theta=0}$ is the standard Marcus theory rate constant.
This reveals that the first-order effect of strong coupling is to raise the
activation barrier, and thus suppress the rate constant.

In Marcus-Levich-Jortner (MLJ) theory the harmonic bath is assumed to consist
of a low frequency part, with reorganisation energy $\lambda_{\mathrm{o}}$,
which can be treated with a high temperature approximation and a single high-
frequency mode, with frequency $\Omega$ and reorganisation energy
$\lambda_{\mathrm{v}}$, which is treated with a low temperature approximation.
Exactly the same renormalisation of the reorganisation energy and free energy
change can be applied in this case, yielding the following modified MLJ theory
rate constant,

$\displaystyle
k_{\mathrm{MLJ},\theta^{*}}\\!=\\!\frac{2\pi\Delta^{2}\sec^{2}(2\theta^{*})}{\hbar\sqrt{{4\pi\tilde{\lambda}_{\mathrm{o},\theta^{*}}k_{\mathrm{B}}T}}}\\!\sum_{\nu=0}^{\infty}\\!\left(e^{-\tilde{S}_{\theta^{*}}}\\!\frac{\tilde{S}_{\theta^{*}}^{\nu}}{\nu!}\\!\right)\exp[-\frac{(\tilde{\lambda}_{\mathrm{o},\theta^{*}}+\Delta\tilde{A}_{\theta^{*}}+\hbar\Omega\nu)^{2}}{4\tilde{\lambda}_{\mathrm{o},\theta^{*}}k_{\mathrm{B}}T}]$
(31)

where
$\tilde{S}_{\theta^{*}}=\tilde{\lambda}_{\mathrm{v},\theta^{*}}/\hbar\Omega$
is the renormalised Huang-Rhys factor. Qualitatively one can see that the
reduction in the Huang-Rhys factor on renormalisation of the diabats by strong
coupling leads to a reduction in tunnelling and therefore a reduction in the
rate constant.

Whilst these approximate theories may be qualitatively useful, it is shown in
Appendix C that the modified Marcus theory is not quantitatively accurate
compared to spin-boson models, although it does generally predict the value of
$\Delta$ at which turnover occurs fairly accurately. This is due to the
neglect of nuclear tunnelling effects in Marcus theory, and as such it is
advisable to use the full quantum theory for the spin-boson rate constant with
a full spectral density for processes deep in the inverted regime with strong
electronic coupling.

##  III Numerical tests

Having presented a theory of strong coupling effects in the Marcus inverted
regime in the previous, I now move on to demonstrate the accuracy of this
method on a series of model systems, including an anharmonic condensed phase
system, for which numerically exact rates can be obtained using the
Hierarchical Equations of Motion (HEOM) method.

###  III.1 Spin-boson models

The first set of model problems considered are spin-boson models, where the
model diabatic potentials are given by Eq. (26). The dynamics of the
electronic variables (and therefore also the rate of transfer from A to B), is
dictated fully by the spectral density,

$\displaystyle\mathcal{J}(\omega)=\frac{\pi}{2}\sum_{\alpha}\frac{c_{\alpha}^{2}}{m_{\alpha}\omega_{\alpha}}\delta(\omega-\omega_{\alpha}).$
(32)

The reorganisation energy $\lambda$ is related to the spectral density by
$\lambda=(1/\pi)\int_{0}^{\infty}\mathcal{J}(\omega)/\omega\differential{\omega}$.
Both Brownian oscillator and Debye bath models have been considered, with
spectral densities given respectively by

$\displaystyle\mathcal{J}_{\mathrm{BO}}(\omega)$
$\displaystyle=\frac{2\lambda\gamma\Omega^{2}\omega}{(\omega^{2}-\Omega^{2})^{2}+\gamma^{2}\omega^{2}},$
(33) $\displaystyle\mathcal{J}_{\mathrm{D}}(\omega)$
$\displaystyle=\frac{2\lambda\omega_{\mathrm{D}}\omega}{\omega_{\mathrm{D}}^{2}+\omega^{2}}.$
(34)

All HEOM calculations (expect those in Fig. 4C, which were taken from Ref. 10)
were performed with the heom-lab code using a Matsubara expansion of the bath
correlation functions with the truncation scheme outlined in Ref 58. with the
low temperature correction scheme outlined in Ref. 59. Further details of the
HEOM calculations are given in Appendix B. In all cases the diabatic coupling
is treated as a constant $\hat{\Delta}=\Delta$ between a value of {0.1} and
around {20} or larger, which corresponds to a range of values typically
encountered for molecular systems, between about {20}cm^-1 and {4000}cm^-1 for
$T=\quantity{298}{K}$. In addition to the spectral density parameters, the
reorganisation energy $\lambda$ and the free energy change
$\epsilon=\epsilon_{\mathrm{A}}-\epsilon_{\mathrm{B}}$ are all varied, to
explore a range of models, from the activationless regime, $\epsilon=\lambda$,
to deep in the inverted regime $\epsilon=5\lambda$.

In Fig. 4 the HEOM rate data are plotted together with the GR rate (Eq. (6))
and the OGR rate (Eq. (22)) as a function of $\beta\Delta$. In general we see
that the OGR approximation correctly captures the Golden Rule limit, where
$\beta\Delta$ is small, and the subsequent turnover of the rate at larger
values of $\beta\Delta$. This is robust to the choice of spectral density,
with the OGR approach working for a highly quantum and underdamped case
$\beta\hbar\Omega=4$, $\beta\hbar\gamma=4$, Fig. 4A–F, a less quantum,
overdamped case, $\beta\hbar\Omega=1$, $\beta\hbar\gamma=16$, Fig. 4G–I and
also for a classical Debye spectral density
$\beta\hbar\omega_{\mathrm{D}}=0.5$, Fig. 4J–L.

In the top row of Fig. 4 (A–C) the ratio of $\epsilon$ to $\lambda$ is held
fixed at 2 and $\lambda$ is systematically increased for an underdamped
Brownian Oscillator spectral density, and we see that the the OGR method
remains accurate across the range of reorganisation energies. This accuracy is
retained as ratio of $\epsilon$ to $\lambda$ is varied at fixed $\lambda$ for
this spectral density, going from the Marcus activation-less case
$\epsilon=\lambda$, to deep in the inverted regime $\epsilon=5\lambda$, as can
be seen in the second row of Fig. 4 (D–F). The accuracy of the OGR method also
holds in the overdamped case and for a Debye spectral density as
$\epsilon/\lambda$ is varied, as is shown in the third and fourth rows of Fig.
4 (G–I and J–L) respectively, although the OGR is less accurate at lower
values of $\epsilon/\lambda$ for the overdamped Brownian oscillator.

In all cases we see that for large values of $\beta\Delta$ the GR rate can
overestimate the exact rate by up to four orders of magnitude, whereas the OGR
has a maximum error of a factor of $\sim 4$, which occurs for the overdamped
Brownian oscillator spectral density at intermediate values of the coupling
and low values of $\epsilon/\lambda$. The turnover behaviour of the rate
constant occurs at larger $\beta\Delta$ deeper in the inverted regime, and the
OGR method also appears more accurate in this regime. The small errors at
intermediate values of $\beta\Delta$ may have their origin in recrossing
effects which are neglected in the OGR method. Despite this the overall near-
quantitative accuracy of the OGR method for these spin-boson models is very
encouraging, and OGR rate is never less accurate than the uncorrected GR rate.
A comparison of these results with the modified Marcus theory Eq. (29) and
Zusman theory[52] are shown in Appendix C, where we see that Zusman theory
fails for these spin-boson models and modified Marccus theory (which neglects
nuclear quantum effects) can only qualitively predict the turnover behaviour.

Figure 5: A) The diabatic potential curves for the A and B states used in the
quartic anharmonic model. B) Rate constants calculated from HEOM for the
quartic anharmonic model using GR and OGR methods with
$\theta=\theta_{\mathrm{A}}$ (the proposed optimal solution) and
$\theta=\theta_{\mathrm{B}}$, as well as the full HEOM results.

###  III.2 Anharmonic model

In order to test the robustness of the OGR method to asymmetric frequencies
and anharmonicity in the diabatic potential energy surfaces, exact and OGR
rates were calculated for a model with quartic reaction coordinate, coupled to
a bath to mimic a condensed phase environment. The model potentials are given
by

$\displaystyle\begin{split}V_{J}(Q,{\bf\it
q})&=\frac{1}{2}M\omega_{J}^{2}\bigg{(}(Q-Q_{0,J})^{2}-\alpha_{J}(Q-Q_{0,J})^{3}\\\
&+\frac{7}{12}(Q-Q_{0,J})^{4}\bigg{)}+V_{\mathrm{sb}}(Q,{\bf\it
q})+\epsilon_{J}\end{split}$ (35) $\displaystyle V_{\mathrm{sb}}(Q,{\bf\it
q})$
$\displaystyle=\sum_{\alpha}\frac{1}{2}m_{\alpha}\omega_{\alpha}^{2}\left({q}_{\alpha}-Q\frac{c_{\alpha}}{m_{\alpha}\omega_{\alpha}^{2}}\right)^{2},$
(36)

Full quantum mechanical rate constants and the OGR rates (using the simplified
version Eq. (23)) can be evaluated for this model with HEOM. Working in units
where $\beta=\hbar=M=1$, the model parameters are set to
$\omega_{\mathrm{A}}=2.5$, $\omega_{\mathrm{B}}=4$,
$\alpha_{\mathrm{A}}=\alpha_{\mathrm{B}}=0.25$, $Q_{0,\mathrm{A}}=0$,
$Q_{0,\mathrm{B}}=0.25$, $\epsilon_{\mathrm{A}}=0$,
$\epsilon_{\mathrm{B}}=-20$. The harmonic bath is set as a Debye bath with
$\lambda=0.5$ and $\omega_{\mathrm{D}}=2$. The diabatic potentials are shown
in Fig. 5A. Further details of the calculations are given in Appendix B.

Figure 6: A) The BODIPY (boron-dipyrromethane)-Anthracene molecule, showing
the orbitals involved in charge recombination from the excited
BD^$\bullet-$-An^$\bullet+$ state to the BD-An ground-state. B) Spectral
densities calculated for this transition using molecular dynamics on both the
BD^$\bullet-$-An^$\bullet+$ and BD-An surfaces. Further details of the
spectral density calculations are given in Ref. 9.

The OGR rate constant, which is optimised at $\theta=\theta_{\mathrm{A}}$,
together with the $\theta=\theta_{\mathrm{B}}$ OGR rate, the GR prediction,
and the full HEOM rates are shown in Fig. 5B as a function of $\beta\Delta$.
In this example the exact rate constant agrees with the GR rate up to about
$\beta\Delta=0.7$, which is followed by a small increase in the rate with
respect to the GR result up to $\beta\Delta=2$. Beyond this value of
$\beta\Delta$ the rate decreases with increasing $\Delta$, dropping to about
$10^{-3}k_{\mathrm{GR}}$ at the largest $\Delta$ considered, $\beta\Delta=10$.
The OGR method still yields accurate rate constants in the weak-coupling limit
$\beta\Delta<1$, and it correctly captures the small increase in the rate
constant up to $\beta\Delta=2$ followed by the decrease in rate for larger
$\Delta$. The largest error in this case is a factor of $\sim\\!3$ error at
$\beta\Delta\approx 7$, where the GR rate overestimates the rate constant by a
factor of 40. The fact that the larger value of the OGR rate at
$\theta=\theta_{\mathrm{B}}$ appears to be more accurate in this example at
larger values of $\beta\Delta$ suggests that the simple minimization of the
rate in the range $\theta_{\mathrm{A}}\leq\theta\leq\theta_{\mathrm{B}}$ may
not always be the best method for optimising $\theta$, although it is still a
considerable improvement over the GR rate. Alternative optimization strategies
will be explored in future work.

##  IV Application: Strong-coupling effects in BODIPY-Anthracene charge
recombination

|  $\ce{BD^{\bullet-}-An^{\bullet+}}$ MD |  $\ce{BD-An}$ MD | Average
---|---|---|---
GR | {2.9d7}s^-1 | {4.4d7}s^-1 | {3.6d7}s^-1
OGR | {0.74d7}s^-1 | {1.2d7}s^-1 | {0.95d7}s^-1
Table 1: GR and OGR rate constants for the BODIPY-Anth charge recombination
calculated with the different spectral densities and their average.

Having established that OGR theory provides a good description of non-
adiabatic rates, particularly for spin-boson models deep in the Marcus
inverted regime, an application of the OGR theory to charge recombination in a
BODIPY-anthracene triplet photosensitizer is considered. This molecule forms a
BODIPY (BD) to anthracene (An) charge transfer state after photo-excitation,
which can either undergo intersystem crossing to the triplet manifold or
recombine to the ground electronic state.[60, 7, 9] Previous work has found
that the $\ce{BD^{\bullet-}-An^{\bullet+}->BD-An}$ charge recombination in
acetonitrile process lies deep in the inverted regime, with $-\Delta A\approx
5\lambda$,[60, 7, 9] and furthermore the good overlap between the BD
$\uppi^{*}$ and An $\uppi$ orbitals means the electronic coupling is very
large at $\quantity{1904}{cm^{-1}}\approx 9.2\ k_{\mathrm{B}}T$ (at
$T=\quantity{298}{K}$).[9] The orbitals involved and the BD-An molecule are
shown in Fig. 6A.

In a previous study[9] a spin-boson mapping was constructed for the charge
recombination process based on molecular dynamics simulations using molecular
mechanics force-fields parameterised for the ground and charge-transfer
states. The Fourier transform of the energy gap correlation function yielded a
spectral density, which can be computed either for dynamics on the BD-An or
$\ce{BD^{\bullet-}-An^{\bullet+}}$ potential energy surface. These spectral
distributions, $\rho(\omega)=\mathcal{J}(\omega)/(\omega\pi\lambda)$, are
given in Fig. 6. Further details of these calculations used to parameterise
the spin-boson models are given in Ref. 9. Using these spectral densities, the
rate constant for the charge recombination process can be estimated using
either the standard GR approach or the OGR theory. A comparison with a
previously constructed model spectral density[9] given in Appendix D suggests
that this problem should be in a regime where OGR theory is reliable. The GR
and OGR rates (with full optimisation, Eq. (22)) for both spectral densities
are given in Table 1, where one can see that for both spectral densities the
rate drops by a factor of 4 on including strong diabatic coupling effects.

The charge recombination rate plays an important role in determining the
triplet quantum yield of BD-An, because charge recombination competes with
intersystem crossing from the singlet $\ce{BD^{\bullet-}-An^{\bullet+}}$ to
the triplet excited state. By taking the charge recombination rate as the
average of the OGR rates for the two spectral densities, the triplet quantum
yield $\Phi_{\mathrm{T}}$ for BODIPY-Anthracene can be calculated. Without
strong coupling effects in the charge recombination, $\Phi_{\mathrm{T}}$ was
previously calculated to be 0.8, however including strong coupling effects
increases this to 0.88, much closer to the experimental value of
0.93[60]–0.96.[7] (Details of this calcualtion are given in Appendix E) This
highlights the importance of strong electronic coupling effects in a real
molecular system. Further details of the other photophysical processes and the
triplet yield calculation are given in Ref. 9.

##  V Conclusions

In this paper I have presented an Optimal Golden Rule theory for calculating
rate constants of non-adiabatic processes in the Marcus inverted regime where
coupling between diabatic electronic states becomes large. The theory is based
on rotating to a new set of diabatic states where the strongly coupled states
are partially mixed, in such a way that the coupling between states is
reduced, enabling pertubation theory to be applied more accurately. This
simple approximation works well for a range of tested spin-boson models and an
anharmonic model. OGR theory can be applied straightforwardly with existing
methods for computing GR rates for anharmonic potential energy surfaces, such
as using spin-boson mappings, instanton approaches and path integral methods.
The simplicity of OGR theory also allows one to obtain a generalization of
Marcus theory and MLJ theory to strong electronic couplings in the inverted
regime, where it was found the OGR theory predicts that strong electronic
coupling effectively renormalises the diabatic coupling, free energy change
and reorganisation energy of the reaction.

Evaluation of the OGR rate is equivalent to the evaluation of a GR rate, so
any method that can be applied with GR rates can also in principle be combined
with OGR theory. This will enable OGR theory to be applied to atomistic models
of reactions, beyond spin-boson models and model potentials considered in this
work. The current simple condition for optimising the mixing angle $\theta$
just requires evaluation of expectation values of quantities in the original A
and B diabatic equilibrium states, which can also be done very
straightforwardly for atomistic models including nuclear quantum effects using
path integral methods.[61] Despite its simplicity, there may exist better
$\theta$ optimisation strategies for obtaining accurate rate constants. For
example it is possible that a more direct minimisation of
$\expectationvalue{\Delta_{+-}^{2}}_{\pm}$ or a minimisation of the initial
value of the flux-flux correlation function
${\Tr}[\hat{F}^{2}e^{-\beta\hat{H}_{\pm}}]\propto
Q_{\pm}\expectationvalue{\Delta_{+-}^{2}}_{\pm}$ (as is done in the NAQI
method[25] or transition state theory) may provide more accurate results for
the OGR method. Furthermore there is no reason that the same value of $\theta$
should be used to evaluate both the approximate flux-flux correlation function
and the reactant partition function, and these could be optimised separately,
although early explorations of this idea have not produced better results than
those presented here. Thorough investigation of these alternative approaches
will be the subject of future work. Although there may still be room for
improvement, OGR theory already enables the accurate calculation of non-
adiabatic rate constants at arbitrary diabatic coupling strengths in the
inverted regime, a regime currently inaccessible to other methods such as the
interpolation formula approach,[10] and it shows considerable promise in
application to atomistic models and real-world molecular systems. OGR theory
also offers a route to generalising theory of other non-adiabatic processes,
such as proton-coupled electron transfer,[62] excitation energy transfer,[63]
both of which are relevant in enzymatic and photosynthetic systems, as well as
electrochemical processes.[64, 4]

### Acknowledgements

I would like to thank Joseph Lawrence for useful discussions as well as David
Limmer and the Limmer group for some early conversations on this work. This
work was supported by “Agence Nationale de la Recherche” through the project
MAPPLE (ANR-22-CE29-0014-01).

### Data availability

All data is provided within the article itself. Code used to perform the HEOM
calculations is available at https://github.com/tomfay/heom-lab and scripts
are available from the corresponding author upon request.

### Conflicts of interest

The author declares no conflicts of interest.

##  Appendix A Derivation of simplified $k_{\mathrm{GR},\theta}$

The coupling between $\ket{+}$ and $\ket{-}$ states involves a term
proportional to $\hat{U}=\hat{H}_{+}-\hat{H}_{-}$. In the static approximation
or classical limit, where the kinetic energy and potential energy operators
commute, the time integral of the propagators becomes,[53]

$\displaystyle\int_{-\infty}^{+\infty}\differential{t}e^{+i\hat{H}_{-}t/\hbar}e^{-i\hat{H}_{-}t/\hbar}\approx
2\pi\hbar\delta(\hat{U}),$ (37)

and therefore it is clear that these terms in the rate constant proportional
to $\hat{U}^{n}$ vanish. We can also show this more generally. First
considering the terms where a single $\hat{U}$ operator appears, we note that
the integrand can be written as a derivative as follows

$\displaystyle\int_{-\infty}^{\infty}$
$\displaystyle\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{+}-i\hat{H}_{+}t/\hbar}\hat{\Delta}e^{+i\hat{H}_{-}t/\hbar}\hat{U}]\differential{t}$
$\displaystyle=i\hbar\int_{-\infty}^{\infty}\derivative{t}\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{+}-i\hat{H}_{+}t/\hbar}\hat{\Delta}e^{+i\hat{H}_{-}t/\hbar}]\differential{t}$
(38)
$\displaystyle=i\hbar\left[\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{+}-i\hat{H}_{+}t/\hbar}\hat{\Delta}e^{+i\hat{H}_{-}t/\hbar}]\right]_{t=-\infty}^{t=\infty}$
(39)

and since this correlation function decays to zero, these contributions are
zero. The same argument applies to the terms where $\hat{\Delta}$ and
$\hat{U}$ are swapped. Applying the same argument twice we find that the terms
involving two factors of $\hat{U}$ are also zero, because

$\displaystyle\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{+}-i\hat{H}_{+}t/\hbar}\hat{U}e^{+i\hat{H}_{-}t/\hbar}\hat{U}]$
$\displaystyle=(i\hbar^{2})\derivative[2]{t}\Tr_{\mathrm{n}}[e^{-\beta\hat{H}_{+}-i\hat{H}_{+}t/\hbar}e^{+i\hat{H}_{-}t/\hbar}].$
(40)

Overall this shows that Eq. (17) simplifies to Eq. (25).

##  Appendix B HEOM calculations

Numerically exact rate constants for the models considered are obtained using
the hierarchical equations of motion (HEOM).[43, 44] The heom-lab code was
used to perform all simulations, using a Matsurbara expansion for constructing
the hierarchy with the low temperature correction from Ref. 59. The HEOM
truncation scheme from Ref. 58 used in all calculations. Up to $M=3$ Matsubara
modes were included explicitly in the calculations with the depth of the
hierarchy set to between $L=40$ and $L=150$ depending on the spectral density
used.

The rate constant was extracted by initialising the system in the
$\hat{\rho}(0)=\outerproduct{\mathrm{A}}{\mathrm{A}}e^{-\beta\hat{H}_{\mathrm{A}}}/Q_{\mathrm{A}}$
state and fitting the population dynamics of state B to
$p_{\mathrm{B}}(t)=(p_{\mathrm{B},\infty}-p_{\mathrm{B},0})(1-e^{-k_{\mathrm{obs}}t})+p_{\mathrm{B},0}$
after initial transient oscillations had decayed. Strictly speaking
$k_{\mathrm{obs}}$ is the sum of forward and backward rates, but in all
examples considered the ratio of the backward to forward rates is $\lesssim
10^{-9}$ so the rate constant is taken to be the fitted value
$k_{\mathrm{obs}}$.

For the anharmonic model the system Hamiltonian was constructed in a Colbert-
Miller DVR basis[65] with 500 basis functions for each diabatic state and then
diagonalised. The $N_{\mathrm{A}}=11$ lowest energy states for diabat A and
$N_{\mathrm{B}}=15$ lowest energy states for diabat B were used to construct
the full system Hamiltonian for the exact rate calculations and the $\Delta=0$
system Hamiltonian for the GR calculations with HEOM.

Figure 7: A comparison of HEOM spin-boson model rates together with OGR and GR
rates and the modified Marcus and Zusman theory rates. All model parameters
are the same as in Fig. 4.

##  Appendix C Comparison of Modified Marcus and Zusman theories

Here a comparison of modified Marcus theory and Zusman theory are shown for
the spin-boson models considered in Sec. III.1. Zusman theory gives the rate
constant as[52, 53, 54]

$\displaystyle
k_{\mathrm{Zus}}=\frac{k_{\mathrm{MT,f}}}{1+\tau_{\mathrm{S,f}}k_{\mathrm{MT,f}}+\tau_{\mathrm{S,b}}k_{\mathrm{MT,b}}}$
(41)

where $k_{\mathrm{MT,f/b}}$ are the forward/backward Marcus theory rate
constants, and $\tau_{\mathrm{S,f/b}}$ are effective solvent relaxation times
given by,

$\displaystyle\begin{split}&\tau_{\mathrm{S,f/b}}=\tau_{\mathrm{S}}\bigg{(}\ln(2)\\\
&+\frac{(\lambda\pm\Delta
A)^{2}}{2k_{\mathrm{B}}T\lambda}{}_{2}F_{2}\left(1,1;\frac{2}{2},2;\frac{(\lambda\pm\Delta
A)^{2}}{4k_{\mathrm{B}}T\lambda}\right)\bigg{)}\end{split}$ (42)

where ${}_{2}F_{2}(a_{1},a_{2};b_{1},b_{2};z)$ is a generalised hypergeometric
function, and $\tau_{\mathrm{S}}$ is the relaxation time of the energy gap
coordinate. Qualitatively Zusman theory predicts that the rate constant
plateaus at large values of $\Delta$, which is clearly not what is observed in
the spin-boson models in Sec. III.1 which show turnover behaviour. For the
Debye spectral density $\tau_{\mathrm{S}}=1/\omega_{\mathrm{D}}$, for the
underdamped Brownian Oscillator $\tau_{\mathrm{S}}=2/\gamma$ and for the
overdamped Brownian Oscillator $\tau_{\mathrm{S}}=\gamma/\Omega^{2}$.

In Fig. 7 we show the spin-boson model results from Fig. 4 together with the
Modified Marcus and Zusman results (note that the modified Marcus rate is not
defined for the activation-less cases: Fig. 4D, G and J). We see that neither
is quantitatively accurate; modified Marcus theory generally over-predicts the
drop in the rate constant at large values of $\Delta$ whilst Zusman theory
simply plateaus at large $\Delta$. Furthermore, since both are classical
theories the neglect of nuclear quantum effects leads to large errors even in
the Golden Rule regime. We see that modified Marcus theory at least roughly
predicts the value of $\Delta$ at which rate turnover occurs relatively
accurately, so it should at least be useful as a guide for when GR theory is
likely to be accurate. The importance of nuclear quantum effects in the
inverted regime with large diabatic couplings clearly means that accounting
for the details of the full spectral density are essential for quantitatively
accurate rate calculations in this regime.

##  Appendix D OGR with the BD-An model spectral density

In Ref. 9 an approximate post-GR correction to the BD-An charge recombination
rate was calculated using a model spectral density, which is a sum of a Debye
spectral density (with $\beta\lambda=10.1459$ and
$\beta\hbar\omega_{\mathrm{D}}=0.1831$) and an underdamped Brownian Oscillator
spectral density (with $\beta\lambda=8.6780$, $\beta\hbar\Omega=6.76$ and
$\beta\hbar\gamma=4$). The model spectral distribution is shown in Fig. 8A.
This model predicted a much smaller reduction in the rate, by only about a
factor of 0.9, relative to the GR rate, comapred to the OGR calculations with
the full spectral density performed here. However this simplified model
overestimates the GR rate compared to the full spectral density by over a
factor of 1000. This is because of the polynomial decaying tails of these
spectral densities which leads to very large tunnelling effects. Despite this
issue with the model, we can use the HEOM results to benchmark OGR theory.
This comparison is shown in Fig. 8B, where we see that OGR theory is very
accurate for this model (note that for BD-An $\beta\Delta\approx 10$). This
suggests that OGR theory should be reliable for the full spetral density as
well.

Figure 8: A) The model spectral distribution for charge recombination in BD-
An. B) A comparison of GR, HEOM and OGR rates with this model as a function of
$\Delta$.

##  Appendix E Triplet yield calculation for BODIPY-Anthracene

The BODIPY-Anthracene singlet charge transfer state can either recombine to
the ground singlet electronic state, or undergo intersystem crossing to yield
the excited triplet state, which can be used as a photosensitizer. This charge
transfer state is in pre-equilibrium with the BD^*-An valence excited state,
so the overall triplet yield $\Phi_{\mathrm{T}}$ is

$\displaystyle\Phi_{\mathrm{T}}=\frac{p_{\mathrm{CT}}k_{\mathrm{ISC,tot}}}{p_{\mathrm{CT}}k_{\mathrm{ISC,tot}}+p_{\mathrm{CT}}k_{\mathrm{CR}}+k_{\mathrm{Fl,tot}}}.$
(43)

${k_{\mathrm{ISC,tot}}}$ is the sum of intersystem crossing rates from the
charge transfer (CT) state, $k_{\mathrm{CR}}$ is the charge recombination rate
back to the ground state (which is recalculated in this work),
$k_{\mathrm{Fl,tot}}$ is the total effective fluorescence rate from both the
CT and BD^*-An states, and $p_{\mathrm{CT}}=e^{-\Delta
A_{\mathrm{CS}}/k_{\mathrm{B}}T}/(1+e^{-\Delta
A_{\mathrm{CS}}/k_{\mathrm{B}}T})$ is the pre-equilibrium population of the CT
state ($\Delta A_{\mathrm{CS}}$ is the free energy change of excited state
charge separation). From Ref. 9
$k_{\mathrm{ISC,tot}}=\quantity{2.05e8}{s^{-1}}$,
$k_{\mathrm{Fl,tot}}=\quantity{1.74e7}{s^{-1}}$ and $p_{\mathrm{CT}}=0.903$.

## References

## References

  * Blumberger [2015] J. Blumberger, “Recent Advances in the Theory and Molecular Simulation of Biological Electron Transfer Reactions,” Chemical Reviews 115, 11191–11238 (2015).
  * Oberhofer, Reuter, and Blumberger [2017] H. Oberhofer, K. Reuter, and J. Blumberger, “Charge Transport in Molecular Materials: An Assessment of Computational Methods,” Chemical Reviews 117, 10319–10357 (2017).
  * Marcus and Sutin [1985] R. Marcus and N. Sutin, “Electron transfers in chemistry and biology,” Biochimica et Biophysica Acta (BBA) - Reviews on Bioenergetics 811, 265–322 (1985).
  * Santos and Schmickler [2022] E. Santos and W. Schmickler, “Models of Electron Transfer at Different Electrode Materials,” Chemical Reviews 122, 10581–10598 (2022).
  * Richter and Fingerhut [2020] M. Richter and B. P. Fingerhut, “Regulatory Impact of the C-Terminal Tail on Charge Transfer Pathways in Drosophila Cryptochrome,” Molecules 25, 4810 (2020).
  * Giannini  _et al._ [2019] S. Giannini, A. Carof, M. Ellis, H. Yang, O. G. Ziogos, S. Ghosh, and J. Blumberger, “Quantum localization and delocalization of charge carriers in organic semiconducting crystals,” Nature Communications 10, 3843 (2019).
  * Wang  _et al._ [2019] Z. Wang, A. A. Sukhanov, A. Toffoletti, F. Sadiq, J. Zhao, A. Barbon, V. K. Voronkova, and B. Dick, “Insights into the Efficient Intersystem Crossing of Bodipy-Anthracene Compact Dyads with Steady-State and Time-Resolved Optical/Magnetic Spectroscopies and Observation of the Delayed Fluorescence,” The Journal of Physical Chemistry C 123, 265–274 (2019).
  * Kubas  _et al._ [2014] A. Kubas, F. Hoffmann, A. Heck, H. Oberhofer, M. Elstner, and J. Blumberger, “Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT, and FODFTB against high-level ab initio calculations,” The Journal of Chemical Physics 140, 104105 (2014).
  * Fay and Limmer [2024] T. P. Fay and D. T. Limmer, “Unraveling the mechanisms of triplet state formation in a heavy-atom free photosensitizer,” Chemical Science , 6726 (2024).
  * Lawrence  _et al._ [2019] J. E. Lawrence, T. Fletcher, L. P. Lindoy, and D. E. Manolopoulos, “On the calculation of quantum mechanical electron transfer rates,” The Journal of Chemical Physics 151, 114119 (2019).
  * Thapa, Fang, and Richardson [2019] M. J. Thapa, W. Fang, and J. O. Richardson, “Nonadiabatic quantum transition-state theory in the golden-rule limit. I. Theory and application to model systems,” The Journal of Chemical Physics 150, 104107 (2019).
  * Heller and Richardson [2020a] E. R. Heller and J. O. Richardson, “Instanton formulation of Fermi’s golden rule in the Marcus inverted regime,” The Journal of Chemical Physics 152, 034106 (2020a).
  * Heller and Richardson [2020b] E. R. Heller and J. O. Richardson, “Semiclassical instanton formulation of Marcus–Levich–Jortner theory,” The Journal of Chemical Physics 152, 244117 (2020b).
  * Lawrence and Manolopoulos [2020a] J. E. Lawrence and D. E. Manolopoulos, “An improved path-integral method for golden-rule rates,” The Journal of Chemical Physics 153, 154113 (2020a).
  * Lawrence and Manolopoulos [2020b] J. E. Lawrence and D. E. Manolopoulos, “Confirming the role of nuclear tunneling in aqueous ferrous–ferric electron transfer,” The Journal of Chemical Physics 153, 154114 (2020b).
  * Ananth [2022] N. Ananth, “Path Integrals for Nonadiabatic Dynamics: Multistate Ring Polymer Molecular Dynamics,” Annual Review of Physical Chemistry 73, 299–322 (2022).
  * Duke and Ananth [2016] J. R. Duke and N. Ananth, “Mean field ring polymer molecular dynamics for electronically nonadiabatic reaction rates,” Faraday Discussions 195, 253–268 (2016).
  * Ranya and Ananth [2020] S. Ranya and N. Ananth, “Multistate ring polymer instantons and nonadiabatic reaction rates,” The Journal of Chemical Physics 152, 114112 (2020).
  * Shi and Geva [2004] Q. Shi and E. Geva, “Nonradiative Electronic Relaxation Rate Constants from Approximations Based on Linearizing the Path-Integral Forward-Backward Action,” The Journal of Physical Chemistry A 108, 6109–6116 (2004).
  * Wolynes [1987] P. G. Wolynes, “Imaginary time path integral Monte Carlo route to rate coefficients for nonadiabatic barrier crossing,” The Journal of Chemical Physics 87, 6559–6561 (1987).
  * Lawrence and Manolopoulos [2018] J. E. Lawrence and D. E. Manolopoulos, “Analytic continuation of Wolynes theory into the Marcus inverted regime,” The Journal of Chemical Physics 148, 102313 (2018).
  * Richardson [2015] J. O. Richardson, “Ring-polymer instanton theory of electron transfer in the nonadiabatic limit,” The Journal of Chemical Physics 143, 134116 (2015).
  * Trenins and Richardson [2022] G. Trenins and J. O. Richardson, “Nonadiabatic instanton rate theory beyond the golden-rule limit,” The Journal of Chemical Physics 156, 174115 (2022).
  * Ansari  _et al._ [2022] I. M. Ansari, E. R. Heller, G. Trenins, and J. O. Richardson, “Instanton theory for Fermi’s golden rule and beyond,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380, 20200378 (2022).
  * Lawrence and Manolopoulos [2020c] J. E. Lawrence and D. E. Manolopoulos, “A general non-adiabatic quantum instanton approximation,” The Journal of Chemical Physics 152, 204117 (2020c).
  * Lawrence and Manolopoulos [2020d] J. E. Lawrence and D. E. Manolopoulos, “Path integral methods for reaction rates in complex systems,” Faraday Discussions 221, 9–29 (2020d).
  * Craig and Manolopoulos [2005] I. R. Craig and D. E. Manolopoulos, “Chemical reaction rates from ring polymer molecular dynamics,” The Journal of Chemical Physics 122, 084106 (2005).
  * Geva, Shi, and Voth [2001] E. Geva, Q. Shi, and G. A. Voth, “Quantum-mechanical reaction rate constants from centroid molecular dynamics simulations,” The Journal of Chemical Physics 115, 9209–9222 (2001).
  * Miller [1975] W. H. Miller, “Semiclassical limit of quantum mechanical transition state theory for nonseparable systems,” The Journal of Chemical Physics 62, 1899–1906 (1975).
  * Richardson [2016] J. O. Richardson, “Derivation of instanton rate theory from first principles,” The Journal of Chemical Physics 144, 114106 (2016).
  * Mannouch and Richardson [2023] J. R. Mannouch and J. O. Richardson, “A mapping approach to surface hopping,” The Journal of Chemical Physics 158, 104111 (2023).
  * Lawrence, Mannouch, and Richardson [2024] J. E. Lawrence, J. R. Mannouch, and J. O. Richardson, “Recovering marcus theory rates and beyond without the need for decoherence corrections: The mapping approach to surface hopping,” The Journal of Physical Chemistry Letters 15, 707–716 (2024).
  * Runeson and Manolopoulos [2023] J. E. Runeson and D. E. Manolopoulos, “A multi-state mapping approach to surface hopping,” The Journal of Chemical Physics 159, 094115 (2023).
  * Runeson, Fay, and Manolopoulos [2024] J. E. Runeson, T. P. Fay, and D. E. Manolopoulos, “Exciton dynamics from the mapping approach to surface hopping: comparison with Förster and Redfield theories,” Physical Chemistry Chemical Physics 26, 4929–4938 (2024).
  * Subotnik  _et al._ [2016] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, “Understanding the Surface Hopping View of Electronic Transitions and Decoherence,” Annual Review of Physical Chemistry 67, 387–417 (2016).
  * Jain and Subotnik [2015] A. Jain and J. E. Subotnik, “Surface hopping, transition state theory, and decoherence. II. Thermal rate constants and detailed balance,” The Journal of Chemical Physics 143, 134107 (2015).
  * Kelly, Brackbill, and Markland [2015] A. Kelly, N. Brackbill, and T. E. Markland, “Accurate nonadiabatic quantum dynamics on the cheap: Making the most of mean field theory with master equations,” The Journal of Chemical Physics 142, 094110 (2015).
  * Pfalzgraff, Kelly, and Markland [2015] W. C. Pfalzgraff, A. Kelly, and T. E. Markland, “Nonadiabatic Dynamics in Atomistic Environments: Harnessing Quantum-Classical Theory with Generalized Quantum Master Equations,” The Journal of Physical Chemistry Letters 6, 4743–4748 (2015).
  * Amati  _et al._ [2022] G. Amati, M. A. C. Saller, A. Kelly, and J. O. Richardson, “Quasiclassical approaches to the generalized quantum master equation,” The Journal of Chemical Physics 157, 234103 (2022).
  * Mulvihill and Geva [2021] E. Mulvihill and E. Geva, “A Road Map to Various Pathways for Calculating the Memory Kernel of the Generalized Quantum Master Equation,” The Journal of Physical Chemistry B 125, 9834–9852 (2021).
  * Wang and Thoss [2003] H. Wang and M. Thoss, “Multilayer formulation of the multiconfiguration time-dependent Hartree theory,” The Journal of Chemical Physics 119, 1289–1299 (2003).
  * Wang, Skinner, and Thoss [2006] H. Wang, D. E. Skinner, and M. Thoss, “Calculation of reactive flux correlation functions for systems in a condensed phase environment: A multilayer multiconfiguration time-dependent Hartree approach,” The Journal of Chemical Physics 125, 174502 (2006).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, “Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath,” Journal of the Physical Society of Japan 58, 101–114 (1989).
  * Tanimura [2020] Y. Tanimura, “Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM),” The Journal of Chemical Physics 153, 020901 (2020).
  * Makri [2014] N. Makri, “Blip decomposition of the path integral: Exponential acceleration of real-time calculations on quantum dissipative systems,” The Journal of Chemical Physics 141, 134117 (2014).
  * Topaler and Makri [1996] M. Topaler and N. Makri, “Path Integral Calculation of Quantum Nonadiabatic Rates in Model Condensed Phase Reactions,” The Journal of Physical Chemistry 100, 4430–4436 (1996).
  * Ziogos and Blumberger [2021] O. G. Ziogos and J. Blumberger, “Ultrafast estimation of electronic couplings for electron transfer between pi-conjugated organic molecules. II,” The Journal of Chemical Physics 155, 244110 (2021).
  * Chandler [2009] D. Chandler,  _Introduction to modern statistical mechanics_ (Oxford Univ. Press, New York, NY, 2009).
  * Richardson and Thoss [2014] J. O. Richardson and M. Thoss, “Non-oscillatory flux correlation functions for efficient nonadiabatic rate theory,” The Journal of Chemical Physics 141, 074106 (2014).
  * Evans and Coalson [1995] D. G. Evans and R. D. Coalson, “Incorporating backflow into a relaxation theory treatment of the dynamics of nonequilibrium nonadiabatic transition processes,” The Journal of Chemical Physics 102, 5658–5668 (1995).
  * Golosov and Reichman [2001] A. A. Golosov and D. R. Reichman, “Reference system master equation approaches to condensed phase charge transfer processes. I. General formulation,” The Journal of Chemical Physics 115, 9848–9861 (2001).
  * Zusman [1988] L. Zusman, “The theory of electron transfer reactions in solvents with two characteristic relaxation times,” Chemical Physics 119, 51–61 (1988).
  * Sparpaglione and Mukamel [1988a] M. Sparpaglione and S. Mukamel, “Dielectric friction and the transition from adiabatic to nonadiabatic electron transfer in condensed phases. II. Application to non-Debye solvents,” The Journal of Chemical Physics 88, 4300–4311 (1988a).
  * Sparpaglione and Mukamel [1988b] M. Sparpaglione and S. Mukamel, “Dielectric friction and the transition from adiabatic to nonadiabatic electron transfer. I. Solvation dynamics in Liouville space,” The Journal of Chemical Physics 88, 3263–3280 (1988b).
  * Hush [1958] N. S. Hush, “Adiabatic Rate Processes at Electrodes. I. Energy-Charge Relationships,” The Journal of Chemical Physics 28, 962–972 (1958).
  * Kubo and Toyozawa [1955] R. Kubo and Y. Toyozawa, “Application of the Method of Generating Function to Radiative and Non-Radiative Transitions of a Trapped Electron in a Crystal,” Progress of Theoretical Physics 13, 160–182 (1955).
  * Marcus [1956] R. A. Marcus, “On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I,” The Journal of Chemical Physics 24, 966–978 (1956).
  * Lindoy, Mandal, and Reichman [2023] L. P. Lindoy, A. Mandal, and D. R. Reichman, “Quantum dynamical effects of vibrational strong coupling in chemical reactivity,” Nature Communications 14, 2733 (2023).
  * Fay [2022] T. P. Fay, “A simple improved low temperature correction for the hierarchical equations of motion,” The Journal of Chemical Physics 157, 054108 (2022).
  * Buck  _et al._ [2019] J. T. Buck, A. M. Boudreau, A. DeCarmine, R. W. Wilson, J. Hampsey, and T. Mani, “Spin-Allowed Transitions Control the Formation of Triplet Excited States in Orthogonal Donor-Acceptor Dyads,” Chem 5, 138–155 (2019).
  * Markland and Ceriotti [2018] T. E. Markland and M. Ceriotti, “Nuclear quantum effects enter the mainstream,” Nature Reviews Chemistry 2, 0109 (2018).
  * Hammes-Schiffer and Stuchebrukhov [2010] S. Hammes-Schiffer and A. A. Stuchebrukhov, “Theory of Coupled Electron and Proton Transfer Reactions,” Chemical Reviews 110, 6939–6960 (2010).
  * Renger [2009] T. Renger, “Theory of excitation energy transfer: from structure to function,” Photosynthesis Research 102, 471–485 (2009).
  * Marcus [1965] R. A. Marcus, “On the Theory of Electron-Transfer Reactions. VI. Unified Treatment for Homogeneous and Electrode Reactions,” The Journal of Chemical Physics 43, 679–701 (1965).
  * Colbert and Miller [1992] D. T. Colbert and W. H. Miller, “A novel discrete variable representation for quantum mechanical reactive scattering via the S -matrix Kohn method,” The Journal of Chemical Physics 96, 1982–1991 (1992).