# An Introduction to Langer’s Theory

Lu Hong
Department of Bioengineering & Therapeutic Sciences
University of California, San Fransicso
San Francisco, CA 94158
<EMAIL_ADDRESS>
###### Abstract

This note provides a pedagogical introduction to Langer’s theory for activated
rate processes in multiple dimensions at the high friction limit, with an
emphasis on the connection between the theory and the property of the backward
committor/splitting probability near the saddle point. The intended audience
is assumed to have some familiarity with linear algebra and statistical
mechanics while knowledge of stochastic processes is not strictly necessary.

Keywords: Kramers’ theory, KLBS theory

###### Contents

  1. 1 A brief review of stochastic process
    1. 1.1 Overdamped Langevin equation
    2. 1.2 Smoluchowski equation
  2. 2 Kramers’ theory at the high friction limit
    1. 2.1 Kramers’ theory via the flux-over-population method
    2. 2.2 Kramers’ theory via MFPT
  3. 3 Langer’s extension to higher dimensions
    1. 3.1 A relation between the steady-state probability and the backward committor
    2. 3.2 An ansatz for the backward committor
    3. 3.3 Properties of the vector normal to the isocommittor planes
      1. 3.3.1 v+ as a generalized eigenvector
      2. 3.3.2 The value of the free constant a
      3. 3.3.3 The magnitude of v+
    4. 3.4 Properties of the steady-state flux vector field
      1. 3.4.1 The flux vector field has a constant direction
      2. 3.4.2 The flux vector field is constant along its streamlines
      3. 3.4.3 Flux surface integrals are constant over dividing surfaces
    5. 3.5 Langer’s rate constant
    6. 3.6 Limitation of Langer’s theory
  4. 4 Simplification of Langer’s theory by Berezhkovskii and Szabo
    1. 4.1 Gaussian surface integral over a plane
    2. 4.2 Projection of Langer’s result to one dimension
  5. 5 Acknowledgement

This note is intended as an introduction to Langer’s theory [1], the
multidimensional extension of Kramers’ theory for activated rate processes in
the high friction limit [2]. The note is organized as follows: Section 1
briefly reviews some elements of the theory of stochastic process, which
prepares for the introduction of Kramers’ theory in Section 2, where the flux-
over-population method is demonstrated in the relatively simple context of a
one-dimensional double-well potential. The majority of this note is devoted to
the discussion of Langer’s theory in Section 3. The derivation begins by
relating the steady-state probability density in the flux-over-population
method to the backward committor, followed by a detailed analysis of the
behavior of the committor near the saddle point. This allows us to derive an
expression for the probability flux vector field, culminating in the
presentation of the multidimensional rate constant in equation (106). Lastly,
Section 4 discusses some more recent development by Berezhkovskii and Szabo
[3] that allows one to project Langer’s result back to one dimension.
Together, the results presented in this note are sometimes referred to as the
KLBS theory.

Given the venerable age of Langer’s theory and the existence of several review
articles [4, 5, 6] and textbook [7] partly devoted to this subject, one might
wonder why this note is necessary. This introduction certainly draws heavily
on the aforementioned works, but differs in two regards. First, almost all
current accounts of Langer’s theory are quite concise. While this feature is
excellent for experts who are in need of a quick refresher, it poses
significant challenges for beginners, because the derivation of the theory can
be quite tedious and involve some mathematical tricks that are not standard
knowledge for a typical reader. One goal of this note, therefore, is to derive
Langer’s theory while erring on the side of presenting too many, rather than
too few, algebraic details.

Second, an interesting feature of Langer’s theory is that the expression for
the rate constant critically depends on several basic geometric features of
the committor/splitting probability near the saddle point. This connection was
not explicitly stated in Langer’s original paper, and is usually not fully
developed in recent accounts of the theory. This is, again, unfortunate, since
this connection provides a simple setting for gaining an intuitive
understanding of the behavior of the committor near an idealized transition
state. Given the prominent role the committor plays in modern rate theories,
such as the transition path theory [8], such intuitions can be valuable in
understanding the more recent development. Therefore, as described above, a
second goal of this note is to make explicit the connection of Langer’s theory
to the committor function and provide a detailed analysis thereof.

The primary intended audience for this note are students and researchers in
chemistry and biophysics interested in condensed-phase simulations of
macromolecules. The reader should be familiar with linear algebra and have
some exposure to classical statistical mechanics and some rate theory (e.g.,
Arrhenius equation); knowledge of stochastic processes is not strictly
necessary although results from the theory will be invoked in some
derivations.

##  1 A brief review of stochastic process

In this section we briefly review some concepts from the theory of stochastic
process relevant for our discussion. Readers who have never studied this
subject should either consult standard references on stochastic processes and
statistical mechanics, or simply take the results presented below as given and
fill in their missing background knowledge later. The presentation in this
section partially follows that in Chapter 15 of [9].

###  1.1 Overdamped Langevin equation

Consider a (closed) system that is in contact with a large heat bath. In a
typical (classical) simulation, the system consists of the macromolecule(s) of
interest as well as some water molecules and counterions. The time evolution
of the full system is determined by the Hamiltonian equations, and the
equilibrium phase space distribution function is microcanonical. To avoid
explicit representation of the bath, we “abstract away” the bath degrees of
freedom, and the remaining system dynamics can be described via the
generalized Langevin equation:

$\displaystyle r_{i}^{\prime}(t)$ $\displaystyle=\frac{p_{i}(t)}{m_{i}}$ (1)
$\displaystyle p_{i}^{\prime}(t)$ $\displaystyle=-\frac{\partial
W(r)}{\partial
r_{i}}-\int_{0}^{t}\sum_{j}\sqrt{m_{i}m_{j}}\gamma_{ij}(t-\tau)r_{j}^{\prime}(\tau)\,d\tau+\xi_{i}(t)$
(2)

Here, $r_{i}$, $p_{i}$, and $m_{i}$ are the position, momentum, and
(renormalized) mass of the $i$th degree of freedom in the remaining system,
and $W(r)$ is the potential of mean force obtained from the full potential
energy function by averaging over the bath degrees of freedom. The generalized
Langevin equation can be derived via either the harmonic bath model or, more
rigorously, via the Mori–Zwanzig theory.

The generalized Langevin equations is reminiscent of Newton’s equation of
motion, but the price we pay for doing away with the bath degrees of freedom
is the difficulty of dealing with a set of stochastic integro-differential
equations. Specifically, compared to Newton’s equations, there are two
additional terms in (2) that provide a coarse-grained model of the effect of
the bath on the system. Here, we have “decomposed” the effect of the bath into
two parts. First, the fluctuation term $\xi_{i}(t)$ represents a random force
(i.e., noise) acting on the system. Although the motion of the bath is fully
deterministic, by ignoring the molecular details we model $\xi_{i}(t)$ as a
random, or stochastic, process. For a system solvated by a dense bath, such as
liquid water, that affects the system dynamics through soft collisions (i.e.,
weak noise), a common model for $\xi_{i}(t)$ is a Gaussian random process with
zero mean.

Second, the dissipation term
$\int_{0}^{t}\gamma_{ij}(t-\tau)r_{j}^{\prime}(\tau)\,d\tau$ is a convolution
integral that acts as a friction force slowing down the system. This integral
is called a memory integral; the term $\gamma_{ij}(t)$ is known as the memory
kernel, or dynamic friction kernel, and it encodes the “memory” of the motion
of the system by the bath. Physically, the memory integral represents the fact
that the bath requires a finite amount of time to respond to fluctuations in
the system and this lag affects the motion of the system. At equilibrium, the
fluctuation term and the dissipation term are related by the second
fluctuation-dissipation theorem,

$\mathbb{E}[\xi_{i}(t)\xi_{j}(t^{\prime})]=k_{B}T\sqrt{m_{i}m_{j}}\gamma_{ij}(|t-t^{\prime}|)$
(3)

Now, we make some further simplifying assumptions to make (2) more
analytically tractable. First, we assume that the bath responds
instantaneously to the motion of the system; i.e., we assume that the memory
integral decays instantaneously and the bath has no memory of the system
history. This is a good model when the (renormalized) mass of the system is
much larger than that of the bath. In this case, the memory kernel becomes
$\gamma_{ij}(t)=2\gamma_{ij}\delta(t)$, where we have defined

$\gamma_{ij}=\int_{0}^{\infty}\gamma_{ij}(t)\,dt$ (4)

as the static friction kernel, or simply the friction coefficient, and
$\delta(t)$ is a Dirac delta function. The fluctuation-dissipation theorem now
reads

$\mathbb{E}[\xi_{i}(t)\xi_{j}(t^{\prime})]=2k_{B}T\sqrt{m_{i}m_{j}}\gamma_{ij}\delta(|t-t^{\prime}|)$
(5)

Stochastic processes with an autocorrelation function of this form are called
white noises, and the terminology reflects the fact that the power spectral
density of the process is a constant over all frequencies. With these
assumptions, the generalized Langevin equation becomes

$\displaystyle r_{i}^{\prime}(t)$ $\displaystyle=\frac{p_{i}(t)}{m_{i}}$ (6)
$\displaystyle p_{i}^{\prime}(t)$ $\displaystyle=-\frac{\partial
W(r)}{\partial
r_{i}}-\sum_{j}\sqrt{m_{i}m_{j}}\gamma_{ij}r_{j}^{\prime}(t)+\xi_{i}(t)$ (7)

which is known as the Langevin equation. Compared to the generalized Langevin
equation, the Langevin equation describes stochastic processes that are
Markovian (i.e., memoryless).

Second, for dense solvent, such as water, the high friction and frequent
collisions with the system leads to, on a short timescale (i.e.,
$t<\gamma^{-1}$), rapid fluctuations in the acceleration $p_{i}^{\prime}(t)$.
However, on a longer timescale the change in the time-averaged velocity will
be small as the effect of collisions cancel out each other. This allows us to
set the acceleration in (7) to zero, which leads to

$\sum_{j}\sqrt{m_{i}m_{j}}\gamma_{ij}r_{j}^{\prime}(t)=-\frac{\partial
W(r)}{\partial r_{i}}+\xi_{i}(t)$ (8)

The motion described by (8) goes by many names, such as diffusion, Brownian
motion, or overdamped Langevin dynamics; the equation itself is known as the
overdamped Langevin equation. In some applications, the cross-correlation
terms in the memory kernel are ignored, in which case (8) takes on a simpler
form

$m_{i}\gamma_{ii}r_{i}^{\prime}(t)=-\frac{\partial W(r)}{\partial
r_{i}}+\xi_{i}(t)$ (9)

###  1.2 Smoluchowski equation

For the purpose of describing the kinetics of activated rate processes at
equilibrium, working directly with (8) is inconvenient, since trajectories
consistent with (8) are individual realizations of the system dynamics, while
we are more interested in the statistics of an ensemble of such realizations.
In other words, we are more interested in $p(r,t)$, the probability density of
the system (strictly speaking, $p(r,t)$ is a conditional probability density
function more appropriated denoted as $p(r,t|r_{0},t_{0})$). The time
evolution of $p(r,t)$ for processes governed by (8) is given by

$\frac{\partial p(r,t)}{\partial t}=-\nabla\cdot J(r,t)\quad\text{with}\quad
J(r,t)=-D\pi(r)\nabla\frac{p(r,t)}{\pi(r)}$ (10)

where $\pi(r)$ is the Boltzmann distribution in the configuration space and
the stationary solution to (10), $J(r,t)$ is the (probability) flux, and $D$
is the (position-independent) diffusion matrix, which we assume to be
symmetric positive definite. Elements of the diffusion matrix are related to
the static friction kernel defined in (4) via the relation

$D_{ij}=\frac{k_{B}T}{\sqrt{m_{i}m_{j}}\gamma_{ij}}$ (11)

Equations such as (10) are known as the Smoluchowski equation, a special case
of a class of partial differential equations known as the Fokker-Planck
equations, or Kolmogorov’s forward equations. It is also common to refer to
the first part of (10) as the continuity equation for the probability density,
and the second part of (10) as a “constitutive” equation.

For the following discussion it is also convenient to rewrite (10) in two
forms. First, we can rewrite (10) component-wise:

$\partial_{t}p(r,t)=\sum_{ij}\partial_{i}D_{ij}\left[\beta\partial_{j}W(r)+\partial_{j}\right]p(r,t)$
(12)

where we have used a short-hand notation $\partial_{t}=\partial/\partial t$
and $\partial_{i}=\partial/\partial r_{i}$. For readers not familiar with the
component-wise notation, they should convince themselves that, e.g., the $i$th
component of a matrix-vector product $Ax$ can be expressed as
$\sum_{j}A_{ij}x_{j}$.

Second, we can rewrite (10) as an equation involving the Fokker-Planck
operator $L^{\dagger}$,

$\frac{\partial p(r,t)}{\partial t}=L^{\dagger}p(r,t)\quad\text{where}\quad
L^{\dagger}=\nabla\cdot D\pi(r)\nabla\pi^{-1}(r)$ (13)

Here, the symbol $\dagger$ indicates that the Fokker-Planck operator
$L^{\dagger}$ is the adjoint of another operator $L$ known as the generator,

$L=\pi^{-1}(r)\nabla\cdot D\pi(r)\nabla$ (14)

The concept of an adjoint operator is a generalization of the Hermitian
transpose of a matrix. As (13) shows, the Fokker-Planck operator dictates the
time evolution of the probability density $p(r,t)$. There is a similar
interpretation for the generator: the generator dictates the time evolution of
(conditional) ensemble averages, or observables, of the form
$u(r,t)=\mathbb{E}[f(r(t))|r(0)=r_{0}]$, for some suitable scalar function $f$
defined on the configuration space; i.e.,

$\frac{\partial u(r,t)}{\partial t}=Lu(r,t)$ (15)

##  2 Kramers’ theory at the high friction limit

In this section we consider Kramers’ theory in the context of a one-
dimensional Brownian particle whose motion is described by the one-dimensional
version of the overdamped Langevin equation (8),

$m\gamma
r^{\prime}(t)=-\frac{dW}{dr}+\xi(t)\quad\text{where}\quad\mathbb{E}[\xi(t)]=0,\mathbb{E}[\xi(t)\xi(t^{\prime})]=2k_{B}Tm\gamma\delta(|t-t^{\prime}|)$
(16)

From (10), it follows that the probability density $p(r,t)$ satisfies the
following Fokker-Planck equation

$\frac{\partial p(r,t)}{\partial t}=-\frac{\partial J(r,t)}{\partial
r}\quad\text{with}\quad
J(r,t)=-D\left[\beta\frac{dW}{dr}+\frac{\partial}{\partial r}\right]p(r,t)$
(17)

where $D=k_{B}T/m\gamma$ is the position-independent diffusion constant. The
stationary solution to (17) is the Boltzmann distribution $\pi(r)=e^{-\beta
W(r)}/Z$, where $Z$ is a configurational partition function.

Here, we take $W(r)$ to be an (asymmetric) double-well potential, where the
minimum of the reactant well $A$ is located at $r_{A}$, the minimum of the
product well $B$ is located at $r_{B}$, and the transition state
$r^{\ddagger}$ is identified as the position of the peak of $W(r)$ between the
two wells; without loss of generality, we assume that
$r_{A}<r^{\ddagger}<r_{B}$. An example of such a potential is shown in Figure
1. Further, we assume that the height of the barrier $W(r^{\ddagger})$ is much
larger than $k_{B}T$, so that there is a separation of timescales between
barrier crossing and within-well equilibration. A related assumption here is
that barrier crossing is much slower than the correlation time of the dynamic
friction kernel so that the static friction kernel approximation can be
justified. Together, these assumptions amount to the situation where a single
slow degree of freedom in the system $r$ is sufficient for describing the
reaction, while the rest of the degrees of freedom relax much faster than the
timescale of barrier crossing along $r$; we should thus interpret “bath” as
also including other degrees of freedom of the macromolecule that are not
explicitly treated, and interpret $r$ not necessarily as a position variable
but more generally as a reaction coordinate (i.e., some function of the
position variables $r_{i}$).

We mention here in passing that Kramers also derived expressions for the rate
constant in the weak and moderate-to-high friction regimes. We will not
discuss these results since they are not as relevant for condensed-phase
simulations.

Figure 1: An example one-dimensional double-well potential of mean force. The
potential is given by
$W(r)=-\ln\left(e^{-2(x+1)^{2}}+e^{-2(x-1)^{2}+1}\right)$.

###  2.1 Kramers’ theory via the flux-over-population method

Kramers derived an expression for the rate constant of the $A\to B$ reaction,
$k_{AB}$, using what is now known as the flux-over-population method. Consider
a hypothetical procedure in which an ensemble of $n$ particles are prepared in
the reactant well $A$ (usually we set $n=1$), where they rapidly reach thermal
equilibrium on a timescale much faster than that of escaping to the product
well $B$. Whenever a particle escapes from the reactant well and reaches the
product well, it is immediately removed and a new one is added to the reactant
well, such that the reactant population is always maintained at $n$. As we
show in the analysis below, the exact positions at which a particle is
considered escaped and at which a new one is inserted are not very important.
As the system reaches a non-equilibrium steady-state, there is a non-zero
probability current $J$, or flux, into the product well, which is the reaction
rate (i.e., the number of transitions from $A$ to $B$ per unit time). We
normalize this reaction rate by the reactant population, which gives the rate
constant,

$k_{AB}=J/n$ (18)

hence the name “flux-over-population”.

Before flushing out the algebraic details, let us pause for a moment and
consider why the flux-over-population procedure is needed. In other words, why
does the calculation of an equilibrium rate constant invoke the flux of a
seemingly contrived non-equilibrium process? A simple answer is that the (net)
flux between any two states (micro- or macro-) in equilibrium is zero, due to
detailed balance, and thus gives no information about kinetics. To see why the
flux-over-population procedure circumvent this problem, let us consider the
behavior of particles at some point near the transition state and see how they
contribute to the flux across that point at equilibrium. At any given moment,
such a particle can be categorized into one of four groups, depending on its
past and future behavior:

  1. 1.

The particle came from $A$ and will move to $B$ before going back to $A$.

  2. 2.

The particle came from $A$ and will go back to $A$ before moving to $B$.

  3. 3.

The particle came from $B$ and will move to $A$ before going back to $B$.

  4. 4.

The particle came from $B$ and will go back to $B$ before moving to $A$.

Behavior described in group 2 and 4 are known as barrier recrossing and does
not contribute to the flux, since each pair of crossing-recrossing cancel out
each other (this is to be contrasted with the procedure in transition-state
theory, where all (re)crossing events are counted towards the total flux
across some dividing surface [10]). Because the net equilibrium flux is zero,
this implies that the flux contributed by group 1 and 3 cancels as well. Note
that, for ergodic dynamics, this categorization is exhaustive and it is not
possible for a particle to cross the barrier once and stay in one of the two
wells forever. By removing particles reaching the product well, the flux-over-
population procedure removes the flux contribution from group 3 (and group 4),
while the only remaining nonzero flux contribution from group 1 stays close to
its equilibrium value at the steady state due to the separation of timescales.
As such, the nonequilibrium steady-state flux $J$ in this procedure gives the
true $A$ to $B$ reaction rate.

At steady state, the population distribution $p_{\text{ss}}(r)$ is time-
independent, and thus by (17) the flux $J$ is now both time- and position-
independent. The two quantities are related via

$J=-D\left[\beta\frac{dW}{dr}+\frac{\partial}{\partial
r}\right]p_{\text{ss}}(r)=-De^{-\beta W(r)}\frac{d}{dr}\left[e^{\beta
W(r)}p_{\text{ss}}(r)\right]$ (19)

Now, divide both sides by $De^{-\beta W(r)}$ and integrate both sides over the
interval $[r_{A},r_{B}]$, which gives

$J\int_{r_{A}}^{r_{B}}D^{-1}e^{\beta W(r)}\,dr=-e^{\beta
W(r)}p_{\text{ss}}(r)\big{|}_{r=r_{A}}^{r_{B}}$ (20)

We now consider the behavior of $p_{\text{ss}}(r)$ near the two boundaries:
$r_{A}$ and $r_{B}$. At $r=r_{B}$, we require $p_{\text{ss}}(r_{B})=0$ to
satisfy the absorbing boundary condition caused by particle removal. At
positions in the reactant well away from the absorbing boundary condition at
$r_{B}$, we assume that the system is close to equilibrium, i.e.,
$p_{\text{ss}}(r)\approx\pi(r)$, because the particle insertion procedure
maintains the reactant population $n$ in the reactant well, and inserted
particles thermalize much faster than the reaction timescale. This observation
has two consequences. First, it implies that
$p_{\text{ss}}(r_{A})=\pi(r_{A})$. Second, it implies that the population in
the reactant well is

$n=\int_{-\infty}^{r^{\ddagger}}p_{\text{ss}}(r)\,dr\approx\int_{-\infty}^{r^{\ddagger}}\pi(r)\,dr$
(21)

Now, apply the boundary conditions at $p_{\text{ss}}(r_{A})$ and
$p_{\text{ss}}(r_{B})$ to (20), then apply the flux-over-population formula
(18), and we arrive at

$\displaystyle k_{AB}$ $\displaystyle=\frac{e^{\beta
W(r_{A})}p_{\text{ss}}(r_{A})-e^{\beta
W(r_{B})}p_{\text{ss}}(r_{B})}{\int_{r_{A}}^{r_{B}}D^{-1}e^{\beta
W(r)}\,dr\int_{-\infty}^{r^{\ddagger}}p_{\text{ss}}(r)\,dr}$ (22)
$\displaystyle=\frac{e^{\beta
W(r_{A})}\pi(r_{A})}{\int_{r_{A}}^{r_{B}}D^{-1}e^{\beta
W(r)}\,dr\int_{-\infty}^{r^{\ddagger}}\pi(r)\,dr}$ (23)
$\displaystyle=\frac{1}{\int_{r_{A}}^{r_{B}}D^{-1}e^{\beta
W(r)}\,dr\int_{-\infty}^{r^{\ddagger}}e^{-\beta W(r)}\,dr}$ (24)

If we assume that $W(r)$ is harmonic around $r^{\ddagger}$ and $r_{A}$, we
write, for $W(r)$ near $r^{\ddagger}$ and $r_{A}$,

$W(r)\approx\Delta
W^{\ddagger}-\frac{1}{2}\kappa^{\ddagger}\left(r-r^{\ddagger}\right)^{2}\quad\text{and}\quad
W(r)\approx\frac{1}{2}\kappa_{A}\left(r-r_{A}\right)^{2}$ (25)

where we have defined $\Delta W^{\ddagger}=W(r^{\ddagger})-W(r_{A})$ and set
$W(r_{A})=0$ (resetting the zero of the potential of mean force has the effect
of changing the partition function, which has no effect on our derivation).
Here, $\kappa_{A}$ and $\kappa^{\ddagger}$ are the force constants of the
harmonic potentials; equivalently, $\kappa_{A}$ and $\kappa^{\ddagger}$ can
also be interpreted as the curvatures at $r_{A}$ and $r^{\ddagger}$,
respectively. Applying (25) to (24) gives the following approximate results

$\displaystyle\int_{r_{A}}^{r_{B}}D^{-1}e^{\beta W(r)}\,dr$
$\displaystyle\approx\int_{-\infty}^{\infty}D^{-1}e^{\beta\left(\Delta
W^{\ddagger}-\kappa^{\ddagger}r^{2}/2\right)}\,dr=D^{-1}\sqrt{\frac{2\pi}{\beta\kappa^{\ddagger}}}e^{\beta\Delta
W^{\ddagger}}$ (26) $\displaystyle\int_{-\infty}^{r^{\ddagger}}e^{-\beta
W(r)}\,dr$
$\displaystyle\approx\int_{-\infty}^{\infty}e^{-\beta\kappa_{A}r^{2}/2}\,dr=\sqrt{\frac{2\pi}{\beta\kappa_{A}}}$
(27)

These results provide a good approximation whenever the region of validity for
the harmonic approximation is large enough that the added probability mass by
taking the integration limit to infinity is negligible.

Substituting (26) and (27) back to (24) gives

$k_{AB}=\frac{\beta D}{2\pi}\sqrt{\kappa_{A}\kappa^{\ddagger}}e^{-\beta\Delta
W^{\ddagger}}$ (28)

This is the rate constant predicted by Kramers’ theory at the high friction
limit. Since the dynamics at high friction is diffusive, it is also known as
the spatial-diffusion-limited rate. The rate constant for the reverse reaction
can be obtained analogously by replacing the reactant well $A$ with the
product well $B$ in the preceding derivation. If the diffusion constant is
position-dependent, the derivation up to (and include) (24) is still valid.
The diffusion constant can be “folded” into the potential of mean force $W(r)$
as

$e^{\beta W(r)}/D(r)=e^{\beta W(r)-\ln D(r)}=e^{\beta\left[W(r)-k_{B}T\ln
D(r)\right]}$ (29)

This defines a new potential surface $W(r)-k_{B}T\ln D(r)$, which may reach
its maximum at a position other than $r^{\ddagger}$.

In the literature it is also common to parameterize the harmonic potential
near $r^{\ddagger}$ and $r_{A}$ as

$W(r)\approx\Delta
W^{\ddagger}-\frac{1}{2}m{\omega^{\ddagger}}^{2}\left(r-r^{\ddagger}\right)^{2}\quad\text{and}\quad
W(r)\approx\frac{1}{2}m\omega_{A}^{2}\left(r-r_{A}\right)^{2}$ (30)

Here, $\omega_{A}$ (and analogously, $\omega^{\ddagger}$) is the angular
frequency of a harmonic oscillator of mass $m$ on the potential surface
$W(r)=m\omega_{A}^{2}r^{2}/2$; $\omega_{A}$ is related to the force
constant/curvature $\kappa_{A}$ by $\kappa_{A}=m{\omega_{A}}^{2}$. Together
with the definition of the diffusion constant $D=1/\beta m\gamma$, (28) can be
written alternatively as

$k_{AB}=\frac{\omega_{A}\omega^{\ddagger}}{2\pi\gamma}e^{-\beta\Delta
W^{\ddagger}}$ (31)

###  2.2 Kramers’ theory via MFPT

There is an alternative approach to deriving (28) using the mean first-passage
time (MFPT). Let $\Omega$ be the reaction coordinate space of the system. For
a subset $B\subseteq\Omega$, the MFPT $m_{B}(r)$ is the average time for the
system initiated at $r$ to reach $B$. It can be shown, using potential theory,
that the MFPT satisfies the following (Dirichlet) boundary value problem,

$\begin{cases}Lm_{B}(r)=-1,&r\in\Omega\backslash B\\\ m_{B}(r)=0,&r\in
B\end{cases}$ (32)

where the notation $\Omega\backslash B$ denotes the complement of $B$ in
$\Omega$. For one-dimensional Brownian dynamics, the generator $L$ defined in
(14) takes on the form

$L=\frac{d}{dr}D\left[\frac{d}{dr}-\beta
W(r)\right]=D\frac{d^{2}}{dr^{2}}-D\beta\frac{dW}{dr}\frac{d}{dr}$ (33)

Unfortunately, a full derivation of (32) from scratch would take too long;
interested readers should consult standard references on Markov processes for
more details. Nevertheless, an “intuitive” justification of (32) goes as
follows: consider a hypothetical procedure where an ensemble of system
trajectories are prepared, all with the same initial condition at position
$r(0)=r_{0}$ outside $B$; we can calculate $m_{B}(r_{0})$ by taking the
average first hitting time to the set $B$ among the trajectories. For each
given trajectory $r(t)$, moving forward in time by a small amount $\delta t$
reduces the first hitting time at $r(t+\delta t)$ by the same amount of
$\delta t$; as such, the time derivative of $m_{B}(r_{0})$ should be $-1$.
Since the generator acting on an observable $m_{B}(r_{0})$ gives its time
derivative, it follows that $Lm_{B}(r_{0})=-1$. On the other hand, if $r_{0}$
is already in $B$, then the first hitting time to $B$ is zero by construction,
hence the boundary condition in (32).

To determine the rate constant for the transition from the reactant well to
the product well, let us consider the MFPT to the product well from some point
$r$ near the reactant well (as we will see, the particular choice of $r$ is
not very important). More specifically, let $\Omega=\mathbb{R}$ and consider
the MFPT for the set $B=[r_{B},\infty)$. The differential equation
$Lm_{B}(r)=-1$ can be solved using an integrating factor,

$\displaystyle
D\frac{d^{2}m_{B}}{dr^{2}}-D\beta\frac{dW}{dr}\frac{dm_{B}}{dr}$
$\displaystyle=-1$ (34) $\displaystyle\frac{d}{dr}\left[e^{-\beta
W(r)}\frac{dm_{B}}{dr}\right]$ $\displaystyle=-D^{-1}e^{-\beta W(r)}$ (35)
$\displaystyle\int_{-\infty}^{t}\frac{d}{ds}\left[e^{-\beta
W(s)}\frac{dm_{B}}{ds}\right]\,ds$
$\displaystyle=-D^{-1}\int_{-\infty}^{t}e^{-\beta W(s)}\,ds$ (36)
$\displaystyle\int_{r}^{r_{B}}\frac{dm_{B}}{dt}\,dt$
$\displaystyle=-\int_{r}^{r_{B}}D^{-1}e^{\beta
W(t)}\int_{-\infty}^{t}e^{-\beta W(s)}\,dsdt$ (37) $\displaystyle m_{B}(r)$
$\displaystyle=\int_{r}^{r_{B}}D^{-1}e^{\beta W(t)}\int_{-\infty}^{t}e^{-\beta
W(s)}\,dsdt$ (38)

It is easy to check that, under the harmonic approximations of (25), the
function $e^{\beta W(t)}e^{-\beta W(s)}$ has a unique global maximum at
$(s=r_{A},t=r^{\ddagger})$ within the domain of integration in (38) (see
Figure 2). Since $e^{\beta W(t)}e^{-\beta W(s)}$ decays exponentially away
from the global minimum, we can simply extend the domain of integration to
$\mathbb{R}^{2}$, in which case the double integral in (38) simplifies to the
integrals considered in (26) and (27). Note that the validity of this
approximation implies that $m_{B}(r)$ is a constant near and to the left of
the reactant well, and that the particular choice for the boundaries of $B$ is
not important as long as its lower bound is sufficiently to the right of
$r^{\ddagger}$.

Taken together, we have shown that the rate constant (28) derived from the
flux-over-population method is equivalent to the inverse of the mean first-
passage time; in fact, this equivalence is exact for a much broader range of
stochastic processes than those considered in this note [11]. Unfortunately,
despite its conceptual simplicity, the MFPT cannot in general be calculated
analytically for problems in more than one dimension and thus we will not
consider this method in the following sections; however, see Section VII. D.
of [4] for an approximate approach to solving for the MFPT in higher
dimensions using asymptotic methods.

Figure 2: The region of integration for the double integral in (38). The
region of integration is shown in green. The maxima of $e^{\beta
W(t)}e^{-\beta W(s)}$ are shown as intersections of the orange lines. Only one
of the two maxima, $(s=r_{A},t=r^{\ddagger})$, is located inside the
integration limit.

##  3 Langer’s extension to higher dimensions

Now let us consider the extension of Kramers’ theory to an $N$-dimensional
potential (of mean force) with the reactant and product wells separated by a
saddle point; recall that the saddle point is a point on the potential energy
surface $W(r)$ where the gradient is zero and the Hessian has exactly one
unstable mode. An example two-dimensional double-well potential is shown in
Figure 3A. Conceptually, the flux-over-population method employed in the one-
dimensional case carries over with few changes. However, as we will see, the
higher-dimensional setting of (10) does pose some algebraic challenges. The
derivation presented here loosely follows the presentation in [1, 4, 5].

Similar to our treatment in one dimension, we assume a harmonic approximation
at $r_{A}$, the bottom of the reactant well $A$,

$W(r)\approx\frac{1}{2}(r-r_{A})^{T}H_{A}(r-r_{A})$ (39)

where $H_{A}$ is the Hessian matrix of second-order derivatives evaluated at
the minimum $r_{A}$, which we assume to be symmetric positive definite; note
that the Hessian is also known as the force constant matrix (to see why, note
that for a harmonic potential $\kappa r^{2}/2$ with force constant $\kappa$,
the second-order derivative at the center $r=0$ is simply $\kappa$). In terms
of notations, in this note we will exclusively use notations of the form
$x^{T}y$ to denote inner products in place of other common notations such as
$x\cdot y$ or $\left<x,y\right>$; readers more familiar with the other
notations should convince themselves that, e.g., the quadratic form in (39)
can be expressed in the angle bracket notation as
$\left<r-r_{A},H_{A}(r-r_{A})\right>/2$.

In addition, we assume that there is a saddle point $r^{\ddagger}$ in between
the reactant and product wells, which is the point with the minimum energy
along the barrier ridge. A harmonic approximation at the saddle point gives

$W(r)\approx\Delta
W^{\ddagger}+\frac{1}{2}(r-r^{\ddagger})^{T}H^{\ddagger}(r-r^{\ddagger})$ (40)

Here, we assume that $H^{\ddagger}$ has exactly one negative eigenvalue
associated with the unstable mode at $r^{\ddagger}$, while the rest of the
eigenvalues are assumed to be strictly positive. In general, it is possible
for $H^{\ddagger}$ to possess one or more eigenvalues of zero, which
correspond to some underlying symmetries of the system. For example, the
Hessian matrix of a system of two one-dimensional particles with positions
$x_{1}$ and $x_{2}$ whose distance is subject to the harmonic restraint
$W(x_{1},x_{2})=\kappa(r_{0}-(x_{1}-x_{2}))^{2}/2$ will have a zero eigenvalue
corresponding to the translational symmetry of the system. Such symmetries can
oftentimes be eliminated (e.g., here by working with the reaction coordinate
$r=x_{1}-x_{2}$) and we shall not consider such cases for the following
analysis. Lastly, we record here for future reference,

$\partial_{i}W(r)=\sum_{j}H^{\ddagger}_{ij}\left(r_{j}-r^{\ddagger}_{j}\right)$
(41)

This result comes from the fact that $H^{\ddagger}$ is symmetric, and, for any
quadratic form $x^{T}Ax$, $\nabla_{x}x^{T}Ax=(A+A^{T})x$.

###  3.1 A relation between the steady-state probability and the backward
committor

Figure 3: An example two-dimensional double-well potential. A) Contour plot of
the double-well potential; the potential is given by
$W(x,y)=\frac{5}{2}(x^{2}-1)^{2}+5y^{2}$ (see also [12]). The reactant well
$A$ and the product well $B$ are defined, somewhat arbitrarily, as disks of
radius $1/4$ centered at the global minima $r_{A}$ and $r_{B}$. The white
dashed vectors show the eigenvectors of $H^{\ddagger}$, while the white solid
vectors show the eigenvectors of $H^{\ddagger}D$, with
$D=\\{\\{4,1\\},\\{1,2\\}\\}$. B) The exact backward committor $q^{-}(x,y)$
and steady-state flux $J(x,y)$. The committor is computed by solving (44)
numerically using the NDSolveValue function with the appropriate
DirichletCondition in Mathematica 12.1 and is shown as a contour plot (with
0.1 contour increment). The steady-state flux is calculated using (76) and
represented as the orange vector field; the lengths of the streamline segments
are proportional to the vector field magnitude. C) The approximate backward
committor and steady-state flux from Langer’s theory. The committor is given
by (74), and the steady-state flux is given by (78) with $\pi(s)$ replaced
with its harmonic approximation at the saddle point. For all calculations
$\beta$ is set to one.

Instead of seeking a steady-state solution $p_{\text{ss}}(r)$ to (10)
directly, we take an indirect approach that relates $p_{\text{ss}}(r)$ to a
function called the backward committor $q^{-}(r)$, which is the probability
that, going back in time, a trajectory starting at $r$ will reach the reactant
well before the product well. As we will see, $p_{\text{ss}}(r)$, a
nonequilibrium quantity, is related to the behavior of $q^{-}(r)$ at
equilibrium by a surprisingly simple equation. Since Brownian dynamics under
detailed balance is time-reversible (in the sense that a time-reversed process
is statistically indistinguishable from the original one), it is more
convenient in this context to think of $q^{-}(r)$ as representing the
probability that, going forward in time, a trajectory starting at $r$ will
reach the reactant well before the product well. The choice of the superscript
notation is intended to distinguish $q^{-}$ from $q^{+}$, the forward
committor, which is the probability that, going forward in time, a trajectory
starting at $r$ will reach the product well before the reactant well; for
time-reversible processes at equilibrium the two quantities are related by
$q^{+}(r)+q^{-}(r)=1$. The committor is also known as splitting probability,
or, in the protein folding literature, $p_{\text{fold}}$. One should be
careful not to confuse the definition of time reversibility in the theory of
stochastic process with the concept of reversibility in Hamiltonian dynamics
(i.e., flipping both time and momentum leaves the process invariant), the
concept of reversibility in thermodynamics, or the concept of reversible
reactions in kinetics.

In this section we show that

$q^{-}(r)=p_{\text{ss}}(r)/\pi(r)$ (42)

First, as in the one-dimensional case, the steady-state distribution
$p_{\text{ss}}(r)$ satisfies the following boundary value problem,

$\begin{cases}L^{\dagger}p_{\text{ss}}(r)=0,&r\in\Omega\backslash(A\cup B)\\\
p_{\text{ss}}(r)=\pi(r),&r\in A\\\ p_{\text{ss}}(r)=0,&r\in B\end{cases}$ (43)

Second, similar to the MFPT, the backward committor also satisfies a boundary
value problem,

$\begin{cases}Lq^{-}(r)=0,&r\in\Omega\backslash(A\cup B)\\\ q^{-}(r)=1,&r\in
A\\\ q^{-}(r)=0,&r\in B\end{cases}$ (44)

In general, the time evolution of the backward committor is actually dictated
by $L^{R}$, the generator for the time-reversed process, but $L^{R}=L$ for
time-reversible processes. Again, we will not attempt to derive (44)
rigorously. An “intuitive” justification is as follows: consider the same
hypothetical procedure as before, where an ensemble of system trajectories are
prepared, all with the same initial condition at position $r(0)=r_{0}$ outside
of $A$ or $B$; we can calculate $q^{-}(r_{0})$ by determining the fraction of
these trajectories that reach $A$ first before $B$. The fate of each
individual realization of the system dynamics is time-independent: a
trajectory either eventually reaches $A$ first before it reaches $B$, or it
doesn’t; therefore, the fraction $q^{-}(r_{0})$ is time-independent, and thus
$Lq^{-}(r_{0})$, which gives the time derivative, is zero. If $r_{0}$ is
already in $A$, then $q^{-}=1$ by construction since it is not possible for
the trajectory to reach $B$ before $A$; similarly, if $r_{0}$ is already in
$B$, then $q^{-}=0$ since the trajectory is already in $B$ before it can reach
$A$.

Comparing (13) with (14), we see that the generator is related to the Fokker-
Planck operator by

$L=\pi(r)^{-1}L^{\dagger}\pi(r)$ (45)

Substituting (45) in (44) gives $\pi(r)^{-1}L^{\dagger}\pi(r)q^{-}(r)=0$; this
implies that $L^{\dagger}\pi(r)q^{-}(r)=0$ since $\pi^{-1}(r)\propto e^{\beta
W(r)}$ is strictly positive for all $r$. Furthermore, $\pi(r)q^{-}(r)=\pi(r)$
for $r\in A$ and $\pi(r)q^{-}(r)=0$ for $r\in B$. Taken together, we see that
$\pi(r)q^{-}(r)$ satisfies the same well-posed boundary value problem (43) as
$p_{\text{ss}}(r)$, which implies that they are the same function, as desired.
We note here that this proof does not work for stochastic processes described
by the Langevin equation outside the overdamped regime.

###  3.2 An ansatz for the backward committor

Given (42), the problem of determining $p_{\text{ss}}(r)$ now reduces to that
of determining $q^{-}(r)$, the backward committor. In this section, we seek an
analytical expression for the backward committor near the saddle point
$r^{\ddagger}$. To make the derivation clearer, we first perform a change of
variable $s=r-r^{\ddagger}$ to shift the origin to the saddle point; this
procedure has no effect on the functional form of any preceding expressions
except for changing the independent variable from $r$ to $s$. At steady state,
use the fact that $p_{\text{ss}}(s)$ is time-independent and
$p_{\text{ss}}(s)=q^{-}(s)\pi(s)$; (12) becomes

$\displaystyle 0$
$\displaystyle=\sum_{ij}\partial_{i}D_{ij}\left[\beta\partial_{j}W(s)+\partial_{j}\right]q^{-}(s)\pi(s)$
(46) $\displaystyle 0$
$\displaystyle=\sum_{ij}\partial_{i}D_{ij}\left[\beta(\partial_{j}W(s))q^{-}(s)\pi(s)+(\partial_{j}q^{-}(s))\pi(s)-\beta\pi
q^{-}(s)\partial_{j}W(s)\right]$ (47) $\displaystyle 0$
$\displaystyle=\sum_{ij}D_{ij}\partial_{i}\left[(\partial_{j}q^{-}(s))\pi(s)\right]$
(48) $\displaystyle 0$
$\displaystyle=\sum_{ij}D_{ij}\left[\partial_{ij}q^{-}(s)-\beta\left(\partial_{i}W(s)\right)\left(\partial_{j}q^{-}(s)\right)\right]$
(49)

In the second equality, we used the fact that
$\partial_{i}\pi(s)=-\beta\pi(s)\partial_{i}W(s)$.

To examine the solution to (49) near the saddle point $r^{\ddagger}$, we first
consider the one dimensional case. Here, (40) becomes $W(s)=\Delta
W^{\ddagger}-\kappa^{\ddagger}s^{2}/2$ as in (25), with $\kappa^{\ddagger}>0$
representing the force constant for the unstable mode, and (49) reduces to

$\frac{d^{2}}{ds^{2}}q^{-}(s)+\beta\kappa^{\ddagger}s\frac{d}{ds}q^{-}(s)=0$
(50)

Equation (50) can be solved by first solving a first-order ODE for
$\nu(s)=dq^{-}/ds$ using an integrating factor and then integrating both sides
of the solution from $0$ to $s$. After some algebra, we arrive at the general
solution in the form of an error function

$q^{-}(s)=A+\frac{B}{\sqrt{2\pi}}\int_{0}^{s\sqrt{\beta\kappa^{\ddagger}}}e^{-u^{2}/2}\,du$
(51)

where $A$ and $B$ are two free constants. Since $q^{-}(s)$ is a backward
committor, we require that $q^{-}\to 1$ as $s\to-\infty$ (i.e., for $s$ in the
reactant well) and $q^{-}\to 0$ as $s\to\infty$ (i.e., for $s$ in the product
well). Substituting these limits into (51) gives $A-B/2=1$ and $A+B/2=0$; it
should be obvious that $A=1/2$ and $B=-1$, which means that the solution to
(50) that satisfies the appropriate boundary conditions is

$q^{-}(s)=\frac{1}{2}-\frac{1}{\sqrt{2\pi}}\int_{0}^{s\sqrt{\beta\kappa^{\ddagger}}}e^{-u^{2}/2}\,du=\frac{1}{\sqrt{2\pi}}\int_{s\sqrt{\beta\kappa^{\ddagger}}}^{\infty}e^{-u^{2}/2}\,du=\sqrt{\frac{\beta\kappa^{\ddagger}}{2\pi}}\int_{s}^{\infty}e^{-\beta\kappa^{\ddagger}u^{2}/2}\,du$
(52)

This solution indicates that the committor drops down from 1 to 0 with a
sigmoidal functional form as $s$ increases. In particular, at $s=0$,
$q^{-}=1/2$.

This analysis of the one-dimensional backward committor provides important
hints about the behavior of the committor in higher dimensions. In this
setting, the set of all points $s$ for which $q^{-}(s)$ has the same value
constitutes an isocommittor surface (or curve, in two dimensions). Near the
transition state, it is not unreasonable to assume that these isocommittor
surfaces can be approximated as planes (or lines, in two dimensions) if
$q^{-}$ is sufficiently smooth. Furthermore, within the region where the
harmonic approximation is valid, this family of isocommittor planes needs to
be parallel, since level surfaces cannot intersect for a well-defined
function.

These considerations motivate us to devise a solution to (49) in terms of
isocommittor planes. In general, a plane passing through the origin can be
parameterized using its normal vector $n$ as $n^{T}s=0$; for planes that do
not cross the origin, they can still be parameterized as $n^{T}(s-s_{0})=0$,
where $s_{0}$ is a displacement vector from some point on the plane to the
origin. In light of this functional form, let us define a vector $v_{+}$ that
is normal to the isocommittor planes (the reason for this notation will become
clear shortly); using this vector, the family of isocommittor planes can be
parameterized as $v_{+}^{T}s=u$, $u\in\mathbb{R}$. We orient $v_{+}$ such that
$q^{-}$ approaches 1 for $u<0$ and approaches 0 for $u>0$.

Taken together, the analysis so far suggests that the solution to (49) has the
form

$q^{-}(s)=\frac{a}{\sqrt{2\pi}}\int_{v_{+}^{T}s}^{\infty}e^{-(au)^{2}/2}\,du$
(53)

where $a>0$ is a free constant to be determined. It is easy to check that with
this solution form, $q^{-}(s)$ is constant on each isocommittor plane,
approaches 1 for $s$ near the reactant well, and approaches 0 for $s$ near to
the product well.

One particular isocommittor plane of interest is the one that crosses the
saddle point; i.e., the plane parameterized by $u=0$ and where $q^{-}=1/2$.
This plane is sometimes known as the stochastic separatrix and often plays a
special role in the analysis of reaction mechanisms: if we define the
transition state (or more accurately, the transition state ensemble) as the
set of configurations where the committor is 1/2, the transition state is the
stochastic saparatrix. In Kramers’ theory, this plane reduces to the peak of
the barrier separating the reactant and product state. Unfortunately, the
simplicity of this equivalence between the transition state and the barrier
peak often does not hold in practice, depending on the choice of the reaction
coordinate and the theoretical framework under consideration; see, e.g., [10,
13].

###  3.3 Properties of the vector normal to the isocommittor planes

The ansatz in (53) still leaves some questions unanswered. Specifically,

  1. 1.

How does the direction of $v_{+}$ relate to the geometry at the saddle point?

  2. 2.

What is the magnitude of $v_{+}$?

  3. 3.

What is the value of the free constant $a>0$?

We address these questions in this section.

####  3.3.1 v+ as a generalized eigenvector

To understand the direction of $v_{+}$, we first argue that $v_{+}$ is
actually an eigenvector of the matrix $H^{\ddagger}D$. In order to show this,
we need to examine the solution of (49) on the $q^{-}=1/2$ plane.

First, using the Leibniz integral rule, we record here the first- and second-
order derivatives of (53) with respect to $s$:

$\partial_{i}q^{-}(s)=-\frac{a}{\sqrt{2\pi}}e^{-\left(av_{+}^{T}s\right)^{2}/2}v_{+i}\quad\text{and}\quad\partial_{ij}q^{-}(s)=\frac{a^{3}}{\sqrt{2\pi}}\left(v_{+}^{T}s\right)e^{-\left(av_{+}^{T}s\right)^{2}/2}v_{+i}v_{+j}$
(54)

On the $q^{-}=1/2$ plane, $v_{+}^{T}s=0$ and (54) reduces to

$\partial_{i}q^{-}(s)=-\frac{a}{\sqrt{2\pi}}v_{+i}\quad\text{and}\quad\partial_{ij}q^{-}(s)=0\quad\text{when}\quad
v_{+}^{T}s=0$ (55)

Substituting in the derivative of $W(s)$ in (41) and the derivatives of
$q^{-}(s)$ in (55), we see that (49) can be further simplified on the
$q^{-}=1/2$ plane as,

$\displaystyle 0$
$\displaystyle=\sum_{ij}D_{ij}\left[-\beta\left(\sum_{k}H^{\ddagger}_{ik}s_{k}\right)\left(-\frac{a}{\sqrt{2\pi}}v_{+j}\right)\right]$
(56) $\displaystyle 0$
$\displaystyle=\sum_{k}s_{k}\left[\sum_{ij}H^{\ddagger}_{ki}D_{ij}v_{+j}\right]$
(57) $\displaystyle 0$ $\displaystyle=s^{T}\left(H^{\ddagger}Dv_{+}\right)$
(58)

In the second equality, we used the fact that $H^{\ddagger}$ is symmetric
(i.e., $H^{\ddagger}_{ik}=H^{\ddagger}_{ki}$). This result indicates that the
vector $H^{\ddagger}Dv_{+}$ is orthogonal to $s$. Since $v_{+}$ is orthogonal
to $s$ confined on the $q^{-1}=1/2$ plane, $H^{\ddagger}Dv_{+}$ is parallel
with $v_{+}$. This shows that $v_{+}$ is an eigenvector of $H^{\ddagger}D$.
Anticipating the analysis in the next section, we denote the corresponding
eigenvalue as $-\lambda_{+}$ (with $\lambda_{+}>0$).

Note that we can write the eigenvector equation
$H^{\ddagger}Dv_{+}=-\lambda_{+}v_{+}$ as

$H^{-\ddagger}v_{+}=-\lambda_{+}^{-1}Dv_{+}$ (59)

where, through an abuse of notation, we have written
$\left(H^{\ddagger}\right)^{-1}=H^{-\ddagger}$. Equations of this form are
called generalized eigenvalue problems whenever $H^{-\ddagger}$ is symmetric
and $D$ is symmetric positive definite. The collection of all such generalized
eigenvectors form a nonsingular matrix $V$ that simultaneously diagonalizes
$H^{\ddagger}$ and $D$, in the sense that

$V^{T}DV=I\quad\text{and}\quad V^{T}H^{-\ddagger}V=\Lambda^{-1}$ (60)

where $\Lambda^{-1}$ is a diagonal matrix containing the corresponding
generalized eigenvalues. The vectors $v_{i}$’s in $V$ are orthogonal with
respect to the inner product induced by $D$; i.e.,

$\left<v_{i},v_{j}\right>_{D}=v_{i}^{T}Dv_{j}=\lambda_{j}v_{i}^{T}H^{-\ddagger}v_{j}=\lambda_{j}\Lambda_{ij}^{-1}=\delta_{ij}$
(61)

where $\delta_{ij}$ is the Kronecker delta function.

####  3.3.2 The value of the free constant a

The magnitudes of $a$ and $v_{+}$ are related; first we determine the value of
$a$ in terms of $v_{+}$. To do so, we consider the solution to (49) outside
the committor 1/2 plane where $v_{+}^{T}s\neq 0$:

$\displaystyle 0$
$\displaystyle=\sum_{ij}D_{ij}\left[\frac{a^{3}}{\sqrt{2\pi}}\left(v_{+}^{T}s\right)e^{-\left(av_{+}^{T}s\right)^{2}/2}v_{+i}v_{+j}-\beta\left(\sum_{k}H^{\ddagger}_{ik}s_{k}\right)\left(-\frac{a}{\sqrt{2\pi}}e^{-(av_{+}^{T}s)^{2}/2}v_{+j}\right)\right]$
(62) $\displaystyle 0$
$\displaystyle=\sum_{ij}D_{ij}\left[a^{2}\left(v_{+}^{T}s\right)v_{+i}v_{+j}+\beta\left(\sum_{k}H^{\ddagger}_{ik}s_{k}\right)v_{+j}\right]$
(63) $\displaystyle 0$
$\displaystyle=a^{2}(v_{+}^{T}s)\sum_{ij}v_{+i}D_{ij}v_{+j}+\beta\sum_{jik}v_{+j}D_{ji}H^{\ddagger}_{ik}s_{k}$
(64) $\displaystyle 0$
$\displaystyle=a^{2}(v_{+}^{T}s)v_{+}^{T}Dv_{+}+\beta(v_{+}^{T}DH^{\ddagger})s$
(65) $\displaystyle 0$
$\displaystyle=a^{2}(v_{+}^{T}s)v_{+}^{T}Dv_{+}-\beta(\lambda_{+}v_{+}^{T})s$
(66) $\displaystyle 0$ $\displaystyle=a^{2}v_{+}^{T}Dv_{+}-\beta\lambda_{+}$
(67) $\displaystyle a$ $\displaystyle=\sqrt{\beta\lambda_{+}/v_{+}^{T}Dv_{+}}$
(68)

Applying (68) to (53) gives

$q^{-}(s)=\sqrt{\frac{\beta\lambda_{+}}{2\pi
v_{+}^{T}Dv_{+}}}\int_{v_{+}^{T}s}^{\infty}\exp\left(-\frac{1}{2}\frac{\beta\lambda^{+}}{v_{+}^{T}Dv_{+}}u^{2}\right)\,du$
(69)

Note that (68) implies that $-\lambda_{+}$ is a negative eigenvalue of
$H^{\ddagger}D$: because $\beta>0$ and $v_{+}^{T}Dv_{+}>0$ (since $D$ is
symmetric positive definite), $\lambda_{+}$ also needs to be positive so that
$a>0$. What does this result mean for the direction represented by the
corresponding eigenvector $v_{+}$? When $D$ can be written in the form of
$\sigma^{2}I$, where $I$ is an identity matrix, the diffusion on the potential
(of mean force) is isotropic (i.e., same diffusivity in every direction). In
this case, $v_{+}$ is simply an eigenvector of $H^{\ddagger}$,

$H^{\ddagger}v_{+}=-\lambda_{+}\sigma^{-2}v_{+}$ (70)

that is, $v_{+}$ is an eigenvector associated with the only negative
eigenvalue of $H^{\ddagger}$ corresponding to the unstable mode separating the
reactant and product well. In the case of anisotropic diffusion,
$H^{\ddagger}D$ still has a single unique unstable mode represented by
$v_{+}$: because $H^{\ddagger}$ is symmetric and $D$ is symmetric positive
definite, a result in linear algebra states that the eigenvalues of
$H^{\ddagger}D$ have the same signs as those of $H^{\ddagger}$. As such,
$v_{+}$ is still an eigenvector associated with the only negative eigenvalue
of $H^{\ddagger}D$ in the anisotropic case. However, the presence of the
matrix $D$ effectively stretches and rotates the dynamics on the surface of
the potential; as such, the unstable mode of $H^{\ddagger}D$ no longer
necessarily points in the same direction as the unstable mode of
$H^{\ddagger}$ (see, e.g., Figure 3A).

####  3.3.3 The magnitude of v+

The magnitude of $v_{+}$ determines how fast the committor $q^{-}(s)$ decays
to zero from the reactant well to the product well. To determine this
quantity, we first note that the derivative of the one-dimensional backward
committor in (52) at the saddle point $s=0$ is given by

$\frac{d}{ds}q^{-}(s)\Big{|}_{s=0}=-\sqrt{\frac{\beta\kappa^{\ddagger}}{2\pi}}e^{-\beta\kappa^{\ddagger}s^{2}/2}\Big{|}_{s=0}=-\sqrt{\frac{\beta\kappa^{\ddagger}}{2\pi}}$
(71)

While the gradient of the multi-dimensional backward committor in (69) at the
saddle point $s=0$ is given by

$\nabla_{s}q^{-}(s)\Big{|}_{s=0}=-\sqrt{\frac{\beta\lambda_{+}}{2\pi
v_{+}^{T}Dv_{+}}}\exp\left(-\frac{1}{2}\frac{\beta\lambda^{+}}{v_{+}^{T}Dv_{+}}(v_{+}^{T}s)^{2}\right)\Big{|}_{s=0}v_{+}=-\sqrt{\frac{\beta\lambda_{+}}{2\pi
v_{+}^{T}Dv_{+}}}v_{+}$ (72)

Comparing these two expressions, we can identify $\kappa^{\ddagger}$ with
$\lambda_{+}$ since both are “eigenvalues” corresponding to the unstable mode
at the saddle point. In order for the absolute value of (71) and the norm of
(72) (induced by $D$) to be equal, we require that

$v_{+}^{T}Dv_{+}=1$ (73)

Recall that the expression $v_{+}^{T}Dv_{+}$ is the inner product of $v_{+}$
with itself induced by $D$ and defines a vector norm via
$\|v_{+}\|_{D}=\sqrt{v_{+}^{T}Dv_{+}}$. Taken together, we see that the
solution to (49) is

$q^{-}(s)=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}\int_{v_{+}^{T}s}^{\infty}e^{-\beta\lambda_{+}u^{2}/2}\,du$
(74)

See Figure 3B and 3C for a comparison of (74) with the exact solution to (44)
for the example double-well potential shown in 3A.

Lastly, we note here that, with the help of (54) (after setting
$a=\sqrt{\beta\lambda_{+}}$), it is fairly straightforward to check that (74)
is in fact an exact solution to (49) under the harmonic approximation; i.e.,
(74) satisfies

$\sum_{ij}D_{ij}\partial_{ij}q^{-}(s)-\beta\sum_{ijk}s_{k}H^{\ddagger}_{ki}D_{ij}\partial_{j}q^{-}(s)=0$
(75)

in addition to satisfying the boundary conditions listed in (44). In other
words, if the potential surface $W(s)$ is exactly harmonic, the isocommittor
surfaces are indeed parallel planes parametrized by $v_{+}^{T}s=u$ for
$u\in\mathbb{R}$.

###  3.4 Properties of the steady-state flux vector field

Now with an expression for the backward committor in hand, we are ready to
write down an expression for the steady-state flux and analyze its behavior
near the saddle point.

####  3.4.1 The flux vector field has a constant direction

Substituting (42) and (74) in the multidimensional Fokker-Planck equation
(10), we see that the steady-state flux is

$\displaystyle J(s)$ $\displaystyle=-D\pi(s)\nabla_{s}q^{-}(s)$ (76)
$\displaystyle=-D\pi(s)\left(-\sqrt{\frac{\beta\lambda_{+}}{2\pi}}e^{-\beta\lambda_{+}(v_{+}^{T}s)^{2}/2}\right)\nabla_{s}(v_{+}^{T}s)$
(77)
$\displaystyle=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}e^{-\beta\lambda_{+}(v_{+}^{T}s)^{2}/2}\pi(s)Dv_{+}$
(78)

One immediate observation here is that near the saddle point, the flux has a
constant direction $Dv_{+}$.

####  3.4.2 The flux vector field is constant along its streamlines

Not only does $J(s)$ have a constant direction $Dv_{+}$, it also has a
constant magnitude along that direction. To see why, we first parameterize $s$
using an orthogonal decomposition

$s(u)=\Delta s+uDv_{+}$ (79)

where $\Delta s$ is a constant that measures the displacement from the origin
(i.e., the saddle point) orthogonal to the direction defined by $Dv_{+}$,
while $u\in\mathbb{R}$ is the independent variable that measures progression
along the $Dv_{+}$ direction. In (78), the only factor dependent on $s$ is
$e^{-\beta\lambda_{+}(v_{+}^{T}s)^{2}/2}\pi(s)$. Using the harmonic
approximation of $\pi(s)$ at near the saddle point in (40), this expression
can be written as

$\exp\left[-\frac{\beta}{2}s^{T}\left(H^{\ddagger}+\lambda_{+}v_{+}v_{+}^{T}\right)s\right]$
(80)

Using the parameterization $s(u)$, the quadratic form in this expression
becomes $(\Delta s+uDv_{+})^{T}(H^{\ddagger}+\lambda_{+}v_{+}v_{+}^{T})(\Delta
s+uDv_{+})$. This expression can be expanded, but some of the terms after the
expansion are constants independent of $u$. Keeping only the $u$-dependent
terms, we get

$u\Delta s^{T}H^{\ddagger}Dv_{+}+u\lambda_{+}\Delta
s^{T}v_{+}v_{+}^{T}Dv_{+}+uv_{+}^{T}D^{T}H^{\ddagger}\Delta s\\\
+u^{2}v_{+}^{T}D^{T}H^{\ddagger}Dv_{+}+u\lambda_{+}v_{+}^{T}D^{T}v_{+}v_{+}^{T}\Delta
s+u^{2}\lambda_{+}v_{+}^{T}D^{T}v_{+}v_{+}^{T}Dv_{+}$ (81)

Using the facts that $H^{\ddagger}$ and $D$ are symmetric,
$H^{\ddagger}Dv_{+}=-\lambda_{+}v_{+}$, and $v_{+}^{T}Dv_{+}=1$, one can show
that every term in this expression cancels, and thus $J(s)$ is constant along
the $Dv_{+}$ direction. Taken together, we can visualize the flux as a flow of
probability density along parallel streamlines near the saddle point, and the
velocity of the probability current along each streamline is constant. See
Figure 3B and 3C for a comparison of the streamlines of (78) with the exact
solution to (76) for the example double-well potential shown in 3A.

####  3.4.3 Flux surface integrals are constant over dividing surfaces

As in the one-dimensional case, we need the total flux from the reactant well
to the product well, which gives the reaction rate. This task is made somewhat
easy by a few considerations. First, we note that, for a double-well
potential, the flux-over-population steady-state flux $J(s)$ has exactly one
source in the reactant well $A$ and one sink in the product well $B$; that is,
$J(s)$ is a divergence-less vector field in $\Omega\backslash(A\cup B)$ at the
steady state (this should be obvious by setting time derivative to zero in
(10)). A consequence of this property is that the total flux from $A$ to $B$
is given by the surface integral of $J(s)$ over a dividing surface between $A$
and $B$ (i.e., a surface that partitions $\Omega$ into two sets, one of which
contains $A$ and the other contains $B$). Importantly, the particular choice
of the dividing surface is irrelevant; if this integral differs on two
dividing surfaces, it would imply that there is a net sink or source located
in the in-between region. Second, since the probability flow is concentrated
near the saddle point, one can restrict the flux surface integral further to
regions close to the saddle point where (78) is valid.

With the help of (108), a result we will prove in Section 4, a particularly
simple choice of dividing surface is the $q^{-}=1/2$ plane. The total flux of
(78) over the stochastic separatrix is

$\displaystyle J$ $\displaystyle=\int\|v_{+}\|^{-1}v_{+}^{T}J(s)\,d\sigma(s)$
(82) $\displaystyle\overset{1}{=}\int v_{+}^{T}J(s)\delta(v_{+}^{T}s)\,ds$
(83)
$\displaystyle=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}v_{+}^{T}Dv_{+}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}\int e^{{}^{-}\beta s^{T}H^{\ddagger}s/2}\delta(v_{+}^{T}s)\,ds$
(84) $\displaystyle=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}\left[(2\pi/\beta)^{1-N}v_{+}^{T}H^{-\ddagger}v_{+}\det
H^{\ddagger}\right]^{-1/2}$ (85)
$\displaystyle\overset{2}{=}\left(\frac{2\pi}{\beta}\right)^{N/2}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}\frac{\beta\lambda_{+}}{2\pi}\left|\det
H^{\ddagger}\right|^{-1/2}$ (86)

Here, equality (1) follows from the relation
$d\sigma(s)=\delta(v_{+}^{T}s)\|\nabla_{s}v_{+}^{T}s\|\,ds$, which is a
consequence of what is known as the coarea formula; note that
$\|\nabla_{s}v_{+}^{T}s\|=\|v_{+}\|$. The volume integrals are over the entire
$\Omega=\mathbb{R}^{N}$ reaction coordinate space of the system; this will be
the assumed domain of integration for the rest of the note unless otherwise
specified. Equality (2) follows from the facts that
$H^{-\ddagger}v_{+}=-\lambda_{+}^{-1}Dv_{+}$ and $\det H^{\ddagger}<0$.

In the rest of this section, we show and confirm that the total flux can be
obtained using a surface integral over an arbitrary plane near the the saddle
point, which gives the same result as (86). Let us parameterize such a plane
by $n^{T}s=\theta$, where $n$ is the normal vector (not necessarily
normalized) and $\theta\in\mathbb{R}$. The only condition we impose on the
plane is that $n$ is not orthogonal to $Dv_{+}$ (i.e., the plane parameterized
by $n^{T}s=\theta$ is not parallel to $Dv_{+}$, the direction of the flux).
Furthermore, without loss of generality we assume that $n$ is oriented such
that $n^{T}Dv_{+}>0$. Using (78) and the coarea formula, the integral of
$J(s)$ over this plane is given by

$\displaystyle J$ $\displaystyle=\int\|n\|^{-1}n^{T}J(s)\,d\sigma(s)$ (87)
$\displaystyle=\int n^{T}J(s)\delta(\theta-n^{T}s)\,ds$ (88)
$\displaystyle=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}n^{T}Dv_{+}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}\int e^{{}^{-}\beta
s^{T}\left(H^{\ddagger}+\lambda_{+}v_{+}v_{+}^{T}\right)s/2}\delta(\theta-n^{T}s)\,ds$
(89)

As with most Gaussian integrals, we seek a change of variable that makes the
volume integral in (89) separable. For a typical Gaussian integral with an
integrand of the form $e^{-x^{T}Ax/2}$, where $A$ is a symmetric matrix, there
is an orthogonal matrix $Q$ that diagonalizes $A$ (i.e., $Q^{T}AQ=M$ for some
diagonal matrix $M$) and enables a change of variable $x=Qy$ that renders the
integral separable. In the current case, recall our discussion of the
generalized eigenvector problem in Section 3.3, where we have defined the
matrix $V$ that simultaneously diagonalizes $H^{\ddagger}$ and $D$ according
to (60). Comparing $V$ and $Q$, we introduce the change of variable $t=V^{T}s$
(as well as $m=V^{T}n$). The Jacobian for this transformation is $\det
V^{-1}$, which is not necessarily one because $V$ may not be orthogonal. After
this transformation, $v_{+}$ relates to a new basis vector $e_{+}$ by
$e_{+}=V^{-1}v_{+}$; to see why, note that together with (60),
$H^{\ddagger}Dv_{+}=-\lambda_{+}v_{+}$ implies

$\Lambda V^{-1}v_{+}=-\lambda_{+}V^{-1}v_{+}$ (90)

that is, $V^{-1}v_{+}$ is an eigenvector of the diagonal matrix $\Lambda$
corresponding to its unique negative eigenvalue $-\lambda_{+}$. Since
$-\lambda_{+}$ has an algebraic multiplicity of one, the eigenvector
$V^{-1}v_{+}$ must have only one nonzero element. Furthermore, the norm of
$V^{-1}v_{+}$ is fixed by
$\|V^{-1}v_{+}\|=v_{+}V^{-T}V^{-1}v_{+}=v_{+}^{T}Dv_{+}=1$. Since $v_{+}$ is a
column vector of $V$, the nonzero element of $e_{+}$ must be positive,
otherwise $Ve_{+}=-v_{+}$. Together, we have shown that $e_{+}$ is a
normalized standard basis vector in the new coordinate system.

With this change of variable, (89) becomes

$\displaystyle J$
$\displaystyle=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}n^{T}Dv_{+}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}\int e^{{}^{-}\beta
t^{T}V^{-1}\left(H^{\ddagger}+\lambda_{+}v_{+}v_{+}^{T}\right)V^{-T}t/2}\delta(\theta-m^{T}V^{-1}V^{-T}t)|\det
V^{-1}|\,dt$ (91)
$\displaystyle=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}n^{T}Dv_{+}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}|\det V^{-1}|\int e^{{}^{-}\beta
t^{T}\left(\Lambda+\lambda_{+}e_{+}e_{+}^{T}\right)t/2}\delta(\theta-m^{T}Dt)\,dt$
(92)

Before further simplifications of (92), a comment about the term
$t^{T}\left(\Lambda+\lambda_{+}e_{+}e_{+}^{T}\right)t$ is in order. In the new
basis, $H^{\ddagger}$ becomes $\Lambda$, while $v_{+}v_{+}^{T}$ becomes
$e_{+}e_{+}^{T}$. In particular, $\Lambda+\lambda_{+}e_{+}e_{+}^{T}$ is a
diagonal matrix whose diagonal elements are the eigenvalues $\lambda_{j}$’s
except for $-\lambda_{+}$, which has been deleted from $\Lambda$ by
$\lambda_{+}e_{+}e_{+}^{T}$. This observation should become especially obvious
using the outer product form of the spectral theorem, which gives
$\Lambda=\sum_{j}\lambda_{j}e_{j}e_{j}^{T}$; in other words, compared to
$\Lambda$ (or $H^{\ddagger}$), $\Lambda+\lambda_{+}e_{+}e_{+}^{T}$ (or
$H^{\ddagger}+\lambda_{+}v_{+}v_{+}^{T}$) is rank-deficient because one of its
eigenspaces corresponding to the unstable mode has been deleted. As a result,
the quadratic form $t^{T}(\Lambda+\lambda_{+}e_{+}e_{+}^{T})t$ sums over all
$\lambda_{j}t_{j}^{2}$ except for $j$ corresponding to
$-\lambda_{+}t_{+}^{2}$. In light of this analysis, we will use notations
$\sum_{j\neq+}$ and $\prod_{j\neq+}$ to denote summation or product over $j$’s
except for $j$ corresponding to $-\lambda_{+}$ and $e_{+}$.

Let us denote the volume integral in (92) as $I$. We will simplify $I$ with
the following strategy. First, we convert the Dirac delta function in (92)
into its integral representation; i.e.,

$\delta(f(x))=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{ikf(x)}\,dk$ (93)

for any sufficiently smooth test function $f(x)$. Second, we integrate over
all $t_{j}$’s for which $j\neq+$, then over $k$, and then over $t_{+}$. This
order of operation ensures that every integral is an analytically tractable
Gaussian integral with a linear term, which can be solved by first completing
the squares:

$\displaystyle I$ $\displaystyle=\frac{1}{2\pi}\int\int e^{{}^{-}\beta
t^{T}\left(\Lambda+\lambda_{+}e_{+}e_{+}^{T}\right)t/2}e^{ik(\theta-m^{T}Dt)}\,dtdk$
(94) $\displaystyle=\frac{1}{2\pi}\int\int
e^{ik\theta}e^{-ik\sum_{l}m_{l}D_{l+}t_{+}}\prod_{j\neq+}\left(\int
e^{{}^{-}\beta\lambda_{j}t_{j}^{2}/2}e^{-ik\sum_{l}m_{l}D_{lj}t_{j}}\,dt_{j}\right)\,dkdt_{+}$
(95) $\displaystyle=\frac{1}{2\pi}\int\int
e^{ik\theta}e^{-ik\sum_{l}m_{l}D_{l+}t_{+}}\prod_{j\neq+}\left(\sqrt{\frac{2\pi}{\beta\lambda_{j}}}e^{-k^{2}\left(\sum_{l}m_{l}D_{lj}\right)^{2}/2\beta\lambda_{j}}\right)\,dkdt_{+}$
(96)
$\displaystyle=\frac{1}{2\pi}\left(\prod_{j\neq+}\sqrt{\frac{2\pi}{\beta\lambda_{j}}}\right)\int
e^{ik\theta}e^{-ik\sum_{l}m_{l}D_{l+}t_{+}}\left(\frac{2\pi\beta}{\sum_{j\neq+}\lambda_{j}^{-1}\left(\sum_{l}m_{l}D_{lj}\right)^{2}}\right)^{1/2}e^{-\frac{\beta}{2}\frac{(\theta-\sum_{l}m_{l}D_{l+}t_{+})^{2}}{\sum_{j\neq+}\lambda_{j}^{-1}\left(\sum_{l}m_{l}D_{lj}\right)^{2}}}\,dt_{+}$
(97)
$\displaystyle=\frac{1}{2\pi}\left(\prod_{j\neq+}\sqrt{\frac{2\pi}{\beta\lambda_{j}}}\right)\left(\frac{2\pi\beta}{\sum_{j\neq+}\lambda_{j}^{-1}\left(\sum_{l}m_{l}D_{lj}\right)^{2}}\right)^{1/2}\left(\frac{2\pi}{\beta}\frac{\sum_{j\neq+}\lambda_{j}^{-1}\left(\sum_{l}m_{l}D_{lj}\right)^{2}}{\left(\sum_{l}m_{l}D_{l+}\right)^{2}}\right)^{1/2}$
(98)
$\displaystyle=\left(\prod_{j\neq+}\sqrt{\frac{2\pi}{\beta\lambda_{j}}}\right)|m^{T}De_{+}|^{-1}$
(99)

Substituting (99) back into (92) gives

$\displaystyle J$
$\displaystyle=\sqrt{\frac{\beta\lambda_{+}}{2\pi}}\left(\prod_{j\neq+}\sqrt{\frac{2\pi}{\beta\lambda_{j}}}\right)Z^{-1}e^{-\beta\Delta
W^{\ddagger}}|\det V^{-1}|\frac{n^{T}Dv_{+}}{|m^{T}De_{+}|}$ (100)
$\displaystyle\overset{1}{=}\left(\frac{2\pi}{\beta}\right)^{N/2}\frac{\beta\lambda_{+}}{2\pi}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}\frac{|\det V^{-1}|}{|\det H^{\ddagger}D|^{1/2}}$ (101)
$\displaystyle\overset{2}{=}\left(\frac{2\pi}{\beta}\right)^{N/2}Z^{-1}e^{-\beta\Delta
W^{\ddagger}}\frac{\beta\lambda_{+}}{2\pi}|\det H^{\ddagger}|^{-1/2}$ (102)

In equality (1), we have used the fact that

$n^{T}Dv_{+}=m^{T}V^{-1}DVe_{+}=m^{T}V^{-1}\left(V^{-T}V^{-1}\right)Ve_{+}=m^{T}\left(V^{-1}V^{-T}\right)e_{+}=m^{T}D^{T}e_{+}=m^{T}De_{+}>0$
(103)

In equality (2), we have used the fact that

$|\det H^{\ddagger}D|=|\det H^{\ddagger}|\det D\quad\text{and}\quad\det D=\det
V^{-T}V^{-1}=\left(\det V^{-1}\right)^{2}$ (104)

both of which follow directly from (60). It should be obvious that (102) is
independent of either $\theta$ or $n$ used to parameterize the dividing plane
and is thus an intrinsic property of the flux. The result in (102) is the same
as (86), as expected. However, the choice of $q^{-}=1/2$ in the derivation of
(86) is special in that the quadratic form
$s^{T}(H^{\ddagger}+\lambda_{+}v_{+}v_{+}^{T})s$ is simply
$s^{T}H^{\ddagger}s$, which is not rank-deficient.

###  3.5 Langer’s rate constant

According to the flux-over-population method, the only quantity left to
determine for computing the rate constant is the steady-state population $n$
in the reactant well. This is given by

$n=Z^{-1}\int
e^{-\beta(r-r_{A})^{T}H_{A}(r-r_{A})/2}\,dr=Z^{-1}\prod_{i}\int_{-\infty}^{\infty}e^{-\beta\mu_{i}t_{i}^{2}/2}\,dt_{i}\\\
=Z^{-1}\left(2\pi/\beta\right)^{N/2}\prod_{i}\mu_{i}^{-1/2}=Z^{-1}\left(2\pi/\beta\right)^{N/2}\left(\det
H_{A}\right)^{-1/2}$ (105)

where $r-r_{A}=Qt$ with $Q$ being an orthogonal matrix that diagonalizes
$H_{A}$, and $\mu_{i}$’s are the eigenvalues of $H_{A}$. Recall that $H_{A}$
is assumed to be symmetric positive definite, and thus $\mu_{i}>0$ for all
$i$’s.

At last, combining (86) and (105), we see that the rate constant is

$k_{AB}=J/n=\frac{\beta\lambda_{+}}{2\pi}\sqrt{\frac{\det H_{A}}{\left|\det
H^{\ddagger}\right|}}e^{-\beta\Delta W^{\ddagger}}\quad\text{where}\quad
H^{\ddagger}Dv_{+}=-\lambda_{+}v_{+}$ (106)

Three brief comments about this result are in order. First, we note here that
it is also common in the literature to define $-\lambda_{+}$ as an eigenvalue
of the matrix $\beta H^{\ddagger}D$, while the magnitude of $v_{+}$ is still
fixed by $v_{+}^{T}Dv_{+}=1$. One can retrace our derivations and see that,
starting from (65), this has the effect of replacing all occurrences of
$\beta\lambda_{+}$ by $\lambda_{+}$, in which case the rate constant now reads

$k_{AB}=\frac{\lambda_{+}}{2\pi}\sqrt{\frac{\det H_{A}}{\left|\det
H^{\ddagger}\right|}}e^{-\beta\Delta W^{\ddagger}}\quad\text{where}\quad\beta
H^{\ddagger}Dv_{+}=-\lambda_{+}v_{+}$ (107)

Second, it is easy to check that the multidimensional result can be reduced to
Kramers’ rate constant (28) in the case of a single reaction coordinate. By
comparing (25) to (39) and (40), we identify $H_{A}$ as $\kappa_{A}$ and
$H^{\ddagger}$ as $-\kappa^{\ddagger}$. Next, since $-\lambda_{+}$ is the
eigenvalue corresponding to the unstable mode of $H^{\ddagger}D$, and the
eigenvalue of a $1\times 1$ “matrix” is the matrix element itself, it follows
that in the one-dimensional case where only the unstable mode is considered,
$-\lambda_{+}$ can be identified as the product of $H^{\ddagger}$ (i.e.,
$-\kappa^{\ddagger}$) and $D$, the one-dimensional diffusion constant. Making
these substitutions reduce (106) to (28).

Lastly, we briefly comment on some computational aspects of the theory. Using
molecular dynamics simulations, the application of (106) requires the
determination of the activation free energy $\Delta W^{\ddagger}$ and the
diffusion matrix $D$. The activation free energy can be obtained through a
variety of well-established equilibrium techniques such as umbrella sampling,
or more recent nonequilibrium techniques based on Jarzynski equality [14].
Some methods for computing the diffusion matrix are described in, e.g., [15,
16, 17, 18, 19].

###  3.6 Limitation of Langer’s theory

Some caution should be exercised in applying Langer’s theory; in this section
we briefly discuss two such pitfalls. First, Langer’s result can fail in some
unexpected cases where the assumption of separation of timescales implied by
the flux-over-population procedure is violated. In the case of highly
anisotropic diffusion, it is possible to arrange the relative positions of the
reactant and product wells such that barrier (re)crossing over the saddle
point can take place on a timescale much shorter than the time required to
relax within either well due to slower within-well diffusion. If the diffusion
anisotropy is extreme, the overall kinetics could become nonexponential,
depending on the initial preparation of the system within the reactant well.
In such cases, Langer’s theory predicts a rate constant that cannot properly
account for the dynamics of within-well equilibration. This issue was first
described in [20], and a corrected analytical expression for the rate constant
in the two-dimensional case was derived in [21]. More generally, anisotropy in
the diffusion matrix and/or potential surface can lead to deviations from
Langer’s theory; some case studies have been documented in, e.g., [22, 23].

Second, for complex systems typically encountered in condensed-phase chemistry
and biophysics, the difficulty with Langer’s theory is not the calculations
entailed by the theory itself, but rather the question of how to select a
small number of reaction coordinates that can provide a “complete” description
of an activated rate process. This question is an active area of research
(see, e.g., [24, 25]), and in practice answering this question is often more
akin to an art that relies on domain-specific knowledge and intuition.
Typically, an incomplete description of the reaction leads to a loss of
Markovianity of the projected stochastic process in the reaction coordinate
space; mathematically, this prevents us from simplifying the generalized
Langevin equation into the Langevin equation (see Section 1.1), which was the
starting point for the derivation of Kramers’ theory. When the memory kernel
can be reasonably approximated, one approach to treating the non-Markovian
dynamics in the reaction coordinate space is the Grote-Hynes theory based on
the stable states picture [26], which is akin to Kramers’ theory but with the
generalized Langevin equation as its starting point; a multidimensional
generalization was developed in [27]. Recent theoretical development has
focused on alternative frameworks such as transition path theory [28],
milestoning [29], weighted ensemble [30], and path sampling methods [31, 32]
that circumvent these difficulties to various extent, usually at the expense
of increased computational costs.

##  4 Simplification of Langer’s theory by Berezhkovskii and Szabo

It is often claimed that the committor is the ideal reaction coordinate to
describe a reaction. In cases where Langer’s theory is adequate, this claim is
supported by a striking result [3], which we will demonstrate in this section;
namely, that projecting the system dynamics along the direction of $v_{+}$ and
then applying Kramers’ theory gives the same rate constant as Langer’s result,
while projections along any other directions either give an overestimation or
outright do not converge. Recall from our extensive discussion in Section 3
that $v_{+}$ is the vector normal to the isocommittor planes near the saddle
point and that $v_{+}$ is the eigenvector of $H^{\ddagger}D$ corresponding to
the unique unstable mode at the saddle point.

###  4.1 Gaussian surface integral over a plane

Before we discuss the main result of this section, we take a detour here and
show that for any vectors $x,x_{0},n\in\mathbb{R}^{N}$ and an invertible,
symmetric matrix $A\in\mathbb{R}^{N\times N}$,

$\int_{\Omega}\delta(\theta-n^{T}x)e^{-\beta(x-x_{0})^{T}A(x-x_{0})/2}\,dx=\frac{e^{-\beta(\theta-n^{T}x_{0})^{2}/2n^{T}A^{-1}n}}{\sqrt{(2\pi/\beta)^{1-N}(n^{T}A^{-1}n)\det
A}}$ (108)

whenever $n^{T}A^{-1}n\det A>0$. This is a Gaussian integral restricted to a
plane parameterized by $n^{T}x=\theta$, with $n$ being the normal vector.

Our strategy for evaluating (108) is somewhat similar to the approach taken in
Section 3.4 for evaluating the surface integral of the steady-state flux
vector; one major difference here is that (108) does not involve a rank-
deficient matrix. Again, we start the derivation by seeking a change of
variable that makes the volume integral in (108) separable. Since $A$ is
symmetric, there exists a diagonal matrix $\Lambda$ and orthogonal matrix $Q$
such that $Q^{T}AQ=\Lambda$. Here, the diagonal elements in $\Lambda$ contain
the eigenvalues of $A$ and the columns of $Q$ are the corresponding
eigenvectors. Let us define the change of variable $x=Qy$, along with
$x_{0}=Qy_{0}$ and $n=Qm$. The Jacobian of this transformation is $\det Q=1$.

Before we apply the change of variable to (108), let us pause for a moment and
consider the meaning of the condition $n^{T}A^{-1}n\det A>0$. Since $A$ is
invertible and symmetric, all the eigenvalues of $A$ are real and nonzero. If
$A$ is furthermore positive definite, then $\det A=\prod\lambda_{j}>0$ and
$n^{T}A^{-1}n>0$ and thus (108) converges for any nonzero
$n\in\mathbb{R}^{N}$. We are more interested in the scenario where $A$ has a
single unstable mode. Let us denote this eigenvector as $v_{j^{\prime}}$ and
the corresponding negative eigenvalue as $\lambda_{j^{\prime}}$. In this case,
$\det A<0$ and thus we require $n^{T}A^{-1}n<0$; this implies that $n$ should
not be too close to being perpendicular to $v_{j^{\prime}}$ (i.e., the
integration should not be over a plane close to being parallel to the unstable
mode). To see why, let us write $n=\sum_{j}m_{j}v_{j}$ using the eigenvectors
$v_{j}$’s in $Q$ as the basis vectors. With this representation,
$n^{T}A^{-1}n=\sum_{j}\lambda_{j}^{-1}m_{j}^{2}$; this sum is negative only if
$m_{j^{\prime}}$ is large enough such that
$|\lambda_{j^{\prime}}^{-1}|m_{j^{\prime}}^{2}>\sum_{j\neq
j^{\prime}}\lambda_{j}^{-1}m_{j}^{2}$. In the following analysis we will prove
(108) in the case of a single unstable mode. For readers not interested in the
algebraic details, the rest of the section can be skipped without loss of
continuity.

Now, let us denote the integral in (108) as $I$. With the change of variable,
we see that

$\displaystyle I$
$\displaystyle=\int\delta(\theta-m^{T}Q^{T}Qy)e^{-\beta(y-y_{0})^{T}\Lambda(y-y_{0})/2}\,dy$
(109) $\displaystyle=\frac{1}{2\pi}e^{-\beta y_{0}^{T}\Lambda y_{0}/2}\int\int
e^{ik(\theta-m^{T}y)}e^{-\beta y^{T}\Lambda y/2}e^{\beta y_{0}^{T}\Lambda
y}\,dydk$ (110) $\displaystyle=\frac{1}{2\pi}e^{-\beta y_{0}^{T}\Lambda
y_{0}/2}\int e^{ik\theta}\int e^{-\beta y^{T}\Lambda y/2}e^{(\beta
y_{0}-ik\Lambda^{-1}m)^{T}\Lambda y}\,dydk$ (111)
$\displaystyle=\frac{1}{2\pi}e^{-\beta y_{0}^{T}\Lambda y_{0}/2}\int
e^{ik\theta}\prod_{j}\int
e^{-\frac{\lambda_{j}}{2\beta}\left[\beta(y_{j}-y_{0j})+ik\lambda_{j}^{-1}m_{j}\right]^{2}+\frac{\lambda_{j}}{2\beta}\left[ik\lambda_{j}^{-1}m_{j}-\beta
y_{0j}\right]^{2}}\,dy_{j}dk$ (112) $\displaystyle=\frac{1}{2\pi}e^{-\beta
y_{0}^{T}\Lambda y_{0}/2}\int
e^{ik\theta}\prod_{j}e^{\frac{\lambda_{j}}{2\beta}\left[ik\lambda_{j}^{-1}m_{j}-\beta
y_{0j}\right]^{2}}\left(\int
e^{-\frac{\lambda_{j}}{2\beta}\left[\beta(y_{j}-y_{0j})+ik\lambda_{j}^{-1}m_{j}\right]^{2}}\,dy_{j}\right)\,dk$
(113)

The second equality follows by using the integral representation of the Dirac
delta function as in (93).

At this point, we split the integrand of $\int dk$ in (113) into the product
of two expressions. The first involves all $j$’s for which $\lambda_{j}>0$,
and the second consists of $j^{\prime}$ for which $\lambda_{j^{\prime}}<0$.
The first group is

$\prod_{j\neq
j^{\prime}}e^{\frac{\lambda_{j}}{2\beta}\left[ik\lambda_{j}^{-1}m_{j}-\beta
y_{0j}\right]^{2}}\int
e^{-\frac{\lambda_{j}}{2\beta}\left[\beta(y_{j}-y_{0j})+ik\lambda_{j}^{-1}m_{j}\right]^{2}}\,dy_{j}=\prod_{j\neq
j^{\prime}}e^{\frac{\lambda_{j}}{2\beta}\left[ik\lambda_{j}^{-1}m_{j}-\beta
y_{0j}\right]^{2}}\sqrt{\frac{2\pi}{\beta\lambda_{j}}}\\\ =\left(\prod_{j\neq
j^{\prime}}\sqrt{\frac{2\pi}{\beta\lambda_{j}}}\right)e^{\frac{\beta}{2}\sum_{j\neq
j^{\prime}}\lambda_{j}y_{0j}^{2}}e^{\sum_{j\neq
j^{\prime}}\left(-k^{2}\lambda_{j}^{-1}m_{j}^{2}/2\beta-ikm_{j}y_{0j}\right)}$
(114)

While the second group is simply

$\int e^{ik\theta-
ikm_{j^{\prime}}y_{j^{\prime}}-\beta\lambda_{j^{\prime}}y_{j^{\prime}}^{2}/2+\beta\lambda_{j^{\prime}}y_{0j^{\prime}}y_{j^{\prime}}}\,dy_{j^{\prime}}$
(115)

The original integral in (108) can now be rewritten as

$\displaystyle I$
$\displaystyle=\frac{1}{2\pi}e^{-\frac{\beta}{2}y_{0}^{T}\Lambda
y_{0}}\left(\prod_{j\neq
j^{\prime}}\sqrt{\frac{2\pi}{\beta\lambda_{j}}}\right)e^{\frac{\beta}{2}\sum_{j\neq
j^{\prime}}\lambda_{j}y_{0j}^{2}}\int\int
I(k)I(y_{j^{\prime}})\,dy_{j^{\prime}}dk$ (116)
$\displaystyle=\frac{1}{2\pi}e^{-\frac{\beta}{2}\lambda_{j^{\prime}}y_{0j^{\prime}}^{2}}(2\pi/\beta)^{N-1}\left(\prod_{j\neq
j^{\prime}}\lambda_{j}^{-1/2}\right)\int I(y_{j^{\prime}})\int
I(k)\,dkdy_{j^{\prime}}$ (117)

with

$\displaystyle I(k)$ $\displaystyle=e^{ik\theta-
ikm_{j^{\prime}}y_{j^{\prime}}+\sum_{j\neq
j^{\prime}}\left(-k^{2}\lambda_{j}^{-1}m_{j}^{2}/2\beta-ikm_{j}y_{0j}\right)}$
(118) $\displaystyle I(y_{j^{\prime}})$
$\displaystyle=e^{-\frac{\beta}{2}\lambda_{j^{\prime}}y_{j^{\prime}}^{2}+\beta\lambda_{j^{\prime}}y_{0j^{\prime}}y_{j^{\prime}}}$
(119)

To further simplify this expression, we first integrate over $k$, which is a
Gaussian integral with a linear term,

$\int I(k)\,dk=\left(\frac{2\pi\beta}{\sum_{j\neq
j^{\prime}}\lambda_{j}^{-1}m_{j}^{2}}\right)^{1/2}e^{-\frac{\beta}{2}\frac{\left(\sum_{j\neq
j^{\prime}}m_{j}y_{0j}+m_{j^{\prime}}y_{j^{\prime}}-\theta\right)^{2}}{\sum_{j\neq
j^{\prime}}\lambda_{j}^{-1}m_{j}^{2}}}$ (120)

Then we integrate over $y_{j^{\prime}}$, which is again a Gaussian integral
with a linear term,

$\displaystyle\int I(y_{j^{\prime}})\int I(k)\,dkdy_{j^{\prime}}$
$\displaystyle=\left(\frac{2\pi\beta}{\sum_{j\neq
j^{\prime}}\lambda_{j}^{-1}m_{j}^{2}}\right)^{1/2}\int
e^{-\frac{\beta}{2}\lambda_{j^{\prime}}y_{j^{\prime}}^{2}+\beta\lambda_{j^{\prime}}y_{0j^{\prime}}y_{j^{\prime}}-\frac{\beta}{2}\frac{\left(\sum_{j\neq
j^{\prime}}m_{j}y_{0j}+m_{j^{\prime}}y_{j^{\prime}}-\theta\right)^{2}}{\sum_{j\neq
j^{\prime}}\lambda_{j}^{-1}m_{j}^{2}}}\,dy_{j^{\prime}}$ (121)
$\displaystyle=\left(\frac{2\pi\beta}{\sum_{j\neq
j^{\prime}}\lambda_{j}^{-1}m_{j}^{2}}\frac{2\pi\lambda_{j^{\prime}}^{-1}\sum_{j\neq
j^{\prime}}\lambda_{j}^{-1}m_{j}^{2}}{\beta\sum_{j}\lambda_{j}^{-1}m_{j}^{2}}\right)^{1/2}e^{\frac{\beta}{2}\lambda_{j^{\prime}}y_{0j^{\prime}}^{2}}e^{-\frac{\beta}{2}\frac{\left(\theta-\sum_{j}m_{j}y_{0j}\right)^{2}}{\sum_{j}\lambda_{j}^{-1}m_{j}^{2}}}$
(122)
$\displaystyle=2\pi\left(\lambda_{j^{\prime}}n^{T}A^{-1}n\right)^{-1/2}e^{\frac{\beta}{2}\lambda_{j^{\prime}}y_{0j^{\prime}}^{2}}e^{-\beta(\theta-n^{T}x_{0})^{2}/2n^{T}A^{-1}n}$
(123)

The last equality follows from the fact that
$n^{T}A^{-1}n=m^{T}\Lambda^{-1}m=\sum_{j}\lambda_{j}^{-1}m_{j}^{2}$.

Taken together, we see that

$\displaystyle I$
$\displaystyle=e^{-\frac{\beta}{2}\lambda_{j^{\prime}}y_{0j^{\prime}}^{2}}(2\pi/\beta)^{N-1}\left(\prod_{j\neq
j^{\prime}}\lambda_{j}^{-1/2}\right)\left(\lambda_{j^{\prime}}n^{T}A^{-1}n\right)^{-1/2}e^{\frac{\beta}{2}\lambda_{j^{\prime}}y_{0j^{\prime}}^{2}}e^{-\beta(\theta-n^{T}x_{0})^{2}/2n^{T}A^{-1}n}$
(124)
$\displaystyle=(2\pi/\beta)^{N-1}\left(\prod_{j}\lambda_{j}^{-1/2}\right)\left(n^{T}A^{-1}n\right)^{-1/2}e^{-\beta(\theta-n^{T}x_{0})^{2}/2n^{T}A^{-1}n}$
(125)
$\displaystyle=\frac{e^{-\beta(\theta-n^{T}x_{0})^{2}/2n^{T}A^{-1}n}}{\sqrt{(2\pi/\beta)^{1-N}(n^{T}A^{-1}n)\det
A}}$ (126)

as desired.

###  4.2 Projection of Langer’s result to one dimension

Let us consider the same $N$-dimensional potential of mean force $W(r)$ as in
Section 3. Here, we are interested in further projecting $W(r)$ to a one-
dimensional potential of mean force. Let us denote this direction by a vector
$n$. Integrating away the degrees of freedom orthogonal to $n$ is equivalent
to performing a surface integral of the probability density $\pi(r)$ over a
family of planes that are orthogonal to $n$. Let us parameterize this family
of planes as $n^{T}r=\theta$ for $\theta\in\mathbb{R}$. After having worked
through our derivation of Langer’s result, it should be obvious to the reader
that the one-dimensional potential of mean force is

$e^{-\beta F(\theta)}=Z^{-1}\int\delta(\theta-n^{T}r)e^{-\beta W(r)}\,dr$
(127)

Adapting Kramers’ result in (24) to $F(\theta)$, the one-dimensional rate
constant is

$k(n)=\left(\int_{\theta_{A}}^{\theta_{B}}\frac{e^{\beta
F(\theta)}}{n^{T}Dn}\,d\theta\int_{-\infty}^{\theta^{\ddagger}}e^{-\beta
F(\theta)}\,d\theta\right)^{-1}$ (128)

where $\theta_{A}$ and $\theta_{B}$ are the positions of free energy minima at
the reactant and product well, respectively, and $\theta^{\ddagger}$ is the
position of the saddle point in the one-dimensional projection. The expression
$n^{T}Dn$ is the one-dimensional diffusion constant along $n$. To understand
this expression, first note that since the diffusion matrix represents a
physical property of the system, the representation of $D$ should change with
coordinate transformations in such a way as to leave the underlying physics
invariant. This makes $D$ a tensor, specifically a second-order contravariant
tensor. Under the transformation $r\mapsto\theta$, the tensor $D_{ij}$
transforms correspondingly to a zeroth-order tensor (i.e., a constant) by

$\sum_{ij}D_{ij}\frac{\partial\theta}{\partial
r_{i}}\frac{\partial\theta}{\partial r_{j}}=\sum_{ij}n_{i}D_{ij}n_{j}=n^{T}Dn$
(129)

as desired.

Using (108) and the harmonic approximations at the saddle point (40), the
first integral in (128) evaluates to

$\displaystyle\int_{\theta_{A}}^{\theta_{B}}\frac{e^{\beta
F(\theta)}}{n^{T}Dn}\,d\theta$ $\displaystyle\approx Z\frac{e^{\beta\Delta
W^{\ddagger}}}{n^{T}Dn}\int\left(\int\delta(\theta-n^{T}r)e^{-\beta(r-r^{\ddagger})^{T}H^{\ddagger}(r-r^{\ddagger})/2}\,dr\right)^{-1}d\theta$
(130) $\displaystyle=Z\frac{e^{\beta\Delta
W^{\ddagger}}}{n^{T}Dn}\sqrt{(2\pi/\beta)^{1-N}(n^{T}H^{-\ddagger}n)\det
H^{\ddagger}}\sqrt{(2\pi/\beta)|n^{T}H^{-\ddagger}n|}$ (131)
$\displaystyle=(2\pi/\beta)^{1-N/2}Ze^{\beta\Delta
W^{\ddagger}}\frac{|n^{T}H^{-\ddagger}n|}{n^{T}Dn}|\det H^{\ddagger}|^{1/2}$
(132)

Note that since $\det H^{\ddagger}<0$, the use of (108) requires that
$n^{T}H^{-\ddagger}n<0$; otherwise the integral diverges as we discussed in
the previous section. Similarly, with the harmonic approximation at the
reactant well (39), the second integral in (128) evaluates to

$\displaystyle\int_{-\infty}^{\theta^{\ddagger}}e^{-\beta F(\theta)}\,d\theta$
$\displaystyle\approx
Z^{-1}\int\int\delta(\theta-n^{T}r)e^{-\beta(r-r_{A})^{T}H_{A}(r-r_{A})/2}\,drd\theta$
(133)
$\displaystyle=Z^{-1}\sqrt{\frac{(2\pi/\beta)(n^{T}H^{-1}_{A}n)}{(2\pi/\beta)^{1-N}(n^{T}H^{-1}_{A}n)\det
H_{A}}}$ (134) $\displaystyle=Z^{-1}(2\pi/\beta)^{N/2}\left(\det
H_{A}\right)^{-1/2}$ (135)

The application of (108) here does not impose any further conditions on $n$,
since $H_{A}$ is symmetric positive definite. Together, the rate constant is

$k(n)=\frac{\beta}{2\pi}\sqrt{\frac{\det H_{A}}{|\det
H^{\ddagger}|}}\frac{n^{T}Dn}{|n^{T}H^{-\ddagger}n|}e^{-\beta\Delta
W^{\ddagger}}$ (136)

The one-dimensional rate constant in (136) is reminiscent of Langer’s result
in (106), except with $\lambda_{+}$ replaced by
$n^{T}Dn/|n^{T}H^{-\ddagger}n|$. In fact,
$n^{T}Dn/|n^{T}H^{-\ddagger}n|=\lambda_{+}$ when $n$ is parallel to $v_{+}$
(this should be clear from the fact that
$H^{-\ddagger}v_{+}=-\lambda_{+}^{-1}Dv_{+}$.) Furthermore $k(v_{+})$ is the
minimum of $k(n)$ (more strictly speaking, any vector proportional to $v_{+}$
will do). To see why, recall the earlier discussion related to the
simultaneous diagonalization of $H^{\ddagger}$ and $D$ by $V$ in (60). With a
change of variable $n=Vm$,

$\frac{n^{T}Dn}{|n^{T}H^{-\ddagger}n|}=\frac{m^{T}V^{T}DVm}{|m^{T}V^{T}H^{-\ddagger}Vm|}=\frac{m^{T}m}{|m^{T}\Lambda^{-1}m|}=\frac{\|m\|^{2}}{|-\lambda_{+}^{-1}m_{+}^{2}+\sum_{j\neq+}\lambda_{j}^{-1}m_{j}^{2}|}$
(137)

The expression $m^{T}\Lambda^{-1}m/m^{T}m$ is an example of a Rayleigh
quotient. In the absence of any constraint, the quotient is bounded between
the largest and smallest eigenvalues of $\Lambda^{-1}$. This bound is
unfortunately not helpful, since we require that $n^{T}H^{-\ddagger}n$ (and
thus $m^{T}\Lambda^{-1}m$) be negative. For any fixed length $\|m\|$, the
expression
$|-\lambda_{+}^{-1}m_{+}^{2}+\sum_{j\neq+}\lambda_{j}^{-1}m_{j}^{2}|$ in the
denominator is maximized under this constraint whenever the positive terms in
the sum $\sum_{j\neq+}\lambda_{j}^{-1}m_{j}^{2}$ are minimized. This is
achieved by setting all $m_{j}=0$ except for $|m_{+}|=\|m\|$; in other words,
$m$ is proportional to $e_{+}$, and thus $n$ is proportional to $v_{+}$ (the
change of variable described here is the same as that in Section 3.4.3, but
one should be careful not to confuse the vectors $n$ and $m$ defined in this
section with those defined in Section 3.4.3). Taken together, we have shown
that $n=v_{+}$ minimizes (136) and the minimum is equivalent to (106), as
desired.

##  5 Acknowledgement

I would like to thank Robert Alberstein for his critical reading of the note,
and I would like to acknowledge support through Chan Zuckerberg Biohub
Investigator funds to Dr. Tanja Kortemme (UCSF).

## References

  * [1] Langer JS (1969) Statistical theory of the decay of metastable states.  Annals of Physics 54(2):258–275.
  * [2] Kramers HA (1940) Brownian motion in a field of force and the diffusion model of chemical reactions.  Physica 7(4):284–304.
  * [3] Berezhkovskii A, Szabo A (2004) One-dimensional reaction coordinates for diffusive activated rate processes in many dimensions.  The Journal of Chemical Physics 122(1):014503.
  * [4] Hänggi P, Talkner P, Borkovec M (1990) Reaction-rate theory: Fifty years after Kramers.  Reviews of Modern Physics 62(2):251–341.
  * [5] Zhou HX (2010) Rate theories for biologists.  Quarterly Reviews of Biophysics 43(2):219–293.
  * [6] Peters B (2016) Reaction Coordinates and Mechanistic Hypothesis Tests.  Annual Review of Physical Chemistry 67(1):669–690.
  * [7] Peters B (2017) Reaction Rate Theory and Rare Events.  (Elsevier), 1st edition.
  * [8] E W, Vanden-Eijnden E (2010) Transition-Path Theory and Path-Finding Algorithms for the Study of Rare Events.  Annual Review of Physical Chemistry 61(1):391–420.
  * [9] Tuckerman ME (2010) Statistical Mechanics: Theory and Molecular Simulation.  (Oxford University Press, Oxford ; New York), Illustrated edition.
  * [10] Vanden-Eijnden E, Tal FA (2005) Transition state theory: Variational formulation, dynamical corrections, and error estimates.  The Journal of Chemical Physics 123(18):184103.
  * [11] Reimann P, Schmid GJ, Hänggi P (1999) Universal equivalence of mean first-passage time and Kramers rate.  Physical Review E 60(1):R1–R4.
  * [12] Metzner P, Schütte C, Vanden-Eijnden E (2006) Illustration of transition path theory on a collection of simple examples.  The Journal of Chemical Physics 125(8):084110.
  * [13] Dellago C, Bolhuis P, Geissler P (2006) Transition Path Sampling Methods in Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1, eds. Ferrario M, Ciccotti G, Binder K.  (Springer, Berlin, Heidelberg), pp. 349–391.
  * [14] Dellago C, Hummer G (2014) Computing Equilibrium Free Energies Using Non-Equilibrium Molecular Dynamics.  Entropy 16(1):41–61.
  * [15] Im W, Roux B (2002) Ions and Counterions in a Biological Channel: A Molecular Dynamics Simulation of OmpF Porin from Escherichia coli in an Explicit Membrane with 1M KCl Aqueous Salt Solution.  Journal of Molecular Biology 319(5):1177–1197.
  * [16] Liu P, Harder E, Berne BJ (2004) On the Calculation of Diffusion Coefficients in Confined Fluids and Interfaces with an Application to the Liquid-Vapor Interface of Water.  The Journal of Physical Chemistry B 108(21):6595–6602.
  * [17] Hummer G (2005) Position-dependent diffusion coefficients and free energies from Bayesian analysis of equilibrium and replica molecular dynamics simulations.  New Journal of Physics 7(1):34.
  * [18] Ma A, Nag A, Dinner AR (2006) Dynamic coupling between coordinates in a model for biomolecular isomerization.  The Journal of Chemical Physics 124(14):144911.
  * [19] Peters B (2009) Competing nucleation pathways in a mixture of oppositely charged colloids: Out-of-equilibrium nucleation revisited.  The Journal of Chemical Physics 131(24):244103.
  * [20] Berezhkovskii AM, Zitserman VY (1989) Anomalous regime for decay of the metastable state: An extension of multidimensional Kramers theory.  Chemical Physics Letters 158(5):369–374.
  * [21] Berezhkovskii AM, Szabo A, Greives N, Zhou HX (2014) Multidimensional reaction rate theory with anisotropic diffusion.  The Journal of Chemical Physics 141(20):204106.
  * [22] Northrup SH, McCammon JA (1983) Saddle-point avoidance in diffusional reactions.  The Journal of Chemical Physics 78(2):987–989.
  * [23] Klosek-Dygas MM, et al. (1989) Diffusion theory of multidimensional activated rate processes: The role of anisotropy.  The Journal of Chemical Physics 90(2):1141–1148.
  * [24] Chen W, Ferguson AL (2018) Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration.  Journal of Computational Chemistry 39(25):2079–2102.
  * [25] Ribeiro JML, Bravo P, Wang Y, Tiwary P (2018) Reweighted autoencoded variational Bayes for enhanced sampling (RAVE).  The Journal of Chemical Physics 149(7):072301.
  * [26] Grote RF, Hynes JT (1980) The stable states picture of chemical reactions. II. Rate constants for condensed and gas phase reaction models.  The Journal of Chemical Physics 73(6):2715–2732.
  * [27] Berezhkovskii AM, Pollak E, Zitserman VY (1992) Activated rate processes: Generalization of the Kramers–Grote–Hynes and Langer theories.  The Journal of Chemical Physics 97(4):2422–2437.
  * [28] Metzner P, Schütte C, Vanden-Eijnden E (2009) Transition Path Theory for Markov Jump Processes.  Multiscale Modeling & Simulation 7(3):1192–1219.
  * [29] Bello-Rivas JM, Elber R (2015) Exact milestoning.  The Journal of Chemical Physics 142(9):094102.
  * [30] Zuckerman DM, Chong LT (2017) Weighted Ensemble Simulation: Review of Methodology, Applications, and Software.  Annual Review of Biophysics 46(1):43–57.
  * [31] Dellago C, Bolhuis PG, Csajka FS, Chandler D (1998) Transition path sampling and the calculation of rate constants.  The Journal of Chemical Physics 108(5):1964–1977.
  * [32] van Erp TS, Bolhuis PG (2005) Elaborating transition interface sampling methods.  Journal of Computational Physics 205(1):157–181.