Comment on “Hyperchaos in constrained Hamiltonian system and its control” by J. Li, H. Wu and F. Mei

The aim of this comment is to show that discovery of hyperchaos in three systems investigated in Li et al. (Nonlinear Dyn 94(3):1703–1720, 2018) is not correct. It is justified both theoretically and numerically. Corrected calculations of Lyapunov exponents and corresponding bifurcation diagram are given. Examples of hyperchaotic Hamiltonian multiple pendulum systems are presented.

hyperchaotic systems in various fields of science. However, it seems that the hyperchaos in Hamiltonian systems still does not attract sufficient attention. The aim of the authors in [1] was to propose a Hamiltonian system with two degrees of freedom and its restricted versions with a holonomic or a non-holonomic constraint exhibiting hyperchaotic phenomena.
We show in this comment that hyperchaos cannot appear neither in Hamiltonian systems with two degrees of freedom nor in its restricted versions with a holonomic or a non-holonomic constraint. This fact follows from the well-known properties of Lyapunov exponents for Hamiltonian systems. Incorrect conclusions of [1] were caused probably by insufficient accuracy of numerical calculations, or by too small number of iterations. Moreover, the initial conditions chosen in [1] for restricted Hamiltonian systems do not satisfy equations of constraints.
Having in mind the importance of the subject and in order to be constructive, we show hyperchaotic phenomena in Hamiltonian systems with more than two degrees of freedom through examples with triple and quartic pendula.
This comment is organized as follows. In Sect. 2, we recall basic properties of Lyapunov exponents that will be used in explanations of results of our corrections made for numerical calculations in [1]. In Sect. 3 will be presented the basic Hamiltonian system considered in this article and with preliminary explanations and description of its dynamics on the base of Poincaré sections. Section 4 contains corrected Lyapunov expo-nents spectra with the corresponding bifurcations diagrams. Section 5 concerns the restricted version of the considered Hamiltonian system with a holonomic or a non-holonomic constraint. It starts from the general explanations concerning the applied restriction procedure. Then in two subsections is considered the analysed Hamiltonian system with a holonomic or a nonholonomic constraint, respectively. Finally, in Sect. 6 we show the true hyperchaotic dynamics in Hamiltonian multiple pendulum systems.

Properties of Lyapunov exponents
Let us consider the general autonomous system of ordinary differential equationṡ Let y = R(x) will be an invertible transformation. The transformed system takes the forṁ where R denotes the Jacobi matrix of the transformation with respect to its argument. If the vector field v has the property R (R −1 (y))v R −1 ( y) = −v( y) (3) or equivalently then we say that system is time reversal because transformation (x, t) → (R(x), −t) preserves the form of system (1). For time-reversal systems, if λ is a Lyapunov exponent, then also −λ is a Lyapunov exponent [2][3][4].
The most known symmetry of this type is R(q, p) = (q, − p) for classical Hamiltonian systems. Benettin et al. in [5] proved that for a Hamiltonian system the Lyapunov exponents obey the pairing rule where N denotes the number of degrees of freedom.
Symmetries and conservation laws have also influence on a spectrum of Lyapunov's exponents. Namely, their presence manifests itself by vanishing of exponents, see Sec. 2.5.6 in [6]. For an autonomous dynamical system, its invariance with respect to time shift, which gives one zero Lyapunov exponent, is the example of such symmetry. Preservation of volume in the phase space implies that sum of all Lyapunov exponents must be equal zero, see, for example, Sec. 2.3 in [7]. Finally, zero exponents may also (non-generically) occur at bifurcation points, where some direction is (linearly) marginally stable and for periodic or quasiperiodic orbits [6,7]. For dissipative systems the sum of all Lyapunov's exponents must be negative. Moreover, the Lyapunov exponents characterize the attractor, see e.g., Sec. 5.3 in [8].
For autonomous Hamiltonian flows, one pair of Lyapunov's exponents is always zero due to time shift continuous symmetry and conservation in time of Hamiltonian. Moreover, the presence of any additional first integral, functionally independent with the Hamiltonian, leads to vanishing of another pair of Lyapunov's exponents. Thus, the determination of a full spectrum of Lyapunov's exponents can be used as an indicator of a number of independent integrals of motion that a considered dynamical system may possess. A detailed explanation can be found in the classical books [6,[9][10][11], while for the rigorous proofs please consult [2,5].
According to the generally accepted terminology for autonomous differential systems, chaos appears when at least one Lyapunov exponent is positive, while hyperchaos is manifested by at least two positive Lyapunov exponents. Thus, the minimal dimension for a dissipative non-time-reversal hyperchaotic system is 4 as, for example, in the first hyperchaotic system due to Rösler [12], and for time-reversal ones, the minimal dimension is 5 [13]. For Hamiltonian systems, hyperchaos can appear when the number of degrees of freedom is at least 3.

Description of the system its dynamics and motivation
General three-parameter family of Hamiltonian system analysed in [1] is given by the Hamiltonian function where (q 1 , q 2 ) are the generalized coordinates and ( p 1 , p 2 ) their corresponding momenta, and θ, η, β ∈ R + . Equations of motion generated by Hamiltonian (6) are the following Since evolution of system (7) takes place in a fourdimensional phase space, we can get a quick insight into its dynamics by making several Poincaré cross sections for fixed values of the parameters, for instance, β = 1, θ = 1, η = 0.01. Points creating patterns presented in Fig. 1 were obtained as traces of intersections of orbits calculated numerically with a properly generic surface of section q 2 = 0 with p 2 > 0, in a three-dimensional hypersurface defined by a constant energy level H = E. Figures are ordered according to increasing energy of the system. They show that, for the chosen values of the parameters, the system possesses chaotic behaviour and it is in general not integrable. As typical in Hamiltonian systems, we can observe a coexistence of three types of motion: periodic (finite numbers of points of orbit), quasi-periodic (closed loops) and chaotic (scattered points). We do not observe an appearance of strange attractors because Hamiltonian systems are divergence-free, and thus, there are no contractions of volume in the whole phase space (see Liouville theorem). Thus, it is unclear why the authors of [1] have used the notion of "hyperchaotic attractor" in the captions of Poincaré maps visible in Figs. 10, 11, 16 of paper [1], which is misleading.
Although the periodic, quasi-periodic and chaotic behaviour of the system can be easily detected in the Poincaré sections, its hyperchaotic behaviour can be recognized only in their Lyapunov exponents spectrum. Indeed, authors of the paper [1] have been searching for hyperchaotic behaviour of system (7) by computing the Lyapunov spectrum for certain values of the parameters. Let us mention just three of them. Namely, for the chosen initial condition and for specified values of the parameters they calculated Lyapunov's exponents spectra using MATLAB, without specifying by which method/ algorithm these exponents were calculated. The results are depicted in Figs. 2, 4 and 5 in [1]. Surprisingly enough, the obtained spectra show that system (7) is hyperchaotic for a wide range values of the parameters. Moreover, for θ = 1, η = 0.01 over the range β ∈ [1.131, 1.735], the system has three positive Lyapunov exponents.
As it is explained in Sect. 2 for a two-degree-offreedom Hamiltonian system, such as system (7), the Lyapunov exponents spectrum is symmetric about zero, i.e.
= {λ, 0, 0, −λ}. Depending on chosen initial conditions, periodic and quasi-periodic motions manifest through all zero Lyapunov exponents, while for chaotic orbit we have λ > 0. From the fundamental reasons, these systems cannot possess a hyperchaotic behaviour as the authors of [1] claimed. In fact, Figs. 2, 4 and 5 in [1] show that accuracy of numerical calculations is not sufficient because two Lyapunov exponents, which should vanish, differ significantly from zero.
In the next section, we give properly calculated Lyapunov exponents spectra of system (7) for the same values of the parameters and initial conditions as the authors of the work [1] have done.

Lyapunov exponents spectra and bifurcations diagrams for system (7)
Algorithms for computing the numerical values of Lyapunov's spectra of continuous dynamical systems have been well established over the years [5,[14][15][16][17][18][19][20]. These different methods are mainly focused on the so-called standard algorithm introduced by Benettin et al. [5,14]. This method is based on the integration of variational equations for n initial conditions with successive applications of the Gram-Schmidt orthonormalization procedure. This procedure can be easily implemented in Mathematica [21]. In this work, we use the standard algorithm with adopting a sufficient amount of k steps, so that the convergence of Lyapunov exponents is  (7) made for the parameters β = 1, θ = 1, η = 0.01, with cross-plain q 2 = 0 and p 2 > 0 assured. We keep the relative and absolute error up to 10 −11 . Figure 2 presents the Lyapunov exponent's spectra of system (7) made for parameters (9) under initial condition (8). As we can notice, these plots are completely different comparing to the ones calculated in [1]. In contrary to [1], computed exponents obey the paring rule (5), as it should be for Hamiltonian system. Of course, there is no hyperchaotic behaviour. Instead of that, we can observe the coexistence of regular and chaotic behaviour depending on values of the control parameter.
From the Lyapunov exponents spectrum, we cannot deduce whether the regular pattern corresponds to periodic or whether to quasi-periodic motion. In order to detect periodic orbits, regular regions (windows) between chaotic layers in chaos, etc., a bifurcation diagram can be constructed. A bifurcation diagram gives insight to the dynamics of a considered system  (7) versus the control parameters with initial condition (8) by plotting dependence a state variable with respect to a change of a control parameter [22].
In Fig. 3 of paper [1], the authors present the bifurcation diagram, by plotting the time series of variable p 1 (t) as a function of a certain control parameter. However, from this diagram we cannot deduce any conclusions since chaotic, quasi-periodic and periodic motions are not distinguishable. In order to get more useful information, we make the following. For given values of constant parameters β = 1, θ = 1, and for initial condition (8), we effectively constituted the bifurcation diagram by sampling points periodically in the same way as the Poincaré sections with surface of sections chosen by q 2 = −1, with gradually increasing the control parameter η ∈ (0, 001, 10). Figure 3a presents the comparison of the largest Lyapunov exponent with the bifurcation diagram. As the first sight, we can observe their very good agreement. From the bifurcation diagram, we are able to select these values of the control parameter for which system possesses chaotic, quasi-periodic and periodic behaviours. For instance, Fig. 3b shows the enlargement of Fig. 3a, which is depicted in colours values of the control parameter η for which system possesses three different types of motions. For better understanding, particular orbits corresponding to these types of motion with marked points of intersections with Poincaré cross-sectional plane are presented in Fig. 4.
It is worth mentioning that in the case of Hamiltonian systems results for bifurcation diagrams and Lyapunov exponents depend strongly on the chosen initial condition. To show that, we present in Fig. 5 the largest Lyapunov exponent with its corresponding bifurcation diagram, computed for the same values of the parameters as Fig. 3 but for different initial conditions. As we can notice, this bifurcation diagram presents completely different structure. First of all, it is symmetric about zero and, moreover, in this diagram can observe the periodic "windows" between completely chaotic layers.

Dynamics of restricted systems
In Section 3 of [1], authors consider Hamiltonian system (7) subjected to constraints. There are many approaches to describe such systems depending on types of constraints and other conditions. Authors of [1] have chosen approach described in [23]. Here, we is restricted by a single constraint ϕ(q, p) = 0 satisfying ∇ p ϕ(q, p) = 0. The constrained Hamilton's equations have the forṁ where term Q(q, p) describes additional 'forces' caused by the constraint. Its explicit form derived in where In the above formula, the gradient of a function is considered as a row vector, while the partial derivative as a column vector. For details, see [23], where also the more general cases with several constraints were considered.
As it is explicitly stated in [23], the constrained dynamics of the system takes place on the manifold Thus, to study properties of the restricted system we have to consider only those solutions (q(t), p(t)) of constrained Hamilton's equations (13) which lie on Σ. To ensure this, it is enough to require that (q(t 0 ), p(t 0 ) ∈ Σ for a certain t 0 . This is because we have the following. ϕ(q, p) is a first integral of restricted system (13).

Lemma 1 Function
But using definition Q given by (14), we easily obtain that The above shows that the total time derivative of ϕ(q, p) vanishes identically, so ϕ(q, p) is a first integral.
Without an additional justification, a study of solutions of constrained Hamilton's equations (13) with initial conditions (q(t 0 ), p(t 0 ) / ∈ Σ, as it was made in [1], does not make big sense because these solutions are related to neither the constrained nor unconstrained dynamics of the system.
For Hamiltonian system (12) restricted by a holonomic constraint ϕ 0 (q) = 0, we define Then, by Lemma 1, we know that ϕ(q, p) is a first integral of constrained system (13). Moreover, the total time derivative of ϕ 0 (q) is The above facts show that 2(n − 1)-dimensional manifold is invariant with respect to the phase flow generated by system (13). That is, if an initial condition (q(t 0 ), p(t 0 )) belongs to Π , then the whole solution (q(t), p(t)) of the constrained Hamilton's equations (13) belongs to Π . Moreover, only these solutions describe the dynamics of constrained system. If we introduce local parameterization of Π , then the constrained Hamilton's equations are canonically mapped into Hamilton's equations in R 2n−2 with respect to a certain symplectic structure; for details, see Proposition 2.1 in [24].

Holonomic constraint
Following [1], we first consider Hamiltonian system (7) with the following holonomic constraint for which, according to (17), we have Then, taking into account formulae (15), we obtain (21) and the constrained Hamilton's equations (13) take the form where The constrained dynamics takes place on the twodimensional manifold Π , which is given by On this manifold, the dynamics is described by solutions of an autonomous Hamiltonian system with one degree of freedom so it is integrable and this precludes any chaotic behaviour.
One can ask whether we can find a chaotic behaviour in the whole phase space of the constrained Hamilton's equations (13) investigating its solutions with arbitrary initial conditions (q 1 , q 2 , p 1 , p 2 ) ∈ R 4 . Numerical experiments show that for arbitrary solution of system (22) all Lyapunov exponents vanish. To understand why it is so, we introduce polar coordinates (r, ϑ) and the corresponding momenta ( p r , p ϑ ) in the standard way q 1 = r cos ϑ, q 2 = r sin ϑ, In these coordinates, the restricted system (22) takes the form Equations for r and p r form the closed subsystem and constraints take the forms Thus, the invariant manifold Π is given by On this manifold, system (23) reduces to the system with one degree of freedoṁ We can obtain the same result using the standard procedure. In fact, first we take the Lagrange function corresponding to the Hamiltonian (6). Then, as local coordinates on the constraint, we can choose the polar angle ϑ by setting q 1 = L cos ϑ and q 2 = L sin ϑ. As a result, we obtain the Lagrange function corresponding to system (6), which is The canonical momentum conjugated with ϑ is p ϑ = ∂ L/∂θ = L 2θ . Neglecting C 0 the corresponding Hamiltonian is exactly (27).
Fact that we obtained completely different results than those in [1] needs an explanation. First of all in the cited paper, the authors do not investigate the constrained Hamilton's equations (22) where and L is considered as a parameter. Although this system coincides with (22) on manifold Π , it has completely different global properties. For instance, function (17) is a first integral of equations (22), but it is not a first integral of (28).
As system (28) has time-reversal symmetry R(q, p) = (q, − p), the Lyapunov exponents obey the paring rule (5), and Lyapunov's exponents spectrum is symmetric about zero: {λ, 0, 0 − λ} with λ > 0. Thus, in system (28) the hyperchaotic behaviours are precluded. To confirm this, we present in Fig. 6 the Lyapunov exponents spectrum of system (28) made for the parameters η = 0.01, L = 1, β ∈ [0, 10] and initial condition (8), which clearly not lie on Π . This spectrum was made for the same values of the parameters and the Fig. 6 Lyapunov exponents spectrum of system (28) versus the control parameter β with initial condition (8). The remaining parameters were specifies as: η = 0.01, and L = 1 initial condition as the spectrum depicted in Fig. 7 in [1], where the authors have found hyperchaos for the entire range of parameter β. As we can notice, Fig. 6 is symmetric about zero, and the Lyapunov exponents obey the parting rule, as it should be.

Non-holonomic constraint
As the second example, we follow [1], where the authors considered Hamiltonian system (7) subjected to a non-holonomic constraint of the form For this constraint, formulae (15) give where I 2 is two-dimensional identity matrix. Thus, constrained system (13) looks as follows: Fig. 7 Lyapunov exponents spectrum and the bifurcation diagram of system (33) versus the control parameter β with crossplain p 2 = 0. The remaining parameters were specified as: θ = 1, η = 1, l = 10, and the initial condition was chosen by (8) where This system has a quadratic first integral of the form Moreover, the divergence of system (31) is In [1], the authors did not investigate the truly constrained Hamiltonian system (31) under constraint (29) but its modification where l was fixed and considered as a control parameter. Again, although this system agrees with the one given in Eq. (31) on constraint p 2 1 + p 2 2 = l 2 , it has very different global properties. For instance, the Lie derivative of function (32) along the vector field (33) is given by so it vanishes only on constraint. In Figs. 12-17, of the work [1] the authors have performed the numerical analysis of their "constrained" system (33) with initial conditions not satisfying the constrained condition p 2 1 (t 0 ) + p 2 2 (t 0 ) − l 2 = 0. They have found certain values of the parameters for which system (33) is hyperchaotic and, moreover, it has three positive Lyapunov exponents. Of course, it is not true and their wrong conclusion was probably caused by not sufficient accuracy of numerical calculations. Since system (33) possesses time-reversal symmetry R(q, p) = (q, − p), the Lyapunov exponents again obey the paring rule. Thus, in system (33) the hyperchaotic behaviours are precluded because only one Lyapunov exponent can be positive. To confirm this, we present in Fig. 7 the Lyapunov exponents spectrum and its corresponding bifurcation diagram of system (33). They were made for the same values of the parameters and the initial condition as Fig. 12 of in [1]. Since system (31) possesses first integral (32), we can reduce its dimension by one. For this purpose, we introduce new variables in the following way In these variables, constrained system (31) reads 3 2 cos ϑ , p = 0. (36) As we can see p = const, and thus, we may consider it as a parameter p = l. Hence, the final form of the reduced, constrained system (31) is given by It has time-reversal symmetry given by R(q 1 , q 2 , ϑ) = (−q 1 , −q 2 , ϑ). In Figs. 8 and 9, the Poincaré cross section, the bifurcation diagram and its corresponding Lyapunov spectrum are visible. As we can notice, for the chosen values of the parameters reduced system (37) still possesses very rich and complex dynamics. However, only one Lyapunov exponent λ > 0, the second by time-reversibility symmetry is −λ and the third vanishes.

Hyperchaos in pendulums systems
As we have already mentioned, the hyperchaotic behaviours in Hamiltonian systems can appear only in systems with at least three degrees of freedom. Multiple pendulums seem to be natural candidates for such systems with many degrees of freedom. In this short section, we will search for the hyperchaotic behaviours of triple and quartic pendulums, which move under the influence of the constant gravity field. For this purpose, we will compute full Lyapunov's exponents spectra and the bifurcations diagrams corresponding to them. Simple triple pendulum is composed of three simple pendulums attached to each other. The system has three degrees of freedom, and it is fully described by three angles ϕ 1 , ϕ 2 , ϕ 3 . Using the notations given in Fig. 10, we can write the corresponding Lagrange function in the form where for simplicity we have introduced Figure 11 presents the Lyapunov exponents spectrum and the bifurcation diagram for the simple triple pendulum system defined by Lagrange function (38). This diagram plots extremal values (amplitudes) of the system variable φ 3 (t) = φ 3,extr , whenφ 3 = p 3 = 0. For values of the parameters m i = l i = g = 1 with the initial condition: (φ i = φ 0 , p i = 0) for i = 1, 2, 3, we show how the change of initial angle of the pen- Fig. 9 Lyapunov exponents spectrum and the bifurcation diagram of system (37) versus the control parameter θ with crossplain ϑ = 0. The remaining parameters and the initial condition were specified as: β = 1, η = 1/100, l = 1, (q 10 , q 20 , ϑ 0 ) = (0, 2, 0) dulums deviation φ 0 ∈ [0, π] affects the dynamics of the system. As one would expect, for small angles the pendulums oscillate near the equilibrium. However, when the angular displacement amplitude of the pendulums is large enough, the regular pattern diverges and the motion finally becomes highly complex. The corresponding Lyapunov exponents spectrum indicates that motion is, in fact, hyperchaotic with two positive exponents.

Flail triple pendulum
Let us now look for the hyperchaotic nature of the flail triple pendulum moving under the influence of the grav- itational force. The geometry of the system is presented in Fig. 12. It consists of three simple pendulums with masses m 1 , m 2 , m 3 and lengths l 1 , l 2 , l 3 . The first pendulum is attached to a fixed point, and to its end mass, the other two pendulums are joined. This system, which was in the absence of the gravity field intensively studied in [25], has the following Lagrange function: (39) Figure 13 presents the Lyapunov exponents spectrum and its corresponding bifurcation diagram versus the initial angle of deviation φ 0 ∈ [0, π]. These plots were made for fixed values of the parameters m i = l i = g = 1 with the initial condition:(φ i (0), p i (0)) = (0, φ 0 , −φ 0 , 0, 0, 0) for i = 1, 2, 3. As in the case of the simple triple pendulum, φ 3,extr states for the extremal values of the system variable φ 3 , whenφ 3 = p 3 = 0. As expected, for small values of φ 0 , the pendulums oscillate near the equilibrium, and we can conclude that for such small values of φ 0 the motion is almost integrable. Nonetheless, for φ 0 ≈ π/2 the motion starts to be more complex and finally becomes chaotic and later hyperchaotic, which is evidenced in Fig. 13. Surprisingly enough, over the range φ 0 ∈ [2.246, 2.293], we can still detect stable solutions localized in the gap between completely hyperchaotic regions. This corresponds to the case when the first pendulum is

Quartic simple pendulum
Let us finally consider the simple quartic pendulum. This system is composed of four simple pendula of masses m i and lengths l i which are connected to each other and move under the influence of the gravity, as shown in Fig. 15. The Lagrange function of this model This system is hyperchaotic with three positive Lyapunov exponents, provided the common initial amplitude φ 0 is large enough. For instance, for m i = l i = g = 1 with the initial condition: (φ i (0) = φ 0 , p i (0) = 0), the motion starts to be hyperchaotic for φ 0 > π/3, as shown in Fig. 16.