# Nonadiabatic instanton rate theory beyond the golden-rule limit

George Trenins<EMAIL_ADDRESS>Jeremy O. Richardson
<EMAIL_ADDRESS>Laboratory of Physical Chemistry, ETH
Zürich, 8093 Zürich, Switzerland

###### Abstract

Fermi’s golden rule describes the leading-order behaviour of the reaction rate
as a function of the diabatic coupling. Its asymptotic (  $\hbar\to 0$ ) limit
is the semiclassical golden-rule instanton rate theory, which rigorously
approximates nuclear quantum effects, lends itself to efficient numerical
computation and gives physical insight into reaction mechanisms. However the
golden rule by itself becomes insufficient as the strength of the diabatic
coupling increases, so higher-order terms must be additionally considered. In
this work we give a first-principles derivation of the next-order term beyond
the golden rule, represented as a sum of three components. Two of them lead to
new instanton pathways that extend the golden-rule case and, among other
factors, account for the effects of recrossing on the full rate. The remaining
component derives from the equilibrium partition function and accounts for
changes in potential energy around the reactant and product wells due to
diabatic coupling. The new semiclassical theory demands little computational
effort beyond a golden-rule instanton calculation. It makes it possible to
rigorously assess the accuracy of the golden-rule approximation and sets the
stage for future work on general semiclassical nonadiabatic rate theories.

††preprint: AIP/123-QED

##  I Introduction

Semiclassical instanton rate theory[1, 2, 3, 4, 5, 6, 7, 8] is becoming a
well-established method for describing chemical reactions on a single
Born–Oppenheimer potential energy surface (PES).[9, 10, 11, 12, 13, 14, 15,
16] The theory accounts for nuclear quantum effects such as tunnelling in an
approximate but rigorous manner and scales favourably with system size. Hence
it can be readily applied to full-dimensional _ab initio_ simulations of
chemical transformations. The favourable scaling arises because (in the
simplest case) the instanton rate derives from a single classical trajectory.
This not only simplifies the calculation compared to the full quantum
treatment, but also offers direct mechanistic insight, as this trajectory
defines the dominant tunnelling pathway, which can easily be visualised.

However, a description in terms of a single Born–Oppenheimer PES is only valid
for systems with strong diabatic coupling $\hat{\Delta}$. Away from this limit
the Born–Oppenheimer approximation breaks down, and the chemical process is
said to be nonadiabatic. For weak diabatic coupling, nonadiabaticity can be
effectively described with perturbation theory, taking the diabatic
Hamiltonian[17] as the reference. In the case of reaction rates, the leading-
order perturbation term (second-order in $\hat{\Delta}$) is the well-known
Fermi’s golden rule (GR).[18, 19] The rule can be formally expressed in terms
of path integrals, whose evaluation by steepest descent leads to semiclassical
golden-rule instanton rate theory.[20, 21, 22, 23, 24] Like its
Born–Oppenheimer counterpart, this semiclassical theory rigorously
approximates the quantum rate using information from (typically) a single
classical trajectory and hence enjoys the same computational advantages. The
theory has been extended[25, 24] into the Marcus inverted regime[26, 27]—often
a considerable methodological challenge for path-integral methods—and has been
used in _ab initio_ calculations of spin crossover rates for thiophosgene,
showing unprecedented accuracy.[28]

That said, Fermi’s golden rule is only accurate if the next term in the
perturbation expansion of the full rate is small by comparison. In other
words, the GR expression should only be used if the diabatic coupling is
sufficiently weak. This assumption is often true for electron-transfer[29, 27,
30, 31] and spin-crossover[32, 33] reactions, however its validity is by no
means guaranteed. One expects to also find systems with intermediate values of
the diabatic coupling, lying in a “grey area” where neither Born–Oppenheimer
nor golden-rule rate theory are valid.[34, 35, 36, 37, 38] This intermediate
regime has always been of great interest, and the search for a corresponding
practical nonadiabatic rate theory is ongoing.[39]

Considerable progress [40, 41, 42, 43, 44] has been made for dissipative
systems that can be mapped onto a spin–boson model,[45] of which electron
transfer in solution is a typical example.[46] Our goal, however, is to
develop an _ab initio_ theory of molecular reactions, for which instanton
methods are particularly well suited. Such reactions usually take place in the
gas phase, so one need not consider solvent effects, particularly the
possibility that the reaction is dominated by diffusive motion along the
reaction coordinate. Instead the key requirement for our theory is that it be
applicable to multidimensional anharmonic potentials.

Early work that has paved the way for such developments includes Landau–Zener
(LZ) theory,[47, 48] which gives the correct classical (high-temperature) rate
constant for a one-dimensional linear system in the GR limit[49, 50] and which
has been used in uniform rate expressions for reduced models of electron
transfer in solution.[41, 42, 46] Zhu–Nakamura theory[51, 52, 53, 54] is a
generalisation of LZ that accounts for nuclear quantum effects and hence
yields (an asymptotic approximation to) the reaction rate, valid at any
temperature or strength of diabatic coupling. Both theories are, however,
limited to one-dimensional systems, and are only strictly rigorous for linear
potentials.[50, 54] Their other limitation is conceptual, as neither provides
the same mechanistic insight as semiclassical instantons, since no stationary-
action pathways are computed.

An alternative approach was proposed in a recent publication by Lawrence  _et
al._ ,[55] who introduce an interpolation formula connecting the golden-rule
and the adiabatic (Born–Oppenheimer) limits. The formula is inspired by the
Zusman equation, originally derived for the rate of electron transfer in the
classical, high-friction limit.[40, 56, 44] In contrast to the Zusman
equation, the interpolation formula can account for nuclear tunnelling and
zero-point energy effects, as it is a function of a golden-rule and an
adiabatic rate that by construction inherits “quantumness” from its input
(originally ring-polymer molecular dynamics[57, 58, 59] and Wolynes
theory,[60] but could equally be the closely related instanton rate
theory[4]). Unlike Landau–Zener and Zhu–Nakamura theories, the interpolation
formula is obtained following a more heuristic approach, wherein lie both its
strength and its weakness. The formula is readily applicable to
multidimensional problems, and its accuracy has been convincingly demonstrated
for the spin–boson model with a Brownian oscillator density. However, since
the inspiration for the interpolation formula—the Zusman equation—is specific
to the spin–boson model, it is unclear how accurately this equation (and hence
the interpolation formula) describe general anharmonic potentials, or how
accurate the approach will prove for systems whose golden-rule and adiabatic
reaction mechanisms are qualitatively different. An extreme example of this is
systems in the Marcus inverted regime, for which the adiabatic rate is
undefined.[55] In light of this we think that the interpolation formula can be
complemented by a first-principles rate theory tackling the intermediate
nonadiabatic regime.

A general instanton solution for arbitrary diabatic coupling strengths has
been proposed by Voth and co-workers.[61, 62, 63, 64, 65] The solution assumes
that the reaction rate is related to the imaginary part of the barrier
partition function. This assumption, known as the
$\mathord{\operatorname{Im}}\mkern 1.0mu\mkern 2.0muF$ premise,[66, 67, 68,
69] has been validated for both adiabatic[70, 6, 8] and golden-rule
instantons[63, 20] by demonstrating that the
$\mathord{\operatorname{Im}}\mkern 1.0mu\mkern 2.0muF$ expressions recover the
semiclassical limits of the corresponding exact rate theories based on the
flux-correlation formalism. The same has not yet been accomplished in the
general nonadiabatic case, and the definitions of (imaginary) barrier
partition functions used by Voth and co-workers are not rigorously justified.
Depending on which definition is used, this can result, for example, in the
wrong golden-rule limit, forcing one to resort to “ _ad hoc_ ”
corrections.[65] Another limitation of this approach is that, like the
interpolation formula, it cannot at present be applied in the Marcus inverted
regime. Hence there remains room for theoretical development.

To this end, we derive a rigorous semiclassical instanton theory of the second
term in the perturbation expansion of the full rate. The term follows
immediately after the golden rule and is fourth-order in the diabatic coupling
$\Delta$. In Sec. II we define the rate problem, develop a perturbation series
for the rate in the diabatic representation and give a brief summary of how
the semiclassical limit of the leading-order term (the GR instanton) can be
derived. In Sec. III we derive the quantum mechanical fourth-order term,
casting it as a sum of three components that are individually amenable to
approximation by steepest descent. Applying the approximation in Sec. IV we
arrive at a semiclassical formula for the new term, whose numerical accuracy
is tested on a one-dimensional predissociation model and a multidimensional
spin–boson system in Sec. V. Our results and outlook on future work are
summarised in Sec. VI.

Because the fourth-order rate comprises three components, and the asymptotic
behaviour of each has to be analysed individually, the amount of mathematical
manipulation in Secs. III and IV is considerable. The salient point is that
the semiclassical instanton expressions in this work derive from a rigorous
flux-correlation formulation of the quantum rate. For the final result one can
proceed directly to Eqs. (4.22)–(4.28), Eqs. (4.39)–(4.45) and Eqs.
(4.47)–(4.50), which define the three components of the fourth-order rate
constant. Despite the lengthy derivations, these expressions are all easily
evaluated and require little computational effort beyond a golden-rule
instanton calculation.

##  II Theoretical background

We consider a system that comprises two diabatic states,[46] $\ket{0}$ and
$\ket{1}$, with corresponding nuclear Hamiltonian operators

$\hat{H}_{n}=\sum_{j=1}^{f}\left[\frac{\hat{p}_{j}^{2}}{2m}+V_{n}(\hat{\mathsf{x}})\right]\\!,$
(2.1)

where $n\in\\{0,1\\}$. Here $f$ is the number of nuclear coordinates, which
have been mass-weighted so that each has the same mass $m$, and $V_{n}$ are
the diabatic potential energy surfaces (PESs). We assume that the system is in
the Marcus normal regime, implying that the diabatic PESs intersect along a
seam for which $\nabla V_{0}\cdot\nabla V_{1}<0$. The states interact via the
diabatic electronic coupling operator,
$\hat{\Delta}=\Delta(\hat{\mathsf{x}})$, which is assumed to vary slowly with
respect to $\mathsf{x}$. The total Hamiltonian expressed in the diabatic basis
is then

$\hat{H}=\hat{H}^{(0)}+\lambda\hat{H}^{(1)}=\left(\\!\begin{array}[]{cc}\hat{H}_{0}&0\\\
0&\hat{H}_{1}\end{array}\\!\right)+\lambda\left(\\!\begin{array}[]{cc}0&\,\hat{\Delta}\\\
\hat{\Delta}&\,0\end{array}\\!\right)\\!,$ (2.2)

where $\lambda$ is a dimensionless parameter used to track the order of the
perturbative terms. At the end of the derivation we set $\lambda=1$.
Identifying $\ket{0}$ and $\ket{1}$ as the reactant and product states
respectively, we introduce the operator $\hat{h}=\ket{1}\\!\bra{1}$ that
projects onto the product state. Its time derivative is the flux operator,

$\hat{F}=\frac{1}{\mathrm{i}\hbar}[\hat{h},\hat{H}]=\frac{1}{\mathrm{i}\hbar}\left(\begin{array}[]{cc}0&\\!-\hat{\Delta}\\\
\hat{\Delta}&\\!\hphantom{-}0\end{array}\right).$ (2.3)

In what follows we give the formally exact expression for the thermal rate
constant associated with the transition from $\ket{0}$ to $\ket{1}$ and
outline the derivation of its $\lambda\to 0$ limit (_i.e._ , Fermi’s golden
rule) and the corresponding semiclassical approximation.

###  II.1 The flux-correlation and Green’s function formalisms

Following Ref. 71, we define the reactant and product partition functions

$Z_{\mathrm{r}}=\operatorname{Tr}\\!\big{[}\mathrm{e}^{-\beta\hat{H}}(1-\hat{h})\big{]},\quad
Z_{\mathrm{p}}=\operatorname{Tr}\\!\big{[}\mathrm{e}^{-\beta\hat{H}}\hat{h}\big{]},$
(2.4)

and the initial reactant population

$\chi_{\mathrm{r}}=\frac{1}{Z_{\mathrm{r}}}\operatorname{Tr}\big{[}\hat{\rho}_{\mathrm{r}}(1-\hat{h})\big{]},$
(2.5)

with the normalisation
$Z_{\mathrm{r}}\equiv\operatorname{Tr}[\hat{\rho}_{\mathrm{r}}].$ Here
$\hat{\rho}_{\mathrm{r}}$ is the symmetric thermal density operator[72, 73]

$\hat{\rho}_{\mathrm{r}}=\mathrm{e}^{-\beta\hat{H}/2}(1-\hat{h})\,\mathrm{e}^{-\beta\hat{H}/2},$
(2.6)

which models the system at thermal equilibrium in the reactant well.[71]
Depending on the system, the initial reaction population may be exactly equal,
or very close to, one. In our theory the deviation of $\chi_{\mathrm{r}}$ from
unity cannot always be neglected, and so we use Eq. (2.21) of Ref. 71 to
define the rate $k(\beta)$ as

$\displaystyle k(\beta)$
$\displaystyle=\frac{\tilde{k}(\beta)}{\widetilde{Z}_{\mathrm{r}}}=\frac{1}{\widetilde{Z}_{\mathrm{r}}}\int_{0}^{\mathrlap{t_{\mathrm{pl}}}}\
\,c(t)\,\mathrm{d}t,$ (2.7a) $\displaystyle c(t)$
$\displaystyle=\frac{1}{\mathrm{i}\hbar}\operatorname{Tr}\left\\{[\hat{H},\hat{\rho}_{\mathrm{r}}]\,\mathrm{e}^{\mathrm{i}\hat{H}t/\hbar}\hat{F}\mathrm{e}^{-\mathrm{i}\hat{H}t/\hbar}\right\\}$
(2.7b) $\displaystyle\rule{0.0pt}{13.99995pt}\widetilde{Z}_{\mathrm{r}}$
$\displaystyle=Z_{\mathrm{r}}\left(\chi_{\mathrm{r}}-[1-\chi_{\mathrm{r}}]Z_{\mathrm{r}}/Z_{\mathrm{p}}\right),$
(2.7c)

where $c(t)$ is the flux correlation function and $\tilde{k}(\beta)$ is a
“reduced” thermal rate introduced for notational convenience. In cases when
$\chi_{\mathrm{r}}=1$ this rate expression reduces to the well-known form[73,
74] with $Z_{\mathrm{r}}$ in place of $\widetilde{Z}_{\mathrm{r}}$.

The expression assumes a separation of time scales, such that the flux
correlation function decays to zero for $t_{\mathrm{pl}}<t\ll
t_{\mathrm{rxn}}$. Here $t_{\mathrm{pl}}$ is some “plateau time” that is
significantly shorter than the reaction time scale $t_{\mathrm{rxn}}$.[17] In
this work we consider the limit of small diabatic coupling, for which
$t_{\mathrm{rxn}}\to\infty$, enabling us to also take the limit
$t_{\mathrm{pl}}\to\infty$.

With the terms in Eq. (2.2) treated as a reference Hamiltonian,
$\hat{H}^{(0)}$, and a perturbation, $\hat{H}^{(1)}$, we can expand the rate
constant $k(\beta)$ as a power series in $\lambda$, namely

$k(\beta)=\sum_{\nu=1}^{\infty}\lambda^{2\nu}k_{2\nu}(\beta).$ (2.8)

Note that in the above expression coefficients of odd powers of $\lambda$ are
identically zero. The leading (second-order) term is obtained by noting that
$\widetilde{Z}_{\mathrm{r}}\sim\widetilde{Z}_{\mathrm{r},0}\equiv\operatorname{Tr}\\!\big{[}\mathrm{e}^{-\beta\hat{H}_{0}}\big{]}$
and
$\mathrm{e}^{-\hat{H}z/\hbar}\sim\hat{\mathcal{K}}_{0}(z)\ket{0}\\!\bra{0}+\hat{\mathcal{K}}_{1}(z)\ket{1}\\!\bra{1}$
as $\lambda\to 0$, where we define the propagator
$\hat{\mathcal{K}}_{n}(z)=\mathrm{e}^{-\hat{H}_{n}z/\hbar}$ and complex time
$z\equiv\tau+\mathrm{i}t$. Here $\mathord{\operatorname{Re}}\mkern
1.0mu(z)=\tau$ corresponds to imaginary time, and
$\mathord{\operatorname{Im}}\mkern 1.0mu(z)=t$ corresponds to real time.
Substituting this into Eq. (2.7) with $t_{\mathrm{pl}}\to\infty$ yields upon
simplification[17, 60, 75]

$\displaystyle k_{2}(\beta)$
$\displaystyle=\frac{1}{\hbar^{2}\widetilde{Z}_{\mathrm{r},0}}\int_{-\infty}^{\infty}\\!c_{2}(z)\,\mathrm{d}t,$
(2.9a) $\displaystyle c_{2}(z)$
$\displaystyle=\operatorname{Tr}\left[\hat{\mathcal{K}}_{0}\big{(}\beta\hbar-z\big{)}\hat{\Delta}\hat{\mathcal{K}}_{1}\big{(}z\big{)}\hat{\Delta}\right].$
(2.9b)

The value of $k_{2}$ does not depend on $\tau$, since both
$\widetilde{Z}_{\mathrm{r},0}$ and the integral are $\tau$-independent. The
latter holds because the integrand $c_{2}(z)$ is an analytic function of $z$,
allowing the integration contour to be shifted along the imaginary-time axis
without affecting the final outcome. This enables one to choose $\tau$ in a
way that simplifies the evaluation of the integral, as discussed later.

Equation (2.9) is precisely the GR limit mentioned previously. To go beyond it
to the next-order contribution, $k_{4}(\beta)$, one could proceed via the
flux-correlation formalism, expanding the propagator in a time-dependent
perturbation series.[76, 77] Doing so results in a $k_{4}$ expressed in terms
of components that have similar magnitudes but opposite signs. Apart from the
numerical difficulties associated with evaluating such an expression
accurately, it is not obvious how it can be effectively approximated using
integration by steepest descent. We will show that a better starting point is
afforded by the Green’s function formalism.[78, 20, 46] Central to this is the
Green’s function operator, which is the Laplace transform of the full
propagator,

$\displaystyle\hat{\mathcal{G}}(E)$ $\displaystyle=\lim_{\eta\to
0^{+}}-\frac{\mathrm{i}}{\hbar}\int_{0}^{\infty}\mathrm{e}^{-\mathrm{i}\hat{H}t/\hbar}\,\mathrm{e}^{\mathrm{i}(E+\mathrm{i}\eta)t/\hbar}\,\mathrm{d}t$
(2.10) $\displaystyle=\lim_{\eta\to 0^{+}}\frac{1}{E+\mathrm{i}\eta-\hat{H}}.$

Its imaginary part can be more simply written as the Fourier transform of the
propagator, and is thus related to the density of states,

$\displaystyle\mathord{\operatorname{Im}}\mkern 1.0mu\hat{\mathcal{G}}(E)$
$\displaystyle=-\frac{1}{2\hbar}\int_{-\infty}^{\infty}\mathrm{e}^{-\mathrm{i}\hat{H}t/\hbar}\,\mathrm{e}^{\mathrm{i}Et/\hbar}\,\mathrm{d}t$
(2.11) $\displaystyle=-\pi\delta(E-\hat{H}).$

This can be used to calculate the cumulative reaction probability at energy
$E$,[73, 78]

$P(E)=2\hbar^{2}\operatorname{Tr}\left[\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}(E)\hat{F}\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}(E)\hat{F}\right],$ (2.12)

which in turn is related to the thermal rate constant in Eq. (2.7a) with
$t_{\mathrm{pl}}\to\infty$ via

$\tilde{k}(\beta)\equiv\int_{0}^{\infty}\\!\\!c(t)\,\mathrm{d}t=\frac{1}{2\pi\hbar}\int_{-\infty}^{\infty}\\!\mathrm{e}^{-\beta
E}P(E)\,\mathrm{d}E.$ (2.13)

One advantage of this formalism is that the perturbation series for the
Green’s function operator, readily obtained from the Dyson equation,[46] has a
particularly simple form,

$\hat{\mathcal{G}}=\sum_{\nu=0}^{\infty}\lambda^{\nu}\left[\hat{\mathcal{G}}^{(0)}\hat{H}^{(1)}\right]^{\\!\nu}\\!\hat{\mathcal{G}}^{(0)},$
(2.14)

where
$\hat{\mathcal{G}}^{(0)}=\hat{\mathcal{G}}_{0}\ket{0}\\!\bra{0}+\hat{\mathcal{G}}_{1}\ket{1}\\!\bra{1}$
is the Green’s function operator for the unperturbed problem, with
$\hat{\mathcal{G}}_{n}$ given by Eq. (2.10) with Hamiltonian $\hat{H}_{n}$.
Substituting this into Eqs. (2.12) and (2.13) leads to a perturbation series
for $\tilde{k}$ that is analogous to Eq. (2.8). More importantly, the
formalism leads to expressions for rate constants that are naturally separated
into terms amenable to semiclassical approximation. As mentioned previously, a
straightforward application of time-dependent perturbation theory in the flux-
correlation picture does not lead to the same helpful separation, even though
the overall expression is mathematically equivalent.

To illustrate how the two formalisms connect, let us re-derive Eq. (2.9)
starting from the perturbation series

$P(E)=\sum_{\nu=1}^{\infty}\lambda^{2\nu}P_{2\nu}(E)$ (2.15)

and considering its $\nu=1$ term,

$P_{2}(E)=4\operatorname{Tr}\left[\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{0}(E)\hat{\Delta}\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{1}(E)\hat{\Delta}\right].$ (2.16)

Taking the Boltzmann average and using the relation in Eq. (2.11) gives the
first term in the perturbation series for $\tilde{k}$,

$\displaystyle\tilde{k}_{2}(\beta)$
$\displaystyle=\frac{1}{2\pi\hbar}\int_{-\infty}^{\infty}\mathrm{e}^{-\beta
E}P_{2}(E)\,\mathrm{d}E$
$\displaystyle=\frac{1}{\hbar^{2}}\int_{-\infty}^{\infty}\frac{\mathrm{d}E}{2\pi\hbar}\int_{-\infty}^{\infty}\\!\\!\mathrm{d}u_{2}\int_{-\infty}^{\infty}\\!\\!\mathrm{d}u_{1}$
(2.17)
$\displaystyle\quad\quad\mathrm{e}^{-E(\beta\hbar-\mathrm{i}(u_{1}+u_{2}))/\hbar}\operatorname{Tr}\left[\hat{\mathcal{K}}_{0}(\mathrm{i}u_{2})\hat{\Delta}\hat{\mathcal{K}}_{1}(\mathrm{i}u_{1})\hat{\Delta}\right].$

The integral over energy can be simplified by the variable transformation
$\bar{t}=u_{1}+u_{2}$, $t=(u_{2}-u_{1})/2$, and the integration contour for
the new variable $\bar{t}$ can be shifted along the imaginary-time axis,
$\bar{t}\to\bar{t}-\mathrm{i}\beta\hbar$, so that the energy integral becomes

$\frac{1}{2\pi\hbar}\int_{-\infty}^{\infty}\mathrm{e}^{\mathrm{i}E\bar{t}/\hbar}\,\mathrm{d}E=\delta(\bar{t}).$
(2.18)

Integrating over $\bar{t}$ leaves

$\tilde{k}_{2}=\frac{1}{\hbar^{2}}\\!\int_{-\infty}^{\infty}\\!\\!\mathrm{d}t\,\operatorname{Tr}\\!\left[\hat{\mathcal{K}}_{0}\big{(}\tfrac{\beta\hbar}{2}-\mathrm{i}t\big{)}\hat{\Delta}\hat{\mathcal{K}}_{1}\big{(}\tfrac{\beta\hbar}{2}+\mathrm{i}t\big{)}\hat{\Delta}\right]\\!.$
(2.19)

At this point we recognize that the remaining integration variable can also be
shifted along the imaginary-time axis, and use
$k_{2}=\tilde{k}_{2}/\widetilde{Z}_{\mathrm{r},0}$ to recover Eq. (2.9). This
approach is generalised in Sec. III to derive the next term in the series,
$k_{4}$.

###  II.2 Semiclassical approximation

We now summarise previous derivations of the semiclassical approximation to
$k_{2}(\beta)$.[8, 24] Evaluating the quantum flux-correlation function is
only computationally feasible for low-dimensional or particularly simple model
systems. Applications to realistic chemical reactions call for approximations
to Eq. (2.9) that make the calculation computationally tractable. The approach
that we use is to expand the trace in terms of position eigenstates and
replace the exact quantum propagator with its semiclassical counterpart,[79,
80, 81]

$\braket{\mathsf{x}_{\mathrm{f}}}{\hat{\mathcal{K}}_{n}(\tau_{n})}{\mathsf{x}_{\mathrm{i}}}\sim\mathcal{K}^{\mathrm{sc}}_{n}(\mathsf{x}_{\mathrm{f}},\mathsf{x}_{\mathrm{i}},\tau_{n})=\sum_{\text{traj.}}\sqrt{\frac{C_{n}}{(2\pi\hbar)^{f}}}\,\mathrm{e}^{-S_{n}/\hbar},$
(2.20)

where ‘$\mathord{\sim}$’ denotes an asymptotic relationship[82] and the sum is
over all classical trajectories that travel from $\mathsf{x}_{\mathrm{i}}$ to
$\mathsf{x}_{\mathrm{f}}$ in imaginary time $\tau_{n}$. $S_{n}\equiv
S_{n}(\mathsf{x}_{\mathrm{f}},\mathsf{x}_{\mathrm{i}},\tau_{n})$ is the
corresponding stationary Euclidean action,

$S_{n}[\mathsf{x}(\mathord{\cdot})]=\int_{0}^{\tau_{n}}\left[\frac{m}{2}\left\lVert\frac{\mathrm{d}{\mathsf{x}}}{\mathrm{d}{u}}\right\rVert^{2}+V_{n}\big{(}\mathsf{x}(u)\big{)}\right]\mathrm{d}u,$
(2.21)

and

$C_{n}=\left\rvert-\frac{\partial^{2}{S_{n}}}{\partial{\mathsf{x}_{\mathrm{i}}}\partial{\mathsf{x}_{\mathrm{f}}}}\right\rvert.$
(2.22)

Hence we can write

$\displaystyle\tilde{k}_{2}\sim\iiint_{-\infty}^{\infty}\\!\\!\mathrm{d}\mathsf{x}^{\prime}\mathrm{d}\mathsf{x}^{\prime\prime}\mathrm{d}t\,\sqrt{\frac{C_{0}C_{1}}{(2\pi\hbar)^{2f}}}\frac{\Delta(\mathsf{x}^{\prime})\Delta(\mathsf{x}^{\prime\prime})}{\hbar^{2}}\,\mathrm{e}^{-S_{2}/\hbar},$
(2.23)

where
$S_{2}(\mathsf{x}^{\prime},\mathsf{x}^{\prime\prime},z)=S_{0}(\mathsf{x}^{\prime},\mathsf{x}^{\prime\prime},\beta\hbar-z)+S_{1}(\mathsf{x}^{\prime\prime},\mathsf{x}^{\prime},z)$.
The semiclassical approximation to the exact propagator in Eq. (2.20) can be
obtained by steepest descent integration of its path-integral
representation.[83] The same technique can be applied to approximate the
remaining integrals in Eq. (2.23). For a one-dimensional example this
typically reads

$\displaystyle\int_{a}^{b}g(x)\mkern
1.0mu\mathrm{e}^{-\phi(x)/\hbar}\,\mathrm{d}x$ (2.24)
$\displaystyle\qquad{}\sim g(\bar{x})\mkern
1.0mu\mathrm{e}^{-\phi(\bar{x})/\hbar}\\!\int_{-\infty}^{\infty}\\!\exp\\!\left[-\frac{\nabla^{2}\phi}{2\hbar}(x-\bar{x})^{2}\right]\mathrm{d}x$
$\displaystyle\qquad{}=\sqrt{\frac{2\pi\hbar}{\nabla^{2}\phi}}\,g(\bar{x})\,\mathrm{e}^{-\phi(\bar{x})/\hbar}$

as $\hbar\to 0$, where $\bar{x}$ is the point where $\phi(x)$ reaches its
minimum value on the interval $(a,b)$, such that the first and second
derivatives of $\phi$ satisfy $\nabla\phi(\bar{x})=0$ and
$\nabla^{2}\phi(\bar{x})>0$. In effect, the steepest-descent approximation
consists in replacing the exponentiated function with its Taylor series
expansion about $\bar{x}$, truncated at the second-order term. All remaining
factors are replaced with the leading-order terms in their series expansions
about the same point, and the integration bounds are extended to $\pm\infty$.
The resulting Gaussian integral can then be evaluated analytically. The
relative error associated with this approximation becomes vanishingly small as
$\hbar\to 0$, provided certain conditions are satisfied. In particular, we can
ignore higher-order terms in the Taylor series expansion of $\phi(x)$ as long
as $\nabla^{2}\phi$ is not itself vanishingly small. Similarly, we can extend
the integration range to span the entire real line provided the stationary
point $\bar{x}$ does not lie infinitesimally close to either of the bounds.
Finally, we can replace $g(x)$ with its value at the stationary point,
provided $g(x)$ does not vary rapidly in its vicinity. Later in the paper we
encounter cases where some of these conditions are not satisfied, at which
point the procedure is modified accordingly.

No such modification is needed in the case of Eq. (2.23), other than a
generalisation to multiple dimensions. The stationary point
$(\mathsf{x}^{\prime},\mathsf{x}^{\prime}\mathrlap{{}^{\prime}},\tau)$
satisfies[20, 8, 24]

$\displaystyle\frac{\partial{S_{2}}}{\partial{\mathsf{x}^{\prime}}}$
$\displaystyle=\mathsf{p}_{0}^{\prime}-\mathsf{p}_{1}^{\prime}=\mathsf{0}$
(2.25a)
$\displaystyle\frac{\partial{S_{2}}}{\partial{\mathsf{x}^{\prime\prime}}}$
$\displaystyle=\mathsf{p}_{1}^{\prime\prime}-\mathsf{p}_{0}^{\prime\prime}=\mathsf{0}$
(2.25b) $\displaystyle\frac{\partial{S_{2}}}{\partial{\tau}}$
$\displaystyle=E_{1}-E_{0}=0$ (2.25c)

where $E_{n}$ are the energies of a classical trajectory comprised of two
parts ($n\in\\{0,1\\}$). The first part corresponds to travelling on $V_{0}$
from $\mathsf{x}^{\prime\prime}$ to $\mathsf{x}^{\prime}$ over an imaginary
time $\beta\hbar-\tau$, with initial and final momenta
$\mathsf{p}_{0}^{\prime\prime}$ and $\mathsf{p}^{\prime}_{0}$. The second part
corresponds to travelling on $V_{1}$ from $\mathsf{x}^{\prime}$ to
$\mathsf{x}^{\prime\prime}$ over an imaginary time $\tau$. The conditions in
Eq. (2.25) follow directly from the definition of $S_{2}$ and the relations

$\frac{\partial{S_{n}}}{\partial{\mathsf{x}_{\mathrm{i}}}}=-\mathsf{p}_{\mathrm{i}},\quad\frac{\partial{S_{n}}}{\partial{\mathsf{x}_{\mathrm{f}}}}=\mathsf{p}_{\mathrm{f}},\quad\frac{\partial{S_{n}}}{\partial{\tau_{n}}}=E_{n}.$
(2.26)

The path that makes the combined $S_{2}$ action stationary is therefore a
periodic classical trajectory in imaginary time with period $\beta\hbar$,
comprised of a reactant and a product segment. Conservation of energy and
momentum imposed by Eq. (2.25) implies that at the hopping points
$\mathsf{x}^{\prime},\,\mathsf{x}^{\prime\prime}$, where the two segments
join, the trajectory is continuous. Together with time-reversal symmetry this
typically imposes
$\mathsf{x}^{\prime}=\mathsf{x}^{\prime\prime}=\mathsf{x}^{\ddagger}$, where
$\mathsf{x}^{\ddagger}$ is a point on the seam along which the reactant and
product potential energies are equal. This stationary path is commonly known
as the golden-rule instanton.

In the Marcus normal regime, which we assume throughout this work, the value
of $\tau$ that extremises the action is in the range
$0<\tau_{\mathrm{g}}<\beta\hbar$. The instanton is a first-order saddle point
of the action and can be found by discretising the trajectory (_i.e._ ,
representing it as a ring polymer) and performing multidimensional
optimisation of the resulting extended classical system.[20, 21, 8, 24] This
is closely related to the analogous procedure for adiabatic instantons,[3, 4,
5] except now one must optimise the action in $\tau$ as well as
$\mathsf{x}^{\prime}$ and $\mathsf{x}^{\prime\prime}$. At the end of the
optimisation one calculates the action $S_{2}$ and its Hessian

$\mathsf{\Sigma}_{2}=\left[\begin{array}[]{*{3}{>{\rule[-0.5em]{0pt}{1.75em}}C{2.5em}}}0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\mathsf{x}^{\prime}}\partial{\mathsf{x}^{\prime}}}$\@add@centering&0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\mathsf{x}^{\prime}}\partial{\mathsf{x}^{\prime\prime}}}$\@add@centering&0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\mathsf{x}^{\prime}}\partial{\tau}}$\@add@centering\cr
0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\mathsf{x}^{\prime\prime}}\partial{\mathsf{x}^{\prime}}}$\@add@centering&0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\mathsf{x}^{\prime\prime}}\partial{\mathsf{x}^{\prime\prime}}}$\@add@centering&0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\mathsf{x}^{\prime\prime}}\partial{\tau}}$\@add@centering\cr
0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\tau}\partial{\mathsf{x}^{\prime}}}$\@add@centering&0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\tau}\partial{\mathsf{x}^{\prime\prime}}}$\@add@centering&0pt]{0.0pt}{17.50002pt}\frac{\partial^{2}{S_{2}}}{\partial{\tau}^{2}}$\@add@centering\end{array}\right]$
(2.27)

at the stationary point
$(\mathsf{x}^{\ddagger},\mathsf{x}^{\ddagger},\tau_{\mathrm{g}})$. The
multidimensional generalisation of Eq. (2.24), along with the Cauchy–Riemann
equations[84] relating the partial derivatives with respect to $\tau$ and $t$,
gives

$\displaystyle\tilde{k}_{2}\sim\tilde{k}_{2}^{\mathrm{sc}}=\sqrt{2\pi\hbar}\mkern
1.0mu\frac{\Delta^{2}}{\hbar^{2}}\sqrt{\frac{C_{0}C_{1}}{-\lvert\mathsf{\Sigma}_{2}\rvert}}\mkern
1.0mu\mathrm{e}^{-S_{2}/\hbar}.$ (2.28)

To complete the derivation, the semiclassical approximation to the reactant
partition function $Z_{\mathrm{r},0}$ is evaluated following the same
approach. This time the trajectory that makes the action stationary is
collapsed at the bottom of the reactant well, and in the absence of
translational or rotational degrees of freedom the expression reduces to[8,
83]

$Z_{\mathrm{r},0}\sim
Z_{\mathrm{r},0}^{\mathrm{sc}}=\mathrm{e}^{-\beta\epsilon_{\mathrm{r}}}\prod_{j=1}^{f}\frac{1}{2}\operatorname{csch}\\!\left(\frac{\beta\hbar\mkern
1.0mu\omega_{\mathrm{r},j}}{2}\right),$ (2.29)

where $\epsilon_{\mathrm{r}}$ is the energy at the minimum of the reactant
well, and $\omega_{\mathrm{r},j}$ is the frequency of its $j$-th vibrational
normal mode. If present, translational and rotational modes can also be
accounted for.[24] Combining Eqs. (2.28) and (2.29) gives the final result,
$k_{2}\sim
k_{2}^{\mathrm{sc}}=\tilde{k}_{2}^{\mathrm{sc}}/Z_{\mathrm{r},0}^{\mathrm{sc}}$.

By representing the instanton trajectory as a ring polymer,[21, 8, 24] the
formula can be readily applied to realistic molecular potentials,[28] for
which the analytic form of the stationary action is not known. In this
representation finding the stationary action is equivalent to a
multidimensional optimisation problem, which can be solved efficiently with
well-established numerical techniques. In what follows we show how the output
of such a GR instanton optimisation can be used, with little additional
computational expense, to calculate the next semiclassical term in the
perturbative expansion of the thermal rate constant.

##  III Exact fourth-order rate expression

Substituting Eq. (2.14) into Eq. (2.12) shows that the fourth-order
contribution to the cumulative reaction probability is

$P_{4}(E)=-P_{{4\mathrm{A}}}(E)+P_{{4\mathrm{B}}}(E),$ (3.1)

where

$\displaystyle P_{{4\mathrm{A}}}(E)$
$\displaystyle=8\operatorname{Tr}\left[\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{0}\hat{\Delta}\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{1}\hat{\Delta}\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{0}\hat{\Delta}\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{1}\hat{\Delta}\right],$ (3.2a) $\displaystyle
P_{{4\mathrm{B}}}(E)$
$\displaystyle=8\operatorname{Tr}\left[\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{0}\hat{\Delta}\mathord{\operatorname{Re}}\mkern
1.0mu\hat{\mathcal{G}}_{1}\hat{\Delta}\mathord{\operatorname{Re}}\mkern
1.0mu\hat{\mathcal{G}}_{0}\hat{\Delta}\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{1}\hat{\Delta}\right].$ (3.2b)

Both terms are comprised of four $\hat{\mathcal{G}}_{n}$ factors, alternating
between the reactant ($n=0$) and product ($n=1$) diabats. This corresponds to
a total of four state changes in a single trace, as opposed to the two changes
in the golden-rule expression [Eq. (2.16)]. One could directly approximate the
terms in Eq. (3.2) with semiclassical methods, arriving at expressions that
have simple physical interpretations in terms of instanton trajectories. As
discussed in Ref. 20, in the forbidden regime the dominant contribution to
$\braket{\mathsf{x}_{\mathrm{f}}}{\mathord{\operatorname{Re}}\mkern
1.0mu\hat{\mathcal{G}}_{n}}{\mathsf{x}_{\mathrm{i}}}$ comes from a trajectory
that goes directly from $\mathsf{x}_{\mathrm{i}}$ to
$\mathsf{x}_{\mathrm{f}}$. On the other hand, the dominant contribution to
$\braket{\mathsf{x}_{\mathrm{f}}}{\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{n}}{\mathsf{x}_{\mathrm{i}}}$ comes from a trajectory
that reaches a turning point where $V_{n}(\mathsf{x})=E$. This feature is
known as a “bounce”. The GR probability in Eq. (2.16) corresponds therefore to
an instanton trajectory that bounces a total of two times: once off the
reactant and once off the product diabat, as shown schematically in Fig. 1(c).
The $P_{{4\mathrm{A}}}$ term corresponds to an instanton with four bounces, as
in Fig. 1(b), and the $P_{{4\mathrm{B}}}$ term to an instanton with two
bounces and two consecutive direct segments on alternating diabats.

Despite the ease of physical interpretation, the semiclassical Green’s
function is not as well behaved as the imaginary-time propagator.[20, 85] For
this reason we convert the microcanonical reaction probabilities into thermal
rate constants, expressed in terms of imaginary-time propagators. Throughout
the derivation it is assumed that we can take the limit
$t_{\mathrm{pl}}\to\infty$, just as in the case of the GR rate constant. The
derivation then amounts to finding the $\mathcal{O}(\lambda^{4})$ term
$\tilde{k}_{4}=\tilde{k}_{{4\mathrm{B}}}-\tilde{k}_{{4\mathrm{A}}}$ in the
expansion of $\tilde{k}(\beta)$ [Eq. (2.13)] and the
$\mathcal{O}(\lambda^{2})$ term $\widetilde{Z}_{\mathrm{r},2}$ of
$\widetilde{Z}_{\mathrm{r}}$ [Eq. (2.7c)]. Together these give the total
fourth-order rate constant

$k_{4}=k_{2}\left(-\frac{\tilde{k}_{{4\mathrm{A}}}}{\tilde{k}_{2}}+\frac{\tilde{k}_{{4\mathrm{B}}}}{\tilde{k}_{2}}-\frac{\widetilde{Z}_{\mathrm{r},2}}{Z_{\mathrm{r},0}}\right).$
(3.3)

Continuing in the same fashion one can obtain analogous expressions for
$k_{6},\,k_{8},$ etc., which we intend to pursue in future work.

###  III.1 _A_ -type term

Following a similar approach to Sec. II.1, we use the Fourier transform
representation of the imaginary Green’s function operator to yield

$\displaystyle\tilde{k}_{{4\mathrm{A}}}=\frac{1}{2\hbar^{4}}\int_{-\infty}^{\infty}\frac{\mathrm{d}E}{2\pi\hbar}\int\mathrm{d}^{4}\bm{u}\,\mathrm{e}^{-E\left(\beta\hbar-\mathrm{i}\sum
u_{i}\right)/\hbar}$ (3.4)
$\displaystyle{}\times\operatorname{Tr}\left[\hat{\mathcal{K}}_{0}(\mathrm{i}u_{4})\hat{\Delta}\hat{\mathcal{K}}_{1}(\mathrm{i}u_{3})\hat{\Delta}\hat{\mathcal{K}}_{0}(\mathrm{i}u_{2})\hat{\Delta}\hat{\mathcal{K}}_{1}(\mathrm{i}u_{1})\hat{\Delta}\right],$

where
$\int\mathrm{d}^{4}\bm{u}\equiv\int_{-\infty}^{\infty}\mathrm{d}u_{1}\cdots\int_{-\infty}^{\infty}\mathrm{d}u_{4}$.
Applying the variable transformation

$\displaystyle\begin{bmatrix}[1.5]u_{1}\\\ u_{2}\\\ u_{3}\\\
u_{4}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}[1.5]0&0&-1&\frac{1}{2}\\\
\frac{1}{2}&-1&0&-\frac{1}{2}\\\ 0&0&1&\frac{1}{2}\\\
\frac{1}{2}&1&0&-\frac{1}{2}\end{bmatrix}\begin{bmatrix}[1.5]\bar{t}\\\
t_{0}\\\ t_{1}\\\ t\end{bmatrix},$ (3.5)

allows us to integrate over $E$ and $\bar{t}$ in the same fashion as before.
The remaining integration variables are shifted into the complex plane, so
that $\mathrm{i}t_{0}\to z_{0}$, $\mathrm{i}t_{1}\to z_{1}$ and
$\mathrm{i}t\to z$, where $z_{n}\equiv\tau_{n}+\mathrm{i}t_{n}$. This leads to
the final expression

$\displaystyle\tilde{k}_{{4\mathrm{A}}}=\frac{1}{2\hbar^{4}}\int_{-\infty}^{\infty}\\!\\!\mathrm{d}t_{0}\int_{-\infty}^{\infty}\\!\\!\mathrm{d}t_{1}\int_{-\infty}^{\infty}\\!\\!\mathrm{d}t\
c_{{4\mathrm{A}}}(z_{0},z_{1},z),$ (3.6a) $\displaystyle
c_{{4\mathrm{A}}}(z_{0},z_{1},z)=\operatorname{Tr}\\!\bigg{[}\hat{\mathcal{K}}_{0}\Big{(}\tfrac{z_{\mathrm{r}}}{2}+z_{0}\Big{)}\hat{\Delta}\,\hat{\mathcal{K}}_{1}\Big{(}\tfrac{z_{\mathrm{p}}}{2}+z_{1}\Big{)}\hat{\Delta}$
$\displaystyle\qquad\qquad\quad{}\times\hat{\mathcal{K}}_{0}\Big{(}\tfrac{z_{\mathrm{r}}}{2}-z_{0}\Big{)}\hat{\Delta}\,\hat{\mathcal{K}}_{1}\Big{(}\tfrac{z_{\mathrm{p}}}{2}-z_{1}\Big{)}\hat{\Delta}\bigg{]},$
(3.6b)

where we introduce the notation $z_{\mathrm{r}}\equiv\beta\hbar-z$ and
$z_{\mathrm{p}}\equiv z$. The expression defines the exact quantum
$\tilde{k}_{{4\mathrm{A}}}$ and will be given a semiclassical treatment in
Sec. IV.1.

###  III.2 _B_ -type term

We rewrite the second component of the fourth-order reaction probability as

$P_{{4\mathrm{B}}}=P_{{4\mathrm{A}}}+8\operatorname{Re}\left\\{\operatorname{Tr}\left[\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{0}\mkern 1.0mu\hat{\Delta}\mkern
1.0mu\hat{\mathcal{G}}_{1}\hat{\Delta}\mkern 1.0mu\hat{\mathcal{G}}_{0}\mkern
1.0mu\hat{\Delta}\mkern 1.0mu\mathord{\operatorname{Im}}\mkern
1.0mu\hat{\mathcal{G}}_{1}\hat{\Delta}\right]\right\\}$ (3.7)

and recast the Green’s function operators as integral transforms according to
Eqs. (2.10) and (2.11), to yield

$\displaystyle\tilde{k}_{{4\mathrm{B}}}=\tilde{k}_{{4\mathrm{A}}}+\operatorname{Re}\bigg{\\{}\\!{-\frac{2}{\hbar^{4}}}\\!\int_{-\infty}^{\infty}\frac{\mathrm{d}E}{2\pi\hbar}\int\\!\mathrm{d}^{4}\bm{u}\
\mathrm{e}^{-E\left(\beta\hbar-\mathrm{i}\sum u_{i}\right)/\hbar}$
$\displaystyle\
{}\times\operatorname{Tr}\left[\hat{\mathcal{K}}_{0}(\mathrm{i}u_{4})\hat{\Delta}\hat{\mathcal{K}}_{1}(\mathrm{i}u_{3})\hat{\Delta}\hat{\mathcal{K}}_{0}(\mathrm{i}u_{2})\hat{\Delta}\hat{\mathcal{K}}_{1}(\mathrm{i}u_{1})\hat{\Delta}\right]\\!\bigg{\\}},$
(3.8)

where the integration ranges are now

$\int\mathrm{d}^{4}\bm{u}\equiv\int_{-\infty}^{\infty}\\!\\!\mathrm{d}u_{1}\int_{0}^{\infty}\\!\mathrm{d}u_{2}\int_{0}^{\infty}\\!\mathrm{d}u_{3}\int_{-\infty}^{\infty}\\!\\!\mathrm{d}u_{4}.$
(3.9)

Under the variable transformation

$\displaystyle\begin{bmatrix}[1.5]u_{1}\\\ u_{2}\\\ u_{3}\\\
u_{4}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}[1.5]0&0&-1&1\\\ 0&1&0&0\\\
0&0&1&0\\\ 1&-1&0&-1\end{bmatrix}\begin{bmatrix}[1.5]\bar{t}\\\ t_{0}\\\
t_{1}\\\ t\end{bmatrix}$ (3.10)

the expression simplifies and can be integrated over $E$ and $\bar{t}$ as
before, so that

$\displaystyle\tilde{k}_{{4\mathrm{B}}}=\tilde{k}_{{4\mathrm{A}}}+\frac{2}{\hbar^{4}}\operatorname{Re}\bigg{\\{}\int_{0}^{\mathrlap{\mathrm{i}\infty}}\
\mathrm{d}z_{0}\\!\int_{0}^{\mathrlap{\mathrm{i}\infty}}\
\mathrm{d}z_{1}\\!\int_{\mathrlap{-\infty}}^{\mathrlap{\infty}}\ \mathrm{d}t\
c_{{4\mathrm{B}}}\bigg{\\}}$ (3.11a) $\displaystyle
c_{{4\mathrm{B}}}(z_{0},z_{1},z)=\begin{aligned}
\operatorname{Tr}\Big{[}&\hat{\mathcal{K}}_{0}(z_{\mathrm{r}}-z_{0})\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{1}(z_{1})\hat{\Delta}\\\
{}\times{}&\hat{\mathcal{K}}_{0}(z_{0})\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{1}(z_{\mathrm{p}}-z_{1})\hat{\Delta}\Big{]}.\end{aligned}$
(3.11b)

The integration contours for $z_{n}$ can be deformed to simplify evaluation.
Introducing

$L_{{4\mathrm{B}}}(z_{0},z_{1})=\frac{2}{\hbar^{4}}\int_{-\infty}^{\infty}\\!c_{{4\mathrm{B}}}(z_{0},z_{1},z)\,\mathrm{d}t,$
(3.12)

we apply Cauchy’s integral theorem[84] iteratively to get

$\displaystyle\int_{0}^{\mathrlap{\mathrm{i}\infty}}\
\,\int_{0}^{\mathrlap{\mathrm{i}\infty}}\ \
L_{{4\mathrm{B}}}(z_{0},z_{1})\,\mathrm{d}z_{0}\mathrm{d}z_{1}={}$ (3.13)
$\displaystyle\left[\int\limits_{0}^{\mathrlap{\tau_{0}^{*}}}\int\limits_{0}^{\mathrlap{\tau_{1}^{*}}}\,+\,\int\limits_{0}^{\mathrlap{\tau_{0}^{*}}}\int\limits_{\tau_{1}^{*}}^{\mathrlap{\sigma_{1}^{*}}}\,+\,\int\limits_{\tau_{0}^{*}}^{\mathrlap{\sigma_{0}^{*}}}\int\limits_{0}^{\mathrlap{\tau_{1}^{*}}}\,+\,\int\limits_{\tau_{0}^{*}}^{\mathrlap{\sigma_{0}^{*}}}\int\limits_{\tau_{1}^{*}}^{\mathrlap{\sigma_{1}^{*}}}\
\ \right]\\!L_{{4\mathrm{B}}}(z_{0},z_{1})\,\mathrm{d}z_{0}\mathrm{d}z_{1},$

where $\sigma^{*}_{n}\equiv\tau^{*}_{n}+\mathrm{i}\infty$. For future
reference we label the integration domains on the right-hand side of Eq.
(3.13) as $\mathcal{A}_{\tau\tau}$, $\mathcal{A}_{\tau t}$,
$\mathcal{A}_{t\tau}$ and $\mathcal{A}_{tt}$. A suitable choice of
$\tau_{n}^{*}$ simplifies the numerical evaluation of the integral by
minimising the oscillations in $L_{{4\mathrm{B}}}(z_{0},z_{1})$. The final
result of this section gives the formally exact definition of the quantum
$\tilde{k}_{{4\mathrm{B}}}$, and will be given a semiclassical treatment in
Sec. IV.2.

At this point it is worth re-emphasizing that we have pursued this particular
route to $\tilde{k}_{4}$ because the $\tilde{k}_{{4\mathrm{A}}}$ and
$\tilde{k}_{{4\mathrm{B}}}$ terms can be well approximated by semiclassical
techniques. As expected from time-dependent perturbation theory,[77] both
terms are expressed as integrals over three-time correlation functions. The
crucial difference is that the application of time-dependent perturbation
theory to Eq. (2.13) splits $\tilde{k}_{4}$ into a different set of
components, not readily amenable to steepest-descent integration.

###  III.3 Partition function

To derive the $\mathcal{O}(\lambda^{2})$ term in the perturbation expansion of
$\widetilde{Z}_{\mathrm{r}}$, we use[86]

$\displaystyle Z_{\mathrm{r}}$
$\displaystyle=\sum_{\nu=0}^{\infty}\lambda^{2\nu}Z_{\mathrm{r},2\nu},$
(3.14a) $\displaystyle Z_{\mathrm{r},2\nu}$
$\displaystyle=\frac{1}{\hbar^{2\nu}}\int_{0}^{\beta\hbar}\\!\\!\mathrm{d}u_{2\nu}\int_{0}^{u_{2\nu}}\\!\\!\mathrm{d}u_{2\nu-1}\cdots\int_{0}^{u_{2}}\\!\\!\mathrm{d}u_{1}$
(3.14b)
$\displaystyle\quad\operatorname{Tr}\Bigg{[}\mathrm{e}^{-\beta\hat{H}_{0}}\\!\prod_{\mu=0}^{\nu-1}\hat{\mathcal{K}}_{0}(r_{\nu-\mu})\hat{\Delta}\mkern
1.0mu\hat{\mathcal{K}}_{1}(s_{\nu-\mu})\hat{\Delta}\Bigg{]},$

where $s_{\mu}\equiv u_{2\mu}-u_{2\mu-1}$, $r_{\mu}\equiv u_{2\mu+1}-u_{2\mu}$
and $u_{2\nu+1}\equiv u_{1}$. Explicitly, the first two terms are

$\displaystyle Z_{\mathrm{r},0}$
$\displaystyle=\operatorname{Tr}\left[\mathrm{e}^{-\beta\hat{H}_{0}}\right],$
(3.15a) $\displaystyle Z_{\mathrm{r},2}$
$\displaystyle=\frac{1}{\hbar^{2}}\int_{0}^{\beta\hbar}\\!\\![\beta\hbar-u]\,c_{2}(u)\,\mathrm{d}u.$
(3.15b)

The integrand in Eq. (3.15b) contains a factor that is precisely the golden-
rule flux correlation function in Eq. (2.9), except it is now integrated over
imaginary time instead of $t$. Analogous results for $Z_{\mathrm{p},2\nu}$ can
be obtained by exchanging the diabatic state labels $0\leftrightarrow 1$,
_e.g._ ,
$Z_{\mathrm{p},0}=\operatorname{Tr}\big{[}\mathrm{e}^{-\beta\hat{H}_{1}}\big{]}$.
We can write a similar perturbation series for the initial reactant
population,

$\chi_{\mathrm{r}}=\sum_{\nu=0}^{\infty}\lambda^{2\nu}\chi_{\mathrm{r},2\nu}$
(3.16)

where $\chi_{\mathrm{r},0}=1$ and

$\displaystyle\chi_{\mathrm{r},2}=\frac{1}{Z_{\mathrm{r},0}}\bigg{\\{}$
$\displaystyle\frac{1}{\hbar^{2}}\int_{0}^{\frac{\beta\hbar}{2}}\\![\beta\hbar-2u]\,c_{2}(u)\,\mathrm{d}u-Z_{\mathrm{r},2}\bigg{\\}}.$
(3.17)

It follows that

$\widetilde{Z}_{\mathrm{r},2}=Z_{\mathrm{r},2}+Z_{\mathrm{r},0}\,\chi_{\mathrm{r},2}\left(1+\frac{Z_{\mathrm{r},0}}{Z_{\mathrm{p},0}}\right),$
(3.18)

which, together with our previous expressions for $\tilde{k}_{{4\mathrm{A}}}$
and $\tilde{k}_{{4\mathrm{B}}}$, gives the overall $k_{4}$ rate. This can be
combined with $k_{2}$ (using straightforward summation or otherwise—see Sec.
V.2) to define the total rate up to fourth order in the diabatic coupling. It
can be shown that under this definition, the total rate satisfies detailed
balance for the forward and backward reactions to the specified order in the
perturbation.

##  IV Fourth-order instanton rate theory

Having split $k_{4}$ into three different terms, we now derive a semiclassical
approximation for each of them in turn. Derivations for
$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$ and
$\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}$ follow the same pattern. First, we
consider the high-temperature limit, which happens to yield the exact
expression for a system of two one-dimensional linear diabatic potentials,
$V_{n}(x)=V^{\ddagger}+\kappa_{n}x$. Then we generalise this result to
arbitrary temperatures and potentials. The remaining term, which describes
changes in the reactant partition function, is obtained directly, without
special consideration of the high-temperature regime. To keep the derivations
as simple as possible we assume here that there are no translational or
rotational normal modes. If such zero-frequency modes are present, they can be
accounted for as described in Refs. 8, 24.

###  IV.1 Semiclassical _A_ -type contribution

In order to derive the semiclassical instanton formulation of
$\tilde{k}_{{4\mathrm{A}}}$, we insert a set of four position-state
resolutions of the identity into Eq. (3.6) and replace each instance of
$\braket{\mathsf{x}_{i}}{\hat{\mathcal{K}}_{n}}{\mathsf{x}_{j}}$ with the
corresponding semiclassical propagator [Eq. (2.20)]. This gives rise to the
combined action

$\displaystyle S_{{4\mathrm{A}}}$ $\displaystyle\equiv
S_{0}\big{(}\mathsf{x}_{1},\mathsf{x}_{4},\tfrac{z_{\mathrm{r}}}{2}+z_{0}\big{)}+S_{1}\big{(}\mathsf{x}_{4},\mathsf{x}_{3},\tfrac{z_{\mathrm{p}}}{2}+z_{1}\big{)}$
$\displaystyle{}+S_{0}\big{(}\mathsf{x}_{3},\mathsf{x}_{2},\tfrac{z_{\mathrm{r}}}{2}-z_{0}\big{)}+S_{1}\big{(}\mathsf{x}_{2},\mathsf{x}_{1},\tfrac{z_{\mathrm{p}}}{2}-z_{1}\big{)}.$
(4.1)

We then look for the path that makes the action stationary. As discussed in
Sec. III, we expect this to be a four-bounce trajectory comprised of two
identical loops, as depicted in Fig. 1(b).

Figure 1:  Instantons for a nonadiabatic reaction involving a pair of diabats
$V_{n}(\mathsf{x})$. The blue and orange segments are classical imaginary-time
trajectories on $V_{0}$ and $V_{1}$ that traverse a path starting and ending
at the hopping point $\mathsf{x}^{\ddagger}$, bouncing once at the turning
point. (a) shows a 4A (top) and a GR (bottom) instanton; (b) and (c) show
expanded diagrams of these two instantons, in the same order. Here the
vertical dimension has no physical meaning and is introduced for clarity. Each
of the two loops of the 4A instanton in (b) is identical to a GR instanton at
twice the temperature.

#### High-temperature limit.

As $\beta\to 0$, the action in Eq. (IV.1) becomes stationary for trajectories
in the immediate vicinity of the minimum-energy crossing point (MECP) between
the diabats. For the sake of convenience, we position it at the origin of our
coordinate system. The compactness of the high-temperature instanton allows us
to approximate the diabatic PESs as[20]

$V_{n}(\mathsf{x})\sim V^{\ddagger}+\kappa_{n}q+\frac{1}{2}\mathsf{Q}^{\mkern
0.0mu\top}\mathsf{H}_{n}\mathsf{Q}.$ (4.2)

The gradients $\mathsf{g}_{n}=\nabla V_{n}(\mathsf{0})$ of the two diabatic
surfaces are antiparallel at the MECP, and we use $\kappa_{n}$ to denote their
signed norms,
$\kappa_{0}=\lVert\mathsf{g}_{0}\rVert,\,\kappa_{1}=-\lVert\mathsf{g}_{1}\rVert$.
We call the Cartesian coordinate aligned with the gradients $q$, and the $f-1$
orthogonal coordinates $\mathsf{Q}$. Each appears in the expansion to its
lowest non-zero order. The results in this section can equally be derived
using a more general expansion, keeping _all_ second-order terms, but such a
derivation is more laborious and has the same $\beta\to 0$ limit. In this
limit, the terms appearing in the semiclassical propagator are

$\displaystyle
C_{n}(\mathsf{x}^{\prime\prime},\mathsf{x}^{\prime},z)\sim\left(\frac{m}{2\pi\hbar
z}\right)^{\\!f/2},$ (4.3a) $\displaystyle
S_{n}(\mathsf{x}^{\prime\prime},\mathsf{x}^{\prime},z)\sim
S_{n}^{\mathrm{lin}}(q^{\prime\prime},q^{\prime},z)+S_{n}^{\mathrm{vib}}(\mathsf{Q}^{\prime\prime},\mathsf{Q}^{\prime},z),$
(4.3b) $\displaystyle
S_{n}^{\mathrm{lin}}(q^{\prime\prime},q^{\prime},z)=V^{\ddagger}z+\frac{mq_{-}^{2}}{2z}+z\kappa_{n}q_{+}-\frac{\kappa_{n}^{2}z^{3}}{24m},$
(4.3c) $\displaystyle
S_{n}^{\mathrm{vib}}(\mathsf{Q}^{\prime\prime},\mathsf{Q}^{\prime},z)\sim\frac{m\lVert\mathsf{Q}_{-}\rVert^{2}}{2z}+\frac{z}{2}\mathsf{Q}_{+}^{\mkern
0.0mu\top}\mathsf{H}_{n}\mathsf{Q}_{+}.$ (4.3d)

Here $q_{-}=q^{\prime\prime}-q^{\prime}$ and
$q_{+}=(q^{\prime}+q^{\prime\prime})/2$, with analogous definitions for
$\mathsf{Q}_{\pm}$. The stationary trajectory follows a path along $q$, with
no lateral displacement ($\mathsf{Q}=\mathsf{0}$). Since the combined action
is quadratic in the positions, integration over these variables can be
performed analytically to yield $c_{{4\mathrm{A}}}(z_{0},z_{1},z)$. Choosing
the real parts of its arguments such that
$\tau=\tfrac{\beta\hbar\kappa_{0}}{\kappa_{0}-\kappa_{1}}$ and
$\tau_{0}=\tau_{1}=0$ makes the three-time correlation function proportional
to $\delta(t)$. By analogy with $L_{{4\mathrm{B}}}$ in Eq. (3.12), we then
define

$\displaystyle
L_{{4\mathrm{A}}}(t_{0},t_{1})=\frac{1}{2\hbar^{4}}\int_{-\infty}^{\infty}c_{{4\mathrm{A}}}(\mathrm{i}t_{0},\mathrm{i}t_{1},\tau+\mathrm{i}t)\,\mathrm{d}t$
(4.4a) $\displaystyle{}\sim\frac{\Delta^{4}}{\hbar^{4}}\sqrt{\frac{\pi
m}{2\beta(\kappa_{0}-\kappa_{1})^{2}}}\,Z^{\ddagger}\mkern
1.0mu\mathrm{e}^{-\beta
V^{\ddagger}}\exp\\!\left[\frac{\beta^{3}\hbar^{2}\kappa_{0}^{2}\kappa_{1}^{2}}{96m(\kappa_{0}-\kappa_{1})^{2}}\right]$
$\displaystyle\qquad\qquad\qquad{}\times\exp\\!\left[-\frac{d_{0}t_{0}^{2}+d_{1}t_{1}^{2}+d_{01}t_{0}^{2}t_{1}^{2}}{2\hbar}\right]$
(4.4b)

as $\beta\to 0$. Here $\Delta$ is evaluated at the MECP, and we define a
transition-state partition function

$\displaystyle
Z^{\ddagger}=\bigg{\lvert}\frac{\beta^{2}\hbar^{2}\widetilde{\mathsf{H}}}{m}\bigg{\rvert}^{-1/2}\
\
\text{with}\qquad\widetilde{\mathsf{H}}=\frac{\kappa_{0}\mathsf{H}_{1}-\kappa_{1}\mathsf{H}_{0}}{\kappa_{0}-\kappa_{1}},$
(4.5)

as well as the coefficients

$d_{n}=\beta\hbar\,\frac{\kappa_{n}^{2}}{4m},\qquad
d_{01}=\frac{1}{\beta\hbar}\frac{(\kappa_{0}-\kappa_{1})^{2}}{m}.$ (4.6)

In general, the asymptotic relation in Eq. (4.4b) only holds at high
temperatures and short times $t_{0}\text{ and }t_{1}$. However the expression
is exact for one-dimensional linear diabats
$V_{n}(x)=V^{\ddagger}+\kappa_{n}x$. It is therefore meaningful to make a
digression into the low-temperature regime ($\beta\to\infty$). In this limit,
the coefficient $d_{01}$ of the quartic term in Eq. (4.4b) becomes vanishingly
small, and so the marginalised time-correlation function
$L_{{4\mathrm{A}}}(t_{0},t_{1})$ becomes well approximated by a Gaussian, as
seen in Fig. 2(a). Neglecting the quartic term and integrating over $t_{n}$
yields the low-temperature limit of $\tilde{k}_{{4\mathrm{A}}}$ for a linear
system,

$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{lt,lin}}=\left[\frac{2\pi
m}{\beta(\kappa_{0}-\kappa_{1})^{2}}\right]^{1/2}\frac{4\pi m\Delta^{4}\mkern
2.0mu\mathrm{e}^{-\beta
V^{\ddagger}}}{\beta\hbar^{4}\lvert\kappa_{0}\kappa_{1}\rvert}\exp\\!\left(\frac{\widetilde{\beta}^{3}}{48}\right)\\!,$
(4.7)

to be compared with the golden-rule rate[20]

$\displaystyle\tilde{k}_{2}^{\mathrm{lin}}=\left[\frac{2\pi
m}{\beta(\kappa_{0}-\kappa_{1})^{2}}\right]^{1/2}\frac{\Delta^{2}\mkern
2.0mu\mathrm{e}^{-\beta
V^{\ddagger}}}{\hbar^{2}}\,\exp\\!\left(\frac{\widetilde{\beta}^{3}}{12}\right),$
(4.8)

where

$\widetilde{\beta}=2\left[\frac{d_{0}d_{1}}{\hbar
d_{01}}\right]^{\mathrlap{\\!1/3}}\
=\beta\left[\frac{\hbar^{2}\kappa_{0}^{2}\kappa_{1}^{2}}{2m(\kappa_{0}-\kappa_{1})^{2}}\right]^{\mathrlap{\\!1/3}}.$
(4.9)

The similarity between Eqs. (4.7) and (4.8) arises because at low temperatures
steepest-descent integration of $c_{{4\mathrm{A}}}$ is a straightforward
generalisation of the procedure for the golden-rule flux correlation function
$c_{2}$. One can therefore expect to see some similarities between the
physical interpretations of these two expressions.

In the golden-rule case, the value of the rate constant $\tilde{k}_{2}$ is
determined by the behaviour of $c_{2}(\tau+\mathrm{i}t)$ in the small-$t$
limit. Resolved in terms of position eigenstates, the function at $t=0$ reads

$c_{2}(\tau)=\iint_{-\infty}^{\infty}\\!\\!\mathrm{d}\mathsf{x}_{\mathrm{i}}\mathrm{d}\mathsf{x}_{\mathrm{f}}\left\lVert\braket{\mathsf{x}_{\mathrm{f}}}{\hat{\mathcal{K}}_{1}\\!\left(\tfrac{\tau_{\mathrm{p}}}{2}\right)\\!\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{0}\\!\left(\tfrac{\tau_{\mathrm{r}}}{2}\right)\\!}{\mathsf{x}_{\mathrm{i}}}\right\rVert^{2}\\!\\!,$
(4.10)

which is dominated by $\mathsf{x}_{\mathrm{i}}$ and $\mathsf{x}_{\mathrm{f}}$
near the GR instanton turning points [see Fig. 1(c)]. Reading the integrand
from right to left, the dominant contribution corresponds to the probability
of the system starting near the reactant turning point, evolving in imaginary
time for $\tau_{\mathrm{r}}/2\equiv(\beta\hbar-\tau)/2$, “switching” from
$V_{0}$ to $V_{1}$, propagating for $\tau_{\mathrm{p}}/2\equiv\tau/2$, and
finally arriving at the product turning point. Imaginary-time evolution
indicates that the transition proceeds via a tunnelling mechanism, with the
dominant tunnelling pathway following the half-instanton from
$\mathsf{x}_{\mathrm{i}}$ to $\mathsf{x}_{\mathrm{f}}$. The dynamics before
and after tunnelling is not explicitly included in the half-instanton but may
be deduced from the locations of the turning points and the corresponding
momenta (which are zero). Given these boundary conditions we conclude that the
thermally activated system initially moves up $V_{0}$ with just enough energy
to reach $\mathsf{x}_{\mathrm{i}}$, at which point it tunnels to
$\mathsf{x}_{\mathrm{f}}$ and continues moving down $V_{1}$. This is precisely
the kind of qualitative picture we were looking for in order to describe the
processes underlying $\tilde{k}_{2}$, and it was possible to deduce from just
the GR half-instanton.

Figure 2:  Plots of $L_{{4\mathrm{A}}}(t_{0},t_{1})/L_{{4\mathrm{A}}}(0,0)$
[Eq. (4.4b)] for a system of linear diabats with
$\kappa_{0}=4,\,\kappa_{1}=-3,\,m=1$ and $\hbar=1$. Panel (a) shows the
function at $\beta=5$, and panel (b) shows the same at $\beta=0.1$, both for
$t_{n}\in[-2,2]$. The darkest blue corresponds to a value of 1, and white
corresponds to 0. The yellow dots highlight some points that make a major
contribution to $\tilde{k}_{{4\mathrm{A}}}$. Half-instanton trajectories
corresponding to these points are plotted directly underneath in panels (c)
and (d). Blue is used for the half-instanton segments that reside on $V_{0}$
and orange is used for $V_{1}$. The dashed black lines depict the system
dynamics on $V_{0}$ before and after the scattering events associated with the
half-instantons.

A similar analysis can be conducted for $\tilde{k}_{{4\mathrm{A}}}$, whose
value at low temperatures is determined by $c_{{4\mathrm{A}}}(z_{0},z_{1},z)$
around $z_{0}=z_{1}=0,$ $z=\tau$ [see Fig. 2(a)]. Expanding
$c_{{4\mathrm{A}}}$ in terms of position eigenstates yields

$\displaystyle c_{{4\mathrm{A}}}(0,0,\tau)=$ (4.11)
$\displaystyle\quad\iint_{-\infty}^{\infty}\\!\\!\mathrm{d}\mathsf{x}_{\mathrm{i}}\mathrm{d}\mathsf{x}_{\mathrm{f}}\left\lVert\braket{\mathsf{x}_{\mathrm{f}}}{\hat{\mathcal{K}}_{0}\\!\left(\tfrac{\tau_{\mathrm{r}}}{4}\right)\\!\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{1}\\!\left(\tfrac{\tau_{\mathrm{p}}}{2}\right)\\!\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{0}\\!\left(\tfrac{\tau_{\mathrm{r}}}{4}\right)\\!}{\mathsf{x}_{\mathrm{i}}}\right\rVert^{2}\\!\\!,$

where the dominant contributions to the integral come from
$\mathsf{x}_{\mathrm{i}}$ and $\mathsf{x}_{\mathrm{f}}$ _both_ in the vicinity
of the reactant turning point. The integrand is once again in the form of a
probability density and describes a double scattering event in which the
system tunnels from the reactant state into the product state, then back to
reactant. This scattering process is associated with the half-instanton
derived from the double-loop trajectory in Fig. 1(b), and is shown
schematically with solid lines in Fig. 2(c). The real-time dynamics
immediately before and after the tunnelling event are indicated on the same
figure with dashed lines.

It is now clear why $\tilde{k}_{{4\mathrm{A}}}$ decreases the full
nonadiabatic rate (recall that it appears in Eq. (3.3) with a minus sign). The
golden-rule expression implicitly assumes that after every reactant-to-product
transition the system remains in the product state. In reality there is a non-
zero probability that the system recrosses back to reactant, which the golden
rule entirely neglects, thus overestimating the full rate constant. The
$\tilde{k}_{{4\mathrm{A}}}$ term accounts for such recrossing events to
leading order in $\hat{\Delta}$, correcting the overestimate. According to our
analysis, this correction is expected to be small at low temperatures, since
the underlying tunnelling mechanism is relatively inefficient compared to the
direct (reactive) transition. For the linear case this can be deduced from the
exponents of Eqs. (4.7) and (4.8), and we will see later that the same applies
to other potentials.

When considering the high-temperature ($\beta\to 0$) limit, our approach has
to be modified, since we can no longer neglect the quartic term in Eq. (4.4b).
The change comes about because the coefficient $d_{01}$ is proportional to
$\beta^{-1}$, whereas the quadratic coefficients $d_{n}$ are linear in $\beta$
[Eq. (4.6)]. This highlights a curious feature of the $c_{{4\mathrm{A}}}$
time-correlation function: unlike its GR counterpart, it cannot be uniformly
approximated as Gaussian in its arguments. Steepest descent integration over
$t$ can always be done in the usual manner, giving us
$L_{{4\mathrm{A}}}^{\mathrm{sc}}(t_{0},t_{1})$, but the remaining two
variables show some unexpected behaviour. At low temperatures,
$L_{{4\mathrm{A}}}^{\mathrm{sc}}$ has an (approximately) Gaussian shape that
broadens with increasing temperature, $T$. Eventually the broadening becomes
noticeably anisotropic and largely confined to the $t_{0}=0$ and $t_{1}=0$
axes. Along all other directions the function shrinks rapidly, so that in the
$\beta\to 0$ limit it assumes the distinctly non-Gaussian four-pointed star
shape in Fig. 2(b). Hence we must modify the steepest-descent procedure
outlined in Sec. II.2 and take the fourth-order term into account. This
results in the asymptotic relation

$\displaystyle\tilde{k}_{{4\mathrm{A}}}(\beta)$
$\displaystyle\sim\frac{\Delta^{4}}{\hbar^{4}}\sqrt{\frac{\pi
m}{2\beta(\kappa_{0}-\kappa_{1})^{2}}}\,Z^{\ddagger}\mkern
1.0mu\mathrm{e}^{-\beta V^{\ddagger}}$ (4.12)
$\displaystyle\quad\quad{}\times\sqrt{\frac{2\pi\hbar}{d_{01}}}\mkern
2.0mu\exp\\!\left(\frac{\widetilde{\beta}^{3}}{48}\right)\mkern-2.0mu\widetilde{K}_{0}\\!\mkern-1.0mu\left(\frac{\widetilde{\beta}^{3}}{32}\right)$

as $\beta\to 0$, where $\widetilde{K}_{0}(y)\equiv\mathrm{e}^{y}K_{0}(y)$ and
$K_{0}$ is a modified Bessel function of the second kind.[87] For one-
dimensional linear diabats this expression is exact at all temperatures. From
the small-argument asymptotic behaviour of $\widetilde{K}_{0}$ it follows that
$\tilde{k}_{{4\mathrm{A}}}$ tends to

$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{ht}}(\beta)=\frac{\pi
m\Delta^{4}Z^{\ddagger}\mkern 1.0mu\mathrm{e}^{-\beta
V^{\ddagger}}}{\hbar^{3}(\kappa_{0}-\kappa_{1})^{2}}\\!\left\\{\ln\\!\left(\frac{64}{\widetilde{\beta}^{3}}\right)-\gamma\right\\}$
(4.13)

as $\beta\to 0$, where $\gamma=0.577\ldots$ is the Euler–Mascheroni constant.
It is absolutely necessary to include the quartic $d_{01}$ term to arrive at
this result.

To understand the change in the asymptotics of $\tilde{k}_{{4\mathrm{A}}}$, we
refer to Fig. 2(b), which shows that $c_{{4\mathrm{A}}}(z_{0},z_{1},z)$ still
has a stationary point at $z_{0}=z_{1}=0$ and $z=\tau$. The corresponding
instanton has the same shape as in the low-temperature regime, but is now
essentially collapsed to a point at the transition state, implying that
nuclear tunnelling does not play a major role. From the same figure we see
that significant contributions to $\tilde{k}_{{4\mathrm{A}}}$ also come from
combinations of arguments where one of $z_{0}$ or $z_{1}$ assumes a non-zero
imaginary value. Without loss of generality we focus on $z_{0}=0$ and
$z_{1}=\mathrm{i}t_{1}$ (the alternative simply interchanges the roles of
$V_{0}$ and $V_{1}$). At this point the $c_{{4\mathrm{A}}}$ correlation
function can be written as

$\displaystyle c_{{4\mathrm{A}}}(0,\mathrm{i}t_{1},\tau)=$ (4.14)
$\displaystyle\iint_{-\infty}^{\infty}\\!\\!\mathrm{d}\mathsf{x}_{\mathrm{i}}\mathrm{d}\mathsf{x}_{\mathrm{f}}\left\lVert\braket{\mathsf{x}_{\mathrm{f}}}{\hat{\mathcal{K}}_{0}\\!\left(\tfrac{\tau_{\mathrm{r}}}{4}\right)\\!\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{1}\\!\left(\tfrac{\tau_{\mathrm{p}}}{2}+\mathrm{i}t_{1}\right)\\!\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{0}\\!\left(\tfrac{\tau_{\mathrm{r}}}{4}\right)\\!}{\mathsf{x}_{\mathrm{i}}}\right\rVert^{2}\\!\\!,$

closely resembling Eq. (4.11). As before, the integrand describes the
probability of an unreactive back-and-forth transition, although now there is
hardly any imaginary-time propagation, since both $\tau_{\mathrm{r}}$ and
$\tau_{\mathrm{p}}$ tend to 0 in the high-temperature limit. The matrix
element in Eq. (4.14) now corresponds to a system starting in the reactant
state and immediately getting scattered into the product state. The product
then follows a real-time trajectory of duration $t_{1}$, at the end of which
it is scattered back into the reactant state. The corresponding half-instanton
follows a path in the complex plane whose imaginary component becomes
vanishingly small as $\beta\to 0$. The real part of the path is shown with
solid lines in Fig. 2(d). On this occasion, the momenta at the end-points of
the half-instanton are non-zero, from which it follows that before and after
the scattering event the system follows the trajectories indicated with dashed
lines.

The change in the asymptotics of $\tilde{k}_{{4\mathrm{A}}}$ is therefore a
consequence of a change in mechanism: at low $T$ the recrossing proceeds
largely via tunnelling, whereas at high $T$ such transitions are mostly due to
high-energy “over-the-barrier” trajectories that overshoot the hopping point
at first and are only scattered on their way back. We note that these
trajectories (emerging naturally from our theory) are precisely what motivates
the Holstein transmission coefficient[88, 50] based on Landau–Zener
theory,[47, 48] which was previously used to construct uniform nonadiabatic
rate expressions.[41, 42, 89] It should then come as no surprise that at high
temperatures, rates based on the Holstein transmission coefficient are
consistent with our Eq. (4.13), as shown in Appendix E. What our new theory
brings to the table is an accurate description of transitions at energies
close to the MECP, whereas LZ only applies in the high-energy limit. As a
result, the Holstein expression can only establish the dependence in Eq.
(4.13) up to an additive constant. Furthermore, this expression only accounts
for over-the-barrier transitions, hence it cannot be used to describe low-
temperature reaction rates, which are dominated by tunnelling under the
barrier. Our new approach can do so rigorously, even for multidimensional
potentials, as shown below.

#### General case.

The high-temperature regime requires additional care because we need to
account for the $\mathcal{O}(t_{0}^{2}t_{1}^{2})$ term in Eq. (4.4a) in order
to obtain the correct semiclassical limit. This can be done with relative ease
for a classical action in its short-time asymptotic form [Eq. (4.3)]. However
deriving the term without relying on the short-time/linear approximation is
not straightforward, since the derivation calls for high-order derivatives of
the action not normally available from standard instanton calculations. At the
same time, our previous discussion indicates that such high-order terms are
negligible at low temperatures, implying that the standard steepest-descent
prescription can be safely followed in this regime. We therefore aim to derive
a semiclassical approximation to Eq. (4.4a) correct to second order in $t_{0}$
and $t_{1}$, expressing all the relevant coefficients in terms of derivatives
of the classical action [Eq. (2.21)]. The quartic term will be based on an
analytically tractable model, yielding a $d_{01}$ coefficient that has the
correct high-temperature limit. It will not necessarily be accurate at low
temperatures, but this should not be an issue provided the coefficient
vanishes sufficiently quickly as $T\to 0$, since in that regime
$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$ is to leading order independent of
$d_{01}$.

To begin, we locate the set of arguments
$\bm{x}=\big{(}\mathsf{x}_{1},\mathsf{x}_{2},\mathsf{x}_{3},\mathsf{x}_{4})$
and $\bm{\tau}=\big{(}\tau_{0},\tau_{1},\tau\big{)}$ that make the combined
action in Eq. (IV.1) stationary. Generally this action will not be known in
closed form, however it can be calculated efficiently using numerical
algorithms.[21, 8, 24] At the stationary point, the derivatives of the action
satisfy

$\frac{\partial{S_{{4\mathrm{A}}}}}{\partial{\bm{x}}}=0,\qquad\frac{\partial{S_{{4\mathrm{A}}}}}{\partial{\bm{\tau}}}=0.$
(4.15)

The trajectory that makes $S_{{4\mathrm{A}}}$ stationary at temperature
$\beta$ is related to the trajectory that makes $S_{2}$ stationary at
temperature $\beta/2$, as shown in Fig. 1. All the information required to
calculate $\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$ can therefore be extracted
from a GR instanton at twice the temperature (that is $\beta/2$). For the
remainder of this section all quantities pertaining to the GR instanton,
namely the hopping point $\mathsf{x}^{\ddagger}$, the imaginary time
$\tau_{\mathrm{g}}$, the stationary action $S_{2}$ and its Hessian
$\mathsf{\Sigma_{2}}$ refer to quantities calculated at this higher
temperature.

Explicitly, the conditions in Eq. (4.15) are satisfied by $\mathsf{x}_{1\ldots
4}=\mathsf{x}^{\ddagger}$, $\tau_{0}=\tau_{1}=0$ and
$\tau=2\tau_{\mathrm{g}}$. It follows that the stationary action
$S_{{4\mathrm{A}}}=2S_{2}$. To relate the derivatives of the two actions, let
us denote $S_{0}\equiv
S_{0}(\mathsf{x}^{\prime\prime},\mathsf{x}^{\prime},\tfrac{\beta\hbar}{2}-\tau_{\mathrm{g}})$
and $S_{1}\equiv
S_{1}(\mathsf{x}^{\prime\prime},\mathsf{x}^{\prime},\tau_{\mathrm{g}})$. We
then define

$\begin{gathered}\mathsf{\Theta}=\frac{\partial^{2}{S_{0}}}{\partial{\mathsf{x}^{\prime}}\partial{\mathsf{x}^{\prime}}}+\frac{\partial^{2}{S_{1}}}{\partial{\mathsf{x}^{\prime\prime}}\partial{\mathsf{x}^{\prime\prime}}},\qquad\mathsf{\Theta}_{n}=\frac{\partial^{2}{S_{n}}}{\partial{\mathsf{x}^{\prime}}\partial{\mathsf{x}^{\prime\prime}}},\\\
\mathsf{w}_{n}=\frac{\partial^{2}{S_{n}}}{\partial{\mathsf{x}^{\prime}}\partial{\tau_{\mathrm{g}}}}=\frac{\partial^{2}{S_{n}}}{\partial{\mathsf{x}^{\prime\prime}}\partial{\tau_{\mathrm{g}}}},\qquad\zeta_{n}=\frac{\partial{{}^{2}S_{n}}}{\partial{\tau_{\mathrm{g}}^{2}}},\end{gathered}$
(4.16)

where $\mathsf{\Theta},\,\mathsf{\Theta}_{n}$ are symmetric
$f\mathord{\times}f$ matrices, $\mathsf{w}_{n}$ are $f$-dimensional column
vectors, $\zeta_{n}$ are scalars, and the derivatives are all evaluated at the
GR stationary point. Thus all of the above quantities are directly available
from a GR instanton calculation and can be evaluated numerically within the
ring-polymer formulation, as described in the Appendix of Ref. 21. It follows
that the Hessian of $S_{{4\mathrm{A}}}$ at its stationary point is the
symmetric matrix111For clarity we only show the lower triangle.

$\mathsf{\Sigma}_{{4\mathrm{A}}}=\frac{\partial^{2}{S_{{4\mathrm{A}}}}}{\partial{\bm{v}\mkern
1.0mu}\partial{\bm{v}}}=\\!\left[\begin{array}[]{*{2}{>{\rule[-0.25em]{0pt}{1.5em}\\!\\!}L{\mycolw}}*{1}{>{\rule[-0.25em]{0pt}{1.5em}}L{0.9em}}
*{3}{>{\rule[-0.25em]{0pt}{1.5em}}L{1.1em}}
>{\rule[-0.25em]{0pt}{1.5em}}L{0.7em}}pt]{0.0pt}{15.00002pt}\\!\\!\
\,2\zeta_{0}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\\!\\!$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright\\\
pt]{0.0pt}{15.00002pt}\\!\\!\phantom{-}0$\@add@raggedright&pt]{0.0pt}{15.00002pt}\\!\\!\
\,2\zeta_{1}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright\\\
pt]{0.0pt}{15.00002pt}\\!\\!\phantom{-}0$\@add@raggedright&pt]{0.0pt}{15.00002pt}\\!\\!\phantom{-}0$\@add@raggedright&pt]{0.0pt}{15.00002pt}\bar{\zeta}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright\\\
pt]{0.0pt}{15.00002pt}\\!\\!{-\mathsf{w}_{0}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\\!\\!{-\mathsf{w}_{1}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\bar{\mathsf{w}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright\\\
pt]{0.0pt}{15.00002pt}\\!\\!\phantom{-}\mathsf{w}_{0}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\\!\\!{-\mathsf{w}_{1}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\bar{\mathsf{w}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}_{\mathrlap{1}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright\\\
pt]{0.0pt}{15.00002pt}\\!\\!\phantom{-}\mathsf{w}_{0}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\\!\\!\phantom{-}\mathsf{w}_{1}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\bar{\mathsf{w}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\,\mathbb{0}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}_{\mathrlap{0}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}$\@add@raggedright&pt]{0.0pt}{15.00002pt}$\@add@raggedright\\\
pt]{0.0pt}{15.00002pt}\\!\\!{-\mathsf{w}_{0}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\\!\\!\phantom{-}\mathsf{w}_{1}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\bar{\mathsf{w}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}_{\mathrlap{0}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\,\mathbb{0}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}_{\mathrlap{1}}$\@add@raggedright&pt]{0.0pt}{15.00002pt}\mathsf{\Theta}$\@add@raggedright\cr\end{array}\right]\\!,$
(4.17)

where $\bm{v}$ is a vector formed by concatenating $\bm{\tau}$ and $\bm{x}$,
$\bar{\mathsf{w}}=\tfrac{1}{2}(\mathsf{w}_{0}+\mathsf{w}_{1})$ and
$\bar{\zeta}=\tfrac{1}{2}(\zeta_{0}+\zeta_{1})$. Next we perform steepest
descent integration over all variables except for $t_{0}$ and $t_{1}$. To this
end we use the relation

$\displaystyle\mathsf{\Sigma}_{{4\mathrm{A}}}$
$\displaystyle\equiv\left[\begin{array}[]{>{\rule[-0.4em]{0pt}{1.6em}}l|l}3.99994pt]{0.0pt}{16.00008pt}\mathsf{M}_{11}&\mathsf{M}_{12}\\\
\hline\cr
3.99994pt]{0.0pt}{16.00008pt}\mathsf{M}_{21}&\mathsf{M}_{22}\end{array}\right]$
(4.18)
$\displaystyle=\left[\begin{array}[]{*{2}{>{\rule[-0.4em]{0pt}{1.6em}$}p{2.2ex}<{$}}}3.99994pt]{0.0pt}{16.00008pt}$\mathbb{1}_{2}&3.99994pt]{0.0pt}{16.00008pt}$\mathsf{B}^{\mkern
0.0mu\top}\\\
3.99994pt]{0.0pt}{16.00008pt}$\mathbb{0}&3.99994pt]{0.0pt}{16.00008pt}$\mathbb{1}_{j}\end{array}\right]\left[\begin{array}[]{>{\rule[-0.4em]{0pt}{1.6em}}C{1em}L{1.5em}}3.99994pt]{0.0pt}{16.00008pt}\mathsf{A}$\@add@centering&\
\mathbb{0}$\@add@raggedright\cr
3.99994pt]{0.0pt}{16.00008pt}\mathbb{0}$\@add@centering&\mathsf{M}_{22}$\<EMAIL_ADDRESS>3.99994pt]{0.0pt}{16.00008pt}$\mathsf{B}&3.99994pt]{0.0pt}{16.00008pt}$\mathbb{1}_{j}\end{array}\right]\\!,$

where $\mathsf{B}=\mathsf{M}_{22}^{-1}\mathsf{M}_{21}$,
$\mathsf{A}\equiv\mathsf{M}_{11}-\mathsf{M}_{12}\mathsf{M}_{22}^{-1}\mathsf{M}_{21}$,
$\mathsf{M}_{11}$ is a $2\mathord{\times}2$ submatrix and the other blocks are
shaped accordingly. Additionally $\mathbb{1}_{k}$ denotes a
$k\mathord{\times}k$ identity matrix and $j=4f+1$. Eq. (4.18) defines a
transformation that casts $\mathsf{\Sigma}_{{4\mathrm{A}}}$ into block-
diagonal form, decoupling $\tau_{n}$ from the other variables. Specifically,

$\begin{gathered}\bm{v}^{\mkern
0.0mu\top}\mathsf{\Sigma}_{{4\mathrm{A}}}\bm{v}=\binom{\tau_{0}}{\tau_{1}}^{\mkern-8.0mu\top}\mathsf{A}\,\binom{\tau_{0}}{\tau_{1}}+\bm{y}^{\mkern
0.0mu\top}\mathsf{M}_{22}\mkern 2.0mu\bm{y}\\\
\bm{y}=\left[\,\mathsf{B}\quad\mathbb{1}_{j}\right]\bm{v},\end{gathered}$
(4.19)

allowing one to integrate over $\bm{y}$ to get a prefactor
$\sqrt{-(2\pi\hbar)^{4f+1}/\lvert\mathsf{M}_{22}\rvert}$ and an effective
action

$S_{{4\mathrm{A}}}^{\mathrm{eff}}=2S_{2}+\frac{1}{2}\binom{z_{0}}{z_{1}}^{\mkern-8.0mu\top}\mathsf{A}\
\binom{z_{0}}{z_{1}}\qquad(\lvert z_{n}\rvert\to 0),$ (4.20)

where $S_{2}$ is evaluated at the stationary point. The second term,
multiplied by $-1/\hbar$, corresponds to the quadratic part of the argument of
the exponential on the second line of Eq. (4.4b). Subsequent calculations are
simplified by the orthogonal transformation

$\mathsf{U}=\left[\begin{array}[]{R{1.2ex}R{1.3em}*{3}{>{\rule[-0.5em]{0pt}{1.7em}}R{2em}}}1$\@add@raggedleft&\mathsf{0}^{\mkern-3.0mu\top}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\mathsf{0}^{\mkern-3.0mu\top}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\mathsf{0}^{\mkern-3.0mu\top}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\mathsf{0}^{\mkern-3.0mu\top}$\@add@raggedleft\cr\mathsf{0}$\@add@raggedleft&\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\frac{1}{2}\mathbb{1}$\@add@raggedleft\cr\mathsf{0}$\@add@raggedleft&\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}-\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}-\frac{1}{2}\mathbb{1}$\@add@raggedleft\cr\mathsf{0}$\@add@raggedleft&\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}-\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}-\frac{1}{2}\mathbb{1}$\@add@raggedleft\cr\mathsf{0}$\@add@raggedleft&\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}-\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}-\frac{1}{2}\mathbb{1}$\@add@raggedleft&0pt]{0.0pt}{16.99998pt}\frac{1}{2}\mathbb{1}$\@add@raggedleft\end{array}\right]$
(4.21)

(here $\mathbb{1}\equiv\mathbb{1}_{f}$), which can be used to show that

$\begin{gathered}\lvert\mathsf{M}_{22}\rvert=\lvert\mathsf{U}^{\mkern
0.0mu\top}\mathsf{M}_{22}\mathsf{U}\rvert=\frac{1}{2}\lvert\mathsf{\Sigma}_{2}\rvert\lvert\mathsf{\Theta_{+}}\rvert\lvert\mathsf{\Theta_{-}}\rvert,\\\
\mathsf{\Theta}_{+}=\mathsf{\Theta}+\mathsf{\Theta}_{0}-\mathsf{\Theta}_{1},\qquad\mathsf{\Theta}_{-}=\mathsf{\Theta}-\mathsf{\Theta}_{0}+\mathsf{\Theta}_{1},\end{gathered}$
(4.22)

with $\mathsf{\Sigma}_{2}$ defined in Eq. (2.27). By the same token,

$\displaystyle\mathsf{A}=-\\!\left[\begin{array}[]{*{2}{>{\rule[-0.2em]{0pt}{1.2em}}C{1em}}}99997pt]{0.0pt}{11.99998pt}d_{0}$\@add@centering&99997pt]{0.0pt}{11.99998pt}0$\@add@centering\cr
99997pt]{0.0pt}{11.99998pt}0$\@add@centering&99997pt]{0.0pt}{11.99998pt}d_{1}$\@add@centering\end{array}\right],$
(4.25) $\displaystyle\rule{0.0pt}{11.99998pt}d_{0}=4\mathsf{w}_{0}^{\mkern
0.0mu\top}\mathsf{\Theta}_{+}^{-1}\mathsf{w}_{0}-2\zeta_{0},\quad
d_{1}=4\mathsf{w}_{1}^{\mkern
0.0mu\top}\mathsf{\Theta}_{-}^{-1}\mathsf{w}_{1}-2\zeta_{1},$

which is the general form of the coefficients $d_{n}$. With this we have all
the information necessary to calculate the $\beta\to\infty$ limit of
$\tilde{k}_{{4\mathrm{A}}}$ for a general, multidimensional potential.

It remains to connect this and the $\beta\to 0$ limit in Eq. (4.12), for which
we require the coefficient $d_{01}$ of the $z_{0}^{2}z_{1}^{2}$ term in the
short-time expansion of $S_{{4\mathrm{A}}}^{\mathrm{eff}}$. This calls for
high-order derivatives of $S_{2}$ not readily available from standard GR
instanton calculations. We can circumvent this by noting that, based on Eq.
(4.6) and on general physical considerations, the ratio $\hbar
d_{01}/d_{0}d_{1}$ tends to zero as $\beta\to\infty$. In this limit the factor
arising from integration over $t_{0}$ and $t_{1}$ behaves as

$\sqrt{\frac{2\pi\hbar}{d_{01}}}\
\widetilde{K}_{0}\\!\left(\frac{d_{0}d_{1}}{4d_{01}\hbar}\right)\sim\frac{2\pi\hbar}{\sqrt{d_{0}d_{1}}}\left(1-\frac{\hbar}{2}\frac{d_{01}}{d_{0}d_{1}}\right)\\!,$
(4.26)

_i.e._ , to leading order it is independent of $d_{01}$. Therefore any
reasonable approximation to $d_{01}$ should lead to a good estimate of
$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$, provided the expression satisfies
$\hbar d_{01}/d_{0}d_{1}\ll 1$ as $\beta\to\infty$, and tends to Eq. (4.6) in
the high-temperature limit. We suggest the following expression, which is
exact for a one-dimensional spin–boson model (that is, Eq. (5.5) with $f=1$):

$\displaystyle
d_{01}=\mathrm{max}\bigg{[}0,\,\tfrac{\omega}{2m}\bigg{\\{}(\kappa_{0}-\kappa_{1})^{2}\operatorname{csch}\\!\left(\tfrac{\beta\hbar\omega}{2}\right)$
(4.27)
$\displaystyle\quad{}+\kappa_{0}\kappa_{1}\operatorname{sech}\\!\left(\tfrac{\beta\hbar\omega}{4}\right)^{\\!2}\left[\left(\tfrac{\kappa_{0}+\kappa_{1}}{\kappa_{0}-\kappa_{1}}\right)^{\mathrlap{\\!2}}+\operatorname{csch}\\!\left(\tfrac{\beta\hbar\omega}{4}\right)^{\\!2}\right]^{-\frac{1}{2}}\\!\bigg{\\}}\bigg{]}.$

Here $\kappa_{n}$ are the signed norms of $\nabla V_{n}$ at the MECP and
$\omega>0$ is a parameter ensuring that $d_{01}$ vanishes sufficiently quickly
as $\beta\to\infty$. A suitable value can be obtained from the curvature along
the direction of steepest descent at the MECP,
$\omega^{2}=m^{-1}\nabla^{2}_{q}\,V_{n}(\mathsf{x}^{\ddagger})$, calculated
for whichever diabatic PES yields the larger value. Unlike Eq. (4.6), this
expression for $d_{01}$ can become negative at sufficiently low temperatures,
tending to zero from below as $\beta\to\infty$. The most straightforward
course of action in such cases is to set $d_{01}$ to zero, which is equivalent
to neglecting contributions from the quartic term.222This is also the more
consistent treatment, since at temperatures low enough that $d_{01}<0$, the
corresponding contribution to the steepest-descent integral is small and no
longer dominates the other quartic terms ($t^{4}$ and $t_{n}^{4}$), which are
neglected. There is therefore no reason to give the $d_{01}$ contribution
privileged treatment and include it into the final result.

More sophisticated approximations to $d_{01}$ could be constructed, but any
modification is expected to only have a significant effect at temperatures
where $\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$ is subdominant to the other
components, namely $\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}$, as we show
later. Thus we arrive at the final expression for the rate constant,

$\displaystyle\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}(\beta)=\left[\tilde{k}_{2}^{\mathrm{sc}}\mkern-1.0mu\big{(}\tfrac{\beta}{2}\big{)}\right]^{2}\\!\\!\sqrt{\frac{\mathllap{-}\lvert\mathsf{\Sigma}_{2}\rvert}{2d_{01}\lvert\mathsf{\Theta}_{+}\rvert\lvert\mathsf{\Theta}_{-}\rvert}}\,\widetilde{K}_{0}\\!\left(\frac{d_{0}d_{1}}{4\mkern
1.0mud_{01}\hbar}\right)\\!,$ (4.28)

where $\tilde{k}_{2}^{\mathrm{sc}}$ is defined in Eq. (2.28). The expression
is a uniform approximation that is expected to hold at any temperature.
Additionally, the $\omega\to 0$ limit of Eq. (4.27) recovers Eq. (4.6), and so
our current prescription for $d_{01}$ yields the exact
$\tilde{k}_{{4\mathrm{A}}}$ for two-state systems comprised of one-dimensional
linear diabats.

###  IV.2 Semiclassical _B_ -type contribution

The considerations in Sec. III.2 allow us to write

$\displaystyle\tilde{k}_{{4\mathrm{B}}}\sim\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}+{}$
(4.29)
$\displaystyle\operatorname{Re}\Bigg{\\{}\\!\left[\iint_{\mathcal{A}_{\mathrlap{\tau\tau}}}+\iint_{\mathcal{A}_{\mathrlap{\tau
t}}}+\iint_{\mathcal{A}_{\mathrlap{t\tau}}}+\iint_{\mathcal{A}_{tt}}\right]L_{{4\mathrm{B}}}^{\mathrm{sc}}(z_{0},z_{1})\,\mathrm{d}z_{0}\mathrm{d}z_{1}\Bigg{\\}},$

where $L_{{4\mathrm{B}}}^{\mathrm{sc}}$ is obtained by substituting
semiclassical propagators into Eq. (3.12) and integrating over positions and
time $t$ by steepest descent. As before, we identify the combined action,

$\displaystyle S_{{4\mathrm{B}}}$ $\displaystyle\equiv
S_{0}\big{(}\mathsf{x}_{1},\mathsf{x}_{4},z_{\mathrm{r}}-z_{0}\big{)}+S_{1}\big{(}\mathsf{x}_{4},\mathsf{x}_{3},z_{1}\big{)}$
$\displaystyle{}+S_{0}\big{(}\mathsf{x}_{3},\mathsf{x}_{2},z_{0}\big{)}+S_{1}\big{(}\mathsf{x}_{2},\mathsf{x}_{1},z_{\mathrm{p}}-z_{1}\big{)},$
(4.30)

which we expect to be minimised by a two-bounce (single-loop) trajectory like
the one shown in Fig. 4. To prove this we need to determine the optimal
integration domains (_i.e._ , the optimal $\tau_{n}^{*}$) which is most easily
done in the high-temperature limit.

#### High-temperature limit.

Once again, we expand the diabatic PESs in a Taylor series about the MECP and
consider the corresponding propagators defined by Eq. (4.3). Integration over
position and time $t$ gives

$\displaystyle
L_{{4\mathrm{B}}}^{\mathrm{sc}}(z_{0},z_{1})\sim\sqrt{\frac{2\pi
m}{\beta\hbar^{2}}}\frac{2\Delta^{4}}{\hbar^{3}\lvert\kappa_{0}-\kappa_{1}\rvert}Z^{\ddagger}\mkern
1.0mu\mathrm{e}^{-\beta V^{\ddagger}}$ (4.31) $\displaystyle\ \
{}\times\exp\left\\{-\frac{1}{\hbar}\left[\frac{(\kappa_{0}-\kappa_{1})^{2}}{2m\beta\hbar}z_{0}^{2}z_{1}^{2}-\frac{\beta\hbar\mkern
1.0mu\kappa_{0}\kappa_{1}}{2m}z_{0}z_{1}\right.\right.$
$\displaystyle\qquad\left.\left.{}-\frac{\kappa_{0}-\kappa_{1}}{2m}(\kappa_{0}z_{0}^{2}z_{1}-\kappa_{1}z_{1}^{2}z_{0})-\frac{\beta^{3}\hbar^{3}\kappa_{0}^{2}\kappa_{1}^{2}}{24m(\kappa_{0}-\kappa_{1})^{2}}\right]\right\\}.$

where all the terms are defined as in Sec. IV.1. From here it can be shown
that the remaining integrals over $z_{n}$ in Eq. (4.29) can be performed most
easily by setting

$\tau_{0}^{*}=\frac{\beta\hbar}{2}\frac{\kappa_{1}}{\kappa_{1}-\kappa_{0}},\qquad\tau_{1}^{*}=\frac{\beta\hbar}{2}\frac{\kappa_{0}}{\kappa_{0}-\kappa_{1}}.$
(4.32)

This choice removes any oscillations in the integrand and allows for further
simplifications, which will be discussed shortly. Note that these values of
$\tau_{n}^{*}$ are related to the stationary imaginary time of a golden-rule
instanton in a linear/high-temperature system at temperature $\beta$, namely
$\tau_{1}^{*}=\tau_{\mathrm{g}}/2$ and
$\tau_{0}^{*}=(\beta\hbar-\tau_{\mathrm{g}})/2$. Under this definition,
contributions from $\mathcal{A}_{\tau t}$ and $\mathcal{A}_{t\tau}$ are purely
imaginary and are therefore discarded. Furthermore, the contribution from
$\mathcal{A}_{tt}$ cancels exactly with
$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$ in Eq. (4.29), meaning that the
semiclassical rate is

$\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}\sim\int_{0}^{\tau_{0}^{*}}\mathrm{d}\tau_{0}\int_{0}^{\tau_{1}^{*}}\mathrm{d}\tau_{1}\,L_{{4\mathrm{B}}}(\tau_{0},\tau_{1})$
(4.33)

as $\beta\to 0$. We observe that similarly to $\tilde{k}_{{4\mathrm{A}}}$, the
$\tilde{k}_{{4\mathrm{B}}}$ rate cannot be approximated with a Gaussian
integral. Also like with $\tilde{k}_{{4\mathrm{A}}}$, the high-temperature
limit of the $\tilde{k}_{{4\mathrm{B}}}$ rate given by Eqs. (4.31)–(4.33) is
exact across all temperatures for the linear system. Hence we may use the
results of this section to inspire the derivation in the general case.

(a)

(b)

Figure 3: Panel (a) shows the scaled integrand of Eq. (4.33),
$L_{{4\mathrm{B}}}(\tau_{0},\tau_{1})/L_{{4\mathrm{B}}}(0,0)$, for linear
diabats with $\kappa_{0}=4$, $\kappa_{1}=-3$, $m=1$, $\hbar=1$. The plotting
ranges are $\tau_{n}\in[0,\tau_{n}^{*}]$ as per Eq. (4.32). The darkest blue
corresponds to $1$, and white corresponds to $0$. For $\beta=4$ we highlight a
representative set of $(\tau_{0},\tau_{1})$ values in yellow, showing the
corresponding minimal action instanton trajectories in panel (b). The $x$-axis
marks the imaginary time along the instanton, and the $y$-axis shows the
position relative to the hopping point $\mathsf{x}^{\ddagger}$. Blue and
orange trajectory segments correspond to imaginary-time propagation along the
reactant and product diabats respectively.

In Fig. 3(a) we plot $L_{{4\mathrm{B}}}$ at different temperatures for
$(\tau_{0},\tau_{1})\in\mathcal{A}_{\tau\tau}$. As $\beta\to 0$, the function
becomes flat across the entire domain of integration; as $\beta\to\infty$, the
function becomes dominated by the edges for which $\tau_{0}$ or $\tau_{1}$ is
$0$. In Fig. 3(b) we show the instanton trajectories that correspond to
different points inside the domain of integration. As anticipated at the
beginning of Sec. III, trajectories that dominate the integral at low
temperatures have two consecutive segments that bounce once, and two
consecutive segments that travel directly between the endpoints. Overall,
these follow the path of a golden-rule instanton, with either a short reactant
segment along a product trajectory or vice versa. As we move on to
trajectories for which neither $\tau_{n}$ is small, the segments become
elongated and the entire instanton more localised around the barrier region.
In the most extreme case (top right trajectory), we arrive at the four-bounce
instanton that corresponds to the $\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$
term in Sec. IV.1. For large $\beta$ these trajectories have considerably
larger action and do not contribute much to the integral in Eq. (4.33). Hence
we expect that the general semiclassical expression for
$\tilde{k}_{{4\mathrm{B}}}$ will be dominated by trajectories for which one of
the direct segments is vanishingly short (one of $\tau_{n}\to 0$). All such
trajectories reduce exactly to a standard (two-bounce) GR instanton at a
temperature $\beta$.

#### General case

In order to derive a general expression for
$\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}$ we need to express the coefficients
in

$\displaystyle L_{{4\mathrm{B}}}^{\mathrm{sc}}(\tau_{0},\tau_{1})\sim
L_{{4\mathrm{B}}}^{\mathrm{sc}}(\tau_{0},0)\exp\\!\left\\{\\!{-\frac{1}{\hbar}}\\!\left[\alpha_{1}\tau_{1}-\frac{\gamma_{1}\tau_{1}^{2}}{2}\right]\\!\right\\}$
(4.34)

(valid as $\tau_{1}\to 0$) in terms of derivatives of the GR action. From Eqs.
(3.11) and (3.12) it follows that

$L_{{4\mathrm{B}}}^{\mathrm{sc}}(\tau_{0},0)=\frac{2\mkern
1.0mu\tilde{k}_{2}^{\mathrm{sc}}(\beta)}{\hbar^{2}}\Delta[\mathsf{x}(\tau_{0})]^{2},$
(4.35)

where $\Delta[\mathsf{x}(\tau_{0})]$ is the diabatic coupling along the
reactant segment of a GR instanton after travelling for an imaginary time
$\tau_{0}$ away from the hopping point. Deriving expressions for $\alpha_{1}$
and $\gamma_{1}$ (both functions of $\tau_{0}$) is somewhat laborious, but the
underlying approach is similar to the one used in Sec. IV.1 to derive $d_{n}$.
Treating $\tau_{0}$ as a parameter, we consider a series expansion of
$S_{{4\mathrm{B}}}$ in the remaining six variables and integrate over
$\left(t,\,\mathsf{x}_{1},\,\mathsf{x}_{2},\,\mathsf{x}_{3},\,\mathsf{x}_{4}\right)$
by steepest descent. The resulting effective action is in the form of the
bracketed expression in Eq. (4.34), and so yields the desired coefficients. To
begin, we transform to a more convenient set of coordinates,

$\displaystyle\mathsf{x}^{\prime}$ $\displaystyle\equiv\mathsf{x}_{1},$
$\displaystyle\quad\mathsf{x}^{\prime\prime}$
$\displaystyle\equiv\mathsf{x}_{2},$ (4.36) $\displaystyle\mathsf{x}$
$\displaystyle\equiv(\mathsf{x}_{3}+\mathsf{x}_{4})/2,$
$\displaystyle\quad\mathsf{x}_{-}$
$\displaystyle\equiv\mathsf{x}_{4}-\mathsf{x}_{3},$

as indicated in Fig. 4. We then define

$\displaystyle S_{0}^{\prime}$ $\displaystyle\equiv
S_{0}(\mathsf{x}\mathrlap{{}^{\prime}},\ \mathsf{x},\
\beta\hbar-\tau-\tau_{0}),$ (4.37a) $\displaystyle S_{0}^{\prime\prime}$
$\displaystyle\equiv S_{0}(\mathsf{x},\
\mathsf{x}\mathrlap{{}^{\prime\prime}},\ \tau_{0}),$ (4.37b) $\displaystyle
S_{1}$ $\displaystyle\equiv S_{1}(\mathsf{x}\mathrlap{{}^{\prime\prime}},\
\mathsf{x}\mathrlap{{}^{\prime}},\ \tau),$ (4.37c)

and introduce the notation $\bm{v}\equiv\left(\tau\ \mathsf{x}^{\prime}\
\mathsf{x}\ \,\mathsf{x}^{\prime\prime}\right)$. Treating $\tau_{0}$ as a
parameter, we expand $S_{{4\mathrm{B}}}$ to second order about $\tau_{1}=0$,
$\mathsf{x}_{-}=\mathsf{0}$, and a $\bm{v}$ such that

$\frac{\partial{}}{\partial{\bm{v}}}\left(S_{0}^{\prime}+S_{0}^{\prime\prime}+S_{1}\right)=\bm{0}.$
(4.38)

Figure 4:  An instanton trajectory with a significant contribution to the
integral in Eq. (4.33). Here $\tau_{\mathrm{r}}=\beta\hbar-\tau_{\mathrm{g}}$,
$\tau_{\mathrm{p}}=\tau_{\mathrm{g}}$ and $\tau_{\mathrm{g}}$ is the
stationary imaginary time of a GR instanton at temperature $\beta$.

We then integrate the (exponential of) the resulting expression over all
variables except for $\tau_{0}$ and $\tau_{1}$, as discussed in Appendix A.
This leads to the desired Taylor series expansion of the effective action,
with coefficients

$\displaystyle\rule[-7.5pt]{0.0pt}{10.00002pt}\alpha_{1}$
$\displaystyle=V_{1}(\mathsf{x})-V_{0}(\mathsf{x}),$ (4.39a)
$\displaystyle\gamma_{1}$
$\displaystyle=\frac{\partial{(V_{1}-V_{0})}}{\partial{\mathsf{x}}}\cdot\left[\mathsf{\Phi}_{0}^{-1}\right]_{\mathsf{xx}}\\!\cdot\frac{\partial{(V_{1}-V_{0})}}{\partial{\mathsf{x}}}$
(4.39b)
$\displaystyle+\frac{\partial{(V_{1}-V_{0})}}{\partial{\mathsf{x}}}\cdot\left[\mathsf{\Phi}_{0}^{-1}\right]_{\mathsf{x}\bm{v}}\\!\cdot\frac{\partial^{2}{(S_{0}^{\prime}+S_{0}^{\prime\prime})}}{\partial{\bm{v}}\partial{\tau_{0}}},$

where $\mathsf{x}$ is the coordinate along the reactant segment of the GR
instanton at time $\tau_{0}$ and

$\mathsf{\Phi}_{0}=\frac{\partial^{2}{(S_{0}^{\prime}+S_{0}^{\prime\prime}+S_{1})}}{\partial{\bm{v}}\partial{\bm{v}}}.$
(4.40)

We use $\left[\mathsf{\Phi}_{0}^{-1}\right]_{\mathsf{x}\bm{v}}$ to denote the
matrix formed from $\mathsf{\Phi}_{0}^{-1}$ by taking the rows corresponding
to $\mathsf{x}$ and all of the columns;
$\left[\mathsf{\Phi}_{0}^{-1}\right]_{\mathsf{x}\mathsf{x}}$ is formed from
the rows _and_ columns corresponding to $\mathsf{x}$. The same considerations
apply to the asymptotic form of $L_{{4\mathrm{B}}}(\tau_{0},\tau_{1})$ as
$\tau_{0}\to 0$, and expressions for $\alpha_{0}$ and $\gamma_{0}$ can be
obtained by swapping the indices $0\leftrightarrow 1$ in the preceding
derivation.

It remains to evaluate the double integral in Eq. (4.33) over the rectangular
domain $\mathcal{A}_{\tau\tau}$, defined by
$\tau^{*}_{1}=\tfrac{1}{2}\tau_{\mathrm{g}}(\beta)$,
$\tau^{*}_{0}=\tfrac{1}{2}[\beta\hbar-\tau_{\mathrm{g}}(\beta)]$. To do this,
we split the domain along the diagonal connecting the bottom left and top
right corners of the rectangles in Fig. 3(a),

$\displaystyle\int_{0}^{\tau^{*}_{0}}\mathrm{d}\tau_{0}\int_{0}^{\tau^{*}_{1}}\mathrm{d}\tau_{1}={}$
(4.41)
$\displaystyle\qquad\qquad\int_{0}^{\tau^{*}_{0}}\mathrm{d}\tau_{0}\int_{0}^{\frac{\tau^{*}_{1}}{\tau^{*}_{0}}\mathrlap{\tau_{0}}}\mathrm{d}\tau_{1}+\int_{0}^{\tau^{*}_{1}}\mathrm{d}\tau_{1}\int_{0}^{\frac{\tau^{*}_{0}}{\tau^{*}_{1}}\mathrlap{\tau_{1}}}\mathrm{d}\tau_{0}.$

We shall focus on the first of the two terms; the result for the second term
follows upon exchanging the labels $0\leftrightarrow 1$. The outer integral
over $\tau_{0}$ is not amenable to steepest descent integration but can
readily be evaluated numerically using information directly available from
standard ring-polymer instanton optimisation. The inner integral is of the
form

$j_{0}^{\mathrm{sc}}(\tau_{0})=\int_{0}^{\frac{\tau^{*}_{1}}{\tau^{*}_{0}}\mathrlap{\tau_{0}}}\mathrm{d}\tau_{1}\,\mathrm{e}^{{-\frac{1}{\hbar}}\\!\left[\alpha_{1}(\tau_{0})\tau_{1}-\tfrac{1}{2}\gamma_{1}(\tau_{0})\tau_{1}^{2}\right]},$
(4.42)

where $\alpha_{1}$ and $\gamma_{1}$ are defined in Eq. (4.39). The standard
procedure (82, see Chapter 6.4 of Ref. ) for calculating this integral by
steepest descent is to replace the integrand with

$\exp\\!\left[-\frac{\alpha_{1}\tau_{1}}{\hbar}+\frac{\gamma_{1}\tau_{1}^{2}}{2\hbar}\right]\sim\mathrm{e}^{-\alpha_{1}\tau_{1}/\hbar}\left[1+\frac{\gamma_{1}\tau_{1}^{2}}{2\hbar}\right]$
(4.43)

and integrate the resulting expression analytically. This yields

$\displaystyle j_{0}^{\mathrm{sc}}(\tau_{0})$
$\displaystyle=\frac{\hbar}{\alpha_{1}}\left[1-\mathrm{e}^{-\phi_{1}}\right]$
(4.44)
$\displaystyle\quad{}+\frac{\hbar^{2}\gamma_{1}}{\alpha^{3}_{1}}\left[1-\mathrm{e}^{-\phi_{1}}\\!\left(\tfrac{\phi_{1}^{2}}{2}+\phi_{1}+1\right)\right],$

with
$\phi_{1}=\tfrac{\alpha_{1}}{\hbar}\mkern-1.0mu\tfrac{\tau_{1}^{*}}{\tau_{0}^{*}}\tau_{0}$,
and hence the semiclassical rate

$\displaystyle\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}(\beta)=\frac{2\tilde{k}_{2}^{\mathrm{sc}}(\beta)}{\hbar^{2}}\bigg{\\{}$
$\displaystyle\int_{0}^{\tau^{*}_{0}}\\!\\!\Delta[\mathsf{x}_{\mathrm{r}}(\tau_{0})]^{2}\,j_{0}^{\mathrm{sc}}(\tau_{0})\,\mathrm{d}\tau_{0}$
(4.45) $\displaystyle{}+$
$\displaystyle\int_{0}^{\tau^{*}_{1}}\\!\\!\Delta[\mathsf{x}_{\mathrm{p}}(\tau_{1})]^{2}\,j_{1}^{\mathrm{sc}}(\tau_{1})\,\mathrm{d}\tau_{1}\bigg{\\}},$

where $\Delta[\mathsf{x}_{\mathrm{r}}(\tau_{0})]$ is the diabatic coupling at
time $\tau_{0}$ along the reactant segment of the GR instanton,
$\mathsf{x}_{\mathrm{r}}(0)$ is the hopping point,
$\mathsf{x}_{\mathrm{r}}(\tau_{0}^{*})$ is the turning point, and analogous
definitions hold for $\mathsf{x}_{\mathrm{p}}$. Despite the lengthy
derivation, Eq. (4.45) is easy to evaluate in practice. All one needs are the
potentials $V_{n}$, their gradient
$\frac{\partial{V_{n}}}{\partial{\mathsf{x}}}$ and the diabatic coupling
$\Delta$ along the GR instanton trajectory at temperature $\beta$, together
with the derivatives of the stationary action directly available from a
standard GR calculation. No reoptimisation is required, and the additional
computational overhead is minimal.

There are three points remaining to be addressed. First, Eq. (4.34) with
$\alpha_{n},\,\gamma_{n}$ defined as in Eq. (4.39) gives the exact
$L_{{4\mathrm{B}}}$ (at all $\tau_{n}$) for linear diabats
$V_{n}=V^{\ddagger}+\kappa_{n}x$. Therefore for that system we need not
perform the expansion in Eq. (4.43) and can instead evaluate Eq. (4.42)
directly, as described in Appendix B. This leads us to the exact
$\tilde{k}_{{4\mathrm{B}}}(\beta)$ for a system of linear diabats, which we
can compare with the asymptotic approximation in Eq. (4.45). It can be shown
that for any combination of parameters the latter tends to the exact result in
both the high- and low-temperature limits, and is never in error by more than
$2.8\%$.

Second, we should stress that the rectangular domains in Fig. 3, which we have
adopted throughout the derivation, are only strictly optimal in the high-
temperature limit (unless the system is linear, in which case the domains are
optimal at any temperature). It is also only in those limits that
contributions from $\mathcal{A}_{\tau t}$ and $\mathcal{A}_{t\tau}$ are purely
imaginary, and therefore discarded, and only there that the integral over
$\mathcal{A}_{tt}$ cancels exactly with
$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$. Nevertheless, any significant
deviations are expected to occur at low temperatures, for which the integral
is entirely dominated by regions of small $\tau_{0}$ or $\tau_{1}$.
Deformation of the optimal $\mathcal{A}_{\tau\tau}$ from the rectangular shape
would only become prominent in parts of the $(\tau_{0},\tau_{1})$ plane that
contribute little to the overall integral. Furthermore, we expect terms
derived from integration over optimal $\mathcal{A}_{\tau t}$,
$\mathcal{A}_{t\tau}$ and $\mathcal{A}_{tt}$, as well as the rate constant
$\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}$, to become subdominant to
$\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}$ in the low-temperature limit.
Therefore any real part acquired by the first two terms, or any imperfect
cancellation between the latter are expected to be negligible, and thus Eq.
(4.45) remains the correct asymptotic limit of the $\tilde{k}_{{4\mathrm{B}}}$
rate constant.

Last comes the question of physical interpretation. Mirroring the discussion
in Sec. IV.1, we consider $c_{{4\mathrm{B}}}(z_{0},z_{1},z)$ at a point that
makes a significant contribution to $\tilde{k}_{{4\mathrm{B}}}$, expressing
the function in the position basis. In what follows we consider the
contribution that corresponds to the instanton in Fig. 4, for which
$z_{0}=\tau_{0}\in(0,\tfrac{\tau_{\mathrm{r}}}{2})$, $z_{1}\to 0$ and
$z=\tau$. This can be written as

$\displaystyle
c_{{4\mathrm{B}}}(\tau_{0},0,\tau)=\iint_{-\infty}^{\infty}\\!\\!\mathrm{d}\mathsf{x}_{\mathrm{i}}\mathrm{d}\mathsf{x}_{\mathrm{f}}\braket{\mathsf{x}_{\mathrm{i}}}{\hat{\mathcal{K}}_{0}\\!\left(\tfrac{\tau_{\mathrm{r}}}{2}\right)\\!\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{1}\\!\left(\tfrac{\tau_{\mathrm{p}}}{2}\right)\\!}{\mathsf{x}_{\mathrm{f}}}\times{}$
$\displaystyle\braket{\mathsf{x}_{\mathrm{f}}}{\hat{\mathcal{K}}_{1}\\!\left(\tfrac{\tau_{\mathrm{p}}}{2}\right)\\!\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{0}(\tau_{0})\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{1}(0)\hat{\Delta}\mkern
2.0mu\hat{\mathcal{K}}_{0}\\!\left(\tfrac{\tau_{\mathrm{r}}}{2}-\tau_{0}\right)\\!}{\mathsf{x}_{\mathrm{i}}}.$
(4.46)

Unlike in Eqs. (4.11) and (4.14), here the integrand is not the square of a
single matrix element and hence cannot be directly identified with the
probability of a scattering event. Instead, the integrand takes the form of a
quantum interference term. The first of the two interfering alternatives is
described by the matrix element on the first line of Eq. (IV.2) and
corresponds to a golden-rule transition, as follows from comparing to Eq.
(4.10). The second alternative, corresponding to the second line of Eq.
(IV.2), is a new type of scattering process that involves an excitation into a
_virtual_ state.

Figure 5:  Schematic depiction of the half-instanton that corresponds to
$c_{{4\mathrm{B}}}(\tau_{0},0,\tau)$ in Eq. (IV.2). The half-instanton
segments residing on $V_{0}$ are shown in blue and those on $V_{1}$ are shown
in orange. The wavy line indicates an excitation into a “virtual” state, which
is infinitesimally short and not energy-conserving. The dashed black lines
depict the system dynamics before and after the scattering event associated
with the half-instanton.

The half-instanton trajectory corresponding to this process is shown in Fig. 5
and is identical to the GR case, except for the infinitesimally short
excitation onto $V_{1}$. We call this excitation virtual because, firstly, the
amount of time spent on $V_{1}$ is vanishingly small and, secondly, the
excitation does not conserve the energy of the half-instanton. As $\tau_{0}$
is varied (_e.g._ , in the first integral in Eq. (4.45)), the excitation moves
from the turning point to the hopping point, “sampling” $V_{1}$ all along the
reactant segment of the half-instanton. When $\tau_{0}\to 0$ and
$\tau_{1}\in(0,\tfrac{\tau_{\mathrm{p}}}{2})$, as in the second integral in
Eq. (4.45), the excitation moves along the product segment and samples $V_{0}$
instead.

It seems reasonable that this type of scattering process should describe how
the _upper_ diabatic state influences the probability of tunnelling from
reactant to product. Casting the system in the adiabatic representation, we
see that the shape of the lower Born–Oppenheimer PES depends on $\Delta$ and
_both_ of the diabats $V_{n}$ all along the tunnelling pathway (not just at
the hopping point). $\tilde{k}_{{4\mathrm{B}}}$ takes this into account and
serves to increase the overall reaction rate. This is consistent with the
notion that a larger $\Delta$ will tend to lower the (adiabatic) activation
energy barrier, increasing the probability of tunnelling and hence the
reaction rate.

###  IV.3 Semiclassical partition function correction

It is often the case that the diabatic coupling $\Delta(\mathsf{x})$ becomes
negligible as $\mathsf{x}$ approaches the minimum of $V_{n}$. The term
considered in this section is then itself negligible and can be omitted
altogether. If the coupling remains considerable around the reactant/product
minimum, it can be accounted for as follows. We begin by rewriting Eq. (3.18)
as

$\displaystyle\hbar^{2}\widetilde{Z}_{\mathrm{r},2}$
$\displaystyle=\int_{0}^{\frac{\beta\hbar}{2}}\\!\\![\beta\hbar-2u]\,c_{2}(u)\,\mathrm{d}u$
(4.47)
$\displaystyle-\bigg{[}\int_{0}^{\frac{\beta\hbar}{2}}\\!\\!u\,c_{2}(u)\,\mathrm{d}u+\int_{0}^{\frac{\beta\hbar}{2}}\\!\\!u\,c_{2}(\beta\hbar-u)\,\mathrm{d}u\bigg{]}\frac{Z_{\mathrm{r},0}}{Z_{\mathrm{p},0}}.$

For imaginary times $u\in[0,\beta\hbar]$, the only stationary point of
$c_{2}(u)$ is a minimum at $u=\tau_{\mathrm{g}}(\beta)$, which corresponds to
the GR instanton. In the absence of maxima, the asymptotic expansion of Eq.
(4.47) is determined by how $c_{2}(u)$ behaves near the edges of the
integration domains, namely $u=0$ and $u=\beta\hbar$. These correspond to
instantons collapsed at the bottom of the reactant and the product well
respectively.

In principle we should also consider $u=\beta\hbar/2$, but this contribution
can be neglected unless the system is in the near-activationless regime or the
temperature is so low that the value of $c_{2}(\beta\hbar/2)$ becomes
comparable to (or exceeds) one of $c_{2}(0)$ or $c_{2}(\beta\hbar)$, which is
possible in asymmetric systems. Other than is such cases, our approximation
holds because the instanton trajectory associated with $c_{2}(\beta\hbar/2)$
is localised in the barrier region. It will thus resemble a GR instanton and
give rise to a relatively large stationary action, making its contribution
negligible compared to those from the reactant and product wells. This is
shown in Fig. 6, where on panel (a) we plot a representative $c_{2}(u)$,
normalised by $\widetilde{Z}_{\mathrm{r},0}$. Panel (b) shows the natural
logarithm of $c_{2}(u)$—in this case simply the stationary action at a given
$u$, multiplied by $-1/\hbar$.

Figure 6:  (a) The exact GR flux-correlation function for a 1D spin–boson
model [Eq. (5.5)] at $\beta=2$, $\epsilon=2$, $d=3$ and all other parameters
equal to $1$. The function is normalised by $\widetilde{Z}_{\mathrm{r},0}$.
Note the change of scale either side of the vertical line at $u=\beta\hbar/2$.
The blue dotted line overlaying the black curve is the asymptotic
approximation as $u\to 0$ [Eq. (4.48)]. The orange dotted line is the
analogous approximation as $u\to\beta\hbar$. (b) The natural logarithm of
$c_{2}/\widetilde{Z}_{\mathrm{r},0}$ (black line) and the corresponding $u\to
0$ and $u\to\beta\hbar$ asymptotic approximations (blue and orange dotted
lines respectively). Note the minimum at $\tau_{\mathrm{g}}/\beta\hbar\approx
0.435$.

Hence to obtain an asymptotic approximation to Eq. (4.47) we require the
asymptotic form of $c_{2}(u)$ as $u\to 0$,

$\displaystyle c_{2}(u)$ $\displaystyle\sim
Z_{\mathrm{r},0}^{\mathrm{sc}}\Delta_{\mathrm{r}}^{2}\,\mathrm{e}^{-\frac{1}{\hbar}(\alpha_{\mathrm{r}}u-\frac{1}{2}\gamma_{\mathrm{r}}u^{2})}$
(4.48) $\displaystyle\sim
Z_{\mathrm{r},0}^{\mathrm{sc}}\Delta_{\mathrm{r}}^{2}\,\mathrm{e}^{-\alpha_{\mathrm{r}}u/\hbar}\left(1+\frac{\gamma_{\mathrm{r}}u^{2}}{2\hbar}\right),$

where $\Delta_{\mathrm{r}}\equiv\Delta(\mathsf{x}_{\mathrm{r}})$, and
$\mathsf{x}_{\mathrm{r}}$ is the minimum of $V_{0}$. The coefficients
$\alpha_{\mathrm{r}}$ and $\gamma_{\mathrm{r}}$ are derived analogously to the
$\alpha_{n}$ and $\gamma_{n}$ in Sec. IV.2,

$\displaystyle\alpha_{\mathrm{r}}$
$\displaystyle=V_{1}\big{(}\mathsf{x}_{\mathrm{r}}\big{)}-V_{0}\big{(}\mathsf{x}_{\mathrm{r}}\big{)}$
(4.49a) $\displaystyle\gamma_{\mathrm{r}}$
$\displaystyle=\sum_{j=1}^{f}\frac{g_{\mathrm{p},j}^{2}\coth\\!\Big{(}\frac{\beta\hbar\omega_{\mathrm{r},\mkern-1.0muj}}{2}\Big{)}}{2m\omega_{\mathrm{r},j}},$
(4.49b)

and $Z_{\mathrm{r},0}^{\mathrm{sc}}$ is defined in Eq. (2.29). As before,
$\omega_{\mathrm{r},j}$ is the frequency of the $j$-th reactant normal mode,
and $g_{\mathrm{p},j}$ is the derivative of the product diabat at
$\mathsf{x}_{\mathrm{r}}$ with respect to that mode. An analogous expansion
can be obtained for $c_{2}(\beta\hbar-u)$ as $u\to 0$ by exchanging the
subscripts $\mathrm{r}\leftrightarrow\mathrm{p}$ and potentials
$V_{0}\leftrightarrow V_{1}$ in Eqs. (4.48) and (4.49). These expressions are
plotted as dotted lines in Fig. 6, alongside the exact $c_{2}(u)$, and are
seen to approximate the exact TCF accurately in regions that contribute
significantly to the integrals over $u$.

Following standard procedure, we substitute the asymptotic expansion of
$c_{2}(u)$ on the second line of Eq. (4.48) into the original integrals, which
yields

$\displaystyle\begin{aligned}
&\frac{1}{\hbar^{2}\widetilde{Z}_{\mathrm{r},0}}\int_{0}^{\frac{\beta\hbar}{2}}\\!\\![\beta\hbar-2u]\,c_{2}(u)\,\mathrm{d}u\sim\frac{1}{\alpha_{\mathrm{r}}}\left[\beta-\frac{2\theta_{\mathrm{r}}}{\alpha_{\mathrm{r}}}\right]\\\
&\qquad{}+\frac{\hbar\gamma_{\mathrm{r}}}{\alpha_{\mathrm{r}}^{2}}\left[-\frac{6\theta_{\mathrm{r}}}{\alpha_{\mathrm{r}}^{2}}+\frac{\beta(3-2\theta_{\mathrm{r}})}{\alpha_{\mathrm{r}}}+\frac{\beta^{2}\mathrm{e}^{-\beta\alpha_{\mathrm{r}}/2}}{4}\right]\end{aligned}$
(4.50a) $\displaystyle\begin{aligned}
&\frac{1}{\hbar^{2}\widetilde{Z}_{\mathrm{s},0}}\int_{0}^{\frac{\beta\hbar}{2}}\\!\\!u\,c_{2}(u_{\mathrm{s}})\,\mathrm{d}u\sim\frac{1}{\alpha_{\mathrm{s}}}\left[\frac{\theta_{\mathrm{s}}}{\alpha_{\mathrm{s}}}-\frac{\beta\mathrm{e}^{-\beta\alpha_{\mathrm{s}}/2}}{2}\right]\\\
&\quad{}+\frac{\hbar\gamma_{\mathrm{s}}}{\alpha_{\mathrm{s}}}\left[\frac{3\theta_{\mathrm{s}}}{\alpha_{\mathrm{s}}^{3}}-\left(\frac{3\beta}{2\alpha_{\mathrm{s}}^{2}}+\frac{3\beta^{2}}{8\alpha_{\mathrm{s}}}+\frac{\beta^{3}}{16}\right)\mathrm{e}^{-\beta\alpha_{\mathrm{s}}/2}\right]\\!,\\!\\!\end{aligned}$
(4.50b)

where $\theta_{\mathrm{s}}=1-\mathrm{e}^{-\beta\alpha_{\mathrm{s}}/2}$, the
subscript ‘s’ is either ‘r’ or ‘p’, $u_{\mathrm{r}}\equiv u$ and
$u_{\mathrm{p}}\equiv\beta\hbar-u$. Substituting Eq. (4.50) into Eq. (4.47)
then gives the semiclassical approximation to $\widetilde{Z}_{\mathrm{r},2}$
that applies when both $V_{0}$ and $V_{1}$ have stable minima. When $V_{1}$ is
unbound, only the first integral in Eq. (4.47) needs to be calculated, as in
this case $Z_{\mathrm{r},0}/Z_{\mathrm{p},0}$ evaluates to zero, and so the
second line of Eq. (4.47) makes no contribution.

The correction discussed here predominantly accounts for changes in the
potential energy at the bottom of the reactant well. In the adiabatic
representation, an increase in $\hat{\Delta}$ lowers the energy of the well,
thus increasing the activation energy barrier. The term in Eq. (4.47) serves
therefore to decrease the overall reaction rate [see Eq. (3.3)]—an effect that
we correctly predict to be negligible if $\hat{\Delta}$ decays to zero in the
vicinity of the PES minima.

Figure 7: The three components of the $k_{4}$ thermal rate, calculated for the
potentials in Eq. (5.1) at constant diabatic coupling $\Delta$. The black
circles are the quantum results. The semiclassical approximations to
$\tilde{k}_{{4\mathrm{A}}}$ on the leftmost panel are calculated according to
Eq. (4.28), with $d_{01}$ given by Eq. (4.27) (red line) or $d_{01}=0$ (blue
line). On the central panel, the semiclassical $\tilde{k}_{{4\mathrm{B}}}$ are
calculated according to Eq. (4.45), in combination with either Eq. (4.44) (red
line) or Eq. (B.1) (blue line). On the rightmost panel, the semiclassical
approximation to $\widetilde{Z}_{\mathrm{r},2}/Z_{\mathrm{r},0}$, in red, is
calculated according to Sec. IV.3. Each panel includes an inset showing the
relative error in the red semiclassical estimate.

##  V Application to model systems

Equation (4.28) for the $\tilde{k}_{{4\mathrm{A}}}$ component of the rate
constant is exact for one-dimensional linear systems. Equation (4.45) for
$\tilde{k}_{{4\mathrm{B}}}$ deviates from the exact result by at most $2.8\%$
for the same type of system. It remains to see how accurate these expressions
are for nonlinear anharmonic potentials and in multiple dimensions. To ensure
high numerical accuracy of our tests, in this section we consider systems for
which the stationary action can either be found using a one-dimensional root-
finding algorithm (Sec. V.1) or is known analytically (Sec. V.2). Since the
semiclassical expressions being tested are based on the golden-rule instanton,
we anticipate no particular difficulties in extending the approach to more
general potentials. We have discussed in other work how to calculate the GR
instanton in such cases using the ring-polymer formalism,[21, 8, 24, 28] and
our fourth-order components can be obtained from the output of these
calculations during post-processing.

###  V.1 Predissociation model

In this example we show that the new theory performs well in anharmonic
potentials, using as our test case the predissociation model considered in
Refs. 92, 93, 25,

$\displaystyle V_{0}(x)$ $\displaystyle=\tfrac{1}{2}m\omega^{2}x^{2}$ (5.1a)
$\displaystyle V_{1}(x)$
$\displaystyle=D_{\mathrm{e}}\mathrm{e}^{-2\alpha(x-\zeta)}-\epsilon,$ (5.1b)

where $m=1$, $\hbar=1$, $\omega=1$, $D_{\mathrm{e}}=2$, $\alpha=0.2$,
$\zeta=5$, $\epsilon=2$. The exact results were calculated by expanding Eqs.
(3.6), (3.11), and the first line of (4.47), in terms of the eigenstates[92]
of $\hat{H}_{0,1}$ and performing the time integrals analytically. Where
applicable, the infinite integration limits for $t$ and $t_{1}$ were replaced
with $t_{\mathrm{max}}=10$, and with $\pi/\omega$ for $t_{0}$. Energies up to
$E_{\mathrm{max}}=75$ were included in the calculation and only those
eigenfunctions of $\hat{H}_{1}$ were considered that had a node at
$x_{\mathrm{max}}=35$, which is equivalent to truncating the space at that
point. $E_{\mathrm{max}}$, $t_{\mathrm{max}}$ and $x_{\mathrm{max}}$ are all
convergence parameters that were confirmed to be sufficiently large for the
range of temperatures considered in this section.

All semiclassical results were derived from the stationary action, which is
directly available[76] for the harmonic reactant potential in Eq. (5.1a), and
can be related to the abbreviated action for the repulsive Morse potential in
Eq. (5.1b), as discussed in Appendix C.

We have taken the diabatic coupling to be constant, $\Delta(x)=\Delta$, and
calculated the quantum and semiclassical GR thermal rates along with their
fourth-order corrections. The three components of the latter are plotted in
Fig. 7 across a set of temperatures ranging from the classical to the deep-
tunnelling regime. For the $\tilde{k}_{{4\mathrm{A}}}$ contribution plotted on
the leftmost panel, we show semiclassical results calculated according to Eq.
(4.28) using two different approximations for $d_{01}$: the red curve uses the
expression in Eq. (4.27), whereas the blue curve sets $d_{01}=0$. The latter
was included to show why it is important to account for the
$t_{0}^{2}t_{1}^{2}$ term in Eq. (4.4b), as failing to do so causes the
semiclassical estimate to diverge from the correct result in the high-
temperature limit. In this example neglecting the quartic term leads to an
overestimation of the $\tilde{k}_{{4\mathrm{A}}}$ component by a factor of 3
at the highest temperature considered in our quantum calculations
($\beta=0.5$).

Figure 8: (a) The golden-rule thermal rate constant $k_{2}$ and (b) the
estimate of its relative deviation from the full rate, $(k-k_{2})/k\sim
k_{4}/k_{2}$, plotted for the systems in Eq. (5.1) at constant diabatic
coupling $\Delta$. The black circles represent the quantum results, and the
red lines are the corresponding semiclassical approximations. Panel (a)
includes an inset showing the error in the semiclassical GR rate relative to
the quantum result, $k_{2}^{\mathrm{sc}}/k_{2}-1$, and panel (b) includes an
inset with the absolute error,
$\Delta^{\\!-2}\left(k_{4}^{\mathrm{sc}}/k_{2}^{\mathrm{sc}}-k_{4}/k_{2}\right)$.

Including the $d_{01}$ term results in a value that is within a few per cent
from the quantum expression down to $\beta=12$. At lower temperatures, the
magnitude of the relative error increases to approximately $20\%$, which is
comparable to the $\approx 25\%$ error developed by the GR instanton rate in
the zero-$T$ limit [Fig. 8(a)]. Further down, the relative error appears to
decrease—an artefact that can be traced down to our quantum results. Their
calculation involves truncating the integration bounds for $t_{0}$ at
$\pm\pi/\omega$ because of recurrences in the $c_{{4\mathrm{A}}}$ flux-
correlation function. At sufficiently low temperatures, this function no
longer decays to zero by the cut-off points, and so the quantum
$\tilde{k}_{{4\mathrm{A}}}$ becomes ill-defined.

The $\tilde{k}_{{4\mathrm{B}}}$ contribution, plotted on the middle panel of
Fig. 7 is also approximated using two different semiclassical expressions: the
red curve follows the general prescription in Eqs. (4.44) and (4.45), and the
blue curve replaces Eq. (4.44) with Eq. (B.1), which would be exact if the
system were linear. The quantum results are calculated according to Eq.
(3.11), with the integration contours for $z_{n}$ deformed as in Eq. (3.13).
We use the $\tau_{\mathrm{g}}$ from the GR instanton to define the deformed
contours and to set the real part of $z$ in Eq. (3.12). Ideally, instead of
using a fixed $\mathord{\operatorname{Re}}\mkern 1.0mu(z)$ for all $z_{n}$,
one would choose an optimal value for each combination $(z_{0},z_{1})$. This
would minimise the oscillations of the integrand in Eq. (3.12) and reduce the
error due to truncating the integration range to
$(-t_{\mathrm{max}},t_{\mathrm{max}})$. However we found this computationally
demanding and used a fixed $\operatorname{Re}(z)=\tau_{\mathrm{g}}$ instead.
We can verify that this is a reasonable approximation by comparing the output
of Eq. (3.11) to the value obtained from only the $\mathcal{A}_{\tau\tau}$
term. The two are expected to be close across the entire temperature range,
and indeed for $\beta\geq 2$ the values are within $10\%$ of each other. For
smaller $\beta$ the magnitudes of $\tilde{k}_{{4\mathrm{A}}}$ and
$\mathcal{A}_{tt}$ are large enough that calculating Eq. (3.11) accurately
poses numerical difficulties. For this reason the $\tilde{k}_{{4\mathrm{B}}}$
plotted at the three highest temperatures includes only the contribution from
$\mathcal{A}_{\tau\tau}$. This change is evident from the small discontinuity
in the relative error shown on the inset.

As anticipated, the two semiclassical approximations are very close at high
temperatures (small $\beta$), but only Eq. (4.45) is well-behaved for large
$\beta$. Eq. (B.1) predicts a rate constant that rapidly diverges from the
quantum result as $\beta\to\infty$. In the same regime, the “good”
approximation develops a constant relative error of around $25\%$. This is not
an error in the fourth-order expression _per se_ , but is rather a consequence
of $\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}$ being proportional to
$\tilde{k}_{2}^{\mathrm{sc}}$, which itself overestimates the GR rate in the
$\beta\to\infty$ limit, as shown in Fig. 8(a). The same observation applies to
the partition function component plotted on the rightmost panel of Fig. 7.
Like $\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}$, it is proportional to the
semiclassical golden-rule rate constant and hence develops a (nearly) constant
relative error in the low-temperature regime, inherited from
$k_{2}^{\mathrm{sc}}$. Otherwise, the quantum and semiclassical results for
this component are in good agreement.

One of the key applications of a fourth-order rate theory is estimating the
error introduced by approximating the full thermal rate with just the second-
order (golden-rule) term. Specifically, we are interested in the relative
error, whose estimate in the small-$\Delta$ limit is $k_{4}/k_{2}$, provided
$k_{4}\neq 0$. The latter is important to bear in mind, since $k_{4}$ contains
both positive and negative components. Their relative magnitudes are
temperature-dependent, such that $-\tilde{k}_{{4\mathrm{A}}}/\tilde{k}_{2}<0$
dominates as $\beta\to 0$ and
$(\tilde{k}_{{4\mathrm{B}}}/\tilde{k}_{2}-\widetilde{Z}_{\mathrm{r},2}/Z_{\mathrm{r},0})>0$
dominates as $\beta\to\infty$. Hence there will be a temperature at which
these expressions cancel exactly. At that point, strictly speaking, one has to
consider the next term in the series, $k_{6}$, in order to estimate the error.

With this in mind, we plot the golden-rule rate constant $k_{2}$ and the ratio
$k_{4}/k_{2}$ in panels (a) and (b) of Fig. 8 respectively, alongside the
corresponding semiclassical approximations. As mentioned previously,
$k_{2}^{\mathrm{sc}}$ develops a constant relative error in the low-
temperature regime, where the rate plateaus. However, because the terms that
dominate $k_{4}^{\mathrm{sc}}$ at such temperatures are themselves
proportional to $k_{2}^{\mathrm{sc}}$, the ratio
$k_{4}^{\mathrm{sc}}/k_{2}^{\mathrm{sc}}$ benefits from near-complete error
cancellation and is a remarkably accurate approximation to the quantum result
across the entire temperature range. Like the quantum result, it
(approximately) plateaus at low temperatures, when only the reactant
vibrational ground state contributes to the process. The semiclassical
expression deviates from quantum by at most $0.06\Delta^{2}$, and the largest
deviations are found for $\beta\leq 2.5$. The errors at these high
temperatures are likely overestimated because of numerical artefacts in the
corresponding quantum calculations, but even the apparent level of accuracy is
more than sufficient to gauge the validity of the GR approximation. Under
favourable circumstances these results should even enable us to correct the GR
approximation, as explored in Sec. V.2.

###  V.2 Spin–boson model

Table 1:  Thermal rate constants for the underdamped spin–boson system in Eq. (5.5), calculated for $\beta=1$, $\hbar=1$, $\Lambda=60$, and $\gamma=\Omega$. The quantum results were calculated by numerical quadrature of the corresponding flux-correlation functions. All results are divided by $\Delta^{2}$.  $\epsilon$ | $0$ | $15$
---|---|---
$\Omega$ | $0.5$ | $2$ | $4$ | $6$ | $0.5$ | $2$ | $4$ | $6$
$k_{2}$ | $7.56\text{\times}{10}^{-08}$ | $1.80\text{\times}{10}^{-07}$ | $1.00\text{\times}{10}^{-06}$ | $5.25\text{\times}{10}^{-06}$ | $5.30\text{\times}{10}^{-05}$ | $1.12\text{\times}{10}^{-04}$ | $4.59\text{\times}{10}^{-04}$ | $1.57\text{\times}{10}^{-03}$
$k_{2}^{\mathrm{sc}}$ | $7.56\text{\times}{10}^{-08}$ | $1.80\text{\times}{10}^{-07}$ | $9.91\text{\times}{10}^{-07}$ | $5.12\text{\times}{10}^{-06}$ | $5.30\text{\times}{10}^{-05}$ | $1.12\text{\times}{10}^{-04}$ | $4.56\text{\times}{10}^{-04}$ | $1.56\text{\times}{10}^{-03}$
$\frac{\tilde{k}_{{4\mathrm{A}}}}{\tilde{k}_{2}}$ | $0.947$ | $0.0744$ | $0.0122$ | $0.0042$ | $0.992$ | $0.0866$ | $0.0181$ | $0.0088$
$\frac{\tilde{k}_{{4\mathrm{A}}}^{\mathrm{sc}}}{\tilde{k}_{2}^{\mathrm{sc}}}$ | $0.859$ | $0.0679$ | $0.0115$ | $0.0042$ | $0.891$ | $0.0778$ | $0.0164$ | $0.0079$
$\frac{\tilde{k}_{{4\mathrm{B}}}}{\tilde{k}_{2}}$ | $0.122$ | $0.0942$ | $0.0637$ | $0.0485$ | $0.114$ | $0.0900$ | $0.0626$ | $0.0487$
$\frac{\tilde{k}_{{4\mathrm{B}}}^{\mathrm{sc}}}{\tilde{k}_{2}^{\mathrm{sc}}}$ | $0.122$ | $0.0956$ | $0.0649$ | $0.0492$ | $0.115$ | $0.0914$ | $0.0639$ | $0.0494$
$\frac{\widetilde{Z}_{\mathrm{r},2}}{Z_{\mathrm{r},0}}$ | $0.0161$ | $0.0162$ | $0.0164$ | $0.0168$ | $0.0225$ | $0.0229$ | $0.0239$ | $0.0255$
$\frac{\widetilde{Z}_{\mathrm{r},2}^{\mathrm{sc}}}{Z_{\mathrm{r},0}^{\mathrm{sc}}}$ | $0.0160$ | $0.0161$ | $0.0164$ | $0.0167$ | $0.0222$ | $0.0225$ | $0.0231$ | $0.0239$
$\frac{k_{4}}{k_{2}}$ | $-0.8417$ | $0.0036$ | $0.0350$ | $0.0274$ | $-0.9000$ | $-0.0194$ | $0.0206$ | $0.0143$
$\frac{k_{4}^{\mathrm{sc}}}{k_{2}^{\mathrm{sc}}}$ | $-0.7527$ | $0.0117$ | $0.0370$ | $0.0284$ | $-0.7986$ | $-0.0088$ | $0.0244$ | $0.0177$

Here our attention turns to the multidimensional spin–boson model. As
mentioned in the introduction, this type of model potential is not our main
target. Even so, the corresponding stationary action is known exactly and can
be evaluated with relative ease (see Appendix D). Furthermore, because the
model has been extensively studied in the past, benchmark quantum results are
available for a broad parameter range.[55]

The diabatic PESs for the spin–boson model in “reaction coordinate” form are
given by

$\displaystyle V_{n}(q,\mathsf{Q})$
$\displaystyle=\frac{1}{2}\Omega^{2}\\!\left(q\pm\sqrt{\frac{\Lambda}{2\Omega^{2}}}\right)^{\mathrlap{\\!\\!2}}+V_{\mathrm{sb}}(q,\mathsf{Q})-\epsilon_{n}$
(5.2a) $\displaystyle V_{\mathrm{sb}}(q,\mathsf{Q})$
$\displaystyle=\sum_{j=1}^{f-1}\frac{1}{2}\omega_{j}^{2}\\!\left(Q_{j}-\frac{c_{j}q}{\omega_{j}^{2}}\right)^{\\!\\!2},$
(5.2b)

where $\mathsf{Q}$ are the bath modes, $q$ is the reaction coordinate and
$\Lambda$ is the Marcus reorganisation energy. The plus and minus signs are
taken for $n=0$ and $1$ respectively, and
$\epsilon_{1}-\epsilon_{0}=\epsilon$. The coupling coefficients $c_{j}$ and
normal-mode frequencies $\omega_{j}$ are all encoded in the spectral density

$J(\omega)=\frac{\pi}{2}\sum_{j=1}^{f}\frac{c_{j}^{2}}{\omega_{j}}\delta(\omega-\omega_{j}).$
(5.3)

Following Ref. 55, we go to the continuum limit ($f\to\infty$) and consider a
purely Ohmic spectral density

$J(\omega)=\gamma\omega,$ (5.4)

where $\gamma$ is the friction coefficient along the reaction coordinate. We
then re-express the diabatic potentials in the conventional spin–boson form,

$V_{n}(\mathsf{x})=\sum_{j=1}^{f}\frac{1}{2}\tilde{\omega}_{j}^{2}\\!\left(x_{j}\pm\frac{\tilde{c}_{j}}{\tilde{\omega}_{j}^{2}}\right)^{\\!\\!2}-\epsilon_{n},$
(5.5)

where the new coordinates $\mathsf{x}$ are related to $(q,\mathsf{Q})$ by an
orthogonal transformation, and $\tilde{\omega}_{j},\tilde{c}_{j}$ derive from
the Brownian oscillator spectral density,[56, 86]

$J(\omega)=\frac{\Lambda}{2}\frac{\gamma\Omega^{2}\omega}{(\omega^{2}-\Omega^{2})^{2}+\gamma^{2}\omega^{2}}.$
(5.6)

The two-state system in Eq. (5.5) has been studied extensively, and some of
the previously derived analytical results are used by us below.

For harmonic potentials, the semiclassical propagator is exact, the
corresponding action $S_{n}$ quadratic in the end-points, and the prefactor
$C_{n}$ independent of positions. Therefore, semiclassical three-time flux-
correlation functions for the spin–boson model (obtained following the
prescriptions in this paper) are identical to their quantum counterparts. As
shown, for example, by Weiss,[86] these are generally of the form
$\mathrm{e}^{-\phi_{4}(z,z_{0},z_{1})/\hbar}$. It follows that the time-
derivatives of $\phi_{4}$ agree exactly with the coefficients $d_{n}$ in Eq.
(4.25), and $\alpha_{n},\gamma_{n}$ in Eq. (4.39). The same considerations
apply to the single-time flux-correlation function from which the GR rate
constant and the partition function correction are derived. In Appendix D we
give the general form of the exact $\phi_{2\nu}(\bm{z})$, where
$\nu=1,\,2,\,\ldots,$ and $\bm{z}$ is a set of $2\nu-1$ complex time
variables. These are used to calculate the quantum $k_{2}$ and $k_{4}$ rate
constants, and it is also from these that we derive all of the coefficients
required for the steepest-descent approximation. The only exception is
$d_{01}$, whose rigorous calculation would in general require high-order
derivatives of the action that are not available from standard GR instanton
calculations. For consistency, we continue using the approximate expression
for $d_{01}$ in Eq. (4.27), with $\omega=\Omega$, even though in this case the
exact value can in principle be obtained.

In Table 1 we list the second- and fourth-order rate constants for a set of
underdamped symmetric ($\epsilon=0$) and asymmetric ($\epsilon=15$) spin–boson
models, with $\beta=1$, $\hbar=1$, $\Lambda=60$, and $\gamma=\Omega$ ranging
from $0.5$ to $6$. With temperature, reorganisation energy and bias kept
constant within each set of systems, rates predicted by classical theories
such as Marcus[27, 29, 30] or Zusman[40, 56, 44] are all of the same order of
magnitude within a given set. Hence any major variation seen in practice is
due to contributions from quantum tunnelling, which becomes more important
with increasing $\Omega$.

All exact results are calculated by numerical integration of the corresponding
time-correlation functions, and $k_{{4\mathrm{B}}}$ only includes the
contribution from the imaginary-time integral over $\mathcal{A}_{\tau\tau}$.
As expected, $k_{2}^{\mathrm{sc}}$ is in close agreement with the exact GR
rate for the entire parameter range. Semiclassical approximations of the
fourth-order rate components are also accurate across the board, lying within
about 10% of quantum.

Notwithstanding the accuracy of its individual components, at $\Omega=2$ the
semiclassical approximation to the _total_ $k_{4}$ rate constant [Eq. (3.3)]
deviates quite significantly from the quantum result. This should come as no
surprise, since the three components contribute to the total fourth-order rate
constant with different signs. One expects therefore to encounter a
combination of parameters for which these contributions cancel exactly. Unless
there is some fortuitous cancellation of errors, semiclassical predictions
around that point are necessarily expected to show a large _relative_ error.
However the key point is that the _absolute_ error remains small. What the
semiclassical theory predicts correctly is that here the fourth-order
dependence of the total rate on $\Delta$ is weak, and that the GR expression
remains accurate up to greater diabatic coupling strengths than would normally
be expected. For additional information one could consider the next term in
the perturbation series, $k_{6}$, which we do not pursue in this work. With
this proviso, the semiclassical results in Table 1 are sufficiently accurate
for establishing whether a system is in the GR regime and estimating the sign
and magnitude of the error made by the GR approximation to the full
nonadiabatic rate.