Applied Mathematics
Vol.5 No.16(2014), Article ID:49878,13 pages DOI:10.4236/am.2014.516251

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

Alberto Vergel1, Rosa M. Benito1, Juan C. Losada1, Florentino Borondo2,3

1Grupo de Sistemas Complejos, Departamento de Física y Mecánica, Escuela Técnica Superior de Ingenieros Agrónomos, Universidad Politécnica de Madrid, Madrid, Spain

2Instituto de Ciencias Matemáticas (ICMAT), Madrid, Spain

3Departamento de Química, Universidad Autónoma de Madrid, Madrid, Spain


Copyright © 2014 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

Received 24 July 2014; revised 26 August 2014; accepted 13 September 2014


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.

Keywords:Geometrodynamical Analysis, Molecular Nonlinear Dynamics, Stability Geometrical Indicator, Riemannian Geometry, Hamiltonian Systems

1. 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] .

For practical applications several chaos indicators are used in literature, such as maximun Lyapunov exponent (LE) [2] , fast Lyapunov indicator [3] [4] , frequency map analysis [5] , small alignment index (SALI) [6] [7] , or the stability geometrical indicator (SGI) [8] .

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 Table1

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 region1 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.

2. Geometrodynamical Perspective

We consider the conservative system described by the Lagrangian2



Table 1. Equivalence between some relevant magnitudes in the Dynamical and Geometrical views of the Classical Mechanics.  

with the kinetic energy tensor and, defined in the domain of the configuration space

. We can make the following identification

being the conjugate momenta in the phase space.

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. with


where is the kinetic energy and the total energy (conserved)

being the differential of length3. 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

where 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 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 the Christoffel symbols from the metric, which are calculated as


Changing ds by dt

in Equation (1.4), we get


that corresponds to trajectories parametrized by physical time.

3. Curvature and Chaos. Jacobi-Levi-Civita Equations (JLC)

The sectional curvature, at the point and tangent plane on the Riemannian manifold at this point (Gaussian curvature in dimension two), determines the velocity of spreading for the geodesics starting at tangent to.

The evolution for this geodesic deviation obeys the JLC equations


with geodesic in, the covariant derivative along the curve and

the Riemann-Christoffel curvature tensor [30] .

In our case, and the geodesic deviation vector field along the fiducial trajectorythen


with the stability tensor, in a local chart and


The matrix is symmetric and defines a selfadjoint mapping relative to as

where, in a local chart and can be expressed as

and the associated quadratic form

with the inner product, that can be diagonalized by solving

with eigenvalues and, and the corresponding eigenvectors and

being mutually orthogonal. These eigenvectors when normalized define a convenient basis obtaining the new convenient moving frame in terms of the canonical local chart of as follows


The equations (1.7) can be easily expressed in this basis being parallel transported along solving the dynamical equations


with we solve (1.11) using the Taylor expansion

. (1.12)

If the vector field is expressed as we can write (1.9)


with the scalar curvature for metric. The evolution for is as much lineal and the stability is determined by the scalar (from now on referred as) along the fiducial curve.

Fixing a constant length for the trajectory, we define the dynamical magnitude Stability Geometric Indicator (SGI) as


with initial conditions. We will use this indicator to compare the mutual relative stabilities for different trajectories.

4. Isomerization System LiNC/LiCN

4.1. 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 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 is the distance from the Li atom to the NC center of mass, is the NC distance, is the angle defined by the R and CN vectors, and are the corresponding conjugate momenta, and

and represent the Li-NC and NC reduced masses, respectively.

The potential energy surface 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 a.u., rad and and LiNC where a.u.,

rad and and a saddle point (hyperbolic point) at a.u., rad and

. The minimun energy path (MEP) connecting the two isomer wells has been plotted as a dashed line superimposed to the potential energy surface in Figure 1. It can be written as the following Fourier expansion:

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

The energy profile associated to the MEP is given in the right panel of Figure 1.

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, , as the sectioning plane [34] . This PSOS is made an area preserving map defining a new set of canonical variables as:

Now the PSOS is defined by the condition and is obtained from the Hamiltonian conservation

corresponding to the negative branch in the calculation of the second-order equation obtained. We represent PSOS values in the interval rad taking into account the symmetry of the molecular system.

4.2. Jacobi Metric

Considering the Hamiltonian of the system (1.15) then4

and the Riemannian manifold is given by and


. (1.17)

The Christoffel symbols for the metric (1.17) are calculated using (1.5), obtaining


Trajectories parametrized by physical time are calculated integrating the equations (1.6) after replacing the corresponding values of (1.18)


4.3. Jacobi-Levi-Civita (JLC) Equations

Geodesic deviation vector field with along the fiducial trajectory

is used in the JLC equation (1.7). Rewriting equation (1.13) for physical time we obtain


where the scalar curvature for the metric determines the local stability for geodesics around the fiducial one at. If is the dimension of the system it can be written as


with the scalar curvature for the space and the Laplace-Beltrami operator both for metric, in our case for



with the Christoffel symbols of the metric,


replacing (1.24) in (1.23)



The expression for the curvature (1.22) of the mechanical manifold needed for solving Equation (1.20) is

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.

5. 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 all along the boundary i.e. and (kinetic energy equals zero) with a value for energy

. Equation (1.22) is singular at these points so the calculations have been done for an energy

lightly higher with value obtaining the results plotted in Figure 2.

Stable trajectories and chaotic ones for the interval rad for the upper region of the boundary (red colour in upper left graph in Figure 2) are represented in Figure 3 labeled from (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 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 locating a chaotic 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. Values of SGI for different trajectories with initial conditions all along the boundary of the Hill’s region (upper left graph) for the energy, , the interval rad of the upper region (red) is zoomed where a stable periodic trajectory (i) and another from the chain of islands (ii) are marked.

Figure 3. Values of SGI for trajectories with initial conditions in rad for the upper region of the boundary (red colour in upper left graph in Figure 2) at the 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 are plotted in the upper graphs.

6. 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

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

Figure 5. Values of SGI for trajectories with initial conditions in rad at the energy. A stable periodic trajectory is labeled with (1) and a chaotic one with (2). Corresponding trajectories in are plotted in the upper graphs.

region (negative curvature) of the boundary and another chaotic trajectory restricted to convex regions of this boundary.

The curvature of the boundary of the Hill’s region considered as the set level can be calculated

We study the convexity for different set levels and plot the points of curvature zero (change of convexity) as the red curve in Figure 7, the non convex points on the lower region of the boundary for the energy are in the interval rad.

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

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. These trajectories correspond to (4) and (i) respectively in Figure 4.

The two trajectories studied with this energy are labeled as (i) (regular) and (4) (chaotic) in Figure 4 where their corresponding Poincaré sections are plotted. These two trajectories in and the zero curvature for the corresponding set level are plotted in Figure 7. The chaotic trajectory only probes the convex region of the boundary while the regular one bounces non convex points.

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 is quite often to understand the potential energy as curving the underlying space (with the kinetic energy tensor as metric), in the particular case of an underlying space of curvature zero (usually) 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 defines an underlying space positively curved with curvature expressed in (1.26). Thus although the potential energy is removed the space is equally positively curved.

Taking into account the stability is ruled by (JLC) Equations (1.13) where includes the underlying curvature we can deduce the corresponding stabilization effect in the dynamics.

7. Conclusions

In this work, we have studied the vibrational dynamics of the isomerizing LiNC/LiCN molecular system described by a two-dimensional Hamiltonian model including a realistic potential function from a geometrodynamical perspective.

We have considered trajectories starting at the boundary of the Hill’s region and characterized the corresponding dynamics. Once defined the geometrical SGI indicator, we have calculated its value for short times along trajectories with initial conditions in the boundary of the Hill’s region, identifying stable periodic trajectories as local minima of that indicator and the width of the corresponding valley as measure of their stability. We have also identified chaotic trajectories as local maxima around the previous valleys just in the extrema of the corresponding stability region.

The importance of the curvature of the Hill’s region for the stability of trajectories has also been studied. We have found that stable trajectories for our system are mainly located at values of rad and rad in the interval, corresponding to the most convex regions of this boundary. In pinciple, this would support the possibility of reducing our system to a classical billiard problem. However, we found a stable periodic trajectory out of this region. This stable periodic trajectory probes non convex points of the Hill’s region. On the other hand we also find a chaotic trajectory inside the supposed “stability region” bouncing uniquely on convex points of the boundary. Accordingly, the study of stability of trajectories in our molecular system cannot be reduced to a classical billiard problem, despite the Gaussian curvature everywhere positive inside the Hill’s region. This result comes to support the hypothesis of the parametric resonance as the main criterion for the onset of chaos in the case of the chaotic trajectories. However, in the case of the regular trajectories by probing non convex regions of the boundary we understand the underlying curvature (Equation (1.26)) to be the responsible for this lack of unstable behavior. This stabilization mechanism, through the scalar curvature of the underlying space, is currently being under active research in our group.


This research was supported by the Ministry of Economy and Competitiveness-Spain under Contracts No. MTM2012-39101 and ICMAT Severo Ochoa SEV-2011-0087.


  1. Borondo, F. and Benito, R.M. (1995) Dynamics and Spectroscopy of Highly Excited Molecules in Frontiers of Chemical Physics. NATO ASI Series C. Kluwer, Dordrecht.
  2. Skokos, Ch. (2010) The Lyapunov Characteristic Exponents and Their Computation. Lecture Notes in Physics, 790, 63.
  3. Lega, E., Guzzo, M. and Froeschlé, C. (2003) Detection of Arnold Diffusion in Hamiltonian Systems. Physica D, 182, 179-187.
  4. Barrio, R. (2005) Sensitivity Tools vs. Poincaré Sections. Chaos, Solitons and Fractals, 25, 711.
  5. Dumas, H.S. and Laskar, J. (1993) Global Dynamics and Long-Time Stability in Hamiltonian Systems via Numerical Frequency Analysis. Physical Review Letters, 70, 2975.
  6. Skokos, Ch., Antonopoulos, Ch., Bountis, T.C. and Vrahatis, M.N. (2004) Detecting Order and Chaos in Hamiltonian Systems by the SALI Method. Journal of Physics A, 37, 6269.
  7. Benitez, P., Losada, J.C., Benito, R.M. and Borondo, F. (2013) Analysis of the Full Vibrational Dynamics of the LiNC/ LiCN Molecular System. Progress and Challenges in Dynamical Systems, 54, 77-88.
  8. Vergel, A., Benito, R.M., Losada, J.C. and Borondo, F. (2014) Geometrical Analysis of the LiCN Vibrational Dynamics. A Stability Geometrical Indicator. Physical Review E, 89, Article ID: 022901.
  9. Arranz, F.J., Borondo, F. and Benito, R.M. (2000) Transition from Order to Chaos in a Floppy Molecule: LiNC/LiCN. Chemical Physics Letters, 317, 451.
  10. Borondo, F., Zembekov, A.A. and Benito, R.M. (1995) Quantum Manifestations of Saddle-Node Bifurcations. Chemical Physics Letters, 246, 421.
  11. Borondo, F., Zembekov, A.A. and Benito, R.M. (1996) Saddle-Node Bifurcations in the LiNC/LiCN Molecular System: Classical Aspects and Quantum Manifestations. Journal of Chemical Physics, 105, 5068.
  12. Borondo, F., Zembekov, A.A. and Benito, R.M. (1997) Semiclassical Quantization of Fragmented Tori: Application to Saddle-Node States of LiNC/LiCN. Journal of Chemical Physics, 107, 7934.
  13. Losada, J.C., Estebaranz, J.M., Benito, R.M. and Borondo, F. (1998) Local Frequency Analysis and the Structure of Classical Phase Space of the LiNC/LiCN Molecular System. Journal of Chemical Physics, 108, 63.
  14. Borondo, F., Losada, J.C. and Benito, R.M. (2001) Scars in Molecular Vibrations and Spectra of LiCN. Foundations of Physics, 31, 147-163.
  15. Losada, J.C., Benito, R.M., Arranz, F.J. and Borondo, F. (2002) Frequency Map Analysis and Scars in Molecular Vibrations. International Journal of Quantum Chemistry, 86, 167-174.
  16. Eisenhart, L.P. (1929) Dynamical Trajectories and Geodesics. Annals of Mathematics, 30, 591.
  17. Cipriani, P. and Bari, M. (1998) Finsler Geometric Local Indicator of Chaos for Single Orbits in the Hénon-Heiles Hamiltonian. Physical Review Letters, 81, 5532.
  18. Horwitz, L., Zion, Y.B., Lewkowics, M., Schiffer, M. and Levitan, J. (2007) Geometry of Hamiltonian Chaos. Physical Review Letters, 98, Article ID: 234301.
  19. Zion, Y.B. and Horwitz, L. (2007) Detecting Order and Chaos in Three-Dimensional Hamiltonian Systems by Geometrical Methods. Physical Review E, 76, Article ID: 046220.
  20. Zion, Y.B. and Horwitz, L. (2008) Applications of Geometrical Criteria for Transition to Hamiltonian Chaos. Physical Review E, 78, Article ID: 036209.
  21. Zion, Y.B. and Horwitz, L. (2010) Controlling Effect of Geometrically Defined Local Structural Changes on Chaotic Hamiltonian Systems. Physical Review E, 81, Article ID: 046217.
  22. Pettini, M. (1993) Geometrical Hints for a Nonperturbative Approach to Hamiltonian dynamics. Physical Review E, 47, 828.
  23. Casetti, L., Livi, R. and Pettini, M. (1995) Gaussian Model for Chaotic Instability of Hamiltonian Flows. Physical Review Letters, 74, 375.
  24. Casetti, L. and Pettini, M. (1993) Analytic Computation of the Strong Stochasticity Threshold in Hamiltonian Dynamics Using Riemannian Geometry. Physical Review E, 48, 4320.
  25. Casetti, L., Clementi, C. and Pettini, M. (1996) Riemannian Theory of Hamiltonian Chaos and Lyapunov Exponents. Physical Review E, 54, 5969.
  26. Bari, M. and Cipriani, P. (1998) Geometry and Chaos on Riemann and Finsler Manifolds. Planetary and Space Science, 46, 1543.
  27. Casetti, L., Pettini, M. and Cohen, E.G.D. (2000) Geometric Approach to Hamiltonian Dynamics and Statistical Mechanics. Physics Reports, 337, 237-241.
  28. Cerruti-Sola, M. and Pettini, M. (1995) Geometric Description of Chaos in Self-Gravitating Systems. Physical Review E, 51, 53.
  29. Cerruti-Sola, M. and Pettini, M. (1996) Geometric Description of Chaos in Two-Degrees-of-Freedom Hamiltonian Systems. Physical Review E, 53, 179.
  30. Carmo, M.P. (1992) Riemannian Geometry. Birkh?user, Boston.
  31. Carmo, M.P. (1976) Differential Geometry of Curves and Surfaces. Prentice Hall, Englewood Cliffs.
  32. Saa, A. (2004) On the Viability of Local Criteria for Chaos. Annals of Physics, 314, 508-516.
  33. Esser, R., Tennyson, J. and Wormer, P.E.S. (1982) An SCF Potential Energy Surface for Lithium Cyanide. Chemical Physics Letters, 89, 223-227.
  34. Benito, R.M., Borondo, F., Kim, J.H., Sumpter, B.G. and Ezra, G.S. (1989) Comparison of Classical and Quantum Phase Space Structure of Nonrigid Molecules, LiCN. Chemical Physics Letters, 161, 60-66.


1Projection of the accessible region in the configuration space on the position coordinates.

2We will use the Einstein’s summation convention for repeated indexes, so


4, i.e. is the inverse matrix of.