Finding NHIM: Identifying High Dimensional Phase Space Structures in Reaction Dynamics using Lagrangian Descriptors

TL;DR

This paper evaluates the effectiveness of Lagrangian descriptors in identifying high-dimensional phase space structures, such as NHIMs, in Hamiltonian systems, by comparing computational results with known analytical solutions.

Contribution

The study demonstrates that Lagrangian descriptors can accurately reveal high-dimensional phase space structures in Hamiltonian systems, validated against known analytical solutions.

Findings

Lagrangian descriptors successfully identify phase space structures.

The method provides clear visualization of high-dimensional structures.

Validation against analytical solutions confirms accuracy.

Abstract

Phase space structures such as dividing surfaces, normally hyperbolic invariant manifolds, their stable and unstable manifolds have been an integral part of computing quantitative results such as transition fraction, stability erosion in multi-stable mechanical systems, and reaction rates in chemical reaction dynamics. Thus, methods that can reveal their geometry in high dimensional phase space (4 or more dimensions) need to be benchmarked by comparing with known results. In this study, we assess the capability of one such method called Lagrangian descriptor for revealing the types of high dimensional phase space structures associated with index-1 saddle in Hamiltonian systems. The Lagrangian descriptor based approach is applied to two and three degree-of-freedom quadratic Hamiltonian systems where the high dimensional phase space structures are known, that is as closed-form analytical expressions. This leads to a direct comparison of features in the Lagrangian descriptor plots and the phase space structures' intersection with an isoenergetic two-dimensional surface and hence provides a validation of the approach.

Figures

Equations

$$\frac{d\mathbf{x}}{dt}=\mathbf{v}(\mathbf{x},t),\mathbf{x}\in {\mathbb{R}}^n,t\in \mathbb{R}$$

$$M_p({\mathbf{x}}_0,t_0,\tau )={\int }_{t_0-\tau }^{t_0+\tau }\sum _{i=1}^n|{\dot{x}}_i(t;{\mathbf{x}}_0){|}^{p}dt$$

$${\mathcal{W}}^s({\mathbf{x}}_0,t_0)$$

$${\mathcal{W}}^u({\mathbf{x}}_0,t_0)$$

$$\mathcal{M}({\mathbf{x}}_0,t_0)$$

$$\mathcal{H}(\mathbf{q},\mathbf{p})=\frac{\lambda }{2}\left(p_1^2-q_1^2\right)+\frac{{\omega }_2}{2}\left(p_2^2+q_2^2\right)+\frac{{\omega }_3}{2}\left(p_3^2+q_3^2\right)$$

$$\mathcal{H}=\sum _{i=1}^3{\mathcal{H}}_i$$

$${\mathcal{H}}_i=\frac{{\omega }_i}{2}\left(p_i^2+q_i^2\right),i\in \{2,3\},{\mathcal{H}}_1=\frac{\lambda }{2}\left(p_1^2-q_1^2\right)$$

$$\begin{array}{rl}{q}_1(t)& =\frac{1}{2}\left[A{e}^{\lambda t}+B{e}^{-\lambda t}\right]\\ p_1(t)& =\frac{1}{2}\left[A{e}^{\lambda t}-B{e}^{-\lambda t}\right]\\ q_i(t)& =q_i^0\cos ({\omega }_{i}t)+p_i^0\sin ({\omega }_{i}t),i\in \{2,3\}\\ p_i(t)& =p_i^0\cos ({\omega }_{i}t)-q_i^0\sin ({\omega }_{i}t),i\in \{2,3\}\end{array}$$

$$A=q_1^0+p_1^0,B=q_1^0-p_1^0$$

$$\mathcal{H}\left({\mathbf{x}}_0\right)={\mathcal{H}}_0=\sum _{i=1}^3{\mathcal{H}}_i^0$$

$${\mathcal{H}}_i^0=\frac{1}{2}\left(p_i^2+{\omega }_i^2{q}_i^2\right),i\in \{1,\dots ,n-1\},{\mathcal{H}}_n^0=\frac{1}{2}\left(p_n^2-{\lambda }^2{q}_n^2\right)$$

$$q_i^2+p_i^2=\frac{2{\mathcal{H}}_i^0}{{\omega }_i}$$

$$\text{NHIM}=\{\left(\mathbf{q},\mathbf{p}\right)\in {\mathbb{R}}^6|{\mathcal{H}}_0=\sum _{i=2}^3\frac{{\omega }_i}{2}\left(p_i^2+q_i^2\right),q_1=p_1=0\}$$

$$\begin{array}{rl}{\mathcal{W}}^u\left(\text{NHIM}\right)& =\{\left(\mathbf{q},\mathbf{p}\right)\in {\mathbb{R}}^6|{\mathcal{H}}_0=\sum _{i=2}^3\frac{{\omega }_i}{2}\left(p_i^2+q_i^2\right),q_1=p_1\}\\ {\mathcal{W}}^s\left(\text{NHIM}\right)& =\{\left(\mathbf{q},\mathbf{p}\right)\in {\mathbb{R}}^6|{\mathcal{H}}_0=\sum _{i=2}^3\frac{{\omega }_i}{2}\left(p_i^2+q_i^2\right),q_1=-p_1\}\end{array}$$

$$M_{\gamma }({\mathbf{x}}_0,t_0,\tau )={\int }_{t_0-\tau }^{t_0+\tau }\sum _{i=1}^3|{\dot{q}}_i(t;{\mathbf{x}}_0){|}^{\gamma }+|{\dot{p}}_i(t;{\mathbf{x}}_0){|}^{\gamma }dt$$

$$M_{\gamma }({\mathbf{x}}_0,t_0,\tau )=M_{\gamma }^h({\mathbf{x}}_0,t_0,\tau )+M_{\gamma }^e({\mathbf{x}}_0,t_0,\tau )$$

$$M_{\gamma }^h({\mathbf{x}}_0,t_0,\tau )={\int }_{t_0-\tau }^{t_0+\tau }|{\dot{q}}_1(t;{\mathbf{x}}_0){|}^{\gamma }+|{\dot{p}}_1(t;{\mathbf{x}}_0){|}^{\gamma }dt$$

$$M_{\gamma }^e({\mathbf{x}}_0,t_0,\tau )={\int }_{t_0-\tau }^{t_0+\tau }\sum _{i=2}^3|{\dot{q}}_i(t;{\mathbf{x}}_0){|}^{\gamma }+|{\dot{p}}_i(t;{\mathbf{x}}_0){|}^{\gamma }dt$$

$$M_{\gamma }^h\left({\mathbf{x}}_0,\tau \right)=M_{\gamma }^{h,q_1}\left({\mathbf{x}}_0,\tau \right)+M_{\gamma }^{h,p_1}\left({\mathbf{x}}_0,\tau \right)={\int }_{-\tau }^{\tau }|{\dot{q}}_1{|}^{\gamma }dt+{\int }_{-\tau }^{\tau }|{\dot{p}}_1{|}^{\gamma }dt$$

$${\dot{q}}_1(t)=\frac{\lambda }{2}\left[A{e}^{\lambda t}-B{e}^{-\lambda t}\right],{\dot{p}}_1(t)=\frac{\lambda }{2}\left[A{e}^{\lambda t}+B{e}^{-\lambda t}\right]$$

$$\begin{array}{rl}{M}_{\gamma }^{h,q_1}({\mathbf{x}}_0,\tau )& ={\int }_{-\tau }^{-{\tau }_0}|{\dot{q}}_1{|}^{\gamma }dt+M_{\gamma }^{h,q_1}({\mathbf{x}}_0,{\tau }_0)+{\int }_{{\tau }_0}^{\tau }|{\dot{q}}_1{|}^{\gamma }dt\\ M_{\gamma }^{h,p_1}({\mathbf{x}}_0,\tau )& ={\int }_{-\tau }^{-{\tau }_0}|{\dot{p}}_1{|}^{\gamma }dt+M_{\gamma }^{h,p_1}({\mathbf{x}}_0,{\tau }_0)+{\int }_{{\tau }_0}^{\tau }|{\dot{p}}_1{|}^{\gamma }dt\end{array}$$

$$\begin{array}{rl}{M}_{\gamma }^{h,q_1}({\mathbf{x}}_0,{\tau }_0)& =M_{\gamma }^{h,q_1}({\mathbf{x}}_0,0)+{\tau }_0\frac{\partial M_{\gamma }^{h,q_1}}{\partial {\tau }_0}({\mathbf{x}}_0,0)+O({\tau }_0^2)=\\ & =2{\tau }_0{\lambda }^{\gamma }|p_1^0{|}^{\gamma }+O({\tau }_0^2)\approx 2{\tau }_0{\lambda }^{\gamma }|p_1^0{|}^{\gamma }\\ M_{\gamma }^{h,p_1}({\mathbf{x}}_0,{\tau }_0)& =M_{\gamma }^{h,p_1}({\mathbf{x}}_0,0)+{\tau }_0\frac{\partial M_{\gamma }^{h,p_1}}{\partial {\tau }_0}({\mathbf{x}}_0,0)+O({\tau }_0^2)=\\ & =2{\tau }_0{\lambda }^{\gamma }|q_1^0{|}^{\gamma }+O({\tau }_0^2)\approx 2{\tau }_0{\lambda }^{\gamma }|q_1^0{|}^{\gamma }\end{array}$$

$$M_{\gamma }^h\left({\mathbf{x}}_0,{\tau }_0\right)\approx 2{\tau }_0{\lambda }^{\gamma }\left(|p_1^0{|}^{\gamma }+|q_1^0{|}^{\gamma }\right)$$

$$\begin{array}{rl}{\int }_{{\tau }_0}^{\tau }|{\dot{q}}_1{|}^{\gamma }dt& =\frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|A{e}^{\lambda t}-B{e}^{-\lambda t}{|}^{\gamma }dt=\frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }{\left|A{e}^{\lambda t}\left(1-\frac{B}{A}{e}^{-2\lambda t}\right)\right|}^{\gamma }dt=\\ & =\frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|A{e}^{\lambda t}{|}^{\gamma }{\left|1-\frac{B}{A}{e}^{-2\lambda t}\right|}^{\gamma }dt=\frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|A{e}^{\lambda t}{|}^{\gamma }+O\left(\frac{|B|}{|A{|}^{1-p}}{e}^{-\gamma \lambda t}\right)dt\end{array}$$

$${\int }_{{\tau }_0}^{\tau }|{\dot{q}}_1{|}^{\gamma }dt\sim \frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|A{e}^{\lambda t}{|}^{\gamma }dt=\frac{|A{|}^{\gamma }{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }{e}^{\gamma \lambda t}dt\sim \frac{|A{|}^{\gamma }{\lambda }^{\gamma -1}}{\gamma 2^{\gamma }}{e}^{\gamma \lambda \tau }$$

$$\begin{array}{rl}{\int }_{-\tau }^{-{\tau }_0}|{\dot{q}}_1{|}^{\gamma }dt& =\frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{-\tau }^{-{\tau }_0}|A{e}^{\lambda t}-B{e}^{-\lambda t}{|}^{\gamma }dt=\frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|B{e}^{\lambda s}{|}^{\gamma }{\left|1-\frac{A}{B}{e}^{-2\lambda s}\right|}^{\gamma }ds\\ & =\frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|B{e}^{\lambda s}{|}^{\gamma }+O\left(\frac{|A|}{|B{|}^{1-p}}{e}^{-\gamma \lambda s}\right)ds\end{array}$$

$${\int }_{-\tau }^{-{\tau }_0}|{\dot{q}}_1{|}^{\gamma }dt\sim \frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|B{e}^{\lambda s}{|}^{\gamma }ds=\frac{|B{|}^{\gamma }{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }{e}^{\gamma \lambda s}ds\sim \frac{|B{|}^{\gamma }{\lambda }^{\gamma -1}}{\gamma 2^{\gamma }}{e}^{\gamma \lambda \tau }$$

$$\begin{array}{rl}{\int }_{{\tau }_0}^{\tau }|{\dot{p}}_1{|}^{\gamma }dt& \sim \frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|A{e}^{\lambda t}{|}^{\gamma }dt\sim \frac{|A{|}^{\gamma }{\lambda }^{\gamma -1}}{\gamma 2^{\gamma }}{e}^{\gamma \lambda \tau }\\ {\int }_{-\tau }^{-{\tau }_0}|{\dot{p}}_1{|}^{\gamma }dt& \sim \frac{{\lambda }^{\gamma }}{2^{\gamma }}{\int }_{{\tau }_0}^{\tau }|B{e}^{\lambda t}{|}^{\gamma }dt\sim \frac{|B{|}^{\gamma }{\lambda }^{\gamma -1}}{\gamma 2^{\gamma }}{e}^{\gamma \lambda \tau }\end{array}$$

$$M_{\gamma }^h\left({\mathbf{x}}_0,\tau \right)\sim \frac{{\lambda }^{\gamma -1}}{\gamma 2^{\gamma -1}}\left(|A{|}^{\gamma }+|B{|}^{\gamma }\right)e^{\gamma \lambda \tau }$$

$$M_{\gamma }^e({\mathbf{x}}_0,\tau )=\sum _{i=2}^3{M}_{\gamma }^{e,i}({\mathbf{x}}_0,\tau )=\sum _{i=2}^3{\int }_{-\tau }^{\tau }|{\dot{q}}_i(t;{\mathbf{x}}_0){|}^{\gamma }+|{\dot{p}}_i(t;{\mathbf{x}}_0){|}^{\gamma }dt$$

Full text

Keywords: Normally hyperbolic invariant manifolds, Stable and unstable manifolds, Lagrangian descriptors, Hamiltonian systems, Chemical reaction dynamics, Phase space transport.

1 Introduction

Identifying invariant manifolds in high dimensional phase space and their anchor, normally hyperbolic invariant manifold (NHIM), is the stepping stone in applying phase space transport methods to a broad array of problems in physics, chemistry, and engineering [45, 22, 24, 8, 9, 17, 18, 44, 43, 46, 47]. These problems are typically formulated in phase space of more than 4 dimensions (that is, $N=2$ or more degrees of freedom) and the geometric approach of computing invariant manifolds that are codimension-1 separatrices requires Poincaré-Birkoff normal form theory. The Poincaré-Birkoff normal form approach for computing NHIMs and their stable and unstable manifolds has only been developed for Hamiltonian systems in the neighborhood of index-k saddle points. This approach in itself is successful and rigorous, but its implementation requires experience in writing or using algebraic manipulation programs. Furthermore, these computations become expensive as the dimensions increase since Poincaré-Birkoff normal form theory is based on Taylor expansion, a high dimensional vector valued polynomial, of the Hamiltonian in a neighborhood of the saddle point to sufficiently high order. Besides, the resulting visualization can be non-obvious leading to cumbersome interpretation and can not work for high dimensions. The succcess and shortcoming of the vizualization approach based on topological methods for four-dimensional space becomes apparent when trying to extend the results in Refs. [30, 2, 23]. This raises the question of detecting the signatures of the high dimensional structures, if they exist, when probed with low dimensional surfaces. If one does manage to detect these structures, how do they manifest on the low dimensional surfaces.

One such method is the Lagrangian descriptors (LDs) developed in Refs. [33, 35, 34, 32]. Lagrangian descriptors are a class of trajectory diagnostic methods that can reveal phase space structures by encoding geometric property of trajectories (such as, phase space arc length, configuration space distance or displacement, cumulative action or kinetic energy) initialised on a two dimensional surface. The method was originally developed in the context of Lagrangian transport studies in fluid dynamics that requires identifying transport barriers which are the invariant manifolds in time dependent 2D fluid flow. Besides, the method applies to both Hamiltonian and non-Hamiltonian systems ([32]) as well as to systems with arbitrary, stochastic and dissipative, time-dependence ([1, 6, 25]). Lagrangian descriptor based detection of transport barriers has also been applied directly to data sets, such as those obtained from satellite observations or numerical simulations ([7, 36, 19, 40]). Furthermore, the method is straightforward to implement computationally and it provides a “high resolution” method for exploring the influence of high dimensional phase space structure on trajectory behaviour. The method of Lagrangian descriptors takes an opposite approach to that of Lyapunov exponent type calculations by emphasizing the initial conditions of trajectories, rather than their advected locations that is involved in calculating normalized rate of divergence. This is achieved by considering a two dimensional section of the full phase space and discretizing with a dense grid of initial conditions. Even though the trajectories wander off in the phase space, as the initial conditions evolve in time, there is no loss in resolution of the two dimensional section. Our objective is to clarify the use of Lagrangian descriptors as a diagnostic on two dimensional sections of high dimensional phase space structures. This diagnostic is also meant to be used as the preliminary step in computing the NHIM, their stable and unstable manifolds using other computational means [27, 3, 13]. In this article, we assess the capability of Lagrangian descriptors to detect the type of high dimensional phase space structures in Hamiltonian systems such as the NHIM, their stable, and unstable manifolds that are used in computing rates of chemical reactions.

Computing chemical reaction rates is based on the fundamental framework of transition state theory as formulated in phase space by Polanyi, Evans, Wigner, Eyring [11, 48, 12, 21]). Further research has shown that constructing a locally recrosssing free, minimal flux orientable surface called a dividing surface (DS) is the phase space structure that provides the correct estimate for reaction rates. This dividing surface of geometry $\mathbb{S}^{2N-2}$ ( a 2N-2 dimensional sphere) is constructed from the NHIM of geometry $\mathbb{S}^{2N-3}$ which forms the equator of the dividing surface, on a constant energy surface [42, 44, 43, 47]. Furthermore, the global dynamics of reactive and non-reactive trajectories is governed by the invariant (stable and unstable) manifolds, of the NHIM, are $\mathbb{S}^{2N-3}\times\mathbb{R}$, and thus act as codimension-1 separatrices on the energy surface. Thus, the NHIM acts as the anchor for the local dynamics via the dividing surface from reactants to products or vice-versa and for the global dynamics via the stable/unstable invariant manifolds. Thus, detecting and constructing the NHIM forms a fundamental step in obtaining quantitative results in reaction dynamics [14, 13]. In the present study, we will focus on detecting and verifying these high dimensional phase space structures.

Recently the applicability of the Lagrangian descriptor based approach to time dependent problems with random and dissipative forcing in chemical reaction dynamics has been shown whereby the transition state trajectory is calculated using the extremal values in Lagrangian descriptor values; see Refs. [4, 5, 6, 26, 27, 28, 29, 15, 41, 39]. The initial conditions for the transition state trajectory is identified by computing the extrema of the Lagrangian descriptor on a two dimensional domain. Comparing extremal and singular features in the Lagrangian descriptor plot with invariant manifolds in three dimensional vector fields, Ref. [20] has also provided numerical evidence for detecting NHIM using the LD on two dimensional surfaces. The authors presented the comparison using a stacked version of a two dimensional linear saddle, the duffing oscillator, and a time-perturbed 3D geophysical model. This has the advantage that the NHIM and the invariant manifolds are known exactly or can be computed with other established methods. For the problems considered therein, the NHIM (a curve embedded in 3D) and its associated invariant manifolds (a 2D surface in 3D) can be visualized in the 3D space since the dimensionality of phase space structures is less than 3. We take this approach to the next logical step by applying this to benchmark problems of two and three degrees-of-freedom Hamiltonian system where the phase space structures are known exactly, that is they can be written as closed-form analytical expressions

This article is outlined as follows. In Section 2, we describe the method of Lagrangian descriptor used in this study and present an analytical result on identifying invariant manifolds using features in the Lagrangian descriptor. In Section 3 and 4, we discuss the benchmark systems and present numerical evidence of our claim using isoenergetic two-dimensional surfaces. In Section 5, we summarize our results on detecting invariant manifolds and discuss related future directions.

2 Theory

2.1 Method of Lagrangian Descriptor

The Lagrangian descriptor (LD) as presented in Ref.[33] is the arc length of a trajectory calculated on a chosen initial time $t_{0}$ and measured for fixed forward and backward integration time, $\tau$. For continuous time dynamical systems, Ref.[32] gives an alternative definition of the LD which is useful for proving rigorous results and can be computed along with the trajectory. It provides a characterization of the notion of singular features of the LD that facilitates a proof for detecting invariant manifolds in certain model situations. In addition, the “additive nature” of this new definition of LD provides an approach for assessing the influence of each degree-of-freedom separately on the Lagrangian descriptor. This property was used in Ref.[10] which showed that a Lagrangian descriptor can be used to detect Lyapunov periodic orbits in the two degrees-of-freedom Hénon-Heiles Hamiltonian system. We will describe this procedure for two and three degrees-of-freedom linear autonomous Hamiltonian systems. We begin by establishing notation in the general setting of a time-dependent vector field where

[TABLE]

where $\mathbf{v}(\mathbf{x},t)\in C^{r}$ ($r\geq 1$) in $\mathbf{x}$ and continuous in time. The definition of LDs depends on the initial condition $\mathbf{x}_{0}=\mathbf{x}(t_{0})$, on the initial time $t_{0}$ (trivial for autonomous systems) and the integration time $\tau$, and the type of norm of the trajectory’s components, and takes the form,

[TABLE]

where $p\in(0,1]$ and $\tau\in\mathbb{R}^{+}$ are freely chosen parameters, and the overdot symbol represents the derivative with respect to time. It is to be noted here that there are three formulations of the function $M_{p}$ in the literature: the arc length of a trajectory in phase space [33], the arc length of a trajectory projected on the configuration space  [27, 26, 28, 29], and the sum of the $p$-norm of the vector field components [31, 32]. Although the latter formulation of the Lagrangian descriptor (2) developed in Ref. [31, 32] does not resemble the arc length, the numerical results using either of these forms have been shown to be in agreement and promise of predictive capability in geophysical flows [7, 36, 19, 40]. The formulation we adopt here is motivated by the fact that this allows for proving rigorous result, which we will discuss in the next section, connecting the singular features and minimum in the LD plots with NHIM and its stable and unstable manifolds. It follows from the result that

[TABLE]

where the stable and unstable manifolds ($\mathcal{W}^{s}(\mathbf{x}_{0},t_{0})$ and $\mathcal{W}^{u}(\mathbf{x}_{0},t_{0})$) denote the invariant manifolds at intial time $t_{0}$ and $\text{\rm argmin}\;(\cdot)$ denotes the argument that minimizes the function $\mathcal{L}^{(\cdot)}(\mathbf{x}_{0},t_{0},\tau)$ in forward and backward time, respectively. In addition, the coordinates of the NHIM at time $t_{0}$ is given by the intersection $\mathcal{W}^{s}(\mathbf{x}_{0},t_{0})$ and $\mathcal{W}^{u}(\mathbf{x}_{0},t_{0})$ of the stable and unstable manifolds, and thus given by

[TABLE]

2.2 Lagrangian descriptor and invariant manifolds

We dedicate this section to prove that the method of Lagrangian descriptor recovers the NHIM and its stable and unstable manifolds for a linear quadratic Hamiltonian system with 3 DoF. We will show that given a sufficiently large integration time, the Lagrangian descriptor value reaches a minimum value at the NHIM and its stable and unstable manifolds, and also that singularities (singular features) in the Lagrangian descriptor values, that is points where the scalar function is non-differentiable, identify the stable and unstable manifolds of the NHIM, and hence also the NHIM. Therefore, this method can be used to detect the NHIM and thus a way to construct the dividing surface that is essential for the computation of reaction rates in chemistry. The arguments we follow to prove this mathematical connection are based on those described in Refs. [32, 10].

Consider the quadratic Hamiltonian with $3$ DoF that has a index-1 saddle at the origin

[TABLE]

Observe that all the DoF that appear in the Hamiltonian $\mathcal{H}$ are uncoupled so we can write:

[TABLE]

where the energy in each DoF is given by

[TABLE]

Given the initial condition $\mathbf{x}_{0}=\mathbf{x}(t_{0})=(\mathbf{q}_{0},\mathbf{p}_{0})\in\mathbb{R}^{6}$ at time $t_{0}=0$, where $\mathbf{q}_{0}=(q_{1}^{0},q_{2}^{0},q_{3}^{0})$ and $\mathbf{p}_{0}=(p_{1}^{0},p_{2}^{0},p_{3}^{0})$, the general solution to Hamilton’s equations is

[TABLE]

where the constants $A$ and $B$ are given by

[TABLE]

This initial condition gives rise to a phase space trajectory with total energy

[TABLE]

which is distributed among all the DoF of the system, and the energy in each DoF is

[TABLE]

Observe that each bath mode $(q_{i},p_{i})$ will give rise to a periodic orbit in the $q_{i}$-$p_{i}$ plane, where $i\in\{2,3\}$. Moreover, this periodic orbit is a cricle

[TABLE]

with radius $R=\sqrt{\dfrac{2\mathcal{H}^{0}_{i}}{\omega_{i}}}$. The isoenergetic NHIM for this Hamiltonian system is:

[TABLE]

which is a $3$-dimensional sphere. The stable and unstable manifolds of the NHIM are

[TABLE]

Take a fixed integration time $\tau>0$ and $\gamma\in(0,1]$, The Lagrangian descriptor is

[TABLE]

Observe that we can decompose this integral into the hyperbolic and elliptic components

[TABLE]

where the hyperbolic part is

[TABLE]

and the elliptic component is given by

[TABLE]

We focus our attention first on the hyperbolic contribution of LDs. Since the dynamical system is autonomous we can take without loss of generality $t_{0}=0$, and we eliminitae the dependence of LDs on the initial time to simplify notation. We can write

[TABLE]

for which we need the time derivative of the position and momentum coordinates

[TABLE]

Observe that it is not possible to calculate the integrals $M_{\gamma}^{h,q_{n}}$ and $M_{\gamma}^{h,p_{n}}$ analytically, so we will accurately approximate their values by means of an asymptotic analysis. Take a small value of $\tau_{0}$, then we can write

[TABLE]

Expanding $M^{h,q_{1}}_{\gamma}(\mathbf{x}_{0},\tau_{0})$ and $M^{h,p_{1}}_{\gamma}(\mathbf{x}_{0},\tau_{0})$ in a Taylor series about $\tau=0$ gives

[TABLE]

This clearly shows that the hyperobolic component of LDs behaves for small $\tau_{0}$ as

[TABLE]

having singularities on the phase space sets $q_{n}^{0}=0$ or $p_{n}^{0}=0$ which do not coincide with the unstable and stable manifolds of the NHIM for the linear Hamiltonian. Nevertheless, as we describe next, when the integration time $\tau$ is sufficiently large, the leading order singularities that appear in LDs align with the true stable and unstalble manifolds of the NHIM.

Take $\tau\gg 1$, we deal first with the integral that involves the time derivative of $q_{1}$

[TABLE]

where we have used the expansion $(1+x)^{p}=1+px^{p-1}+O(x^{2})$ for small $x$. Then the asymptotic behavior of the integral is

[TABLE]

Analogously, for negative times

[TABLE]

Then the asymptotic behavior of the integral is

[TABLE]

If we develop a similar asymptotic argument for the integral that involves the time derivative of $p_{1}$ it can be easily shown that

[TABLE]

In summary, the asymptotic behavior as $\tau\gg 1$ of the hyperbolic component of LDs is

[TABLE]

Therefore, we have show that $M_{\gamma}^{h}$ grows exponentially with $\tau$ and also that the leading order singularities in $M_{\gamma}^{h}$ occur when $|A|=0$, that is, when $p^{0}_{1}=-q^{0}_{1}$, which corresponds to initial conditions on the stable manifold of the NHIM, or in the case where $|B|=0$, that is, $p^{0}_{1}=q^{0}_{1}$, corresponding to initial conditions on the unstable manifold of the NHIM. Moreover, $M_{\gamma}^{h}$ is non-differentiable at the NHIM, since it is given by the intersection of the stable and unstable manifolds.

To finish the proof we focus now on the elliptic component of LDs. We have that

[TABLE]

Since all the bath modes follow the same type of dynamics, we focus on analyzing the contribution of the bath mode $(q_{i},p_{i})$ where $i\in\{2,3\}$ to the elliptic component of LDs, that is, $M^{e,i}_{\gamma}(\mathbf{x}_{0},\tau)$. We recall that

[TABLE]

Then we have to evaluate the integral

[TABLE]

We deal first with the integral

[TABLE]

Using that $|\dot{q}_{i}|$ is periodic with period $T_{i}=\pi/\omega_{i}$, we can always write $\tau=NT_{i}+r$ for some integer $N$ and $r\in[0,T_{i}]$. With this decomposition we can write

[TABLE]

Since the first integral covers one whole period, without loss of generality we can compute its value by starting on the initial condition $p_{i}^{0}=0$ and $q_{i}^{0}=\sqrt{2\mathcal{H}_{i}^{0}/\omega_{i}}$. We can do this because the bath mode $(q_{i},p_{i})$ follows a circle in the $q_{i}$-$p_{i}$ plane with that radius. This yields

[TABLE]

where $B$ is the beta function and $\Gamma$ is the gamma function. Threfore we have shown that

[TABLE]

Followng this argument, we can also compute the integral

[TABLE]

by using the decomposition

[TABLE]

The integral over one period in this case, choosing $q_{i}^{0}=0$ and $p_{i}^{0}=\sqrt{2\mathcal{H}_{i}^{0}/\omega_{i}}$, is

[TABLE]

This gives:

[TABLE]

Consequently, the elliptic component of LDs corresponding to the bath mode $(q_{i},p_{i})$ is

[TABLE]

which clearly shows that it is a linearly increasing function of $\tau$, in contrast to the hyperbolic component of LDs that grows exponentially with the integration time. If we divide the elliptic part of LDs by $2\tau$ (the integration period of the trajectory) and we let $\tau\to\infty$ we get that

[TABLE]

where $R_{i}$ is the radius of the circle described by the bath coordinate $(q_{i},p_{i})$. So the time average of the elliptic part of LDs converges to a value that depends on the energy $\mathcal{H}_{i}^{0}$ of the initial condition $(q_{i}^{0},p_{i}^{0})$, and therefore in the limit its value is constant on the periodic orbit described by the bath mode on the $q_{i}$-$p_{i}$ plane. This result is of interest because it allows us to use LDs as a tool to recover phase space KAM tori, given that the dynamical system under study satisfies the conditions of the Ergodic Partition Theorem [37]. Therefore, if one analyzes thr time average of LDs on a specific initial condition in phase space $\mathbf{x}_{0}$, and its value converges, then this initial copndition would lie in an invariant phase space set cosisting of points that share the same time average value. Therefore, the contours of the time average of LDs, when it converges, identify invariant phase space sets. It is also important to highlight that the limit value to which $M^{e,i}$ converges is directly related to the limit value to which the classical definition of LDs, representing the arclength of the trajectory of the bath mode, converges. Indeed:

[TABLE]

From Eqs. (18) and (19) it follows that

[TABLE]

To conclude the proof we will show how LDs attains a local minimum at the phase space points corresponding to the stable and unstable manifolds of the NHIM and a global minimum at the NHIM. Given an energy of the system $\mathcal{H}_{0}$ above that of the origin, we know that a family of NHIM parametrized by the energy bifurcates from the rank-1 saddle. The phase space points that lie on the stable (or unstable) manifold, that is $q_{1}^{0}=-p_{1}^{0}$ ($q_{1}^{0}=p_{1}^{0}$), contribute to the hyperbolic component of LDs described in Eq. (15) with $|A|=0$ ($|B|=0$) so that the Lagrangian descriptor has a local minimum at the manifold. Moreover, if the initial condition is on the NHIM, that is $q_{1}^{0}=p_{1}^{0}=0$, then $|A|=|B|=0$ and consequently $M^{h}_{\gamma}(\mathbf{x}_{0},\tau)=0$. Furthermore, all the energy of the system for these points concentrates on the bath modes that evolve periodically, which implies that $M^{e}_{\gamma}(\mathbf{x}_{0},\tau)>0$. As a result, Lagrangian descriptors attain a global minimum at these points. We would like to conclude by noting that due to Moser’s theorem (generalization of the Liapunov’s subcenter theorem) on index-1 saddle in a Hamiltonian system, this result on Lagrangian descriptor for a linear system extends to the local neighborhood of an index-1 saddle equilibrium point in the full nonlinear system [38].

3 Two Degrees of Freedom

In this section we use the Lagrangian descriptor method to identify the NHIM and its stable and unstable manifolds for a 2 DoF separable quadratic Hamiltonian with an index-1 saddle at the origin. The advantage of using this model Hamiltonian is that we can compare the analytical expressions for the phase space structures with features in LD plots.

3.1 Decoupled quadratic Hamiltonian

We consider a quadratic Hamiltonian for a two degree-of-freedom system given by

[TABLE]

with the corresponding vector field

[TABLE]

The equilibirum point is located at $(0,0,0,0)$ on the zero total energy surface. It is trivial to check that the eigenvalues of the linearized system around the equilibrium point are $\pm\lambda$ and $\pm i\omega_{2}$, and hence the equilibrium point is of saddle $\times$ center type. In this form, the Hamiltonian (20) is decoupled into the “reactive” mode given by $H_{r}$ and the “bath” mode given by $H_{b}$, hence it will be referred to as separable quadratic Hamiltonian (SQH). This representation lends to discussing the distribution of total energy between the two modes for the phase space structures in uncoupled coordinates. In this form, a chemical reaction is said to have occurred when the $q_{1}$ coordinate of a trajectory changes sign and thus, an isoenergetic, $H(q_{1},p_{1},q_{2},p_{2})=h$, dividing surface (DS) can be defined by $q_{1}=0$ hypersurface. The constant energy defines a three dimensional surface in the four dimensional phase space and is given by

[TABLE]

The dividing surface, $q_{1}=0$, for a constant energy is

[TABLE]

which is a two dimensional surface, or precisely of geometry $\mathbb{S}^{2}$, that is a 2-sphere on the three dimensional energy surface. Thus it is codimension-1 and partitions the energy surface into reactant $p_{1}-q_{1}>0$ and product $p_{1}-q_{1}<0$ regions by the forward and backward “reaction” dividing surfaces as shown in Ref. [42] and given by

[TABLE]

The forward and backward DS is joined at $p_{1}=0$, thus

[TABLE]

which is of geometry $\mathbb{S}^{1}$, that is a circle centered at the origin and of radius $\sqrt{2h/\omega_{2}}$ in the $(q_{2},p_{2})$ plane. This is a normally hyperbolic invariant manifold (NHIM) associated with the index-1 saddle and parametrized by total energy $H(q_{1},p_{1},q_{2},p_{2})=h$ [46]. Invariance follows from the vector field (21) since when $q_{1}=p_{1}=0$, $\dot{q}_{1}=\dot{p}_{1}=0$. Thus $q_{1}$ and $p_{1}$ always remain zero, and trajectories with these initial conditions remain on the NHIM, that is, $q_{1}=p_{1}=0$ is invariant. It is normally hyperbolic since the directions normal to the NHIM, that is, the $(q_{1},p_{1})$ surface, have linear saddle like dynamics. For a two degree-of-freedom system, this NHIM is more commonly referred in the literature as an unstable periodic orbit.

In order to understand the relationship between the NHIM and the index 1 saddle point, we note that for $H_{r}=0$ and $H_{b}=0$ the NHIM reduces to the point $(q_{1},p_{1},q_{2},p_{2})=(0,0,0,0)$, which is the index-1 saddle point on the energy surface $H_{r}+H_{b}=0$. Therefore, as the total energy is increased from [math], with $H_{b}$ increasing from zero, the NHIM “grows” from the index-1 saddle point on the zero energy surface into an invariant 1-sphere. This shows how the “influence” of the index-1 saddle point is carried to higher energy sufaces on which the saddle point does not exist.

The stable and unstable manifolds of the NHIM (25) are given by

[TABLE]

which is two dimensional surface and have geometry $\mathbb{R}\times\mathbb{S}^{1}$ for a fixed energy. Thus, the codimension-1 geometry of the manifolds partition the phase space into “reactive” and “non-reactive” trajectories as shown in Fig. 1. The detailed visualization is only possible for the two degrees-of-freedom system which is also the starting point for testing the Lagrangian descriptor based approach for detecting the phase space structures.

3.1.1 Detecting the unstable periodic orbit and its manifolds

To identify NHIM and its invariant manifolds, we compute the Lagrangian descriptor in a square domain of size 2 units around the origin and discretize the coordinates of the two-dimensional surface. Next, we pick a constant value for one of the two remaining coordinates, and use the total energy equation to solve for the fourth coordinate. Due to the form of the Hamiltonian (20), obtaining the coordinate from the constant energy condition reduces simply to solving a quadratic equation.

Isoenergetic two-dimensional surface parametrized by $(q_{1},p_{1})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{1},p_{1})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{1}p_{1}}^{+}$ with the NHIM (25) becomes

[TABLE]

Thus, the NHIM is located at the origin $(0,0)$ and marked by a red cross in the LD plot (Fig. 2). Furthermore, the intersection of the two-dimensional surface with the unstable (26) and stable manifolds (27) is given by

[TABLE]

which are one-dimensional for a fixed energy, and represent lines passing through the origin shown as dashed red (unstable) and white (stable) lines, respectively, in Fig. 2. The only points of local minima in the LD plot (Fig. 2) also lie along the lines passing through the origin and correspond to the manifolds of the NHIM.

Isoenergetic two-dimensional surface parametrized by $(q_{2},p_{2})$ — Next, on the fixed energy surface, $H(q_{1},p_{1},q_{2},p_{2})=h$, we compute the Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{2},p_{2})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{2}p_{2}}^{+}$ with the NHIM (25) is given by

[TABLE]

which represents a circle of radius $\sqrt{2h/\omega_{2}}$. The radius is $\approx 0.632$ for $h=0.2,\omega_{2}=1.0$ and is marked as a dashdot red circle in Fig. 2. The points on the circle are also locations of minima as shown by the one-dimensional slices at constant $p_{2}$.

Next, intersection of the two-dimensional $U_{q_{2}p_{2}}^{+}$ with the stable (27) and unstable manifolds (26) is given by

[TABLE]

which are circles of same radius $\sqrt{2h/\omega_{2}}$. These manifolds are identified along the boundary of the LD plot in Fig. 2 along which the minima also occurs.

Isoenergetic two-dimensinal surface parametrized by $(q_{1},q_{2})$ — Next, on the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2})=h$, we compute the Lagrangian descriptor on a two-dimensional surface by defining

[TABLE]

where

[TABLE]

The intersection of the NHIM (25) with the two-dimensional surface $U_{q_{1}q_{2}}^{+}$ is

[TABLE]

which represent points on the line $q_{1}=0$ which is marked by the dashdot line in Fig. 2.

Next, intersection of the two-dimensional $U_{q_{1}q_{2}}^{+}$ with the stable (27) and unstable manifolds (26) is given by

[TABLE]

which denote lines parallel to the $q_{2}$ axis and marked by the dashed red (unstable) and white (stable) lines in Fig. 2. These manifolds are identified by minima in the Lagrangian descriptor plot in Fig. 2 as also shown in the one-dimensional slice for constant $q_{2}$.

Now the question is how do these phase space structures manifest when the system is non-separable into “reactant” and “product” coordinates. We will answer this using a linear symplectic transformation of the separable quadratic Hamiltonian and then compute the LD on different 2D surfaces for a fixed energy.

3.2 Coupled quadratic Hamiltonian

To couple the coordinates in the separable quadratic Hamiltonian (20), we introduce a linear transformation, $C$, such that it satisfies the symplectic condition

[TABLE]

where $\mathcal{J}$ is the 2N $\times$ 2N matrix and $I_{N}$ denotes the $N\times N$ identity matrix. The symplectic transformation $C$ acts on the non-separable (coupled) coordinates $(x,y,p_{x},p_{y})$ to give the decoupled coordinates $(q_{1},q_{2},p_{1},p_{2})$:

[TABLE]

since, $C$ satisfies $C^{T}=C^{-1}=-C$ and some examples are given in Appendix B.

Let us consider the symplectic transformation

[TABLE]

The coupled (non-separable) coordinates become

[TABLE]

This set of transformed coordinates applied to (20) gives the non-separable quadratic Hamiltonian

[TABLE]

which gives the Hamiltonian vector field

[TABLE]

where the equilibrium point is at $(0,0,0,0)$, and its total energy is [math]. The Jacobian at this equilibrium point has eigenvalues of type $\lambda,-\lambda,i\omega_{2},-i\omega_{2}$ (as shown in the Appendix A) and is of saddle $\times$ center type, that is index-1.

In the decoupled (separable) quadratic Hamiltonian, the dividing surface is given by $q_{1}=0$ Eqn. (23) which in the transformed coordinates becomes $p_{x}=0$. The dividing surface on the fixed energy surface $\mathcal{H}=h$ is given by

[TABLE]

In the decoupled coordinates, the NHIM is defined by $q_{1}=0,p_{1}=0$ on a fixed energy surface, which gives $-x+p_{x}+p_{y}=0$ that is $p_{y}=x$ in the coupled coordinates, and is given by

[TABLE]

Next, in the coupled coordinates, the stable and unstable manifolds of the NHIM are given by

[TABLE]

Invertible linear symplectic transformations are clearly $C^{\infty}$ diffeomorphisms. They not only preserve the Hamiltonian structure, but they preserve the phase space structures relevant to reaction dynamics. In particular, Lyapunov exponents are preserved under symplectic diffeomorphisms which implies that normal hyperbolicity is preserved. The no-recrossing property of the dividing surface is a result of transversality of the Hamiltonian vector field to the dividing surface, and such transversality properties are preserved under diffeomorphisms. Hence the NHIM (50) also has a geometry $\mathbb{S}^{1}$ and the invariant manifolds (51) and (52) are also of geometry $\mathbb{R}\times\mathbb{S}^{1}$.

3.2.1 Detecting the unstable periodic orbit and its manifolds

Now we illustrate the procedure for detecting the unstable periodic orbit, its stable and unstable manifolds by computing the Lagrangian descriptor on isoenergetic two-dimensional surfaces parametrized by a pair of coordinates of the separable quadratic Hamiltonian (47).

Isoenergetic two-dimensional surface parametrized by $(x,p_{x})$ — On a fixed energy surface, $\mathcal{H}(x,p_{x},y,p_{y})=h$, we consider a two-dimensional surface parametrized by $(x,p_{x})$ coordinates by defining

[TABLE]

Combining the dividing surface condition $p_{x}=0$, NHIM condition $p_{y}=x$, the directionality condition $\dot{y}(x,y,p_{x},p_{y};h)>0$, gives $x>0$. Thus, for a fixed energy, the intersection of the one-dimensional NHIM (50) with the isoenergetic two-dimensional surface (53) is given by

[TABLE]

which is a point on the line $p_{x}=0$ and $x=\sqrt{2h/\omega_{2}}$, and is marked by the red cross in Fig. 3. This point is also identified by the minimum value of the Lagrangian descriptor at this coordinate (as shown along the one-dimensional slice at constant $p_{x}=0$) and agrees with (54) to within the grid resolution used for the discretization of the two-dimensional surface. Since the NHIM is the intersection of stable and unstable manifolds, its coordinate is also a minimum and singular point in the LD values.

Next, the intersection of the unstable and stable manifolds (52) with the isoenergetic two-dimensional surface $U_{xp_{x}}^{+}$ is given by

[TABLE]

where the inequalities are derived from the directionality condition $\dot{y}>0$ for the surface. This is detected by the minima in the LD plot in Fig. 3 and highlighted by the dashed red (unstable) and dashed white (stable) curves. Points on these manifolds are also picked up by the one-dimensional slices along the constant $p_{x}=-0.60,0.30$.

Isoenergetic two-dimensional surface parametrized by $(x,y)$ — On a fixed energy surface $\mathcal{H}(x,p_{x},y,p_{y})=h$, we compute the Lagrangian descriptor on a two-dimensional surface parametrized by $(x,y)$ coordinates by defining

[TABLE]

Thus, the intersection of the NHIM (50) with this isoenergetic two-dimensional surface (57) is given by

[TABLE]

which represents a one-dimensional curve and is marked as a dashdot red line in Fig. 3. The NHIM is also identified by the minima in Lagrangian descriptor values along one-dimensional sections at constant $y$ in Fig. 3.

Next, the intersection of the unstable and stable manifolds (52) with the isoenergetic two-dimensional surface $U_{xy}^{+}$ is given by

[TABLE]

where each manifold represents the same one-dimensional curve as the NHIM and are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 3. These manifolds are also identified by points of minima in the Lagrangian descriptor values as shown by one-dimensional sections at constant $y$ values.

Isoenergetic two-dimensional surface parametrized by $(y,p_{y})$ — On the fixed energy surface $\mathcal{H}(x,p_{x},y,p_{y})=h$, we take the Lagrangian descriptor slice by defining a two-dimensional surface

[TABLE]

Thus, the intersection of the NHIM (50) with this isoenergetic two-dimensional surface (61) is given by

[TABLE]

which is a point $y=\sqrt{2h/\omega_{2}}$ on the $p_{y}=0$ line and shown as a red cross in Fig.3. This point is also identified by the minima in the Lagrangian descriptor values as shown along one-dimensional slice in Fig. 3.

Next, the intersection of the unstable (51) and stable manifolds (52) with the isoenergetic surface $U_{yp_{y}}^{+}$ manifest as

[TABLE]

which are shown as dashed red (unstable) and dashed white (stable) curves in the Fig. 3. These curves are also identified by the points of minima in the Lagrangian descriptor values as evident by the one-dimensional slices in Fig. 3.

4 Three Degrees of Freedom

In this section we use the Lagrangian Descriptor method to identify the NHIM and its stable and unstable manifolds for a 3 DoF separable quadratic Hamiltonian with an index-1 saddle point at the origin. The advantage of using this model Hamiltonian is that we can compare the analytical expressions for the phase space structures present with features in LD plots.

4.1 Decoupled quadratic Hamiltonian

We consider the three degrees of freedom benchmark example in Ref.[47] which is a linear system (quadratic Hamiltonian) with an equilibrium point of saddle-center-center equilibrium type at the origin. The Hamiltonian of this system is of the form

[TABLE]

with the corresponding Hamiltonian vector field given by:

[TABLE]

Since the total Hamiltonian decouples into the independent Hamiltonians $H_{r}$, $H_{b_{1}}$ and $H_{b_{2}}$ we can analyze the phase portraits for each separately. In the language of chemical reaction dynamics $H_{r}$ corresponds to the “reactive mode”, $H_{b_{1}}$ and $H_{b_{2}}$ are “bath modes”. The equilibrium point $(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=(0,0,0,0,0,0)$ is an index-1 saddle point for the full three DoF system on the zero (total) energy surface.

In this system a trajectory “reacts” when its $q_{1}$ coordinate changes sign. Therefore the surface $q_{1}=0$ is a dividing surface (DS) for trajectories, separating “reactants” from “products” (i.e. before and after reaction). This DS is a five dimensional surface in the six dimensional phase space. We discuss its geometrical structure, in both the phase space and in a fixed energy surface, in more detail.

First, we consider the “energetics” of the reaction. In order for $q_{1}$ to change sign we must have $H_{r}>0$. Also, it is clear from the form of $H_{b_{1}}$ and $H_{b_{2}}$ that $H_{b_{1}}\geqslant 0$ and $H_{b_{2}}\geqslant 0$. Therefore, reaction requires $H=H_{r}+H_{b_{1}}+H_{b_{2}}>0$.

The energy surface in which reaction occurs is given by

[TABLE]

The intersection of the DS, $q_{1}=0$, with this energy surface is given by:

[TABLE]

This is the isoenergetic DS. It is a 4-sphere in the six dimensional $(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})$ phase space. It has two hemispheres that are DSs for the forward and backward reactions, respectively:

[TABLE]

[TABLE]

Thes two hemispheres “meet” at $p_{1}=0$:

[TABLE]

which is a normally hyperbolic invariant 3 sphere. invariance follows from (66) since for $q_{1}=p_{1}=0$ then $\dot{q}_{1}=\dot{p}_{1}=0$. Thus $q_{1}$ and $p_{1}$ always remain zero, and trajectories with these initial conditions remain on (71), i.e., $q_{1}=p_{1}=0$ is invariant. It is normally hyperbolic since the directions normal to (71), i.e. $q_{1}-p_{1}$, are linearized saddle like dynamics.

In order to understand the relationship between the NHIM and the index-1 saddle point we note that for $H_{r}=0,\,H_{b_{1}}=0$ and $H_{b_{2}}=0$ the NHIM reduces to the point $(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=(0,0,0,0,0,0)$, which is the index-1 saddle point on the energy surface $H_{r}+H_{b_{1}}+H_{b_{2}}=0$. Therefore as the total energy is increased from [math], with $H_{b_{1}}$ increasing from zero and/or $H_{b_{2}}$ increasing from zero, we see that the NHIM “grows” from the index-1 saddle point on the zero energy surface into an invariant 3 sphere. This shows how the “influence” of the index-1 saddle point is carried to higher energy sufaces on which the saddle point does not exist.

The stable and unstable manifolds of the NHIM are given the following equations:

[TABLE]

and they are four dimensional on a fixed five dimensional energy surface. The topology of the manifolds is the product of a line ($q_{1}=p_{1}$ or $q_{1}=-p_{1}$) with a 3 sphere ($\frac{\omega_{2}}{2}\left(p_{2}^{2}+q_{2}^{2}\right)+\frac{\omega_{3}}{2}\left(p_{3}^{2}+q_{3}^{2}\right)=H_{b_{1}}+H_{b_{2}}>0$), that is $\mathbb{R}\times\mathbb{S}^{3}$ and sometimes such geometry is referred to as “spherical cylinders”. Since the lines $q_{1}=p_{1}$ and $q_{1}=-p_{1}$ correspond to the contour $H_{r}=0$, in the six dimensional phase space, these manifolds have energy $H=H_{r}+H_{b_{1}}+H_{b_{2}}=0+H_{b_{1}}+H_{b_{2}}>0$.

To understand the dynamics governed by the invariant manifolds of the NHIM, let us choose an initial condition $(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})$ on $\mathcal{W}^{u}\left(\mathcal{M}(h)\right)$. Then as $t\rightarrow+\infty$ the $(q_{2},p_{2},q_{3},p_{3})$ components of the trajectory with this initial condition evolve quasiperiodically and the $(q_{1},p_{1})$ components grow at an exponential rate. Similarly, if we choose an initial condition on $\mathcal{W}^{s}\left(\mathcal{M}(h)\right)$, then as $t\rightarrow+\infty$ the $(q_{2},p_{2},q_{3},p_{3})$ components of the trajectory evolve quasiperiodically and the $(q_{1},p_{1})$ components decay to zero at an exponential rate as $t\rightarrow\infty$. In other words, trajectories starting on $\mathcal{W}^{u}\left(\mathcal{M}(h)\right)$ decay at an exponential rate to the NHIM as $t\rightarrow-\infty$ and trajectories starting on $\mathcal{W}^{s}\left(\mathcal{M}(h)\right)$ decay at an exponential rate to the NHIM as $t\rightarrow+\infty$.

4.1.1 Detecting NHIM and its manifolds

We consider isoenergetic two-dimensional surfaces parametrized by two coordinates and compute the Lagrangian descriptor in a square domain of size 2 units around the origin. We discretize the coordinates of the two dimensional surface and pick constant values for three of the four remaining coordinates, and use the total energy equation to solve for the sixth coordinate. Due to the form of the Hamiltonian (65), obtaining the coordinate from the constant energy condition is simply solving a quadratic equation.

Isoenergetic two-dimensional surface parametrized by $(q_{1},p_{1})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{1},p_{1})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{1}p_{1}}^{+}$ with the NHIM (71) becomes

[TABLE]

Thus, the NHIM is located at the origin $(0,0)$ and marked by a red cross in the LD plot (Fig. 4).

Next, the intersection of the two-dimensional surface with the unstable (72) and stable manifolds (73) is given by

[TABLE]

which are one-dimensional for a fixed energy, and represent lines passing through the origin shown as dashed red (unstable) and white (stable) lines, respectively, in Fig. 4. The only points of local minima in the LD plot (Fig. 4) also lie along the lines passing through the origin and correspond to the manifolds of the NHIM.

Isoenergetic two-dimensional surface parametrized by $(q_{2},p_{2})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{2},p_{2})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{2}p_{2}}^{+}$ with the NHIM (71) becomes

[TABLE]

Thus, the NHIM is the circle of radius $\sqrt{2h/\omega_{2}}$ and marked by a dashed line in the LD plot (Fig. 4).

Next, the intersection of the two-dimensional surface with the unstable (72) and stable manifolds (73) is given by

[TABLE]

which are one-dimensional for a fixed energy, and marked by dashed red (unstable) and white (stable) lines, respectively, in Fig. 4. The only points of local minima in the LD plot (Fig. 4) also lie along the lines passing through the origin and correspond to the manifolds of the NHIM.

Isoenergetic two-dimensional surface parametrized by $(q_{3},p_{3})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{3},p_{3})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{3}p_{3}}^{+}$ with the NHIM (71) becomes

[TABLE]

Thus, the NHIM is the circle of radius $\sqrt{2h/\omega_{3}}$ and marked by a dashed line in the LD plot (Fig. 4).

Next, the intersection of the two-dimensional surface with the unstable (72) and stable manifolds (73) is given by

[TABLE]

which are one-dimensional for a fixed energy and represent circles of radius $\sqrt{2h/\omega_{3}}$, and marked by dashed red (unstable) and white (stable) lines, respectively, in Fig. 4. The only points of minima and singularity in the LD plot (Fig. 4) is along the circle and thus identify the manifolds of the NHIM.

Isoenergetic two-dimensional surface parametrized by $(q_{1},q_{2})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{1},q_{2})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{1}q_{2}}^{+}$ with the NHIM (71) is given by

[TABLE]

Thus, on the isoenergetic two-dimensional surface $U_{q_{1}q_{2}}^{+}$, the NHIM is the line $q_{1}=0$ and marked by a dashdot line in the LD plot (Fig. 5).

Next, the intersection of the unstable (72) and stable manifolds (73) with two-dimensional surface is given by

[TABLE]

which are one-dimensional for a fixed energy, and marked by dashed red (unstable) and white (stable) lines, respectively, in Fig. 5. The only points of minima in the LD plot (Fig. 5) also lie along this line at $q_{1}=0$ and identify the manifolds of the NHIM.

Isoenergetic two-dimensional surface parametrized by $(q_{2},q_{3})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{2},q_{3})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{2}q_{3}}^{+}$ with the NHIM (71) is given by

[TABLE]

Thus, on the isoenergetic two-dimensional surface $U_{q_{2}q_{3}}^{+}$, the NHIM is an ellipse of semi-major axis $\sqrt{2h/\omega_{2}}$ and semi-minor axis $\sqrt{2h/\omega_{3}}$ if $\omega_{2}<\omega_{3}$, and vice-versa, otherwise. In our case, $\omega_{2}=\omega_{3}$, the intersection of the NHIM with the isoenergetic two-dimensional surface becomes a circle and marked by a dashdot line in the LD plot (Fig. 5).

Next, the intersection of the unstable (72) and stable manifolds (73) with the isoenergetic two-dimensional surface is given by

[TABLE]

which have the same geometry as the NHIM on the isoenergetic two-dimensional surface, $U_{q_{2}q_{3}}$, and marked by dashed red (unstable) and white (stable) lines, respectively, in Fig. 5. The only points of minima in the LD plot (Fig. 5) also lie along this circle and identify the manifolds of the NHIM.

Isoenergetic two-dimensional surface parametrized by $(q_{3},q_{1})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{3},q_{1})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{3}q_{1}}^{+}$ with the NHIM (71) is given by

[TABLE]

which represents two points $q_{3}=\pm\sqrt{2h/\omega_{3}}$ on the line $q_{1}=0$, and marked by a red cross in Fig. 5. These points are also identified by the minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $q_{1}$.

Next, the intersection of the unstable (72) and the stable manifold (73) with the isoenergetic two-dimensional surface is given by

[TABLE]

where each manifold represent lines parallel to $q_{1}$ axis. The unstable manifold lies on the $q_{1}>0$ plane and the stable manifold lies on the $q_{1}<0$ plane, and are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 5. These manifolds are again identified by points of minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $q_{1}$.

4.2 Coupled quadratic Hamiltonian

To couple the coordinates in the separable quadratic Hamiltonian (65), we introduce a linear transformation such that it satisfies the symplectic condition (43).

Let us consider the symplectic transformation

[TABLE]

The change of coordinates is given by

[TABLE]

This transformation applied to (65) gives the Hamiltonian in the transformed coordinates

[TABLE]

which gives the vector field

[TABLE]

where the equilibrium point is at $(0,0,0,0,0,0)$ and its total energy is [math]. The Jacobian at this equilibrium point has eigenvalues $\lambda,-\lambda,i\omega_{2},-i\omega_{2},i\omega_{3},-i\omega_{3}$ (as shown in the Appendix A) and is saddle $\times$ center $\times$ center, thus index-1. The energy surface is a five dimensional surface given by (106) in the six dimensional phase space.

In the decoupled (separable) system (65), the dividing surface is defined by $q_{1}=0$ which becomes $p_{x}=0$ in the coupled (non-separable) coordinates. The dividing surface for a fixed energy $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$ is given by

[TABLE]

On this dividing surface, the NHIM is defined by $p_{1}=0$ in the separable coordinates, which gives $-x+p_{x}+p_{y}+p_{z}=0$ that is $p_{y}+p_{z}=x$ in the non-separable coordinates, and can be expressed as

[TABLE]

Next, in the coupled coordinates, the stable and unstable manifolds of the NHIM are given by

[TABLE]

As noted earlier, invertible linear symplectic transformations preserve the normal hyperbolicity and the transversality of the Hamiltonian vector field to the dividing surface. Hence the NHIM (109) still has a geometry $\mathbb{S}^{3}$ and the invariant manifolds (110) and (111) are also $\mathbb{R}\times\mathbb{S}^{3}$.

4.2.1 Detecting NHIM and its manifolds

Now we illustrate the procedure for detecting NHIM and its stable and unstable manifolds using features in Lagrangian descriptor on isoenergetic two-dimensional surfaces for the non-separable quadratic Hamiltonian (47).

Isoenergetic two-dimensional surface parametrized by $(x,p_{x})$ — On the constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$ (106), we compute the Lagrangian descriptor on a two-dimensional surface parametrized by $(x,p_{x})$ coordinates by defining

[TABLE]

Thus, on the five-dimensional energy surface, the intersection of the three-dimensional NHIM (109) with the two-dimensional surface is zero-dimensional and given by

[TABLE]

which is a point on the line $p_{x}=0$, and marked by a red cross in Fig. 6. This is also identified by the location of minima in the Lagrangian descriptor values as shown by the one-dimensional slice in Fig. 6 and agrees with $x=\sqrt{2h/\left(\omega_{2}+2\omega_{3}\right)}$ to within the grid resolution.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (112) becomes

[TABLE]

which are one-dimensional and are shown as a dashed red line (unstable) and a dashed white line (stable) in Fig. 6. These lines are also identified by the minima in the Lagrangian descriptor values as shown by one-dimensional slice at constant $p_{x}$.

Isoenergetic two-dimensional surface parametrized by $(y,p_{y})$ — On a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$, we compute Lagrangian descriptor on a two-dimensional surface by defining

[TABLE]

The intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which is a point on the isoenergetic two-dimensional surface (116), and shown as red cross in Fig. 6. This point is also identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices for constant $y$.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (116) becomes

[TABLE]

which are lines on the two-dimensional isoenergetic surface (116) and are also identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices at constant $y$.

Isoenergetic two-dimensional surface parametrized by $(z,p_{z})$ — On a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$ (106), we compute Lagrangian descriptor on a two-dimensional surface by defining

[TABLE]

The intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which is a point and shown as red cross in Fig. 6. This point is also identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices for constant $z$.

Next, the intersection of the four-dimensional unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (120) becomes

[TABLE]

which are one-dimensional and shown as dashed red (unstable) and dashed white (stable) lines in Fig. 6. These invariant manifolds are also identified by minima in the Lagrangian descriptor values as shown by the one-dimensional slice at constant $z$.

Isoenergetic two-dimensional surface parametrized by $(x,p_{z})$ — On a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$, we compute the Lagrangian descriptor on a two-dimensional surface by defining

[TABLE]

The intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which is also identified by the minima in the Lagrangian descriptor values as shown in the one-dimensional slices at constant $p_{z}$.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (124) becomes

[TABLE]

which project as the same one-dimensional curve as the NHIM and identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices at constant $p_{z}$ in Fig. 7.

Isoenergetic two-dimensional surface parametrized by $(y,p_{z})$ — Next, on a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$, we compute Lagrangian descriptor on a two-dimensional surface by defining

[TABLE]

Thus, the intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which represents a portion of the ellipse and shown as the dashdot red curve in Fig. 7. As shown by the one-dimensional slices for constant $p_{z}$, Eqn. (129) is identified by the minima in the Lagrangian descriptor values.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (128) becomes

[TABLE]

which are shown as red dashed line (unstable) and white dashed line (stable) in Fig. 7. We note here that the intersection of the manifolds with the isoenergetic two-dimensional surface (128) overlays on the NHIM’s intersection (129). These curves are also identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slice for constant $p_{z}$.

Isoenergetic two-dimensional surface parametrized by $(x,z)$ — Next, on a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$ (106), we compute the Lagrangian descriptor on a two-dimensional surface by defining

[TABLE]

The intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which represents a portion of an ellipse and is identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices for constant $z$ in Fig. 7.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (128) becomes

[TABLE]

which are also be identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slice for constant $z$ in Fig. 7.

5 Summary and Outlook

In this article, we assessed Lagrangian descriptor (LD) based detection of high dimensional phase space structures, that is normally hyperbolic invariant manifolds, their unstable and stable manifolds, in the two and three degrees-of-freedom quadratic Hamiltonian systems with index-1 saddle. This is done using a sytematic comparison of numerical results with analytical expressions for a comprehensive set of coordinates; see Appendix C for more pairs of coordinates that support the results herein.. It is to be noted that in $N$ DoF system, there are $2N(2N-1)(N-1)$ two-dimensional sections, so the detection strategy should include inspecting a few pairs of coordinates. Based on our investigation, a judicious choice of coordinate for the two-dimensional surface is to use one of the reaction coordinates in phase space, that is $q_{1}$ or $p_{1}$, and the other one from the remaining $2N-1$ coordinates.

The closed-form analytical expressions of the NHIM and their manifolds are precisely identified by the minima and singular features in the Lagrangian descriptor values. This provides further numerical evidence for the use of LDs in detecting high dimensional phase space structures with a simple computation that can be implemented along with trajectory integration. So at least for the form of Hamiltonian considered here, one can rely with certainity on using Lagrangian descriptor to detect the NHIM and its stable and unstable manifolds. Furthermore, this detection approach has the potential to be combined with machine learning type methods [16]. Future work on this approach will include specific nonlinear systems that are inspired by applications in celestial mechanics, ship dynamics, structural mechanics, and chemical reaction dynamics.

Acknowledgements

We acknowledge the support of EPSRC Grant No.  EP/P021123/1 and ONR Grant No. N00014-01-1-0769. We would like to thank Vladimír Krajňák for useful discussions and feedback.

Appendix

A Eigenvalues of linear symplectic transformed system

The Jacobian of the vector field (48) near the equilibrium point is given by

[TABLE]

which gives the characteristic polynomial

[TABLE]

which has solutions $\lambda,-\lambda,i\omega_{2},-i\omega_{2}$.

The Jacobian of the vector field (107) near the equilibrium point is given by

[TABLE]

which gives the characteristic polynomial

[TABLE]

which has solutions $\lambda,-\lambda,i\omega_{2},-i\omega_{2},i\omega_{3},-i\omega_{3}$.

B Examples of symplectic transformations

Along with the transformations $C:(x,y,p_{x},p_{y})\rightarrow(q_{1},q_{2},p_{1},p_{2})$ used for demonstrating the procedure of detecting NHIM, its unstable, and stable manifolds, we have verified the following transformations for symplectic condition.

Two degrees of freedom

[TABLE]

Three degrees of freedom

The following $C:(x,y,z,p_{x},p_{y},p_{z})\rightarrow(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})$ have been verfied to be symplectic.

[TABLE]

C Lagrangian Descriptor on isoenergetic two-dimensional surfaces

C.1 Decoupled quadratic Hamiltonian: 2 DoF

Isoenergetic two-dimensional surface parametrized by $(p_{1},p_{2})$ — Similarly, on the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2})=h$, we take the Lagrangian descriptor slice in $(q_{2},p_{2})$ coordinates by defining the two dimensional surface

[TABLE]

where

[TABLE]

The intersection of the NHIM (25) with the two-dimensional surface $U_{p_{1}p_{2}}^{+}$ is

[TABLE]

which represents two points on the line $p_{1}=0$. Their location is $p_{2}=\pm\sqrt{2h/\omega_{2}}$ which is $\approx 0.632$ for $h=0.2,\omega_{2}=1.0$ and is marked as red crosses in Fig. 8. These points are also locations of minima in LD values as shown by the one-dimensional slices at constant $p_{1}$.

Next, intersection of the two-dimensional $U_{p_{1}p_{2}}^{+}$ with the stable (27) and unstable manifolds (26) is given by

[TABLE]

where each manifold represent lines parallel to $p_{1}$ axis for $p_{2}=\pm\sqrt{\frac{2h}{\omega_{2}}}$. The unstable manifold lies on the $p_{1}>0$ plane and stable manifold lies on the $p_{1}<0$ plane, and are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 8.

Isoenergetic two-dimensional surface parametrized by $(q_{1},p_{2})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2})=h$, we take the LD slice by defining the two dimensional surface

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{1}p_{2}}^{+}$ with the NHIM (25) is given by

[TABLE]

which gives two points on the $q_{1}=0$ line with $p_{2}$ satisfying $p_{2}=\pm\sqrt{2h/\omega_{2}}$. This is shown by the dashdot line in Fig. 8.

Next, the intersection of the two-dimensional $U_{q_{1}p_{2}}^{+}$ with the stable (27) and unstable manifolds (26) is given by

[TABLE]

where each manifold represents lines parallel to the $q_{1}$ axis. The unstable manifold lies on the $q_{1}>0$ plane and the stable manifold lies on the $q_{1}<0$ plane, and are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 8. These manifolds are identified by points of minima located on the $q_{1}=0$ line.

Isoenergetic two-dimensional surface parametrized by $(q_{2},p_{1})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2})=h$, we take the LD slice by defining the two dimensional surface

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{2}p_{1}}^{+}$ with the NHIM (25) is given by

[TABLE]

which gives points on the $p_{1}=0$ line with $q_{2}$ satisfying $q_{2}=\pm\sqrt{2h/\omega_{2}-p_{2}^{2}}$. This is shown as a dashdot red line in Fig. 8 and is also the line with minima as shown by the one-dimensional slice at constant $q_{2}$.

Next, the intersection of the two-dimensional surface $U_{q_{2}p_{1}}^{+}$ with the stable (27) and unstable manifolds (26) is given by

[TABLE]

where each manifold represents a line parallel to the $q_{1}$ axis. The unstable manifold lies on the $q_{1}>0$ plane and stable manifold lies on the $q_{1}<0$ plane, and are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 8.

C.2 Coupled quadratic Hamiltonian: 2 DoF

Isoenergetic two-dimensional surface parametrized by $(p_{x},p_{y})$ — On the fixed energy surface $\mathcal{H}(x,p_{x},y,p_{y})=h$, we compute the Lagrangian descriptor by defining a two-dimensional surface

[TABLE]

Thus, the intersection of the NHIM (50) with this isoenergetic two-dimensional surface (61) is given by

[TABLE]

which represents a point at the origin on the two-dimensional surface, $U_{p_{x}p_{y}}^{+}$, and shown as a red cross in Fig. 9. This point is also identified by the minima in the Lagrangian descriptor values as evident by the one-dimensional slices at constant $p_{y}$.

Next, the intersection of the unstable (51) and stable manifolds (52) with the isoenergetic surface $U_{yp_{y}}^{+}$ manifest as

[TABLE]

which are shown as dashed red (unstable) and dashed white (stable) curves in the Fig. 9. These curves are also identified by the points of minima in the Lagrangian descriptor values as evident by the one-dimensional sections in Fig. 9.

Isoenergetic two-dimensional surface parametrized by $(x,p_{y})$ — On the fixed energy surface $\mathcal{H}(x,p_{x},y,p_{y})=h$, we compute the Lagrangian descriptor by defining a two-dimensional surface

[TABLE]

Thus, the intersection of the NHIM (50) with this isoenergetic two-dimensional surface (161) is given by

[TABLE]

which represents a point at the origin on the two-dimensional surface, $U_{xp_{y}}^{+}$, and shown as a red cross in Fig. 9. This point is also identified by the minima in the Lagrangian descriptor values as evident by the one-dimensional slices at constant $p_{y}$.

Next, the intersection of the unstable (51) and stable manifolds (52) with the isoenergetic surface $U_{xp_{y}}^{+}$ manifest as

[TABLE]

which are shown as dashed red (unstable) and dashed white (stable) curves in the Fig. 9. These curves are also identified by the points of minima in the Lagrangian descriptor values as evident by the one-dimensional sections in Fig. 9.

Isoenergetic two-dimensional surface parametrized by $(y,p_{x})$ — On the fixed energy surface $\mathcal{H}=h$, we compute the Lagrangian descriptor by defining a two-dimensional surface

[TABLE]

Thus, the intersection of NHIM (50) with this isoenergetic two-dimensional surface (165) is given by

[TABLE]

which represents a point at the origin on the two-dimensional surface, $U_{yp_{x}}^{+}$, and shown as a red cross in Fig. 9. This point is also identified by the minima in the Lagrangian descriptor values as evident by the one-dimensional slices at constant $p_{x}$.

Next, the intersection of the unstable (51) and stable manifolds (52) with the isoenergetic surface $U_{xp_{y}}^{+}$ manifest as

[TABLE]

which are shown as dashed red (unstable) and dashed white (stable) curves in the Fig. 9. These curves are also identified by the points of minima and singularities in the Lagrangian descriptor values as evident by the one-dimensional sections in Fig. 9.

C.3 Decoupled quadratic Hamiltonian: 3 DoF

Isoenergetic two-dimensional surface parametrized by $(p_{2},p_{3})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(p_{2},p_{3})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{p_{2}p_{3}}^{+}$ with the NHIM (71) is given by

[TABLE]

which represents points on an ellipse, and marked by a dashed line in Fig. 10. These points are also identified by the minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $p_{3}$.

Next, the intersection of the unstable (72) and stable manifolds (73) with the isoenergetic two-dimensional surface $U_{p_{2}p_{3}}^{+}$ is given by

[TABLE]

where both the manifolds represent an ellipse and are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 10. These manifolds are again identified by points of minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $p_{3}$.

Isoenergetic two-dimensional surface parametrized by $(q_{1},p_{3})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{1},p_{3})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{1}p_{3}}^{+}$ with the NHIM (71) is given by

[TABLE]

which represents points on the line $q_{1}=0$, and marked by a dashed line in Fig. 10. These points are also identified by the minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $p_{3}$.

Next, the intersection of the unstable (72) and stable manifolds (73) with the isoenergetic two-dimensional surface $U_{q_{1}p_{3}}^{+}$ is given by

[TABLE]

where both manifolds represent points on the line $q_{1}=0$ which are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 10. These manifolds are again identified by points of minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $p_{3}$.

Isoenergetic two-dimensional surface parametrized by $(q_{3},p_{2})$ — On the constant energy surface, $H(q_{1},p_{1},q_{2},p_{2},q_{3},p_{3})=h$, we compute Lagrangian descriptor on a two-dimensional surface parametrized by $(q_{3},p_{2})$ coordinates by defining

[TABLE]

where

[TABLE]

The intersection of the two-dimensional surface $U_{q_{3}p_{2}}^{+}$ with the NHIM (71) is given by

[TABLE]

which represents points on an ellipse, and marked by dashed line in Fig. 10. These points are also identified by the minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $p_{2}$.

Next, the intersection of the unstable (72) and stable manifolds (73) with the isoenergetic two-dimensional surface $U_{q_{3}p_{2}}^{+}$ is given by

[TABLE]

where both manifolds are ellipses and are marked by dashed red (unstable) and dashed white (stable) lines in the Fig. 10. These manifolds are again identified by points of minima in the Lagrangian descriptor values as shown by one-dimensional slices at constant $p_{2}$.

C.4 Coupled quadratic Hamiltonian: 3 DoF

Isoenergetic two-dimensional surface parametrized by $(y,z)$ — Next, on a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$, we compute the Lagrangian descriptor on a two-dimensional surface by defining

[TABLE]

The intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which is an ellipse as shown in Fig. 11 and is also identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices for constant $y$.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (128) becomes

[TABLE]

which are also be identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slice for constant $y$ in Fig. 11.

Isoenergetic two-dimensional surface parametrized by $(z,p_{y})$ — Next, on a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$ (106), we consider a two-dimensional surface by defining

[TABLE]

Thus, the intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which represents a portion of an ellipse and is identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices for constant $p_{y}$ in Fig. 11.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (188) becomes

[TABLE]

which are identified by the minima in the Lagrangian descriptor values in Fig. 11 as evident by the one-dimensional slice at constant $p_{y}$.

Isoenergetic two-dimensional surface parametrized by $(p_{y},p_{z})$ — Next, on a constant energy surface, $\mathcal{H}(x,p_{x},y,p_{y},z,p_{z})=h$, we compute the Lagrangin descriptor on a two-dimensional surface by defining

[TABLE]

The intersection of the NHIM (109) with the two-dimensional surface is given by

[TABLE]

which are two points and shown as the red crosses in Fig. 11. These are also identified by the minima in the Lagrangian descriptor values as shown by the one-dimensional slices for constant $p_{z}$.

Next, the intersection of the unstable (110) and stable (111) manifolds with the isoenergetic two-dimensional surface (192) becomes

[TABLE]

which are also identified by the minima in the Lagrangian descriptor values and shown by the one-dimensional slice for constant $p_{z}$ in Fig. 11.