# Reactive Islands for Three Degrees-of-Freedom Hamiltonian Systems

Vladimír Krajňák<EMAIL_ADDRESS>corresponding author School of
Mathematics, University of Bristol,
Fry Building, Woodland Road, Bristol, BS8 1UG, United Kingdom.  Víctor J.
García-Garrido<EMAIL_ADDRESS>Departamento de Física y Matemáticas,
Universidad de Alcalá,
Madrid, 28871, Spain.  Stephen Wiggins<EMAIL_ADDRESS>School of
Mathematics, University of Bristol,
Fry Building, Woodland Road, Bristol, BS8 1UG, United Kingdom.

###### Abstract

We develop the geometrical, analytical, and computational framework for
reactive island theory for three degrees-of-freedom time-independent
Hamiltonian systems. In this setting, the dynamics occurs in a 5-dimensional
energy surface in phase space and is governed by four-dimensional stable and
unstable manifolds of a three-dimensional normally hyperbolic invariant
sphere. The stable and unstable manifolds have the geometrical structure of
spherinders and we provide the means to investigate the ways in which these
spherinders and their intersections determine the dynamical evolution of
trajectories. This geometrical picture is realized through the computational
technique of Lagrangian descriptors. In a set of trajectories, Lagrangian
descriptors allow us to identify the ones closest to a stable or unstable
manifold. Using an approximation of the manifold on a surface of section we
are able to calculate the flux between two regions of the energy surface.

Keywords: Phase space, Hamiltonian Systems, Stable and Unstable Manifolds,
Normally Hyperbolic Invariant Manifolds, Reactive Islands, Spherinders,
Lagrangian descriptors.

##  I Introduction

The deterministic evolution of a system between two stable states through an
intermediate unstable state is a fundamental setting for an important form of
dynamical evolution that informs our way of thinking of “transition
phenomena”. This dynamical situation is common to many diverse fields in
science and engineering, such as celestial mechanics[1, 2], structural
mechanics[3, 4, 5], and chemical reaction dynamics. This setting is natural
for studies in the latter area and will be the focus of the applications that
we consider.

The mathematical setting that we will consider is that of two and three
degrees-of-freedom (DoF) autonomous Hamiltonian systems (hence, four and six
dimensional phase spaces, respectively). The Hamiltonians have the form of the
sum of potential energy and kinetic energy, where the potential energy is a
function of the position variables (“configuration space”) and the kinetic
energy is a function of the momentum variables (variables conjugate to the
position variables). The topography of the potential energy function, or
“potential energy surface” (PES) plays an important role in the reaction
problem that we study. Potential wells (i.e. minima of the PES) are identified
with stable states, or “reactants” and “products”, and they are separated by
saddle points on the PES. There are different types of saddle points
characterized by their index [6]. In this paper we will be concerned with
index-one saddle points as this is the most common saddle point considered in
this type of transition phenomena. Trajectories passing between potential
wells are dynamical objects that are most naturally studied in phase space,
i.e., the space of configuration and momentum variables. The situation is very
different for two and three DoF systems, and we will describe each separately.

###  I.1 Two Degrees-of-Freedom Hamiltonian Systems

The geometry and dynamics of the phase space in a neighborhood of an index-one
saddle was first worked out by Conley for two DoF systems in his studies of
the restricted three body problem [7, 8], and we briefly describe that
situation. The Lyapunov subcenter theorem [9, 10] states that for a range of
energies above that of the index-one saddle (the precise range is not
specified in the theorem) there exists an unstable periodic orbit (UPO) having
two dimensional stable and unstable manifolds. In the three dimensional energy
surface the stable and unstable manifolds have the geometry of cylinders
(“tubes”) and these stable and unstable cylinders mediate transitions between
the wells as we now describe.

As a result of the invariance of the stable and unstable manifolds, and their
cylindrical structure, trajectories that start within the tubes must remain in
the tubes for all positive and negative time. All trajectories inside the
stable cylinder approach the UPO in positive time, pass through the regions
bounded by the UPO, and exit the region through the unstable cylinder. The
stable cylinder (and the unstable cylinder) has two “branches” joined at the
UPO. Considering a neighborhood of the UPO in the energy surface, one branch
is in the reactant region and the other is in the product region.

All trajectories that react must lie in a branch of the stable cylinder.
Trajectories starting in the branch of the stable cylinder lying in the
reactant region correspond to forward reacting trajectories, i.e. trajectories
that evolve from reactants to products in forward time. Trajectories starting
in the branch of the stable cylinder lying in the product region correspond to
backward reacting trajectories, i.e. trajectories that evolve from products to
reactants in forward time. Hence the stable cylinders form the pathways for
reaction in phase space. They can be used to monitor reacting trajectories in
the following manner. For definiteness we consider the branch of the stable
cylinder in the reactant region. Similar reasoning applies to the branch of
the stable cylinder in the product region.

A dividing surface (DS) can be constructed to monitor, or “count”, the
trajectories passing from reactants to products. In order for trajectories to
be counted only once in their evolution from reactants to products they must
cross the DS only once. A DS with this property is said to have the “no-
recrossing” property. Seminal work of Pollak, Pechukas, and Child [11, 12, 13,
14, 15] showed how to use this UPO to construct a DS with these properties, a
so-called “periodic orbit dividing surface”, or PODS. The DS is topologically
a 2 sphere. The equator of the 2 sphere is the UPO which divides that 2 sphere
into two halves. One half is intersected by trajectories contained in the
branch of the stable cylinder in the reactant region (forward reaction) and
the other half is intersected by trajectories contained in the branch of the
stable cylinder in the product region (backward reaction).

The two dimensional invariant cylinders can exhibit a complicated geometrical
structure in the three dimensional energy surface. We describe a lower
dimensional method for monitoring their evolution towards reaction as follows.
In the region of the energy surface corresponding to reactants (a similar
construction can be carried out for the product region) we construct a two
dimensional Poincaré section, i.e. a two dimensional surface where the
Hamiltonian vector field is everywhere transverse to the surface. The Poincaré
section is constructed in such a way that the stable cylinder intersects it in
a topological circle [16, 14]. The region bounded by this topological circle
is referred to as a “reactive island” [17]. The return map of the Poincaré
section into itself is the map that associates to a point its first return to
the Poincaré section under the flow generated by the Hamiltonian vector field.
The inverse of this return map of the Poincaré section into itself is the map
that associates to a point its first return to the Poincaré section under
negative time, i.e. the point “where it came from”.

If one considers the preimages of this reactive island by letting it evolve
backwards under the inverse of the Poincaré map one obtains a reactive island
on the Poincaré section that returns to the original reactive island, and then
reacts by crossing the DS for positive time evolution. By repeating this
construction we obtain a temporal ordered sequence of reactive islands which
sequentially map to each other before reacting [18]. We would like to remark
that pathological intersections of the cylinders with the Poincaré section can
occur, and are discussed in [19].

We illustrate in Fig. 1 the concept of reactive islands for 2 DoF Hamiltonian
systems on a suitable Poincaré section. The intricate geometry of the reactive
islands, obtained when the stable and unstable manifolds originating from an
unstable periodic orbit intersect transversely an adequately chosen phase
space section, is revealed using the method of Lagrangian descriptors
calculated on that surface of section. This technique provides a clear and
highly-detailed visualization of the reactive island picture, offering an
understanding of the nonlinear interactions between the manifolds that exist
in the phase space of the Hamiltonian system.

Consider for example a double well Hamiltonian defined by a potential energy
surface, where the two wells are separated by an index-1 saddle that lies
between them. In the diagram displayed in Fig. 1, we have represented the
boundary of the energy surface where motion takes place in green, and $q_{1}$
is the reactive DoF of the system. We have depicted in magenta a reactive
trajectory that evolves in forward time from the ‘reactant’ to the ‘product’
well regions of the PES, by moving across the phase space bottleneck that
exists in the neighborhood of the UPO associated to the index-1 saddle. This
trajectory starts inside the spherical cylinder corresponding to the stable
manifold on the right side of the phase space, and once it crosses the
dividing surface associated to the UPO, it enters the left well region
following the unstable manifold. We show the reactive islands formed by the
tube manifolds of the system by computing LDs on the Poincaré sections
$\Sigma^{-}_{q_{1}=a}$ and $\Sigma^{+}_{q_{1}=b}$ taken on both well regions,
where the plus and minus sign notation represents the directionality of the
section, that is, $\dot{q}_{1}>0$ or $\dot{q}_{1}<0$ respectively. Moreover,
on the section $\Sigma^{-}_{q_{1}=a}$ LDs have been calculated using a longer
integration time in order to unveil the successive sequence of reactive
islands that arise when the manifolds intersect the section for several times.
The deformation of the successive reactive islands observed in the Poincaré
section is a fingerprint of the chaotic behavior of this system.

Figure 1: Reactive islands as revealed by the method of Lagrangian descriptors
on the phase space of a coupled double well Hamiltonian system with 2 DoF. The
cylindrical geometry of the stable and unstable manifolds emanating from the
unstable periodic orbit in the bottleneck region is highlighted when LDs are
applied on the two-dimensional Poincaré sections $\Sigma^{-}_{q_{1}=a}$ and
$\Sigma^{+}_{q_{1}=b}$. A reactive trajectory is depicted in magenta and the
intersections of its path with the phase space planes (taking into account the
directionality of the sections) is marked with cyan dots. Arrows are included
in the picture to indicate forward time evolution, and the green surface
represents the projection of the boundary of the energy surface.

The flux of trajectories across the phase space bottleneck region, where the
UPO associated to the index-1 saddle point is located, is a key ingredient for
the computation of the classical reaction rate. The flux through the DS,
$\varphi_{\text{DS}}$, can be calculated by means of applying Stokes theorem
using the action integral [20], that is:

$\varphi_{\text{DS}}=\int_{\text{DS}}d\mathbf{p}\wedge
d\mathbf{q}=\oint_{\text{UPO}}\mathbf{p}\cdot d\mathbf{q}$ (1)

where $\mathbf{q}$ represents the configuration space variables and
$\mathbf{p}$ the canonically conjugate momenta. Following the Poincaré-Cartan
theorem [21], we have:

$\oint_{\text{UPO}}\mathbf{p}\cdot d\mathbf{q}=\oint_{\Gamma}\mathbf{p}\cdot
d\mathbf{q},$ (2)

for any curve $\Gamma$ that results from a transverse intersection of the
stable or unstable manifolds of the UPO with a Poincaré section. Hence, the
flux across the dividing surface in the phase space bottleneck is the same as
the flux across any of the reactive islands.

A theory of reaction dynamics for two DoF systems based on the geometry of
these stable and unstable cylinders and “reactive islands” [17] was developed
in the late 1980’s and 90’s in [19, 17, 18, 22, 23, 24, 25, 26], although as
remarked in [18], the ideas of fluxes of trajectories in tubes in phase space
and their circular cross sections were considered in [16, 14]. The phase space
structures and geometrical considerations that we have described above are
very specific to two DoF Hamiltonian systems [18]. However, in essentially all
of the papers noted above the need and interest for an analogous theory is
described as “future work”. However, such future work has not appeared. There
are two reasons for this. One is that new types of dynamical objects are
required, and the other is that a numerically efficient method for generating
these structures is required. In recent years progress has been made along
both lines that enable us to realize a reactive island theory for three DoF
systems, which we now describe.

###  I.2 Three Degrees-of-Freedom Hamiltonian Systems

Three DoF Hamiltonian systems have a 6-dimensional phase space and dynamics is
constrained to a 5-dimensional energy surface. An UPO has 3-dimensional stable
and unstable manifolds in the 5-dimensional energy surface. Hence, they are
not codimension one, i.e. one less dimension than the energy surface, and
therefore they cannot trap trajectories or act as barriers in the phase space.

It has been shown that an appropriate generalization of an UPO for 3 DoF
Hamiltonian systems is a normally hyperbolic invariant 3-sphere having
4-dimensional stable and unstable manifolds in the 5-dimensional energy
surface [27, 28, 29, 30, 31, 32, 33]. The stable and unstable manifolds have
the geometry of a spherical cylinder, or “spherinder”, see e.g.
https://en.wikipedia.org/wiki/Spherinder. The 3-sphere can be used to
construct a DS with the same properties as the PODS for 2 DoF Hamiltonian
systems that we described earlier.

However, the main issue is how would this approach be implemented in specific
systems. In particular, how can one determine the existence and compute the
normally hyperbolic invariant 3-sphere, compute its 4-dimensional stable and
unstable manifolds (the “spherinders”), and determine, and understand, the
geometry of their intersection with a 4-dimensional Poincaré section? These
intersection sets are 3-dimensional on the Poincaré section and are the analog
of the reactive islands for the 2 DoF Hamiltonian systems discussed earlier.

For 2 DoF Hamiltonian systems there is a rich body of knowledge concerning
computing UPOs and their stable and unstable manifolds. However, the
generalization of these 2 DoF techniques to the 3 DoF setup described above is
simply not practical. For example, suppose we had a realization of the
normally hyperbolic invariant 3-sphere adequate for computational purposes
(this is a big step that we are just assuming, for the moment). The next step
would be to compute its 4-dimensional stable and unstable manifolds. The
analog of the situation for 2 DoF would be to generate a 4-dimensional mesh of
initial conditions near the linearized unstable directions of the sphere and
propagate them in time in the 5-dimensional energy surface. However, as it is
propagated gaps would quickly develop in the mesh and repeated remeshing in
space after intervals of time would be required. This has been carried out for
the computation of 2-dimensional stable and unstable manifolds in 3
dimensions, but not in higher dimensions.

We note that Poincaré maps “work” for 2 DoF systems because they are
2-dimensional and it is relatively easy to compute many trajectories that
return to the Poincaré section and give rise to structures that can be related
to dynamical concepts. For 3 DoF systems the Poincaré section is 4-dimensional
and if sufficient trajectories returned it is not obvious how one could
extract “recognizable” structures in this 4-dimensional space. Of course, one
could attempt to extract structures in the 4-dimensional Poincaré section by
considering 2-dimensional slices of the Poincaré section. However, the
probability that a trajectory with initial conditions on this 2-dimensional
slice returning to this slice in 5 dimensions is zero. For completeness we add
that 2-tori around stable periodic orbits were investigated on 4-dimensional
Poincaré sections in 3 DoF systems using color as the fourth dimension [34]
and on “thickened” 3-dimensional slices in 4-dimensional maps [35]. It is not
clear how these approaches could be used for 3-dimensional objects.

The relatively new trajectory diagnostic known as “Lagrangian descriptors” is
giving very promising results for uncovering phase space structures in high-
dimensional systems. The method was originally developed in the context of
Lagrangian transport studies in fluid dynamics [36, 37, 38], but the wide
applicability of the method has recently been recognized in chemistry [39, 40,
41, 42, 43, 44, 45, 46, 47]. The method is straightforward to implement
computationally [48], the interpretation is clear, and it provides a “high
resolution” method for exploring high-dimensional phase space with low-
dimensional slices [49, 50].

Phase space structures consist of trajectories. But as discussed above, for
high-dimensional phase space trajectory approaches for revealing phase space
structure are problematic and prone to issues of interpretation since a
tightly grouped set of initial conditions may result in trajectories that
become “lost” with respect to each other in phase space. The method of
Lagrangian descriptors turns this problem on its head by emphasizing the
initial conditions of trajectories, rather than the precise location of their
futures and pasts, after a specified amount of time. In fact, Lagrangian
descriptors have already been used to reveal reactive islands in 2 DoF
Hamiltonian systems [51, 52, 53, 45]. A low-dimensional “slice” of phase space
can be selected and sampled with a grid of initial conditions of high
resolution. Since the phase space structure is encoded in the initial
conditions of the trajectories, no resolution is lost as the trajectories
evolve in time. In this paper we will show how this can be done in order to
realize the setting of reactive islands for 3 DoF Hamiltonian systems. We will
show how the reactive islands can be extracted from Lagrangian descriptor
values and we will be able to calculate flux by means of quantifying the
volume enclosed by them.

Recall that, as we have discussed above, reactive islands in 3 DoF systems
have the topology of $S^{3}$. To facilitate the work we will decompose this
high-dimensional structures into a family of two-dimensional sphere-like
objects, $S^{2}$, embedded in 3-dimensional slices, that we triangulate and
accumulate their contributions into the flux integral. Furthermore, we will
show how to calculate the volume enclosed by the intersection of two reactive
islands, a crucial step towards determining reaction rates.

Detecting these structures in this high-dimensional setup, and using them to
compute flux, is a challenging task that we address in this work. We provide a
procedure for systematically carrying out these computations, and illustrate
how this algorithm can be applied to a simple benchmark Hamiltonian system
with 3 DoF, a double well potential which describes a reactive DoF
quadratically coupled to two harmonic oscillators.

This paper is outlined as follows. Section II introduces the 3 DoF double well
system that we will use throughout this work as a simple model to study the
concept of reactive islands, and also for analyzing the computation of flux.
In Sec. III this Hamiltonian model is extended further by adding a quadratic
coupling between the reactive and bath (i.e. harmonic oscillator) degrees-of-
freedom. Section IV focuses on describing how to select an adequate Poincaré
surface of section so that the reactive island approach for calculating flux
works. Section V is devoted to explaining how the method of Lagrangian
descriptors can be used to extract the high-dimensional reactive islands (and
hence visualize their geometry) formed by the intersections of the stable and
unstable manifolds in the system. This step is crucial in order to develop an
algorithm that allows us to determine fluxes for the uncoupled and coupled
system scenarios from the reactive islands detected by means of Lagrangian
descriptors. The computation of flux for the uncoupled Hamiltonian is
described in Sec. VI, whereas Sec. VII does the same for the coupled system.
Finally, in Sec. VIII we present the conclusions of this work and discuss some
ideas for future research lines where the methodology we have developed for
identifying high-dimensional reactive islands using Lagrangian descriptors
could provide relevant insights and lead to interesting results.

##  II Double Well Potential Hamiltonian with Three Degrees-of-Freedom

In this section we present a benchmark 3 DoF example that allows analytical
calculation of the phase space structures that play the central role for a
reactive island theory in 3 DoF Hamiltonian systems. It is an uncoupled 3 DoF
Hamiltonian system consisting on a 1 DoF double well PES and two harmonic
oscillators. This is the natural generalization of the motivational example in
[18] and follows the spirit of [33]. In particular, we consider the following
Hamiltonian system:

$\mathcal{H}(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})=\dfrac{p_{1}^{2}}{2}+\dfrac{q_{1}^{4}}{4}-\dfrac{q_{1}^{2}}{2}+\dfrac{\omega_{2}}{2}\left(q_{2}^{2}+p_{2}^{2}\right)+\dfrac{\omega_{3}}{2}\left(q_{3}^{2}+p_{3}^{2}\right)\;,$
(3)

where $\omega_{2},\,\omega_{3}>0$ are the angular frequencies of the harmonic
oscillators. Notice that we have written the contributions of the harmonic
oscillators to the Hamiltonian differently from the classical expression
containing the squares of the frequencies. In fact, both expressions are
equivalent by a canonical transformation. This approach is more convenient to
use and is similar to the convention followed in chemistry when working with
mass-weighted coordinates [54]. The corresponding Hamilton’s equations of
motion that describe the dynamics of the system are as follows:

$\begin{cases}\dot{q}_{1}=\dfrac{\partial\mathcal{H}}{\partial
p_{1}}=p_{1}\\\\[11.38092pt] \dot{q}_{2}=\dfrac{\partial\mathcal{H}}{\partial
p_{2}}=\omega_{2}p_{2}\\\\[11.38092pt]
\dot{q}_{3}=\dfrac{\partial\mathcal{H}}{\partial
p_{3}}=\omega_{3}p_{3}\\\\[11.38092pt]
\dot{p}_{1}=-\dfrac{\partial\mathcal{H}}{\partial
q_{1}}=q_{1}-q_{1}^{3}\\\\[11.38092pt]
\dot{p}_{2}=-\dfrac{\partial\mathcal{H}}{\partial
q_{2}}=-\omega_{2}q_{2}\\\\[11.38092pt]
\dot{p}_{3}=-\dfrac{\partial\mathcal{H}}{\partial
q_{3}}=-\omega_{3}q_{3}\end{cases}$ (4)

We will construct reactive islands for a fixed energy. Hence, the
5-dimensional energy surface is given by:

$\mathcal{E}_{E}=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\dfrac{p_{1}^{2}}{2}+\dfrac{q_{1}^{4}}{4}-\dfrac{q_{1}^{2}}{2}+\dfrac{\omega_{2}}{2}\left(q_{2}^{2}+p_{2}^{2}\right)+\dfrac{\omega_{3}}{2}\left(q_{3}^{2}+p_{3}^{2}\right)=E\right\\}\,.$
(5)

The following results are proven in [33]. The isoenergetic normally hyperbolic
invariant manifold (NHIM) is given by:

$\mathcal{N}_{E}=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\dfrac{\omega_{2}}{2}\left(q_{2}^{2}+p_{2}^{2}\right)+\dfrac{\omega_{3}}{2}\left(q_{3}^{2}+p_{3}^{2}\right)=E\,,\,q_{1}=p_{1}=0\right\\}\,.$
(6)

and it has the form of a 3-sphere in the 5-dimensional energy surface. In this
example the stable and unstable manifolds (“spherinders”) of the NHIM are
coincident and are given by:

$\mathcal{W}^{s,u}\left(\mathcal{N}_{E}\right)=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\dfrac{\omega_{2}}{2}(q_{2}^{2}+p_{2}^{2})+\dfrac{\omega_{3}}{2}(q_{3}^{2}+p_{3}^{2})=E\;,\;\dfrac{p_{1}^{2}}{2}+\dfrac{q_{1}^{4}}{4}-\frac{q_{1}^{2}}{2}=0\;,\;q_{1}\neq
0\right\\}\,.$ (7)

They are 4-dimensional in the 5-dimensional energy surface. The dividing
surface ($q_{1}=0$) is given by:

$\mathcal{D}_{E}=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\dfrac{1}{2}p_{1}^{2}+\dfrac{\omega_{2}}{2}\left(q_{2}^{2}+p_{2}^{2}\right)+\dfrac{\omega_{3}}{2}\left(q_{3}^{2}+p_{3}^{2}\right)=E\;,\;q_{1}=0\right\\}\,.$
(8)

and has the form of a 4-sphere in the 5-dimensional energy surface. Its
equator is given by $p_{1}=0$, which is the NHIM.

###  II.1 Intersection of $\mathcal{W}^{s,u}\left(\mathcal{N}_{E}\right)$ with
a hyperplane

As discussed in the introduction, reactive islands are the intersection of the
stable manifold of the normally hyperbolic invariant 3-sphere (“spherinder”)
with a 4-dimensional Poincaré section in the energy surface. In particular, we
will show that in an energy surface, a transverse hyperplane intersects
$\mathcal{W}^{s,u}\left(\mathcal{N}_{E}\right)$ in a 3-sphere. We define a
hyperplane that intersects $\mathcal{W}^{s,u}\left(\mathcal{N}_{E}\right)$
transversely by:

$\Sigma_{0.5}=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\mathcal{H}(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})=E\,,\,q_{1}=-0.5\,,\,\dot{q}_{1}>0\right\\}\,.$
(9)

then we have:

$\Sigma_{0.5}\cap\mathcal{W}^{s,u}\left(\mathcal{N}_{E}\right)=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\dfrac{\omega_{2}}{2}(q_{2}^{2}+p_{2}^{2})+\dfrac{\omega_{3}}{2}(q_{3}^{2}+p_{3}^{2})=E\;,\;q_{1}=-0.5\;,\;p_{1}=\sqrt{\dfrac{7}{32}}\right\\}.$
(10)

which is a 3-sphere, the higher-dimensional analogue of the reactive islands
shown on $\Sigma^{-}_{q_{1}=a}$ and $\Sigma^{+}_{q_{1}=b}$ in Fig. 1. In Fig.
2 we show several $p_{3}=\text{const}>0$ slices of this 3-sphere.

Figure 2: Section of the stable invariant manifold
$\Sigma_{0.5}\cap\mathcal{W}^{s}$ of the uncoupled system with $E=0.1$ and
$\omega_{2}=\omega_{3}=1$. The section consists of several
$p_{3}=\text{const}>0$ slices ($p_{3}<0$ is related by symmetry) of which the
first six are shown. The analytically calculated manifold in yellow is
overlayed with green points calculated numerically as detailed in Sec. V.

##  III The Double Well Potential with Quadratic Coupling

We now introduce a quadratic coupling term between the reactive DoF, $q_{1}$,
and each of the bath DoF, $q_{2}$ and $q_{3}$, for the model system described
in Sec. II. The coupled Hamiltonian has the form:

$\mathcal{H}(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})=\dfrac{p_{1}^{2}}{2}+\dfrac{q_{1}^{4}}{4}-\dfrac{q_{1}^{2}}{2}+\dfrac{\omega_{2}}{2}\left(q_{2}^{2}+p_{2}^{2}\right)+\dfrac{\omega_{3}}{2}\left(q_{3}^{2}+p_{3}^{2}\right)+\dfrac{\varepsilon}{2}\left(q_{1}-q_{2}\right)^{2}+\dfrac{\varepsilon}{2}\left(q_{1}-q_{3}\right)^{2}\;,$
(11)

where the degrees of freedom are coupled quadratically with the coupling
parameter $\varepsilon>0$. This system gives rise to the following Hamilton’s
equations:

$\begin{cases}\dot{q}_{1}=\dfrac{\partial\mathcal{H}}{\partial
p_{1}}=p_{1},\\\\[11.38092pt] \dot{q}_{2}=\dfrac{\partial\mathcal{H}}{\partial
p_{2}}=\omega_{2}p_{2}\\\\[11.38092pt]
\dot{q}_{3}=\dfrac{\partial\mathcal{H}}{\partial
p_{3}}=\omega_{3}p_{3}\\\\[11.38092pt]
\dot{p}_{1}=-\dfrac{\partial\mathcal{H}}{\partial
q_{1}}=q_{1}-q_{1}^{3}+\varepsilon\left(q_{2}-q_{1}\right)+\varepsilon\left(q_{3}-q_{1}\right)\\\\[11.38092pt]
\dot{p}_{2}=-\dfrac{\partial\mathcal{H}}{\partial
q_{2}}=-\omega_{2}q_{2}+\varepsilon\left(q_{1}-q_{2}\right)\\\\[11.38092pt]
\dot{p}_{3}=-\dfrac{\partial\mathcal{H}}{\partial
q_{3}}=-\omega_{3}q_{3}+\varepsilon\left(q_{1}-q_{3}\right)\end{cases}\;.$
(12)

While the NHIM $\mathcal{N}$ associated with the index-1 saddle is still a
$3$-sphere, its location and shape depends on the coupling strength
$\varepsilon$. It can be approximated by a nearly invariant linear combination
of $q_{1},q_{2},q_{3}$ on the energy surface $\mathcal{H}=E$. For the sake of
simplicity we derive this linear combination for $\omega_{2}=\omega_{3}=1$ as

$cq_{1}-\varepsilon(q_{2}+q_{3})=0\,,$ (13)

where

$c=\dfrac{1}{2}\left(\varepsilon-2-\sqrt{\left(\varepsilon-2\right)^{2}+8\varepsilon^{2}}\right),$
(14)

follows from the near-invariance under the equations of motion linearized
around $q_{1}=0$ (or the origin):

$cp_{1}-\varepsilon(p_{2}+p_{3})=0\,,\quad
c\dot{p}_{1}-\varepsilon(\dot{p}_{2}+\dot{p}_{3})=0.$ (15)

We define the dividing surface constructed using the approximated NHIM as:

$\mathcal{D}_{E}=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\mathcal{H}(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})=E\;,\;cq_{1}-\varepsilon(q_{2}+q_{3})=0\right\\}\,.$
(16)

In all of the following calculations for the coupled system we consider
$\varepsilon=0.1$ and $\omega_{2}=\omega_{3}=1$.

##  IV Choice of surface of section and tangencies

For the reactive island approach to work correctly, it is crucial that the
hyperplane $\Sigma$ is defined transversely to the investigated invariant
manifold. It is well known that discontinuities occur when a trajectory on an
invariant manifold is tangent to the hyperplane in systems with $2$ DoF [18],
and the only difference with 3 DoF systems involves the dimensionality of the
tangent set, which can vary from $0$ to $2$.

Figure 3: Transverse intersection of $\mathcal{W}^{s}$ (left) and non-
transverse intersection of $\mathcal{W}^{u}$ (right) with $\Sigma_{0.5}$ and
$\left\\{p_{3}\approx 1/3\right\\}$ in the coupled system with
$\varepsilon=0.1$. The hole in the non-transverse intersection is highlighted
in red.

While a transverse intersection of $\Sigma$ and $\mathcal{W}^{s,u}$ is a
3-sphere, see Eq. (10), a non-transverse intersection causes this sphere to
have holes as we show in Fig. 3. As a consequence, a portion of the invariant
manifold and the energy surface volume it encloses do not reach $\Sigma$,
which leads to errors of qualitative and quantitative nature. We remark that
one can always estimate the likelihood of tangent intersections by considering
the Hamiltonian vector field on the 3-sphere, for example how close $p_{1}$ is
to $0$ in case of $\Sigma$.

##  V Approximating Reactive Islands With Lagrangian Descriptors

To approximate invariant manifolds, we employ the method of Lagrangian
descriptors (LDs). This technique has its origins in fluid dynamics, more
precisely in Lagrangian transport phenomena [36], and has since become a
popular tool for exploring the phase space of dynamical systems. For a
detailed exposition on the method, with examples on how it can be used to
reveal the phase space of Hamiltonian systems, refer to [48].

A Lagrangian descriptor is a non-negative function $f$ of the phase space
variables that is integrated along trajectories of a system for a fixed amount
of time $\tau$ forward and backward, assigning every initial condition
$(\mathbf{q}_{0},\mathbf{p}_{0})$ at time $t=t_{0}$ a value given by:

$L(\mathbf{q}_{0},\mathbf{p}_{0},t_{0},\tau)=\int^{t_{0}+\tau}_{t_{0}-\tau}f(\mathbf{q},\mathbf{p};\mathbf{q}_{0},\mathbf{p}_{0})\,dt\,.$
(17)

Among the common choices for the function $f$ are arclength, $L^{p}$ norm
($p\geq 1$) or semi-norm ($p<1$) and action [36, 37, 38]. If $f$ is chosen
appropriately, initial conditions leading to qualitatively different dynamics,
for example oscillation and rotation, attain different ranges of values for
the Lagrangian descriptor. Abrupt transitions in the values of $L$ correspond
to (codimension-$1$) invariant manifolds that form the boundary of areas of
qualitatively different dynamics.

Here we follow [45] to directly approximate specific invariant manifolds using
some knowledge of NHIMs in the systems. We are interested in invariant
manifolds of codimension-$1$ in the energy surface $\mathcal{H}=E$, that are
asymptotic to $\mathcal{N}_{E}$, which is of codimension-$2$ in the energy
surface and codimension-$3$ in phase space. Any codimension-$3$ object can be
written as the intersection of three hypersurfaces, namely the energy surface
and two other distinct hypersurfaces $f_{1}(q,p)=0$ and $f_{2}(q,p)=0$. For
example, in the uncoupled system in Sec. II, these can be $q_{1}=0$ and
$p_{1}=0$. Due to the invariance of $\mathcal{N}_{E}$, for any initial
condition $(q,p)$ on $\mathcal{N}_{E}$ we have that:

$\int^{\pm\infty}_{0}|\dot{f_{i}}(q,p)|\,dt=0\;,\quad i=1,2\;.$ (18)

For initial conditions on the stable invariant manifold
$\mathcal{W}^{s}(\mathcal{N}_{E})$ (unstable invariant manifold
$\mathcal{W}^{u}(\mathcal{N}_{E})$) $f_{i}\rightarrow 0$ as
$t\rightarrow\infty$ ($t\rightarrow-\infty$) and the integral in Eq. (18) does
not vanish, but attains a local minimum. Provided the integration time $\tau$
is sufficiently large, a local minimum is attained for finite values. For
having the most accurate result, it is important to select initial conditions
in such a way, that trajectories on the manifolds converge to
$\mathcal{N}_{E}$ without any preceding complicated evolution. In the system
with $2$ DoF considered in [45], accurate results were obtained this way when
$\tau$ was similar to the period of the periodic orbits. We remark that it
suffices to integrate with lower precision than that used to approximate NHIMs
in the system.

To find points on the invariant manifolds for the uncoupled double well
system, see Fig. 2, we uniformly sample the surface $\Sigma_{0.5}$. Since
$\mathcal{W}^{s}(\mathcal{N}_{E})$ and $\mathcal{W}^{u}(\mathcal{N}_{E})$
coincide, we only need to approximate one of them. As argued above, it is
important that trajectories immediately approach $\mathcal{N}_{E}$, we
therefore calculate Lagrangian descriptor values for the stable manifold.

For the coupled double well with $\varepsilon=0.1$ considered here,
$\mathcal{W}^{u}$ does not intersect $\Sigma_{0.5}$ transversely, see Sec. IV
and Fig. 3. For simplicity, in the coupled system we always use the hyperplane

$\Sigma_{0.7}=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\mathcal{H}(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})=E\,,\,q_{1}=-0.7\,,\,\dot{q}_{1}>0\right\\}$
(19)

In general, a surface constructed using the stable object associated with the
local minimum at $q_{1}=\pm\sqrt{1-2\varepsilon}$ will intersect
$\mathcal{W}^{s,u}$ transversely. From $\Sigma_{0.7}$, we integrate
trajectories until $q_{1}=-0.5,\,\dot{q}_{1}>0$ in forward time and
$q_{1}=-0.5,\,\dot{q}_{1}<0$ in backward time and subsequently we calculate
Lagrangian descriptor values for a fixed time interval starting at
$q_{1}=-0.5$. This way we omit the contribution of the dynamics near the local
minimum from Lagrangian descriptor values, making them comparable to values
obtained in the uncoupled system. We remark that higher order reactive islands
can be calculated analogously using trajectories that return to $\Sigma_{0.7}$
the corresponding number of times.

Points closer to the invariant manifolds attain lower LD values. In the
uncoupled system, it is possible to explicitly calculate the values of LDs
along invariant manifolds, we then only considered points whose LD did not
exceed, for example, twice the actual value. In general a suitable cut-off
value can be determined using the lowest attained LD value or from a histogram
of calculated LD values.

The resulting cloud of points approximates the intersection of the invariant
manifold with $\Sigma_{0.5}$ or $\Sigma_{0.7}$, a $3$-sphere, and we perform
calculations on its triangulation. For a more understandable presentation, we
cut the sphere in (disjoint) slices of the form $p_{3}=\text{const}$,
resulting in a continuum of $2$-spheres embedded in $3$-dimensional space. We
triangulate each $2$-sphere separately using the function
`geometry.TriangleMesh.create_from_point_cloud_poisson` of `open3d` [55],
which reconstructs a smooth surface approximated by oriented triangles as
shown in Fig. 3. This Poisson reconstruction implementation creates a convex
hull, with the parameter _depth_ balancing smoothness (low values) and level
of detail (high value). Throughout this work we use $depth=5$, where the
default is value $8$. The resulting mesh may be simplified to an arbitrary
number of triangles. We verify the connectedness of the resulting triangulated
surface using the built in method of mesh objects `is_watertight` and its
Euler characteristic ($2$ for a sphere) using `euler_poincare_characteristic`.
Software for triangulation of higher-dimensional manifolds exists or is in
development and should be used for treatment of higher-dimensional systems.

##  VI Calculating flux

If we select the Poincaré section:

$\Sigma=\left\\{(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})\in\mathbb{R}^{6}\;\Big{|}\;\mathcal{H}(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})=E\,,\,q_{1}=\text{const}\,,\,\dot{q}_{1}>0\right\\}\,,$
(20)

where $p_{1}$ is implicitly given by energy, the volume of the
(unidirectional) flux of energy surface volume enclosed by $\mathcal{W}^{s}$
across a section $\Sigma$ is defined by

$\int_{\mathcal{W}^{s}\cap\Sigma}p_{2}\;dq_{2}\wedge dp_{3}\wedge dq_{3},$
(21)

and analogously for $\mathcal{W}^{u}$. We remark that his flux is a discrete
time flux associated with the return map, that is a volume per iteration of a
return map, which is different to the volume per unit of time $t$.

When cut in finitely many $p_{3}=\text{const}$ slices, integral (21) is
approximated by

$\sum_{i}(p_{3_{i}}-p_{3_{i-1}})\,\int_{S_{i}}p_{2}\,dq_{2}\wedge dq_{3}\,,$
(22)

where $S_{i}=\mathcal{W}^{s}\cap\Sigma\cap\\{p_{3}=p_{3_{i}}\\}$ and upon
triangulation of $S_{i}$ this becomes

$\sum_{i}(p_{3_{i}}-p_{3_{i-1}})\sum\limits_{\Delta}\int\limits_{\Delta}p_{2}\,dq_{2}\wedge
dq_{3}.$ (23)

Given the planarity of individual triangles and adjustable triangulation
density, approximating $p_{2}$ on the triangle by a constant value
$\bar{p}_{2}$, we can further simplify the integral

$\int\limits_{\Delta}p_{2}\,dq_{2}\wedge
dq_{3}\approx\bar{p}_{2}\int\limits_{\Delta}dq_{2}\wedge dq_{3},$ (24)

where $\int\limits_{\Delta}dq_{2}\wedge dq_{3}$ is the area of the projection
of the triangle onto the $q_{2},q_{3}$-plane.

For the uncoupled double well, the integral in Eq. (21) is equal to

$\int_{\mathcal{N}_{E}}dp_{2}\,\wedge dq_{2}\wedge dp_{3}\wedge
dq_{3}=\text{vol}\left(\mathcal{N}_{E}\right)=\dfrac{2\pi^{2}E^{2}}{\omega_{2}\omega_{3}},$
(25)

which for $E=0.1$ and $\lambda=\omega_{2}=\omega_{3}=1$ is approximately
$0.1974$. Table 1 summarises the flux across the first reactive island
estimated by our algorithm for a $60\times 60\times 60\times 10$ uniform grid
in $\\{0\leq q_{2},p_{2},q_{3},p_{3}\leq c_{0.5}\\}$ (negative values are
obtained via symmetry) on $\Sigma_{0.5}$, where $c_{0.5}=\sqrt{2E+7/32}$ is
the maximal value attainable by any of the harmonic oscillators for $q_{1}=\pm
0.5$. Following the Poincaré-Cartan theorem [21] the flux across the reactive
island is the same as (25). We used the Lagrangian descriptor

$\int^{\tau}_{0}|\dot{q}_{1}|\,dt\,,$ (26)

where $\tau$ is chosen close to the period $2\pi$ of orbits on the NHIM.

| cut-off
---|---
|  | 0.55 | 0.6 | 0.8 | 1.0 | 1.5
$\tau$ | 5 | 0.1982 | 0.1983 | 0.1983 | 0.2006 | 0.2373
6 | 0.1976 | 0.1976 | 0.1977 | 0.1977 | 0.1989
7 | - | - | 0.1972 | 0.1972 | 0.1972
| 8 | - | - | 0.1960 | 0.1971 | 0.1972
Table 1: Estimated flux across the first reactive island in the uncoupled
system using LD in Eq. (26) for selected $\tau$ and cut-off. “-” stands for
unreliable estimate due to holes in the triangulated surface and deviating
Euler characteristics.

##  VII Flux in the coupled double well

In the uncoupled system, $\mathcal{W}^{s}$ and $\mathcal{W}^{u}$ coincide,
leading all enclosed trajectories back and forth between the wells. In the
coupled system, $\mathcal{W}^{s}$ and $\mathcal{W}^{u}$ do not coincide (see
Fig. 4). The analogous structure for systems with 2 DoF has been described in
[56] as a partial barrier, which leads some trajectories from one well to the
other, while trapping other trajectories in the wells.

Figure 4: The stable $\mathcal{W}^{s}$ (green) and unstable $\mathcal{W}^{u}$
(purple) manifolds in the coupled system no longer coincide.

Figure 4 shows a single $p_{3}=\text{const}$ slice of the first forward and
backward reactive islands on $\Sigma_{0.7}$. Recall that all trajectories
inside the first forward reactive island (enclosed by $\mathcal{W}^{s}$) are
carried across the dividing surface $\mathcal{D}_{E}$ defined in Eq. (16) in
forward time without returning to $\Sigma_{0.7}$. This includes trajectories
on $\mathcal{W}^{u}$. Following [57], we can determine which part of the first
forward reactive island lies inside the first forward reactive island by
integrating trajectories on $\mathcal{W}^{u}$ forward in time - ones that
reach $\mathcal{D}_{E}$ before $\Sigma_{0.7}$ are inside, ones that reach
$\Sigma_{0.7}$ first are outside. To minimise computational cost, crossings of
$\mathcal{D}_{E}$ and $\Sigma_{0.7}$ are monitored during the Lagrangian
descriptor calculation. Should the integration time be too short for a
significant number of trajectories to reach either surface, we suggest
considering these as lying inside the reactive island.

We proceed in the same way for trajectories on $\mathcal{W}^{s}$ in backward
time. This allows us to identify pieces of $\mathcal{W}^{s}$ and
$\mathcal{W}^{u}$ that bound the intersection of the first forward reactive
island with the first backward reactive island. We remark that the union of
the two reactive islands can be determined analogously. The result is
displayed in Fig. 5. We approximate reactive islands using:

$\int^{\tau}_{0}\Big{|}\dot{q}_{1}-\dfrac{\varepsilon}{c}(\dot{q}_{2}+\dot{q}_{3})\Big{|}\,dt\,,$
(27)

where $c$ is defined in Eq. (14), and a 4-dimensional grid of size $60\times
60\times 60\times 10$ that covers the region

$\\{|q_{2}|,|p_{2}|,|q_{3}|,|p_{3}|\leq c_{0.7}\\},$ (28)

on $\Sigma_{0.7}$, where

$c_{0.7}=\sqrt{2E-\dfrac{(0.7)^{4}}{2}+\left(1-2\varepsilon\right)\cdot(0.7)^{2}}$
(29)

is the maximal value attainable by any of the harmonic oscillators for
$q_{1}=\pm 0.7$. The estimated fluxes are shown in Tab. 2 and 3, which we
compare to a Monte-Carlo estimate of $0.1779$ for the flux across either
reactive island and $0.1102$ for their intersection obtained using $10^{6}$
trajectories launched from $\mathcal{D}_{E}$. It is possible to estimate the
fluxes on a coarser grid, we were able to obtain estimates within $\pm 10\%$
on a uniform $30\times 30\times 30\times 8$ grid, but we consider these
unreliable due to holes in the triangulated surface and deviating Euler
characteristics.

Figure 5: 4D section $q_{1}=-0.5,p_{1}>0$ of the intersection of $\mathcal{W}^{s}$ and $\mathcal{W}^{u}$ on the energy surface $E=0.1$ consisting of several $p_{3}=\text{const}>0$ slices ($p_{3}<0$ is related by symmetry). The parameters are $\omega_{2}=\omega_{3}=1$ and $\varepsilon=0.1$. | cut-off
---|---
|  | 0.6 | 0.8 | 1.0 | 1.1 | 1.2
$\tau$ | 5 | 0.1820 | 0.1843 | - | - | -
6 | 0.1813 | 0.1803 | 0.1811 | - | -
7 | - | 0.1807 | 0.1799 | 0.1797 | -
| 8 | - | - | 0.1814 | 0.1809 | -
Table 2: Estimated flux across the first forward reactive island in the coupled system ($\varepsilon=0.1$) using LD (27) for selected $\tau$ and cut-off. Estimates for flux across the first backward reactive island are within $10^{-3}$. “-” stands for unreliable estimate due to holes in the triangulated surface and deviating Euler characteristics. | cut-off
---|---
|  | 0.6 | 0.8 | 1.0 | 1.1 | 1.2
$\tau$ | 5 | 0.1147 | 0.1155 | - | - | -
6 | 0.1152 | 0.1133 | 0.1144 | - | -
7 | - | 0.1139 | 0.1131 | 0.1130 | -
| 8 | - | - | 0.1146 | 0.1132 | -
Table 3: Estimated flux across the intersection of the first forward and
backward reactive islands in the coupled system ($\varepsilon=0.1$) using LD
(27) for selected $\tau$ and cut-off. “-” stands for unreliable estimate due
to holes in the triangulated surface and deviating Euler characteristics.

When the Poisson reconstruction returns a spherical surface without holes, the
agreement with the Monte-Carlo estimate is within $5\%$ and our results are
robust under parameter variations. This depends mostly on the density of the
grid, as $16$ times as many points (via symmetries) in the uncoupled system
returned an accurate results in a broader range of parameters.

On top of the flux estimates we have carried out, our approach allows for a
qualitative investigation of the intersection of $\mathcal{W}^{s}$ and
$\mathcal{W}^{u}$, which is invaluable for the study of transport and
bifurcation analysis.

##  VIII Conclusions and future work

In this paper we have developed the geometrical, analytical, and computational
framework for reactive island theory for three degrees-of-freedom time-
independent Hamiltonian systems. In the past, reactive island theory has been
successfully applied to studies of chemical reaction dynamics from the phase
space point of view in 2 DoF systems, including the computation of rates in
chemical reactions. This work advances reactive island theory to 3 DoF
systems, providing a foundational step towards achieving this goal for high-
dimensional systems with 3 or more DoF.

For 2 DoF systems, the dividing surface anchored at the unstable periodic
orbit provides the frontier between two dynamically distinct regions (e.g.
reactants and products). The 2-dimensional stable and unstable invariant
manifolds of the unstable periodic orbit have the geometry of two-dimensional
cylinders (or “tubes”) and they act as phase space conduits that determine
whether or not trajectories cross between these dynamically distinct regions.
The intersections of these manifolds with an adequately selected isoenergetic
two-dimensional Poincaré section transverse to the tubes gives rise to the
reactive island picture. From this approach, one can calculate the flux of
trajectories from one region to the other by means of Stokes theorem.

The geometrical structure of higher-dimensional systems based on normally
hyperbolic invariant manifolds and their stable and unstable invariant
manifolds has been known for some time, but required a substantial conceptual
leap in the computational realisation of these structures. Reactive islands
theory in 3 DoF systems is a departure from 2-dimensional Poincaré sections
and traditional ways of finding stable and unstable invariant manifolds.
Instead, recent extensions of the method of Lagrangian descriptors provide the
ability to computationally realize and probe these high-dimensional structures
using low-dimensional representations. Using these advances, we have shown how
this approach can be used to realize a reactive island theory for the
benchmark example of a two-well potential energy surface defining a 3 DoF
Hamiltonian system.

While the visualisation of the reactive islands and the triangulation
procedure we employed only works for 3 DoF systems, these limitations can be
overcome to treat higher dimensional systems. Reliance on the visual
inspection of high-dimensional structures can and will have to be replaced by
numerical approaches, such as the tests of connectedness and Euler
characteristic presented in this work. Triangulation procedures for higher
dimensional manifolds exist or are being developed by computational geometers.
Therefore our framework can be used in higher dimensional systems and only
relies on the knowledge of the relevant normally hyperbolic invariant
manifolds.

Our work provides an entry point into high-dimensional studies of the geometry
of phase space reaction dynamics. Future work might involve the development of
kinetic reaction rate models along the lines of earlier work on reactive
island theory as well as studies of realistic molecular systems in higher
dimensions.

## Acknowledgments

The authors would like to acknowledge the financial support provided by the
EPSRC Grant No. EP/P021123/1 and the Office of Naval Research Grant No.
N00014-01-1-0769.

## References

  * Koon  _et al._ [2000] W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, Chaos 10, 427 (2000).
  * Gómez  _et al._ [2001] G. Gómez, W. S. Koon, M. W. Lo, J. E. Marsden, J. Masdemont, and S. D. Ross,  _Invariant manifolds, the spatial three-body problem and space mission design_ , 109 (American Astronautical Society, 2001).
  * Zhong  _et al._ [2018] J. Zhong, L. N. Virgin, and S. D. Ross, International Journal of Mechanical Sciences 149, 413 (2018).
  * Ross  _et al._ [2018] S. D. Ross, A. E. BozorgMagham, S. Naik, and L. N. Virgin, Phys Rev E 98, 052214 (2018).
  * Zhong and Ross [2020] J. Zhong and S. D. Ross, Communications in Nonlinear Science and Numerical Simulation 82, 105033 (2020).
  * Agaoglou  _et al._ [2019] M. Agaoglou, B. Aguilar-Sanjuan, V. J. García-Garrido, R. García-Meseguer, F. González-Montoya, M. Katsanikas, V. Krajňák, S. Naik, and S. Wiggins,  _Chemical Reactions: A Journey into Phase Space_ (Zenodo, 2019).
  * Conley [1969] C. C. Conley, Journal of Differential Equations 5, 136 (1969).
  * Conley [1968] C. C. Conley, SIAM Journal on Applied Mathematics 16, 732 (1968).
  * Moser [1958] J. Moser, Communications on Pure and Applied Mathematics 11, 257 (1958).
  * Kelley [1967] A. Kelley, Journal of mathematical analysis and applications 18, 472 (1967).
  * Pollak and Pechukas [1978] E. Pollak and P. Pechukas, The Journal of Chemical Physics 69, 1218 (1978).
  * Pechukas and Pollak [1977] P. Pechukas and E. Pollak, The Journal of Chemical Physics 67, 5976 (1977).
  * Pechukas and Pollak [1979] P. Pechukas and E. Pollak, The Journal of Chemical Physics 71, 2062 (1979).
  * Pollak and Child [1980] E. Pollak and M. S. Child, J. Chem. Phys. 73, 4373 (1980).
  * Pollak  _et al._ [1980] P. Pollak, M. S. Child, and P. Pechukas, J. Chem. Phys. 72, 1669 (1980).
  * De Vogelaere and Boudart [1955] R. De Vogelaere and M. Boudart, J. Chem. Phys. 23, 1236 (1955).
  * Marston and De Leon [1989] C. C. Marston and N. De Leon, The Journal of Chemical Physics 91, 3392 (1989).
  * Ozorio de Almeida  _et al._ [1990] A. M. Ozorio de Almeida, N. De Leon, M. A. Mehta, and C. C. Marston, Physica D 46, 265 (1990).
  * De Leon and Marston [1989] N. De Leon and C. C. Marston, The Journal of Chemical Physics 91, 3405 (1989).
  * Mackay [1990] R. Mackay, Physics Letters A 145, 425 (1990).
  * Arnold [1978] V. I. Arnold,  _Mathematical Methods of Classical Mechanics_ (Springer-Verlag New York, 1978).
  * De Leon  _et al._ [1991a] N. De Leon, M. A. Mehta, and R. Q. Topper, The Journal of Chemical Physics 94, 8310 (1991a).
  * De Leon  _et al._ [1991b] N. De Leon, M. A. Mehta, and R. Q. Topper, The Journal of Chemical Physics 94, 8329 (1991b).
  * De Leon [1992] N. De Leon, The Journal of Chemical Physics 96, 285 (1992).
  * Mehta [1990] M. A. Mehta,  _Classical and quantum dynamics of phase space cylindrical manifolds_ , Ph.D. thesis, Yale University (1990).
  * Fair  _et al._ [1995] J. R. Fair, K. R. Wright, and J. S. Hutchinson, The Journal of Physical Chemistry 99, 14707 (1995).
  * Wiggins [1990] S. Wiggins, Physica D 44, 471 (1990).
  * Komatsuzaki and Nagoaka [1997] T. Komatsuzaki and M. Nagoaka, Chem. Phys. Lett. 265, 91 (1997).
  * Komatsuzaki and Berry [2000] T. Komatsuzaki and R. S. Berry, J. Mol. Struct. 506, 55 (2000).
  * Wiggins  _et al._ [2001] S. Wiggins, L. Wiesenfeld, C. Jaffé, and T. Uzer, Phys Rev Lett 86, 5478 (2001).
  * Uzer  _et al._ [2002] T. Uzer, C. Jaffé, J. Palacián, P. Yanguas, and S. Wiggins, Nonlinearity 15, 957 (2002).
  * Komatsuzaki and Berry [2002] T. Komatsuzaki and R. S. Berry, J. Phys. Chem. A 106, 10945 (2002).
  * Wiggins [2016] S. Wiggins, Regul Chaotic Dyn 21, 621 (2016).
  * Patsis and Zachilas [1994] P. Patsis and L. Zachilas, Int. J. Bifurcat. Chaos 04, 1399 (1994).
  * Lange  _et al._ [2014] S. Lange, M. Richter, F. Onken, A. Bäcker, and R. Ketzmerick, Chaos 24, 024409 (2014).
  * Jiménez Madrid and Mancho [2009] J. A. Jiménez Madrid and A. M. Mancho, Chaos 19, 013111 (2009).
  * Mancho  _et al._ [2013] A. M. Mancho, S. Wiggins, J. Curbelo, and C. Mendoza, Communications in Nonlinear Science and Numerical Simulation 18, 3530 (2013).
  * Lopesino  _et al._ [2017] C. Lopesino, F. Balibrea-Iniesta, V. J. García-Garrido, S. Wiggins, and A. M. Mancho, Int J Bifurc Chaos 27, 1730001 (2017).
  * Craven and Hernandez [2015] G. T. Craven and R. Hernandez, Phys Rev Lett 115, 148301 (2015).
  * Craven and Hernandez [2016] G. T. Craven and R. Hernandez, Phys Chem Chem Phys 18, 4008 (2016).
  * Craven  _et al._ [2017] G. T. Craven, A. Junginger, and R. Hernandez, Phys Rev E 96, 022222 (2017).
  * Junginger  _et al._ [2017] A. Junginger, L. Duvenbeck, M. Feldmaier, J. Main, G. Wunner, and R. Hernandez, J. Chem. Phys. 147, 064101 (2017).
  * Feldmaier  _et al._ [2017] M. Feldmaier, A. Junginger, J. Main, G. Wunner, and R. Hernandez, Chem. Phys. Lett. 687, 194 (2017).
  * Revuelta  _et al._ [2019] F. Revuelta, R. M. Benito, and F. Borondo, Phys Rev E 99, 032221 (2019).
  * Krajňák  _et al._ [2020] V. Krajňák, G. S. Ezra, and S. Wiggins, Int. J. Bifurcat. Chaos 30, 2050076 (2020).
  * Montoya and Wiggins [2020a] F. G. Montoya and S. Wiggins, Journal of Physics A: Mathematical and Theoretical 53, 235702 (2020a).
  * Montoya and Wiggins [2020b] F. G. Montoya and S. Wiggins, Physical Review E 102, 062203 (2020b).
  * 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, and S. Wiggins,  _Lagrangian Descriptors: Discovery and Quantification of Phase Space Structure and Transport_ (zenodo: 10.5281/zenodo.3958985, 2020).
  * Naik  _et al._ [2019] S. Naik, V. J. García-Garrido, and S. Wiggins, Communications in Nonlinear Science and Numerical Simulation 79, 104907 (2019).
  * Naik and Wiggins [2019] S. Naik and S. Wiggins, Phys Rev E 100, 022204 (2019).
  * Patra and Keshavamurthy [2015] S. Patra and S. Keshavamurthy, Chemical Physics Letters 634, 1 (2015).
  * Patra and Keshavamurthy [2018] S. Patra and S. Keshavamurthy, Phys. Chem. Chem. Phys. 20, 4970 (2018).
  * Naik and Wiggins [2020] S. Naik and S. Wiggins, Phys. Chem. Chem. Phys. 22, 17890 (2020).
  * Wilson  _et al._ [1980] E. Wilson, J. Decius, and P. Cross,  _Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra_ , Dover Books on Chemistry Series (Dover Publications, 1980).
  * Zhou  _et al._ [2018] Q.-Y. Zhou, J. Park, and V. Koltun, arXiv:1801.09847 (2018).
  * MacKay  _et al._ [1984] R. S. MacKay, J. D. Meiss, and I. C. Percival, Physica D 13, 55 (1984).
  * Krajňák and Waalkens [2020] V. Krajňák and H. Waalkens, Journal of mathematical chemistry 58, 292–339 (2020).