# Elementary exposition of realizing phase space structures relevant to chemical reaction dynamics

In this article, we review the analytical and numerical approaches for computing the phase space structures in two degrees-of-freedom Hamiltonian systems that arise in chemical reactions. In particular, these phase space structures are the unstable periodic orbit associated with an index-1 saddle, the periodic orbit dividing surface, and the stable and unstable invariant manifolds of the unstable periodic orbit. We review the approaches in the context of a two degrees-of-freedom Hamiltonian with a quartic potential coupled with a quadratic potential. We derive the analytical form of the phase space structures for the integrable case and visualize their geometry on the three dimensional energy surface. We then investigate the bifurcation of the dividing surface due to the changes in the parameters of the potential energy. We also review the numerical method of turning point and present its new modification called the turning point based on configuration difference for computing the unstable periodic orbit in two degrees-of-freedom systems. These methods are implemented in the open-source python package, UPOsHam <cit.>.

Comments

There are no comments yet.

## Authors

• 1 publication
• 4 publications
• 3 publications
07/18/2021

### Support vector machines for learning reactive islands

We develop a machine learning framework that can be applied to data sets...
03/08/2021

### Self-learning Machines based on Hamiltonian Echo Backpropagation

A physical self-learning machine can be defined as a nonlinear dynamical...
02/25/2022

### Connecting Gaits in Energetically Conservative Legged Systems

In this work, we present a nonlinear dynamics perspective on generating ...
09/23/2019

### Computing the degrees of freedom of rank-regularized estimators and cousins

Estimating a low rank matrix from its linear measurements is a problem o...
10/01/2018

### Raster Grid Pathology and Other Ptychographic Ambiguities

Ptychography with an unknown probe and object is analyzed for the raster...
02/04/2020

### On the generation of periodic discrete structures with identical two-point correlation

Strategies for the generation of periodic discrete structures with ident...
10/01/2018

### Raster Grid Pathology and Other Ambiguities in Blind Ptychography

Blind ptychography is a phase retrieval method using multiple diffractio...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Brief background.  Hamiltonian systems [Meyer2009, wiggins2003applied] appear in a wide range of areas in natural sciences and engineering, for example, chemical reactions [NHIMinchem], celestial mechanics [koon_heteroclinic_2000, Jaffe2002, Dellnitz2005, Waalkens2005Escape], ship dynamics [naik_geometry_2017], and fluid dynamics [Mendoza2010]

. Dynamical systems is the study of continuous (or discrete) time evolution of initial conditions governed by ordinary differential equations (ODEs), partial differential equations, or function iterations, and a Hamiltonian system is an energy conservative system. For a

degrees-of-freedom system, the Hamiltonian dynamical system is described by a set of position coordinates and momenta which gives a dimensional space called the phase space where the initial conditions evolve in time. The time evolution of the initial conditions is governed by an autonomous or non-autonomous ordinary differential equation with respect to time called the Hamilton’s equations. Each initial condition is represented by a point in the phase space and the initial condition’s past and future location in the phase space is obtained by evolving the initial condition in backward and forward time, respectively. One of the phase space structures is the equilibrium point of a Hamiltonian dynamical system, which is also called as the critical point (where the gradient vanishes) of the potential energy surface, which exists at only one value of the total energy [Wiggins2017book]. In two DOF systems, equilibrium points can be of different types and in this review, we will focus on the center saddle (called an index-1 saddle) which corresponds to a top of the potential energy barrier and the center center which corresponds to a potential energy well. The type of the equilibrium point is deduced by linear stability analysis of the equilibrium point. One important point to consider is that the equilibrium points have zero momentum

(this follows from the definition) and exist at one energy but their influence at non-zero momentum which corresponds to higher energies is affected by the periodic orbit associated with the equilibrium point. The periodic orbit associated with an index-1 saddle equilibrium point is classified as unstable and its generalization to high-dimensional (

) phase space is called a Normally Hyperbolic Invariant Manifold (NHIM) [Wiggins1990, Wiggins1994]. The unstable periodic orbit (and its general form, NHIM) has been shown to be the appropriate phase space structure to construct a dividing surface that partitions the energy surface into the reactant and product regions and which is free of local recrossings [waalkens2004direct, Waalkens2005]. A phase space dividing surface leads to more accurate calculation of the reaction flux and rate constant in a chemical reaction [Ezra2009]. In this review, we illustrate the geometry of a phase space dividing surface by varying parameters of a potential energy surface. Thus, a dividing surface is the phase space structure (it is a dynamical object since momenta is taken into account) used in reaction rate calculations. In a sense this provides the local pictures of crossing the energetic barrier located by the saddle equilibrium point. However, another interest in studying chemical reactions is to test arbitrary configurations of a molecule for reactivity (potential of reaction in an environment). In this case, the global dynamics of the configurations needs to be understood and the appropriate phase space structures are the stable and unstable invariant manifolds of the unstable periodic orbit (NHIM in -DOF systems) associated with the saddle equilibrium point. The stable and the unstable manifolds are codimension-1 in the energy surface and are impenetrable barriers [wiggins_impenetrable_2001] that partition trajectories with different dynamical fates.

Hamiltonian mechanics in the context of chemical reactions.  In the context of chemistry, the configuration space is described by the bond lengths between atoms in a molecule and denoted by the position coordinates, , in a degrees-of-freedom system. In this setup, one has to define what reaction means from a mathematical point of view. For example, a dissociation reaction involving occurs when one or more bonds break which translates to one or more position coordinates becoming unbounded when they evolve in time. These position coordinates are called the reaction coordinates. The position coordinates are also involved in describing the potential energy in a chemical reaction. Thus, the configuration space perspective is concerned with the features and shape of a potential energy function. However, the time evolution of the position coordinates, , depends on the momenta, , due to the kinetic energy of the bonds which arises due to the molecular motion. The combination of the potential energy and kinetic energy defines the Hamiltonian while the space in which the coordinates and momenta evolve is called the phase space [agaoglou_chemical_2019] of the Hamiltonian system. Thus, the configuration space and phase space of a molecular system are inherently related and the appropriate setting for understanding reactions is the phase space.

In this review, we illustrate the geometry of the aforementioned phase space structures by using a two DOF model Hamiltonian where all the relevant structures can be calculated analytically. It is to be noted that linear Hamiltonian systems are considered as the first step [NHIMinchem, Geomodels] in these theory since only one equilibrium point, either a saddle center or center center, is possible. The natural next step is to consider a model nonlinear Hamiltonian by adding higher order terms in the potential energy function and thus allows for multiple equilibrium points. We consider a Hamiltonian (Eqn. (1)) with a quartic term in the potential energy function, , and use this Hamiltonian to discuss the associated phase space structures along with the visualizations. In particular, we obtain analytical expressions for the NHIM which is an unstable periodic orbit associated with the index-1 saddle equilibrium point which makes the visualization feasible in the three-dimensional energy surface of the two DOF Hamiltonian (Eqn. (1)). We demonstrate the bifurcation of the phase space DS as the parameters in the potential energy are varied and which can not happen in linear Hamiltonian systems. Next, we demonstrate numerical methods that need to be adopted when analytical expressions of the UPOs can not be obtained which is what happens in most realistic models of chemical reactions. Some of these numerical methods include differential correction and numerical continuation [naik_finding_2019b, Koon2011], multiple shooting [farantos_pomult_1998], turning point [Pollak1980], and turning point based on configuration difference [Lyu2020]. We describe the turning point methods for the model Hamiltonian in this article.

The article is outlined as follows. First, we briefly describe the two degrees-of-freedom (DOF) model Hamiltonian in section 2 along with the basic dynamical systems analysis. Then, we obtain the phase space structures for different combinations of the signs of the parameters in the potential energy in section 3.1. We discuss the resulting bifurcations in section 3.2 and the phase space structures are visualized in section 3.3. We present the results obtained using the turning point methods in section 3.4. Then, we summarize and present our outlook in section 4.

## 2 Two DOF model: Uncoupled quartic Hamiltonian

The Hamiltonian structure is given by

 H=p212−αq212+βq414H1+ω2(p22+q22)H2, α, β are some parameters. (1)

The Hamilton’s equations of are given by:

 ˙q1 =∂H∂p1=p1, (2) ˙p1 =−∂H∂q1=αq1−βq31, (3) ˙q2 =∂H∂p2=ωp2, (4) ˙p2 =−∂H∂q2=−ωq2. (5)

Note that () stands for derivative with respect to time .

Now we want to classify all equilibrium points and their stability associated to . If we let

 ˙q1 =∂H∂p1=p1=0, ˙p1 =−∂H∂q1=αq1−βq31=0, ˙q2 =∂H∂p2=ωp2=0, ˙p2 =−∂H∂q2=−ωq2=0.

we can see that equilibrium points are

 (¯q1,¯p1,¯q2,¯p2)=(0,0,0,0),(±√α/β,0,0,0) (6)

Note that is an equilibrium point for all and are equilibrium points if and only if . However we need to consider the case separately as the fraction does not make sense when

. The eigenvalues of the Jacobian associated to

are and the eigenvalues of the Jacobian J associated to are the solutions of the characteristic equation:

 det(J−λI)=det(−λ1α−3βq21−λ)=λ2−(α−3βq21)=0 (7)

where

 J=(01α−3βq210) (8)

Rearranging we have .

## 3 Results and Discussions

### 3.1 Phase space structures

##### 1. degenerate case α=β=0.

First, let’s consider the trivial case when . Our becomes and (2), (3) become

 ˙q1 =∂H∂p1=p1, ˙p1 =−∂H∂q1=0.

which imply that is a constant, i.e. phase diagram is consisted of horizontal lines in plane.

##### 2. degenerate case α≠0,β=0.

In this case, our Hamiltonian becomes

 H1=p212−αq212

we only have one equilibrium point and the eigenvalues of the Jacobian are . It’s exactly the same to what we have discussed before.

1. If , the eigenvalues of the Jacobian are purely imaginary, i.e. the origin is a centre.

2. If , the eigenvalues of the Jacobian are real and of opposite signs, i.e. the origin is a saddle. Therefore is a "saddle-centre" equilibrium point. The full energy surface is given by

 H=p212−αq212+ω2(p22+q22)=H1+H2>0,H1>0,H2≥0

The DS can be taken as and it is nonisoenergetic since it’s not consisted of trajectories in the phase space. The DS has two hemispheres which correspond to forward and backward reactions.

 (q1,p1,q2,p2)∈R4:p1=+√2√H1+H2−ω2(p22+q22)>0, forward DS,
 (q1,p1,q2,p2)∈R4:p1=−√2√H1+H2−ω2(p22+q22)<0, backward DS.

The forward and backward DS meet at :

 S1NHIM(h)=Σ(h)∩{q1=0,p1=0}={(q1,p1,q2,p2)∈R4:ω(p22+q22)=2h}, q1=0, p1=0.

The NHIM is given by

 ω2(p22+q22)=H1+H2≥0

the stable and unstable manifolds are given by

 Wu(NHIM) ={(q1,p1,q2,p2)|p1=√αq1, ω2(p22+q22)=H2}, Ws(NHIM) ={(q1,p1,q2,p2)|p1=−√αq1, ω2(p22+q22)=H2}.
##### 3. degenerate case α=0,β≠0.

If , the Hamiltonian becomes

 H1=p212+βq414

we only have one equilibrium point, which is the origin. The eigenvalues of the Jacobian are i.e. the equilibrium point is non hyperbolic. Therefore is nonhyperbolic. We can still draw the phase diagrams in plane depending on the signs of .

1. If , we can plot the potential energy function which is a "" shape and the phase diagram is similar to the one which has centre stability.

2. If , similarly we can plot the potential energy function which is a "" shape and the phase diagram is similar to the one which has saddle stability. The invariant manifolds associated to the origin are given by

 p21=−β2q41

We can write explicitly as

 Wc((0,0)) ={(q1,p1)|p1=√−β2q21}, Wc((0,0)) ={(q1,p1)|p1=−√−β2q21}.
##### 4. Nonzero α,β.

Now we consider . The Hamiltonian is

 H1=p212−αq212+βq414 (9)

Recall that the eigenvalues of the Jacobian are .

1. If , we only have 1 equilibrium point as in (6) is negative. are purely imaginary, i.e. the origin is a centre.

2. If , we only have 1 equilibrium point for the same reason. are real and of opposite signs, i.e. the origin is a saddle. Therefore is a "saddle-centre" equilibrium point. The full energy surface is given by

 H=p212−αq212+βq414+ω2(p22+q22)=H1+H2>0,H1>0,H2≥0. (10)

The DS can be taken as . The DS is nonisoenergetic since it’s not consisted of trajectories in the phase space. The DS has two hemispheres which correspond to forward and backward reactions.

 (q1,p1,q2,p2)∈R4:p1=+√2√H1+H2−ω2(p22+q22)>0, forward DS, (11)
 (q1,p1,q2,p2)∈R4:p1=−√2√H1+H2−ω2(p22+q22)<0. backward DS. (12)

The forward and backward DS meet at :

 S1NHIM(h)=Σ(h)∩{q1=0,p1=0}={(q1,p1,q2,p2)∈R4:ω(p22+q22)=2h}, q1=0, p1=0. (13)

The NHIM is given by

 ω2(p22+q22)=H1+H2≥0 (14)

the stable and unstable manifolds are given by

 p21=αq21−βq412 (15)

We can write explicitly as

 Wsf(NHIM) ={(q1,p1,q2,p2)|p1=√αq21−βq412, ω2(p22+q22)=H2}, q1<0, p1>0, (16) Wsb(NHIM) ={(q1,p1,q2,p2)|p1=−√αq21−βq412, ω2(p22+q22)=H2}, q1>0, p1<0, (17) Wuf(NHIM) ={(q1,p1,q2,p2)|p1=√αq21−βq412, ω2(p22+q22)=H2}, q1>0, p1>0, (18) Wub(NHIM) ={(q1,p1,q2,p2)|p1=−√αq21−βq412, ω2(p22+q22)=H2}, q1<0, p1<0. (19)
3. If , we have 3 equilibrium points . The stability of these equilibrium points are determined by . When , we have , i.e. the origin is a saddle. when , we have , i.e. are centres. Note that the equilibrium point exists on the energy surface and exist on the energy surface . The phase space structure of is given in Fig. 1. Note also that is a "saddle-centre" equilibrium point.

In order for reaction to occur, we must have . Therefore the full energy surface is given by

 H=p212−αq212+βq414+ω2(p22+q22)=H1+H2>0, H1>0,H2≥0. (20)

The DS can be taken as . The DS is nonisoenergetic since it’s not consisted of trajectories in the phase space. The DS has two hemispheres which correspond to forward and backward reactions.

 (q1,p1,q2,p2)∈R4:p1=+√2√H1+H2−ω2(p22+q22)>0, forward DS, (21)
 (q1,p1,q2,p2)∈R4:p1=−√2√H1+H2−ω2(p22+q22)<0, backward DS. (22)

The forward and backward DS meet at :

 S1NHIM(h)=Σ(h)∩{q1=0,p1=0}={(q1,p1,q2,p2)∈R4:ω(p22+q22)=2h}, q1=0, p1=0. (23)

The NHIM is given by

 ω2(p22+q22)=H1+H2≥0 (24)

the stable and unstable manifolds are given by

 Wsf(NHIM) ={(q1,p1,q2,p2)|p1=√αq21−βq412, ω2(p22+q22)=H2}, q1<0, p1>0, (25) Wsb(NHIM) ={(q1,p1,q2,p2)|p1=−√αq21−βq412, ω2(p22+q22)=H2}, q1>0, p1<0, (26) Wuf(NHIM) ={(q1,p1,q2,p2)|p1=√αq21−βq412, ω2(p22+q22)=H2}, q1>0, p1>0, (27) Wub(NHIM) ={(q1,p1,q2,p2)|p1=−√αq21−βq412, ω2(p22+q22)=H2}, q1<0, p1<0. (28)
4. If , we have 3 equilibrium points . The stability of these equilibrium points are given by . When , we have , i.e. the origin is a centre. when , we have , i.e. are saddles. Note that the equilibrium point exists on the energy surface and exist on the energy surface , i.e. the saddles exist on the energy surface which has more energy than the origin. The phase space structure of is given in Fig. 2. Note also that are "saddle-centre" equilibrium points.

In order for reaction to occur, we must have . Therefore the full energy surface is given by

 H=p212−αq212+βq414+ω2(p22+q22)=H1+H2>−α24β, H1>−α24β>0,H2≥0. (29)

The DS can be taken as or . The DS is nonisoenergetic since it’s not consisted of trajectories in the phase space. Each DS has two hemispheres which correspond to forward and backward reactions.

 (q1,p1,q2,p2)∈R4:p1=+√2√H1+H2−ω2(p22+q22)+α24β>0, q1=√α/β, (30) forward DS(q1=√α/β), (31)
 (q1,p1,q2,p2)∈R4:p1=−√2√H1+H2−ω2(p22+q22)+α24β<0, q1=√α/β, (32) backward DS(q1=√α/β), (33)
 (q1,p1,q2,p2)∈R4:p1=+√2√H1+H2−ω2(p22+q22)+α24β>0, q1=−√α/β, (34) forward DS(q1=−√α/β), (35)
 (q1,p1,q2,p2)∈R4:p1=−√2√H1+H2−ω2(p22+q22)+α24β<0, q1=−√α/β, (36) backward DS(q1=−√α/β). (37)

The forward and backward DS meet at or :

 S1NHIM(h)=Σ(h)∩{q1=±√α/β,p1=0}={(q1,p1,q2,p2)∈R4:ω(p22+q22)=2h+α22β}, (38) q1=±√α/β,p1=0. (39)

The NHIM is given by

 ω2(p22+q22)=H1+H2≥−α24β (40)

the invariant manifolds of the NHIMs are given by

 −α24β=p212−αq212+βq414 (41)

the stable and unstable manifolds are given by

 Ws(NHIM with q1=√α/β,q2=0) ={(q1,p1,q2,p2)|p1=√αq21−βq412−α22β, ω2(p22+q22)=H2}, q1<0, (42) Wu(NHIM with q1=√α/β,q2=0) ={(q1,p1,q2,p2)|p1=−√αq21−βq412−α22β, ω2(p22+q22)=H2}, q1<0, (43) Ws(NHIM with q1=−√α/β,q2=0) ={(q1,p1,q2,p2)|p1=−√αq21−βq412−α22β, ω2(p22+q22)=H2}, q1>0, (44) Wu(NHIM with q1=−√α/β,q2=0) ={(q1,p1,q2,p2)|p1=√αq21−βq412−α22β, ω2(p22+q22)=H2}, q1>0. (45)

### 3.2 Bifurcation analysis

We have discussed the number and stability of the equilibrium points for different signs of and therefore we can plot the Bifurcation diagram of the Hamiltonian in plane in Fig. 3.

When , , we have 1 centre equilibrium point, the origin. If we fix and decrease the value of , we still have 1 equilibrium point until . As becomes negative, the origin remains a centre and we have 2 new saddle equilibrium points at . In this case, our parameters and . As the number of equilibrium points change from 1 to 3, we have Hamiltonian Pitchfork Bifurcation [Wiggins2017book].

Now we fix and increase the value , we have 3 equilibrium points until . As evolves from negative to positive, disappear and the origin is non hyperbolic when and becomes a hyperbolic saddle point when . When , , we only have 1 equilibrium point and therefore we have Hamiltonian Pitchfork Bifurcation.

If we fix and increase from negative to positive, we have 3 equilibrium points. The origin remains a saddle and two centres occur at . In this case, , we still have Hamiltonian Pitchfork Bifurcation.

Finally, if we fix and let decrease from positive to negative. We have 3 equilibrium points until . The origin loses its hyperbolicity at . When , , we have 1 centre equilibrium point again. Hamiltonian Pitchfork Bifurcation occurs for the same reason as above.

### 3.3 Visualisation of the phase space structures

In this section we use Matlab to simulate various phase structures associated to the Hamiltonian discussed before. We are only interested in the case and and details of other cases can be found in  [Geomodels, NHIMinchem]. The method we use to visualise the phase space structures is model I in  [Geomodels] where is the energy surface with total energy . The sign of is postive on and the sign of is negative on . Given an initial condition , the trajectory evolves on an invariant manifold:

 ΛI,J(h)={(q1,p1,q2,p2)∈R4:H1=(p01)22−α(q01)