Geometrodynamical Analysis to Characterize the Dynamics and Stability of a Molecular System through the Boundary of the Hill ’ s Region

In this paper we study the dynamics and stability of a two-dimensional model for the vibrations of the LiCN molecule making use of the Riemannian geometry via the Jacobi-Levi-Civita equations applied to the Jacobi metric. The Stability Geometrical Indicator for short times is calculated to locate regular and chaotic trajectories as the relative extrema of this indicator. Only trajectories with initial conditions at the boundary of the Hill’s region are considered to characterize the dynamics of the system. The importance of the curvature of this boundary for the stability of trajectories bouncing on it is also discussed.


Introduction
The phase space structure for the vibrational dynamics of molecular systems with strong chaotic behavior is a very interesting topic in chemical processes [1].
In this paper we study the nonlinear Hamiltonian dynamics and stability of the vibrational dynamics of the isomerizing LiNC/LiCN molecular system [9]- [15] from a geometrodynamical perspective using the SGI indicator, by rewriting the dynamics in a geometric language, where a Riemannian manifold is defined by adding a particular metric to the corresponding configuration space.To this end different metrics are used in literature.Among them, the Eisenhart [16], Finsler [17], Horwitz [18]- [21] or Jacobi [22]- [29] are worth mentioning.We consider the Jacobi metric obtained from the Maupertuis' Principle (1744), where geodesics [30] [31] (natural motions) are identified with trajectories of the system, and the stability can be analyzed through the study of the curvature of this Riemannian manifold via the Jacobi-Levi-Civita equations (JLC) [31].Once the dynamical system is geometrized, physical magnitudes can be identified with geometrical ones, as indicated in Table 1.
This approach has been applied successfully in the calculation of Lyapunov exponents in Hamiltonian systems with many degrees of freedom [25], after simplifying the JLC equations by coarse graining.We can define via the JLC equations the SGI indicator along trajectories.This indicator has been recently used to study the stability of trajectories in this molecular system [8].We will use SGI for short times to locate stable periodic trajectories (with its corresponding chains of islands) and chaotic ones as local minima or maxima respectively of this indicator.
The boundary of the Hill's region 1 for systems with Gaussian curvature everywhere positive (locally stable) has been considered as a very important set to generate unstabilities because of its curvature [32].As in the classical billiard problem non convex points (negative curvature) of the boundary are considered defocusing points related with unstable behavior of trajectories reaching these points.Thus trajectories probing the non convex points of the boundary are always chaotic and those trajectories reaching uniquely convex points of the boundary are regular as indicated in [32].In this paper we show that this statement is not true in general showing that although our system is positively curved everywhere we can locate chaotic trajectories restricted to convex parts of the boundary of the Hill's region and regular trajectories reaching non convex points.
The organization of the paper is as follows.In Sections 2 and 3 the geometric approach is described, relating chaos with JLC equations.In Section 4 the model LiNC/LiCN and its geometrization are presented.In Section 5 we present the results for trajectories with initial conditions at the boundary of the Hill's region using the SGI.In Section 6 the importance of the boundary in order to study the stability of trajectories is discussed.Finally, the conclusions are summarized in Section 7.

Geometrodynamical Perspective
We consider the conservative system described by the Lagrangian 2 ( ) ( ) ( ) Table 1.Equivalence between some relevant magnitudes in the Dynamical and Geometrical views of the Classical Mechanics.

Conserved quantities Symmetries, killing fields
Instability Curvature 1 Projection of the accessible region in the configuration space on the position coordinates. 2We will use the Einstein's summation convention for repeated indexes, so a q the kinetic energy tensor and d d , defined in the domain of the configuration space ( ) , i j q q  .We can make the following identification p the conjugate momenta in the phase space ( ) , i j q p .The Hamilton's principle of least action can be connected with the Maupertuis' principle which establishes as the natural motions those whose asynchronous trajectories with fixed end points are extrema for the action integral  , i.e.
where W is the kinetic energy and E the total energy (conserved) being ds the differential of length 3 .So the trajectories for the system corresponds to geodesics for the Riemannian manifold endowed with the metric that in coordinates can be written as called Jacobi metric.The configuration space for the system is the previous Riemannian manifold ( ) ( ) , , , n q q q = q  represent the local coordinates and being referred to as mechanical manifold.It can be shown that the trajectory followed by a particle affected by the potential energy function ( ) V q in the Euclidean space is equivalent to the motion of a free particle in the restricted mechanical manifold . Thus there is a correspondence between dynamics in an Euclidean space (zero curvature) and kinematics in a Riemannian manifold.Geodesics are autoparallel curves, so they are described through these equations with ∇ the affine connection compatible with the metric, in coordinates and i jk Γ the Christoffel symbols from the metric ij g , which are calculated as ( ) Changing ds by dt ( ) ( ) arc-length parameter.s ≡ that corresponds to trajectories parametrized by physical time t .

Curvature and Chaos. Jacobi-Levi-Civita Equations (JLC)
The sectional curvature ( ) at the point p M ∈ and tangent plane σ on the Riemannian manifold M at this point (Gaussian curvature in dimension two), determines the velocity of spreading for the geodesics starting at p tangent to σ .
The evolution for this geodesic deviation obeys the JLC equations the covariant derivative along the curve ( ) the Riemann-Christoffel curvature tensor [30].
In our case, ( ) the geodesic deviation vector field along the fiducial trajectory ( ) with i k S the stability tensor, in a local chart ( ) The matrix ( ) with , the inner product, that can be diagonalized by solving 0, with eigenvalues 1 0 λ = and 2 1212 R g λ = , and the corresponding eigenvectors ( )    being mutually orthogonal.These eigenvectors when normalized define a convenient basis , , = e e u u u u obtaining the new convenient moving frame in terms of the canonical local chart of M as follows q q e e e q q     (1.10) The Equations (1.7) can be easily expressed in this basis being parallel transported along ( ) .11) with ( ) ( ) = e e e we solve (1.11) using the Taylor expansion (1.12) If the vector field ( ) we can write (1.9) ( ) with R the scalar curvature for metric ij g .The evolution for 1 J is as much lineal and the stability is determined by the scalar 2 J (from now on referred as J ) along the fiducial curve.Fixing a constant length s for the trajectory, we define the dynamical magnitude Stability Geometric Indicator (SGI) as ( ) ( ) with initial conditions 0 x .We will use this indicator to compare the mutual relative stabilities for different trajectories.

The Model
We study the two degrees of freedom model for the isomerizing molecular system LiNC/LiCN, where the bond NC is kept frozen at its equilibrium distance 2.186 e r = a.u.Studying the problem in the phase space, the classical vibrational (with zero total angular momentum) Hamiltonian for the dynamical system in Jacobi coordinates is given by where R is the distance from the Li atom to the NC center of mass, e r is the NC distance, θ is the angle de- fined by the R and CN vectors, R p and p θ are the corresponding conjugate momenta, and The potential energy surface ( ) V R θ is defined by an expansion in Legendre polynomials ( ) where parameters are taken from the literature [33].As can be seen in Figure 1 it presents two potential wells (elliptic points), corresponding each one to a stable isomer with linear configuration: LiCN where

Jacobi Metric
Considering the Hamiltonian of the system (1.15) then 4 ( ) and the Riemannian manifold ( ) is the inverse matrix of ( ) jk a .
( ) ( ) The Christoffel symbols for the metric (1.17) are calculated using (1.5), obtaining Trajectories parametrized by physical time t are calculated integrating the Equations (1.6) after replacing the corresponding values of i jk

Jacobi-Levi-Civita (JLC) Equations
Geodesic deviation vector field ( ) along the fiducial trajectory ( ) s γ is used in the JLC Equation (1.7).Rewriting Equation (1.13) for physical time t we obtain where the scalar curvature R for the metric ( ) g R θ determines the local stability for geodesics around the fiducial one at ( ) , R θ .If N is the dimension of the system it can be written as ( ) with R  the scalar curvature for the space and 2 ∇ the Laplace-Beltrami operator both for metric ( ) . 6 .
The expression for the curvature (1.22) of the mechanical manifold ( ) This scalar curvature (Gaussian curvature) is everywhere positive, i.e. every trajectory is locally stable.The differential equations system to integrate numerically consists on (1.19) and (1.20) simultaneously.

Locating Stable and Chaotic Trajectories with SGI at the Boundary of the Hill's Region
In order to characterize the dynamics of the system the selection of initial conditions to integrate trajectories is very important.All the trajectories we have integrated reach the boundary of the Hill's region at some time.This fact motivate to study the dynamics in our molecular system through trajectories with initial conditions at this boundary.
We have calculated the value of SGI for short times of trajectories with initial conditions ( ) , R θ all along the boundary ( ) R p = and 0 p θ = (kinetic energy equals zero) with a value for energy  1) to (8).
The stable periodic trajectories correspond with local minima values of SGI in particular the points (1), ( 2), (3), ( 5) and (7), and the chaotic trajectories with ( 4) and (8).The trajectory represented with (7) is the most stable trajectory of the system related with the KAM torus with the most irrational relation of frequencies, being this torus the last one broken in the transition from order to chaos.The corresponding trajectories in ( ) , R θ are represented in the upper graphs of Figure 3.
Chain of islands are represented as those forming the corresponding valleys in the global SGI diagram, the stability of the resonant periodic orbit in the minimun is indicated by the width of the valley as illustrated in the Poincaré Sections in Figure 4.
In order to locate a chaotic trajectory from the global SGI diagram we search for hyperbolic points located between the chains of islands (Poincaré Birkchoff Theorem) so we take initial conditions at the extrema of the previous valleys represented by local maxima labeled by ( 4) and (8) in Figure 3.The Poincaré sections for all these trajectories are illustrated in Figure 4.
The indicator SGI is calculated along trajectories with initial conditions at the boundary for energy   .Stable periodic trajectories are marked as (1), ( 2), (3), ( 5), ( 7) and chaotic trajectories as ( 4) and ( 8).The quasiperiodic trajectory ( 6) is the most stable trajectory of the system.Corresponding trajectories in ( ) , R θ are plotted in the upper graphs.

SGI and the Curvature of the Boundary of the Hill's Region
The boundary of the Hill's region and particularly the sign of its curvature is often related with the stability of trajectories bouncing on it.In systems with Gaussian curvature everywhere positive inside the Hill's region is quite natural consider the non convex (negative curvature) regions of the boundary as source of unstabilities for trajectories bouncing on it.In this sense we could say that regular trajectories should probe or touch only convex regions (positive curvature) of the boundary of the Hill's region as well as chaotic trajectories should probe the non convex region at least once [32].
We show that this statement is not true in general finding a stable periodic trajectory bouncing on non convex  The curvature κ of the boundary of the Hill's region considered as the set level ( ) ( ) We study the convexity for different set levels   The behavior of the chaotic trajectory comes to support the well known idea of the parametric resonance as fundamental reason for the onset of chaos.But the regular trajectory motivates the idea of an additional stabilization mechanism.
When defining the Jacobi metric ( ) , ij g R θ is quite often to understand the potential energy as curving the underlying space (with the kinetic energy tensor ( ) , ij a R θ as metric), in the particular case of an underlying space of curvature zero (usually ij ij a δ = ) the first term in the right hand side of Equation (1.22) is zero and the curvature is caused by the energy potential.In our molecular system the kinetic energy tensor ( ) , ij a R θ defines an underlying space positively curved with curvature R  expressed in (1.26).Thus although the potential energy is removed the space is equally positively curved.
kSis symmetric and defines a selfadjoint mapping A γ relative to ij g as : -NC and NC reduced masses, respectively.
minimun energy path (MEP) connecting the two isomer wells has been plotted as a dashed line superimposed to the potential energy surface in Figure1.It can be written as the following Fourier expansion:

Figure 1 .
Figure 1.(Left) Potential energy surface for the LiNC/LiCN isomerizing system.Contours lines have been plotted every 1 1000 cm − .The minimum energy path connecting the two isomers is shown as a dashed line.(Right) Energy profile along the minimum energy path.

+ 1 .
The energy profile associated to the MEP is given in the right panel of Figure The dynamics is monitored by PSOS (Poincaré surface of section).In order to obtain the most complete dynamical information we compute them taking the minimun energy path, ( ) e R θ , as the sectioning plane[34].This PSOS is made an area preserving map defining a new set of canonical variables ( )Nowthe PSOS is defined by the condition 0 ρ = and P ρ is obtained from the Hamiltonian conservation ( ) negative branch in the calculation of the second-order equation obtained.We represent PSOS values in the interval ( ) 0, π ψ ∈ rad taking into account the symmetry of the molecular system.

obtaining the results plotted in Figure 2 .
1.22) is singular at these points so the calculations have been done for an energy E ˆ lightly higher with value Stable trajectories and chaotic ones for the interval ( )1.7, 2.1 θ ∈ rad for the upper region of the boundary (red colour in upper left graph in Figure 2) are represented in Figure 3 labeled from ( trajectory and a periodic one represented by the maximum and minimum points in the graph for SGI in Figure 5.The Poincaré sections are illustrated in Figure 6.

Figure 2 .
Figure 2. Values of SGI for different trajectories with initial conditions all along the boundary of the Hill's region (upper left graph) for the energy, 1 1500.5 cm E − = , the interval

Figure 4 .
Figure 4. Poincaré sections for all trajectories represented in Figure 2 and Figure 3 for the energy 1 1500.5 cm E − = .The black curve represents the boundary of the accessible region in the phase space.

1 1500
points of curvature zero (change of convexity) as the red curve in Figure7, the non convex points on the lower region of the boundary for the energy

Figure 6 .
Figure 6.Poincaré sections for the two trajectories represented in Figure 5 at the energy 1 500 cm E − = .The black curve represents the boundary of the accessible region in the phase space.

Figure 7 . 1 1500 4 .
Figure 7.A chaotic trajectory (a) restricted to convex points of the boundary of the Hill's region and a stable periodic trajectory (b) reaching the boundary at non convex points at the energy 1 1500.5 cm E − = .These trajectories correspond to (4) and (i) respectively in Figure 4.