Tilting and Squeezing: Phase space geometry of Hamiltonian saddle-node bifurcation and its influence on chemical reaction dynamics

TL;DR

This paper investigates how a Hamiltonian saddle-node bifurcation alters phase space structures that govern chemical reaction dynamics, using Lagrangian descriptors to analyze invariant manifolds in high-dimensional systems.

Contribution

It introduces a phase space analysis framework for Hamiltonian saddle-node bifurcations, extending from one to two degrees of freedom, and examines the impact on reaction dynamics.

Findings

Identification of invariant manifolds via Lagrangian descriptors.

Description of phase space changes due to bifurcation.

Analysis of NHIMs and their manifolds in reaction dynamics.

Abstract

In this article we present the influence of a Hamiltonian saddle-node bifurcation on the high-dimensional phase space structures that mediate reaction dynamics. To achieve this goal, we identify the phase space invariant manifolds using Lagrangian descriptors, which is a trajectory-based diagnostic suitable for the construction of a complete ``phase space tomography'' by means of analyzing dynamics on low-dimensional slices. First, we build a Hamiltonian system with one degree-of-freedom (DoF) that models reaction, and study the effect of adding a parameter to the potential energy function that controls the depth of the well. Then, we extend this framework to a saddle-node bifurcation for a two DoF Hamiltonian, constructed by coupling a harmonic oscillator, i.e. a bath mode, to the other reactive DoF in the system. For this problem, we describe the phase space structures associated with the rank-1 saddle equilibrium point in the bottleneck region, which is a Normally Hyperbolic Invariant Manifold (NHIM) and its stable and unstable manifolds. Finally, we address the qualitative changes in the reaction dynamics of the Hamiltonian system due to changes in the well depth of the potential energy surface that gives rise to the saddle-node bifurcation.

Equations

$$H(u,v)=\frac{1}{2}{v}^2-\mu u+\frac{1}{3}{u}^3,$$

$$\{\begin{array}{rl}\dot{u}& =\frac{\partial H}{\partial v}=v\\ \dot{v}& =-\frac{\partial H}{\partial u}=\mu -u^2\end{array},$$

$$\mathbb{J}(u,v)=\left(\begin{array}{cc}\frac{{\partial }^{2}H}{\partial u\partial v}& \frac{{\partial }^{2}H}{\partial v^2}\\ -\frac{{\partial }^{2}H}{\partial u^2}& -\frac{{\partial }^{2}H}{\partial v\partial u}\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -2u& 0\end{array}\right).$$

$$\mathbb{J}(-\sqrt{\mu },0)=\left(\begin{array}{cc}0& 1\\ 2\sqrt{\mu }& 0\end{array}\right),$$

$$\mathbb{J}(\sqrt{\mu },0)=\left(\begin{array}{cc}0& 1\\ -2\sqrt{\mu }& 0\end{array}\right);,$$

$$\{\begin{array}{rl}u& =q-\sqrt{\mu }\\ v& =p\end{array},\mu \in {\mathbb{R}}^{+}\cup \{0\},$$

$$\{\begin{array}{rl}\dot{q}& =p\\ \dot{p}& =2\sqrt{\mu }q-q^2\end{array},$$

$$H(q,p)=\frac{1}{2}{p}^2-\sqrt{\mu }{q}^2+\frac{1}{3}{q}^3.$$

$$\mathbb{J}(q,p)=\left(\begin{array}{cc}\frac{{\partial }^{2}H}{\partial q\partial p}& \frac{{\partial }^{2}H}{\partial p^2}\\ -\frac{{\partial }^{2}H}{\partial q^2}& -\frac{{\partial }^{2}H}{\partial p\partial q}\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 2\sqrt{\mu }-2q& 0\end{array}\right)$$

$$\mathbb{J}(0,0)=\left(\begin{array}{cc}0& 1\\ 2\sqrt{\mu }& 0\end{array}\right),\mathbb{J}(2\sqrt{\mu },0)=\left(\begin{array}{cc}0& 1\\ -2\sqrt{\mu }& 0\end{array}\right)$$

$$H(q,p)=\frac{1}{2}{p}^2-\sqrt{\mu }{q}^2+\frac{\alpha }{3}{q}^3\equiv T(p)+V(q),$$

$$\{\begin{array}{l}\dot{q}=\frac{\partial H}{\partial p}=p\\ \dot{p}=-\frac{\partial H}{\partial q}=2\sqrt{\mu }q-\alpha q^2\end{array}.$$

$$J(q,p)=\left(\begin{array}{cc}0& 1\\ 2\sqrt{\mu }-2\alpha q& 0\end{array}\right).$$

$$J(0,0)=\left(\begin{array}{cc}0& 1\\ 2\sqrt{\mu }& 0\end{array}\right),J\left(2\sqrt{\mu }/\alpha ,0\right)=\left(\begin{array}{cc}0& 1\\ -2\sqrt{\mu }& 0\end{array}\right).$$

$$V(q)=-\sqrt{\mu }{q}^2+\frac{\alpha }{3}{q}^3,$$

$$\mathcal{F}\equiv V(0)-V\left(\frac{2\sqrt{\mu }}{\alpha }\right)=-\frac{4\sqrt{{\mu }^3}}{3{\alpha }^2}.$$

$$\mathcal{F}=\frac{\sqrt{\mu }}{3}{\mathcal{D}}^2,$$

$$\frac{d\mathcal{F}}{d\mathcal{D}}=\frac{2\sqrt{\mu }}{3}\mathcal{D}=\frac{1}{3}{\lambda }_0^2\mathcal{D}.$$

$$H(q,x,p,p_x)=\frac{1}{2}\left(p^2+p_x^2\right)-\sqrt{\mu }{q}^2+\frac{\alpha }{3}{q}^3+\frac{{\omega }^2}{2}{x}^2+\frac{\epsilon }{2}{\left(x-q\right)}^2,$$

$$T(p,p_x)=\frac{1}{2}\left(p^2+p_x^2\right),V(q,x)=-\sqrt{\mu }{q}^2+\frac{\alpha }{3}{q}^3+\frac{{\omega }^2}{2}{x}^2+\frac{\epsilon }{2}{\left(x-q\right)}^2$$

$$\{\begin{array}{rl}\dot{q}& =\frac{\partial H}{\partial p}=p\\ \dot{x}& =\frac{\partial H}{\partial p_x}=p_x\\ \dot{p}& =-\frac{\partial H}{\partial q}=2\sqrt{\mu }q-\alpha q^2+\epsilon (x-q)\\ {\dot{p}}_x& =-\frac{\partial H}{\partial x}=-{\omega }^{2}x+\epsilon (q-x)\end{array}$$

$$q_e=\frac{2\sqrt{\mu }}{\alpha }-\frac{{\omega }^2\epsilon }{\alpha ({\omega }^2+\epsilon )},x_e=\frac{\epsilon }{{\omega }^2+\epsilon }{q}_e$$

$$H({\mathbf{x}}_1^e)=0,H({\mathbf{x}}_2^e)=\left(-2\sqrt{\mu }+\frac{{\omega }^4\epsilon -2{\omega }^2{\epsilon }^2}{{\left({\omega }^2+\epsilon \right)}^2}\right)\frac{q_e^2}{6}$$

$${\epsilon }_c=\frac{2\sqrt{\mu }{\omega }^2}{{\omega }^2-2\sqrt{\mu }},$$

$$\mathbb{J}(q,x,p,p_x)=\left(\begin{array}{cccc}\frac{{\partial }^{2}H}{\partial q\partial p}& \frac{{\partial }^{2}H}{\partial x\partial p}& \frac{{\partial }^{2}H}{\partial p^2}& \frac{{\partial }^{2}H}{\partial p_x\partial p}\\ \frac{{\partial }^{2}H}{\partial q\partial p_x}& \frac{{\partial }^{2}H}{\partial x\partial p_x}& \frac{{\partial }^{2}H}{\partial p\partial p_x}& \frac{{\partial }^{2}H}{\partial p_x^2}\\ -\frac{{\partial }^{2}H}{\partial q^2}& -\frac{{\partial }^{2}H}{\partial x\partial q}& -\frac{{\partial }^{2}H}{\partial p\partial q}& -\frac{{\partial }^{2}H}{\partial p_x\partial q}\\ -\frac{{\partial }^{2}H}{\partial q\partial x}& -\frac{{\partial }^{2}H}{\partial x^2}& -\frac{{\partial }^{2}H}{\partial p\partial x}& -\frac{{\partial }^{2}H}{\partial p_x\partial x}\end{array}\right)=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 2\sqrt{\mu }-2\alpha q-\epsilon & \epsilon & 0& 0\\ \epsilon & -{\omega }^2-\epsilon & 0& 0\end{array}\right)$$

$$\mathbb{J}({\mathbf{x}}_1^e)=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 2\sqrt{\mu }-\epsilon & \epsilon & 0& 0\\ \epsilon & -{\omega }^2-\epsilon & 0& 0\end{array}\right)$$

$$\det \left(\mathbb{J}({\mathbf{x}}_1^e)-\beta \mathbb{I}\right)=\left({\beta }^2-2\sqrt{\mu }+\epsilon \right)\left({\beta }^2+{\omega }^2+\epsilon \right)-{\epsilon }^2=0$$

$${\lambda }_0=\sqrt[4]{4\mu },{\omega }_0=\omega$$

$$\left({\beta }^2-{\lambda }_0^2+\epsilon \right)\left({\beta }^2+{\omega }_0^2+\epsilon \right)-{\epsilon }^2={\beta }^4+\left({\omega }_0^2-{\lambda }_0^2+2\epsilon \right){\beta }^2+\left({\omega }_0^2-{\lambda }_0^2\right)\epsilon -{\lambda }_0^2{\omega }_0^2=0$$

$$\xi =\frac{{\lambda }_0^2-{\omega }_0^2}{2}-\epsilon \pm \sqrt{{\left(\frac{{\lambda }_0^2+{\omega }_0^2}{2}\right)}^2+{\epsilon }^2}.$$

Full text

Keywords: Saddle-Node bifurcation, Hamiltonian systems, Phase space structure, Lagrangian descriptors, Chemical reaction dynamics.

1. I Introduction

2. II Hamiltonian Saddle-Node Bifurcation Models

3. II.1 One Degree-of-Freedom Hamiltonian

4. II.2 Two degrees-of-freedom Hamiltonian

5. III Results and Discussion: Influence of well-depth on reaction dynamics

6. III.1 One Degree-of-Freedom Hamiltonian

7. III.2 Two degree-of-freedom Hamiltonian

8. IV Conclusions

9. A Lagrangian Descriptors

10. B Visualization of phase space structure for the quadratic normal form Hamiltonian

11. C Computation of NHIM and its invariant manifolds for the 2 DoF system

I Introduction

Background.

Molecular reaction dynamics is concerned with the breaking and forming of bonds between the atoms that make up a molecule. The description of such mechanisms requires analyzing how the bonds change in time, that is, it is a dynamical mechanism requiring a phase space description. The phase space will have coordinates that correspond to the bond configurations with a corresponding canonical momentum for each bond configuration coordinate. Thus, the dynamics of chemical bonds can be formulated in the language of Hamiltonian dynamics and its resulting phase space geometry. In recent years there has been significant developments in the mathematical description and analysis of chemical reaction dynamics in phase space Komatsuzaki and Berry (2000); Uzer et al. (2002); Komatsuzaki and Berry (2002).

Our objective is to describe the change in a bond that signals the occurrence of a chemical reaction in terms of a unique characteristic of a trajectory of Hamilton’s equations. This characteristic requires an understanding of the geometry of the phase space of the molecule in a way that enables us to divide an appropriate volume of phase space into a region corresponding to “reactants” and a region corresponding to “products”. The passage from reactants to products, that is “reaction”, occurs when a trajectory crosses the “dividing surface” (DS) between the reactants and products. In this framework for understanding chemical reactions, the flux through such a dividing surface would be related to the reaction rate, and therefore the construction of DS between such regions is of interest in reaction dynamics. The description of the regions of reactants and products can often be inferred from the nature of the development of the coordinates used in the mathematical model of the chemical reaction. Once the model is developed, the DS must be constructed in the context of this model. In relating the flux through the DS to a reaction rate it is desirable that trajectories crossing the DS from the reactant side proceed to the product side before possible recrossing the DS back to the reactant region. Thus, we require the DS to have the “no-recrossing” propertyMacKay (1990); Waalkens and Wiggins (2004). Hamiltonian dynamics conserves energy and therefore the trajectories evolve on a fixed energy surface of one less dimension (referred to as codimension one) than the phase space. The DS is required to be codimension one in the energy surface. Hence, we require a DS between reactants and products to be codimension one in the energy surface and to have the no-recrossing property. We note that Wigner had already described these properties for a DS in phase space much earlierWigner (1938, 1939), and a review of classical and quantum versions of DS constructed in phase space can be found in Waalkens et al. (2008).

Traditionally, the construction of DS was initially focused on critical points of the potential energy surface (PES), that is, in the configuration space describing the molecular system. Critical points on the PES do have significance in phase space; they are the equilibrium points for zero momentum. But they continue to have influence for nonzero momentum for a range of energies above the energy of the equilibrium point. The precise manner of this dynamical influence has only been understood recently and we will describe this shortly Komatsuzaki and Nagoaka (1997); Komatsuzaki and Berry (2000); Waalkens and Wiggins (2004). The construction of a DS separating the phase space into two parts, reactants and products, has been a focus from the dynamical systems point of view in recent years. However, the lack of a firm theoretical basis for the construction of such surfaces for molecular systems with three and more degree-of-freedom (DoF) has until recently been a major obstacle in the development of the theory. In phase space, that is for nonzero momentum, the role of the saddle point is played by an invariant manifold of saddle stability type, the normally hyperbolic invariant manifold (NHIM) (see §: II.2) Wiggins (1988, 1990, 2013). In order to fully appreciate the NHIM and its role in reaction rate theory, it is useful to begin with a precursor concept – the periodic orbit dividing surface or PODS. For systems with two DoF described by a natural Hamiltonian, kinetic plus potential energy, the problem of constructing the DS in phase space was solved during the 1970s by McLafferty, Pechukas and Pollak Pechukas and McLafferty (1973); Pechukas and Pollak (1977); Pollak and Pechukas (1978); Pechukas and Pollak (1979). They demonstrated that the DS at a specific energy is related to an invariant phase space structure, an unstable periodic orbit (UPO). The UPO defines (it is the boundary of) the bottleneck in phase space through which the reaction occurs and the DS which intersects trajectories evolving from reactants to products can be shown to have the geometry of a hemisphere in phase space whose boundary is the unstable PO Wiggins et al. (2001a); Waalkens and Wiggins (2004). The same construction can be carried out for a DS intersecting trajectories crossing from products to reactants and these two hemispheres form a sphere for which the UPO is the equator. Generalisation of this construction of DS to high dimensional systems has been a central question in reaction dynamics and has only received a satisfactory answer in recent years Wiggins et al. (2001a); Uzer et al. (2002). The key difficulty concerns the high dimensional analogue of the unstable PO used in the two DoF system for the construction of the DS. This difficulty is resolved by considering the NHIM, which has the appropriate dimensionality for anchoring the dividing surface in phase space.

Results from dynamical systems theory show that transport in phase space is controlled by high dimensional manifolds, NHIMs, which are the natural generalisation of the UPO of the two DoF case Wiggins (1990). Normal hyperbolicity of these invariant manifolds means that their stability, in a precise sense, is of saddle type in the transverse direction, which implies that they possess stable and unstable invariant manifolds that are impenetrable barriers and mediate transport in phase space. These invariant manifolds of the NHIM are structurally stable, that is, stable under perturbation Wiggins (2013). For two DoF systems, the NHIM is an unstable PO, and for an $n>2$ DoF system at a fixed energy, the NHIM has the topology of a $(2n-3)$-dimensional sphere and is the equator of a $(2n-2)$-dimensional sphere which constitutes the DS. The DS can be used to divide the $(2n-1)$-dimensional energy surface into two parts, reactants and productsGillilan and Ezra (1991); Komatsuzaki and Nagaoka (1996); Komatsuzaki and Nagoaka (1997); Komatsuzaki and Berry (2000, 2002). An elementary description of the role of the NHIM in reaction dynamics is given in Wiggins (2016a). Fundamental theorems assure the existence of the phase space structures —NHIM and its invariant manifolds —for a range of energies above that of the saddle Wiggins (2013). However, the precise extent of this range, as well as the nature and consequences of any bifurcations of the phase space structures that might occur as energy is increased, is not known and is a topic of continuing researchLi et al. (2009); Inarrea et al. (2011); Allahem and Bartsch (2012); Mauguière et al. (2013); MacKay and Strub (2014, 2015).

Motivation.

In this article, we investigate the changes in geometry of the phase space structures in a normal form Hamiltonian that undergoes saddle-node bifurcation. In general, a Hamiltonian system with a saddle equilibrium point whose eigenvalues depend on the potential well depth parameter will exhibit bifurcation when the potential energy barrier height decreases. If the potential energy surface has a single well, then the result is a collision of the stable and saddle equilibrium points as we illustrate in Figs. 3 and 5. This leads to drastic changes in the geometry of phase space structures — NHIM and its invariant manifolds — which can then be used to account for corrections to Kramers’ reaction rate as barrier height decreases Hathcock and Sethna (2019). Furthermore, quantifying rates of crossing low or vanishing barrier is significant for experimental study of single bond dynamics of molecules and control of micro and nano-electromechanical devices Hathcock and Sethna (2019); Husson et al. (2009); Miller and Shaw (2012); Herbert and Bouchet (2017). We would also like to point out that the questions addressed in this article are further motivated by the work of Borondo and co-authors Borondo et al. (1995, 1996); Revuelta et al. (2019) who noted the significance of a saddle-node bifurcation in the isomerization of LiCN/LiNC molecule.

Conservative dynamics on an open potential well has received considerable attention because the geometry of phase space structures explains the intricate fractal structure of ionization rates Mitchell et al. (2003a, b, 2004a). Furthermore, the discrepancies in observed and predicted ionization rates in atomic systems has also been explained by accounting for the topology of these phase space structures. These have been connected with the breakdown of the ergodic assumption that is the basis for using ionization and dissociation rate formulae De Leon and Berne (1981). This rich literature on chaotic escape of electrons from atoms sets a precedent for delineating the changes in phase space structures due to saddle-node bifurcation that can be expected in open potential wells where the eigenvalues of the saddle equilibrium point depend on the well depth parameter Mitchell et al. (2004b, a); Mitchell and Ilan (2009); Mitchell and Delos (2007); Wang et al. (2010).

Approach in this article.

We first develop a Hamiltonian model that exhibits saddle-node bifurcation when the well depth is decreased and obtain a relationship between the parameters of the model for the bifurcation in the 2 DoF system. Then, the NHIM (an unstable PO in the 2 DoF system) is computed using differential correction and numerical continuation, and its invariant manifolds are computed using globalization. The role of these invariant manifolds in reaction dynamics and the implications of the change in their geometry by varying well depth is discussed. We compare the results with another numerical method that can be used for identifying high dimensional phase space (4 or more dimensions) structures called Lagrangian descriptors (LDs). Lagrangian descriptors is a trajectory diagnostic method for revealing invariant manifolds that mediate phase space transport. The method was originally developed in the context of Lagrangian transport studies in fluid dynamics Madrid and Mancho (2009), but the wide applicability of the method has recently been recognized as a tool to construct the dividing surface via identification of NHIM Craven and Hernandez (2016, 2015); Craven et al. (2017); Revuelta et al. (2019); Junginger et al. (2017); Feldmaier et al. (2017); Patra and Keshavamurthy (2018). The method is straightforward to implement computationally, and it provides a high resolution method for locating high dimensional invariant manifolds using low dimensional slices of the high dimensional phase space Demian and Wiggins (2017); Naik et al. (2019); Naik and Wiggins (2019). Therefore, this technique provides us with a first step towards realizing a complete “phase space tomography” of high dimensional invariant manifolds. It also applies to both Hamiltonian and non-Hamiltonian systems Lopesino et al. (2017) as well as to stochastic dynamical systems Balibrea-Iniesta et al. (2016). Moreover, Lagrangian descriptors can be applied directly to data sets obtained from the numerical solution of geophysical models and satelite observations Garcia-Garrido et al. (2016); Ramos et al. (2018); Balibrea-Iniesta et al. (2019). We provide more details on our use of Lagrangian descriptors in Appendix A.

This paper is outlined as follows. Section II is devoted to the description of the Hamiltonian models with one and two DoF that we will use in our study of saddle-node bifurcation phenomena in phase space. In Section III we explore the implications of these model problems for reaction dynamics, and discuss our results by analyzing the phase space geometrical structures responsible for the escaping dynamics from the PES and their relevance for chemical reactions. Finally, in Section IV we summarize the conclusions of this work. The reader can find all the relevant information related to the method of Lagrangian Descriptors and its implementation details in Appendix A. Furthermore, Appendix C describes the numerical techniques used (differential correction and manifold globalization) for the computation of the NHIM and its stable and unstable manifolds associated to the rank-1 saddle point in the 2 DoF Hamiltonian system that characterizes the saddle-node phase space bifurcation.

II Hamiltonian Saddle-Node Bifurcation Models

In this section we introduce the normal form Hamiltonian that exhibit saddle-node bifurcation and their implications for reaction dynamics. First, we describe the normal form Hamiltonian system with one DoF, and parametrize this model by adding a parameter that controls the potential well depth. Next, we extend this model to a two DoF Hamiltonian system by coupling a bath mode (a harmonic oscillator) to the reaction DoF.

II.1 One Degree-of-Freedom Hamiltonian

The normal form for the one DoF Hamiltonian that undergoes a saddle-node bifurcation in phase spaceWiggins (2017) can be written as

[TABLE]

where $\mu\in\mathbb{R}$ is the bifurcation parameter. Hamilton’s equations are given by

[TABLE]

and the equilibria of this system are located at $(u,v)=(\pm\sqrt{\mu},0)$. We observe that the dynamical system has two different equilibrium points for $\mu>0$, which approach each other as the bifurcation parameter goes to zero. When $\mu=0$, both equilibria “collide” into one equilibrium point at $(0,0)$. For $\mu<0$, the resulting dynamical system has no equilibria.

In order to study the stability of the equilibria, we linearize the vector field and calculate the Jacobian

[TABLE]

Evaluating the Jacobian at the equilibrium point $(-\sqrt{\mu},0)$, we get

[TABLE]

which has eigenvalues $\pm\sqrt[4]{4\mu}$, and hence the equilibrium point $(-\sqrt{\mu},0)$ is a saddle. The Jacobian at $(\sqrt{\mu},0)$ is

[TABLE]

which has eigenvalues $\pm\sqrt[4]{4\,|\mu|}\,i$,yielding thst the equilibrium point $(\sqrt{\mu},0)$ is a center. For $\mu=0$, zero is the only eigenvalue of the Jacobian matrix, and therefore the linearization does not provide enough information, and one needs to include higher order terms to determine the stability of the equilibrium point $(0,0)$.

We compute LDs with the goal of detecting the invariant manifolds of the dynamical system given in Eq. (2). As shown by the phase portraits in Fig. 1, trajectories of (2) can escape off to infinity in finite time depending on the total energy of the system and the location of the initial conditions. Therefore, this issue needs to be accounted for in the computation of LDs as explained in Appendix A, and the approach we will follow to resolve this problem has also been adopted in the chemical reactions literature Craven and Hernandez (2015); Craven et al. (2017). In particular we will use the $p$-norm in the LD definition (59) with $p=1/2$. Therefore, we evolve initial conditions forward and backward in time for $\tau=8$ and whenever a trajectory leaves the domain defined by the circle of radius $15$ about the origin, we stop the numerical integration of that trajectory. In Fig. 1 we show the LD contour maps for three different values of the bifurcation parameter $\mu=-0.25,\,0,\,0.25$ and compare the results obtained by means of LDs with the corresponding phase portraits. We can clearly see in Fig. 1 how the method succesfully recovers all the relevant phase space structures. The stable and unstable manifolds are highlighted by the singularities present in the LD contour map. Furthermore, another interesting aspect to highlight from Fig. 1D is that, for $\mu=0$, the LD values in the neighborhood of the origin seem to indicate that a bifurcation is going to take place in the phase space structure, as one can observe that a ‘ghost’ center structure is about to be created close to the cusp at the origin.

Fixing the saddle point at the origin.

In the Hamiltonian saddle-node model (1), the equilibria move as the bifurcation parameter $\mu$ is varied. In order to simplify the bifurcation analysis necessary to address the implications for chemical reaction dynamics, we fix the saddle point at the origin. Moreover, this form of non-moving saddle will become useful when discussing the 2 DoF model of the saddle-node bifurcation. This transformation will also facilitate the construction of the dividing surface from the NHIM associated with the ‘non-moving’ rank-1 saddle.

Let us consider the linear change of coordinates

[TABLE]

and substitute this transformation into Eq. (2) to yield the following dynamical system in the new coordinates

[TABLE]

with corresponding Hamiltonian

[TABLE]

The equilibrium points of Eq. (7) are $(0,0)$ and $(2\sqrt{\mu},0)$. Their stability is characterized by the eigenvalues of the Jacobian matrix

[TABLE]

evaluated at each equilibrium point

[TABLE]

with eigenvalues $\pm\sqrt[4]{4\mu}$ and $\pm\sqrt[4]{4\,|\mu|}\,i$, respectively. Hence, the origin is a saddle and $(2\sqrt{\mu},0)$ is a center equilibrium point. This is also supported by the result that eigenvalues do not change under translational change of coordinates.

In Fig. 2, we show the variable integration time LD contour maps for Eq. (7) at two different values of the bifurcation parameter, $\mu=0$ and $\mu=0.25$. These plots illustrate how the stable and unstable manifolds of the saddle equilibrium point are identified by the singular features, points where LDs are non-differentiable, in the LD contour maps. These singular features can be visualized by taking one dimensional slices of the LD contour map, where jump discontinuities mark the initial conditions on an invariant manifold.

Introducing a parameter to control “depth” of the PES.

The depth of the well on the PES is controlled by the cubic term that appears in the Hamiltonian (8). Therefore, we introduce a parameter that allows us to vary the amplitude of this term, and consequently the strength of the nonlinearity in the vector field (7). Thus, the Hamiltonian becomes

[TABLE]

where $\mu\in\mathbb{R}^{+}\cup\{0\}$ and $\alpha\in\mathbb{R}^{+}$ are the two parameters, $T$ is the kinetic energy of the DoF and $V$ its potential energy. Hamilton’s equations are given by

[TABLE]

The equilibria are $(0,0)$ and $\left(2\sqrt{\mu}/\alpha,0\right)$, and the Jacobian of the vector field is given by

[TABLE]

We determine the stability of the equilibria by evaluating the Jacobian which gives

[TABLE]

The eigenvalues of the Jacobian at the equilibrium point $(0,0)$ are $\pm\sqrt[4]{4\mu}$ with corresponding eigenvectors $(1,\pm\sqrt[4]{4\mu})$. The eigenvalues for the Jacobian at the equilibrium point $\left(2\sqrt{\mu}/\alpha,0\right)$ are $\pm\sqrt[4]{4\mu}\,i$ with eigenvectors $(1,\pm\sqrt[4]{4\mu}\,i)$. Clearly, the eigenvalues of the Hamiltonian (8) are perserved, and $(0,0)$ is a saddle, while $\left(2\sqrt{\mu}/\alpha,0\right)$ is a center. We note here that the eigenvalues and eigenvectors of the Jacobian evaluated at the equilibrium points only depend on $\mu$.

“Depth” of the well on the PES, referred to as well-depth from here on, is determined by the difference between the potential energy of the saddle and the potential energy of the center (minimum of the well) equilibrium points. The potential energy function is given by

[TABLE]

and this difference is given by

[TABLE]

Hence, for a fixed $\mu$, the well-depth is increased by decreasing $\alpha$. In fact, if we denote $\mathcal{D}=2\sqrt{\mu}/\alpha$ as the distance between the saddle and the center equilibrium point, we get

[TABLE]

therefore, as $\alpha$ increases, the well-depth approaches zero faster than the distance between the saddle and the center equlibrium points. Furthermore, the rate at which the well-depth changes as we change the distance between equilibria is given by

[TABLE]

which is proportional to the product of the square of the eigefrequency associated to the saddle equilibrium and the distance between equilibria.

In Fig. 3(a) we illustrate the well-depth together with the distance between the saddle and the center. In Fig. 3(b) shows how the potential energy function changes, for a fixed value $\mu=1$, as the well-depth parameter $\alpha$ (dstrength of the nonlinearity) is increased.

Recall that, as we discussed in the introduction, in the context of chemical reaction dynamics the escape from a potential well problem can be identified for instance with dissociation or fragmentation reactions, where a chemical transformation takes place if a bond of a molecule A breaks up, giving rise to two products B and C. In this situation, the equilibrium conformation of the given molecule A is represented by a potential well in a PES, and dissociation into B and C takes place if the system has sufficient energy to croos the potential barrier that separates bounded (vibration) from unbounded (bond breakup) motion. This setting can be modeled for instance by the potential energy function in Eq. (15), and is illustrated in Fig. 4(a), where a locally cubic potential energy function about the potential barrier describes this phenomenon, which could be the result of a saddle-node bifurcation that has occurred in the phase space of the Hamiltonian system. In Fig. 4(b) we show the phase portrait for the dynamical system (12) which shows the equilibrium points and the homoclinic orbit. Fig. 4(c-d) display the changes in the phase space structure using the variable integration time LDs for different values of the well-depth parameter $\alpha$.

II.2 Two degrees-of-freedom Hamiltonian

In this section we introduce the normal form for the two DoF Hamiltonian that undergoes a saddle-node bifurcation by extending the one DoF Hamiltonian discussed above. To do so, we add another degree of freedom in the form of a harmonic oscillator with mass $m=1$ and frequency $\omega$, known as a bath mode in the terminology of chemical reaction dynamics, and parametrize the coupling of the reaction and bath modes. Thus, the Hamiltonian becomes

[TABLE]

where $\alpha>0$ is the well-depth parameter and measures the strength of the nonlinearity, $\omega>0$ is the frequency of the harmonic oscillator or the bath mode, $\mu\geq 0$ is the bifurcation parameter, and $\varepsilon\geqslant 0$ is the strength of the coupling between the reaction and the bath mode. Identifying this Hamiltonian’s kinetic energy, $T(p,p_{x})$, potential energy, $V(q,x)$, and assuming $H(q,x,p,p_{x})=T(p,p_{x})+V(q,x)$, we get that

[TABLE]

The corresponding Hamilton’s equations are given by:

[TABLE]

For this Hamiltonian system, the phase space is four dimensional, and since energy is conserved, the trajectories evolve on a three dimensional energy surface. The equilibria for this system are located at $\mathbf{x}_{1}^{e}=(0,0,0,0)$ and $\mathbf{x}_{2}^{e}=\left(q_{e},x_{e},0,0\right)$ where

[TABLE]

The energy of the system at the equilibrium points is:

[TABLE]

One can easily show that there exists a critical value of the coupling strength (also interpreted as a perturbation to the one DoF model) given by

[TABLE]

for which the Hamiltonian (21) has only one equilibrium point at the origin. We note here that this critical value is independent of the nonlinearity strength parameter $\alpha$ and requires $\omega^{2}>2\sqrt{\mu}$ to be satisfied. Moreover, it is characterized by a functional relationship between the squares of the eigenfrequencies of the reactive and bath modes of the uncoupled system ($\varepsilon=0$). We will analyze the influence of the perturbation strength $\varepsilon$ on the geometry of the phase space structures up until this critical condition in Section III.

We illustrate in Fig. 5 how the geometry of the PES (20) and the equipotentials in configuration space change as we vary the coupling strength (perturbation parameter) $\varepsilon$. For this visualization, we have used the following values of the model parameters $\mu=0.25$, $\alpha=2$, and $\omega=1.25$. The reason for doing so is that they satisfy the condition $\omega^{2}>2\sqrt{\mu}$ and consequently a critical value of the perturbation parameter exists, given by $\varepsilon_{c}=25/9$, at which the two equilibrium points ‘collide’ at the origin resulting in a saddle-node bifurcation. The dynamics after this collision is beyond the scope of this study and therefore we will only focus on the description of the system dynamics for different values of the coupling strength when $\varepsilon<\varepsilon_{c}$. We observe in Fig. 5(a,b) that the location of the well when the DoF are uncoupled $(\varepsilon=0)$ lies on the $q$-axis. As $\varepsilon$ is increased, Fig. 5(c,d) show that its effect is to tilt the PES with respect to the configuration space plane, and the position of the center equilibrium point moves off the $q$-axis towards the origin. Finally, the situation for which the coupling strength reaches the critical value $\varepsilon_{c}$ is shown in Fig. 5(e,f) when the collision has happened and there is only one equilibrium point at the origin.

Linear stability analysis of the equilibrium point at the origin.

We show here that the equilibrium point $\mathbf{x}_{1}^{e}$ is a rank-1 saddle on the PES and that $\mathbf{x}_{2}^{e}$ is a center corresponding to the bottom of the potential well. This linear stability analysis will be useful later for computing the rank-1 saddle’s NHIM, which carries its influence to higher energies, and also to determine the NHIM’s invariant stable and unstable manifolds that act as codimension one impenetrable barriers on the constant energy surface. The global geometry of the invariant manifolds is paramount for quantifying reaction rates such as escape rate from potential well.

We determine the stability of the equilibria by linearizing Eq. (21), which gives the Jacobian matrix

[TABLE]

The stability of $\mathbf{x}_{1}^{e}=(0,0,0,0)$ is given by the eigenvalues of the Jacobian

[TABLE]

which has the characteristic equation

[TABLE]

where $\beta$ are the eigenvalues and $\mathbb{I}$ the identity matrix.

When the reaction and bath modes are decoupled, that is $\varepsilon=0$ in (27), the eigenvalues are given by $\pm\lambda_{0}$ and $\pm\omega_{0}\,i$, where

[TABLE]

Therefore, the origin is a rank-1 saddle equilibrium point, since the linearized system has exactly one pair of real eigenvalues, $\lambda_{0}$ and $-\lambda_{0}$ and the saddle plane is spanned by their eigenvectors. We know from the Moser’s generalization of the Lyapunov Subcenter TheoremWiggins (2013, 2003) that when the energy of the system is above that of the rank-1 saddle, there is a two dimensional plane spanned by the eigenvectors of $\pm\omega_{0}\,i$, known as the center invariant manifold. This invariant manifold with normal hyperbolicity is referred to as the normally hyperbolic invariant manifold (NHIM) and has the topology of a $(2N-3)$-sphere, that is $S^{2N-3}$. In the two DoF setting, a NHIM is simply an unstable periodic orbit whose geometry is topologically equivalent to a circle, i.e. $S^{1}$.

As the energy of the system is increased above the energy of the rank-1 saddle equilibrium point, a bottleneck opens in the phase space connecting dynamically different phase space regions and allowing trajectories to move between them. This phenomenon results in a phase space transport mechanism. This framework of understanding chemical reactions is realized by computing the stable and unstable manifolds associated with the unstable periodic orbit. These invariant manifolds have cylindrical geometry, that is $\mathbb{R}\times\S^{1}$, and their global behavior is referred to as tube dynamics. The cylindrical manifolds on the constant energy surface are impenetrable barriers (since they are two dimensional on the three dimensional energy surface) separating the reactive and non-reactive trajectories in the phase space. Thus, they determine the initial conditions that will pass through the bottleneck in some future time (or had passed through it in the past) during their evolution Wiggins et al. (2001b). The unstable periodic orbit provides us with the scaffolding to construct the dividing surface that separates the trapped motion in the well region of the PES and the escape to infinity of particle trajectories through the phase space bottleneck. The local (linearized) dynamics mediated by these phase space structures is shown in Appendix B. We will resume the discussion on computing the NHIM and its invariant manifolds in §:III.2.

When $\varepsilon\neq 0$, the reaction and bath modes are coupled and the Hamilton’s equations (21) are non-integrable. This means we can expect chaotic trajectories to appear that might lead to reaction by escaping the potential well. However, the geometry of the NHIM and its invariant manifolds still governs this reacting and non-reacting behavior. To simplify the eigenvalue problem arising in the stability analysis, we observe that Eq. (27) can be rewritten in terms of the eigenvalues of the uncoupled system as

[TABLE]

Introducing $\xi=\beta^{2}$ the solutions are

[TABLE]

We note here that the two possible values of $\xi$ have opposite signs for all values $\varepsilon>0$ of the coupling strength when $\omega_{0}^{2}\leq\lambda_{0}^{2}$, and for $0<\varepsilon<\varepsilon_{c}$ whenever the condition $\omega_{0}^{2}>\lambda_{0}^{2}$ is satisfied, where the critical coupling strength $\varepsilon_{c}$ is given by Eq. (24). If we denote $\xi_{1,\varepsilon}$ and $\xi_{2,\varepsilon}$ as the positive and negative roots respectively in Eq. (30), then the eigenvalues of the Jacobian matrix are $\pm\lambda_{\varepsilon}$ and $\pm\omega_{\varepsilon}i$ where

[TABLE]

This shows that in this case the equilibrium point at the origin is a rank-1 saddle. In particular, when the coupling between the reaction and the bath mode is weak, that is $\varepsilon\ll 1$, we have

[TABLE]

which confirms that the rank-1 saddle structure of the equilibrium point at the origin persists under small perturbations.

We finish the linear stability analysis of the equilibrium point at the origin with the computation of the eigenvectors. Let us denote the eigenvector of the matrix $J(\mathbf{x}_{1}^{e})$ by $\mathbf{v}=\left[v_{1},v_{2},v_{3},v_{4}\right]^{T}$, then the eigenvalue $\beta$ satisfies

[TABLE]

The eigenvectors associated with the real eigenvalues $\beta=\pm\lambda_{\varepsilon}$ correspond to the tangent directions of the unstable and stable manifolds, respectively, of the NHIM at the rank-1 saddle, and are given by

[TABLE]

The eigenvectors associated with the complex eigenvalues $\beta=\pm\omega_{\varepsilon}i$ span the center subspace of the NHIM at the rank-1 saddle and are given by

[TABLE]

where $\omega_{\varepsilon}$ is the magnitude of the complex eigenvalue in (31). Therefore, we can write the solution to the linearized system at the rank-1 saddle as

[TABLE]

where $C_{1},C_{2}\in\mathbb{R}$ and $\eta=\eta_{1}+\eta_{2}\,i\in\mathbb{C}$ are constants to be determined from an initial condition. We will use this general solution to the linear system to select an initial guess to search for the NHIM using the differential correction and numerical continuation method Koon et al. (2011).

For the sake of completeness, we show next the linear stability analysis for the equilibrium point $\mathbf{x}_{2}^{e}$. The Jacobian is given by

[TABLE]

The characteristic equation $\det\left(J(\mathbf{x}_{2}^{e})-\beta\,\mathbb{I}\right)=0$ becomes

[TABLE]

Introducing $\xi=\beta^{2}$ the solutions are

[TABLE]

It can be easily shown that both solutions are negative, and consequently $\mathbb{J}(\mathbf{x}_{2}^{e})$ has two pairs of complex eigenvalues. This implies that the equilibrium point $\mathbf{x}_{2}^{e}$ is a center.

III Results and Discussion: Influence of well-depth on reaction dynamics

In this section we discuss the implications of changing the well-depth in terms of the geometry of the invariant manifolds that mediate reacting trajectories (or escape from the potential well). This is done by identifying the invariant manifolds using Lagrangian descriptors and a numerical continuation with globalization method (see Appendix C). The geometry of these invariant manifolds and character of changes affected by varying the potential well-depth are involved in the computation of reaction fraction (or escape rates).

III.1 One Degree-of-Freedom Hamiltonian

Here we return to the 1 DoF Hamiltonian (11) where the reaction is defined as the change in the sign of the configuration coordinate $q$. For the rank-1 saddle equilibrium point at $(0,0)$, the energy is $H(0,0)=H_{c}=0$ which we call the critical energy, and that of the center equilibrium point at $(2\sqrt{\mu}/\alpha,0)$, at the bottom of the potential well, is $H(2\sqrt{\mu}/\alpha,0)=H_{w}=-4\sqrt{\mu^{3}}/3\alpha^{2}$. In this model system, the phase space is two-dimensional and since the total energy is conserved, the trajectories evolve on the one-dimensional constant energy lines or isoenergetic contours of the Hamiltonian in the 2D phase space. Therefore, the system is completely integrable, and for a fixed energy $H_{0}$ the trajectories evolve on the one-dimensional curve

[TABLE]

The value of total energy $H_{0}$ with respect to $H_{c}$ and $H_{w}$ characterizes the nature of trajectories as follows We discuss the nature of the trajectories depending on the energy $H_{0}$ of the system:

• •

Case $H_{0}<H_{w}$: The initial conditions that satisfy $q<-\sqrt{\mu}/\alpha$, which lie to the left of the potential energy barrier at the origin. Initial conditions in this configuration space initially climb the potential and at some point their velocity reverses direction and they roll down the potential going to infinity. In the phase space, they come from infinity and fly-by the rank-1 saddle on their way towards infinity. They can not cross the barrier and thus do not lead to reaction shown by the blue curve in Fig. 6.

• •

Case $H_{w}\leq H_{0}<0$: Now two types of trajectories are possible. For the initial conditions that satisfy $-\sqrt{\mu}/\alpha\leq q<0$, trajectories will show same fly-by behavior as in the previous case. If initial conditions satisfy $0<q_{0}<3\sqrt{\mu}/\alpha$, then they are confined in the potential well on the right of the potential barrier. Since they are bounded by the homoclinic orbit, trajectories will never escape the well, and thus do not lead to reaction and shown as green curves in Fig. 6.

• •

Case $H_{0}=H_{c}=0$: When the total energy of the initial conditions is equal to the potential energy of the barrier, resulting trajectories approach the barrier asymptotically in forward and backward time. The initial conditions that start on the left of the potential barrier and asymptotically approach it in forward and backward time form pieces of the stable and unstable manifolds of the saddle equilibrium point. The initial conditions that start on the right of the potential barrier also approach the barrier as time goes to infinity. These trajectories combine to form the homoclinic orbit and shown as a black curves in Fig. 6.

• •

Case $H_{0}>H_{c}=0$: The energy of the system is above that of the barrier so trajectories can escape from the well, that is the configuration space coordinate $q$ can change sign and lead to reaction. The dividing surface that trajectories must cross once in a given directionWiggins (2016a) (locally no recrossing), that is when escaping from the potential well, is given by

[TABLE]

We note that in this case the dividing surface for a fixed energy and given direction of crossing is a point, or has a geometry of $\mathbb{S}^{0}$, thus it can partition the one-dimensional (isoenergetic) curves into “reactant” and “‘product” regions. To further characterize the reaction dynamics, we resort to the linear stability analysis near the saddle equilibrium point at the origin. The eigenvalues of the saddle equilibrium point are $\pm\sqrt[4]{4\mu}$ and the corresponding eigenvectors $\left(\pm 1/\sqrt[4]{4\mu},1\right)$. Therefore, the origin is a NHIM since the dynamics is hyperbolic in directions normal to it Wiggins (2016b). Furthermore, we remark that the stability of the NHIM only depends on the $\mu$ parameter. The eigenvectors tell us how the unstable and stable manifolds of the NHIM are oriented and can be used to numerically globalize the linear approximation. However, we can compute the stable and unstable manifolds of the saddle equilibrium point analytically by noting that they lie on the zero level curve of the Hamiltonian (total energy). Thus, the stable and unstable manifolds are given by

[TABLE]

where $\Gamma$ is the homoclinic orbit

[TABLE]

and

[TABLE]

We note here that the geometry of the invariant manifolds is $\mathbb{R}\times\mathbb{S}^{0}$, that is codimension-1 in $\mathbb{R}^{2}$, and hence can form the impenetrable barrier between the reactive and non-reactive trajectories as shown in Fig. 6.

III.2 Two degree-of-freedom Hamiltonian

We return to the 2 DoF Hamiltonian (19) which introduces the bath degree-of-freedom as a harmonic oscillator into the reaction dynamics. We show the changes in the geometry of the invariant manifolds due to the changes in the well-depth which leads to a saddle-node bifurcation. This is done using Lagrangian descriptors along with numerical continuation and globalization for computing the NHIM and its invariant manifolds.

Let us consider a fixed energy $H_{0}$ and since the model has 2 DoF we know that the dynamics is on a three-dimensional energy surface given by

[TABLE]

The projection of the energy surface onto the $(q,x)$ configuration space is the region of energetically possible motion for a fixed energy $H_{0}$, and is given by

[TABLE]

This projection denotes configurations with positive kinetic energy and has been known in classical mechanics as the Hill’s region. The boundary of $M(H_{0})$ is defined as the locus of points in the $(q,x)$ plane where the kinetic energy is zero, that is $(H_{0}-V(x,y))=0$, and is called the zero velocity curve. The trajectories are only able to move on the side of the curve where the kinetic energy is positive, shown as white regions in Fig. 7(b-d).

To identify the invariant manifolds, we take two dimensional slices of the energy surface and determine the intersection of the invariant manifolds with these low-dimensional slices. In particular, we calculate Lagrangian desriptor and compare with qualitative understandings from Poincaré surface-of-section. The isoenergetic surfaces-of-section are

[TABLE]

where $q_{e}$ is the configuration space coordinate of the equilibrium point at the bottom of the well on the PES (22).

We return to the Hamiltonian (19) where the “reaction” and “bath” DoF are uncoupled, that is $\varepsilon=0$. Therefore, the system is integrable and the trajectories are regular, thus the separable Hamiltonian is

[TABLE]

where $H_{r}$ is the Hamiltonian for the reaction and $H_{b}$ is the Hamiltonian for the bath DoF.

Given a fixed total energy of the system $H_{0}$, the necessary condition for the reaction to take place is when the total energy is above that of the barrier of the PES located at the origin, which is zero. For $H_{0}\leq 0$ the energy surface divides phase space into two disconnected regions as illustrated in Fig. 7 so that we have bounded motion in the potential well region. In Fig. 8, we compare the bounded trajectories (regular dynamics) of the system for $H_{0}=0$ by computing LDs and Poincaré section on the surface-of-section, $\mathcal{U}_{qp}^{+}$ (47). Both methods clearly recover, the trajectories on tori, as known for integrable Hamiltonian systems Meyer et al. (2009), that fills the energy surface. We observe that the high values (white regions in 8) in the LD contour map recovers the quasiperiodic trajectories, which is a consequence of the relationship between the convergence of time averages of LDs with the Ergodic Partition Theorem as explained in Appendix A.

When the energy of the system is above the barrier, that is $H_{0}>0$, the topology of the energy surface changes and a phase space bottleneck opens up in the barrier region allowing the escape from the well as shown in Fig. 7(d). Thus, reaction can take place by crossing the bottleneck and passing from the $q>0$ to $q<0$ or vice versa. Therefore, $q=0$ is a natural choice for defining a dividing surface (DS) separating reactants (bounded motion in the potential well region) from products (escape to infinity) in phase space or vice versa. Hence, the isoenergetic DS is given by

[TABLE]

which has the geometry of a 2-sphere, that is $\mathbb{S}^{2}$ in $\mathbb{R}^{3}$. To be precise, it is an ellipsoid with semi-major axis $\sqrt{2H_{0}}$ in the $p$, $p_{x}$ axis, and $\sqrt{2H_{0}}/\omega$ in the $x$ axis. This ellipsoid has two hemispheres, known as the forward and backward DS with the form

[TABLE]

Forward “reaction” occurs when trajectories cross $\mathcal{D}_{f}$ and back “reaction” when trajectories cross $\mathcal{D}_{b}$. The forward and backward DS is joined at the equator along the normally hyperbolic invariant manifold given by

[TABLE]

which has the topology of $S^{1}$. To be precise, it is an ellipse with semiaxis $\sqrt{2H_{0}}$ in the $p_{x}$ direction and $\sqrt{2H_{0}}/\omega$ in the $x$ direction. As discussed earlier, for a 2 DoF Hamiltonian the NHIM is an unstable periodic orbit which extends the influence of the rank-1 saddle equilibrium point in the bottleneck of the PES (a configuration space concept) into phase space. We note that for $H_{0}>0$ the NHIM is the correct phase space structure that anchors the barriers to the reaction and carry the effect of the saddle equilibrium point to a range of energies as given by Moser’s generalization of Lyapunov Subcenter Manifold Theorem Wiggins (2003). The stable and unstable manifolds of the unstable periodic orbit are

[TABLE]

where $\Gamma$ is the homoclinic orbit

[TABLE]

and the left branches are

[TABLE]

We note that the stable and unstable invariant manifolds have the structure of a cartesian product of a curve in the $(q,p)$ saddle space and an ellipse in the $(x,p_{x})$ center space, and thus become cylindrical (or tube) manifolds. We show the energy surface in Fig. 7(a) at 3 energy values for $\epsilon=0$ where the bottleneck only opens for $H_{0}>0$. The NHIM computed using differential correction and continuation, and its invariant manifolds computed using globalization are shown in Fig. 9(a) for $H_{0}=0.05$ and $\epsilon=0$.

When $\varepsilon\neq 0$, the “reaction” and “bath” modes are coupled and for $\varepsilon$ small, we can think about the resulting dynamics as a perturbation of the uncoupled case . Given an energy of the system below that of the barrier, that is $H_{0}\leq 0$, the phase space bottleneck is closed and trajectories are trapped in the potential well region. However, due to the perturbation, the unstable regular motion on the tori in the uncoupled system is destroyed as given by the KAM theorem. Therefore, chaotic motion arises in some regions of the phase space due to a non-zero coupling of the “reaction” and “bath” modes. We illustrate the system’s behavior for the energy $H_{0}=0$ (energy of the rank-1 saddle) in Fig. 10 using LDs and Poincaré maps on the surfaces of section $\mathcal{U}^{+}_{qp}$ (47) and $\mathcal{U}^{+}_{xp_{x}}$ (48). We observe that there is a distinct correlation between the qualitative dynamics revealed by the Lagrangian descriptor (LD) contour maps and Poincaré sections. That is, chaotic regions of phase space which appear as a sea of points in the Poincaré section and hide the underlying structures of stable and unstable manifolds are completely resolved by LD contour maps where the tangled geometry of the manifolds is revealed by the points where the function attains a local minimum as has also been shown in Naik and Wiggins (2019). The capability of LDs to identify the phase space structures relevant in chemical reaction dynamics can also be found in recent literature Craven and Hernandez (2016); Demian and Wiggins (2017); Naik et al. (2019). In addition, the generation of Poincaré section relies on tracking the crossing of a 2D surface which can not be guaranteed in high dimensional phase space, while LDs just accumulate a positive scalar quantity along the trajectory and thus have potential to reveal high dimensional phase space structures.

Now let us consider the dynamics when the total energy is above that of the rank-1 saddle, that is $H_{0}>0$, and in particular we will set $H_{0}=0.05$ in our analysis. In this situation, a phase space bottleneck opens up in the barrier region, as shown in Fig. 7. In order to describe the structures that mediate reaction dynamics, that is the NHIM (or unstable periodic orbit in this case) and its stable and unstable manifolds, we use differential correction and continuation along with globalization as described in Appendix C. In Fig. 9(c-d), we have shown the cylindrical manifolds along with the energy surface and the unstable periodic orbit for a fixed energy. In order to recover the homoclinic tangle geometry of the invariant manifolds, we calculate LDs on the surfaces of section $\mathcal{U}^{+}_{qp}$ and $\mathcal{U}^{+}_{xp_{x}}$, and compare with the direct numerical construction of these invariant manifolds. This LD based diagnostic is similar to performing a “phase space tomography” of the high dimensional phase space structures using a low dimensional slice. In Fig. 11 we show the computation of variable time LD for an integration time $\tau=10$ on the slice $\mathcal{U}^{+}_{qp}$. We observe that LD clearly identify, by means of points which attain a minimum, the location of the stable and unstable manifolds and the NHIM at the intersection of the invariant manifolds. Since we are using a small integration time of $\tau=10$ to compute LDs in Fig. 11, the complete geometry of the homoclinic tangle is not fully revealed.

Therefore, in order to recover a more complete and intricate dynamical picture of the homoclinic tangle, the integration time to compute LDs has to be increased. This is shown in Fig. 12(b), where $\tau=30$ and we observe that the regular motion obtained in the middle of the Poincaré section displayed in Fig. 12(a) corresponds to trajectories that remain trapped in the potential well region and never escape. The trapped dynamics on the tori are also captured by the LD contour map shown in 12(b) which also reveals the homoclinic tangle of the stable and unstable manifolds and the resulting lobe dynamics Beigie and Wiggins (1992). We note that in all these computations we are using the variable time definition of LDs, since the open potential surface causes trajectories to escape to infinity through the bottleneck in finite time and resulting in NaN values in the LD contour map. This issue is discussed in Appendix A and illustrated in Fig. 14, and would hide the important underlying phase space structures making the interpretation of results difficult. In addition, when we analyze the dynamics using Poincaré sections, we can not ensure that trajectories return to the surface of section when escaping to infinity is possible. This will result in blank regions in the Poincaré sections as shown in Figs. 12(a) and (c). However, we note that trapped trajectories in the potential well corresponding to regular (motion on the tori) and chaotic motion are highlighted as expected in the Poincaré sections. These trajectories are non-reactive and will remain so until they satisfy the sufficient condition for reaction which is entering the cylindrical manifolds.

Identifying the regions inside the cylindrical (tube) stable and unstable manifolds can be done using globalization or using the LD contour map on appropriate isoenergetic surfaces. We compute LD contour maps to recover the intersections of these tube manifolds with the surface of section (48). The intersection of the invariant manifolds with the isoenergetic surface of section yields topological ellipses, known in the chemical reactions literature as reactive islands De Leon et al. (1991a, b); De Leon (1992); Patra and Keshavamurthy (2018) and is involved in calculating reaction rates/fraction. In order to illustrate the capability of LDs to recover the reactive island structure, we have compared the LD contour map and manifold intersection with the surface of section $\mathcal{U}^{+}_{xp_{x}}$ obtained using globalization in Fig. 13. We have shown the first intersection of the stable and unstable manifold with the surface of section as a blue/red curve superimposed on the LD contour map. We observe that the minima in the LD contour maps and the manifolds intersections are in agreement, thus verifying the LD based identification of reaction islands. We also observe that successive fold and resulting intersection of the tube manifolds with the surface of section are also revealed by the points with minima in the LD values. We note here that the first intersection of the manifolds encloses a large phase space volume in the potential well which indicates that a large portion of the potential well escapes to infinity through the bottleneck. This is despite the fact that we have chosen a very small value for the energy of the system, $H_{0}=0.05$, compared to the energy of the barrier. This is a consequence of using a high coupling strength, $\varepsilon=0.25$, which makes the vase-like shaped potential well region small and narrow as the center equilibrium point at the bottom of the potential well approaches the rank-1 saddle equilibrium point located at the origin as $\varepsilon$ is increased. Movies illustrating the change in both the shape of the energy surface and the equipotentials in configurations space as we vary the coupling strength can be found at here and here, respectively. We can see that the effect of increasing the coupling strength from zero is to tilt and squeeze the vase-like shape of the energy surface that qualitatively increases the number of reactive trajectories. This action of tilting and squeezing of the vase-like container (boundary defined by the total energy surface) is to pour out its contents (the reactive trajectories) and hence the increase in reaction fraction. A quantitative investigation of this phenomenon is current work in progress and beyond the scope of this article.

IV Conclusions

In this article we have presented new results, by means of using a normal form Hamiltonian that models a saddle-node bifurcation, to investigate the influence on phase space dynamics of decreasing the potential well depth of a PES. The decrease in potential well depth causes the saddle and center equilibrium points to collide —an effect that takes place in the configuration space. The resulting saddle-node bifurcation, when the well depth is decreased, manifests in the phase space as tilting and squeezing of the energy surface. This change in geometry of the energy surface due to decreasing well depth leads to more reacting (or escaping) trajectories for a given energy. This is supported by the increase in the area of the first intersection of the invariant manifolds with an appropriate surface of section. The trajectory diagnostic method of Lagrangian descriptors also identifies these changes in the geometry of the phase space structures. This method also provides an approach for revealing the influence of potential well depth on high dimensional phase space structures. Detecting the qualitative changes in the geometry of the invariant manifolds using Lagrangian descriptors gets us closer to the goal of achieving a complete high-dimensional “phase space tomography” for realistic molecular systems.

Acknowledgments

We acknowledge the support of EPSRC Grant No. EP/P021123/1 and Office of Naval Research Grant No. N00014-01-1-0769. The authors would like to acknowledge the EPSRC sponsored Chemistry and Mathematics in Phase Space (CHAMPS) workshop held March 19, 2019 at the University of Bristol where they had the opportunity to discuss with Florentino Borondo about the role of the saddle-node bifurcation in the isomerization of lithium cyanide.

Appendix A Lagrangian Descriptors

The computational tool that we use in this work to explore the template of geometrical structures governing phase space transport is the method of Lagrangian descriptors (LDs). This mathematical technique is a trajectory-based diagnostic that was originally developed in the context of Lagrangian transport studies in geophysical fluid dynamicsMadrid and Mancho (2009); Mancho et al. (2013). The fundamental idea behind this methodology is to integrate a positive scalar function along particle trajectories of a dynamical system of general time-dependence in the form:

[TABLE]

where the vector field $\mathbf{v}(\mathbf{x},t)\in C^{r}(r\geq 1)$ in $x$ and continuous in time.

Lagrangian Descriptors were first introduced in Madrid and Mancho (2009); Mancho et al. (2013) by means of a scalar function, referred to as the function $M$, to identify distinguished hyperbolic trajectories, i.e. moving saddles, of a dynamical system with general time dependence. In this original approach, the function $M$ was based on the computation of the arclength of a trajectory starting at an initial condition $\mathbf{x}(t_{0})=\mathbf{x}_{0}$ as it evolves forward and backward for a specified time period $\tau>0$. In this context, LDs were defined as:

[TABLE]

where $\|\cdot\|$ stands for the Euclidean distance and $\mathbf{v}(\mathbf{x},t)$ is the vector field of the dynamical system defined in Eq. (56). The connection between the function $M$ and invariant manifolds in phase space has been primarily demonstrated through numerical simulation experiments. Hyperbolic trajectories and their stable and unstable manifolds are revealed by the function $M$ through sharp changes, which we call singular features, in the values of the $M$ field where the gradient becomes very large and changes abruptly. Moreover, it is important to highlight that if the function $M$ is broken into forward and backward integration:

[TABLE]

where we have that:

[TABLE]

then $M^{(f)}$ detects the stable manifolds in phase space and $M^{(b)}$ does the same for unstable manifolds. Other structures such as tori-like invariant manifolds are highlighted due to the relation of the time averages of the function $M$ with the Ergodic Partition Theory Mezic and Wiggins (1999). This can be done by defining the time average:

[TABLE]

and analyzing its convergence as $\tau\to\infty$. A detailed description of the detection and visualization of phase space structures with the function $M$ can be found in Mancho et al. (2013); Lopesino et al. (2017). It is important to note here that in the context of chemical reaction dynamics, this definition of LDs would measure phase space arclength of trajectories. In particular, in the Transition State Theory literature Craven and Hernandez (2015); Junginger et al. (2016); Craven et al. (2017); Junginger et al. (2017), configuration space arclength has been used for the computation of LDs.

Recently, a rigorous mathematical foundation for Lagrangian descriptors has been established in Lopesino et al. (2017) by means of introducing an alternative definition based on the $p$-norm of the components of the vector field that define the dynamical system in Eq. (56). Consider the scalar function:

[TABLE]

For this LD one can mathematically prove in certain model problems that hyperbolic points and their stable and unstable manifolds are detected as singularities of the $M_{p}$ field, that is, points in which the function is non-differentiable. Therefore, this alternative definition provides a characterization for the concept of singular features. Moreover, tori-like invariant structures are also related to time averages of $M_{p}$ by means of the Ergodic Partition Theory. An important aspect to highlight from the alternative definition given in Eq. (59) is that it allows to decompose the phase space analysis by separating the integral into the different DoFs of the system under study, making it possible to isolate and assess their elliptic and hyperbolic dynamical contributions separately. This property has been shown to be relevant for the detection of unstable periodic orbits in Hamiltonian system with two DoF for the classical Hénon-Heiles systemDemian and Wiggins (2017). Furthermore, this property has been used recently to address the Barbanis system with three DoFNaik and Wiggins (2019), and also to mathematically prove that for the normal form of a Hamiltonian system with three DoFs with an rank-1 saddle, the NHIM and its stable and unstable manifolds are located at the singularities, which are also local minima, of the LD field. Hence, it is shown that if we decompose Eq. (59) into forward and backward integration as we also did for Eq. (57), for a sufficiently large integration time $\tau$ we have:

[TABLE]

where $\mathcal{W}^{u}$ and $\mathcal{W}^{s}$ are, respectively, the unstable and stable manifolds calculated at time $t_{0}$ and argmin denotes the phase space coordinates $\mathbf{x}_{0}$ that minimize the function $M_{p}$. In addition, the NHIM at time $t_{0}$ can be calculated as the intersection of the stable and unstable manifolds:

[TABLE]

At this stage, we would like to point out other relevant properties reated to the method of Lagrangian descriptors. First, the integration parameter $\tau$ that appears in the definition of both $M$ and $M_{p}$ plays a critical role in the identification of the underlying dynamical structures in phase space. For applications, the value of $\tau$ is chosen so that the dynamical history of particle trajectories covers all significant timescales for the problem under study. Therefore, there is no general rule on how to choose the appropriate value for $\tau$ in order to achieve the successful visualization of phase space structures. Consequently, all dynamical systems have to be investigated on a case by case basis by means of trial and error simulations. Observe that small values of $\tau$ would yield blurry phase space structures, which are barely recognizable, since this regime resembles the Eulerian (instantaneous) description of the flow. On the other hand, large $\tau$ values result in a richer and more complex geometrical description of phase space structures, because LDs resolve them in great detail, making the task of interpreting transport from the obained picture difficult. So there is always a compromise between these two situations. Another interesting feature of LDs is that its computational implementation is straightforward, even for high dimensions, and that the method can be parallelized easily to run on high performance computers, because it deals with initial conditions separately.

But probably the most important capability of LDs in all its variants is that it provides us with a high-resolution methodology to explore and visualize high-dimensional phase space dynamics. The successful results obtained by applying this tool to uncover the phase space geometry are in part vindicated from its naive approach of emphasizing initial conditions of particle trajectories rather than focusing on their precise location and long-term evolution. As a result of this crucial point, it offers tremendous advantages for the analysis of high-dimensional phase space, where the evolution of ensembles of initial conditions may yield trajectories that become disperse and “lost” with respect to each other, making the interpretation of phase space structures problematic and challenging. The goal of identifying phase space structures for high-dimensional systems with LDs is achieved by analyzing the behavior of initial conditions over low-dimensional phase space slices Demian and Wiggins (2017); García-Garrido et al. (2018); Naik et al. (2019); Naik and Wiggins (2019). In fact, any low-dimensional surface can be selected as a probe, which can be sampled with arbitrary high resolution by defining an adequate grid of initial conditions. Therefore, no resolution will be lost in the analysis, as the trajectories corresponding to these initial conditions evolve in time, since phase space structure is encoded in the initial conditions of the trajectories themselves.

In recent studiesJunginger et al. (2017); Naik and Wiggins (2019) it has been highlighted that computing fixed-time LDs, that is, integrating all initial conditions chosen on a phase space surface for the same integration time $\tau$, could give rise to some issues. First, it can obscure the detection of the NHIM, which is a crucial step, for instance in Transition State Theory, in order to determine chemical reaction rates by analyzing flux across the dividing surface constructed from the NHIM. This difficulty appears as a result of bounded trajectories recrossing the barrier region that surrounds the NHIM, due to the Poincaré recurrence theorem. Consequently, multiple minima and singularities occur in the LD plots, which makes the location of the true NHIM trajectories a challenging task. Another issue, which takes place in the saddle-node bifurcation problem that we study in this paper, is that some of the trajectories that escape the potential well of the PES can go to infinity in finite time. The trajectories that show this behavior will give NaN values in the LD scalar field, hiding some regions of the phase space, and therefore obscuring the detection of invariant manifolds. In order to illustrate this problem, we have calculated the fixed-time $M_{p}$ function given in Eq. (59) with $p=1/2$ for the 1 DoF Hamiltonian described by Eq. (11) using the model parameters $\mu=0.25$ and $\alpha=2$. The results are shown in Fig. 14 for different values of the integration time. In Fig. 14A, which corresponds to the value $\tau=3$, we can clearly see how LDs detect the hyperbolic fixed point at the origin and the elliptic point at $(1/2,0)$. Moreover, the stable and unstable manifolds that originate from the saddle point are barely visible, as well as the homoclinic orbit. This is a consequence of the integration time being small. However, flat regions on the left part of the phase space, corresponding to initial conditions that have escaped to infinity in finite time, start to appear as a result of NaN values in the LD. This effect is emphasized further in Fig. 14B, where we have calculated fixed-time LDs using $\tau=5$. As we increase the integration time more trajectories escape to infinity and, consequently, a larger region of phase space disappears in the LD picture.

In order to circumvent these issues, we will apply in this work the approach that has been recently adopted in the literature Junginger et al. (2017); Naik and Wiggins (2019) known as variable time Lagrangian Descriptors. In this methodology, LDs are calculated for a given initial condition until the trajectory leaves a certain barrier (or saddle) region $\mathcal{R}$ defined in the phase space. Therefore, the total integration time in this strategy depends on the initial condition itself, that is $\tau(\mathbf{x}_{0})$. In this formulation, the $p$-norm definition of LDs has the form:

[TABLE]

and, for a fixed integration time $\tau_{0}$, the total integration time is defined as:

[TABLE]

where $t^{+}$ and $t^{-}$ are the times for which the trajecory leaves the barrier region $\mathcal{R}$ in forward and backward time, respectively. Since the origin is an rank-1 sddle for the Hamiltonian models that we are using to analyze saddle-node bifurcations, we choose for the barrier region:

[TABLE]

To conclude, it is important to point out here that if the selected barrier region is large enough, the variable time LD definition given above in Eq. (62) will approach the fixed-time LD definition in Eq. (62). Normally hyperbolic invariant manifolds and their stable and unstable manifolds will be captured by the phase space points for which the LD is non-differentiable and tori-like structures can be determined from time averages of LDs. Moreover, if the barrier region is very small, the NHIM and its invariant stable and unstable manifolds will appear as local maxima in the LD field, while for large barrier regions the local minimum behavior given in Eqs. (60) and (61) is recovered. Consequently, the variable integration time LD provides us with a suitable methodology to study the phase space geometrical structures that characterize the type of saddle-node problems that we discuss in this work, since it avoids the issue of trahectories escaping from the potential well of the PES to infinity in finite time. A detailed analysis on the theoretical background on the variable integration time Lagrangian Descriptors technique will be carried out in future work.

Appendix B Visualization of phase space structure for the quadratic normal form Hamiltonian

Consider the normal form for a quadratic Hamiltonian system with two DoF given by

[TABLE]

where Hamilton’s equations are

[TABLE]

The equilibirum point is located at $(0,0,0,0)$ and has zero energy. It is trivial to check that the eigenvalues of the linearized system about the equilibrium point are $\pm\lambda$ and $\pm\omega\,i$, and hence the equilibrium point is of saddle$\times$center stability type, which is known as a rank-1 saddle. In this form, the Hamiltonian (65) is decoupled into the “reactive” mode given by $H_{r}$ and the “bath” mode given by $H_{b}$. For this reason, it is known as a separable quadratic Hamiltonian (SQH). This representation allows us to address the phase space structures and discuss the distribution of the total energy of the system between the two modes 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=H_{0}$ dividing surface (DS) can be defined by the $q_{1}=0$ hypersurface. The constant energy defines a three-dimensional energy surface in the four dimensional phase space given by

[TABLE]

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

[TABLE]

which is a two dimensional surface, and has the geometry $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. Waalkens and Wiggins (2004), which are given by

[TABLE]

The forward and backward DS join at $p_{1}=0$ giving

[TABLE]

which is of geometry $S^{1}$, that is a circle centered at the origin with radius $\sqrt{2H_{0}/\omega}$ in the $(q_{2},p_{2})$ plane. This is a normally hyperbolic invariant manifold (NHIM) associated with the rank-1 saddle and parametrized by total energy $H(q_{1},q_{2},p_{1},p_{2})=H_{0}$ Wiggins (2013). Invariance follows from the vector field (66), since when we have that $q_{1}=p_{1}=0$, this gives $\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})$ plane, have linear saddle-like dynamics. For a two DoF system, this NHIM is more commonly referred to in the literature as an unstable periodic orbit and is shown in Fig. 15.

In order to understand the relationship between the NHIM and the rank-1 saddle point, we note that for $H_{r}=H_{b}=0$ the NHIM reduces to the point $(q_{1},q_{2},p_{1},p_{2})=(0,0,0,0)$, which is the rank-1 saddle point on the energy surface $H_{0}=0$. Moreover, recall that this rank-1 saddle equilibrium point is a configuration space concept and is located on the potential energy surface. Therefore, as the total energy of the system is increased from [math], with the bath mode energy $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 rank-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 (70) are given by

[TABLE]

which are two-dimensional surfaces and have the geometry of $\mathbb{R}\times 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. 15.

Appendix C Computation of NHIM and its invariant manifolds for the 2 DoF system

In this aapendix we describe the steps followed in order to calculate the NHIM and its stable and unstable invariant manifolds associated with the rank-1 saddle located at the origin for the Hamiltonian system with 2 DoF given by Eq. (19), which models the saddle-node bifurcation phenomena in phase space.

Step 1: Select an excess energy above the critical value. The Lyapunov subcenter theoremWiggins (2014) tells us that when the energy of the Hamiltonian system is raised above that of the equilibrium point at the origin of saddle$\times$center (rank-1 saddle) stability, which is known as critical energy, a family of NHIM with geometry $S^{2N-3}$ bifurcates, where $N$ indicates the number of DoF. In our problem, $N=2$ so the topology of the NHIM is $S^{1}$, that is, the NHIM is an unstable periodic orbit (UPO). So we pick a value for the total energy $H_{0}$ greater than the critical value $H_{c}=0$, which gives an excess energy $\Delta H=H_{0}-H_{c}=H_{0}$. This total energy is also the energy of trajectories on the invariant manifolds which partition the phase space into nonreactive and reactive trajectories at the same energy. The excess energy can be an arbitrary value, up to the value at which the NHIM bifurcates.

Step 2: Obtain the NHIM at the selected excess energy. As we have discussed, for the 2 DoF problem the NHIM is an UPO which touches the equipotential contour corresponding to the selected energy at the bottleneck region. However, due to the saddle-like dynamics in the transverse directions, any numerical error in the computation of the periodic orbit will get exponentially amplified with each time step, which will eventually destroy the periodic nature of the trajectory, since it will fail to come back to its starting point.

In order to generate the desired unstable periodic orbit corresponding to a selected excess energy, we consider a procedure that starts with a small initial condition (“seed”) obtained from the linearized equations of motion near the rank-1 saddle, and uses differential correction and numerical continuationKoon et al. (2011); Naik and Ross (2017); Ross et al. (2018) on that initial guess. The result is an unstable periodic orbit at target energy $H_{0}$ of period $T$ which will be close to $2\pi/\omega_{\varepsilon}$, where $\pm\omega_{\varepsilon}i$ is the pair of imaginary eigenvalues of the linearization about the rank-1 saddle point. In order to choose an initial guess for the search of the UPO we can use the linearized equations of motion about the rank-1 saddle equilibrium point $\mathbf{x}_{1}^{e}=\mathbf{0}$. The Jacobian matrix at the origin is given by Eq. (26) and the general solution to the linearized dynamical system can be written as:

[TABLE]

where $C_{1},C_{2}\in\mathbb{R}$ and $\eta=\eta_{1}+\eta_{2}\,i\in\mathbb{C}$. The real eigenvalues $\pm\lambda_{\varepsilon}$ with corresponding eigenvectors $\mathbf{v}_{+},\mathbf{v}_{-}$ are described in Eqs. (31) and (34) and the pair of complex eigenvalues $\pm\omega_{\varepsilon}i$ with corresponding eigenvectors $\mathbf{w}_{+},\mathbf{w}_{-}$ are described in Eqs. (31) and (35). The idea is to use the complex eigenvalues and their corresponding eigenvectors to obtain a starting guess as an initial condition to initialize the method that searches for the UPO. To do so, we choose a small amplitude, $A_{x}\approx 2\times 10^{-5}$, periodic orbit in the center manifold of the linearized system by selecting in Eq. (73) the values $\eta=-A_{x}/2$ (this eliminates the factor 2 in the formula), $C_{1}=C_{2}=0$ and $t=0$. Thus, the initial guess is:

[TABLE]

which has a period of $T_{0,g}=2\pi/\omega_{\varepsilon}$.

Step 3: Differential correction of the initial guess. The differential correction procedure that we apply in this step to the initial guess will only yield a good approximation to the true UPO of the nonlinear Hamiltonian system whenever $A_{x}\ll 1$, which corresponds to a much smaller excess energy $\Delta H$ than the one we originally selected in step 1. The reason for this is that the construction of the initial guess is based on the linear approximation near the rank-1 saddle equilibrium point. The convergence criterion that we use is based on the basic property that a periodic orbit returns to the starting point after a given period $T$, which is in fact the period of the PO. If $\bar{\mathbf{x}}_{0}=\bar{\mathbf{x}}_{\rm po}(0)$ is a true initial condition on the PO $\mathbf{x}_{po}$ of period $T$, the convergence is checked using the condition:

[TABLE]

for some tolerance $\epsilon\ll 1$. In this approach, we hold the configuration coordinate $q$ constant, while applying correction to the $x$ configuration coordinate of the initial guess. While the $p_{x}$ momentum coordinate is used as a stopping criterion for the differential correction procedure, the $p$ momentum coordinate is used as a terminating event (crossing the $p=0$ plane in phase space) for the integration. It is to be noted that this combination of coordinates is suitable for the structure of the initial guess at hand, so for other problems or other forms of the initial guess, would require some permutation of the phase space coordinates to achieve a stable implementation.

Let us denote the flow map of a differential equation $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})$ with initial condition $\mathbf{x}(t_{0})=\mathbf{x}_{0}$ by $\phi(t;\mathbf{x}_{0})$. Thus, the displacement of a reference trajectory $\bar{\mathbf{x}}(t)$ after a time $\delta t$ becomes:

[TABLE]

Thus, a Taylor expansion of the displacement after $t_{1}+\delta t_{1}$ gives:

[TABLE]

where the first term on the right hand side is the state transition matrix, $\mathbf{\Phi}(t_{1},t_{0})$, evaluated along the reference trajectory initialized at $t_{0}$. The state transition matrix along a trajectory is obtained from the numerical solution of the variational equations along with the Hamilton’s equation of motion Parker and Chua (1989). Suppose that we want to land at a point $\mathbf{x}_{\rm d}$ (this would be the starting initial condition for a periodic orbit), after an integration time interval $t_{1}$, and starting from the initial guess $\bar{\mathbf{x}}_{\rm 0,g}$, then we have:

[TABLE]

where the error $\delta\bar{\mathbf{x}}_{1}$ is the applied first order correction obtained from the state transition matrix evaluated along the trajectory with initial condition $\bar{\mathbf{x}}_{\rm 0,g}$ and integrated for $t_{1}$ time units. For the rank-1 saddle equilibrium point under consideration, we initialize the guess as obtained in Eq. (74) with $\bar{\mathbf{x}}_{\rm 0,g}$. Using numerical integration, the initial condition is integrated until a $p=0$ event occurs with a high specified tolerance (typically $10^{-14}$). This results in $\bar{\mathbf{x}}(t_{1})$, which for the guess initial condition denotes the half-period location $t_{1}=T_{0,g}/2$. Then, we evaluate the state transition matrix $\mathbf{\Phi}(t_{1},0)$ for the trajectory obtained from the guess intial condition. Now, this is used to correct the initial value of $x_{0,g}$ while keeping $q_{0,g}$ constant and iterating until $p_{x}=0$. Since the $x$ configuration coordinate is kept constant, the first order correction is given by:

[TABLE]

where $\Phi_{ij}$ is the $(i,j)^{th}$ entry of $\mathbf{\Phi}(t_{1},0)$ and the acceleration terms come from the equations of motion evaluated at the crossing $t=t_{1}$ when $p_{1}=\delta p_{1}=0$. Thus, we obtain the first order correction $\delta x_{0}$ as:

[TABLE]

which is iterated until $|p_{x_{1}}|=|\delta p_{x_{1}}|<\epsilon$ for some tolerance $\epsilon$, since we want the final point of the periodic orbit to be of the form $\bar{\mathbf{x}}_{t_{1}}=(q_{1},x_{1},0,0)$. We remark that in all this process, differential correction assumes that the guess periodic orbit has a small error (for example in this system, of the order of $10^{-2}$) and can be corrected using first order form of the correction terms. If, however, large a corrective step is used, the half-orbit overshoots between successive steps leading to failure in converging to a closed orbit. Once the appropriate conditions are chosen, differential correction generates a family of periodic orbits and takes 2-3 iterations per unstable periodic orbit.

Step 4: Numerical continuation to the UPO at the selected excess energy. The procedure described above yields an initial condition for an unstable periodic orbit from an initial guess. Since our initial guess came from the linearization near the rank-1 saddle equilibrium pointgiven by (73), we can use this procedure for small amplitudes of order $2\times 10^{-5}$. We remark that this procedure is based on computations presented for this problem and will vary for a different nonlinear system. This small amplitude corresponds to small excess energy, typically of the order $10^{-5}$, and to obtain the unstable periodic orbit of arbitrarily large amplitude, we resort to numerical continuation for generating a family of periodic orbits that reach the selected excess energy.

This procedure starts with the initial conditions for two nearby unstable periodic orbits of small amplitude to obtain an initial guess for the next periodic orbit. The initial guess obtained by a simple extrapolation can then be corrected using differential correction. To this end, we proceed as follows. Suppose we find two nearby small amplitude unstable periodic orbits with initial conditions $\bar{\mathbf{x}}_{0}^{(1)}$ and $\bar{\mathbf{x}}_{0}^{(2)}$, correct to within the tolerance $d_{\rm tol}$ obtained using the differential correction procedure described above. We can then generate a family of periodic orbits with increasing amplitudes around $\bar{\mathbf{x}}_{\rm eq}$ as:

[TABLE]

A linear extrapolation to an initial guess of slightly larger amplitude, $\bar{\mathbf{x}}_{0}^{(3)}$ is given by:

[TABLE]

Thus, we can use differential correction on this guess initial condition to compute an accurate solution $\bar{\mathbf{x}}_{0}^{(3)}$ from the initial guess $\bar{\mathbf{x}}_{\rm 0,g}^{(3)}$ and repeat as an iterative step to generate a family of unstable periodic orbits.

Next, to compute an unstable periodic orbit at the selected excess energy, we track the energy of each unstable periodic orbit in the family until we have two solutions, $\bar{\mathbf{x}}_{0}^{\rm(k)}$ and $\bar{\mathbf{x}}_{0}^{\rm(k+1)}$, whose energy brackets the selected excess energy $\Delta H$. Then, we can resort to combining a bisection type method with differential correction on the two periodic orbits until we converge to the desired periodic orbit to within a specified tolerance. Thus, the result is an unstable periodic orbit at the selected energy $H_{0}$ and is specified by the initial condition $\bar{\mathbf{x}}_{\rm po}(0)$ and time period $T_{\rm po}$.

Step 5: Computation of invariant manifolds of the NHIM. We find the global approximation to the unstable and stable manifolds of the periodic orbit from the eigenvectors of the monodromy matrix. The local linear approximation of the unstable (or stable) manifolds (initial conditions displaced along the saddle space eigenvectors) is integrated using the full nonlinear equations of motion to produce the global approximation of the unstable (or stable) manifolds. This procedure is known as globalization of the manifolds and we proceed as follows.

First, the state transition matrix $\Phi(t)$ along the periodic orbit with initial condition $X_{0}$ can be obtained by numerical integration of the variational equations along with the equations of motion from $t=0$ to $t=T_{\rm po}$. This is gives the monodromy matrix $M=\Phi(T_{\rm po})$ and its eigenvalues are obtained. For Hamiltonian systemsMeyer et al. (2009), it is known that the eigenvalues of $M$ are of the form:

[TABLE]

The eigenvector associated with eigenvalue $\lambda_{1}$ is in the unstable direction, the eigenvector associated with eigenvalue $\lambda_{2}$ is in the stable direction. Let $e^{s}(\bar{\mathbf{x}}_{\rm po}(0))$ denote the normalized stable eigenvector, and $e^{u}(\bar{\mathbf{x}}_{\rm po}(0))$ denote the normalized unstable eigenvector. We can compute the invariant manifolds by initializing along these eigenvectors as:

[TABLE]

for the stable/unstable manifold of the desired periodic orbit. Here the small displacement from $\mathbf{X}_{0}$ is denoted by $\epsilon>0$ and its magnitude is taken to be small enough for the validity of the linearization, yet not so small that the time of flight becomes too large due to asymptotic nature of the stable and unstable manifolds. It has been suggestedKoon et al. (2011) that typical values of $\epsilon$ around $10^{-9}$ correspond to nondimensional position displacements of magnitude around $10^{-6}$.

By numerically integrating the unstable vector forwards in time, using both $\epsilon$ and $-\epsilon$, for the forward and backward branches respectively, we generate trajectories shadowing the two branches, $\mathcal{W}^{u}_{+}$ and $\mathcal{W}^{u}_{-}$. Similarly, by integrating the stable vector backwards in time, using both $\epsilon$ and $-\epsilon$, for forward and backward branches respectively, we generate trajectories shadowing the stable manifold, $\mathcal{W}^{s}_{+,-}$. For the manifold at $\mathbf{X}(t)$, one can simply use the state transition matrix to transport the eigenvectors from $\mathbf{X}_{0}$ to $\mathbf{X}(t)$

[TABLE]

It is to be noted that since the state transition matrix does not preserve the norm, and hence the resulting vector must be normalized. The globalized invariant manifolds associated with rank-1 saddles are known as Conley-McGehee tubesMarsden and Ross (2006). These tubes form the impenetrable phase space conduits (codimension-1 barriers) for the trajectories to react by crossing the $q=0$ dividing surface.

In summary, the computation of invariant manifolds of the unstable periodic orbit associated with the rank-1 saddle begins with the linearized equations of motion. This is obtained after a coordinate transformation to the rank-1 saddle equilibrium point and a Taylor expansion of the equations of motion. Keeping the first order terms in this expansion, we obtain the eigenvalues and eigenvectors of the linearized system. The eigenvectors corresponding to the center subspace provide the starting guess for computing the unstable periodic orbits of small excess energy, $\Delta H<<1$. This iterative procedure performs a small correction to the starting guess for initial condition based on the terminal condition of the periodic orbit until a desired tolerance is satisfied. This procedure is known as differential correction and generates initial condition for an unstable periodic orbit at small excess energy. Next, a numerical continuation procedure is adopted to follow the small energy (small amplitude) periodic orbit out to high excess energy. Once the unstable periodic orbit is obtained, the globalization of its invariant manifolds is done using the initial conditions along the eigenvectors of the state transition matrix computed along the unstable periodic orbit. The result of the steps described above is shown in Fig. 9.