Non-affine geometrization can lead to nonphysical instabilities

Non-affine geometrization can lead to nonph ysical instabilities Eduardo Cuerv o-Reyes 1, 2 , ∗ and Ramis Mo v assagh 3, 4 , † 1 Swiss F e der al L ab or atories of Materials Scienc e and T e chnolo gy, ¨ Ub erlandstr asse 129, CH-8600 D ¨ ub endorf Switzerland 2 Swiss F e der al Institute of T e chnolo gy (ETH), CH-8093 Z¨ urich, Switzerland 3 Dep artment of Mathematics, Northe astern University, Boston, MA, 02115 4 Dep artment of Mathematics, Massachusetts Institute of T e chnolo gy, Cambridge, MA, 02139 (Dated: June 21, 2021) Geometrization of dynamics consists of representing tra jectories by geo desics on a configuration space with a suitably defined metric. Previously , efforts were made to show that the analysis of dynamical stabilit y can also be carried out within geometrical framew orks, b y measuring the broadening rate of a bundle of geo desics. Tw o known formalisms are via Jacobi and Eisenhart metrics. W e find that this geometrical analysis measures the actual stability when the length of any geo desic is prop ortional to the corresponding time interv al. W e prov e that the Jacobi metric is not alw ays an appropriate parametrization by showing that it predicts chaotic b eha vior for a system of harmonic oscillators. F urthermore, we show, by explicit calculation, that the corresp ondence b et w een dynamical- and geometrical-spread is ill-defined for the Jacobi metric. W e find that the Eisenhart dynamics corresp onds to the actual tangen t dynamics and is therefore an appropriate geometrization scheme. I. ST ABILITY OF HAMIL TONIAN SYSTEMS A. Ly apunov exponent Man y ph ysical systems are w ell represen ted b y the time ev olution of the co ordinates ( q i ( t ), i = 1 . . . N ) of parti- cles with inertia matrix a ij ( q ), moving under the influ- ence of the p oten tial V ( q ). The q i ( t ), and the momenta p i ( t ) = a ij ˙ q j ( t ), where repeated indices are summed o v er, are found by integrating Hamilton’s equations ˙ p i = − ∂ H ∂ q i , ˙ q i = ∂ H ∂ p i , (1) with the Hamiltonian H ( q , p ) = 1 2 a ij ( q ) p i p j + V ( q ) , (2) where a ik a kj = δ i j . Since, in practical situations, ini- tial conditions ( q ( t 0 ) = q 0 , and p ( t 0 ) = p 0 ) are sub- ject to uncertain ties, one often needs to kno w the exten t to which a pair of tra jectories, which are infinitesimally close at some time, will remain close at later times; i.e., the extent to whic h the dynamics is stable. The stability of dynamical systems is also relev an t in the con text of statistical metho ds, as it has strong implications on the ergo dicit y of isolated systems. A Hamiltonian system in general may hav e stable or unstable tra jectories, depending on the region of the phase space where the dynamics takes place [1–3]. The degree of instability of a set of tra jectories can b e quan- ∗ Electronic address: eduardo.cuervoreyes@empa.ch † Electronic address: ramis.mov@gmail.com tified b y the Ly apunov exp onen t λ ≡ lim t →∞ 1 t ln  || x ( t ) || || x (0) ||  , (3) where || x ( t ) || is the norm of the 2 N -dimensional v aria- tion of the phase-space tra jectory . Note that there are in principle 2 N ly apunov exp onen ts. Since it is sufficient for detecting c haos, we will only consider the largest and simply refer to it as “the” lyapuno v exp onen t. The v ari- ation x ( t ) is obtained from the solution of the linearized v ariation of the equations of motion (the tangen t dynam- ics) [4–7] ∆ ˙ p i = −  ∆ q j ∂ ∂ q j + ∆ p j ∂ ∂ p j  ∂ H ( q , p ) ∂ q i , (4a) ∆ ˙ q i =  ∆ q j ∂ ∂ q j + ∆ p j ∂ ∂ p j  ∂ H ( q , p ) ∂ p i . (4b) λ quantifies the rate of separation of infinitesimally close tra jectories, and λ ≤ 0 in the stable regimes with equal- it y holding for conserv ativ e systems. T o find regions of the phase space corresponding to (un)stable dynamics, an exhaustive computation of tra jectories, and the cor- resp onding tangent dynamics, is in general necessary . It is, therefore, desirable to ha v e a metho d which could pre- dict the stability , from static “geometrical” prop erties. B. Geometrization of dynamics Geometrization of dynamics consists of represen ting the physical tra jectories in time b y geo desics on a mani- fold with a suitable metric. Geometrization has the the- oretical upshot in that it provides an alternative frame- w ork for analyzing the dynamics. Geometrization do es not necessarily pro vide a computationally efficien t adv an- tage ov er the standard methods, y et it has pro v ed useful 2 in quantifying stability from the curv ature of the corre- sp onding metric (see for example [9]). Inspired by the general theory of relativity , Eisenhart [8] prop osed a geometrical description of classical dynam- ics, for systems which can b e describ ed by the least ac- tion principle. Eisenhart’s metric is defined on the en- larged configuration space-time with N + 2 co ordinates x ν = { q 0 ≡ t, q i , q N +1 : i = 1 . . . N } , by the differential arc length ( ds ) ds 2 = − 2 V ( q )( dq 0 ) 2 + a ij dq i dq j + 2 dq 0 dq N +1 . (5) F or an y geo desic, the extra coordinate has the solution q N +1 = κ 2 2 t + C 0 − Z t 0 L dt 0 , (6) where L is the Lagrangian, and C 0 and κ are arbitrary real constan ts. With this metric, physical motions sat- isfy an affine p ar ametrization ds 2 = κ 2 ( dt ) 2 (i.e., the arc length of any geo desic is prop ortional to the time elapsed). C. Geometrical v ersion of the stability analysis It has b een prop osed that stability can b e quantified based on a lyapuno v exp onen t in a geometrical frame- w ork [10]. Analogous to ξ T ≡ ∆ q , one defines the vector field of geo desic spread ξ i G ( s ) ≡ ∆ q i ( s ) (7) as a v ariation of geodesics at constan t arc length. ξ G sat- isfies the linearized Jacobi-Levi-Civita (JLC) equations D 2 ξ i G ds 2 + R i j km dq j ds ξ k G dq m ds = 0 , (8) where D ds , and R i j km are the co v arian t deriv ative along geo desics, and the Riemann-Christoffel curv ature tensor [11] 1 , resp ectiv ely . The stabilit y is quantified b y the ge- ometrical ly apuno v exponent [10] λ G ≡ lim s →∞ 1 s ln  || ξ G ( s ) || || ξ G (0) ||  . (9) 1 In terms of the Christoffel sym b ols Γ k lj = 1 2 g km  g lm,j + g mj,l − g lj,m  respectively . The cov ariant differentiation and the curv ature tensor are defined as Dξ i = dξ i + Γ i lj ξ l dq j , R i j lk = ∂ Γ i j k ∂ q l − ∂ Γ i j l ∂ q k + Γ m j k Γ i ml − Γ m j l Γ i mk . Eq.(8) can b e rewritten as d 2 ξ k ds 2 + 2Γ k lj dq l ds dξ j ds + Γ k lm,j dq l ds dq m ds ξ j = 0 . This definition formally excludes the v ariations with resp ect to the momen ta of the initial conditions. Never- theless, in practice most of the dynamical features per- taining to stability can b e explored by v ariations in co- ordinates alone. With Eisenhart’s metric, the spacial comp onen ts of Eq. (8) are equiv alent to Eqs. (4) [12]; therefore, λ G ≡ λ . In w ords, the Ly apuno v exponent obtained with Eisenhart’s metric is equal to the one obtained with the tangen t dy- namics. F urthermore, P ettini et.al. [13, 14] hav e sho wn that the (in)stability of some Hamiltonian systems can b e quantified by means of the Ricci curv ature, rendering unnecessary the tedious in tegration of the equations of motion. The formalism based on Eisenhart’s ( N + 2)- dimensional space also has the adv antage of being appli- cable to systems with time-dep enden t Hamiltonians. Sev eral works in the last tw o decades ha v e b een dedi- cated to sho w that, in the N -dimensional configuration- space, the geometrical-stabilit y analysis can b e also car- ried out employing the kinetic energy metric, also known as the Jac obi metric [15–17] ( g J ) ij ≡ 2[ E − V ( q )] a ij ( q ) . (10) The interest in Jacobi’s metric was partly due to the fact that the resulting scalar curv ature, Ricci’s curv ature, and sectional curv atures are positive, in many systems of in- terest [12, 16, 17]. This seemed to supp ort the h yp othesis that negative curv ature is not the fundamen tal source of instabilit y , and that the non-negativ e oscillating curv a- ture leads the systems to chaos through parametric res- onance [10]. In order to test the v alidit y of the geometrical approach within Jacobi’s metric, λ J G has b een computed for some mo del Hamiltonians and these results hav e b een com- pared with those from tangen t dynamics. The evidence has not b een conclusive. JCL equations within the Ja- cobi metric are generally cumbersome; in fact, they are tangen t dynamics equations with some extra terms ha v- ing no clear ph ysical interpretation [15]. V ery few exam- ples of exact numerical in tegration of Eq.(8) hav e b een presen ted [12, 17], and there is no intuitiv e understand- ing of the results [12, 15]. An approximate version of Eq.(8) was used for stability analysis of a large system of self-gra vitating particles, and a surprising suppression of c haos was observ ed, with increasing num b er of particles [16]. This unexp ected b eha vior w as assumed to originate from the mathematical approximations that w ere made for obtaining an equation for the dynamics of the scalar || ξ G || [16, 17]. Szydlo wski et.al. [18, 19] hav e already p oin ted out some inadequacies of the Jacobi metric. A central argu- men t has b een that the curv ature tensor b ecomes singular at the b oundaries of the configuration space. There, the kinetic energy v anishes, and there is an infinite n um ber of geo desics, none of which corresp onds to a ph ysical tra jec- tory . One may susp ect that for many degrees of freedom this framew ork could giv e meaningful results, since the probabilit y of reac hing the b oundaries (i.e., the system 3 coming to a full halt) is very lo w. How ever, this has not b een made rigorously . Another limitation of Jacobi’s metric is that, since it dep ends critically on the total energy ( E ), only those v ariations that do not c hange E m ust be considered. In this paper w e sho w that the geo desic spread and the tangen t dynamics fields are equiv alent when the geodesic length is only a function of time (if any pair of geo desics, whic h co v er equal time interv als, ha ve equal arc length). W e apply the formalism to a system at the onset of c haos, and to a trivially stable system (i.e., system of harmonic oscillators). In these cases w e found that Jacobi’s metric predicts an unstable dynamics for stable systems, ev en for N  2, where the kinetic energy do es not v anish. The geo desic-spread within the Jacobi metric suffers from non-ph ysical parametric resonance, and that this is re- lated to the fluctuations of the kinetic ener gy . As such, endo wing the configuration space with Ja- cobi metric do es not alwa ys provide an appropriate mea- sure for the calculation of the Lyapuno v exp onen t, and the consequent analysis of dynamical instability . The man uscript is organized as follows. In section I I the rela- tionship b et ween tangent dynamic and geo desic spread, and the differences in the corresp onding equations are deriv ed and discussed. In section I II stabilit y analysis by means of Jacobi metric is compared to the results from tangen t dynamics for a t wo-dimensional system. In sec- tion IV a similar comparison is made employing a stable and analytically soluble system (harmonic oscillators). General conclusions are dra wn in section V. I I. RELA TIONSHIP BETWEEN ξ T AND ξ G A v ariation (∆ q i ) of the co ordinates q i (measured be- t ween tw o tra jectories) can b e written as ∆ q i = ξ i T + ˙ q i ∆ t, (11) where ξ i T is the v ariation at a fixed time, and the second term accounts for a time v ariation. Let the arc length of a geodesic be parametrized b y s = Z T 0 F ( q , ˙ q ) dt. (12) The difference in the arc le ngth of tw o tra jectories, up to first order in co ordinates and time v ariations, is ∆ s =  F − ˙ q i ∂ F ∂ ˙ q i  ∆ t | T + ∂ F ∂ ˙ q i ∆ q i | T − ∂ F ∂ ˙ q i ∆ q i | 0 + Z T 0  ∂ F ∂ q i − d dt ∂ F ∂ ˙ q i  ξ i T dt . (13) The r.h.s. in tegral v anishes along any geo desic b ecause the terms in paren thesis are the equations of motion. Re- mem b er that the latter are obtained from the extremal condition (∆ s = 0) with resp ect to v ariations that keep the boundary conditions unchanged. Th us, the ∆ s be- t ween t w o ph ysical tra jectories with slightly differen t ini- tial conditions is ∆ s =  F − ˙ q i ∂ F ∂ ˙ q i  ∆ t | T + ∂ F ∂ ˙ q i ∆ q i | T − ∂ F ∂ ˙ q i ∆ q i | 0 ; (14) and it simplifies to different forms depending on the v ariational formalism. In geometrical formalisms with tra jectory-indep enden t arc length (such as Eisenhart’s approac h), the difference in the arc length of tw o neigh- b oring geo desics reduces to ∆ s = F ∆ t , (15) where F turns out to be constant. As such, v ariations at constant t corresp ond to v ariations at constan t s . In v ariational formalisms where s is minimal with resp ect to v ariations with unconstrained time (e.g., within Jacobi metric), the term m ultiplying ∆ t | T v anishes; i.e., F − ˙ q i ∂ F ∂ ˙ q i ≡ 0 . (16) Remem b er that in the case of the Jacobi metric F = L + E , and v ariations are tak en at constant energy E , with unconstrained time. In such cases, s is tra jectory dep enden t, and the difference in the arc-length of neigh- b oring geo desics is ∆ s = ∂ F ∂ ˙ q i ∆ q i | T − ∂ F ∂ ˙ q i ∆ q i | 0 . (17) Consequen tly , a geo desic spread at constant s , ξ i G ≡ ∆ q i | ∆ s =0 , is only p ossible if ∂ F ∂ ˙ q i ξ i G is constant in time, or equiv alently if d dt  ∂ F ∂ ˙ q i ξ i G  = 0 . (18) If w e assume that Eq.(18) can b e fulfilled, then, one can m ultiply Eq.(11) b y ∂ F ∂ ˙ q i and find the corresp onding time mismatc h ∆ t b et w een the tw o geo desics with equal arc, and we arrive at the following transformation b et w een ξ G and ξ T : ξ G = M ξ T + C  ∂ F ∂ ˙ q i ˙ q i  − 1 ˙ q , (19) 4 where C = ∂ F ∂ ˙ q i ξ i G , and M j i = δ j i − ∂ F ∂ ˙ q i ˙ q j / ∂ F ∂ ˙ q k ˙ q k . (20) It is easy to sho w that M 2 = M , and T r( M ) = N − 1; th us, M is a pro jector of rank ( N − 1). The tangen- tial comp onen t of ξ G is not w ell-defined in the ab o ve relations. A t first sight, one could think that a suitable c hoice of co ordinates might isolate this comp onen t. How- ev er, the geodesic spread at constant arc length turns out to b e ill-defined. Eq.(18) cannot b e fulfilled and the con- stan t C do es not exist. T o show this, we pro ceed as follo ws. Carrying out the differen tiation in Eq.(18), it takes the form ξ i G d dt  ∂ F ∂ ˙ q i  + ∂ F ∂ ˙ q i ˙ ξ i G = 0 ; (21) whic h, together with the equations of motion, give ξ i G ∂ F ∂ q i + ∂ F ∂ ˙ q i ˙ ξ i G = 0 . (22) A t the same time, since Eq.(16) is an identit y that has to b e fulfilled by any geo desic, its linearized v ariation must b e identically zero; i.e., 0 ≡ ∆  F − ˙ q i ∂ F ∂ ˙ q i  = ∂ F ∂ q i ξ i G − ˙ q i ∆  ∂ F ∂ ˙ q i  . (23) Satisfying both Eq.(23) and Eq.(22) w ould require that ∆  ˙ q i ∂ F ∂ ˙ q i  = 0 . (24) Eq.(24) is a very strong condition whic h, within the Ja- cobi metric, means that the v ariations at constant arc length m ust also leav e the kinetic energy unc hanged. F or the sak e of simplicity , and since it does not affect the gen- eralit y of the analysis, let us rep eat the abov e deriv ation for the Jacobi metric assuming that the elements of the inertia matrix ( a ij ) are constant. The fact that the Jacobi metric depends critically on E requires that we restrict the v ariations to a constan t en- ergy h ypersurface. This is actually equiv alent to Eq.(23); and from this condition, the v ariations must satisfy a ij ˙ q i ˙ ξ j G + V ,k ξ k G = 0 . (25) On the other hand, Eq.(22) gives a ij ˙ q i ˙ ξ j G − V ,k ξ k G = 0 . (26) The tw o constraints (Eq.(25) and Eq.(26)) require that ξ lea ves b oth total energy and kinetic energy unchanged at any time. In other w ords, ˙ ξ has to b e orthogonal to the momen tum, and ξ orthogonal to the force. W e will sho w now that these t wo conditions are incompatible. But first, let us show that the tangent dynamics do es allo w ξ to remain in the constant-energy h ypersurface. The equations of the tangent dynamics, for constan t a ij , are the set of equations ¨ ξ n T + a nk V ,kl ξ l T = 0 . (27) The set Eq.(25) constrains the v ariation to a constant energy h ypersurface. The dynamics of ξ must allow this to hold ov er time; therefore the time deriv ative of Eq. (25) must b e iden tically zero. By p erforming this deriv a- tiv e, and taking Hamilton’s equations in to accoun t, one obtains p n  ¨ ξ n T + a nk V ,kl ξ l T  = 0 . (28) Because of Eq.(27), all co efficien ts m ultiplying the mo- men ta in Eq. (28) are zero, and therefore the ab o ve con- dition is satisfied for any ph ysical tra jectory . Within the Jacobi metric, this is not the case. The corresp onding JLC equations are ¨ ξ n G + a nk V ,kl ξ l G = − 1 T  a nm V ,m { a ij ˙ q i ˙ ξ j G + V ,l ξ l G } − ˙ q n { V ,il ˙ q i ξ l G + V ,j ˙ ξ j G + 1 T V ,i ˙ q i V ,l ξ l G }  , (29) where T = E − V . If w e could guarantee that the v aria- tions ( ξ G and ˙ ξ G ) do not alter the total energy , the first term in curly braces would v anish, and the equations w ould reduce to ¨ ξ n G + a nk V ,kl ξ l G = ˙ q n T  V ,il ˙ q i ξ l G + V ,j ˙ ξ j G + 1 T V ,i ˙ q i V ,l ξ l G  . (30) Comparing Eqs.(27) with Eqs.(30), w e see that the latter con tains extra terms (the r.h.s.). These terms can b e 5 rewritten as ˙ q n d dt (( E − V ) − 1 ∆ V ). Th us, the differential equations for ξ G w ould b ecome iden tical to the equations of the tangent dynamics only when d dt V ,k ξ k G E − V ! ≡ 0 . (31) Ho wev er, Eq.(31) is not protected b y the dynamics de- riv ed from Eqs.(30). Namely , this equality cannot b e de- riv ed from the equations of motion. As a result, even if w e choose initial conditions for ξ G and ˙ ξ G , which satisfy δ T = δ E = 0, Eq.(31) will not hold o ver time. The work of Sospedra and co-work ers[12] apparently sho ws that, in contrast to our results, the comp onen t of ξ G in the direction of the tra jectory can b e decoupled from the system of equations, and it do es not acceler- ate. How ev er, their result is a consequence of replacing the co v arian t deriv ative of a quan tity that is not a true scalar (the pro jections of ξ G on a given basis) by a simple deriv ative. W e will show numerically in the next sections that the extra terms in Eq. (30) in troduce non-ph ysical instabili- ties. But let us first make a simple heuristic analysis to get the feeling of the problem. As before, we take the tangen t dynamics as a reference of v alidit y to compare against. When the energy of a system is close enough to its v alue at a minimum of the potential, the motion is re- stricted to small p eriodic oscillations around this equi- librium p oin t; the dynamics is stable. W e can expand the potential up to second order in q . Then, all V ,kl are constan t, and the tangent dynamics (from Eq.(27)) giv es oscillatory solutions for ξ T (therefore, λ = 0), in corre- sp ondence to a dynamically stable system. If we consider no w tra jectories with higher energy , higher orders in the q dep endence of V m ust b e taken in to account. This ma y driv e the system tow ards chaotic b ehavior. V ,kl are q -dependent, and oscillate at harmonics of the frequen- cies of the system. They may pro duce parametric reso- nance in the tangent dynamics, and result in λ > 0 [20]. Th us, parametric resonance in ξ -dynamics is relev an t at the on-set of c haos [10]. F or ev en higher energies, the tra- jectories may cross hyperb olic p oin ts b et w een the minima of the p oten tial; the system may visit most regions of the phase space, and the tra jectories are no longer p erio dic. No w w e repeat the analysis using Jacobi metric. Because of the first order deriv ativ es ( V ,j ), and the ˙ q -dependence of the the r.h.s. of Eqs.(30), the co effi- cien ts in these linear equations are time-dep enden t, even when V ( q ) is quadratic in q . These terms oscillate at frequencies that are harmonics of the fundamen tal mo des of the system, and so they may create parametric reso- nance [20], resulting in p ositive λ G ’s for a stable system. A t the actual onset of c haos, when higher orders in the q -dependence of V are imp ortan t, real parametric reso- nance app ears from the second term on the l.h.s.. It is not p ossible to discriminate the false exp onen tial diver- gence from the physical one, when the exp onen t is com- puted from the ev olution of ξ G according to Eqs.(30). A p ositiv e exp onen t ma y b e obtained in b oth, stable and unstable regimes. I II. A TWO-DIMENSIONAL EXAMPLE. PHYSICAL AND UNPHYSICAL INST ABILITIES The first representativ e example that we tak e is a tw o- dimensional system describ ed by the Hamiltonian H = P 2 R 2 µ 1 +  1 2 µ 1 R 2 + 1 2 µ 2 r 2 e  P 2 θ + V ( R , θ ) . (32) This represen ts the energy of a particle mo ving with re- sp ect to the center of mass of a rigid dimer of length r e = 3 . 0271 ˚ A , when the total angular momen tum is J = 0; it has been the sub ject of inv estigation due to its in terest for molecular dynamics [12, 21–26]. R and θ are p olar coordinates, with p olar axis running along the dimer. The interaction b et w een the parti- cle and eac h member of the dimer is represented by a Morse p oten tial V ( r ) = D  1 − e − α ( r − d )  2 , with param- eters D = 40 . 75 cm − 1 , α = 1 . 56 ˚ A − 1 , and d = 4 . 36 ˚ A [12]. The resulting potential surface has tw o equiv a- len t minima at ( R = p r 2 e / 4 + d 2 ; θ = ± π / 2). These minima are connected by tw o equiv alen t trenches go- ing around the dimer, with saddle p oin ts at ( R = d + 1 α ln(cosh( αr e )) − 1 α ln(cosh( αr e / 2)); θ = 0 , π ). The rel- ativ e height E saddle − E min = 2 D [1 + coth 2 ( αr e / 2)] − 1 is ≈ 40 . 66 cm − 1 . The reduced masses are taken as µ 1 = 18 am u, and µ 2 = 64 amu [12]. In the follo wing, all relev ant quan tities are expressed in pow ers of cm. Energy , linear momen ta, and λ are in cm − 1 , time and distances are in cm and angular momen tum is dimensionless. FIG. 1: P oincare surface section for ∆ E = 6 . 5 W e hav e performed an exhaustive integration of tra- jectories, and we ha ve analyzed P oincare surface sections (PSS) for different v alues of the total energy . W e present in Fig.1 the PSS θ - P θ , for ∆ E = E − E min = 6 . 5 cm − 1 . A t this energy , there are wide regions that con tain sta- ble limit-cycles, and there are also regions of unstable dynamics. W e chose initial conditions at the edge of a 6 stable region, and computed the exp onen ts ( λ G , and λ ), for ∆ E = 6 . 5 cm − 1 , and for several energies ab o v e and b elo w. This should allo w us to observe how the on-set of chaos is detected by the geometrical exp onen t, and to compare it with the results from tangen t dynamics. FIG. 2: The Lyapuno v exp onen t obtained from the tangent dynamics of tra jectories at the edge betw een stable and un- stable regions in Fig.1, for sev eral v alues of the total energy . FIG. 3: The conv ergence of the Lyapuno v exp onen t for dif- feren t energies using Eisenhart metric (tangent dynamics). W e plot the Lyapuno v exponent obtained with the Eisenhart metric (actually , tangent dynamics) v ersus the total energy , in Fig.2. The on-set of c haos seems to tak e place around ∆ E = 4 . 3 cm − 1 , as indicated b y p ositiv e v alues of λ . F or a reason that will b ecome clear so on, it is interesting to lo ok at the early-time evolution of Y ( t ) = 1 t log( || ξ ( t ) || || ξ (0) || ) for different energies, ranging from v ery stable (∆ E = 0 . 5 cm − 1 ) to unstable (∆ E = 15 . 5 cm − 1 ) regimes. Fig.3 sho ws the stabilization of Y around the corresp onding limiting v alue. W e intended to cal- culate the Lyapuno v exp onen t using Jacobi metric, and w e faced serious memory ov erflo ws, due to a explosiv e gro wth of k ξ G k . Its time evolution is exponential even for tra jectories in the stable regime (e.g., ∆ E = 2 . 5 cm − 1 ), as it can b e seen in Fig.4. W e extracted from the Y FIG. 4: Kinetic Energy and “conv ergence” of the Ly apuno v exp onen t for ∆ E =2.5 cm − 1 using Jacobi metric. FIG. 5: Zo oming in on Fig.4: Kinetic Energy and the corre- sp onding “kicks” in ξ for ∆ E =2.5 cm − 1 . dynamics a very short time window whic h is shown in Fig.5. In contrast to the results from tangent dynam- ics, one can see the ”kicks” pro duced in the evolution of Y , associated to the oscillations of the kinetic energy . Y con verges to about 130 cm − 1 ; in con tradiction to the v anishing exp onen t obtained with Eisenhart metric, and also in contradiction with PSS calculated in the corre- sp onding energy range. These results ga ve us the idea to lo ok at a mo del where one could quantify the oscillations of the kinetic energy with a minimal error and at minimal cost, in order to searc h for a correlation betw een the positive λ G and the fluctuations of T . This is the sub ject of the next section. W e like to comment here that during the review of this manuscript w e became aw are 2 of other results ob- tained b y Motter and co-work ers, which are in goo d agreemen t with ours [27–29]. They found that the Lya- puno v exp onen t could be a reliable indicator of stabil- 2 W e thank Prof. R. Montgomery for ha ving kindly suggested these v aluable references within the review of our manuscript. 7 it y as long as the time-reparametrization do es not cre- ate singularities in the inv arian t measure. In turn, sin- gular reparametrizations can shift the Ly apunov exp o- nen t. The Jacobi metric is actually an example of singu- lar reparametrizations. In a v ery recen t work [30] it has b een also sho wn that, within the Jacobi metric, there are tra jectories whic h fail to minimize the arc length. IV. λ G FR OM JACOBI METRIC VS FLUCTUA TIONS OF THE KINETIC ENER GY. Here w e use a paradigm of stabilit y: a system of in- dep enden t harmonic oscillators. W e apply the stability analysis using Jacobi metric; any evidence of c haos can then b e understo od as a failure of the metho d. W e in- v estigate whether there is a relation b et ween λ G and the fluctuation of the kinetic energy . The time-independent Hamiltonian H ( q , p ) = 1 2  δ ij p i p j + ω 2 δ ij q i q j  , (33) with p i = δ ik ˙ q k is our basic mo del. The solutions to the equations of motion are q k ( t ) = C k cos( ω t + θ k ) , (34) where C k and θ k dep end on the initial conditions, and w e c hose them as C k = 1, and θ k = k 2 π f N , k = 1 , ..., N . As b efore, N stands for the num ber of degrees of freedom. The phases, θ k , are homogeneously distributed on a frac- tion ( f ) of 2 π (phase circle). These settings allow us to find an analytical expression for the fluctuation √ σ , and to tune it v arying the v alues of f and N . Using Eq.(34), w e obtain for T , and for its normalized v ariance ( σ ) T = N  ω C 2  2  1 − √ 2 σ cos(2 ω t + 2 π f N + 1 N )  , (35) σ ≡ h T 2 i − h T i 2 h T i 2 =  sin(2 π f ) √ 2 N sin(2 π f / N )  2 . (36) σ decreases with increasing N , having the limit σ → | sin(2 π f ) 2 π f √ 2 | as N → ∞ . When all of the oscillators are in phase ( f = 0), σ tak es its maximum v alue and T b ecomes zero ev ery ∆ t = π ω . JLC equations with Eisenhart metric (and the equa- tions of the tangen t dynamics) take the form d 2 ξ k T dt 2 + ω 2 ξ k T = 0 , (37) whereas with the Jacobi metric, d 2 ξ k G dt 2 + ω 2 ξ k G + ω 2 T δ lj "  q k ˙ q l − q l ˙ q k  dξ j G dt +  ω 2 q k q l − ˙ q k ˙ q l − ω 2 T δ im ˙ q i q m ˙ q k q l  ξ j G # = 0 . (38) Eq.(37) gives λ = 0 for any initial condition. On the other hand, when tra jectories (34) are substituted in Eq.(38), the latter tak e the form d 2 ξ k G dt 2 + ω 2 ξ k G + ω I k j dξ j G dt + ω 2  J k j + K k j  ξ j G = 0 , (39) with the couplings I k j , J k j and K k j giv en b y I k j = ω 2 C 2 T sin( θ k − θ j ) , (40a) J k j = ω 2 C 2 T cos(2 ω t + θ k + θ j ) , (40b) K k j = − ω 4 C 4 2 T 2 sin(2 ω t + 2 π f N + 1 N ) × (40c) [sin(2 ω t + θ k + θ j ) − sin( θ k − θ j )] . Eqs.(39) hav e a t ypical structure from whic h parametric resonance ma y arise [20]. The basic frequency of the os- cillators is one half of the frequency at which the terms of the “restoring force matrix” oscillate (which is the ideal condition for a first order parametric resonance). Ho wev er, since the comp onen ts ξ k G ha ve b ecome coupled through the terms I k j , J k j and K k j , one could still doubt whether a normal mo de resonates with the fluctuating parameters. Despite the simplicit y of the c hosen system, the an- alytical integration of these equations seems quite c hal- lenging. Therefore, we lea v e the analytical ev aluation of the Ly apuno v exp onen t for future w ork. W e obtained λ G from the numerical in tegration of Eq.(39), using several v alues of f , and N . In Fig.6 we sho w λ G vs √ σ . Note that, although we v aried N and f indep enden tly , and σ dep ends on both of these parameters, w e obtained a smo oth curve for λ G vs. √ σ . This confirms our h yp oth- esis that, within Jacobi metric, λ G gro ws with σ , and it do es not measure the actual stability of the physical sys- tem. Actually , un-physical parametric resonance is not the only manifestation of the failure of this metho dology . The scalar curv ature ( K ) of the manifold is K = N − 1 8( E − V ) 3  4( E − V ) ∇ 2 V − ( N − 6) |∇ V | 2  , (41) whic h takes negative v alues in some regions of the acces- 8 FIG. 6: λ G vs. √ σ for N = 2 + j 2 with j = 1 , . . . , 14 and f = 0 . 05 , 0 . 1 , . . . , 0 . 45. sible configuration space, for an y system with N > 6. F or our simple example of harmonic oscillators, the curv ature can be rewritten as K = ω 2 ( N − 1) 4( E − V ) 3 [2 N ( E − V ) − ( N − 6) V ] , (42) whic h, for N m uch greater than 6, will be negativ e at any p oin t with V > 2 E / 3. Eq.(41) mak es eviden t that the sign of K is not an indicator of the (in)stability of the underlying dynamical system. W e would lik e to finish this discussion sho wing a special (and illustrative) case of the ab o v e example, where the differential equations for ξ G are simple; the solutions are quite telling b y mere insp ection. Let us take t w o identical one dimensional simple har- monic oscillators with coordinates x and y , resp ectiv ely . T ak e the initial conditions such that they both hav e non- v anishing amplitudes, and a phase difference of π 2 . With this choice, the kinetic energy never v anishes. This is formally equiv alen t to a single tw o dimensional oscillator with a central p oten tial V ( r ) = 1 2 ω 2 r 2 . The orbits are ellipses; a circle corresp onding to the special case where b oth v ariables oscillate with equal amplitude. The cir- cular tra jectory has a constant kinetic energy . This ex- ample also allows a transparent w a y to control the initial conditions of ξ G , so to keep the total energy unchanged (whic h has b een said to b e a k ey to obtain a meaningful stabilit y analysis). The tra jectories can b e parametrized in polar coordinates as r 2 ( t ) = R 2 + ∆ 2 cos(2 ω t ) , (43) ˙ θ 2 ( t ) = r − 4 ω 2 ( R 4 − ∆ 4 ) , (44) with the total energy E and the angular momen tum L giv en b y E = R 2 ω 2 , (45) L 2 = ω 2 ( R 4 − ∆ 4 ) . (46) R 2 , and ∆ 2 are the a verage of r 2 ( t ), and its oscilla- tion amplitude, respe ctiv ely . The latter is directly re- lated to the fluctuations of the kinetic energy , since T = ω 2 2 [ R 2 − ∆ 2 cos(2 ω t )]. The constants satisfy ∆ 2 ≤ R 2 , with ∆ = 0 corresp onding to the circular orbit, and ∆ = R to the undesired one dimensional case. By making use of conserv ation la ws, one can reduce the tw o coupled equations of ξ G to a single v ariable problem. Eisenhart’s metric, which is equiv alent the tangen t dy- namics, giv es ¨ ξ T +  V ,rr + 3 L 2 r 4  ξ T = 0 . (47) One can show, b y differen tiation of Newton’s equation, that the solutions of eq.(47) are prop ortional to ˙ r ( t ), and are stable, as exp ected. F or the harmonic p oten tial, one can also show that ξ T ∝  r − (2 L 2 /E ) r − 1  . On the other hand, in p olar co ordinates, the com- p onen ts of the metric tensor for the Jacobi metric are g rr = 2[ E − V ], g rθ = 0, and g θθ = 2[ E − V ] r 2 . The differen tial equation for the radial comp onen t of ξ is ¨ ξ G +  V ,rr + 3 L 2 r 4  ξ G =  V ,rr  2 − L 2 r 2 ( E − V )  (48) + V 2 ,r  1 E − V − 2 L 2 r 2 ( E − V ) 2  + V ,r 3 L 2 r 3 ( E − V )  ξ G . The extra terms on the righ t hav e no ph ysical meaning. The first thing to note is that the r.h.s. of eq.(48) v an- ishes for tra jectories with ∆ = 0. Thus, we observ e here again that the equation from the geometrical formalism with Jacobi metric is equiv alent to equation from the tangen t dynamics only for tra jectories with constan t ki- netic energy . The reader might argue that although the equations are different, the stability analysis might still b e equiv alen t. T o sho w that this is not the case, it is enough to ev aluate eq.(48) for 2 ω t = (2 n + 1) π 2 , to ob- tain ¨ ξ G + ω 2  4 − 7 ∆ 4 R 4  ξ G = 0 , for t = (2 n + 1) π 2 . (49) The effective restoring force is negative in time inter- v als con taining 2 ωt = (2 n + 1) π 2 , for an y tra jectory with ∆ 4 R 4 > 4 7 . The length of these “explosiv e” interv als is larger for larger v alues of ∆. Then, one finds again an infinite n umber of stable tra jectories for which the ge- ometrical approach with Jacobi metrics predicts an un- stable b eha vior. This example shows that not only the effectiv e frequencies b ecome time dep enden t, they take imaginary v alues, perio dically . In summary , the computation of the Lyapuno v exp o- nen t within the kinetic energy metric has m ultiple short- comings. In addition to the resulting equations b eing m uch more in v olved than those corresponding to the tan- gen t dynamics, small kinetic energies make these equa- tions numerically unstable. Moreov er, the metho d by definition, is restricted to conserv ative systems. W e also 9 find that the sign of the curv ature, and in general, the sign of any elemen t of the curv ature tensor, is not a go od measure of stabilit y . Lastly , the geode sic spread generally fails to describ e the stability of the underlying ph ysical system, by introducing un-physical parametric resonance and negativ e restoring forces. V. CONCLUSIONS In this paper, we hav e sho wn that the vector field of geo desic spread ξ G , and the tangen t dynamic vector field ξ T are equiv alent when the arc length measured along an y geo desic is proportional to the time in terv al. This is the case with the Eisenhart metric. When this is not fulfilled, as in the Jacobi metric, the geo desic spread is ill-defined. In the Jacobi metric, the equations of motion satisfied by ξ G con tain extra terms, that do not app ear in ξ T . These terms do not seem to hav e a clear ph ysi- cal meaning, and can b e resp onsible for the non-physical parametric resonance seen. F urthermore, they cause un- stable modes (imaginary frequencies) in stable systems. The Ly apuno v exp onen t calculated within the geomet- rical formalism with Jacobi metric is equiv alen t to the tangen t dynamics if the total kinetic ener gy is c onserve d , whic h is unrealistic for interacting systems. The time ev olution of the geo desic spread is not compatible with the mo v emen t in a constan t-energy h yp ersurface. Using tw o represen tativ e examples, we demonstrated that the geometrical Lyapuno v exp onen t, calculated with Jacobi metric, is correlated with the fluctuation of the kinetic energy , irresp ectiv e of the actual dynamical sta- bilit y of the system. This implies that, in a statistical system of N interacting particles, where the fluctuations in the kinetic energy decrease as N increases, the formal- ism based on Jacobi metric would falsely predict that the system b ecomes “less” chaotic as N → ∞ . This, w e b eliev e, is the reason why an unexp ected reduction of c haos with the increase of N was seen in simulations of a self-gra vitating system [16]. W e conclude that great care m ust be taken when using the Jacobi metric to derive stability results, esp ecially in computing the Lyapuno v exponents. In this resp ect the Eisenhart metric, or equiv alen tly , the tangen t dynamics remain as a simpler metho d that do not suffer from sin- gularities. It w ould be in teresting to study whether other non-affine parameterizations hav e similar drawbac ks. [1] H. P oincar´ e, Les M ´ etho des Nouvel les de la M´ echanique Celeste Blanchard, Paris, vol. 3 (1987). [2] V. I. Arnold, Mathematic al methods for classical me chan- ics (MIR, Moscow, 1976). [3] N. S. Krylov, Works on the F undations of Statistic al Me- chanics (Princeton Univ. Press, Princeton, 1979) [4] M. T o da, Phys. Lett. A 48 , 335 (1974). [5] P . Brumer and J. W. Duff, J. Chem. Phys. 65 , 3566 (1976). [6] J. W. Duff and P . Brumer, J. Chem. Phys . 67 , 4898 (1977). [7] R. Koslo v and S. Rice, J. Chem. Phys. 74 , 1947 (1981). [8] L. P . Eisenhart, Ann. Math. 30 , 591 (1929) [9] R. Montgomery , Erg. Theory and Dyn. Systems, v. 25; 921-947 (2005) [10] M. P ettini, Phys. Rev. E 47 828-987 (1993). [11] T. Levi-Civita, Ann. Math. 97, 291, (1926). [12] R. Sosp edra-Alfonso, L. V elazquez, and J. Rubay o- Soneira, Chem. Phys. Lett. 375 , 261 (2003). [13] L. Casseti, R. Livi and M. Pettini, Ph ys. Rev. Lett. 74 375-378 (1995) [14] L. Caiani, L. Casetti, C. Clementi and M. Pettini, Phys. Rev. Lett. 79 4361-4364 (1997). [15] M. Cerruti-Sola, R. F ranzosi and M. Pettini, Phys. Rev. E 56 , 4872-4875 (1997). [16] M. Cerruti-Sola and M. Pettini, Phys. Rev. E 51 53-64 (1995). [17] M. Cerruti-Sola and M. P ettini, Phys. Rev. E 53 179-188 (1996). [18] M. Szydlowski, M. Heller, W. Sasin, J. Math. Ph ys. 37 1 345, (1996). [19] M. Szydlo wski, J. Math. Ph ys. 40 7 345, (1999). [20] L. D. Landau and E. M. Lifshitz Me chanics (1976). [21] O. Roncero, J. A. Beswick, N. Halberstadt, P . Villarreal and G. Delgado-Barrio, J. Chem. Phys. 92 , 3348 (1990). [22] D. M. Willb erg, M. Gutmann, J. J. Breen and A. H. Zew ail, J. Chem. Ph ys. 96 , 198 (1992). [23] M. Gutmann, D. M. Willberg and A. H. Zewail, J. Chem. Ph ys. 97 , 8037 (1992). [24] M. Gutmann, D. M. Willberg and A. H. Zewail, J. Chem. Ph ys. 97 , 8048 (1992). [25] M. Grueb ele and A. H. Zewail, J. Chem. Ph ys. 98 , 883 (1993). [26] J. Rubay o-Soneira, A. Garc ´ ıa-V ela, G. Delgado-Barrio and P . Villarreal, Chem. Phys. Lett. 243 , 236 (1995). [27] A. E. Motter, Phys. Rev. Lett. 91 , 231101 (2003). [28] A. E. Motter and A. Saa, Ph ys. Rev. Lett. 102 , 184101 (2009). [29] K. Gelfert and A. E. Motter, arXiv:1010.1791 (2010). [30] R. Mon tgomery , arXiv: 1407.7163 (2014).