# Influence of low frequency modes on dynamical concertedness in double proton
transfer dynamics

Priyanka Pandey<EMAIL_ADDRESS>Department of Chemistry, Indian Institute
of Technology
Kanpur, Uttar Pradesh 208016, India  Shibabrat Naik<EMAIL_ADDRESS>School of Mathematics, University of Bristol
Fry Building, Woodland Road, Bristol BS8 1UG, United Kingdom  Srihari
Keshavamurthy<EMAIL_ADDRESS>Department of Chemistry, Indian Institute of
Technology
Kanpur, Uttar Pradesh 208016, India

###### Abstract

We analyze the classical phase space dynamics of a three degree of freedom
Hamiltonian that models multiple bond breaking and forming reactions. The
model Hamiltonian, inspired from studies on double proton transfer reactions,
allows for exploring the dynamical consequences of higher index saddles on
multidimensional potential energy surfaces. Studies have shown that coupling
of low frequency transverse modes to the reaction coordinate can significantly
influence the reaction mechanism, concerted or sequential, as inferred from a
reduced dimensional analysis. Using the notion of dynamically concerted and
sequential pathways, we provide insights into the role of the transverse modes
by studying the delay times between the formation of two bonds. The delay time
distribution, used extensively in earlier studies, is placed on a firm
dynamical footing by correlating it with the phase space manifolds, determined
using the technique of Lagrangian descriptors. We establish the utility of
Lagrangian descriptors in identifying the phase space manifolds responsible
for the dynamically concerted and dynamically sequential pathways.

###### keywords:

Double proton transfer , Sequential and concerted mechanisms , Delay time
distributions , Lagrangian descriptor , Phase space structures , Higher index
saddles

##  1 Introduction

The theory of nonlinear dynamical systems is a natural framework for
understanding chemical reactions. There are several reasons for such a claim,
but two of them are key. Firstly, breaking of a bond is only possible if the
vibrations are modeled as nonlinear oscillators. Secondly, the canonical
paradigm of associating an energized molecule with many such nonlinear
oscillators that are coupled together [1, 2, 3, 4] leads to a rich and complex
dynamical behaviour that necessitates a phase space perspective for proper
analysis and interpretation [5, 6]. Indeed the central notion of a transition
state is best understood as a dynamical bottleneck that is formed by certain
invariant manifolds in the phase space [7, 8, 9, 10, 11, 12, 13]. Reaction
rates can then be associated with fluxes through appropriate bottlenecks [14,
15, 16]. Thus, a local dynamical perspective on Transition state theory (TST)
has provided fresh insights into the usefulness and limitations of TST in the
microcanonical [17], canonical [18] and more general [19, 20, 21, 22, 23]
settings.

Apart from the rates, there is yet another important aspect of a chemical
reaction that is enshrined in the TS - the mechanism. In fact, identifying the
correct TS is essentially equivalent to a knowledge of the mechanism of the
reaction. For reactions involving a single TS (elementary reactions) one
therefore associates a single mechanism that leads to the transformation of
the reactants to products. However, several reactions are associated with
potential energy surface (PES) that exhibit novel features [24] like extended
flat regions (calderas or more generally entropic intermediates), ambimodal
TS, valley ridge inflection points, several distinct saddle points (multiple
TSs), and saddle points with more than one unstable direction (higher index
saddles). It is now clear that the existence of such features on the PES can
lead to significant dynamical effects. Examples include dynamical matching
[25, 26, 27, 28, 29, 30], nonstatistical branching ratios [31, 32, 33], energy
dependent product selectivity [34, 35, 36], and switching of reaction
mechanisms [37, 38, 39]. Consequently, there is an increased focus now on
trajectory-based analysis of complex reactions.

In the current work we are interested in understanding the dynamics of
reactions that involve breaking and forming of multiple bonds. Here one
invariably has to face up to a fundamental and essential mechanistic question:
is the process occurring sequentially or in a concerted fashion? In this
regard, the Diels-Alder and the double proton transfer (DPT) reactions have
provided a rich arena to explore the role of dynamics in determining the
correct reaction mechanism [40, 41, 42, 43, 44, 39, 45, 46, 47, 48]. From a
fundamental point of view the possibility of more than one distinct pathway is
linked with the presence of several distinct TSs. Although traditionally one
associates TS with a index-$1$ saddle point on the multidimensional PES,
several studies indicate that the dynamical influence of higher index saddles
on the PES can also be a deciding factor in identifying the dominant mechansim
[49, 50, 51, 52, 53]. For instance, for energies above the index-$2$ saddle
one can have a time-dependent switching between the concerted and sequential
pathways. Recently, it was shown [37] that such a dynamical mechanism switch
is an inherently classical phenomenon. Moreover, owing to the mixed regular-
chaotic nature of the classical phase space, initial quantum wavepackets that
are centered at specific regions of the classical phase space can undergo
strikingly different mechanism-switching dynamics [37]. Interestingly, the
switching timescale is typically of the order of a bond stretching time period
and hence ultrafast. These observations therefore raise questions on the
utility of a purely non-dynamical classification of the mechanism as concerted
or sequential. Such concerns have been raised by Carpenter in his early work
on the dynamic matching phenomenon wherein he emphasizes the “hazards
associated with partitioning of mechanisms into stepwise and concerted
categories” based purely on the features on the static potential energy
surface [26]. More recently, Houk and coworkers introduced quantitative
measures for classifying the mechanism as dynamically concerted or sequential.
Thus, for a given trajectory, if the time delay between the formation of the
first bond and the second is shorter than a specified timescale then that
particular trajectory is classified as dynamically concerted [43]. As a
consequence the central quantity of interest is the distribution of the delay
times associated with an appropriate ensemble of trajectories. Depending on
the nature of the delay time distributions one can identify the mechanism as
dynamically concerted or sequential.

Note that the approach of Houk and coworkers [43] implicitly invokes the
dynamics in the full classical phase space. Understandably, a detailed phase
space analysis of the ab initio molecular dynamics based studies of reactions
like the Diels-Alder is far from easy. At the same time, there is no denying
the fact that rationalizing the dynamics based on the phase space structures
is expected to yield rich dividends in terms of our ability to predict rather
than simply observe or compute. Thus, from a nonlinear dynamical systems point
of view it is natural to expect that the delay time distributions are
intimately linked to the disposition of the stable and unstable phase space
manifolds that lead to transport from the reactant to the product regions.
However, identifying, let alone computing, such manifolds in very high
dimensions is not feasible at the present moment. A crucial question then is
this: can reduced dimensional models capture enough of the essential dynamics
to allow for at least qualitative predictions? The answer, as apparent from
the several studies utilizing “minimal” models, is yes. For example,
significant dynamical insights into the phenomenon of roaming and dynamic
matching have come from phase space studies on the low dimensional model
systems [54, 55, 56, 57, 58, 25, 26, 27, 28, 29]. Nevertheless, the detailed
study of a electrocyclic ring opening reaction by Kramer et al. highlights the
central issues in this regard [59]. A comparison of the direct dynamics
calculations (in a $36$-dimensional phase space) with the reduced two-
dimensional model dynamics for the same reaction revealed that the two do
share dynamical similarities. However, they make an important point - in the
event that large amplitude modes, which would be considered as ”spectator”
modes in the reduced dimensional treatment, couple to the reaction coordinate,
the dynamics may be more complicated then what would be predicted by the
reduced dimensional models. Note that one can associate large amplitude modes
with low frequency vibrations and in a molecule with symmetry the various low
frequency modes can couple to the reactive mode in different ways due to the
symmetry constraints. Thus, apart from leading to a more complicated dynamics,
can the coupling of specific low frequency modes alter the inferred reduced
dimensional mechanism itself?

Since the present study focuses on the DPT reaction in a specific class of
molecules, we mention a few examples from earlier studies that highlight the
importance of the low frequency modes. In their extensive review of multiple
proton transfer dynamics, Smedarchina et al. [60] have argued for the
importance of the coupling of low frequency skeletal vibrations to the proton
transfer modes. Furthermore, in a path integral molecular dynamics simulation
Yoshikawa et al. have shown that the low frequency out-of-plane vibration can
suppress the concerted pathway in porphycene molecule [61, 62]. Another
example comes from the Car-Parrinello molecular dynamics study of porphycene
by Walewski et al., where it was observed that excitations of selective low
frequency modes, and combinations thereof, tend to enhance or suppress the
different mechanisms [63]. It is also relevant to point out the classical ab
initio molecular dynamics study of DPT by Ushiyama and Takatsuka where, apart
from hints to the importance of delay time distributions, the crucial role of
skeletal vibrations to the second proton transfer was emphasized [45].

Clearly, and as discussed in detail in sec. 2, for a DPT reaction with two
reactive modes, coupling of even one low frequency transverse mode results in
a system with three degrees of freedom. In this work we investigate the
classical dynamics of such a model system with the aim of explicitly
correlating the delay time distributions with the appropriate phase space
manifolds. In particular, as mentioned above, we investigate the influence on
the delay time distributions due to the coupling of a third mode belonging to
a specific symmetry class. The results show that while high frequency modes do
not significantly change the fraction of concerted trajectories, the low
frequency modes can substantially reduce the fraction. An explanation of our
delay time results in terms of the phase space structures is given by
computing the appropriate manifolds using the technique of Lagrangian
descriptors(LD). We show that the the LD maps faithfully capture the changes
in the delay time distributions with varying coupling strength of the third
mode. In sec. 2 and A we motivate the model Hamiltonian used in our study. The
influence of the third mode on the delay time distributions are presented in
sec. 3.3, followed by the results of the delay time distributions for varying
frequencies and coupling strengths. In sec. 3.4, aided by the stability
analysis of the linearized flow in B, the dynamical trajectory observations
are correlated with the LD-based determination of the relevant phase space
manifolds. Finally, sec. 4 concludes with a brief summary and future outlook.

##  2 Model Hamiltonian

Figure 1: (a) The porphycene molecule with $N=38$ atoms. The nitrogen, carbon,
and hydrogen atoms are indicated in blue, grey, and white colors respectively.
(b) Schematic for the double proton transfer model (indicated by arrows in
both the panels) corresponding to the dashed square region shown in (a). The
coordinates $(x_{1},x_{2})$ correspond to the two proton transfer events. Each
proton transfer subsystem is described by an appropriate double well potential
with minima at $\pm\Delta x_{0}$. Interaction between the two subsystems is
mediated by the coupling constant $G$ and $D$ at the leading order. See the
main text and A for details.

To study the double proton transfer reaction, we consider a model three degree
of freedom Hamiltonian motivated by the models introduced by Smedarchina et
al. [64, 65] in their extensive studies. In Fig. 1 we show a schematic for the
DPT process. The two proton transfer events occurring in the molecule (labeled
as subsystem $1$ and $2$ in the figure) are described by one dimensional
coordinates $x_{1}$ and $x_{2}$ with associated masses $m_{1}$ and $m_{2}$,
which are taken to be equal to the proton mass $m_{H}$. As shown in detail in
A, an appropriate model two degrees of freedom dimensionless Hamiltonian for
the coupled proton transfer is conveniently expressed in terms of the
coordinates
$(X_{s},X_{a})\equiv(\sqrt{M}x_{s},\sqrt{M}x_{a})=(\sqrt{M}(x_{1}+x_{2})/2,\sqrt{M}(x_{1}-x_{2})/2)$
with $M=m_{1}+m_{2}=2m_{H}$. The Hamiltonian is of the form

$H({\bf X},{\bf P})=\frac{1}{2}(P_{s}^{2}+P_{a}^{2})+U(X_{s},X_{a})$ (1)

with $(P_{s},P_{a})$ being the momenta conjugate to $(X_{s},X_{a})$. The two
dimensional potential energy surface is given by

$U({\bf X})=\bar{\alpha}_{s}\left[X_{s}^{2}-(\Delta
X_{s})^{2}\right]^{2}+\bar{\alpha}_{a}\left[X_{a}^{2}-(\Delta
X_{a})^{2}\right]^{2}+2\bar{R}X_{s}^{2}X_{a}^{2}+{\cal U}(G,D)$ (2)

In the above we have denoted $({\bf X},{\bf
P})\equiv(X_{s},X_{a},P_{s},P_{a})$ with the parameters
$\bar{\alpha}_{s}=\bar{\alpha}_{a}\equiv(1-D)/M^{2}$ and
$\bar{R}\equiv(3+D)/M^{2}$. The various minima on the PES are given in terms
of the quantities

$\Delta X_{s,a}=\sqrt{\frac{M(1\pm G)}{1-D}}$ (3)

and the constant energy shift is denoted as

${\cal U}(G,D)=1-\bar{\alpha}_{s}(\Delta X_{s})^{4}-\bar{\alpha}_{a}(\Delta
X_{a})^{4}$ (4)

The parameters $G$ and $D$ are specific to a given system (molecule) and
correspond to the coupling of the two proton transfer coordinates. We refer
the reader to the A for a detailed derivation of the above Hamiltonian along
with the relevant mass, length, and time scales. In Fig. 2 the two dimensional
PES are shown in the two different sets of coordinates. Note that, in general,
the number and type of critical points on the PES [65] depend on the values of
$G$ and $D$, mimicking a wide variety of dynamical systems. For the values of
interest to us in the current work the PES exhibits a total of nine critical
points which, as seen in Fig. 2, include four minima, four index-$1$ saddles
and one index-$2$ saddle.

Figure 2: (Left panel) A contour plot of the model potential energy surface in
the local ($x_{1},x_{2}$) coordinates. (Right panel) A three dimensional plot
of the potential energy surface in the mass-weighted ($X_{s},X_{a}$)
coordinates. The parameter values are $G=0.2$, and $D=0.15$. Different regions
of the potential energy surface are shown inside the dotted space of the
contour plot where $R$, $P$, and $I_{1,2}$ represent the reactant, product and
intermediate regions respectively. The central box represents the index-$2$
saddle region. Note that the PES also exhibits four index-$1$ saddles, which
are shown as $\times$.

In the molecular context the Hamiltonian in Eqn. 1 captures the dynamics
corresponding to the key reactive degrees of freedom. Thus, for an $N$-atom
molecule of interest with $(3N-6)$ vibrational degrees of freedom Eqn. 1
accounts for two of the degrees of freedom. However, the remaining $(3N-8)$
modes that are transverse to the reactive modes typically do couple to the
${\bf X}$ degrees of freedom to varying extents. Moreover, if the molecule of
interest has a certain point group symmetry then the various transverse modes
are constrained to couple to ${\bf X}$ with specific functional form of the
coupling potentials. For instance, in the context of DPT all the $(3N-8)$
modes denoted by ${\bf Y}$ couple at leading order via the potential

$U_{\rm coup}({\bf
X},Y_{k})=\frac{1}{2}\omega_{kY}^{2}\left[Y_{k}-\frac{\lambda}{\omega_{kY}^{2}}g({\bf
X})\right]^{2}$ (5)

with $g({\bf X})=X_{s,a},X_{s}X_{a}$, and $X_{s,a}^{2}$ depending on the
symmetry class to which the $Y_{k}$-mode belongs. In Eqn. 5 the $Y_{k}$ mode
is modeled as a harmonic oscillator with the frequency of the mode denoted by
$\omega_{kY}$ and $\lambda$ being a measure of the coupling strength. It is
interesting to note that with the above form of coupling one can still think
of the ${\bf Y}$-modes as providing a “bath”, albeit a structured one. This is
in contrast to the usual system-bath models wherein all the “bath” ${\bf
Y}$-modes would couple bilinearly with a specified spectral density. Clearly,
the dynamical implications of a structured and a non-structured bath are
expected to be quite different for the reaction process.

As mentioned in the introduction, several studies have indicated the
importance of including the additional modes since they can have a significant
effect on the mechanism inferred from an analysis of the reduced dimensional
Hamiltonian. Thus, certain symmetry modes tend to enhance a specific mechanism
(concerted or sequential) whereas certain other symmetry modes act in an
opposite manner [63]. Therefore, in order to rationalize the observed DPT
rates one minimally needs to include two of the ${\bf Y}$-modes, resulting in
a four degree of freedom Hamiltonian. However, understanding global phase
space transport and linking it to the dynamical influence of the coupled ${\bf
Y}$-modes in this case is a challenging task. Instead, here we address a
simpler yet nontrivial question - can the coupling of a ${\bf Y}$-mode with a
given symmetry significantly influence the reaction mechanism as inferred from
the low dimensional system in Eqn. 1? And, if so, what is the dynamical origin
of such a modulation? As noted above in Eqn. 5, there are several choices for
the model Hamiltonian according to the $g({\bf X})$ of interest. In this work
we focus on the specific three degree of freedom Hamiltonian

$\displaystyle H({\bf X},{\bf P},Y,P_{Y})$ $\displaystyle=$
$\displaystyle\frac{1}{2}(P_{s}^{2}+P_{a}^{2}+P_{Y}^{2})+U({\bf
X})+\frac{1}{2}\omega_{Y}^{2}\left[Y-\frac{\lambda}{\omega_{Y}^{2}}X_{s}^{2}-\frac{\lambda^{\prime}}{\omega_{Y}^{2}}X_{a}^{2}\right]^{2}$
(6) $\displaystyle\equiv$
$\displaystyle\frac{1}{2}(P_{s}^{2}+P_{a}^{2}+P_{Y}^{2})+V({\bf X},Y)$ (7)

with $U({\bf X})$ being the potential in Eqn. 2 and we continue to adopt the
mass-weighted coordinate representation. Note that the above form of coupling
corresponds to the so called $a_{g}$-symmetry $Y$-mode and the importance of
this coupling to DPT has been noted in several earlier studies. For example,
in the $N=38$ atom porphycene molecule one has a total of $3N-8=106$ ${\bf
Y}$-modes. Among these modes the ones with $a_{g}$ symmetry have substantial
projection [66] onto the two reactive modes ${\bf X}$.

Figure 3: Potential energy visualized as equipotential surfaces for three
energies. The green surface denotes $E<E_{s}$, red surface denotes $E=E_{s}$,
and the blue surface denotes $E>E_{s}$. The energy of the index-2 saddle is
$E_{s}$ and denoted by the $+$, while the index-1 saddles are shown as
$\times$. The parameter values are $G=0.2$, $D=0.15$, $\omega_{Y}=0.3$,
$\lambda=0.3$ and $\lambda^{\prime}=0.1$ for the total energy, $E=1.1$ to
visualize the equipotential surfaces.

For the purpose of the current study we choose $0<G<1/2$ and $|D|<2G$ which
yields a total of nine critical points. The details associated with the
critical points are given in Table 1. Note that there is no restriction on the
sign of $D$. At the same time, for a given $G$, the dynamics corresponding to
positive or negative $D$ can be sufficiently different. For the rest of the
paper we fix the values $G=0.2$ and $D=0.15$. Moreover, we fix the total
energy at $E=1.1$, which is slightly above the index-$2$ saddle energy (cf.
Table 1). Consequently, both the concerted and sequential pathways (examples
can be seen in Fig. 4) from reactant to product are available classically. A
key objective of the current study is to relate the phase space dynamics of
the Hamiltonian in Eqn. 7 with the mechanism of DPT. In particular, we intend
to assess the influence of the transverse $Y$-mode with both low and high
frequencies $\omega_{Y}$ over a range of the couplings
$(\lambda,\lambda^{\prime})$. In Fig. 3 a representation of the PES
$V(\mathbf{X},Y)$ is shown for three values of the energy,
$E<E_{s},E=E_{s},E>E_{s}$, where $E_{s}$ is the energy of the index-2 saddle
at the origin.

Configuration space coordinates | Total energy | Linear stability | Description
---|---|---|---
$\left(\pm\Delta X_{s},0,\dfrac{\lambda}{\omega_{Y}^{2}}(\Delta X_{s})^{2}\right)$ | $1-\bar{\alpha}_{s}(\Delta X_{s})^{4}$ | C-C-C | reactant and product $R,P$
$\left(0,\pm\Delta X_{a},\dfrac{\lambda^{\prime}}{\omega_{Y}^{2}}(\Delta X_{a})^{2}\right)$ | $1-\bar{\alpha}_{a}(\Delta X_{a})^{4}$ | C-C-C | intermediates $I_{1,2}$
$\left(X_{s}^{\ddagger},X_{a}^{\ddagger},Y^{\ddagger}\right)$ | $H^{\ddagger}$ | S-C-C | index-$1$ saddles
$(0,0,0)$ | $E_{s}=1$ | S-S-C | index-$2$ saddle
Table 1: Equilibrium points, their respective energies and linear phase space
stability. The stability types are denoted by $C$ for center and $S$ for
saddle. For an explicit expression for
$\left(X_{s}^{\ddagger},X_{a}^{\ddagger},Y^{\ddagger}\right)$ and the
associated energy $H^{\ddagger}$ see B.

##  3 Results and Discussions

###  3.1 Computational preliminaries: defining the initial ensemble, relevant
regions, and delay time

In order to understand the influence of the $(\lambda,\lambda^{\prime})$
couplings and the frequency $\omega_{Y}$ of the transverse $Y$-mode on the
dynamics of the Hamiltonian in Eqn. 7, we compute the delay time distribution.
The concept of delay time distribution is motivated by the dynamical studies
of Diels–Alder reactions by Houk and coworkers [43]. As the name suggests, the
delay time corresponds to the time difference between the transfer of the
first proton and the subsequent transfer of the second proton. By definition,
the delay time is zero for a pure concerted pathway that proceeds directly
from the reactant $R$ to the product $P$ via the index-$2$ saddle, and without
visiting the intermediate regions. On the other hand, a sequential pathway
from $R$ to $P$ via any or both the intermediate regions $I_{1,2}$ yields a
finite value for the delay time. Therefore, demarcating different regions in
the configuration space (i.e., reactant, intermediate, product, and index-2
saddle) is necessary in order to compute the delay time and related measures.
Consequently, in Fig. 2, the definitions of the various regions of the PES
used in this study are shown as boxes of various sizes. We have divided the
PES into five regions. The product and the reactant regions are centered at
the minima ($\pm\Delta X_{s},0$), and the region boundaries are at ($\pm
0.1\Delta{X_{s}},\pm 0.15$) from the minima. Similar regions are defined
around the intermediate minima ($0,\pm\Delta{X_{a}}$). However, for the
intermediate regions, the region boundaries are at ($\pm 0.40,\pm
0.33\Delta{X_{a}}$) from the minima. Since the height and the spread of the
wells containing the global and local minima vary, therefore we are using
different box sizes. The index-$2$ saddle region is centered at (0,0) with the
boundaries at ($\pm 0.05,\pm 0.05$). Note that there are other ways to assign
the different regions on the PES and, naturally, the quantitative delay time
distributions will be sensitive to the specific choice. However, within
reasonable definitions of the regions, at a qualitative level the results are
not expected to be significantly different i.e., the key trends with varying
parameters are preserved.

Figure 4: Examples trajectories with different delay times projected onto the
$(X_{s},X_{a})$ space. The various region definitions are shown as dashed
boxes. (Left column) Region definitions used in the present work. (Middle
column) Increased $R$ and $P$ region sizes. (Right column) Decreased $I_{1}$
and $I_{2}$ region sizes. All trajectories start at the index-$2$ saddle
(solid circle) with a total energy $E=1.1$. The forward and backward time
propagation from the index-$2$ saddle are shown in grey and black
respectively. The delay time assigned in each case is indicated in the
respective panels. Note that the four index-$1$ saddles are indicated by
$\times$. The example with zero delay time is a “pure” concerted trajectory.

For our calculations, we choose the initial values of position coordinates
$X^{(0)}_{s}$ and $X^{(0)}_{a}$ randomly from the index-$2$ saddle region.
This choice of the initial ensemble is to focus on the influence of the
index-$2$ saddle on the DPT. In addition, we randomly choose the initial
momentum $P^{(0)}_{s}>0$ and fix the initial value of the third mode
coordinate $Y^{(0)}$ and its conjugate momentum $P^{(0)}_{Y}$ at $(0,0)$.
Finally, the initial momentum $P^{(0)}_{a}$ is obtained by the energy
conservation condition i.e., $H({\bf X}^{(0)},{\bf
P}^{(0)},Y^{(0)}=0,P^{(0)}_{Y}=0)=E=1.1$. The specific total energy value,
fixed for the rest of the study, corresponds to being just above the index-$2$
saddle energy. Furthermore, note that the choice $P^{(0)}_{s}>0$ corresponds
to the trajectories at the index-$2$ saddle having momentum in the direction
of the product $P$. It is well known that the product selectivity of a
chemical reaction in a trajectory calculation is strongly associated with the
momentum distribution at the TS [67, 25]. Thus, although different initial
momentum distributions at the index-$2$ saddle can lead to quantitatively
different results, we believe that the qualitative insights are fairly robust.
A total of $10^{4}$ trajectories were initiated from the index-$2$ saddle
region and propagated both in the forward and backward direction until they
reach the product and the reactant regions respectively. Trajectories,
propagated up to a final time $t_{f}=300$, are deemed to be reactive if they
form the product in the forward direction and the reactant in the backward
direction. Depending on the path a reactive trajectory takes and the
associated delay time, we can characterize them as concerted or sequential
trajectories. In Fig. 4 we show examples of a concerted and several sequential
trajectories.

For the delay time computation we adopt the following strategy. The first
instance when the forward time trajectory enters either of the intermediate
regions $I_{1,2}$ is noted as $t_{I}$. This event corresponds to the transfer
of either one of the proton (cf. Fig. 1). Subsequently, the time at which the
trajectory enters the defined product region $P$ is denoted as $t_{P}$. This
time corresponds to the second proton transfer and hence
$\Delta\tau=t_{P}-t_{I}$ is associated as the delay time for the specific
trajectory. Two points are important to note at this stage. First, between
$t_{I}$ and $t_{P}$ the trajectory may visit the reactant region $R$ or visit
the intermediate regions several times. Second, the $\Delta\tau$ as defined is
sensitive to the extent of the different regions. Thus, variations in the
region sizes can change the $\Delta\tau$ for a given trajectory. Examples for
the same are shown in Fig. 4 and it is clear that some of the large
$\Delta\tau$ can become considerably smaller or a short $\Delta\tau$
sequential trajectory can turn into a concerted trajectory. However, the
effect of such variations on the distribution $P(\Delta\tau)$ shown in Fig. 6
is not expected to be significant. In particular, the qualitative trends seen
in Fig. 5 and Fig. 6 are robust to small variations in the region sizes.

Figure 5: The fraction of concerted trajectories
($f_{conc}=N_{conc}/N_{reac}$) as a function of the coupling parameters. Here
$N_{reac}$ is the total number of reactive trajectories i.e., ones that start
at the index-$2$ saddle and go to the product P and reactant R regions in
forward and backward time propagation. The total integration time is
$t_{f}=300$. (a) Variation with $\lambda$ for fixed $\lambda^{\prime}=0$ (b)
Variation with $\lambda^{\prime}$ for fixed $\lambda=0$ (c) Variation with
$\lambda^{\prime}$ for fixed $\lambda=0.3$. The frequency $\omega_{Y}$ of the
transverse $Y$-mode are shown in the legend. It is important to note that
$N_{reac}$ varies with the $(\lambda,\lambda^{\prime},\omega_{Y})$ parameters.

###  3.2 Effect of the third degree of freedom on the concerted pathways

Before discussing our results for the delay time distributions, in Fig. 5 we
show the influence of the third degree of freedom coupling in Eqn. 7 on the
“pure” (as opposed to dynamically) concerted mechanism. In particular, Fig.
5(a) and (b) show the fraction of concerted trajectories $f_{conc}$ upon
coupling only the $X_{s}-Y$ modes ($\lambda^{\prime}=0)$ and the $X_{a}-Y$
modes ($\lambda=0$) respectively. Note that the Hamiltonian in Eqn. 7 involves
both the couplings and hence the results in Fig. 5(a) and (b) are a bit
artificial. Nevertheless, such an analysis allows for dissecting, and a better
understanding, of the results for the actual $a_{g}$-symmetry coupling form.
It is clear from Fig. 5(a) that, apart from the initial oscillatory
nature111We remark here that the $f_{conc}$ in Fig. 5(a) exhibits peaks at
certain values of $\lambda\equiv\lambda_{p}$. Interestingly, these peaks seem
to occur when the reactant well harmonic frequencies $\omega_{s}(\lambda)$ and
$\omega_{a}$ become degenerate. An approximate estimate is
$\lambda_{p}\approx\omega_{Y}\left[(\omega_{a}/2\Delta
X_{s})^{2}-2\bar{\alpha}_{s}\right]^{1/2}$. For the parameters of interest,
$\lambda_{p}\approx 0.47\omega_{Y}$. Note that the degeneracy is driven by the
$X_{s}-Y$ coupling and hence an explcit three degrees of freedom effect. At
the present moment we do not have a dynamical insight into this observation.,
low $\omega_{Y}$ tend to drastically reduce $f_{conc}$ for increasing
$\lambda$, while relatively larger $\omega_{Y}$ lead to a slight increase. On
the other hand, the results in Fig. 5(b) indicate that $f_{conc}$ increases
moderately upon increasing the $X_{a}-Y$ mode coupling strengths. Therefore,
the full $a_{g}$-symmetry coupling case with $\lambda,\lambda^{\prime}\neq 0$
should encode the subtle competition between the two couplings. This is
confirmed in Fig. 5(c) where, as an example, the variation in $f_{conc}$ with
$\lambda^{\prime}$ for a fixed value of $\lambda=0.3$ is shown. Interestingly,
now the high $\omega_{Y}$ cases show very little variation over a significant
range of the $X_{a}-Y$ coupling strengths. In contrast, for $\omega_{Y}\leq
0.4$ the results are more complex with $f_{conc}$ increasing with
$\lambda^{\prime}$ and the oscillations seen in Fig. 5(a) being absent. As
expected, for $\lambda^{\prime}>\lambda=0.3$ one observes $f_{conc}$
increasing substantially. Nevertheless, it is evident from Fig. 5(c) that even
for relatively large $\lambda^{\prime}$ values the $\omega_{Y}=0.2$ and $0.3$
cases have considerably lower $f_{conc}$ in comparison to the uncoupled case.
We remark that these results agree with the general expectation that coupling
of the large amplitude (low frequency) modes to the reaction coordinate can
lead to dynamical behaviours that are vastly different from the dynamics of
reduced dimensional systems.

Note that Fig. 5 pertains to the pure concerted pathways and hence, by
definition, zero delay times. Based on the discussions in the introduction, a
useful perspective is to focus on the fraction of dynamically concerted
trajectories. Thus, although Fig. 5 indicates that low values of $\omega_{Y}$
lead to a reduced $f_{conc}$, is it possible that most of the trajectories are
still dynamically concerted for a reasonable choice for the delay time cutoff
$\Delta\tau_{c}$. In other words, if the distribution of delay times
$P(\Delta\tau)$ associated with the initial ensemble of trajectories in Fig. 5
is strongly peaked for $\Delta\tau\leq\Delta\tau_{c}$ then the mechanism would
be labeled as dynamically concerted. Consequently, as argued by Black et al.
[43], the significant lowering of $f_{conc}$ for small $\omega_{Y}$ values
observed in Fig. 5 need not really imply a major change in the reaction
mechanism. Therefore, to ascertain if this indeed is the case we now turn our
attention to the computation of the delay time distributions.

###  3.3 Delay time distributions: importance of the low frequency transverse
modes

From the discussions above, it is clear that in order to analyze
$P(\Delta\tau)$ results for our model system it is essential to define the
cutoff $\Delta\tau_{c}$. One possible choice for this cutoff time is related
to the lifetime of a TS according to the Eyring equation [68, 69]. This
timescale is set by the prefactor of TST i.e., $\Delta\tau_{c}\sim h/k_{B}T$
with $h,k_{B}$, and $T$ being the Planck constant, Boltzmann constant and
temperature respectively. However, as we are dealing with a index-$2$ saddle
and the dynamics is at zero temperature, we choose $\Delta\tau_{c}$ based on
the timescales associated with the reactant oscillations [43] or the unstable
motion at the index-$2$ saddle. Such criteria have been invoked before in
several studies [70]. For our model system and parameters of interest, as
shown in B, the unstable frequencies $\Omega_{s}^{*}\sim 1.6$ and
$\Omega_{a}^{*}\sim 1.3$ at the index-$2$ saddle. These frequencies,
independent of the couplings $(\lambda,\lambda^{\prime})$ due to the form of
the Hamiltonian, translate roughly to a timescale $T^{*}_{s,a}\sim 5$. On the
other hand, of the two harmonic frequencies around the reactant minimum, only
$\Omega_{s}$ depends on $(\lambda,\omega_{Y})$ and varies from $\sim 5.5$
($\omega_{Y}=0.2$) to $\sim 3.0$ ($\omega_{Y}=0.5$), while $\Omega_{a}\sim
2.7$ stays fixed. Consequently, the harmonic timescales associated with the
proton transfer modes at the reactant minimum is about $T_{h}\sim 2$. In this
work, we therefore choose the conservative estimate $\Delta\tau_{c}\sim 5$ for
discussing the delay time results.

Figure 6: Delay time distributions $P(\Delta\tau)$ for fixed $\lambda=0.3$ and
varying $\lambda^{\prime}$. The frequency $\omega_{Y}$ of the transverse
Y-mode is indicated in each case. The histograms in blue, red, and green
correspond to $\lambda^{\prime}$ value $0.1$, $0.3$ and $0.5$ respectively.
The orange vertical bar at $\Delta\tau_{c}\sim 5$ indicates the timescale
associated with the unstable directions at the index-$2$ saddle. Note that the
axis scales are different for each case and for comparison in each case the
vertical dashed line corresponds to the harmonic period $2\pi/\omega_{Y}$
associated with the transverse mode.

The results of the delay time computations are shown in Fig. 6 for fixed
$\lambda=0.3$ as histograms222Note that as discussed previously and shown in
Fig. 4, changing the region sizes will lead to some reshuffling of the counts,
particularly for those with very large delay times. Nevertheless, the small to
moderate time counts and their observed shifts should be robust.. In each
panel of Fig. 6 the transverse mode frequency $\omega_{Y}$ is fixed and the
delay time distributions for three values of $\lambda^{\prime}$ are shown.
Note that the both $f_{conc}$, shown in Fig. 5, and $P(\Delta\tau)$ are
computed using same initial ensemble. In addition, the parameters used for
generating Fig. 6 are fairly representative of other parameter sets as well.
For the value of $\omega_{Y}=0.4$ and $0.5$ Fig. 6 (bottom panels) shows that
a large fraction of the distribution is concentrated for
$\Delta\tau\leq\Delta\tau_{c}$, implying that the mechanism is dynamically
concerted. In contrast, for the case of $\omega_{Y}=0.2$ it is clear that the
mechanism is sequential for $\lambda^{\prime}=0.1$ and $0.3$ with the
emergence of dynamically concerted behaviour for larger coupling strengths.
However, the fact that there are substantial peaks for $\Delta\tau\gg 10$ does
hint at a fairly complex reaction dynamics. Clearly, the most complex
variations in the distribution are seen in Fig. 6 for the $\omega_{Y}=0.3$
case. Here, despite the general trend of the onset of dynamical concerted
behaviour with increasing $\lambda^{\prime}$, even for the largest coupling a
substantial fraction of the trajectories exhibit dynamically sequential
mechanism. Given the opposing trends in $f_{conc}$ observed in Fig. 5(a) and
(b), one perhaps anticipates the $\lambda\sim\lambda^{\prime}$ case for lower
values of $\omega_{Y}$ to be in a sort of “crossover” region.

It is worthwhile pointing out the following interpretation of the delay time
distribution results presented here. In a given molecular system,
characterized by the parameters $G$ and $D$, the multitude of $a_{g}$-symmetry
modes couple with a range of $\omega_{Y},\lambda$, and $\lambda^{\prime}$
values. The results in Fig. 6 then suggest that in the full multidimensional
system whether the mechanism is dynamically concerted or sequential depends
rather sensitively on the set of ratios
$[(\lambda/\lambda^{\prime})_{1},(\lambda/\lambda^{\prime})_{2},\ldots,(\lambda/\lambda^{\prime})_{n_{ag}}]$,
where $n_{ag}$ is the total number of $a_{g}$-symmetry modes in a specific
molecule. Clearly, similar criteria should exist for other transverse modes
belonging to different symmetry classes. At the moment there is not much known
about the dynamical competition between two or more low frequency modes with
different symmetries. Nevertheless, Fig. 6 does provide a clue as to why any a
priori decision on the mechanism based solely on the static PES features is
bound to be problematic. To this end, in the following section we provide
further support by establishing a link between the phase space manifolds and
the delay time distributions.

###  3.4 Phase space viewpoint: Lagrangian descriptors are correlated with
delay time distributions

A crucial observation, as shown in B, is that the linear analysis of the
index-2 saddle equilibrium does not shed any light on the changes in the
fraction of concerted trajectories with coupling strengths shown in Fig. 5.
The eigenvalues of the linearized system at the index-2 saddle is independent
of the coupling strengths, $\lambda,\lambda^{\prime}$. This implies that the
competition of concerted vs sequential pathways for initial conditions
launched from the vicinity of the index-2 saddle is inherently mediated by the
global phase space structures. More so, these are global invariant manifolds
in the phase space and transport initial conditions between intermediate and
product wells. Thus, differentiating which initial conditions have low delay
time, that is dynamically concerted, and high delay time, that is dynamically
sequential.

Figure 7: (Top panel) Forward LD map for $\tau=10$ and (bottom panel) delay
time map in the ($P_{s},X_{s})$ space at total energy of $E=1.1$. The
parameter values are $\lambda=0.3$, $\lambda^{\prime}=0.1$. The frequency
$\omega_{Y}$ of the transverse Y-mode is indicated above each column. The
color scales associated with LD and delay time map are indicated in the
respective panels. Note that the initial conditions of “pure” concerted and
dynamically concerted trajectories are shown in blue. The green color
indicates the initial condition of trajectories with delay time $10$ or
greater in the delay time map. The empty (white color) spaces in the delay
time map correspond to initial conditions that are non-reactive up to the
final time of integration.

In this study, we use Lagrangian descriptors [71, 72, 73] (see C for details
on the method) to identify the changes in the phase space structures with the
changes in the coupling strengths and frequency of the third mode. In the case
of three or more degrees of freedom systems, this method has been used to
detect invariant manifolds and reactive islands [74, 75, 76], discovering
structure in the nuclear phase space in nonadiabatic quantum dynamics [77],
while there is an increasing number of analysis for one and two degrees of
freedom system with and without dissipation and time dependence. We refer the
reader to the references in the open-source book on Lagrangian descriptors
[73]. However, for three degrees of freedom systems with multiple saddles with
varying indices, the use of LD has not been studied carefully and we present
some preliminary discussion of this method.

For the three degrees of freedom system, we define the two dimensional section
on the five dimensional energy surface

$\Sigma_{X_{s}P_{s}}^{+}=\left\\{(X_{s},X_{a},Y,P_{s},P_{a},P_{y})\in\mathbb{R}^{6}\,|\,X_{a}=0,Y=0,P_{y}=0,\dot{X}_{a}>0\right\\}$
(8)

to inspect changes in the phase space structures with changes in the coupling
strength and frequency of the third mode. We compute the Lagrangian descriptor
for initial conditions on the reactant side ($X_{s}\leq 0$) of the section
(Eqn. (8)). The integration time used if $10$ time units which is almost
double the cut-off time for dynamical concerted ($\Delta\tau_{c}\sim 5$)
pathway. While most studies using the LD method, and supported by theoretical
arguments, tend to choose high integration time, we found that the short
integration time of $10$ time units gave sufficient time for the structures to
form and did not generate the many intricate stretching and folding of the
global invariant manifolds. We compare the LD contour map with the delay time
map in Figs. 7-9 to show, for the first time, a striking correspondence
between the invariant manifolds and delay time distribution. First, we observe
that a direct correspondence in the contours of delay times and LD values
across all the coupling strengths and frequency of the third mode with a cut-
off time for dynamically concerted behaviour of $\Delta\tau_{c}\sim 5$ time
units. The invariant manifolds identified in the LD contour map correspond to
initial conditions with high delay time. However, the regions bounded by the
invariant manifolds have two distinct delay times, that is either below
$\Delta\tau<5$ or $\Delta\tau>10$. It implies that the invariant manifolds
partition the initial conditions into dynamically concerted or sequential
mechanisms. However, it is unclear as to which invariant manifolds can be
unambiguously tied to a given mechanism. In order to discern which invariant
manifolds mediate dynamically concerted and dynamically sequential mechanisms
for energies above the index-2 saddle, one needs to evolve an ensemble of
trajectories inside regions bounded by the invariant manifolds. This needs to
be paired with a computation of the normally hyperbolic invariant manifolds
[78] and it’s associated global invariant manifolds. The geometry of the
normally hyperbolic invariant manifolds (3-sphere) associated with the index-2
saddle and its stable and unstable invariant manifolds (spherical cylinders or
with geometry $\mathbb{S}^{2}\times\mathbb{R}$) is still an area of continued
interest [79, 80, 37] and we expect their structure and stability in the
parameter space will shed light on the precise phase space mechanism of the
competition between the concerted and sequential pathways.

Figure 8: Same as in Fig. 7 with parameter values $\lambda=0.3$ and
$\lambda^{\prime}=0.3$. Figure 9: Same as in Fig. 7 with parameter values
$\lambda=0.3$ and $\lambda^{\prime}=0.5$.

##  4 Conclusion and Outlook

In this work we have studied the classical dynamics of a three degrees of
freedom Hamiltonian which models the double proton transfer reaction in a
particular class of molecules. However, the analysis and techniques presented
here are expected to be relevant for other types of systems which involve
breaking and forming of multiple bonds. The key points that emerge from our
study are as follows:

  1. 1.

Coupling of additional low frequency modes to the reactive modes can lead to a
change in the reaction mechanism inferred from lower dimensional studies. In
particular, in the context of the DPT reaction studied in this work, it is
seen that even a single low frequency mode can substantially change the
fraction of reactive trajectories that proceed along the concerted pathway. It
would be instructive to construct the phase space dividing surface for the
index-$2$ saddle, along the lines of the earlier work by Collins, Ezra, and
Wiggins [80], to gain further insights into the modulation of the fraction of
concerted trajectories.

  2. 2.

Inspired by several earlier studies on various reactions that involve multiple
bond formation, we have explored the utility of classifying reactions as
dynamically concerted or sequential. To this end a simple, but dynamical,
measure involving the time delay between the formation of two bonds was used.
We suggest that the distributions of delay times provides much more
information then the fraction of concerted trajectories and its use is
motivated by the work of [43]. More importantly, we show a direct
correspondence between the delay times and the phase space invariant manifolds
for the coupling strengths and frequency of the third mode. This observation,
therefore, places the earlier studies on a firm dynamical basis.

  3. 3.

We have shown that the technique of Lagrangian descriptors can be invoked to
map out the relevant invariant manifolds in high dimensional phase space. In
particular, although not explored further in the current work, we observe that
the LDs do encode the manifolds responsible for both the dynamically concerted
and the dynamically sequential reactive pathways. Further studies on the
changes in the structure of the index-$2$ normally hyperbolic invariant
manifolds (NHIM) [8, 81] with the transverse mode frequency and connection to
the LD maps will be the focus of our future work.

Several issues arise in the context of our, admittedly preliminary, study and
we briefly mention a few. Firstly, are the concerted and sequential pathways
uncorrelated? One way to address this is to compute the so called gap time
distribution [82, 83, 17, 84, 85] for the model system and the possible
connections to the delay time distributions. Such a connection, along with an
unambiguous disentangling of the phase space invariant manifolds for the two
mechanisms, will then allow for decomposing the rate of the reaction in terms
of “concerted rates” and “sequential rates”. Secondly, the extent of
intramolecular vibrational energy redistribution (IVR) [4, 5, 86] amongst the
modes needs to be brought out clearly. Since the model has three degrees of
freedom, it would be relevant to map out the Arnold web structure [87, 6, 88]
in the intermediate wells and correlate with the residence time distributions
[89]. Such insights from the IVR dynamics may lead to the identification of
the “trigger” modes of the molecule that ultimately result in a concerted or
sequential mechanism. Finally, our entire study is classical and raises the
question of whether the quantum dynamics also allows for a dynamically
concerted or sequential classification. The issue is subtle since, apart from
tunneling which is relevant at low temperatures and correlated dynamics due to
quantum entanglement [65], even a proper definition of the delay time may pose
difficulties.

We conclude by noting that the Hamiltonian in eq. 7 leads to a very rich and
complex dynamics. Our study here has explored only a thin “slice” of the vast
parameter range. We hope that a more detailed phase space analysis and
classical-quantum correspondence study of the model presented herein will lead
to further insights into the dynamical implications of high index saddles on
reaction mechanisms.

##  5 Acknowledgement

We thank Antonio Fernández-Ramos and Zorka Smedarchina for clarifying certain
aspects of the potential energy surface used in their various studies of DPT.
SK’s research is supported by the Science and Engineering Research Board
(SERB) India (project no. EMR/006246). SN acknowledges the support of EPSRC
Grant No. EP/P021123/1 (CHAMPS project). PP thanks the IIT Kanpur for graduate
fellowship and the High Performance Computing Facility at IIT Kanpur for
computing resources.

##  Appendix A Double proton transfer Hamiltonian: scaled $2$D model

For a single proton transfer process, labeled as subsystem $1$ in Fig. 1, the
standard model corresponds to a quartic double well oscillator with the
Hamiltonian

$\bar{H}(\bar{x},\bar{p})=\frac{1}{2\bar{m}}\bar{p}^{2}+U(\bar{x})$ (9)

and the potential energy function

$U(\bar{x})=-a\bar{x}^{2}+b\bar{x}^{4}+U_{0}$ (10)

with $\bar{a},\bar{b}>0$ and $U_{0}$ being the barrier height for the single
barrier proton transfer. The critical points of the potential are determined
as $\bar{x}_{c}=0,\pm\Delta x_{0}$ with $\Delta x_{0}\equiv(a/2b)^{1/2}$. The
point $x_{c}=0$ corresponds to the maximum with $U(0)\equiv U_{0}=a^{2}/4b$
while $x_{c}=\pm\Delta x_{0}$ are the two minima with $U(x_{c}=\pm\Delta
x_{0})=0$. Thus, $U_{0}$ is the barrier height. As indicated in Fig. 1, the
distance between the two minima is equal to $2\Delta x_{0}$.

We introduce scaled variables as follows. The coordinate is scaled by $\Delta
x_{0}$ as $\bar{x}=x\Delta x_{0}$, with $x$ being dimensionless. Thus, the
potential energy transforms as

$\displaystyle U(x)$ $\displaystyle=$ $\displaystyle-a(\Delta
x_{0})^{2}x^{2}+b(\Delta x_{0})^{4}x^{4}+U_{0}$ (11) $\displaystyle=$
$\displaystyle U_{0}\left[x^{2}-1\right]^{2}$ (12)

Consequently, the Hamiltonian can be written down as

$\bar{H}(x,P)=\frac{(\Delta
x_{0})^{2}}{2\bar{m}}P^{2}+U_{0}\left[x^{2}-1\right]^{2}$ (13)

where $P\equiv\bar{m}\dot{x}$. Measuring mass in units of the proton mass
$m_{H}$ we have $\bar{m}=m_{H}m$ and $P=\bar{m}\dot{x}=m_{H}(m\dot{x})\equiv
m_{H}\tilde{p}$. Finally, scaling the energy by $2U_{0}$ and time as
$t=\alpha\tau$ with $\alpha=[(\Delta x_{0})^{2}m_{H}/2U_{0}]^{1/2}$ we obtain
the transformed dimensionless Hamiltonian

$H(x,p)=\frac{1}{2m}p^{2}+\frac{1}{2}\left[x^{2}-1\right]^{2}$ (14)

with the identification $p\equiv m(dx/d\tau)$.

We now consider the double proton transfer scenario shown in Fig. 1 wherein
the system has two such equivalent protons tunneling sites. A two degree of
freedom Hamiltonian can then be expressed in terms of the two proton
coordinates ${\bf x}=(x_{1},x_{2})$ and their corresponding conjugate momenta
${\bf p}=(p_{1},p_{2})$

$H(\mathbf{x},\mathbf{p})=H_{0}(\mathbf{x},\mathbf{p})+U_{\rm
coup}(\mathbf{x})$ (15)

where, the zeroth-order Hamiltonian is generalized from Eqn. 14 and of the
form

$H_{0}(\mathbf{x},\mathbf{p})=\sum_{j=1,2}\left[\frac{1}{2m_{\rm
j}}p_{j}^{2}+U_{0}({x_{j}})\right]$ (16)

with

$U_{0}(x_{j})=\frac{1}{2}\left[{x_{j}}^{2}-1\right]^{2}$ (17)

The zeroth-order form is appropriate in the limit that the two protons being
transferred are not correlated. However, typically, the two proton motions are
coupled and general symmetry-based arguments indicate that the correct form of
the coupling potential is given by

$U_{\rm
coup}(\mathbf{x})=-2G{x_{1}}{x_{2}}-Dx_{1}^{2}x_{2}^{2}-C(x_{1}^{3}x_{2}+x_{1}x_{2}^{3})+\ldots$
(18)

Thus, in principle there are couplings of all order between the two modes.
However, as has been noted earlier, from a perturbative perspective the first
two leading order terms in the above expansion for $U_{\rm coup}({\bf x})$ are
sufficient to capture most of the essential dynamical features of the system.
Therefore, in what follows we take $G,D\neq 0$ and ignore the higher order
terms.

At this stage it is useful to switch from the local coordinates used above to
the normal mode coordinates $(x_{s},x_{a})$ with the transformation
$(x_{1},x_{2})=(x_{s}+x_{a},x_{s}-x_{a})$. The Hamiltonian in this new
representation is given by

$H(x_{s},x_{a},p_{s},p_{a})=\frac{1}{2M}(p_{s}^{2}+p_{a}^{2})+\left(\frac{\delta
m}{M^{2}}\right)p_{s}p_{a}+U(x_{s},x_{a})$ (19)

with $M\equiv m_{1}+m_{2}$ and $\delta m=m_{1}-m_{2}$. Note the presence of
the momentum coupling term. This term vanishes when we are looking at the
symmetric $m_{1}=m_{2}$ cases, as in the present work. However, when
considering singly substituted isotope case like $m_{1}=m_{H}$ and
$m_{2}=m_{D}$, for instance, then $\delta m\neq 0$. The potential thus
transforms into

$U(x_{s},x_{a})=\alpha_{s}\left[x_{s}^{2}-(\Delta
x_{s})^{2}\right]^{2}+\alpha_{a}\left[x_{a}^{2}-(\Delta
x_{a})^{2}\right]^{2}+2Rx_{s}^{2}x_{a}^{2}+{\cal U}(G,D)$ (20)

where we have denoted $R=3+D$ and $\alpha_{s}=\alpha_{a}\equiv\alpha=1-D$ with

$\Delta x_{s,a}&=\sqrt{\frac{1\pm G}{1-D}}\\\ {\cal
U}(G,D)=1-\alpha_{s}(\Delta x_{s})^{4}-\alpha_{a}(\Delta x_{a})^{4}$ (21)

As a final transformation, and preparation for the three degree of freedom
Hamiltonian of interest to the current work, we transform to mass-weighted
coordinates via $X_{s,a}\rightarrow\sqrt{M}x_{s,a}$. We thus obtain the
Hamiltonian

$H({\bf X},{\bf P})=\frac{1}{2}(P_{s}^{2}+P_{a}^{2})+\left(\frac{\delta
m}{M}\right)P_{s}P_{a}+U(X_{s},X_{a})$ (22)

with the potential energy term

$U({\bf X})=\bar{\alpha}_{s}\left[X_{s}^{2}-(\Delta
X_{s})^{2}\right]^{2}+\bar{\alpha}_{a}\left[X_{a}^{2}-(\Delta
X_{a})^{2}\right]^{2}+2\bar{R}X_{s}^{2}X_{a}^{2}+{\cal U}(G,D)$ (23)

In the above we have denoted $({\bf X},{\bf
P})\equiv(X_{s},X_{a},P_{s},P_{a})$ with the parameters $\Delta
X_{s,a}=\sqrt{M}\Delta x_{s,a}$,
$\bar{\alpha}_{s}=\bar{\alpha}_{a}\equiv\alpha/M^{2}$,and $\bar{R}\equiv
R/M^{2}$. The above Hamiltonian with $\delta m=0$ corresponds to Eqn. 1 in the
main article.

##  Appendix B Hamiltonian vector field and linear stability of equilibria

The Hamiltonian vector field is given by

$\displaystyle\dot{X}_{s}=$ $\displaystyle~{}P_{s}+\dfrac{\delta m}{M}P_{a}$
(24) $\displaystyle\dot{X}_{a}=$ $\displaystyle~{}P_{a}+\dfrac{\delta
m}{M}P_{s}$ $\displaystyle\dot{Y}=$ $\displaystyle~{}P_{y}$
$\displaystyle\dot{P}_{s}=$
$\displaystyle~{}-4\left(\bar{\alpha_{s}}(X_{s}^{3}-X_{s}(\Delta
X_{s})^{2})+\bar{R}X_{s}X_{a}^{2}\right)+2\lambda
X_{s}\left(Y-\dfrac{\lambda}{\omega_{Y}^{2}}X_{s}^{2}-\dfrac{\lambda^{\prime}}{\omega_{Y}^{2}}X_{a}^{2}\right)$
$\displaystyle\dot{P}_{a}=$
$\displaystyle~{}-4\left(\bar{\alpha_{a}}(X_{a}^{3}-X_{a}(\Delta
X_{a})^{2})+\bar{R}X_{a}X_{s}^{2}\right)+2\lambda^{\prime}X_{a}\left(Y-\dfrac{\lambda}{\omega_{Y}^{2}}X_{s}^{2}-\dfrac{\lambda^{\prime}}{\omega_{Y}^{2}}X_{a}^{2}\right)$
$\displaystyle\dot{P}_{y}=$
$\displaystyle~{}-\omega_{Y}^{2}\left(Y-\dfrac{\lambda}{\omega_{Y}^{2}}X_{s}^{2}-\dfrac{\lambda^{\prime}}{\omega_{Y}^{2}}X_{a}^{2}\right)$

This vector field has same number of equilibria as the two degrees of freedom
model except that each equilibria also has a third coordinate. The equilibria,
total energies, and their linear stability is summarised in the Table. 1,
where the locations of the index-$1$ saddles and their energy are given by

$\displaystyle\left(X_{s}^{\ddagger},X_{a}^{\ddagger},Y^{\ddagger}\right)$
$\displaystyle=\left(\pm\sqrt{\dfrac{\bar{\alpha_{s}}\bar{\alpha_{a}}(\Delta
X_{s})^{2}-\bar{\alpha_{a}}\bar{R}(\Delta
X_{a})^{2}}{(\bar{\alpha_{s}}\bar{\alpha_{a}}-\bar{R}^{2})}},\pm\sqrt{\dfrac{\bar{\alpha_{s}}\bar{\alpha_{a}}(\Delta
X_{a})^{2}-\bar{\alpha_{s}}\bar{R}(\Delta
X_{s})^{2}}{(\bar{\alpha_{s}}\bar{\alpha_{a}}-\bar{R}^{2})}},\dfrac{\lambda}{\omega_{Y}^{2}}(X_{s}^{\ddagger})^{2}+\dfrac{\lambda^{\prime}}{\omega_{Y}^{2}}(X_{a}^{\ddagger})^{2}\right)$
(25) $\displaystyle H^{\ddagger}=$
$\displaystyle\dfrac{1}{M^{2}(\bar{R}^{2}-\bar{\alpha_{a}}\bar{\alpha_{s}})}\left[\bar{R}^{2}\bar{\alpha_{s}}(\Delta
X_{s})^{4}+\bar{R}^{2}\bar{\alpha_{a}}(\Delta X_{a})^{4})-2\bar{R}(\Delta
X_{s}\Delta
X_{a})^{2}\bar{\alpha_{s}}\bar{\alpha_{a}}\right]+\mathcal{U}(G,D).$ (26)

The linear stability of the equilibria in Table 1 is determined by the
eigenvalue problem associated with the linearized vector field,
$\mathbb{J}\mathbf{v}=\beta\mathbf{v}$, where the Jacobian $\mathbb{J}$ is
given by Eqn. (27), $\beta,\mathbf{v}$ are the eigenvalues and eigenvectors of
the Jacobian at the equilibrium point. For the reactant, product, and
intermediate wells, the eigenvalues are of the form
$\pm\omega_{c1}^{e},\pm\omega_{c2}^{e},\pm\omega_{c3}^{e}$. For the saddles at
energy $H^{\ddagger}$, the eigenvalues are of the form
$\pm\lambda^{e},\pm\omega_{s1}^{e},\pm\omega_{s2}^{e}$ which makes these
index-1 saddles. For the saddle at energy $E_{s}=1$, the eigenvalues are of
the form $\pm\lambda_{s1}^{e},\pm\lambda_{s2}^{e},\pm\omega_{s3}^{e}$ which
makes this an index-2 saddle.

The Jacobian of the vector field,
$\mathbb{J}(X_{s},X_{a},Y,P_{s},P_{a},P_{y})$, is given by

$\displaystyle\begin{pmatrix}0&0&0&1&\dfrac{\delta m}{M}&0\\\
0&0&0&\dfrac{\delta m}{M}&1&0\\\ 0&0&0&0&0&1\\\
-\dfrac{\partial^{2}V}{\partial X_{s}^{2}}&-\dfrac{\partial^{2}V}{\partial
X_{a}\partial X_{s}}&-\dfrac{\partial^{2}V}{\partial Y\partial X_{s}}&0&0&0\\\
-\dfrac{\partial^{2}V}{\partial X_{s}\partial
X_{a}}&-\dfrac{\partial^{2}V}{\partial
X_{a}^{2}}&-\dfrac{\partial^{2}V}{\partial Y\partial X_{a}}&0&0&0\\\
-\dfrac{\partial^{2}V}{\partial X_{s}\partial
Y}&-\dfrac{\partial^{2}V}{\partial X_{a}\partial
Y}&-\dfrac{\partial^{2}V}{\partial Y^{2}}&0&0&0\end{pmatrix},$ (27)

where

$\displaystyle-\dfrac{\partial^{2}V}{\partial X_{s}^{2}}=$
$\displaystyle-4\left[\bar{\alpha}_{s}(3X_{s}^{2}-(\Delta
X_{s})^{2})+\bar{R}X_{a}^{2}\right]+2\lambda\left(Y-\dfrac{\lambda}{\omega_{Y}^{2}}3X_{s}^{2}-\dfrac{\lambda^{\prime}}{\omega_{Y}^{2}}X_{a}^{2}\right)$
(28) $\displaystyle-\dfrac{\partial^{2}V}{\partial X_{a}^{2}}=$
$\displaystyle-4\left[\bar{\alpha}_{a}(3X_{a}^{2}-(\Delta
X_{a})^{2})+\bar{R}X_{s}^{2}\right]+2\lambda^{\prime}\left(Y-\dfrac{\lambda}{\omega_{Y}^{2}}X_{s}^{2}-\dfrac{\lambda^{\prime}}{\omega_{Y}^{2}}3X_{a}^{2}\right)$
(29) $\displaystyle-\dfrac{\partial^{2}V}{\partial X_{a}\partial X_{s}}=$
$\displaystyle-\dfrac{\partial^{2}V}{\partial X_{s}\partial
X_{a}}=-4\left(2\bar{R}+\dfrac{\lambda\lambda^{\prime}}{\omega_{Y}^{2}}\right)X_{s}X_{a}$
(30) $\displaystyle-\dfrac{\partial^{2}V}{\partial Y\partial X_{s}}=$
$\displaystyle-\dfrac{\partial^{2}V}{\partial X_{s}\partial
Y}=\phantom{-}2\lambda X_{s}$ (31)
$\displaystyle-\dfrac{\partial^{2}V}{\partial Y\partial X_{a}}=$
$\displaystyle-\dfrac{\partial^{2}V}{\partial X_{a}\partial
Y}=\phantom{-}2\lambda^{\prime}X_{a}$ (32)
$\displaystyle-\dfrac{\partial^{2}V}{\partial Y^{2}}=$
$\displaystyle\phantom{-}-\omega_{Y}^{2}$ (33)

Figure 10: Eigenvalues variation with $\lambda^{\prime},\lambda,\omega_{Y}$.
(Top row) $\lambda=0.0$ (Middle row) $\lambda=0.3$ (Bottom row)
$\lambda^{\prime}=0.0$. Left, middle, and right columns correspond to the
index-1 saddles, intermediate wells, and product (or reactant) well,
respectively. The magnitude of the eigenvalues at $\omega_{Y}=0.2$ and
$\omega_{Y}=0.6$ are shown as continuous line and dash-dot line, respectively,
with the shaded region between the lines denoting the variation for
$0.2\leq\omega_{Y}\leq 0.6$. Other parameters are $\delta
M=0,M=2,G=0.20,D=0.15$.

We track the changes in the linear stability of the index-1 and index-2
saddles with changes in the coupling strength by tracking the eigenvalues of
the linearized vector field (Jacobian (27)) evaluated at the equilibrium
points as the coupling strength is continuously varied. Note that in this work
we have $\delta m=0$.

The eigenvalues at the phase space point $(0,0,0,0,0,0)$ are given by

$\left[\pm 2\sqrt{\bar{\alpha}_{s}}\Delta X_{s},\pm
2\sqrt{\bar{\alpha}_{a}}\Delta X_{a},\pm i\omega_{Y}\right]$ (34)

which has the structure of an index-2 saddle and is only dependent on the
parameters $D,G,C,\omega_{Y}$. This supports the parametric study in this work
where the index-2 saddle maintains its linear stability as we vary
$\lambda,\lambda^{\prime}$. Further, this also points to the fact that mere
linear (local in the neighborhood of the equilibrium point) analysis will not
reflect the dynamical mechanism due to the influence of the coupling
parameters on the fraction of concerted trajectories. For the remaining
equilibria, we show the changes in the magnitude of the eigenvalues along
coupling parameters $\lambda^{\prime}$, $\lambda$ and third mode frequency
$\omega_{Y}$ for $\lambda=0.3$, $\lambda=0.0$ and $\lambda^{\prime}=0.0$ in
Fig. 10. We observe that the eigenvalues of the intermediate wells are
independent of $\lambda$ while the eigenvalues of the product (or reactant)
well are independent of $\lambda^{\prime}$. In general, there are no critical
changes in the eigenvalues as the coupling strengths and frequency of the
third mode are varied. The stability type of the equilibria stays the same.

##  Appendix C Lagrangian descriptor: method to reveal the invariant manifolds

We briefly describe the method of Lagrangian descriptors, which reveals
regions with qualitatively distinct dynamical behavior by showing the
intersection of the invariant manifolds with the two dimensional section. For
a general time-dependent dynamical system given by

$\dfrac{d\mathbf{x}}{dt}=\mathbf{f}(\mathbf{x},t)\;,\quad\mathbf{x}\in\mathbb{R}^{n}\;,\;t\in\mathbb{R}\;,$
(35)

where the vector field $\mathbf{f}(\mathbf{x},t)$ is assumed to be
sufficiently smooth both in space and time. The vector field $\mathbf{f}$ can
be prescribed by an analytical model or given from numerical simulations as a
discrete spatio-temporal data set. For instance, the vector field could
represent the velocity field of oceanic or atmospheric currents obtained from
satellite measurements or from the numerical solution of geophysical models.
For any initial condition $\mathbf{x}(t_{0})=\mathbf{x}_{0}$, the system of
first order nonlinear differential equations (given in Eqn. (35)) has a unique
solution represented by the trajectory that starts from that initial point
$\mathbf{x}_{0}$ at time $t_{0}$.

In this study, we adopt the LD definition

$\mathcal{L}_{p}(\mathbf{x}_{0},t_{0},\tau)=\int^{t_{0}+\tau}_{t_{0}-\tau}\,\sum_{k=1}^{n}|f_{k}(\mathbf{x}(t;\mathbf{x}_{0}),t)|^{p}\;dt\;,\quad
p\in(0,1]$ (36)

where $f_{k}$ is the $k-$the component of the vector field, Eqn. (35) and use
$p=1/2$. We note that the integral can be split into its forward and backward
time parts to detect the intersection of stable and unstable manifolds
separately. This relates to finding the escape and entry channels into the
potential well. In this study, we keep the forward part of the integral given
by

$\begin{split}\mathcal{L}_{p}^{f}(\mathbf{x}_{0},t_{0},\tau)&=\int^{t_{0}+\tau}_{t_{0}}\sum_{k=1}^{n}|f_{k}(\mathbf{x}(t;\mathbf{x}_{0}),t)|^{p}\;dt\end{split}$
(37)

Although this definition of LD does not have an intuitive physical
interpretation as that of the arclength definition [71], it allows for a
rigorous proof that the “singular features” (non-differentiable points) in the
LD contour map identify intersections with stable and unstable invariant
manifolds [72]. Another important aspect of what is known in LD literature as
the $p$-(quasi)norm is that degrees of freedom with relevance in
escape/transition (reaction) dynamics can be decomposed and computed. This
definition was used to show that the method can be used to successfully detect
NHIMs and their stable and unstable manifolds in Hénon-Heiles Hamiltonian [90,
91]. For this system, where both fixed (or variable) integration time is used,
it has also been shown that the LD scalar field attains a minimum (or maximum)
value along with singularity at the intersections of the stable and unstable
manifolds, and given by

$\mathcal{W}^{s}(\mathbf{x}_{0},t_{0})=\textrm{argmin
}\mathcal{L}_{p}^{f}(\mathbf{x}_{0},t_{0},\tau)\;,$ (38)

where $\mathcal{W}^{s}(\mathbf{x}_{0},t_{0})$ are the stable manifolds
calculated at time $t_{0}$ and argmin denotes the phase space coordinates on
the two dimensional section that minimize the scalar field,
$\mathcal{L}_{p}^{f}(\mathbf{x}_{0},t_{0},\tau)$, over the integration time,
$\tau$. Thus, the scalar field plotted as a contour map identifies the
intersection of the stable manifold with a two dimensional section. This
ability of LD contour map to partition trajectories with different phase space
geometry is shown in the right panel of Fig. 11 as singular values of LD
identify the intersection of the manifolds with the chosen section.

Figure 11: Lagrangian descriptor (forward) on the section 8 identifying the
invariant manifolds by the singular points with minima of the contour map
shown by the one dimensional slices on the right. Other parameters are
$\lambda=\lambda^{\prime}=0,\omega_{Y}=0.2$ and $E=1.1$.

## References

  * Uzer and Miller [1991] T. Uzer, W. Miller,  Theories of intramolecular vibrational energy transfer,  Phys. Rep. 199 (1991) 73–146. doi:10.1016/0370-1573(91)90140-H.
  * Logan and Wolynes [1990] D. E. Logan, P. G. Wolynes,  Quantum localization and energy flow in many-dimensional Fermi resonant systems,  J. Chem. Phys. 93 (1990) 4994–5012. doi:10.1063/1.458637.
  * Gruebele [2003] M. Gruebele,  Mechanism and control of molecular energy flow: a modeling perspective,  Theor. Chem. Acc. 109 (2003) 53–63. doi:10.1007/s00214-002-0394-2.
  * Leitner [2015] D. M. Leitner,  Quantum ergodicity and energy flow in molecules,  Adv. Phys. 64 (2015) 445–517. doi:10.1080/00018732.2015.1109817.
  * Keshavamurthy [2013] S. Keshavamurthy,  Scaling perspective on intramolecular vibrational energy flow: Analogies, insights, and challenges,  Adv. Chem. Phys 153 (2013) 43–110. doi:10.1002/9781118571767.ch2.
  * Karmakar and Keshavamurthy [2020] S. Karmakar, S. Keshavamurthy,  Intramolecular vibrational energy redistribution and the quantum ergodicity transition: a phase space perspective,  Phys. Chem. Chem. Phys. 22 (2020) 11139–11173. doi:10.1039/D0CP01413C.
  * Waalkens and Wiggins [2010] H. Waalkens, S. Wiggins,  Geometrical models of the phase space structures governing reaction dynamics,  Regul. Chaotic Dyn. 15 (2010) 1–39. doi:10.1134/S1560354710010016.
  * Wiggins [2016] S. Wiggins,  The role of normally hyperbolic invariant manifolds (NHIMS) in the context of the phase space setting for chemical reaction dynamics,  Regul. Chaotic Dyn. 21 (2016) 621–638. doi:10.1134/S1560354716060034.
  * Pollak and Pechukas [1978] E. Pollak, P. Pechukas,  Transition states, trapped trajectories, and classical bound states embedded in the continuum,  J. Chem. Phys. 69 (1978) 1218–1226. doi:10.1063/1.436658.
  * Pechukas and Pollak [1979] P. Pechukas, E. Pollak,  Classical transition state theory is exact if the transition state is unique,  J. Chem. Phys. 71 (1979) 2062–2068. doi:10.1063/1.438575.
  * Waalkens et al. [2007] H. Waalkens, R. Schubert, S. Wiggins,  Wigner’s dynamical transition state theory in phase space: classical and quantum,  Nonlinearity 21 (2007) R1–R118. doi:10.1088/0951-7715/21/1/r01.
  * Uzer et al. [2002] T. Uzer, C. Jaffé, J. Palacián, P. Yanguas, S. Wiggins,  The geometry of reaction dynamics,  Nonlinearity 15 (2002) 957. doi:10.1088/0951-7715/15/4/301.
  * Jaffé et al. [2005] C. Jaffé, S. Kawai, J. Palacián, Y. Patricia, T. Uzer,  A new look at the transition state: Wigner’s dynamical perspective revisited,  Adv. Chem. Phys 130 (2005) 171–216. doi:10.1002/0471712531.ch3.
  * Pechukas [1976] P. Pechukas,  Statistical approximations in collision theory,  in: Dynamics of molecular collisions, Springer, 1976, pp. 269–322. doi:10.1007/978-1-4757-0644-4_6.
  * Miller [1998] W. H. Miller,  Spiers memorial lecture quantum and semiclassical theory of chemical reaction rates,  Faraday Discuss. 110 (1998) 1–21. doi:10.1039/A805196H.
  * Pollak and Talkner [2005] E. Pollak, P. Talkner,  Reaction rate theory: What it was, where is it today, and where is it going?,  Chaos 15 (2005) 026116\. doi:10.1063/1.1858782.
  * Ezra et al. [2009] G. S. Ezra, H. Waalkens, S. Wiggins,  Microcanonical rates, gap times, and phase space dividing surfaces,  J. Chem. Phys. 130 (2009) 164118. doi:10.1063/1.3119365.
  * Collins et al. [2010] P. Collins, G. S. Ezra, S. Wiggins,  Phase space structure and dynamics for the Hamiltonian isokinetic thermostat,  J. Chem. Phys. 133 (2010) 014105. doi:10.1063/1.3455712.
  * Bartsch et al. [2005] T. Bartsch, R. Hernandez, T. Uzer,  Transition state in a noisy environment,  Phys. Rev. Lett. 95 (2005) 058301. doi:10.1103/PhysRevLett.95.058301.
  * Kawai and Komatsuzaki [2009a] S. Kawai, T. Komatsuzaki,  Dynamic pathways to mediate reactions buried in thermal fluctuations. I. Time-dependent normal form theory for multidimensional Langevin equation,  J. Chem. Phys. 131 (2009a) 224505. doi:10.1063/1.3268621.
  * Kawai and Komatsuzaki [2009b] S. Kawai, T. Komatsuzaki,  Dynamic pathways to mediate reactions buried in thermal fluctuations. II. Numerical illustrations using a model system,  J. Chem. Phys. 131 (2009b) 224506. doi:10.1063/1.3268622.
  * Feldmaier et al. [2020] M. Feldmaier, J. Reiff, R. M. Benito, F. Borondo, J. Main, R. Hernandez,  Influence of external driving on decays in the geometry of the LiCN isomerization,  J. Chem. Phys. 153 (2020) 084115. doi:10.1063/5.0015509.
  * Çiftçi and Waalkens [2013] Ü. Çiftçi, H. Waalkens,  Reaction dynamics through kinetic transition states,  Phys. Rev. Lett. 110 (2013) 233201. doi:10.1103/PhysRevLett.110.233201.
  * Tantillo [2018] D. J. Tantillo, Applied Theoretical Organic Chemistry, World Scientific, 2018.
  * Carpenter [1985] B. K. Carpenter,  Trajectories through an intermediate at a fourfold branch point. implications for the stereochemistry of biradical reactions,  J. Am. Chem. Soc. 107 (1985) 5730–5732. doi:10.1021/ja00306a021.
  * Carpenter [1995] B. K. Carpenter,  Dynamic matching: The cause of inversion of configuration in the [1, 3] sigmatropic migration?,  J. Am. Chem. Soc. 117 (1995) 6336–6344. doi:10.1021/ja00128a024.
  * Collins et al. [2014] P. Collins, Z. C. Kramer, B. K. Carpenter, G. S. Ezra, S. Wiggins,  Nonstatistical dynamics on the caldera,  J. Chem. Phys. 141 (2014) 034111. doi:10.1063/1.4889780.
  * Katsanikas and Wiggins [2018] M. Katsanikas, S. Wiggins,  Phase space structure and transport in a caldera potential energy surface,  Int. J. Bifurc. Chaos 28 (2018) 1830042. doi:10.1142/S0218127418300422.
  * Katsanikas et al. [2020] M. Katsanikas, V. J. García-Garrido, S. Wiggins,  The dynamical matching mechanism in phase space for caldera-type potential energy surfaces,  Chem. Phys. Lett. 743 (2020) 137199. doi:10.1016/j.cplett.2020.137199.
  * Geng et al. [2021] Y. Geng, M. Katsanikas, M. Agaoglou, S. Wiggins,  The influence of a pitchfork bifurcation of the critical points of a symmetric caldera potential energy surface on dynamical matching,  Chem. Phys. Lett. 768 (2021) 138397. doi:10.1016/j.cplett.2021.138397.
  * Rehbein and Carpenter [2011] J. Rehbein, B. K. Carpenter,  Do we fully understand what controls chemical selectivity?,  Phys. Chem. Chem. Phys. 13 (2011) 20906–20922. doi:10.1039/C1CP22565K.
  * Collins et al. [2013] P. Collins, B. K. Carpenter, G. S. Ezra, S. Wiggins,  Nonstatistical dynamics on potentials exhibiting reaction path bifurcations and valley-ridge inflection points,  J. Chem. Phys. 139 (2013) 154108. doi:10.1063/1.4825155.
  * Thomas et al. [2008] J. B. Thomas, J. R. Waas, M. Harmata, D. A. Singleton,  Control elements in dynamically determined selectivity on a bifurcating surface,  J. Am. Chem. Soc. 130 (2008) 14544–14555. doi:10.1021/ja802577v.
  * Kurouchi and Singleton [2018] H. Kurouchi, D. A. Singleton,  Labelling and determination of the energy in reactive intermediates in solution enabled by energy-dependent reaction selectivity,  Nat. Chem. 10 (2018) 237–241. doi:10.1038/nchem.2907.
  * Kurouchi et al. [2016] H. Kurouchi, I. L. Andujar-De Sanctis, D. A. Singleton,  Controlling selectivity by controlling energy partitioning in a thermal reaction in solution,  J. Am. Chem. Soc. 138 (2016) 14534–14537. doi:10.1021/jacs.6b09052.
  * Quijano and Singleton [2011] L. M. M. Quijano, D. A. Singleton,  Competition between reaction and intramolecular energy redistribution in solution: observation and nature of nonstatistical dynamics in the ozonolysis of vinyl ethers,  J. Am. Chem. Soc. 133 (2011) 13824–13827. doi:10.1021/ja2043497.
  * Pandey et al. [2021] P. Pandey, S. Naik, S. Keshavamurthy,  Classical and quantum dynamical manifestations of index-2 saddles: Concerted versus sequential reaction mechanisms,  Regul. Chaotic Dyn. 26 (2021) 165–182. doi:10.1134/S1560354721020052.
  * Teramoto et al. [2011] H. Teramoto, M. Toda, T. Komatsuzaki,  Dynamical switching of a reaction coordinate to carry the system through to a different product state at high energies,  Phys. Rev. Lett. 106 (2011) 054101. doi:10.1103/PhysRevLett.106.054101.
  * Accardi et al. [2010] A. Accardi, I. Barth, O. Kühn, J. Manz,  From synchronous to sequential double proton transfer: Quantum dynamics simulations for the model porphine,  J. Phys. Chem. A 114 (2010) 11252–11262. doi:10.1021/jp103435d.
  * Goldstein et al. [1996] E. Goldstein, B. Beno, K. Houk,  Density functional theory prediction of the relative energies and isotope effects for the concerted and stepwise mechanisms of the Diels- Alder reaction of butadiene and ethylene,  J. Am. Chem. Soc. 118 (1996) 6036–6043. doi:10.1021/ja9601494.
  * Pham and Houk [2014] H. V. Pham, K. Houk,  Diels-Alder reactions of allene with benzene and butadiene: Concerted, stepwise, and ambimodal transition states,  J. Org. Chem. 79 (2014) 8968–8976. doi:10.1021/jo502041f.
  * Houk et al. [2021] K. N. Houk, F. Liu, Z. Yang, J. I. Seeman,  Evolution of the Diels-Alder reaction mechanism since the 1930s: Woodward, Houk with Woodward, and the influence of computational chemistry on understanding cycloadditions,  Angew. Chem. Int. Ed. 60 (2021) 12660–12681. doi:10.1002/anie.202001654.
  * Black et al. [2012] K. Black, P. Liu, L. Xu, C. Doubleday, K. N. Houk,  Dynamics, transition states, and timing of bond formation in Diels-Alder reactions,  Proc. Natl. Acad. Sci. U.S.A 109 (2012) 12860–12865. doi:10.1073/pnas.1209316109.
  * Takeuchi and Tahara [2007] S. Takeuchi, T. Tahara,  The answer to concerted versus step-wise controversy for the double proton transfer mechanism of 7-azaindole dimer in solution,  Proc. Natl. Acad. Sci. U.S.A 104 (2007) 5285–5290. doi:10.1073/pnas.0610141104.
  * Ushiyama and Takatsuka [2001] H. Ushiyama, K. Takatsuka,  Successive mechanism of double-proton transfer in formic acid dimer: A classical study,  J. Chem. Phys. 115 (2001) 5903–5912. doi:10.1063/1.1398090.
  * Homayoon et al. [2014] Z. Homayoon, J. M. Bowman, F. A. Evangelista,  Calculations of mode-specific tunneling of double-hydrogen transfer in porphycene agree with and illuminate experiment,  J. Phys. Chem. Lett. 5 (2014) 2723–2727. doi:10.1021/jz501482v.
  * Abdel-Latif and Kühn [2011] M. K. Abdel-Latif, O. Kühn,  Laser control of double proton transfer in porphycenes: towards an ultrafast switch for photonic molecular wires,  Theor. Chem. Acc. 128 (2011) 307–316. doi:10.1007/s00214-010-0847-y.
  * Abdel-Latif and Kühn [2010] M. K. Abdel-Latif, O. Kühn,  Infrared laser driven double proton transfer. an optimal control theory study,  Chem. Phys. 368 (2010) 76–82. doi:10.1016/j.chemphys.2009.12.021.
  * Nagahata et al. [2013] Y. Nagahata, H. Teramoto, C.-B. Li, S. Kawai, T. Komatsuzaki,  Reactivity boundaries for chemical reactions associated with higher-index and multiple saddles,  Phys. Rev. E 88 (2013) 042923. doi:10.1103/PhysRevE.88.042923.
  * Pradhan and Lourderaj [2019] R. Pradhan, U. Lourderaj,  Can reactions follow non-traditional second-order saddle pathways avoiding transition states?,  Phys. Chem. Chem. Phys. 21 (2019) 12837–12842. doi:10.1039/C9CP02431J.
  * Lu et al. [2014] Z. Lu, Y. C. Chang, Q.-Z. Yin, C. Y. Ng, W. M. Jackson,  Evidence for direct molecular oxygen production in $CO_{2}$ photodissociation,  Science 346 (2014) 61–64. doi:10.1126/science.1257156.
  * Quapp and Maria Bofill [2015] W. Quapp, J. Maria Bofill,  Embedding of the saddle point of index two on the PES of the ring opening of cyclobutene,  Int. J. Quantum Chem. 115 (2015) 1635–1649. doi:10.1002/qua.24996.
  * Rashmi et al. [2021] R. Rashmi, K. Yadav, U. Lourderaj, M. Paranjothy,  Second-order saddle dynamics in isomerization reaction,  Regul. Chaotic Dyn. 26 (2021) 119–130. doi:10.1134/S1560354721020027.
  * Mauguière et al. [2014] F. A. Mauguière, P. Collins, G. S. Ezra, S. C. Farantos, S. Wiggins,  Multiple transition states and roaming in ion–molecule reactions: a phase space perspective,  Chem. Phys. Lett. 592 (2014) 282–287. doi:10.1016/j.cplett.2013.12.051.
  * Mauguière et al. [2015] F. A. Mauguière, P. Collins, Z. C. Kramer, B. K. Carpenter, G. S. Ezra, S. C. Farantos, S. Wiggins,  Phase space structures explain hydrogen atom roaming in formaldehyde decomposition,  J. Phys. Chem. Lett. 6 (2015) 4123–4128. doi:10.1021/acs.jpclett.5b01930.
  * Mauguière et al. [2016] F. A. Mauguière, P. Collins, Z. C. Kramer, B. K. Carpenter, G. S. Ezra, S. C. Farantos, S. Wiggins,  Phase space barriers and dividing surfaces in the absence of critical points of the potential energy: Application to roaming in ozone,  J. Chem. Phys. 144 (2016) 054107. doi:10.1063/1.4940798.
  * Mauguière et al. [2017] F. A. Mauguière, P. Collins, Z. C. Kramer, B. K. Carpenter, G. S. Ezra, S. C. Farantos, S. Wiggins,  Roaming: A phase space perspective,  Annu. Rev. Phys. Chem. 68 (2017) 499–524. doi:10.1146/annurev-physchem-052516-050613.
  * Montoya and Wiggins [2020] F. G. Montoya, S. Wiggins,  Revealing roaming on the double morse potential energy surface with Lagrangian descriptors,  J. Phys. A: Math. Theor. 53 (2020) 235702. doi:10.1088/1751-8121/ab8b75.
  * Kramer et al. [2015] Z. C. Kramer, B. K. Carpenter, G. S. Ezra, S. Wiggins,  Reaction path bifurcation in an electrocyclic reaction: ring-opening of the cyclopropyl radical,  J. Phys. Chem. A 119 (2015) 6611–6630. doi:10.1021/acs.jpca.5b02834.
  * Smedarchina et al. [2006] Z. Smedarchina, W. Siebrand, A. Fernández-Ramos,  Multiple proton transfer: From stepwise to concerted,  Hydrogen-Transfer Reactions (2006) 895–945. doi:10.1002/9783527611546.ch29.
  * Yoshikawa et al. [2010] T. Yoshikawa, S. Sugawara, T. Takayanagi, M. Shiga, M. Tachikawa,  Theoretical study on the mechanism of double proton transfer in porphycene by path-integral molecular dynamics simulations,  Chem. Phys. Lett. 496 (2010) 14–19. doi:10.1016/j.cplett.2010.07.009.
  * Yoshikawa et al. [2012] T. Yoshikawa, S. Sugawara, T. Takayanagi, M. Shiga, M. Tachikawa,  Quantum tautomerization in porphycene and its isotopomers: Path-integral molecular dynamics simulations,  Chem. Phys. 394 (2012) 46–51. doi:10.1016/j.chemphys.2011.12.007.
  * Walewski et al. [2010] Ł. Walewski, J. Waluk, B. Lesyng,  Car- Parrinello molecular dynamics study of the intramolecular vibrational mode-sensitive double proton-transfer mechanisms in porphycene,  J. Phys. Chem. A 114 (2010) 2313–2318. doi:10.1021/jp907754r.
  * Smedarchina et al. [2007] Z. Smedarchina, W. Siebrand, A. Fernández-Ramos,  Correlated double-proton transfer. I. Theory,  J. Chem. Phys. 127 (2007) 174513. doi:10.1063/1.2785186.
  * Smedarchina et al. [2018] Z. Smedarchina, W. Siebrand, A. Fernández-Ramos,  Entanglement and co-tunneling of two equivalent protons in hydrogen bond pairs,  J. Chem. Phys. 148 (2018) 102307. doi:10.1063/1.5000681.
  * Smedarchina et al. [2014] Z. Smedarchina, W. Siebrand, A. Fernández-Ramos,  Tunneling splitting in double-proton transfer: Direct diagonalization results for porphycene,  J. Chem. Phys. 141 (2014) 174312. doi:10.1063/1.4900717.
  * Wang et al. [2009] Z. Wang, J. S. Hirschi, D. A. Singleton,  Recrossing and dynamic matching effects on selectivity in a Diels-Alder reaction,  Angew. Chem. Int. Ed. 48 (2009) 9156–9159. doi:10.1002/anie.200903293.
  * Horn et al. [1996] B. Horn, J. Herek, A. Zewail,  Retro-Diels- Alder femtosecond reaction dynamics,  J. Am. Chem. Soc. 118 (1996) 8755–8756. doi:10.1021/ja9620696.
  * Xu et al. [2010] L. Xu, C. E. Doubleday, K. Houk,  Dynamics of 1, 3-dipolar cycloadditions: energy partitioning of reactants and quantitation of synchronicity,  J. Am. Chem. Soc. 132 (2010) 3029–3037. doi:10.1021/ja909372f.
  * Yang et al. [2019] Z. Yang, C. S. Jamieson, X.-S. Xue, M. Garcia-Borràs, T. Benton, X. Dong, F. Liu, K. Houk,  Mechanisms and dynamics of reactions involving entropic intermediates,  Trends Chem. 1 (2019) 22–34. doi:10.1016/j.trechm.2019.01.009.
  * Mancho et al. [2013] A. M. Mancho, S. Wiggins, J. Curbelo, C. Mendoza,  Lagrangian descriptors: A method for revealing phase space structures of general time dependent dynamical systems,  Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 3530–3557. doi:10.1016/j.cnsns.2013.05.002.
  * Lopesino et al. [2017] C. Lopesino, F. Balibrea-Iniesta, V. J. García-Garrido, S. Wiggins, A. M. Mancho,  A Theoretical Framework for Lagrangian Descriptors,  Int. J. Bifurcat. Chaos 27 (2017) 1730001. doi:10.1142/S0218127417300014.
  * Agaoglou et al. [2020] M. Agaoglou, B. Aguilar-Sanjuan, V. J. García Garrido, F. González-Montoya, M. Katsanikas, V. Krajňák, S. Naik, S. Wiggins, Lagrangian Descriptors: Discovery and Quantification of Phase Space Structure and Transport, Zenodo, 2020. URL: https://champsproject.github.io/lagrangian_descriptors/. doi:10.5281/zenodo.3958985, EPSRC Grant Number: EP/P021123/1.
  * Naik et al. [2019] S. Naik, V. J. García-Garrido, S. Wiggins,  Finding NHIM: Identifying high dimensional phase space structures in reaction dynamics using Lagrangian descriptors,  Commun. Nonlinear Sci. Numer. Simul. 79 (2019) 104907. doi:10.1016/j.cnsns.2019.104907.
  * Naik and Wiggins [2019] S. Naik, S. Wiggins,  Finding normally hyperbolic invariant manifolds in two and three degrees of freedom with Hénon-Heiles-type potential,  Phys. Rev. E 100 (2019) 022204. doi:10.1103/PhysRevE.100.022204.
  * Naik and Wiggins [2020] S. Naik, S. Wiggins,  Detecting reactive islands in a system-bath model of isomerization,  Phys. Chem. Chem. Phys. 22 (2020) 17890–17912. doi:10.1039/D0CP01362E.
  * Eklund and Ananth [2021] E. C. Eklund, N. Ananth,  Investigating the Stability and Accuracy of a Classical Mapping Variable Hamiltonian for Nonadiabatic Quantum Dynamics,  Regul. Chaotic Dyn. 26 (2021) 131–146. doi:10.1134/S1560354721020039.
  * Wiggins [2014] S. Wiggins, Normally hyperbolic invariant manifolds in dynamical systems., Springer, 2014\.
  * Ezra and Wiggins [2009] G. S. Ezra, S. Wiggins,  Phase-space geometry and reaction dynamics near index 2 saddles,  J. Phys. A: Math. Theor. 42 (2009) 205101. doi:10.1088/1751-8113/42/20/205101.
  * Collins et al. [2011] P. Collins, G. S. Ezra, S. Wiggins,  Index k saddles and dividing surfaces in phase space with applications to isomerization dynamics,  J. Chem Phys. 134 (2011) 244105. doi:10.1063/1.3602465.
  * Ezra and Wiggins [2009] G. S. Ezra, S. Wiggins,  Phase-space geometry and reaction dynamics near index 2 saddles,  J. Phys. A Math. Theor. 42 (2009) 205101. doi:10.1088/1751-8113/42/20/205101.
  * Slater [1956] N. B. Slater,  New formulation of gaseous unimolecular dissociation rates,  J. Chem. Phys. 24 (1956) 1256–1257. doi:10.1063/1.1742756.
  * Slater [1959] N. B. Slater, Theory of unimolecular reactions, Cornell University Press, 1959\.
  * Thiele [1962] E. Thiele,  Comparison of the classical theories of unimolecular reactions,  J. Chem Phys. 36 (1962) 1466–1472. doi:10.1063/1.1732765.
  * Thiele [1963] E. Thiele,  Comparison of the classical theories of unimolecular reactions. II. A model calculation,  J. Chem Phys. 38 (1963) 1959–1966. doi:10.1063/1.1733903.
  * Nesbitt and Field [1996] D. J. Nesbitt, R. W. Field,  Vibrational energy flow in highly excited molecules: Role of intramolecular vibrational redistribution,  J. Phys. Chem. 100 (1996) 12735–12756. doi:10.1021/jp960698w.
  * Karmakar and Keshavamurthy [2018] S. Karmakar, S. Keshavamurthy,  Relevance of the resonance junctions on the Arnold web to dynamical tunneling and eigenstate delocalization,  J. Chem Phys. 122 (2018) 8636–8649. doi:10.1021/acs.jpca.8b08626.
  * Karmakar et al. [2020] S. Karmakar, P. K. Yadav, S. Keshavamurthy,  Stable chaos and delayed onset of statisticality in unimolecular dissociation reactions,  Commun. Chem. 3 (2020) 1–11. doi:10.1038/s42004-019-0252-y.
  * Shojiguchi et al. [2007] A. Shojiguchi, C.-B. Li, T. Komatsuzaki, M. Toda,  Fractional behavior in multidimensional Hamiltonian systems describing reactions,  Phys. Rev. E 76 (2007) 056205. doi:10.1103/PhysRevE.76.056205.
  * Demian and Wiggins [2017] A. S. Demian, S. Wiggins,  Detection of periodic orbits in Hamiltonian systems using Lagrangian descriptors,  Int. J. Bifurcat. Chaos 27 (2017) 1750225. doi:10.1142/S021812741750225X.
  * Naik and Wiggins [2019] S. Naik, S. Wiggins,  Finding normally hyperbolic invariant manifolds in two and three degrees of freedom with Hénon-Heiles type potential,  Phys. Rev. E 100 (2019) 022204. doi:10.1103/PhysRevE.100.022204.