Dynamics and integrability of the swinging Atwood machine generalisations

This paper studies the dynamics and integrability of two generalisations of a 3D Swinging Atwood’s Machine with additional Coulomb’s interactions and Hooke’s law of elasticity. The complexity of these systems is presented with the help of Poincaré cross sections, phase-parametric diagrams and Lyapunov exponents spectrums. Amazingly, such systems possess both chaotic and integrable dynamics. For the integrable cases we find additional first integrals and we construct general solutions written in terms of elliptic functions. Moreover, we present bifurcation diagrams for the integrable cases and we find resonance curves, which give families of periodic orbits of the systems. In the absence of the gravity, both models are super-integrable.


Introduction
An amazing property of the laws of nature is the fact that analytical formulas can be used to describe simple models that are built to illustrate these properties. Such models are said to be solvable. This property is closely related with conservation laws. In the case of models described by ordinary differential equationṡ the conservation laws are first integrals. A first integral of equations of motion (1) is a non-constant scalar function F(x) : R n → R, which does not change during time evolution of the system, i.e., F(x(t)) = F(x(0)) for an arbitrary t. If we are able to find a large enough number of independent first integrals, then the system is called integrable and its solutions can be found. However, the absence of first integrals indicates that the behaviour of our system is complicated or chaotic. If we know only a few first integrals (fewer than the number necessary for integrability), we can reduce the dimension of the system. For most systems of a physical origin, the bestknown examples of first integrals are the energy, linear momentum, and angular momentum of the system. These properties are directly related to conservations of certain physical quantities. Moreover, even in a case of simple models, the universal conservation laws are not sufficient to describe them analytically. For instance, to solve the heavy top problem, it is not enough to know the laws of conservation of energy, and the conservation of the projection of the angular momentum into vertical directions. Additional first integrals exist only for special inertia tensor and the location of the fixed point in the body. Moreover, distinguishing exceptional cases when the heavy top problem is integrable is not at all obvious.
A second system we discuss which is integrable for a certain set of values of parameters is the classical Swinging Atwood's Machine (SAM). This is an ordinary Atwood's machine in which one of the weights can swing in a plane. This model, which is the system of two degrees of freedom, was intensively studied by Tufillaro and co-workers, see [1][2][3][4][5][6][7]. They performed a numerical analysis via the Poincaé sections and they showed that for most values of the mass ratio M/m > 1 SAM is not integrable displaying highly chaotic behavior. However, as they noted, for the mass ratio M/m = 4n 2 −1 where n ∈ Z, the system behaves less chaotically and motion looks regular. Indeed, after introducing the parabolic coordinates, Tufillaro in [8] showed that Hamilton-Jacobi equation for the classical SAM system for M/m = 3 is separable and the equations of motion can be solved in terms of elliptic functions, which is quite unusual for systems with two degrees of freedom.
For some time, it was unclear whether SAM system could have other integrable cases. However, after five years of research, thanks to new criterion for the integrability of homogeneous Hamiltonian systems obtained by Yoshida in [9], the non-integrability of the system for M/m = μ p , where μ p = p( p +1)/( p( p + 1) − 4) for every integer p > 2 was proved [5]. This result was later completed for the case p = 2 in [10] by Martínez and Simó. Currently, the SAM is considered to be a paradigm of a simple, mechanical system that exhibits chaotic as well as regular properties. Indeed, one can find numerous web pages and video clips devoted to constructions of real SAM models and reviews of its dynamics. However, it seems that a study of the dynamics and integrability of certain generalisations of the SAM model do not attract sufficient attention. To the best of our knowledge, only two modifications of the classical SAM system exist, both of which were considered in [11][12][13][14]. Both systems, however, manifest complicated and, in fact, not integrable dynamics for all values of the parameters.
The only integrable generalisation of the SAM model that was found was the one from a threedimensional swinging Atwood's machine [15]. As in the 2D version, the problem is integrable when the mass ratio M/m = 3. However, what is particularly interesting is the fact that the additional first integral is a polynomial of the fourth degree with respect to momenta. Unlike the 2D case with quadratic first integrals, integrable Hamiltonian systems with first integrals of higher degree do not necessary possess separable solutions of Hamilton-Jacobi equations. Such systems are absolutely unique and only a few examples are known. This is why in our paper we focus on dynamics and integrability studies of two generalized classes of 3D SAM systems with additional Coulomb's interactions between two point charges and Hooke's law of elasticity. We show that, in general, both models are not integrable and chaotic. However, for specific values of parameters, the models admit additional first integrals. We use these first integrals to construct general solutions of equations of motion.
The paper is organised as follows: A description of the systems and their complex dynamics analysis is given in Sects. 2 and 3 respectively. For this purpose, we calculate the Poincaré sections, phase-parametric diagrams and Lyapunov exponents spectrums. Moreover, we present the integrable cases that were found and we construct phase-parametric diagrams of the integrable systems as well as the general solutions of equations of motion. The composition of both models and its integrability is discussed in Sect. 4. Final conclusions are drawn in Sect. 5.

Presentation of the system and its equations of motion
In Fig. 1, the geometry of system under consideration is presented. As in the classical SAM, the model consists of two masses M and m linked by an inextensible string of length l. The mass m is allowed to swing in 3dimensional space, while the second mass M plays the role of a counterweight and it is constrained to move only along the vertical direction. The constant gravitational field acts on the masses, so their gravitational potential energy is V g = g(mz − M X). The mass m swings in 3-dimensional space, while M is constrained to move along the vertical direction. The whole system moves under the influence of the potential (10). The Hamiltonian function corresponding to this model is defined in (9) We propose the following modification of the classical SAM. Namely, in the centre of the left pulley, we put an electric static point charge Q. The mass m has a charge q, so the Coulomb potential of these two charges is V c = −γ / x 2 + y 2 + z 2 , where γ is proportional to Q q. We assume that the mass M is electrically neutral.
The Lagrange function of the system is as follows The holonomic constraint imposed on the system is Thus, we introduce the spherical coordinates according to the constraint In these coordinates, the Lagrange function (2) takes the form where Since ϕ is the cyclic variable, i.e., ∂ L/∂ϕ = 0, the corresponding generalized momentum is the first integral of motion. The existence of this first integral corresponds to the axial symmetry of the system. Performing the Legendre transformation we obtain the Hamiltonian function The system has two degrees of freedom and the effective potential reads In the case when Π = 0, the pendulum with mass m oscillates on a plane, while for Π = 0, the pendulum can perform spacial oscillations.
Hamilton's equations of motion generated by Hamiltonian (9) form a system of four first-order ordinary differential equations The time dependence of cyclic variable ϕ can be determined by

Poincaré cross-sections
We present the complex dynamics of the system under consideration with the help of numerical analysis.  The plots indicate the complex dynamics of the system. Each colour corresponds to distinct initial condition Cross section plane was specified as r = 1, with p r > 0. The plots present chaotic dynamics of the system. Each colour corresponds to distinct initial condition system (9) with gradually increased values of energy. These plots were constructed for two sets of constant parameters The section plane was specified as r = 1 and a crosssection is marked on the figure if p r > 0. Figure 2 was made for a zero level of the cyclic first integral Π = 0, which implies that mass m moves in the vertical plane. The Poincaré sections visible in Fig. 3 were constructed for a non-zero value of the cyclic first integral, i.e., Π = 1/4, which corresponds to the case that mass m performs spacial oscillations. Figures 2a and 3a present Poincaré sections for the lowest energies. As is typical in Hamiltonian systems, for energies close to the energy minimum, the patterns are very regular. In the central part, the stable particular periodic solution is visible. It is surrounded by closed loops related to the quasi-periodic motion of the system. They are just intersection of invariant tori. It is noticeable that the amplitude of the oscillations of mass m is very small. We may conclude that for small values of the energy, the system is very close to an integrable one.
The picture becomes more complex when we increase the value of the energy. In Fig. 2b and in Fig. 3b, we can notice that the invariant tori are visibly deformed and some of them decayed into sequences of stable and unstable periodic solutions. For instance, at the prescribed energy levels, the stable particular solutions depicted in the centre of Figs. 2a and 3a have been bifurcated, which led to the emergence of stable periodic solutions enclosed by the remains of the separatrix. Moreover, we observe the first signs of chaotic behaviour of the system visible in the plane as scattered points appearing in the regions where the separatrices where located.
The last pair of Poincaré sections, visible in Figs. 2c and 3c, show highly chaotic stages of the system dynamics. Indeed, almost all invariant tori located in the central parts decayed into the global chaos. For a not integrable Hamiltonian system, these random-looking points correspond to the fact that trajectories can freely wander over larger regions of a phase space. Trajectories (in contrary to those in an integrable system) are no longer confined to surfaces of a set of nested tori, but they begin to move outside the tori. We may conclude that the Hamiltonian (9), for the prescribed values of the parameters, is highly not integrable.

Spectrum of Lyapunov exponents
The Poincaré sections visible in Figs. 2 and 3 show the beauty of coexistence of three types of motion: periodic (finite numbers of points of orbit in the cross-section plane), quasi-periodic, and chaotic (scattered points). However, from these plots, we cannot deduce any information about the strength of chaos. For this purpose, the Lyapunov exponents can be computed. The use of Lyapunov exponents is an essential tool for quantifying chaos in dynamical systems. It measures the exponential divergence of two close trajectories (orbits) in a phase space. According to the chaos theory, for autonomous differential systems, chaos appears when at least one Lyapunov exponent is positive [16,17].
In this paper, we use the standard algorithm introduced by Benettin et al. [18]. This method is based on the successive integrations of variational equations for n initial conditions together with applications of the Gram-Schmidt orthonormalization procedure, see also [19]. This procedure can be effectively implemented in Mathematica [20]. We adopt a sufficient amount of k steps, so that the convergence of Lyapunov exponents is assured. Figure 4 presents the Lyapunov exponents spectrum of the system (11) made for constant values of the parameters (14) with the initial conditions r 0 = 1.05,ṙ 0 =θ 0 = 0, where ϑ 0 ∈ (0, π), (15) is the control parameter. Here, {λ, λ 2 , λ 3 , λ 4 } denotes the full spectrum, where λ ≡ λ 1 is the largest Lyapunov exponent. Figure 4 presents how the change of the initial swing angle of the pendulum (with zero initial velocities) affects the dynamics of the system. We can observe the coexistence of regular and chaotic behaviour depending on the values of the control parameter. Specifically, regions of the spectrum with positive values of λ are responsible for chaotic motion of the system, while regions with λ ≈ 0 indicate regular (non-chaotic) behaviour. First, we can notice that the spectrum obeys the pairing rule, i.e., Lyapunov exponents appear in additive inverse pairs. This is due to the time reversibility of the Hamiltonian vector field (11). Moreover, the existence of the first integral, which is the conservation of the energy H = E, implies that one pair of Lyapunov exponents is zero [18,[21][22][23]. Hence, Fig. 4 Lyapunov exponents spectrum of system (11) versus the initial swinging angle ϑ 0 ∈ (0, π) with the initial condition (r 0 = 1.05, ϑ 0 = ϑ 0 ,ṙ 0 =θ 0 = 0). The constant parameters where chosen as: M = 2, m = 1, g = 1, γ = −0.5, Π = 0.25. Here, {λ, λ 2 , λ 3 , λ 4 } denotes the full spectrum, where λ ≡ λ 1 is the largest Lyapunov exponent. The intervals with positive values of λ are responsible for the chaotic motion of the system, while the regions with λ ≈ 0 indicate regular (non-chaotic) behaviour for the Hamiltonian system (11), which is the system of two degrees of freedom, the possible spectrum is {λ, 0, 0, −λ}. In the case considered, the maximum pick value of λ is at point ϑ 0 ≈ 2.9, where λ ≈ 0.05.

Phase-parametric diagram
The spectrum of Lyapunov exponents visible in Fig. 4 gives us an insight into the system dynamics by specifying the intervals of the control parameter ϑ 0 where the chaotic motion takes place. However, from this figure, we are not able to classify whether the regular patterns with λ ≈ 0 are responsible for periodic or quasi-periodic motion. For this purpose, a phase-parametric diagram is very useful. For Hamiltonian systems, the phase-parametric diagram enables the detection of periodic orbits and their families, routes to chaos and "periodic windows" between chaotic layers   by plotting a state variable as a function of a certain control parameter [24,25]. Figure 5 presents the phase-parametric diagram of the system (11) constructed for the same values of the parameters (14) and initial conditions (15) as the Lyapunov exponents spectrum visible in Fig. 4. For the calculated phase-parametric diagram, we display the extreme values of ϑ(t) versus the initial swing angle ϑ 0 ∈ (0, π). That is, for a given initial condition, we numerically integrate equations of motion (11), and we constituted the diagram by sampling points ϑ(t ) periodically whenθ(t ) = 0. In order to make the analysis more effective, we combine the phase-parametric diagram with the maximal Lyapunov exponent λ. The colour scale is proportional to value of λ.
At first sight, we observe very good agreement of the phase-parametric diagram with λ. Moreover, while not being visible via Lyapunov exponents, we can distinguish from the regular regimes periodic and quasiperiodic orbits. This is due to the fact that we plot the amplitudes of ϑ(t) not its time series as was done for another model in [26]. Thanks to that, we are able to select from the computed phase-parametric diagram the values of the control parameter, ϑ 0 , for which the system possesses periodic, quasi-periodic and chaotic behaviours. The exemplary values of ϑ 0 corresponding to these types of motion are marked on the x-axis of Fig. 5b. For better understanding, we also plot these particular orbits with marked points of intersections planeθ(t) = 0 in Fig. 6. Figure 5 shows, moreover, that for the initial range of ϑ 0 the system performs regular, quasi-periodic oscillations up to the point ϑ 0 ≈ 0.125, where the periodic solution appeared. However, we claim that this periodic solution is unstable because in the right-neighbourhood of this point: at ϑ 0 ≈ 0.13, to be precise, we observe the chaotic behaviour in terms of the large pick with λ ≈ 0.038. Further, over the range ϑ 0 ∈ (0.13, 2.02], the motion is regular with mostly quasi-periodic orbits. As ϑ 0 is increased, the regular pattern diverges up to the appearance of global chaos, where the orbits' maxima densely occupy the large admissible region. However, as is evidenced in Fig. 5c, even for very large values of the initial amplitude ϑ 0 , we still detect stable solutions localized in the gaps (windows) between completely chaotic regions. Finally, Fig. 7 shows the regular and chaotic orbits of system (11) plotted in 3D Euclidean space. These plots confirm that the motion of the pendulum is spatial. Respective values of ϑ 0 were taken from the phase-parametric diagram (see the coloured labels in Fig. 5b and also Fig. 6).

The separable case
From the Poincaré sections, the phase-parametric diagram and Lyapunov exponents spectrum, we can conclude so far that the system seems to be not integrable. However, numerical analysis was performed for certain fixed values of the parameters. For other values, the dynamics can be completely different and the system can be even integrable. For instance, Fig. 8 presents the Poincaré sections of the system obtained for values of parameters at energy level E = 10. As we can notice, even for such high value of the energy, the dynamics of the system is very regular. Moreover, the Lyapunov exponent spectrum and the phase-parametric diagram (see Figs. 9 and 10) also indicate the presence of an additional first integral, functionally independent with the Hamiltonian. This is because another pair of Lyapunov exponents tends to zero for all values of the initial swing angle ϑ ∈ (0, π) of the pendulum. Hence, for parameters (16), the numerical analysis excludes chaotic behaviours of the system and, in fact, suggests its integrability. Surprisingly, more detailed numerical analysis (not included) shows that for the system seems to be integrable. However, to ensure the integrability, we have to find an additional first integral functionally independent with the Hamiltonian. It seems that for this purpose, the Hamilton-Jacobi theory could be very useful. The main idea of this theory is to introduce an appropriate canonical transformation after which Hamilton's equations become separable. However, for study the system considered, we perform a separation of variables directly. We consider the case M = 3m, and (without loss of generality) we set m = 1. Then, by means of the following canonical change of the variables the Hamiltonian (9) reads As we can notice, for Π = 0 the system separates. Namely, from equality we separate terms which depend only either on (u, p u ) or (v, p v ). Then, we find that is the additional first integral, which is functionally independent with the Hamiltonian (19). Hence, the SAM model with Coulomb's interactions is integrable at zero level of the cyclic integral, Π = 0.
Performing the inverse change of variables (18) to (20) we get the desired constant of motion written in the polar coordinates. The explicit form of the integrable system becomes It is evident that the integrable system (21) is the generalisation of the classical SAM model by the additional term in the potential coming from the Coulomb interaction. We can easily recover the classical SAM by putting γ = 0. Moreover, in the absence of gravitational potential, the system is obviously superintegrable with the additional cyclic first integral, J = p ϑ .

Bifurcation diagram of the integrable system
In further consideration we put Π = 0 and to simplify analysis, we chose units in such a way that g = 1. Let us consider a common level of the first integrals (19)-(20), namely It is easy to notice that on the common level (22), the following equalities are fulfilled. Notice also, that Eq. (23) is invariant under the change These values of (h, f ) belong to an image of the so-called momentum map A set of critical points of the momentum is defined by It is an invariant set completely filled with critical trajectories. By condition (27) set S is a union of two disjoint subsets In a typical situation, S 0 consists of isolated fixed points, while S 1 includes families of curves, which are either periodic or asymptotic to unstable equilibria. We refer to such orbits as critical orbits in order to distinguish them from periodic ones that lie in regular level sets with resonant characteristic frequencies.
A bifurcation diagram of an integrable system is defined as a region of possible motions Φ(M) on a plane of first integrals together with the images of the critical sets (29) depicted on the same plane.
Let us assume that P ∈ S 0 . Condition rank dΦ( P) = 0 implies that gradients ∇ H ( P) and ∇ F( P) restricted to M vanish. Hence, one may conclude that S 0 coincides with the set of equilibria of the Hamiltonian vector field. Indeed, set S 0 is defined by the equations Using the relation (23), we can show that in terms of parameters (γ , f, h) ∈ R 3 and variables (u, v) ∈ R 2 image Σ 0 of set S 0 on a plane of first integrals (h, f ) consists of two families of distinct fixed points P 1 , P 2 , such that with u 2 = v 2 = h/4, with v = 0 and u 2 = h/4. In order to find set S 1 of critical periodic orbits σ i ∈ Σ 1 , we need to specify the conditions when rank dΦ(P) = 1 for P ∈ S 1 . This indicates that gradients ∇ H ( P) and ∇ F( P) are parallel at P. It transpires that the only possible critical periodic solutions contained in S 1 are determined by two different families of equations Eliminating the momenta by relation (23), we obtain We remark that one can obtain the mirror family of the above equations by a change of variables defined in (24). Taking into account that the coordinates (u, v) and momenta (23) must be real, we can easily determine admissible range of critical (bifurcation) curves σ i ∈ Σ 1 on a plane of first integrals (h, f ). These curves are parametrized by γ , such that where v 2 = h/4 and −u 0 ≤ u ≤ u 0 , where v = 0 and −u 0 ≤ u ≤ u 0 . Figure 11 shows the 3-dimensional bifurcation diagram of an integrable system (19) by presenting the region of possible motions Φ(M) depicted in a (γ , h, f )-space. For better understanding, we also present in Fig. 12 the projections of Fig. 11 in a plane of first integrals (h, f ) for γ = ±1. As expected, according to invariance (24), the bifurcation diagram is symmetric with respect to the f -axis. Moreover, we can notice that Φ(M) consists of four admissible regions separated by bifurcation curves ±σ i . It is worth emphasizing that in each connected component of the level set M h, f has the same topological type. However, if constants h and f pass though the bifurcations lines, then the topological type of M h, f may change and the bifurcation of invariant tori takes place. In order to specify all topological types of M h, f over regions I-IV, the Fomenko theorem of bifurcation for Liouville tori can be adopted, see Chpt. 2, and Chpt. 6 in [27]. A more detailed description of algebraic and geometric aspects of two-dimensional integrable Hamiltonian systems can be found in [28].
For our purpose, let us show that inside the regions separated by the bifurcation lines, the values of h and f correspond to distinct dynamical behaviours of the system. First, from the bifurcation diagram, we specify values of first integrals (h, f ), and thus we can find initial conditions for (u, v, p u , p v ) from the admissible regions where trajectories are confined in M h, f . Then, by making the inverse change of variables (18), we find initial conditions for (r, ϑ, p r , p ϑ ). Figure 13 shows three exemplary solutions of the system plotted in (x, y)-plane with (h, f ) taken from admissible regions of the bifurcation diagram. Despite the fact that in all these cases the motion is quasi-periodic, the structures of the orbits and the behaviour of the swinging mass are quite different.
As we have noted previously, at the bifurcation lines, where invariant tori decay and level sets M h, f becomes degenerated, exceptional families of particular periodic solutions may appear. In Fig. 14 we present the periodic orbits of the system corresponding to values of (h, f ), which lie at the bifurcation curve σ 1 . As in the ordinary swinging Atwood machine model, we observe characteristic "loops" and "smile" trajectories. As it transpires, orbits (not presented) lying at the bifurcation curve σ 2 correspond to singular orbits approaching into the fixed point (x, y) = (0, 0). Surprisingly, we found a family of periodic but not critical solutions (rank dΦ(x) = 2), which lie at zero level of f . Figure 15 shows exemplary periodic orbits with various values of h. It will be clear shortly that zero value of f is, in some sense, also exceptional.

Derivation of a general solution
To derive explicit solutions, we perform a time transformation t → τ by setting dt dτ = 4(u 2 + v 2 ).
Then, the dynamics of the system at level H = h is described by Hamilton's equations generated by the following Hamiltonian restricted to its zero level. As from (23), we obtain where prime denotes the derivative with respect to τ and The right-hand sides of (43) are polynomials of degree four. Thus, their solutions can be expressed in terms of the elliptic functions. An explicit form of these solutions depends of the parameters' values. In general, we can represent them in terms of the Weierstrass ℘function or in terms of the Jacobi elliptic functions. We choose the latter option because it is simpler and easier for interpretation. We assume that the motion takes place within Region IV of the bifurcation diagram, i.e., we restrict In this region we can set where the real positive parameters (a, b, c, d) are defined in the following way Thus, the motion takes place within the region Next, we have to perform inversion of the following integrals Taking into account equalities (46), we can rewrite (49) by Both integrals are of the same type. Making a substitution x = A cn(y, k), for appropriate constants A and k, we immediately obtain the desired inversion formulae. This is the standard procedure described in [29]. We fix u = a and v = c, then we obtain where The function cn(x, k) has period 4K (k) where K (k) is the complete elliptic integral of the first kind. Hence, the function u(τ ) has period For the case considered, we can find an explicit relation between t and τ . Direct integration of (40) gives  where and E(x, k) = x 0 dn(s, k)ds denotes the Jacobi epsilon function.
Let us now consider motion of the system restricted to Region II of the bifurcation diagram. Thus, we assume The formula for p u remains unchanged. However, inside this region, the polynomial P − (v 2 ) has four real roots, so we set where c and d are defined by Thus, we have to invert the following integral It can be done easily -it is enough to make a substitution The If (h, f ) belongs to Region I, both polynomials P ± (x 2 ) have four real roots. This region consists of two parts. The first one is given by and the second one by In this case v(τ ) is given by formula (59). To find u(τ ) we set Then, a and b are defined by Proceeding like in the previous case, we obtain that Hence, the real period of u(τ ) is Obtaining solutions (u(τ ), v(τ )) and their periods, we are able to construct resonance curves corresponding to periodic solutions lying on regular tori. The ratio of periods of functions (u(τ ), v(τ )), defined by is called the rotation number. If R n:m is rational, then the orbit (u(τ ), v(τ )) is periodic, otherwise it is quasiperiodic. Hence, we can effectively find families of periodic solutions within the indicated regions by sampling R n:m on the plane of first integrals (h, f ). Figure 16 presents the contour plots of such families. We can read from it the initial conditions for which motion of the system is periodic with an appropriate rotational number R n:m . For instance, it is now clear that previously mentioned exceptional curves (see Fig. 15) lying at level f = 0 are, in fact, periodic orbits that lie in regular level sets with a resonant characteristic frequency R 1:1 . Moreover, in Figs. 17, 18 and 19, we present examples of trajectories with regions of possible motion plotted in the Cartesian (x, y)-plane given by numerical computations with initial values related to small ratio resonances. It is worth emphasizing that, within the regions indicated, every initial condition (u 0 , v 0 ), which can be easily translated to (x 0 , y 0 ), corresponds to an orbit with distinct shape but the ratio R n:m is the same for prescribed values of first integrals (h, f ). For instance, Fig. 18 presents orbits for R 1:2 , R 1:3 and R 2:3 , with two different initial condition each.

3D integrable case with a quartic first integral
In terms of the parabolic-cylinder coordinates (18), the Hamiltonian (9) becomes separable at zero level of first integral (7). So, we put Π = 0 and, in this case, the system admits first integral (20), which is the quadratic polynomial with respect to the momenta. However, Poincaré sections show in Fig. 8b that the system can be also integrable for Π = 1/4. In fact, more detailed analysis (not included) suggests integrability of the system for arbitrary Π ∈ R. So, let us try to confirm integrability by constructing a complementary first integral of a higher degree. The Hamiltonian (9) is an even function with respect to the momenta, and also equations of motion (11) are invariant by the time reflection t → t = −t, i.e., Hence, without loss of generality, we can restrict our search for even or odd degree polynomial first integrals. Moreover, as noted and successfully used in papers [15,[30][31][32], to build a complementary first integral, we can apply the so-called order doubling. For this purpose, the quadratic integral I 2D (21) can be adopted. We postulate that the additional first integral, quartic in the momenta, is following where a i = a i (r, ϑ) are coefficients to be determined. Calculating the time derivative of I 3D , and using equations of motion (11), we receive a polynomial with respect to the momenta. Thus, the equationİ 3D = 0 gives rise a system of first-order, nonlinear partial differential equations for unknown functions a i (r, ϑ). The solution of this system gives the final form of the first integral we look for. After, tedious symbolic calculations, it can be written as follows where Π, g, γ ∈ R. It is clear that at zero level of cyclic first integral (7), Π = 0, the quartic first integral I 3D reduces to a quadratic one I 3D = I 2D I 2D , where I 2D is defined in (21). The Hamiltonian (68) defines the new three-parameter family of an integrable system, which is the generalisation of the model studied by Elmandouh in [15]. Surprisingly, the addition of the term related to Coulomb interactions V c = −γ /r , does not destroy the structure of integrability. Let us summarize the main results of this section by the following theorem. (9) is completely integrable, for arbitrary g, γ, Π ∈ R.

Theorem 1 If M = 3m, then the 3D-SAM with static Coulomb interactions governed by Hamiltonian
The integrable system has for its explicit form defined by (68).

Presentation of the system
In Fig. 20 the geometry of the second system under consideration is presented. As in the classical case, the masses M and m are connected to each other by one string. The mass M is constrained to move only along the vertical direction, while the second mass m is allowed to move in 3D space. However, we propose the following modification. We introduce, two springs with spring constants k 1 , k 2 ∈ R + and one additional string (coloured red). Spring k 1 with natural length d is rigidly restrained. It stretches vertically and at its end we have a massless ring to which the red string is attached. This string is thread on two massless pulleys and at its end, the mass M is hung. Spring k 2 connects mass m with spring k 1 . Potential energies coming from Hooke's interactions between springs k 1 and k 2 are given by The Hamiltonian function of the system, written in spherical coordinates (r, ϑ, ϕ), is given as follows where the effective potential reads V (r, ϑ) = Π 2 2mr 2 sin 2 ϑ + gr(M − m cos ϑ) As for the first model, Π represents a value of cyclic integral J = p ϕ = Π , associated to cyclic coordinate ϕ, while s = l − d ≥ 0 has been introduced for simplicity. Hamilton's equations generated by Hamiltonian (70) are as follows where cyclic variable ϕ can be determined by relation (12).

Numerical analysis
Figures 21 and 22 present Poincaré sections of system (72). We observe the beautiful coexistence of periodic, quasi-periodic and chaotic motion. We see, in particular, the magnification of the Poincaré section presenting families of high-order resonance orbits. Hence, the ratio of the masses M/m = 3 does not imply integrability of the system for arbitrary remaining parameters which was the case in the first model. Here, values of the remaining parameters Π, s, k 1 , k 2 matter. Moreover, as evidenced in Fig. 23, due to the appearance of the harmonic interactions, the model considered for the mass ratio M/m < 1 still possesses very complex dynamics. This is particularly visible in the phaseparametric diagram presented in Fig. 24. The diagram shows the extreme values of ϑ(t) versus an initial swing angle ϑ 0 ∈ (0, π) of the pendulum. The colour scale is proportional to value of the maximal Lyapunov exponent λ. As we can notice, this figure presents completely different dynamics compared to the classical SAM model for which assumption M/m > 1 is needed for oscillations, see [4]. We observe the rich structure of high-order resonance periodic solutions coexisting with chaotic orbits in their neighbourhoods. Since a value of the cyclic integral is different from zero (we put Π = 1), the pendulum with mass m performs spacial oscillations. Some of the regular orbits plotted in 3D Euclidean space are depicted in Fig. 25. The initial conditions were taken from the phase-parametric diagram (see the coloured labels in Fig. 24).

Separable case
Numerical analysis shows the complex and, in fact, highly chaotic behaviour of the system for a wide range values of the parameters. Surprisingly, however, a more detailed analysis (not included) suggests that for the dynamics is very regular and it suggests the integrability of the system. So, let us try to ensure its integrability by constructing an additional first integral. We start by making a canonical change of dependent variables from the polar coordinates into the paraboliccylinder ones (18). For values of parameters (73) and (without a loss of generality) with m = 1 and k 2 = 3k, the Hamiltonian (70) transforms to the form It is obvious that at zero level of the cyclic integral, i.e., p ϕ = Π = 0, variables u, v are separated. Hence, we put Π = 0. Now, from equality we separate terms which depend only either on (u, p u ) or (v, p v ). Then, we find that is the additional first integral. Hence, we conclude that a SAM system with additional Hooke's potential (69) is integrable at zero level of cyclic integral Π = 0. We perform an inverse change of variables (18) on (76) to obtain a desired constant of motion written in the polar coordinates. The explicit form of the integrable system is as follows System (77) is the new integrable model. It is the generalization of the classical 2D SAM model. This generalisation is given by the appearance of the additional term in the potential which comes from the Hooke interactions between the springs.

Bifurcation diagram of the integrable system
We chose units in such a way that g = 1. Assuming Π = 0, we consider a common level of first integrals (74) and (76), specifically Values of (h, f ) belong to an image of the momentum map From equations H = 2h and F = 2 f we find that Notice, that these equations are invariant under change (24). Hence, we can restrict our consideration for f ≥ 0. In order to build a bifurcation diagram we processed as in the first model. Let S = S 0 ∪ S 1 be a set of The section plane was specified as r = 5, with p r > 0. Central region of the section is responsible for highly chaotic stages of the system dynamics. Each colour corresponds to distinct initial condition critical points P of the momentum map (79), such that rank dΦ( P) < 2. If rank dΦ( P) = 0, then S 0 contains isolated fixed points. Otherwise, if rank dΦ( P) = 1, then S 1 consists of families of periodic or asymptotic curves. It can be easily verified that set S 0 = ∅. Hence, the Hamiltonian vector field governed by Hamiltonian function (74) does not possesses any equilibrium point. To find set S 1 of critical periodic orbits, we have to find points P for which gradients ∇ H (P) and ∇ F(P) are Fig. 24 The phase-parametric diagram of system (72) versus the initial swing angle ϑ 0 . The diagram is combined together with the maximal Lyapunov exponent λ, where colour scale is proportional to the value of λ. Parameters where chosen as M = 1, m = 3, g = 1, s = 5, k 1 = 2, k 2 = 0.5, Π = 1 with initial condition (r 0 = 3.78, ϑ 0 = ϑ 0 ,ṙ 0 =θ 0 = 0). Here ϑ(t ) denotes extreme values of ϑ(t) whenθ(t ) = 0 for some t . Periodic motion diverges as ϑ 0 increases and finally becomes chaotic. We can also detect the stable "windows" between completely chaotic regions. Coloured labels correspond to the respective initial conditions for which periodic orbits are plotted in Fig. 25 linearly dependent. Direct calculations show two distinct families of equations exist, which can define possible critical periodic solutions contained in S 1 . They have the form Then, we eliminate momenta ( p u , p v ) using formulae defined in (80). Hence, we receive We remark a mirror family of the above equations exists according to invariance (24). Taking into account that the coordinates (u, v) and momenta (80) must be real, we can easily find critical curves λ i ∈ Σ = Φ(S 1 ), on a plane of first integrals (h, f ). They are as follows where v 2 = (−1+ √ 1 + 3h k)/3k and −u 0 ≤ u ≤ u 0 , where v = 0 and −u 0 ≤ u ≤ u 0 . Figure 26 shows a 3D bifurcation diagram of an integrable system (74) by presenting the region of possible motion depicted in a (k, h, f )-space. For greater detail, see Fig. 28, which presents the projection of the bifur-cation diagram to a plane of first integrals (h, f ). Due to invariance (24), the bifurcation diagram is symmetric about f -axis. Moreover, we can notice that Φ(M) consists of two admissible regions separated by critical curves ±σ i . At these curves, the bifurcations occur and exceptional families of periodic solutions appeared, see Fig. 28.

Derivation of a general solution
Next, we perform time a transformation t → τ by setting dt dτ Then, from (42) and (80), we obtain where Direct integrations of these equations lead to hyperelliptic integrals so a change of variable has to be made to reduce them to elliptic integrals. Let us set Then, we obtain The above equation defines the Weierstrass elliptic function w(τ ) = ℘ (τ, g + 2 , g + 3 ), with invariants g + 2 = 2 10 3 Thus, from Eq. (88) we have In the similar way, we obtain where As we already have noted, the case with f = 0 is exceptional because at this level, invariant tori bifurcate. Figure 28b-c presents two different closed trajectories. As its turns out, for a small initial angle ϑ 0 ,    Fig. 28 The loop, the teardrop and the heart orbits of the system plotted in (x, y)-plane. In the all cases, the initial conditions were specified as u 0 = 1, v 0 = −1, which implies r 0 = 1 and ϑ 0 = 0 we get the teardrop orbit, while for a sufficiently large ϑ 0 the orbit is called the heart orbit. This unexpected feature was observed in the ordinary SAM system by Tufillaro [6]. In our model, however, we have the additional term in the potential related with the Hooke interactions.
Let us make a deeper insight into the structure of an explicit form of solutions. We put f = 0. Then, from Eq. (86), we have to perform inversion of the following integral Direct elementary integration (by substitution) gives where After some algebraic manipulation, we find that a general solution of the system for f = 0 is defined by In analogous way we find that Finally, one can go through to the real-time parametrization by integrating relation (85) which gives Let us now find an equation of orbits in (r, ϑ)-plane. Instead of using explicit formulae for u and v, it seems easier way to eliminate dτ from Eq. (86). By separating the variables, the orbit equation becomes which can be easily integrated to give Here, c is a constant of the integration. Taking the exponents of both sides, we find the implicit orbit equation Taking the inverse of transformation (18), namely we get an implicit equation of orbits in the polar coordinates as .
(104) Fig. 29 The phase-parametric diagram of super-integrable system (109) for Π = 1 at variations of the initial angle ϑ 0 . The initial conditions was chosen as (r 0 = 1, ϑ 0 = ϑ 0 ,ṙ 0 =θ 0 = 0). Here, ϑ(t ) denotes the extreme values of ϑ whenθ(t ) = 0 for some t . The diagram indicates the periodic orbits of the superintegrable system without any regimes responsible for quasiperiodic motion. See Fig. 30, which shows exemplary periodic orbits plotted in 3D Euclidean space angle ϑ 0 ∈ (0, π) of the pendulum. We obtain elegant curves corresponding to periodic motion of the system even for large values of ϑ 0 . Compared to the integrable system (see Fig. 10), there are no quasi-periodic orbits at all. Since Π = 0, the pendulum with mass m performs spacial oscillations. Some of the particular periodic orbits, plotted in 3D Euclidean space, are visible in Fig. 30. The initial conditions were taken from the phase-parametric diagram. As we can notice, each orbit is closed and we obtain 3D versions of the heart orbit and the loop orbit visible in Fig. 30a and in Fig. 30c. Let us summarize the results of this section in the following theorem.

Theorem 2
If M = 3m, s = 0 and 3k 1 = 4k 2 , then the generalisation of 3D-SAM governed by Hamiltonian (70) is completely integrable. In the absence of gravity, g = 0, the system is maximally superintegrable with three functionally independent first integrals. The integrable and the super-integrable systems have for their explicit forms defined in (108) and in (109).

Composition of the models
Finally, let us now combine both integrable models into one and let us study its integrability. Specifically, we consider the following Hamiltonian The question of what can be added to an integrable Hamiltonian function without destroying its integrability has a long history. KAM theory says that most of perturbations destroy integrability. However, it does not give any conditions for possible integrable perturbations. As a matter of fact, most modifications of an integrable potential gives rise a non-integrable one, see e.g. [35][36][37][38][39]. Surprisingly enough, the Hamiltonian (111) is integrable and an additional first integral is a quartic polynomial with respect to the momenta. An explicit form of the integrable system is as follows J 3D = J 2D + Π 2 sin(ϑ/2) r sin 2 ϑ 2 + Π 2 p 2 r 16 sin 2 (ϑ/2) + 2gr + 4kr 2 , where J 2D = r −1 sin ϑ 2 p 2 ϑ − 1 2 cos ϑ 2 p r p ϑ − gr 2 cos ϑ 2 sin ϑ −γ sin ϑ 2 − 2kr 3 cos ϑ 2 sin ϑ, is the first integral of the system at zero level of p ϕ = Π = 0. This is the new, not trivial, integrable model with physical origin.

Conclusions and open problems
This paper has discussed a fundamental problem of dynamical systems which is the distinction of integrable systems from non integrable ones. We focused our theoretical analysis on two generalized classes of 3D SAM system with additional Coulomb's and Hooke's interactions. For most values of the parameters, both models exhibit highly complex dynamics, which was presented via the Poincaré sections, the Fig. 30 Periodic solutions of super-integrable system (109) plotted in 3D Euclidean space for Π = 1. The presented orbits are 3-dimensional generalisations of the heart and the loop orbits visible in Fig. 30a and in Fig. 30c phase-parametric diagrams and the Lyapunov exponents spectrums. However, surprisingly, we have also found integrable cases: for the first model, the ratio of the masses M/m = 3 was sufficient for integrability, while for the second model, the additional condition for the ratio of the spring constant k 1 /k 2 = 4/3 was needed for integrability. We constructed the bifurcation diagrams for the integrable systems and we found explicit forms of general solutions of the equations of motion. These solutions are expressed in terms of the elliptic Jacobi functions. Moreover, we showed that in the absence of the gravitational potential, both models are maximally super-integrable This is an exceptional feature for Hamiltonian systems with more than one degree of freedom. All in all, the results presented in this manuscript were obtained with powerful tools, whose application seems to be of importance and usefulness to the study of other oscillatory systems. Moreover, the systems considered in this paper can be treated as unique theoretical models of Hamiltonian systems with two degrees of freedom which are chaotic, integrable, and superintegrable for certain values of the parameters. Therefore, such systems provide excellent examples for the teaching of Lagrangian and Hamiltonian mechanics.
Nevertheless, there is a number open problems related to the analysis of these models. For instance, the question of existence of other integrable cases for different values of the parameters remains open. For this purpose the Morales-Ramis theory can be adopted [40,41]. This theory enables a distinction between values of the parameters for which a considered model is suspected to be integrable. Morales-Ramis theory has already been successfully applied to various, important physical systems, see for instance [11,[42][43][44][45][46][47][48].
As second problem results from a lack of research on the dynamics and integrability the Atwood machine with two swinging masses. This model, which is the system of three degrees of freedom, was preliminary considered in [13]. However, its dynamics and its integrability analysis still does not attract sufficient attention and it is interesting whether this model possesses an additional first integral or even if it integrable for the mass ratio M/m = 3.
Finally, the physical construction of the models considered is now under way. We plan to obtain experimental results concerning their dynamics and compare these results with numerical and analytical results obtained in this paper. Therefore, our future work arguably will complete the above theoretical results. statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.