Data-driven reconstruction of chaotic dynamical equations: the H\'enon-Heiles type system

In this study, the classical two-dimensional potential $V_N=\frac{1}{2}\,m\,\omega^2\,r^2 + \frac{1}{N}\,r^N\,\sin(N\,\theta)$, $N \in {\mathbb Z}^+$, is considered. At $N=1,2$, the system is superintegrable and integrable, respectively, whereas for $N>2$ it exhibits a richer chaotic dynamics. For instance, at $N=3$ it coincides with the H\'enon-Heiles system. The periodic, quasi-periodic and chaotic motions are systematically characterized employing time series, Poincar\'e sections, symmetry lines and the largest Lyapunov exponent as a function of the energy $E$ and the parameter $N$. Concrete results for the lowest cases $N=3,4$ are presented in complete detail. This model is used as a benchmark system to estimate the accuracy of the Sparse Identification of Nonlinear Dynamical Systems (SINDy) method, a data-driven algorithm which reconstructs the underlying governing dynamical equations. We pay special attention at the transition from regular motion to chaos and how this influences the precision of the algorithm. In particular, it is shown that SINDy is a robust and stable tool possessing the ability to generate non-trivial approximate analytical expressions for periodic trajectories as well.


I. INTRODUCTION
The interest on the Hénon-Heiles potential [1] was originally motivated by the astronomy community to investigate the 3D motion of a moving star in the gravitational field of an axis-symmetric galaxy.From the two evident integrals, namely the energy and the angular momentum around the symmetry axis, astronomers pondered the question on the existence of an additional conserved quantity I in the Liouville sense (the so called third integral ) that would make the system integrable.This searching, without any knowledge of the physical origin of I or about its possible mathematical expression.It turned out that the system is chaotic and no third integral exists.Accordingly, several analytical and numerical systematic studies on galactic potentials were conducted in the 1960s (see, for instance Contopoulos [2]) whilst recent references on the closely related classical and quantum chaos can be found in [3], [4].
Exploiting the conservation of angular momentum, Hénon and Heiles eventually obtained a reduced Hamiltonian with two degrees of freedom only.This reduced problem is called the Hénon-Heiles system (HHS).The corresponding classical Hamiltonian is of the form: where x is the altitude and y denotes the radius (the distance to the aforementioned symmetry axis), p x and p y are the corresponding canonical momenta and g > 0, the coupling constant, is a real parameter.The total energy * maqj@fata.unam.mxH = E of the above 2D system is the single conserved quantity of the problem.The potential of HHS can be viewed as a two-dimensional isotropic harmonic oscillator perturbed with nonlinear cubic terms.Also, it is worth mentioning that the Hamiltonian (1) possesses a dihedral symmetry, i.e., a D 3 symmetry.As indicated in [1], the Hamiltonian ( 1) is one of the simplest models to capture the relevant dynamics of non-integrable systems.In particular, it incorporates the physical feature that stars are trapped by an attractive potential, and those stars with sufficient energy may escape from the galaxy.
A straightforward numerical analysis [1] of HHS revealed the existence of a critical value E c of energy separating the bounded trajectories from the unbounded ones.For small values E ≪ E c , the cubic terms in the potential are negligible and the system displays a regular behaviour mostly.As the energy increases, these terms become relevant and a mixed combination of regular and chaotic motion occurs.Eventually, at E ∼ E c the chaotic regions dominate the phase space landscape.
Various generalizations of the HHS have been studied in the literature either by considering the coefficients of the cubic terms, in the potential, as free parameters or by adding higher order terms.There is a plethora of works about the integrability or non-integrability of Hénon-Heiles type systems with two parameters; see, for example, Llibre and Jiménez [5] and references included therein.Remarkably, the integrability of a generalized HHS that introduces two parameters in the quadratic terms of the potential was demonstrated by Grammaticos et al. [6], see also [7,8].
Needless to say that several applications of Hénon-Heiles type systems appear in different fields of physics.For example, an interesting theoretical study of a generalized 3D Hénon-Heiles potential was applied to model the dynamics of ions in a 3D axialy symmetric Penning trap [9], [10].
In the present study we analyze a one-parametric generalization of (1) with the following properties: (I) the potential is the sum of an isotropic 2D harmonic oscillator plus a polynomial function of degree N = 1, 2, . . ., in the variables (x, y), (II) the system admits a 2D discrete D N symmetry (group of symmetries of a regular polygon with N sides), (III) there exists a critical value of energy that separates the bounded and unbounded motion, and (IV) mixed regions in the phase space of regular and chaotic dynamics appears for N ≥ 3.
The goal of this work is two-fold.Firstly, we aim to investigate the rich dynamics of the system as a function of N and the energy E. To this end, we will employ time series, Poincaré sections, symmetry lines and the largest characteristic Lyapunov exponents.For instance, in the computation of the symmetry lines we exploit the discrete rotational symmetry of the Hamiltonian.Secondly, we will take this model as a benchmark system to test qualitative and quantitatively the algorithm of Machine Learning SINDy which reconstructs the governing dynamical equations of a given system using only the time series data.In particular, we are interested in the influence of the degree of chaoticity on the accuracy provided by SINDy.

II. THE MODEL
We consider a classical non-relativistic system with two degrees of freedom on the plane.The Hamiltonian is given by with scalar potential here r = x 2 + y 2 , θ = tan −1 (y/x) are polar coordinates, p r and p θ are their canonical conjugate momentum variables, respectively, m is the mass of the particle, ω > 0 and N > 0 is a positive integer number.The phase space is 4-dimensional.
The Hamiltonian (2) possesses a discrete rotational dihedral symmetry θ → θ + 2 π N , and for even N the reflection symmetry r → −r is present formally.Furthermore, the generalized Hénon-Heiles systems (2) are reversible.Since H N is a quadratic homogeneous function in the momenta, the Hamilton's equations remain invariant under the time and momenta inversions (r, θ, p r , p θ , t) → (r, θ, −p r , −p θ , −t) .
At N = 1, 2, the system is superintegrable and integrable, respectively, whereas for any N > 2 it exhibits a chaotic motion.In particular, at N = 3 with m = ω = 1 it coincides with the celebrated Hénon-Heiles potential (1) at g = 1.
In Cartesian coordinates, x = r cos θ and y = r sin θ, the Hamiltonian (2) reads where again p x and p y are the corresponding canonical momentum variables and P N (x, y) in ( 5) is an homogeneous polynomial of order N (N > 1) in coordinates x, y.Explicitly, for the lowest values of N the potential (6) The Hamiltonian (5) defines a fourth-order autonomous system with variables (x, y, p x , p y ).It is a dynamical system with polynomial equations of motion.
It is worth noting that under the rescaling r → α r accompanied by the transformations (7) the dynamics of the Hamiltonian (2) is scale-invariant.
For any N = 1, 2, 3, . . ., the time evolution of the variables x and y takes the form: thus, they are linear while for the momentum coordinates (p x , p y ) the dynamics, ṗx = − ∂ ∂ x H N and ṗy = − ∂ ∂ y H N , is N −dependent and non-linear starting from N = 3.
In the case N = 1, making the canonical change of variables x → x, y → y − 1 m ω 2 with momenta p x , p y unchanged, we arrive to the 2D isotropic harmonic oscillator.Thus, the Hamiltonian H N =1 , the angular momentum L z = x p y − y p x as well as S xy = p x p y + m 2 ω 2 x y are three algebraically independent integrals of motion.The system H N =1 and its quantum counterpart are maximally superintegrable [11,12].Similarly, it can be shown that the case N = 2 corresponds to an integrable system.For higher values of N , the equations of motion become nonlinear.In particular, at N = 3, we obtain: For N = 4 we have the equations: At N = 5: Hereafter, we put m = 1 and ω = 1.For fixed N , there exists a critical value (escape energy) above which the energy surface H N = E is unbounded, i.e., the trajectories of the system transit from bounded to unbounded motion at E = E c .This critical energy is a monotonously increasing function of N .It runs from − 1 2 at N = 1 to 1 2 for N → ∞.In Figure 1, for the cases N = 3, 4, 5, 6 the level curves, on the plane (x, y), of the potential V N (x, y) (3) are depicted.The dihedral D N symmetry dictates the regular form of these curves.

III. BOUNDED MOTION
In this section, for the cases N = 3 and N = 4 we explore the bounded dynamics of the Hamiltonian H N (5) as a function of the energy E.

A. Dynamics in configuration space
In the case of finite motion, the boundaries in configuration space are defined by the equation For the values N = 3 (top panel of Fig. 2) and N = 4 (bottom panel of Fig. 2), these boundaries exhibit the 2D dihedral D N symmetry, as expected.We also present in Fig. 2 examples of generic periodic, quasi-periodic and chaotic trajectories occurring for 0 < E < E c .These trajectories were calculated numerically.

B. Time series
Now, for the cases N = 3 and N = 4, we display time series (trajectories) obtained by solving numerically the Hamilton's equations of motion.In particular, we plot the data of the dynamical conjugate variables y(t) and p y (t).For each value of energy E = 3 4 E c , 99 100 E c , we present three representative trajectories of the regular (periodic), quasi-periodic and chaotic motion, whereas for E = 1 4 E c only periodic and quasi-periodic ones are considered (at this small value of energy there is no significant presence of chaos).The corresponding results for N = 3 and N = 4 are shown in Fig. 3 and 4, respectively.A straightforward Fourier decomposition of the aforementioned time series confirms the nature of these trajectories.For completeness, in Fig. 5 we show the associated dynamical motion in the phase space plane (y, p y ).

IV. POINCAR É SECTIONS AND LYAPUNOV EXPONENTS
In this section, we present the oriented Poincaré sections on the (y, p y ) plane, considering m = 1 and ω = 1, for the cases N = 3, 4 as a function of energy E. These sections are obtained from the intersection of orbits, associated with specific initial conditions within the state space, with a lower-dimensional subspace (x = 0, y, p x , p y ).These sections are transversal planes to the flow of the system.Interestingly, the Poincaré section can be seen as a discretized version of the dynamical system, For higher values of N , the equations of motion become nonlinear.In particular, at N = 3, we obtain: For N = 4 we have the equations: At N = 5: Hereafter, we put m = 1 and ω = 1.For fixed N , there exists a critical value (escape energy) above which the energy surface H N = E is unbounded, i.e., the trajectories of the system transit from bounded to unbounded motion at E = E c .This critical energy is a monotonously increasing function of N .It runs from In Figure 1, for the cases N = 3, 4, 5, 6 the level curves, on the plane (x, y), of the potential V N (x, y) (3) are depicted.The dihedral D N symmetry dictates the regular form of these curves.

III. BOUNDED MOTION
In this section, for the cases N = 3 and N = 4 we explore the bounded dynamics of the Hamiltonian H N (5) as a function of the energy E.

A. Dynamics in configuration space
In the case of finite motion, the boundaries in configuration space are defined by the equation For the values N = 3 (top panel of Fig. 2) and N = 4 (bottom panel of Fig. 2), these boundaries exhibit the 2D dihedral D N symmetry, as expected.We also present in Fig. 2 examples of generic periodic, quasi-periodic and chaotic trajectories occurring for 0 < E < E c .These trajectories were calculated numerically.

B. Time series
Now, for the cases N = 3 and N = 4, we display time series (trajectories) obtained by solving numerically the Hamilton's equations of motion.In particular, we plot the data of the dynamical conjugate variables y(t) and p y (t).For each value of energy E = 3 4 E c , 99 100 E c , we present three representative trajectories of the regular (periodic), quasi-periodic and chaotic motion, whereas for E = 1  4 E c only periodic and quasi-periodic ones are considered (at this small value of energy there is no significant presence of chaos).The corresponding results for N = 3 and N = 4 are shown in Fig. 3 and 4, respectively.A straightforward Fourier decomposition of the aforementioned time series confirms the nature of these trajectories.For completeness, in Fig. 5 we show the associated dynamical motion in the phase space plane (y, p y ).

IV. POINCAR É SECTIONS AND LYAPUNOV EXPONENTS
In this section, we present the oriented Poincaré sections on the (y, p y ) plane, considering m = 1 and ω = 1, for the cases N = 3, 4 as a function of energy E. These sections are obtained from the intersection of orbits, associated with specific initial conditions within the state space, with a lower-dimensional subspace (x = 0, y, p x , p y ).These sections are transversal planes to the flow of the system.Interestingly, the Poincaré section can be seen as a discretized version of the dynamical system,  retaining many characteristics of the original continuous system but operating in a reduced state space [13][14][15].Complementarily, we analyze the trajectories' nature by calculating the largest Lyapunov exponents.The Lyapunov exponent provides a measure of the average exponential divergence of states that began infinitesimally close.By observing the evolution of nearby trajectories over time, the Lyapunov exponent offers valuable insights into the system's sensitivity to initial conditions, enabling the characterization of chaotic behavior in dynamical systems [16,17].The largest Lyapunov exponents can be obtained by computing the double limit: where δ ⃗ X(t) = Φ • δ ⃗ X(t 0 ) with ⃗ X = [x, y, p x , p y ] T and δ ⃗ X(t 0 ) the perturbed initial state.Φ is the fundamental matrix whose elements are calculated by solving Φ = JΦ and Hamilton equations, simultaneously.Here, J denotes the Jacobian matrix.
Figure 6 and 7 show oriented Poincaré sections and largest Lyapunov exponents for the Hamiltonians H N =3 and H N =4 , respectively, considering three different energies: To draw the Poincaré sections, we develop numerical simulations for 70 different (random) initial conditions with a simulation time of 7000 seconds.In contrast, for each value of the energy, Lyapunov exponents plots were constructed by integrating numerically 15000 random initial conditions during 7000 seconds.Computations were performed in Mathematica 12.1.
In Figure 6(a),(d), the Poincaré sections present regular dynamics for low energy, E = 1 4 E c .This is confirmed by Lyapunov exponents extremely close to zero.However, for the cases E = 3 4 E c (Fig. 6(b),(e)), Poincare sections show regions with scattered points whose Lyapunov exponents plot exhibit values different from zero, unveiling the presence of chaos in the system.Importantly, some islands of periodic behaviors are clearly visible which are denoted by dark regions.Despite the pre- retaining many characteristics of the original continuous system but operating in a reduced state space [13][14][15].Complementarily, we analyze the trajectories' nature by calculating the largest Lyapunov exponents.The Lyapunov exponent provides a measure of the average exponential divergence of states that began infinitesimally close.By observing the evolution of nearby trajectories over time, the Lyapunov exponent offers valuable insights into the system's sensitivity to initial conditions, enabling the characterization of chaotic behavior in dynamical systems [16,17].The largest Lyapunov exponents can be obtained by computing the double limit: where δ ⃗ X(t) = Φ • δ ⃗ X(t 0 ) with ⃗ X = [x, y, p x , p y ] T and δ ⃗ X(t 0 ) the perturbed initial state.Φ is the fundamental matrix whose elements are calculated by solving Φ = JΦ and Hamilton equations, simultaneously.Here, J denotes the Jacobian matrix.Here , whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x 0 , p y 0 ) are those used for the bottom panel of Fig. 2.

A. Symmetry lines
Additionally, we show the symmetry lines (defined in [18]) for the cases of N = 3 and N = 4 .These symmetry lines are used for obtaining the initial conditions of periodic orbits.
Moreover, for even values of N , they also possess the symmetry of coordinates and time inversion: For these reversible systems, the Poincaré map T can be expressed as the composition of two involutions [18]: where T, I 0 , I 1 : Σ → Σ, here Σ denotes a Poincaré section.Introducing I j = T j I 0 , the infinite set of transformations T j , I i : Σ → Σ, with i, j ∈ Z, forms an infinite discrete group.The fixed points of I i are referred to as the symmetry lines: Of particular significance are the fundamental lines Γ 0 and Γ 1 , they lead to the property This implies that Γ 0 generates all symmetry lines with values of the energy.Here E = 1 4 E c in (a), E = 3 4 E c in (b), whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x 0 , p y 0 ) are those used for the top panel in Fig. 2.Here , whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x 0 , p y 0 ) are those used for the bottom panel of Fig. 2.
, whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x0 , p y0 ) are those used for the bottom panel of Fig. 2. even indices: Γ 2i = T i Γ 0 whilst Γ 1 produces all symmetry lines with odd indices: Γ 2i+1 = T i Γ 1 .It can be proven that the intersections of the symmetry lines Γ i ∩Γ j correspond to periodic orbits, with the period being a divisor of the integer |i − j|.
For the Henon-Heiles system (2) with N = 3, the involution I 0 is used in the Σ plane, along with its fixed points: For the case of N = 4, the relevant symmetry is as follows: These symmetry transformations are employed to obtain initial conditions for periodic orbits.

V. RECONSTRUCTION OF THE HAMILTON'S EQUATIONS FROM TIME SERIES
This section explores the process of reconstructing the Hamilton's equations based only on the time series data computed in Section III B. First of all, we briefly describe the key elements of the sparse identification of nonlinear dynamical systems (SINDy) algorithm [19][20][21][22][23].
The goal is to estimate quantitatively the accuracy of this data-driven algorithm.

A. SINDy Algorithm
Let X ∈ R n×4 be a matrix in which each column contains n samples of the time series for each of the four dynamical variables x, y, p x and p y .Moreover, Ẋ ∈ R n×4 denotes a matrix whose columns are numerical derivatives of each time series data.With this in mind, we write the following matrix equation Here, Θ(x, y, p x , p y ) ∈ R 1×m is the vector of candidate functions that could potentially constitute the recon-

A. Symmetry lines
Additionally, we show the symmetry lines (defined in [18]) for the cases of N = 3 and N = 4 .These symmetry lines are used for obtaining the initial conditions of periodic orbits.
Moreover, for even values of N , they also possess the symmetry of coordinates and time inversion: (q, p, t) → (−q, p, −t) .
For these reversible systems, the Poincaré map T can be expressed as the composition of two involutions [18]: where T, I 0 , I 1 : Σ → Σ, here Σ denotes a Poincaré section.Introducing I j = T j I 0 , the infinite set of transformations T j , I i : Σ → Σ, with i, j ∈ Z, forms an infinite discrete group.The fixed points of I i are referred to as the symmetry lines: Of particular significance are the fundamental lines Γ 0 and Γ 1 , they lead to the property This implies that Γ 0 generates all symmetry lines with even indices: Γ 2i = T i Γ 0 whilst Γ 1 produces all sym- structed Hamilton's equations and Θ(X) ∈ R n×m represents the evaluation of the vector of candidate functions for the time series X.On the other hand, the entries of Ξ ∈ R m×4 will be the coefficients of such candidate functions and therefore the unknowns.Since the matrix equation ( 22) is overdetermined, we aim to find a solution such that Ξ minimizes || Ẋ − Θ(X) Ξ|| where || • || denotes the Euclidean norm.
It is important to note that this least squares solution may not necessarily be sparse, meaning that the majority of its elements are not zero.The work in [19] proposed SINDy which is a method to find a sparse solution for the system identification problem, whose results suggest that, with an appropriate choice of candidate functions, a least squares sparse solution will result in a parsimonious model that aligns with the governing equations.Finally, given a sparsity solution obtained via SINDy, the reconstructed model is presented as follows: Our work uses the Python library PySINDy [24,25] to reconstruct the Hamilton's equations from the time series data computed in Section III B. For a more comprehensive understanding of the SINDy algorithm, we direct the reader to [19].

B. Results and Discussion
As previously mentioned in subsection V A, concerning the reconstruction of state equations from the time series, it is necessary to choose candidate functions that will comprise the reconstructed model.A physical and mathematically based selection is fundamental to obtain optimal results.Since in this case we know the exact model ( 5), for each state equation ẋ, ẏ, ṗx , ṗy , a suitable option is to take a polynomial function of degree K > N − 1 in the variables x, y, p x , p y whose C-coefficients will be the decision variables (the unknowns).For example, the most general polynomial ansatz of order K is of the form being the C-constants real parameters.Similar expressions for the remaining three state equations are employed.
In order to determine the accuracy of SINDy, let us 7 metry lines with odd indices: Γ 2i+1 = T i Γ 1 .It can be proven that the intersections of the symmetry lines Γ i ∩Γ j correspond to periodic orbits, with the period being a divisor of the integer |i − j|.
For the Henon-Heiles system (2) with N = 3, the involution I 0 is used in the Σ plane, along with its fixed points: For the case of N = 4, the relevant symmetry is as follows: These symmetry transformations are employed to ob-tain initial conditions for periodic orbits.

V. RECONSTRUCTION OF THE HAMILTON'S EQUATIONS FROM TIME SERIES
This section explores the process of reconstructing the Hamilton's equations based only on the time series data computed in Section III B. First of all, we briefly describe the key elements of the sparse identification of nonlinear dynamical systems (SINDy) algorithm [19][20][21][22][23]. .The (marked) bigger points in green, black and red indicate initial conditions for periodic, quasi-periodic and chaotic trajectories, respectively; they were used to construct Figs. 2-5 .
introduce the error parameter ∆C τ (σ) as follows here, C exact τ (σ) represents the coefficient of the monomial σ in the state equation τ .Note that the constant (25), which measures the precision of SINDy, depends on the value dt of the subinterval of time used in the numerical integration of the equations of motion.In general, for sufficiently small values of dt and taking solely the set of data corresponding to a single trajectory, the algorithm of SINDy reproduces exactly the same functional form of the original Hamilton's equations.Naturally, as we increase the value of dt, extra terms in the reconstructed equations occur.We define a critical time step dt c as the smallest value of dt for which an extra term appears in any of the state equations.Summarizing, N is fixed and for each value of energy E =1 4 E c , 3 4 E c , 99 100 E c , we select a trajectory to generate the data (a numerical time series) that (together with polynomial candidate functions) will be used as an input to reconstruct via SINDy the (state) Hamilton's equations.Afterwards, we calculate the corresponding value of ∆C τ (σ) as a function of the step dt.This process is repeated in the same manner for a periodic, quasi-8 The goal is to estimate quantitatively the accuracy of this data-driven algorithm.

A. SINDy Algorithm
Let X ∈ R n×4 be a matrix in which each column contains n samples of the time series for each of the four dynamical variables x, y, p x and p y .Moreover, Ẋ ∈ R n×4 denotes a matrix whose columns are numerical derivatives of each time series data.With this in mind, we write the following matrix equation Here, Θ(x, y, p x , p y ) ∈ R 1×m is the vector of candidate functions that could potentially constitute the reconstructed Hamilton's equations and Θ(X) ∈ R n×m repre-sents the evaluation of the vector of candidate functions for the time series X.On the other hand, the entries of Ξ ∈ R m×4 will be the coefficients of such candidate functions and therefore the unknowns.Since the matrix equation ( 22) is overdetermined, we aim to find a solution such that Ξ minimizes || Ẋ − Θ(X) Ξ|| where || • || denotes the Euclidean norm.
It is important to note that this least squares solution may not necessarily be sparse, meaning that the majority of its elements are not zero.The work in [19] proposed SINDy which is a method to find a sparse solution for the system identification problem, whose results suggest that, with an appropriate choice of candidate functions, a least squares sparse solution will result in a parsimonious model that aligns with the governing equations.Finally, given a sparsity solution obtained via SINDy, the recon- Concrete results, taking K = N and the threshold Γ = 0.05, for the cases N = 3 and N = 4 are presented in Figs. 9, 11 and Figs. 12, 13, respectively.In general, most of the curves ∆C τ (σ) vs dt are smooth increasing positive functions vanishing as dt → 0. This implies that, regardless of the nature of motion (periodic or chaotic), the SINDy algorithm recovers the correct underlying dynamical equations of the system where a single trajectory was employed in the reconstruction.
Essentially, at fixed (N, E, dt), the parameter ∆C τ (σ) takes its minimum value for a periodic trajectory whereas for a chaotic one it reaches its maximum.Nevertheless, if dt is sufficiently small ≈ 10 −4 , the presence of chaos does not influence the accuracy provided by SINDy.
given by: ẋSINDy = p x ẏSINDy = p y ṗxSINDy = 0.00322 y − 1.00132 x − 0.75788 x y 2 + 2.81860 y p x p y − 1.88935 y x 2 − 1.19774 y p 2 x − 1.16865 x p x p y − 0.25909 x 3 + 0.49646 x p 2 x ṗySINDy = − y + 2.48528 y x 2 − 0.00042 x p x p y − 1.51478 x 3 + 0.00012 x p 2 x , (26) here the coefficients were computed with an accuracy of ≈ 3−4 significant digits.Since we know the exact model, the rhs of the above equations can be decomposed into the sum of the exact part plus extra terms, namely Clearly, no extra terms occur for ẋSINDy and ẏSINDy .However, in the case of ṗxSINDy and ṗySINDy they read given by: ẋSINDy = p x ẏSINDy = p y ṗx SINDy = 0.00322 y − 1.00132 x − 0.75788 x y 2 + 2.81860 y p x p y − 1.88935 y x 2 − 1.19774 y p 2 x − 1.16865 x p x p y − 0.25909 x 3 + 0.49646 x p 2 x ṗy SINDy = − y + 2.48528 y x 2 − 0.00042 x p x p y − 1.51478 x 3 + 0.00012 x p 2 x , (26) here the coefficients were computed with an accuracy of ≈ 3−4 significant digits.Since we know the exact model, the rhs of the above equations can be decomposed into the sum of the exact part plus extra terms, namely Clearly, no extra terms occur for ẋSINDy and ẏSINDy .However, in the case of ṗx SINDy and ṗy SINDy they read S ṗx ,extra = 0.00322 y − 0.00132 x − 0.75788 x y 2 − y 3   27), the exact and extra parts, obtained by SINDy.These terms were evaluated over the periodic trajectory displayed in Fig. 2 (d).The variation of the extra terms is of order ≤ 10 −7 .Now, evaluated over the original numerical trajectory with initial conditions of Fig. 2 (d), the contribution of S ṗx ,extra and S ṗy ,extra above is effectively zero!, see Fig. 8.In particular, from the second equation in (28) neglecting the terms with coefficients of order ≈ 10 −4 , we arrive to the approximate analytical relation y = 0.41 x , which corresponds with the periodic trajectory in configuration space displayed in Fig. 2 (d).Hence, for arbitrary small values of dt, SINDy can generate spurious terms in the reconstruction of the governing equations of motion.Relative Error C ( )   Relative Error C ( )   Relative Error C ( )     Relative Error C ( )     For instance, the periodic trajectory y = y(x) occurring at N = 3 in Fig. 2 (d) does not admits an accurate representation in the polynomial basis (candidate functions) employed by SINDy.Thus, in that case the reconstruction of the dynamical equations is indeed correct (no extra terms appear).
Remark: For a given Hamiltonian system, the SINDy algorithm together with an educated election of the candidate functions (a basis) can be used to find approximate analytical expressions of periodic trajectories.Remark: For a given Hamiltonian system, the SINDy algorithm together with an educated election of the candidate functions (a basis) can be used to find approximate analytical expressions of periodic trajectories.For instance, the periodic trajectory y = y(x) occurring at N = 3 in Fig. 2 (d) does not admits an accurate representation in the polynomial basis (candidate functions) employed by SINDy.Thus, in that case the reconstruction of the dynamical equations is indeed correct (no extra terms appear).
Remark: For a given Hamiltonian system, the SINDy algorithm together with an educated election of the candidate functions (a basis) can be used to find approximate analytical expressions of periodic trajectories.

VI. CONCLUSIONS
In this study, we defined a Hénon-Heiles like potential V N = 1 2 m ω 2 r 2 + 1 N r N sin(N θ), N ∈ Z + .For arbitrary N > 2, in Cartesian coordinates (x, y), it corresponds to a polynomial function of degree N possessing a D N dihedral symmetry.For E < E c , the motion is finite and the system displays a chaotic behavior.A systematic characterization of the dynamics of the system as a function of N and the energy E was carried out.In particular, the explicit calculation of symmetry lines allowed us to compute periodic trajectories with high accuracy.In the second part of the paper, this model was employed to test SINDy which is a method to find the underlying governing (Hamilton's) equations using only the time series data of the trajectories obtained numerically.In general, for sufficiently small time-step of the data, SINDy provides convergent results (even for chaotic trajectories) to the exact functional form of the dynamical equations.Otherwise, the chaotic behaviour decreases the accuracy of the method.Interestingly, it was pointed out that when the corresponding series time data admits an accurate finite expansion in the corresponding basis of candidate functions, then this data-driven algorithm cannot recover the equations of motion correctly.Notably, for a given Hamiltonian, that is to say, when a reconstruction of the governing equations is not needed, SINDy can also be exploited to obtain approximate analytical expressions for the periodic trajectories which is, in general, a non-trivial task.For future work, among the interesting open questions we can mention the quantum counterpart of the model as well as the extension to the case with rational N .We plan to study the use of the SINDy algorithm and other methods alike to find symmetries (integrals of motion) of generic Hamiltonian systems as well.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Figure 1 :
Figure 1: Level curves, on the plane (x, y), of the potential V N (3) as a function of N .The corresponding values of the critical energy E crit = N −2 2 N are displayed as well.The value of the combination m ω 2 = 1 is used in the calculations.

Figure 1 :
Figure 1: Level curves, on the plane (x, y), of the potential V N (3) as a function of N .The corresponding values of the critical energy E crit = N −2 2 N are displayed as well.The value of the combination m ω 2 = 1 is used in the calculations.

Figure 6
Figure 6 and 7 show oriented Poincaré sections and largest Lyapunov exponents for the Hamiltonians H N =3 and H N =4 , respectively, considering three different energies: E = 1 4 E c , E = 3 4 E c and E = 99 100 E c .To draw the Poincaré sections, we develop numerical simulations for 70 different (random) initial conditions with a simulation time of 7000 seconds.In contrast, for each value of the energy, Lyapunov exponents plots were constructed by integrating numerically 15000 random initial conditions during 7000 seconds.Computations were performed in Mathematica 12.1.InFigure6(a),(d), the Poincaré sections present regular dynamics for low energy, E = 1 4 E c .This is confirmed by Lyapunov exponents extremely close to zero.However, for the cases E =3  4 E c (Fig.6(b),(e)), Poincare sections show regions with scattered points whose Lyapunov exponents plot exhibit values different from zero, unveiling the presence of chaos in the system.Impor-

Figure 3 :Figure 4 :
Figure 3: Case N = 3: time series of the dynamical variables y(t) and p y (t) of H N =3 (5) for three representative trajectories of the periodic (green dashed), quasi-periodic (black solid) and chaotic (red dotted) motion at different values of the energy.Here E = 1 4 E c in (a), E = 3 4 E c in (b), whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x 0 , p y 0 ) are those used for the top panel in Fig.2.

Figure 3 :
Figure 3: Case N = 3: time series of the dynamical variables y(t) and p y (t) of H N =3 (5) for three representative trajectories of the periodic (green dashed), quasi-periodic (black solid) and chaotic (red dotted) motion at different values of the energy.Here E = 1 4 E c in (a), E = 3 4 E c in (b), whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x0 , p y0 ) are those used for the top panel in Fig. 2.

Figure 4 :
Figure 4: Case N = 4: time series y(t) and p y (t) of H N =4 (5) for three representative trajectories of the periodic (green dashed), quasi-periodic (black solid) and chaotic (red dotted) motion at different values of the energy.HereE = 1 4 E c in (a), E = 3 4 E c in (b), whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x 0 , p y 0 ) are those used for the bottom panel of Fig.2.

Figure 4 :
Figure 4: Case N = 4: time series y(t) and p y (t) of H N =4 (5) for three representative trajectories of the periodic (green dashed), quasi-periodic (black solid) and chaotic (red dotted) motion at different values of the energy.Here E = 1 4 E c in (a), E = 3 4 E c in (b), whereas E = 99 100 E c for (c).The corresponding initial conditions (x 0 , y 0 , p x0 , p y0 ) are those used for the bottom panel of Fig. 2.

Figure 5 :
Figure 5: Cases N = 3 (top) and N = 4 (bottom): numerical trajectories p y = p y (y), corresponding to those displayed in Fig. 2, for representative orbits of the periodic (green dashed), quasi-periodic (black solid) and chaotic (red dotted) motion.Here E = 1 4 E c in (a) and (d), E = 3 4 E c in (b) and (e), while E = 99 100 E c for (c) and (f).

Figure 5 :
Figure 5: Cases N = 3 (top) and N = 4 (bottom): numerical trajectories p y = p y (y), corresponding to those displayed in Fig. 2, for representative orbits of the periodic (green dashed), quasi-periodic (black solid) and chaotic (red dotted) motion.Here E = 1 4 E c in (a) and (d), E = 3 4 E c in (b) and (e), while E = 99 100 E c for (c) and (f).

Figure 6 :
Figure 6: Case N = 3: (a)-(c) oriented Poincaré sections and symmetry lines Γ 0 and Γ 2 (16), and (d)-(f) the largest Lyapunov exponent on the plane (y, p y ), for the Hamiltonian H N =3 (5) with m = 1, ω = 1 at different values of the energy E = 1 4 E c , 3 4 E c , 99 100 E c .In this case, the critical energy takes the value E c = 1 6 .The (marked) bigger points in green, black and red indicate initial conditions for periodic, quasi-periodic and chaotic trajectories, respectively; they were used to construct Figs. 2-5 .

Figure 6 :
Figure 6: Case N = 3: (a)-(c) oriented Poincaré sections and symmetry lines Γ 0 and Γ 2 (16), and (d)-(f) the largest Lyapunov exponent on the plane (y, p y ), for the Hamiltonian H N =3 (5) with m = 1, ω = 1 at different values of the energy E = 1 4 E c , 3 4 E c , 99 100 E c .In this case, the critical energy takes the value E c = 1 6 .The (marked) bigger points in green, black and red indicate initial conditions for periodic, quasi-periodic and chaotic trajectories, respectively; they were used to construct Figs. 2-5 .

For
the sake of clarifying our notation, let us examine the scenario where N = 3.The exact state equation for the derivative of p x with respect to time is ṗx = −x − 2 y x.Thus, τ = ṗx , σ = x, x y, and the non-zero coefficients are C exact ṗx (x) = −1 and C exact ṗx (x y) = −2 only.Similarly, C SINDy τ (σ) will denote the coefficient of the monomial σ in the (approximate) state equation for τ constructed by the SINDy algorithm.

Figure 7 :
Figure 7: Case N = 4: (a)-(c) oriented Poincaré sections and symmetry lines Γ 0 and Γ 2 (16), and (d)-(f) largest Lyapunov exponents on the plane (y, p y ), for the Hamiltonian H N =4 (5) with m = 1, ω = 1 at different values of the energy E = 1 4 E c , 3 4 E c , 99 100 E c .In this case, the critical energy is E c = 1 4 , see text.The (marked) bigger points in green, black and red indicate initial conditions for periodic, quasi-periodic and chaotic trajectories, respectively; they were used in Figs.2-5

Figure 7 :
Figure 7: Case N = 4: (a)-(c) oriented Poincaré sections and symmetry lines Γ 0 and Γ 2 (16), and (d)-(f) largest Lyapunov exponents on the plane (y, p y ), for the Hamiltonian H N =4 (5) with m = 1, ω = 1 at different values of the energy E = 1 4 E c , 3 4 E c , 99 100 E c .In this case, the critical energy is E c = 1 4 , see text.The (marked) bigger points in green, black and red indicate initial conditions for periodic, quasi-periodic and chaotic trajectories, respectively; they were used in Figs.2-5

Figure 8 :
Figure 8: Case N = 4: the time series of the two terms in (27), the exact and extra parts, obtained by SINDy.These terms were evaluated over the periodic trajectory displayed in Fig. 2 (d).The variation of the extra terms is of order ≤ 10 −7 .

3 C
x (p x ) C y (p y ) C px (x) C px (xy) C py (y) C py (x 2 ) C py (y 2 )

1e 5 CFigure 9 :
Figure 9: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory and (b) quasi-periodic case.The critical time integration increments (see text) are dt c = 0.01763 and dt c = 0.02628, respectively.These trajectories correspond to the energy value E = 1 4 E c .

Figure 8 : 2 + 2 . 2 x− 1 . 3 + 2 . 2 x− 1 .
Figure 8: Case N = 4: the time series of the two terms in (27), the exact and extra parts, obtained by SINDy.These terms were evaluated over the periodic trajectory displayed in Fig. 2 (d).The variation of the extra terms is of order ≤ 10 −7 .Now, evaluated over the original numerical trajectory with initial conditions of Fig. 2 (d), the contribution of S ṗx,extra and S ṗy,extra above is effectively zero!, see Fig. 8.In particular, from the second equation in (28) neglecting

Figure 8 : 3 C
Figure 8: Case N = 4: the time series of the two terms in (27), the exact and extra parts, obtained by SINDy.These terms were evaluated over the periodic trajectory displayed in Fig. 2 (d).The variation of the extra terms is of order ≤ 10 −7 .Now, evaluated over the original numerical trajectory with initial conditions of Fig. 2 (d), the contribution of S ṗx,extra and S ṗy,extra above is effectively zero!, see Fig. 8.In particular, from the second equation in (28) neglecting the terms with coefficients of order ≈ 10 −4 , we arrive to

Figure 9 :
Figure 9: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory and (b) quasi-periodic case.The critical time integration increments (see text) are dt c = 0.01763 and dt c = 0.02628, respectively.These trajectories correspond to the energy value E = 1 4 E c .

Figure 9 : 3 C
Figure 9: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory and (b) quasi-periodic case.The critical time integration increments (see text) are dt c = 0.01763 and dt c = 0.02628, respectively.These trajectories correspond to the energy value E = 1 4 E c .

Figure 10 :
Figure 10: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.The critical time increments are dt c = 0.01925, dt c = 0.02726 and dt c = 0.13303, respectively.These trajectories correspond to the same energy value E = 3 4 E c .

Figure 11 :
Figure 11: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.The critical time increments are dt c = 0.02120, dt c = 0.02406 and dt c = 0.08443, respectively.These trajectories correspond to the same energy value E = 99 100 E c .

Figure 10 : 3 C
Figure 10: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.The critical time increments are dt c = 0.01925, dt c = 0.02726 and dt c = 0.13303, respectively.These trajectories correspond to the same energy value E = 3 4 E c .

Figure 10 :
Figure 10: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.The critical time increments are dt c = 0.01925, dt c = 0.02726 and dt c = 0.13303, respectively.These trajectories correspond to the same energy value E = 3 4 E c .

Figure 11 :
Figure 11: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.The critical time increments are dt c = 0.02120, dt c = 0.02406 and dt c = 0.08443, respectively.These trajectories correspond to the same energy value E = 99 100 E c .

Figure 11 :
Figure 11: Case N = 3.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.The critical time increments are dt c = 0.02120, dt c = 0.02406 and dt c = 0.08443, respectively.These trajectories correspond to the same energy value E = 99 100 E c .

3. N = 4 12
For instance, the periodic trajectory y = y(x) occurring at N = 3 in Fig.2 (d) does not admits an accurate representation in the polynomial basis (candidate functions) employed by SINDy.Thus, in that case the reconstruction of the dynamical equations is indeed correct (no extra terms appear).

Figure 12 :
Figure 12: Case N = 4.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory and (b) quasi-periodic case.These trajectories correspond to the energy value E = 1 4 E c .The critical time increment for the periodic trajectory is zero, but for the quasi-periodic case is given by dt c = 0.0305.

Figure 13 :
Figure 13: Case N = 4.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.For (a) the critical time increment is zero whereas for (b) and (c), they are dt c = 0.0136 and dt c = 0.01005, respectively.These three trajectories correspond to the energy value E = 3 4 E c .

Figure 12 :
Figure 12: Case N = 4.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory and (b) quasi-periodic case.These trajectories correspond to the energy value E = 1 4 E c .The critical time increment for the periodic trajectory is zero, but for the quasi-periodic case is given by dt c = 0.0305.

Figure 12 :
Figure 12: Case N = 4.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory and (b) quasi-periodic case.These trajectories correspond to the energy value E = 1 4 E c .The critical time increment for the periodic trajectory is zero, but for the quasi-periodic case is given by dt c = 0.0305.

Figure 13 :
Figure 13: Case N = 4.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.For (a) the critical time increment is zero whereas for (b) and (c), they are dt c = 0.0136 and dt c = 0.01005, respectively.These three trajectories correspond to the energy value E = 3 4 E c .

Figure 13 :
Figure 13: Case N = 4.The relative error ∆C τ (σ), see (25), of the coefficients obtained by SINDy for the reconstructed Hamilton's equations; (a) periodic trajectory, (b) quasi-periodic and (c) chaotic motion.For (a) the critical time increment is zero whereas for (b) and (c), they are dt c = 0.0136 and dt c = 0.01005, respectively.These three trajectories correspond to the energy value E = 3 4 E c .