Coherent Riemannian-geometric description of Hamiltonian order and chaos with Jacobi metric

TL;DR

This paper defends the use of the Jacobi metric in geometrizing Hamiltonian dynamics, showing it can reliably describe chaos and stability without inconsistencies, countering recent criticisms.

Contribution

It demonstrates that the Jacobi metric provides a consistent geometric framework for analyzing Hamiltonian chaos, clarifying misconceptions about its limitations.

Findings

Jacobi metric accurately describes chaos in Hamiltonian systems.

Observed instabilities with Jacobi metric are artefacts, not fundamental flaws.

Supports geometric approach as a valid method for stability analysis.

Abstract

By identifying Hamiltonian flows with geodesic flows of suitably chosen Riemannian manifolds, it is possible to explain the origin of chaos in classical Newtonian dynamics and to quantify its strength. There are several possibilities to geometrize Newtonian dynamics under the action of conservative potentials and the hitherto investigated ones provide consistent results. However, it has been recently argued that endowing configuration space with the Jacobi metric is inappropriate to consistently describe the stability/instability properties of Newtonian dynamics because of the non-affine parametrization of the arc length with physical time. To the contrary, in the present paper, it is shown that there is no such inconsistency and that the observed instabilities in the case of integrable systems using the Jacobi metric are artefacts.

Equations

$$L(\dot{q},q)=\frac{1}{2}{\overline{g}}_{ij}(q){\dot{q}}^i{\dot{q}}^j-V(q),$$

$${\dot{q}}^i(t)=\frac{\mathrm{d}{q}^i}{\mathrm{d}t}(t)$$

$$S[\mathbf{q}(t)]={\int }_{t_0}^{t_1}L(\dot{\mathbf{q}}(t),\mathbf{q}(t))\mathrm{d}t$$

$$\delta S[\gamma ]=0\text{with}\gamma (t_0)=a\text{and}\gamma (t_1)=b.$$

$${\overline{\nabla }}_{\dot{\gamma }}\dot{\gamma }=-\overline{\mathrm{grad}}V$$

$$\delta l[\tilde{\gamma }]=0\text{with}l[\tilde{\gamma }(s)]={\int }_{s_0}^{s_1}\sqrt{{\tilde{g}}_{ij}\frac{\mathrm{d}{q}^i}{\mathrm{d}s}\frac{\mathrm{d}{q}^j}{\mathrm{d}s}}\mathrm{d}s\tilde{\gamma }(s_0)=a\tilde{\gamma }(s_1)=b.$$

$$H(\dot{q},q)=\left({\dot{q}}^i\frac{\partial L}{\partial {\dot{q}}^i}-L\right)=\frac{1}{2}{\overline{g}}_{ij}(q){\dot{q}}^i{\dot{q}}^j+V(q)$$

$$W=\frac{1}{2}{g}_{ij}{\dot{q}}^i{\dot{q}}^j=H(\dot{q},q)-V(q)$$

$$W(q(t;E))=E-V(q)$$

$$\begin{array}{rl}{S}_E[q(t;E)]={\int }_{t_0}^{t_1}L(\dot{q}(t;E),q(t;E))\mathrm{d}t& ={\int }_{t_0}^{t_1}\left[E+2W(q(t;E))\right]\mathrm{d}t\\ & =(t_1-t_0)E+{\int }_{t_0}^{t_1}2W(q(t;E))\mathrm{d}t.\end{array}$$

$$\begin{array}{rl}{S}_E[q(t;E)]& ={\int }_{t_0}^{t_1}2W(q(t;E))\mathrm{d}t\\ & ={\int }_{t_0}^{t_1}\sqrt{2W(q(t;E))}\sqrt{2W(q(t;E))}\mathrm{d}t={\int }_{t_0}^{t_1}\sqrt{2\left[E-V(q)\right]}\sqrt{{\bar{g}}_{ij}(q){\dot{q}}^i{\dot{q}}^j}\mathrm{d}t=\\ & ={\int }_{t_0}^{t_1}\sqrt{2\left[E-V(q)\right]{\bar{g}}_{ij}{\dot{q}}^i{\dot{q}}^j}\mathrm{d}t={\int }_{t_0}^{t_1}\sqrt{g_{ij}{\dot{q}}^i{\dot{q}}^j}\mathrm{d}t={\int }_{t_0}^{t_1}\sqrt{{\left(\frac{\mathrm{d}s}{\mathrm{d}t}\right)}^2}\mathrm{d}t\\ & ={\int }_{s_0}^{s_1}\mathrm{d}s\end{array}$$

$$g_{ij}:=2W(q){\overline{g}}_{ij},W(q):=E-V(q)$$

$$\mathrm{d}{s}^2=g_{ij}\mathrm{d}{q}^i\mathrm{d}{q}^j=2\left[E-V(q)\right]{\overline{g}}_{ij}\mathrm{d}{q}^i\mathrm{d}{q}^j=4[W(q){]}^2\mathrm{d}{t}^2.$$

$$d{s}^2=4W^2(q)d{t}^2.$$

$$g_J(\dot{q}(s),\dot{q}(s))=2W\frac{d{q}^i}{ds}{\delta }_{ij}\frac{d{q}^j}{ds}=1$$

$$\frac{d{q}^i}{dt}{g}_{ij}\frac{d{q}^j}{dt}=2W.$$

$$Y^i:=\frac{\mathrm{d}{q}^i}{\mathrm{d}s}$$

$$X^i:=\frac{\mathrm{d}{\tilde{q}}^i}{\mathrm{d}t}=\frac{\mathrm{d}s}{\mathrm{d}t}\frac{\mathrm{d}{q}^i}{\mathrm{d}s}=2W{Y}^i$$

$${\nabla }_{Y}Y=0$$

$${\nabla }_{\alpha X}\left(\alpha X\right)=\alpha \left({\nabla }_X\alpha \right)+{\alpha }^{2}X=0$$

$${\nabla }_{X}X=-{\nabla }_X(\log \alpha )X=X{\nabla }_X\left(\log 2W\right).$$

$$\frac{{\mathrm{d}}^2{q}^i}{\mathrm{d}{t}^2}+{\Gamma }_{jk}^i\frac{\mathrm{d}{q}^j}{dt}\frac{\mathrm{d}{q}^k}{dt}={\delta }_k^i\frac{\mathrm{d}{q}^k}{\mathrm{d}t}\frac{\mathrm{d}{q}^j}{\mathrm{d}t}{\partial }_j\log (2W).$$

$$\begin{array}{rl}{\Gamma }_{jk}^i=& {\overline{\Gamma }}_{jk}^i+\frac{1}{2}\left[{\delta }_j^i{\partial }_k\log (2W)+{\delta }_k^i{\partial }_j\log (2W)-g_{jk}{\overline{\mathrm{grad}}}^i\log (2W)\right]\end{array}$$

$${\overline{\nabla }}_{X}X=\frac{{\mathrm{d}}^2{q}^i}{\mathrm{d}{t}^2}+{\overline{\Gamma }}_{jk}^i\frac{\mathrm{d}{q}^j}{dt}\frac{\mathrm{d}{q}^k}{dt}=\frac{1}{2}{\overline{\mathrm{grad}}}^i(2W)=-{\overline{\mathrm{grad}}}^{i}V$$

$$\frac{{\nabla }^2{J}^i}{d{s}^2}+R_{jkl}^i\frac{d{q}^j}{ds}{J}^k\frac{d{q}^l}{ds}=0$$

$$R_{jkl}^i={\partial }_k{\Gamma }_{jl}^i-{\partial }_l{\Gamma }_{jk}^i+{\Gamma }_{lj}^m{\Gamma }_{km}^i-{\Gamma }_{kj}^m{\Gamma }_{ml}^i$$

$$\begin{array}{rl}\frac{{\nabla }^2{X}^i}{d{s}^2}& =\frac{1}{2W}\frac{\nabla }{dt}\left(\frac{1}{2W}\frac{\nabla X^i}{dt}\right)\\ & =\frac{1}{4W^2}\frac{{\nabla }^2{X}^i}{d{t}^2}+\frac{1}{4W}\frac{\nabla X^i}{dt}\frac{d}{dt}\left(\frac{1}{W}\right)\end{array}$$

$$\begin{array}{r}\frac{{\nabla }^2{J}^i}{d{s}^2}+R_{jkl}^i\frac{d{q}^j}{ds}{J}^k\frac{d{q}^l}{ds}=0\mapsto \frac{{\nabla }^2{J}^i}{d{t}^2}+R_{jkl}^i\frac{d{q}^j}{dt}{J}^k\frac{d{q}^l}{dt}=\frac{\nabla J^i}{dt}\frac{d}{dt}\ln W.\end{array}$$

$$\begin{array}{r}{\nabla }_{\mathbf{\xi }}^{2}J+R(J,\mathbf{\xi })\mathbf{\xi }=0.\end{array}$$

$$\begin{array}{rl}{\nabla }_{\mathbf{\xi }}J& =Y\\ {\nabla }_{\mathbf{\xi }}Y& =R_{\mathbf{\xi }}J.\end{array}$$

Full text

Hamiltonian Chaos, Differential Geometry

I Introduction

Elementary tools of Riemannian differential geometry can be successfully used to explain the origin of chaos in Hamiltonian flows or, equivalently, in Newtonian dynamics. Natural motions of Hamiltonian systems can be viewed as geodesics of the configuration-space manifold $\displaystyle M$ equipped with the Riemannian metric $\displaystyle g_{J}$, known as the Jacobi metric (or kinetic energy metric). The stability/instability properties of such geodesics can be investigated by means of the Jacobi–Levi-Civita (JLC) equation for geodesic spread. It has been shown that chaos in physical geodesic flows does not stem from hyperbolicity of $\displaystyle M$: phase space trajectories/geodesics are destabilized both by regions of negative curvature and by parametric instability caused by positive curvature varying along the geodesics marco ; physrep ; book . Another remarkable fact is that the JLC equation written for a geometrization of Hamiltonian systems in an enlarged configuration space-time endowed with a metric due to Eisenhart eisenhart yields the standard tangent dynamics equation commonly used in numerical computations of the largest Lyapunov exponent (LLE). These two different geometric framework have been proven to give the same information about order and chaos and about the strength of chaos as well. This has been checked in the case of two-degrees of freedom Hamiltonian systems, for the Hénon-Heiles model and for two coupled quartic oscillators, respectively cerruti1996geometric ; rick , and in the case of a large number of degrees of freedom (from 150 up to 1000) cerruti1997lyapunov . It has been found that the JLC equation stemming from the Jacobi metric gives exactly the same quantitative results of the tangent dynamics equation.

This notwithstanding, in Ref. cuervo2015non it has been argued that the non-affine parametrization of the arc-length with time in configuration space endowed with the Jacobi metric leads to nonphysical instabilities. More precisely, the JLC equation written for the Jacobi metric seems to give chaos also for a system of harmonic oscillators. Moreover, such alleged non-physical instabilities are found to be stronger in systems with few degrees of freedom. Such an argument seems to radically exclude the use of Jacobi metric in configuration space to consistently investigate Hamiltonian chaos; in fact, the non-affine parametrization of the arc length with respect to physical time is an unavoidable consequence of the way Jacobi metric is derived from Maupertuis’ least action principle. Some mathematical works montgomery2014s ; giambo2014morse ; giambo2015normal , partly motivated by the results reported in Ref.cuervo2015non , have investigated the behaviour of the geodesics in configuration space endowed with the Jacobi metric near the so called Hill’s boundaries, i.e. the regions in configuration space where $\displaystyle E=V(\mathbf{q})$ and where the Jacobi metric is singular ($\displaystyle g_{J}=0$). In particular, all these works emphasized the phenomenon of geodesic reflection near Hill’s boundaries and the relevance of these reflections to characterize periodic orbits, with, among the others, an original proposal dating back to Ref. seifert1948periodische . Finally, other authors yamaguchi2001geometric have tackled the geometrization of Hamiltonian dynamical systems by lifting the Jacobi metric from configuration space $\displaystyle M$ to its cotangent bundle $\displaystyle T^{*}M$, i.e. to the phase space. In this framework JLC equations are rewritten as a system of first order linear differential equations on the tangent bundle $\displaystyle TT^{*}M$ of phase space: this allows to identify the adequate degrees of freedom to compute the Geometric Largest Lyapunov Exponent(GLLE). Within this framework, the GLLE (JLC equations) and LLE (tangent dynamics) are found in very good agreement to characterize the chaotic regime (strong/weak chaos) of the Hénon-Heiles model. All these studies, together with the manifest contradiction among the outcomes of Ref.cuervo2015non and those of Refs. cerruti1996geometric ; rick have motivated the present work. The present paper is organized as follows. In Section II we briefly discuss some aspects of the construction of the Jacobi metric, with emphasis on the consequences of non-affine parametrization of the arc length with time for the JLC equation, and on the presence of boundaries where the metric is singular. In Section III, the JLC equation is rewritten by introducing a parallel transported frame for which explicit expressions are then given in Subsections III.1 and III.3 for $\displaystyle N=2$ and $\displaystyle N=3$, respectively. Then in Subsections III.2 and III.4 the results of the corresponding numerical simulations are reported for two and three (resonant) harmonic oscillators, respectively, using the most critical values of the parameters which, according to Ref.cuervo2015non , should yield non physical instabilities. It is shown that this is not the case. In Section IV, we discuss how the measure concentration phenomenon, which takes place at a large number of degrees of freedom, completely removes any problem making unnecessary even resorting to a parallel transported frame. Some conclusions are drawn in Section V.

II Effects of Non-Affine Parametrization of the arc length with Jacobi Metric

Among the different possibilities of rephrasing Newtonian dynamics in geometric terms, as reported in physrep ; book , the Jacobi metric in configuration space leads to the mathematically richest structure [$\displaystyle(M,g_{J})$ is geodesically complete in the sense of the Hopf-Rinow theorem]. Let us consider a thorough investigation of the geometrization of Newtonian dynamics by means of $\displaystyle(M,g_{J})$ for systems described by a Lagrangian of the form

[TABLE]

where

[TABLE]

$\displaystyle\overline{g}_{ij}$ are the components of the kinetic energy metric on configuration space $\displaystyle M$ with the associated Levi-Civita connection $\displaystyle\overline{\nabla}$ specified by the Christoffel coefficients $\displaystyle\overline{\Gamma}^{i}_{jk}$, and $\displaystyle V(q)$ is the potential energy. It is well known that, given a chart $\displaystyle({U},\phi)$ on the configuration space $\displaystyle M$ and a curve $\displaystyle\gamma:I\subset\mathbb{R}\longrightarrow M$, the natural motions $\displaystyle\phi(\gamma(t))=\boldsymbol{\gamma}(t)=\boldsymbol{q}(t)$ are the class of curves that make stationary the action functional

[TABLE]

on the class of curves with $\displaystyle\boldsymbol{q}(t_{0})=\boldsymbol{a}$ and $\displaystyle\boldsymbol{q}(t_{1})=\boldsymbol{b}$ fixed, i.e.

[TABLE]

The Newton equations are derived from the Euler-Lagrange equations, i.e.

[TABLE]

where $\displaystyle\overline{\mathrm{grad}}^{i}f(q)=\overline{g}^{ik}\partial_{k}f$ is the gradient. Natural motions $\displaystyle\gamma(t)$ belong to a class of curves of configuration space satisfying the ”physical” variational principle (4), therefore these curves can be identified with geodesics of configuration space which also satisfy a variational principle but in this case of a ”geometrical” kind. In fact, geodesics $\displaystyle\tilde{\gamma}(s)$ are curves of a Riemannian manifold $\displaystyle(M,\tilde{g})$ endowed with a metric $\displaystyle\tilde{g}$ that makes stationary the length functional between two fixed points, i.e.

[TABLE]

One possible way to provide an identification of natural motion with geodesics is provided by the introduction of the Jacobi metric on a subspace of configuration space. Let us consider the total energy function

[TABLE]

which is obviously a conserved quantity along the natural motions, i.e. $\displaystyle d{H}(\dot{q}(t),q(t))/dt=0$. Since the Lagrangian is a homogeneous function in $\displaystyle\dot{q}^{i}$ it follows that $\displaystyle H(\dot{q},q)=2W-L$ where

[TABLE]

is the kinetic energy. If we consider a class of isoenergetic trajectories $\displaystyle q(t;E)$ in configuration space, that is having the same total energy value $\displaystyle H(\dot{q}(t;E),q(t;E))=E$, as the kinetic energy is non negative, the trajectories of the system in configuration space are confined in the region $\displaystyle\mathcal{M}_{V<E}=\left\{q\in M|V(q)<E\right\}$. Moreover, for isoenergetic trajectories the kinetic energy can be expressed as a function of the coordinates, i.e.

[TABLE]

thus the action functional of Eq. (3) can be rewritten in the form

[TABLE]

as we are interested in the variational principle and $\displaystyle t_{0},t_{1},E$ are fixed quantities, the first term in the last equality can be neglected. The integral in Eq.(10) can be interpreted as a length integral in configuration space, in fact

[TABLE]

where the new metric

[TABLE]

called also Jacobi metric, has been introduced with the associated the arc-length element

[TABLE]

A central point of the following discussion is related to the non-affine parametrization of the arc-length $\displaystyle s$ with respect to the physical time $\displaystyle t$, which is clearly a necessary consequence of the construction of Jacobi metric. Moreover, we observe that the Jacobi metric is related to kinetic energy metric $\displaystyle\overline{g}_{ij}$ through a conformal rescaling via a factor proportional to the kinetic energy $\displaystyle[E-V(q)]$ preserving the signature of the metric only in the interior of $\displaystyle\mathcal{M}_{V\leq E}$, i.e. $\displaystyle\mathcal{M}_{V<E}=\mathcal{M}_{E}=\left\{q\in M|V(q)<E\right\}$, the so called Hill’s region. By endowing the region $\displaystyle\mathcal{M}_{E}$ with the Jacobi metric $\displaystyle g$ the natural motions with fixed energy $\displaystyle E$ are the same as geodesics $\displaystyle\gamma(s)$ of the manifold $\displaystyle(M_{E},g)$.

This approach has remarkable consequences, as it is discussed in the following. In fact, the geometric description of Newtonian dynamics, identifying the solutions of Newton equations with the geodesics of suitable Riemannian manifolds, provides a powerful conceptual and mathematical framework to study the stability/instability of dynamics in terms of the stability properties of a geodesic flow, described by the geodesic spread equation relating stability/instability with geometry.

The standard observable to define the presence of dynamical chaos and to measure its strength is the largest Lyapunov exponent, and, as we will see in the following, in the geometrical framework a Geometrical Lyapunov exponent can be defined.

Now, a key point in derivation of Jacobi metric from Maupertuis’ principle is to set the relation between the physical time $\displaystyle t$ and the arc-length $\displaystyle s$ as

[TABLE]

It follows that a generic geodesic parametrized by the arc-length $\displaystyle s$, i.e. $\displaystyle\dot{\boldsymbol{q}}(s)=\{\dot{q}^{i}(s)\}_{i\in[1,n]}$ has a unit velocity

[TABLE]

while if parametrized with respect to the physical time

[TABLE]

So, in general, the physical time is a non-affine parametrization of the geodesics of Jacobi metric. We derive in what follows the effect of the reparametrization (14) on the geodesic equation. Let us introduce the vector fields

[TABLE]

and

[TABLE]

defined along the geodesic $\displaystyle q^{i}(s)=q^{i}(s(t))=\tilde{q}^{i}(t)$. Using the definition of geodesic for the vector field $\displaystyle Y$

[TABLE]

it is possible to derive the equations for the vector field $\displaystyle X=\alpha Y$ (with $\displaystyle\alpha=(2W)^{-1}$)

[TABLE]

that implies

[TABLE]

In a natural coordinate system $\displaystyle\{q^{i}\}_{i=1,...,N}$ the equation (20) reads

[TABLE]

The Christoffel symbols in Jacobi metric take the form

[TABLE]

where $\displaystyle\overline{\Gamma}^{i}_{jk}$ are Christoffel symbols of the kinetic energy metric $\displaystyle g_{ij}$ and $\displaystyle\overline{\mathrm{grad}}^{i}f:=g^{ik}\partial_{i}f$ is the gradient with respect to kinetic energy metric. Substituting (22) in (21) and using (15) we obtain

[TABLE]

i.e. the Newton’s equations of the dynamical system. As already mentioned above, dynamical chaos can now be investigated by means of the equation for the geodesic spread describing the stability of a geodesic flow of a Riemannian manifold. The geodesic spread is measured by a vector field $\displaystyle J$ which locally gives the distance between nearby geodesics. This vector field evolves along a reference geodesic according to the Jacobi-Levi Civita equation which, in components, reads

[TABLE]

where

[TABLE]

is the Riemann curvature tensor. In Jacobi metric we have

[TABLE]

that substituted into (24) yields to

[TABLE]

Every non-affine parametrization in Jacobi metric can be reduced to a relation between the physical time and the arc-length parameter of the form $\displaystyle t\mapsto ds=f(q)dt$, where $\displaystyle f(q)$ is a function of the point and generates a term proportional to the derivative of $\displaystyle f(q)$ into each equations. Thus, the JLC equation written for the Jacobi metric is not invariant for time reparametrization. We can rephrase the main point raised by the work in Ref.cuervo2015non as attributing to the right hand-side of Eq.(27) the origin of chaotic-like instabilities even for integrable systems, that is, the origin of non-physical artifacts. In Ref. cuervo2015non it has been surmised that the larger the fluctuations of kinetic energy $\displaystyle W$ the more dramatic the occurrence of non-physical instabilities stemming from Eq.(27). Remarkably, the geometrization of Newtonian dynamics in an enlarged configuration space-time equipped with the Eisenhart metric tensor marco yields the JLC equation in the form of the Tangent Dynamics Equation which is commonly used to compute the Largest Lyapunov Exponent (LLE) marco ; book , therefore the authors of Ref. cuervo2015non claim that this is the only consistent geometrization of Newtonian dynamics to investigate chaos, while the Jacobi metric would be unsuitable for the same task.

II.1 Geometrical Lyapunov exponent

We conclude the present Section by giving a definition of a Geometrical Lyapunov exponent. Of course, the starting point is the JLC which, in intrinsic notation, is

[TABLE]

By defining $\displaystyle Y=\nabla_{\boldsymbol{\xi}}J$ and $\displaystyle R_{\boldsymbol{\xi}}J:=R(\boldsymbol{\xi},J)\boldsymbol{\xi}$ the above equation becomes

[TABLE]

Then, by putting

[TABLE]

the JLC equation reads

[TABLE]

A Geometrical Lyapunov exponent, in analogy with the definition of the standard Lyapunov exponent, can be defined after having expressed as a function of physical time the solution of Eq. (30) as

[TABLE]

where the norm of $\displaystyle\mathcal{J}$ is

[TABLE]

Let us note that throughout the literature, marco , book , cerruti1996geometric and cerruti1997lyapunov , the geometrical Lyapunov exponent has been expressed as a function of physical time $\displaystyle t$, whereas in cuervo2015non the Geometrical Lyapunov exponent is defined as a function of the arc-length. As a final comment, let us remark that the existence of many different frameworks to rephrase Newtonian dynamics in geometric terms book can lead, a-priori, to different quantitative evaluations of the strength of chaos, to the contrary, the use of different geometric frameworks must lead to the same qualitative description of the stability/instability properties of the dynamics. For example, the transition from weak to strong chaos must be at least qualitatively and coherently reproduced in any geometric framework, as is for example actually shown for high dimensional Hamiltonian flows in Ref. cerruti1997lyapunov .

III Parallel Transported Frame for a system of $\displaystyle N$ particles in Jacobi Manifold

Let us now work out the JLC equation for a parallel transported orthonormal frame along a reference geodesic. The advantage of this representation with respect to the use of standard local coordinates is that by making parallel transported frames to ”incorporate” the geodesics reflection, when they approach the Hill’s boundaries in configuration space, eliminates a source of artefacts in the numerical solution of the JLC equation. In fact, the sharp reflection of geodesics close to the Hill’s boundary of a mechanical manifold would require a prohibitively high numerical precision to avoid the introduction of an error amplification mimicking chaos even for integrable systems. To the contrary, with respect to a parallel transported frame - ”incorporating” the geodesic reflection - nearby geodesics are no longer affected by the fake error amplification due to the reflection. As a matter of fact, we will show that the solutions of the JLC equation - written for the Jacobi metric - have the correct physical meaning also in the ”pathological” cases where unphysical instabilities were found in Ref.cuervo2015non .

The parallel transported frame is built by requiring that the covariant derivative $\displaystyle\nabla_{\boldsymbol{\xi}}X$ of all vectors $\displaystyle X$ with respect to the geodesic flow $\displaystyle\boldsymbol{\xi}$ is orthogonal to $\displaystyle\boldsymbol{\xi}$. Let us introduce a reference frame $\displaystyle e_{i}:=\partial/\partial q^{i}$ on the tangent bundle $\displaystyle TM$ and the corresponding dual frame $\displaystyle\theta^{i}:=dq^{i}$ such that $\displaystyle\theta^{i}(e_{j})=\delta^{i}_{j}$. The Jacobi metric tensor is the multi-linear map $\displaystyle\boldsymbol{g}_{J}:TM\times TM\longrightarrow{\mathbb{R}}$ such that, if written with respect to the natural basis, is given by

[TABLE]

the components of which are

[TABLE]

with the differential arc-length

[TABLE]

To build a parallel transported frame we consider a geodesic $\displaystyle\gamma:I\subset{\mathbb{R}}\longrightarrow M$, with tangent vector field $\displaystyle\boldsymbol{\xi}:I\longrightarrow TM$ and the Jacobi vector field $\displaystyle J:I\longrightarrow TM$ such that respect to the natural basis $\displaystyle\{e_{i}\}_{i\in[1,n]}$ they are written as

[TABLE]

where the Jacobi field, by definition, verifies

[TABLE]

We shall show that there exist reference frames parallel transported along any geodesic. Such systems exist if an orthogonal tensor field $\displaystyle\boldsymbol{\Omega}:I\longrightarrow SO_{N}(\gamma(I))$ is defined at each point along the geodesic flow. This tensor field allows to define the required frame where the basis on $\displaystyle TM$ is given by $\displaystyle\{E_{i}\}_{i=i,\cdots,N}$ and dual frame on $\displaystyle T^{*}M$ $\displaystyle\{\Theta^{i}\}_{i=i,\cdots,N}$ and, by definition of dual frame we have $\displaystyle\Theta^{i}(E_{j})=\delta^{i}_{j}$. Therefore, the parallel transported frame is defined by

[TABLE]

In the Jacobi-Levi Civita equation, written in (34), we can define the Ricci tensor along the geodesic flow $\displaystyle Ric(\boldsymbol{\xi}):\gamma(I)\longrightarrow\mathcal{T}^{1}_{1}(\gamma(I))$ by

[TABLE]

A canonical isomorphism exists between the rank-two tensor $\displaystyle Ric(\boldsymbol{\xi})$ and a symmetric matrix $\displaystyle N\times N$, that we shall denote with $\displaystyle\boldsymbol{R}$. For every $\displaystyle\omega\in T^{*}M$ and $\displaystyle X\in TM$, such a matrix is given by

[TABLE]

where $\displaystyle Y_{\omega}:=\boldsymbol{g}^{-1}\omega$. We note the following symmetry properties of $\displaystyle Ric(\boldsymbol{\xi})$

[TABLE]

and by redefining the indices $\displaystyle j\rightleftharpoons i$ and $\displaystyle l\rightleftharpoons k$ the symmetry of $\displaystyle Ric(\boldsymbol{\xi})$ is evident. Hence the associated matrix can be diagonalized at each point along the flow. Let note that the symmetry of Riemann tensor entails

[TABLE]

Then, the tangent vector to the geodesic $\displaystyle\boldsymbol{\xi}$ is an eigenvector of the matrix $\displaystyle\boldsymbol{R}$ associated to $\displaystyle Ric(\boldsymbol{\xi})$ with vanishing eigenvalue. A posteriori, one observes that for harmonic oscillators - geometrized through the Jacobi metric - the eigenbasis of the matrix $\displaystyle\boldsymbol{R}$ coincides with the parallel transported frame. In general, this is not true and we have two different orthogonal matrices, that one which diagonalises $\displaystyle\boldsymbol{R}$ and another one which transforms the natural basis into the parallel transported frame. Fortunately, in the present case, given the natural frame $\displaystyle\{e_{i}\}_{i\in[1,N]}$ on the tangent bundle $\displaystyle TM$, there exists an orthogonal transformation $\displaystyle\boldsymbol{\Omega}:I\longrightarrow SO_{N}(\gamma(I))$ such that

[TABLE]

namely, such that it transforms the natural basis into the parallel transported frame $\displaystyle\{E_{i}\}_{i\in[1,N]}$ and, moreover, such that it diagonalises the matrix $\displaystyle\boldsymbol{R}$:

[TABLE]

This allows to write the Jacobi field with respect to such a basis and thus

[TABLE]

Let us consider the JLC equation (34) and proceed to substitute $\displaystyle J\mapsto\tilde{J}$ then

[TABLE]

where now the repeated indices do not stand for summation. In this way, we obtain $\displaystyle N-1$ second order differential equations in the unknown functions $\displaystyle\{\tilde{J}^{k}\}_{i\in[1,N-1]}$, and the $\displaystyle N$-th equation is that for the vector tangent to the reference geodesic, equation which corresponds to the geodesic equation thus, by denoting with $\displaystyle\tilde{J}^{1}$ the function for this equation, we have

[TABLE]

It is interesting to remark that $\displaystyle\{\Lambda_{k}\}$ are just sectional curvatures, i.e the principal directions of curvature identified by the vectors $\displaystyle\{E_{k}\}_{k\in[1,N]}$, i.e

[TABLE]

Now, passing from the arc-length parameter $\displaystyle s$ to the physical time $\displaystyle t$ the equations (44) become

[TABLE]

whence

[TABLE]

The second equation in LABEL:before_canonical_way can be immediately solved by setting $\displaystyle h(t):=\dot{J}^{1}(t)$ and then tackling the differential equation

[TABLE]

giving

[TABLE]

The other equations can be written in canonical form, namely, by redefining the function $\displaystyle\tilde{J}$ as follows

[TABLE]

and substituting it into (LABEL:before_canonical_way) we obtain

[TABLE]

By choosing $\displaystyle A(t)$ such that the coefficient of $\displaystyle\dot{f}(t)$ vanishes, we get the following conditions for $\displaystyle A(t)$

[TABLE]

thus giving

[TABLE]

The equation (49) becomes

[TABLE]

To apply this procedure to the equations (44), we set

[TABLE]

and for every $\displaystyle k\in[2,N]$ we obtain the final form for the components of the Jacobi-Levi Civita equation

[TABLE]

where the functions in the second term of the l.h.s. can be interpreted as time dependent frequencies

[TABLE]

where

[TABLE]

III.1 Parallel Transported Frame for a system of $\displaystyle 2$ harmonic oscillators

In Refs. cerruti1996geometric ; rick it has been shown that for the Hénon-Heiles model and for two coupled quartic oscillators, respectively, the geometrization through the Jacobi metric perfectly discriminates between ordered and chaotic motions by investigating the stability/instability of geodesics through the JLC equation expressed in a parallel transported frame. In this section, we are going to show that for two harmonic oscillators (of course an integrable system) the norm of the geodesic separation vector remains bounded, in spite of the fluctuations of kinetic energy which, according to the claim of Ref.cuervo2015non , should have entailed apparent instability of the regular motions of this system. The Hamiltonian of this system is

[TABLE]

and the associated Jacobi metric, having set $\displaystyle W(q):=E-V(q_{1},q_{2})$, is

[TABLE]

The Ricci tensor along the geodesic flow is

[TABLE]

The eigenvectors of this matrix are

[TABLE]

with the corresponding eigenvalues:

[TABLE]

With these eigenvectors the matrix for the basis transformation is simply obtained in the form

[TABLE]

so that the JLC equation for the components of the parallel transported Jacobi vector field $\displaystyle\tilde{J}$ is written as

[TABLE]

which, written for the physical time $\displaystyle t$ and with the notations of the preceding Section, read

[TABLE]

where

[TABLE]

III.2 Numerical results

In Ref.cuervo2015non , it has been claimed that there exists a set of initial conditions for which the JLC equations written for the Jacobi metric lead to unstable solutions even in the case of two harmonic oscillators because of the affine parametrization of the arc-length with time and the consequent fluctuations of the kinetic energy. By starting from the general solutions for a system of two harmonic oscillators given by

[TABLE]

where $\displaystyle A_{1},A_{2}$ and $\displaystyle\phi_{1},\phi_{2}$ are determined by the initial conditions, the authors of Ref.cuervo2015non used polar coordinates to rewrite the previous equations in the compact form

[TABLE]

where

[TABLE]

then they have reported that every initial condition fulfilling the condition

[TABLE]

yields unstable solutions. We have adopted the same initial conditions for two harmonic oscillators and representing the solutions of the JLC equations with respect to a parallel transported frame, already used in cerruti1996geometric , and we have found that the equations (64) with the frequencies (LABEL:frequency_jacobimetric_2dimension) display stable solutions. We have solved the equations by using two different conditions both fulfilling Eq. (68). These are

[TABLE]

Let us consider the condition (69). This is obtained by plugging through Eqs.(66) the initial condition

[TABLE]

into the frequency (LABEL:frequency_jacobimetric_2dimension). The time variation of $\displaystyle\tilde{\omega}_{2}(t)$ is reported in Figure 1.

Correspondingly, the time dependence of the Geometric Lyapunov Exponent, $\displaystyle\lambda(t)$, reported in Figure 2, clearly displays the typical time dependence found for regular trajectories, that is, $\displaystyle\lambda(t)\sim 1/t$.

Coming now to the second condition given in (70), again this is implemented by plugging through Eqs.(66) the initial condition

[TABLE]

into the frequency (LABEL:frequency_jacobimetric_2dimension).

The time variation of $\displaystyle\tilde{\omega}_{2}(t)$ is now reported in Figure 3. The corresponding time variation of the Geometric Lyapunov Exponent, $\displaystyle\lambda(t)$, is now reported in Figure 4. Again it is found that $\displaystyle\lambda(t)$ decays as $\displaystyle 1/t$, as expected for regular motions.

III.3 Parallel Transported Frame for a system of $\displaystyle 3$ harmonic oscillators

A relevant step forward is now obtained by considering three degrees of freedom because now the geodesic separation vector has two nontrivial components (since the component parallel to the velocity vector does not accelerate). Therefore, we consider $\displaystyle 3$ harmonic oscillators described by the Hamiltonian

[TABLE]

and the corresponding Jacobi metric, again with $\displaystyle W(q):=E-V(q_{1},q_{2},q_{3})$, is

[TABLE]

The eigenvectors of the Ricci tensor $\displaystyle Ric(\boldsymbol{\xi},\boldsymbol{\xi})$ are

[TABLE]

with the associated eigenvalues

[TABLE]

respectively.

With these eigenvectors the matrix for the basis transformation now is

[TABLE]

In order to make the parallel transported reference frame orthonormal, we have to orthogonalise the above eigenvectors with respect to the Jacobi metric, thus for $\displaystyle i=1,2,3$ we have

[TABLE]

namely

[TABLE]

These factors are

[TABLE]

The reference frame is composed by the normalized vectors

[TABLE]

and by representing the Jacobi vector field as $\displaystyle J=J^{i}E_{i}$, we write the Jacobi-Levi Civita equations for the three components as

[TABLE]

Finally, passing to the physical time and by using Eq. (53) we get

[TABLE]

with, again,

[TABLE]

The results reported in the following are worked out by numerically integrating the following equations

[TABLE]

with the following expressions for $\displaystyle\omega_{(2)}(t)$ and $\displaystyle\omega_{(3)}(t)$  :

[TABLE]

III.4 Numerical results

In Ref. cuervo2015non , the alleged definitive argument to rule out the use of Jacobi metric to consistently describe the stability/instability of Hamiltonian dynamics was given by considering $\displaystyle N$ decoupled harmonic oscillators. The claim was that non-vanishing fluctuations of kinetic energy (due to the non affine parametrization of the arc length with time) entail parametric resonance in the JLC equation mimicking chaos for an integrable system. The authors considered the solutions of this system in the form

[TABLE]

where $\displaystyle\theta_{k}=k\frac{2\pi f}{N}$, with the phases $\displaystyle\theta_{k}$ distributed on a fraction $\displaystyle f$ of the interval $\displaystyle 2\pi$. It has been reported that the smaller $\displaystyle f$ the larger fluctuation of kinetic energy and the larger the Lyapunov exponent. The fluctuation of kinetic energy $\displaystyle\sqrt{\sigma}$ is given in cuervo2015non by

[TABLE]

Notice that in the $\displaystyle N\to\infty$, the kinetic energy fluctuation magnitude $\displaystyle\sqrt{\sigma}$ has a non-vanishing value, so that the authors claim that any dimension this basic integrable system would display non-physical instabilities. In the next Section we will argue against this claim on the basis of an argument related with the concentration of measure at high dimension.

As shown in the preceding Section, and before in Refs.cerruti1996geometric ; rick , a consistent description of order and chaos is obtained using the Jacobi metric and writing the JLC equation for a parallel transported frame which is quite simple to be found for $\displaystyle N=2$, while the first non trivial extension is given by the $\displaystyle N=3$ case. The JLC equation for the three dimensional case is given by Eqs.(83). There are three principal directions of curvature, namely, the sectional curvatures [Eqs. (76)], identified by the planes generated by the velocity vector along a geodesic, $\displaystyle E_{1}=\boldsymbol{\xi}$, and the parallel transported basis vectors $\displaystyle E_{k}$ with $\displaystyle k=2,3$. These sectional curvatures coincide with the eigenvalues of the operator $\displaystyle Ric(\boldsymbol{\xi},\boldsymbol{\xi})=R(\cdot,\boldsymbol{\xi})\boldsymbol{\xi}$; one of these is obviously zero because $\displaystyle Ric(\boldsymbol{\xi},\boldsymbol{\xi})\boldsymbol{\xi}=0$ while the other two are given by $\displaystyle\Lambda_{k}=\boldsymbol{g}(E_{k},Ric(\boldsymbol{\xi},\boldsymbol{\xi})E_{k})$.

In Figure 6 of Ref.cuervo2015non , the largest value of $\displaystyle\lambda$ corresponds to a kinetic energy fluctuation level $\displaystyle\sqrt{\sigma}\simeq 0.7$ which can be obtained with different values of $\displaystyle N$ and $\displaystyle f$. The Geometrical Lyapunov exponent and the norm of the Jacobi vector field have been worked out by numerically integrating Eqs.(83) for the following cases: $\displaystyle N=3,\ f=0.05$ corresponding to $\displaystyle\sqrt{\sigma(0.05)}=0.6968$; for $\displaystyle N=3,\ f=0.1$ corresponding to $\displaystyle\sqrt{\sigma(0.05)}=0.6663$; and for $\displaystyle N=3,\ f=0.45$ corresponding to $\displaystyle\sqrt{\sigma(0.45)}=0.0900$. The outcomes are reported in Figures 6 and 6, in Figures 8 and 8, and in Figures 10 and 10, respectively.

Equations (83) have been integrated with a fourth-order Runge-Kutta algorithm along the $\displaystyle q_{k}(t)$ given by (85) and setting the phases $\displaystyle\theta_{k}$ uniformly distributed on a fraction $\displaystyle f$ of the interval $\displaystyle 2\pi$. It is well evident that the norm of the Jacobi geodesic separation vector is always bounded, coherently with the decay with $\displaystyle 1/t$ of the running value of $\displaystyle\lambda(t)$. The JLC equation written for the Jacobi metric and with a parallel transported frame provides the correct result: no instability of the trajectories of an integrable system is found, contrary to the claim of Ref.cuervo2015non .

IV Concentration of measure of the volume occupied by accessible configurations with Jacobi metric

We know from the work in Ref.cerruti1997lyapunov that the Jacobi-Levi Civita equation, written in the natural reference frame for a large number of degrees of freedom, appropriately works by producing Geometrical Lyapunov exponents in both qualitative and quantitative agreement with the standard Lyapunov exponents.

In view of the above discussed problems due to the bouncing of phase trajectories/geodesics on the Hill’s boundary, let us see why at large $\displaystyle N$ (large meaning just a few tens) the non physical divergencies are not found in spite of the use of the natural reference frame, in place of the parallel transported one, and in spite of the presence of kinetic energy fluctuations.

Consider a system composed by a large number of harmonic oscillators, denote by $\displaystyle\boldsymbol{Q}=\{q_{1},\cdots,q_{N}\}$ and $\displaystyle\boldsymbol{P}=\{p_{1},\cdots,p_{N}\}$ the conjugate momenta, with $\displaystyle\kappa=1$ for simplicity, the Hamiltonian is

[TABLE]

and the Jacobi metric in the Hill’s region $\displaystyle M_{E}$ is

[TABLE]

The associated Riemannian volume form is

[TABLE]

where the determinant of the metric is

[TABLE]

Therefore, the total volume of $\displaystyle M_{E}$ is given by the following integral

[TABLE]

It is worth noticing a remarkable property of the Jacobi metric, its associated volume is proportional (up to an $\displaystyle N$-dependent factor) to the microcanonical ensemble measure

[TABLE]

where $\displaystyle\mathcal{M}_{E}=\left\{(\boldsymbol{P},\boldsymbol{Q})\in T^{*}M_{E}\,|H(\boldsymbol{P},\boldsymbol{Q})=E\right\}$. Due to the quadratic form of the kinetic energy $\displaystyle 2W=R^{2}$, the integral (92) can be rewritten as

[TABLE]

where $\displaystyle\mathcal{A}(\mathbb{S}^{N-1})$ is the area of the unitary sphere. Let us now consider the volume of $\displaystyle M_{E}$ given by integral $\displaystyle\mathcal{V}_{N}(E)$ and show that it concentrates around an $\displaystyle N-1$-dimensional manifold $\displaystyle\Sigma_{\bar{V}}=\left\{\boldsymbol{Q}\in M_{E}\,|V(\boldsymbol{Q})=\bar{V}\right\}$, in other words, the overwhelming contribution to the volume integral is given by microscopic configurations far from the boundary of $\displaystyle M_{E}$. Although it could be questionable to provide a statistical argument for an integrable system for which ergodic hypothesis does not hold, the statistical averaging is intended as an averaging over all the possible configurations compatible with the constraint $\displaystyle H(\boldsymbol{P},\boldsymbol{Q})=E$. It is convenient to rewrite the integral (91) with the volume element expressed in spherical coordinates

[TABLE]

where $\displaystyle Q=\|\boldsymbol{Q}\|_{{\mathbb{R}}^{N}}=\sqrt{2V(Q)}$, because integrating over the angular variables one obtains

[TABLE]

where $\displaystyle x=V/E$ is the relative value of the potential energy with respect to the total energy and

[TABLE]

As we are interested in the limit of large $\displaystyle N$, we can apply the Laplace approximation to evaluate the previous integral, i.e. we consider the Taylor expansion around the minimum with respect to $\displaystyle x$ of $\displaystyle F$ in the interval $\displaystyle(0,1)$,

[TABLE]

where $\displaystyle\sigma^{2}=\left(F^{\prime\prime}(\bar{x})N\right)^{-1}>0$. For a generic value of $\displaystyle N$, the solution of $\displaystyle F^{\prime}(x)=0$ is

[TABLE]

which is actually a minimum since

[TABLE]

This means that the largest part of the volume $\displaystyle\mathcal{V}_{N}(E)$ is concentrated around the hypersurface at constant potential energy $\displaystyle\Sigma_{\bar{x}E}$, a result very close to what is expected from the virial theorem note .

From (99) and $\displaystyle\sigma^{2}=\left(F^{\prime\prime}(\bar{x})N\right)^{-1}$, it follows that the largest part of the volume ($\displaystyle\sim 99.7\%$) is concentrated around $\displaystyle\Sigma_{\bar{x}E}$ in the interval $\displaystyle[\bar{x}-3\sigma,\bar{x}+3\sigma]$ with

[TABLE]

This statistical argument shows that in the case of a large number of degrees of freedom the volume of the manifold is concentrated around a submanifold constant energy hypersurface $\displaystyle V=E/2$, far from the boundary of the Hill region where Jacobi metric is singular.

For example, with $\displaystyle N=100$ kinetic energy fluctuations of absolute value of $\displaystyle 10\%$ occur with a probability of $\displaystyle 68\%$ while the probability of configurations hitting the boundary ($\displaystyle E-V=0$) is $\displaystyle\sim 5.52\times 10^{-88}$.

V Discussion

Even though the point raised in Ref.cuervo2015non is interesting, the conclusion put forward by the authors is incorrect. The fluctuations of kinetic energy along a trajectory/geodesic of the Jacobi metric associated with an integrable system, like a collection of harmonic oscillators, are by no means responsible for the activation of parametric instability mimicking a chaotic behaviour. When the number of degrees of freedom of a Hamiltonian system is small, the associated geodesics can often approach the boundary $\displaystyle\partial M_{E}=\{q\in M|V(q)=E\}$ of the mechanical manifold $\displaystyle M_{E}$, and, in so doing, the geodesics bounce on $\displaystyle\partial M_{E}$. The sharp reflection of the geodesics on the so-called Hill’s boundaries montgomery2014s ; giambo2014morse ; giambo2015normal ; seifert1948periodische are at the origin of numerical instabilities which in principle could be perhaps avoided by a prohibitively high precision of the integration algorithm for the Jacobi–Levi-Civita equation describing the geodesic spread. However, throughout this paper we have shown that this problem can be fixed by choosing a parallel transported coordinate system. The stability/instability of geodesics is an intrinsic property thus in principle independent of the choice of the coordinate system, however, not all the coordinate systems are necessarily equivalent from the point of view of their numerical implementation and reliability of the corresponding outcomes. And, in fact, the sharp reflection of the geodesics by the boundaries $\displaystyle\partial M_{E}$ is accounted for by a sudden reflection of the coordinate axes of the parallel transported frames thus separating the true geometric origin of stability/instability of geodesics from the source of numerical artefacts related with their peculiar shape.

We have then shown that when the number of degrees of freedom increases, then the probability of approaching the boundary of the corresponding mechanical manifold $\displaystyle M_{E}$ gets lower and lower and, even if at finite $\displaystyle N$ the kinetic energy fluctuates it does not affect the strength of chaos measured through the outcomes of the JLC equation written for both the Jacobi and Einsenhart metrics which are in perfect agreement, as shown in Ref.cerruti1997lyapunov . Thus already for a few tens of degrees of freedom the JLC equation for $\displaystyle(M_{E},g_{J})$ written in natural chart cerruti1997lyapunov ; book

[TABLE]

can be safely used, at most with the exclusion of a zero measure set of initial conditions. For very weakly coupled harmonic oscillators and $\displaystyle N=128$ these equations give $\displaystyle\lambda$ as small as $\displaystyle 2\times 10^{-6}$ .

In conclusion, the study of order and chaos of Hamiltonian flows - identified as geodesic flows of the Jacobi metric in configuration space - is legitimate and coherent, although not unique.