Trajectories of a ball moving inside a spherical cavity using first integrals of the governing nonlinear system

Analytical study of ball vibration absorber behavior is presented in the paper. The dynamics of trajectories of a heavy ball moving without slipping inside a spherical cavity are analyzed. Following our previous work, where a similar system was investigated through various numerical simulations, research of the dynamic properties of a sphere moving in a spherical cavity was carried out by methods of analytical dynamics. The strategy of analytical investigation enabled definition of a set of special and limit cases which designate individual domains of regular trajectories. In order to avoid any mutual interaction between the domains along a particular trajectory movement, energy dissipation at the contact of the ball and the cavity has been ignored, as has any kinematic excitation due to cavity movement. A governing system was derived using the Lagrangian formalism and complemented by appropriate non-holonomic constraints of the Pfaff type. The three first integrals are defined, enabling the evaluation of trajectory types with respect to system parameters, the initial amount of total energy, the angular momentum of the ball and its initial spin velocity. The neighborhoods of the limit trajectories and their dynamic stability are assessed. Limit and transition special cases are investigated along with their individual elements. The analytical means of investigation enabled the performance of broad parametric studies. Good agreement was found when comparing the results achieved by the analytical procedures in this paper with those obtained by means of numerical simulations, as they followed from the Lagrangian approach and the Appell–Gibbs function presented in previous papers.


Introduction
The motivation for this work comes from the results of two recent papers by the authors which examined the behavior of a ball-type absorber using numerical simulations [36,37]. In these engineering-oriented articles, a number of particular phenomena were identified, the causes of which had to be estimated from numerical results alone. However, these phenomena can have a major impact on the stability of the whole structure including the absorber. The presented analysis based on the first integrals, which is made possible by the assumption of a fixed cavity, enables a detailed analysis of the movement of the sphere inside the cavity. Therefore, this work does not intend to primarily examine the design details and effectiveness of the ball-type absorber. Instead, it focuses on the theoretical explanation and interpretation of phenomena observed during numerical and physical experiments, and on a clear separation of individual cases corresponding to the natural properties of the system.
The family of passive tuned mass vibration absorbers is well established in the engineering literature; see the exhaustive review paper [13] which reflects the situation until 2017. Conventional passive absorbers are of the pendulum type, mostly based on the autoparametric principle; see, for instance, [16,30,31,45]. The basic principle of ball-type absorbers comes from the rolling movement of a metallic ball of radius r inside a metallic spherical cavity of radius R > r (Fig. 1), possibly with a rubber lining. Such a simple device is practically maintenance free; this is a topic which is gaining increasing popularity [10]. Its vertical dimension can be relatively small, and thus, it can be used in cases where a pendulum absorber is inapplicable due to lack of vertical space or difficult maintenance. Not surprisingly, this type of absorber is growing in popularity in connection with wind turbines, e.g., [8].
The first papers dealing with the theory and practical aspects of ball absorbers were published during the past decades, see [41,42], and are based on engineering approaches. However, to the best of the authors' knowledge, the first papers that dealt with this issue from the perspective of rational (analytical) dynamics were published only some years ago; see [38,39]. They presented a basic nonlinear model in 2D together with its numerical evaluation, and a report on its practical application, including some results of long-term in situ measurements. The 2D approach is satisfactory when Fig. 1 Ball vibration absorber in a dynamic testing laboratory the absorber is limited to acting in a specific direction, e.g., in bridges, systems of multiple one-directionally acting elements, etc. However, structures like masts or towers require equipment which works simultaneously in both horizontal directions, and hence, a 3D model must be considered.
Modeling of a homogeneous sphere rolling on a perfectly rough surface has a long tradition in classical mechanics. Several 3D approaches are available which respect the spatial and strongly nonlinear character of the system. The system is non-holonomic with linear constraints in the first derivatives with respect to time. The classical setting of several particular cases, including the one that this paper is concerned with, is considered by Routh [43]. He and other authors of popular monographs use the general Lagrangian methodology [4,32,40], which is by far the most popular approach in classical mechanics; it offers many advantages which are actively applied in practice and research, including the rolling sphere problems [17]. For example, a study regarding the rolling of a ball over a spherical surface based on homogeneous Lagrange equations has recently been published [18]. There the author analyzes the free and undamped movement of a ball in the vertical plane using phase portraits for different initial conditions together with generalizations regarding vibro-impact dynamics.
The approach based on the Appell-Gibbs function, even though not frequently used, has proven very effective in subsequent numerical simulations. For details and some specific attributes of this approach, see, e.g., papers [11,14,26,28,46,48]. Recently, the rolling of a ball over a curved surface was dealt with on an abstract basis using the Lie group theoretical methods, e.g., [19]. Borisov et al. [7] use a similar abstract methodology and present a more complete analysis of the Routh solution for the solid of revolution. They derived new integrals for the ball rolling on nonsymmetrical surfaces of the second order. Jurdjevic and Zimmerman [20] extended the problem to a hyperbolic analogue in which the spheres are replaced by the hyperboloids, and rolling is taken in an isometric sense in either Euclidean or Riemann geometry. A numerical analysis of the ball rolling in a spherical recess, also based on the Appell-Gibbs function, was studied numerically by Legeza [24].
Another derived topic with high publication activity regards the rolling of a so-called Chaplygin sphere-a dynamically non-symmetric non-homogeneous spherewhich represents one of the best known integrable systems of classical non-holonomic mechanics [5], either with or without considering a spin [22]. There is no doubt that devices based on similar effects find their use in absorbing unwanted vibrations. The usage of nonhomogeneous spheres, hemispheres or semi-elliptic spheres would allow the absorber to be fine-tuned for a precisely limited nonlinear damping effect or multidirectional damping [27,28]. These topics are already popular in the engineering literature, however, still only partially treated analytically. For example, a prospective tuned vibration absorber based on the nested ball principle has been vaguely described in a patent proposal [33]. The theoretical analysis of this setup is partially addressed in an abstract way by Borisov [6], and a case with a semi-elliptic cavity was described by Legeza [25]. It is also worth noting that the overview [13] does not mention any text concerning the dynamic stability analysis of structures equipped with vibration absorbers. This is despite the fact that the subject is mentioned among the topics that require considerable attention.
The above-mentioned abstract solutions offer valuable tools of investigation. From the engineering perspective, however, the actual trajectories that can be encountered in a particular vibration absorbing device are interesting because they determine its efficiency. The authors addressed two possible strategies for this purpose: (1) the Lagrangian formalism in 2D [38] and (2) the Appell-Gibbs function in 3D; see [36,37]. In each case, a governing differential system was composed, and the solution itself was conducted numerically with subsequent analysis of the extensive data set. Important physical properties of the "ball-spherical cavity" system were discovered in these papers, and many particular trajectory types were identified numerically, however, without proper analytical explanation. This indicates that the problem should also be investigated at an analytical level using energy, momentum and angular momentum balance principles related to relevant first integrals. This strategy makes a qualitative investigation of the system behavior possible at least in the setting of the homogeneous problem. A similar technique was also adopted for investigating a nonlinear spherical pendulum; see [35]. In general, such an approach can significantly improve insight into the internal character of the system and increase the possibilities for practical applications. It enables a systematic identification of limit trajectories and the definition of relevant categories. The existence of limit trajectories depends on the system parameters and the initial conditions setting. The possibility of delimiting the domains of parameters can significantly contribute to an analysis of a system's reliability and its lifetime period. Moreover, the obtained results can serve in constructing effective forms of the Lyapunov function intended for particular cases; see [9,23] and other resources in the domain of dynamic stability.
Finally, it is worth mentioning that many of the algebraic manipulations required to derive or verify formulas in this paper were done using the Wolfram Mathematica Package [47].
The paper is organized as follows. First, after this introduction, the governing system and three of the first integrals based on the Lagrangian approach are derived. A particular characteristic function is then introduced as a basic tool for further classification of the response trajectories. Next, the particular "separation circle" trajectory is introduced, which separates two main trajectory groups. Sections 4 and 5 describe settings with and without consideration of the initial spin of the ball. A particular type of an almost planar rocking is then analyzed in detail in Sect. 6. Finally, the last section concludes.

Lagrangian system and first integrals
The rolling without slipping of a ball on a surface is a non-holonomic problem because constraints relat-ing to the mutual movement of a ball and a surface include velocity components. When putting together an expression for the kinetic and potential energies T , V , and external forces Q = [Q 1 , . . . , Q n ], the relevant Lagrangian equations should be written as follows: . . , n, and with non-holonomic constraints: where q j ( j = 1, . . . , n) are generalized coordinates. Symbols B m j and B m are generally functions of q j . Explicit time can be usually omitted as constraints are considered scleronomous as a rule if external kinematic excitation is absent and only initial conditions provide energy to the system. Provided that kinematic excitation works, the respective parameters B m = 0 and should then be considered as functions of time. Symbols λ m (m = 1, . . . , l) are Lagrangian multipliers which are used to add non-holonomic conditions to the Hamiltonian functional. Therefore, the system Eqs.
(1), (2) includes n + l unknowns q j , λ m ; see, for instance, the popular monographs [2,3,12,15]. The non-holonomic constraints Eq. (2) are formulated as linear functions of velocitiesq j ; this has been shown to be satisfactory concerning the problems considered. For more details and generalization, see [4,44]. In order to arithmetize the mathematical model, we need to introduce three adequate coordinate systems. In accordance with Fig. 2, the fixed Cartesian coordinates (x, y, z) are obvious. Their origin is in the "Southern Pole of the Cavity" (SPC) denoted A. The position of the ball center is described in standard spherical coordinates with their origin being in the center of the cavity and α, γ denoting the polar and azimuthal angles, respectively. The origin of the moving coordinates is located in the center of the moving ball, so that it lies on the concentric sphere with radius = R−r . Moving axis p follows a tangent of the concentric sphere meridian in the vertical plane (x , z), axis q is always horizontal, and axis n has the direction of the upward directed normal at the contact of both bodies; see Fig. 2. Rotation of the ball with respect to the moving coordinates is represented by the Euler angles ϕ, θ, ψ. Components of the angular velocity vector ω = [ω p , ω q , ω n ] T in moving coordinates are positive as corresponds to the usual convention of the right hand rule [40]. The velocities of the ball center with respect to global coordinates In further text, only the components of vectors v, ω will be used. Nevertheless, their relation to the velocity componentsq j used in Eqs. (1), (2) is obvious.
The basic formulae for kinetic and potential energies with respect to moving coordinates read where m, g represent the mass of the ball and gravitational acceleration, respectively. In order to relate the angular velocity vector ω and the velocity components expressed in Euler angles ϕ, θ, ψ, we write Because the spherical cavity has a constant curvature 1/ at every point regardless of direction, the relation between angles α, γ and velocities v p , v q can be expressed simply, see Fig. 2, where the expression v n = 0 represents one of the three contact constraints. For the same reasons, the Pfaff contact conditions of perfect rolling without slipping can be easily expressed. They specify the relations between angular velocities ω p , ω q , ω n and angles α, γ or displacements v p , v q . With respect to Eq. (5), the contact conditions can be reformulated as follows: Taking into account Eqs. (5) and (6), the expressions for energies from Eq. (3) can be rewritten in the form With reference to the problem definition, no external excitation or energy dissipation is assumed. Hence, the internal energy introduced to the system by nonhomogeneous initial conditions has the form: Here, the spin of the moving sphere is given in the form ω n =ψ +γ cos α with respect to Eq. (4c) and the geometric properties of the cavity. Every conservative Lagrangian system (in the sense of energy balance) possesses at least one first integral, which can be considered a multiple of the total energy of the system. Therefore, with respect to Eq. (8), it can be formulated as follows: (9) where the following notations were adopted: The dimensionless constant μ reflects the ratio between the ball radius and the distance between the centers of the ball and the cavity; it is, in fact, closely related to the ratio of both radii. ω 2 0 denotes the squared natural frequency of the rocking of the homogeneous ball in the spherical cavity, and E relates to the internal energy of the system.
Due to the transparent structure of the problem, although it is non-holonomic, it can be rewritten in three degrees of freedom only, and no procedure via Lagrangian multipliers has to be applied. Moreover, since no explicit external excitation is present, both terms on the right-hand side of Eq. (1) identically vanish. Therefore, after some modifications, we obtain three Lagrangian equations for the three unknowns α, γ and ω n : It is obvious that Eqs. (11b,c) are the first integrals, because γ and ω n are cyclic coordinates. Therefore, we can write Parameters E 0 and H 0 represent the energy and angular momentum of the system, respectively; they are given by Eqs. (8,12a) or (9, 12a), where the initial conditions α,α,γ , ω n are substituted. Equations (12) are the second and third first integrals of the system. The first one represents the conservation of the system's angular momentum with respect to axis z at the constant level H , while Eq. (12b) shows the spin velocity (proportional to S) of the ball with respect to the normal n. We can see that the spin velocity is constant throughout the whole period of the system movement, although it interacts with the other angular velocity components of the ball through relations Eq. (4). Note that this very special character of the spin is a direct consequence of the spherical shape of both the cavity and the ball. Any other combination of two bodies in contact would lead to a variable spin velocity. The general form of Eq. (11c) would encompass an additional term dependent on differences in the principal curvatures of the two bodies. This difference is zero for spherical surfaces.

Characteristic equation
The three first integrals Eqs. (9,12) represent certain invariants and an alternative description of the system behavior; they provide a possibility for a transparent introduction of energy and movement through initial parameters. Thus, they enable the investigation of particular states of the system, an analytical formulation of various characteristics of trajectories and, consequently, much greater insight into the nature of system parameters than one based solely on a numerical integration of the original differential system.
Let δ denote the height of the ball above the bottom of the cavity, i.e., the z coordinate of the center of the sphere: Ifγ is eliminated from Eq. (9) using Eq. (12a), one obtains after some manipulatioṅ In further text, we term f (δ) defined by Eq. (14) the "Characteristic Function" (CF), and f (δ) = 0 the "Characteristic Equation" (CE). Symbols E and H in Eq. (14) represent the measure of energy introduced into the ball at the instant t = 0 by means of the initial conditions: δ c -the initial height of the ball,γ c -"Initial Horizontal Velocity" (IHV), and ω n -"Initial Spin Velocity" (ISV); see Eqs. (9, 12a). Based on these initial conditions, the energy in the system may be quantified using the following relations: Function f (δ) is a cubic parabola which attains positive or negative values based on system parameters μ, ω 0 , state variables α,α,γ , ω n together with their initial values α c ,α c ,γ 0 , ω n0 (hidden in E, H ) and with respect to independent variable δ. The general form of the CF is obvious in Fig. 3. The interval δ ∈ (0, 2) spans the whole diameter 2R of the cavity from the SPC to the NPC (NPC-"Northern Pole of the Cavity").
Obviously, it holds that The cubic polynomial Eq. (14a) is physically meaningful only for values where f (δ) > 0, asδ is considered to be real. The first derivative of f (δ) may be formally written as a general quadratic polynomial: Quadratic equation Eq. (17) has two real roots because its discriminant is always positive: The inequality Eq. (18) is fulfilled trivially for ω n H > 0. Otherwise, introduction of initial conditions δ c ,γ 0 and ω n into E, H , defined by Eqs. (9, 12a), implies that Eq. (18) is valid for 0 ≤ δ c < 2. Therefore, the cubic parabola Eq. (14a) has two extremes. The analogous procedure confirms that C > 0. Thus, and both the extremes are situated on the positive semiaxis: 0 ≤ δ e1 ≤ δ e2 ; the first one is positive, and the second is negative: Summarizing the above contemplation, one can conclude that the CE Eq. (14a) has three real roots satisfying the following conditions: The first two roots are physically meaningful, as they delimit an interval on axis δ whereδ 2 ≥ 0, which is a necessary condition for the energy accumulated in the system to be real. For geometrical reasons, the values δ > δ 3 > 2 do not represent a physically meaningful state, althoughδ 2 ≥ 0 there as well. Note that zero and the coinciding roots can occur. As we will see later, they represent physically important cases.

Definition and relevance of the Separation Circle
Previous studies by the authors [36,37] have shown that the most important separation limit between trajectory types (or groups) that start at a certain point is a trajectory running at constant angular velocity Γ along a parallel of the cavity. This trajectory will be called the "Separation Circle" (SC), because it separates qualitatively different trajectories. The SC can be characterized as follows: The ball is pulled up along a meridian of the cavity to a certain level δ c to the "Starting Point of the Trajectory" (SPT), and subsequently, a horizontal impulse is applied to it. The intensity of the impulse is so high that it sets the ball onto a horizontal circular trajectory at the vertical level δ c . This condition can be used reversely to determine the necessary IHVγ c = Γ or v qc , the values of which represent the constant angular or tangential velocities relevant to movement along the SC. The initial spin velocity ω n is considered zero in the first step. The spin velocity can also be set such that it compensates forγ c when different from Γ in order to maintain the SC trajectory. Let us evaluate the notion of the SC from the viewpoint of the CE discussed above. Figure 4a shows parabola f (δ) for a given height δ c and three different values of the IHV. Intervals where f (δ) ≥ 0 for δ ∈ 0, 2 represent possible heights at which the trajectory of the ball in the cavity can occur. In the case of the SC trajectory, i.e., for IHVγ c = Γ , one double root δ 1 = δ 2 = δ c of f (δ) = 0 occurs; this case provides an active area of zero width, see the dotted curve c in Fig. 4. The dependence of the width of the active area on the initial horizontal velocityγ c is illustrated in Fig. 4b.
Forγ c < Γ , the active area spans between δ 1 and δ 2 which coincide with δ c . An initial velocity higher than Γ leads to a trajectory within the spherical strip above the δ c = δ 1 boundary and goes up to δ 2 . In such a case, it can occur that δ 2 > 1, which means that the ball passes into the upper hemisphere of the cavity. The limit case is reached when the IHV approaches an infinite value. Then, δ 1 , δ 2 are symmetrically distributed with respect to δ = 1. The upper boundary of the active strip is represented by the SC mirrored in the upper hemisphere. The trajectory becomes planar again although the plane is slanted passing the SPT and the center of the cavity.
The graph in Fig. 4b also shows the effect of nonzero initial spin ω n . The dashed curve 2 shows the case when ω n < 0, and the dot-dashed curve 3 corresponds to positive ω n . Each of the curves also has a second branch above δ = 2 that have no physical meaning.
The classification strategy based on the SC was intuitively adopted by the authors in their earlier studies, e.g., [36,37]. Actually it appears that this classification well describes all possible trajectories of a ball rolling inside a spherical cavity and starting from a certain point. Whatever the orientation and intensity of the initial impulse are, the movement of the ball takes place b Schema of the dependence of δ on the initial horizontal angular velocityγ c /Γ for three values of ω n : solid curve 1 : ω n = 0, dashed curve 2 : ω n < 0, dot-dashed curve 3 : ω n > 0. Verticals ac correspond to curves ac in plot a) within a uniquely defined spherical strip delimited by the two lower roots δ 1 , δ 2 of the CE, Eq. (14), which are related to the energy contained in the system.

Position of the separation circle
Let us inspect Eq. (11a). Provided the ball follows the SC at a height given by initial condition α c < π/2, i.e., δ c < 1, its angular velocity is constant,γ 0 = Γ , and vertical accelerationα vanishes, so that we can write This quadratic equation has two real roots: When zero spin velocity is assumed, ω n = 0, these relations simplify greatly: The IHV v qc or angular velocity Γ satisfying Eq. (23) produces the SC for the given polar angle α c . The roots represent two opposite directions with respect to the trajectory initial point. Their ratio depends on the sign of the spin velocity ω n . For zero initial spin, the image is symmetrical in the horizontal plane with respect to the initial point of the trajectory. The schematic plots in Fig. 5 demonstrate the dependence of velocity v qc on both the initial height (given by parameter α c ) and the ratio of the sphere to the cavity. Graph (a): fixed height α c = π/4, plot (b): fixed ratio r/R. In this latter case, the dependence starts from zero for the SPC (α c = 0) and tends to infinity for the "Equator of the Cavity" (EQC) (α c = π/2). The bold-black curves in plots (a) and (b) represent the spin-free state (ω n = 0), while the color curves show the influence of the initial spin. The relation between the velocity v qc and the initial spin ω n , which maintains the trajectory in the SC, may be deduced from Eq. (22) and is illustrated in picture (c). A decrease in horizontal velocity is compensated by a negative spin, whereas an increase in horizontal velocity implies a proportional increase in positive spin. One-sided limits of ω n for v qc → 0 + and v qc → ∞ exist and equal ±∞.

Dynamic stability of the separation circle
We now examine the neighborhood near the SC. We revisit Eq. (11a) and also the first integral Eq. (12a). The vertical position of the ball on the SC is given by the angle α c , and the horizontal angular velocity Γ is determined by Eq. (23). These are the nominal values which are subjected to small perturbations generated by dispersion in the initial conditions setting. Thus, we reformulate the initial conditions aṡ where ζ and η are small values. Substituting perturbed initial conditions Eq. (25) into Eq. (11a), a linearized equation for ζ(t) can be deduced. Disregarding the higher-order terms ζ 2 , η 2 , η· ζ , one obtains The first line of Eq. (26) vanishes due to Eq. (11a), and the coefficient of the term with cos α c on the second line disappears because of Eq. (22). Hence, it holds thaẗ This equation is solvable for zero initial conditions in a closed form: At the level of the first approximation, supposing that a small increase in initial tangential velocity is considered, it is obvious that the trajectory lies above the SC within the narrow spherical strip δ ∈ (δ c = δ 1 , δ 2 ), where δ 1 < δ 2 . Similarly, decreasing the IHV, we get a trajectory below the SC within the limits 0 < δ 1 < δ 2 = δ c . The width of the strip in both cases is |η|/Ω v .
The explicit form of perturbation Eq. (28) represents a form of harmonic ripples which oscillate above or below the separation circle. The period of perturbations is One loop around the SC takes T loop ≈ 2π/(Γ + η). Therefore, the number of small waves during one loop is which, in general, is not a rational number, and the trajectory alternating above or below the SC does not pass the SPT. We should be aware that the estimates Eq. (25) and also the results Eq. (28) are applicable if 0 < α c < π/2. Indeed, for small α c , the values α c and ζ are commeasurable as are the values Γ, η, and classification according to the powers of a small parameter becomes invalid. Of course, the perturbations η, ζ should also remain small in order for linear approximation to be justified.
Finally, one can conclude that the SC (perhaps with the exception of the SPC neighborhood) is dynamically stable, and a small initial perturbation does not cause a receding of this trajectory from the original SC.
Let us show an example of a time history of the trajectory corresponding to the SC; see

Roots of the CE
Let us examine the case where δ c = δ 1 ≤ δ 2 and ω n = 0, i.e., the ball is rolling above the SC and no spin is assumed,γ 0 > Γ 0 ; see also curve (a) in Fig. 4. Setting the position for the SPT on a meridian of the cavity at a certain level characterized by polar angle 0 < α c < π/2 or equivalently by height δ 1 ∈ (0, 1) (in the lower hemisphere of the cavity) in fact means that the lowest root δ 1 = δ c of the CE is fixed. Considering that one root is known, we can rewrite Eq. (14a) in the form of a partial decomposition with respect to the root factors: The detailed form of coefficients K , L , M is derived from the full CE, Eq. (14a), where a zero ISV was substituted. The energy E, Eq. (9), and the angular momentum H , Eq. (12a), contained in Eq. (14a), were included in the initial parameter values for the SPT.
The SC is regarded as the lower boundary of the spherical strip on the cavity surface within which a par-ticular trajectory runs. The upper boundary is then the lower of the remaining roots δ 2 , δ 3 . They can be calculated from the quadratic equation, which is outlined in Eq. (31): The discriminant D is always positive, Eq. (32c), which fact confirms that Eq. (14a) has three real roots. Furthermore, 0 < D < L 2 and, consequently, all roots are positive and fulfill conditions 0 < δ c = δ 1 < δ 2 < 2 and δ 3 > 2. Of course, root δ 3 is geometrically out of scope, and it is no longer considered; see Sect. 2.2.

Circumferential periodicity of trajectories
In order to outline the basic character of a trajectory occurring between boundaries δ c = δ 1 , δ 2 , we inspect the expression for the circumferential velocityγ following from Eq. (14b). In general, it can be assumed that the trajectory is periodical in the vertical direction, where the ratio of the period to the SC lengths is not a rational number. Therefore, one round along the SC will not contain a whole number of periods. The angular momentum H is always positive, see Eq. (12a), and the second term of the numerator vanishes since ω n = 0. The denominator in this fraction is also positive, because δ(2 − δ) > 0 for δ ∈ (0, 2). Therefore,γ > 0 regardless of the settings for the initial conditions at the SPT. Moreover, it can be supposed that the variability ofγ during one period for given initial settings will not be dramatic. Indeed, it is obvious that the angular momentum H is constant during one period. Consequently, the horizontal angular velocity at the point where the trajectory touches the upper boundary of the strip, δ 2 , follows from Eq. (12), where ω n = 0, and it holds thaṫ where α c ,γ 0 or α 2 ,γ 2 are values of the respective parameters at the initial point (δ c = δ 1 ) or at the touching point on the upper boundary (δ 2 ). Hence, it can be written: The horizontal velocity is slightly lower at the upper apex due to an increase in potential energy. At the same time, some qualitative confirmation of this fact follows from the constant angular momentum H and δ, which changes more or less monotonously. This implies that no singular points emerge during one vertical period, and that the angular velocity along the trajectory is mildly variable. Thus, the trajectory has the shape of a simple periodic curve without any turnabout points reversing velocityγ .

Vertical periodicity of trajectories
Here, we assess the length and duration of one vertical period of the trajectory. The IHV leads to a trajectory that is symmetrical with respect to both points where it touches the lower boundary (δ c = δ 1 , where the IHV is applied) and the upper boundary (δ 2 ); see Fig. 7. Thus, it is sufficient to examine only the first half of the period for δ increasing between δ 1 and δ 2 . We revisit both relations in Eq. (14). Assuming ω n = 0, the following differential system can be written: The first equation isγ independent, and, therefore, it can be solved as the first step. The pair ±δ could be put into the second equation to obtain γ (t) by means of integration between δ c and δ 2 . However, it comes to light that both points on the strip boundaries δ c , δ 2 are singular and represent bifurcation points. In addition, the relevant Jacobi matrix is also singular and, consequently, does not enable us to predict the principal directions in the point neighborhood. Three solutions start from the SPT (including the constant δ c ), and all of them have the zero derivative. Therefore, neither an analytical nor a numerical stable solution can be deduced from these points. This difficulty can be overcome by differentiating Eq. (35a) with respect to time.
After reducing byδ, the modified system reads The denominator in Eq. (36b) is always positive because δ ∈ (δ c , δ 2 ). Then, Eq. (36) can be solved for initial conditions γ (0) = 0, δ(0) = δ c ,δ(0) = 0 sequentially putting partial results of Eq. (36a) into Eq. (36b). Finally, the time is eliminated, and one obtains δ as a function of γ . Some samples of the trajectory are plotted in Fig. 7, in which a set of five cases Notice the shape of the high curves in this figure. They correspond to the dominating first term in the right-hand side of Eq. (36a) whenγ 0 is large. Omitting the second term whenγ 0 Γ 0 , we obtain a linear equation of harmonic oscillation with an adequate constant right-hand side. The resulting equation forγ 0 → ∞ is analytically resolvable; its complicated resulting expression δ(γ ) is 2π periodic with the extremal values δ c and 2 − δ c at γ = 2kπ and 2(k + 1)π, k ∈ Z, respectively. The shape of the limiting function is indicated by the dashed black curve segments. This limit case ofγ 0 → ∞ will be discussed separately in Sect. 4.2. On the other hand, γ changes in steps from 0 to π whenever δ(t) passes the SPC foṙ γ 0 = 0, which makes γ T = π . In a general case, for given initial height δ c , the vertical component of the trajectory δ(γ ) increases its angular period and amplitude as the initial angular velocitẏ γ 0 increases. This can be observed in Fig. 8; the blue curves above the SPT,γ 0 > Γ 0 are regarded in this section. For a fixed Γ 0 and increasing IHV, the angular period γ T slowly rises and approaches the full angle foṙ γ 0 → ∞; see picture (a). Picture (b) demonstrates the rising of the strip's upper level δ 2 with increasing IHV. The parameter δ 2 approaches the horizontal asymptote 2 − δ c , symmetrically placed with respect to the EQC, in the particular case from 0.3 to 1.7. When δ c changes its value, the dependences shown in Fig. 8 remain valid, but the strip (δ c , 2 − δ c ) decreases/increases its width.

Discussion
Using the above analytical results, we can outline some detailed trajectory properties with respect to IHV and the height of the SPT above the SPC. A typical trajectory shape is plotted in Fig. 9. Comparing Figs. 6 and 9, we can see that the trajectory time history is still a simple periodic curve synchronous in both horizontal coordinates with a slight modulation. This modulation depends on both values of δ c andγ 0 . It relates to the difference between frequencies of γ (t) and δ(t) and vanishes whenγ 0 = Γ 0 orγ 0 → ∞.
The frequency of vertical displacement u cz in Fig. 9a is related to that of horizontal displacements via the Example of the trajectory above the SC with no ISV (ω n = 0): a time history, b top view, c axonometric demonstration angular period γ T which depends on the δ 1 = δ c and δ 2 boundaries of the spherical strip. No visible influence of higher harmonics is observed, although a very light quasi-periodic character can be noticed, as it follows from the geometric character of the system. The shape in the top view is helical and resembles the form of a prolate-type hypotrochoid close to a hypotrochoid form of prolate types. As no sharp apexes or loop multiple points are detected, we can conclude that no basic cycloid or curtate trochoid is approached.

General shape of the trajectory
Let us discuss the trajectories above the SC that emerge when the IHV distinctly exceeds the velocity Γ 0 or v qc , and simultaneously, no ISV is applied, ω n = 0. The CF for an infinite IHV has a form corresponding to Fig. 3, which means that zero points δ 1 , δ 2 are symmetrically distributed with respect to the point δ = 1. A lower IHV shifts δ 2 to the left, and the position of δ 1 remains the same. The trajectories are again concentrated within the spherical strip delimited by roots 0 < δ 1 = δ c < δ 2 < 2 of the characteristic equation Eq. (14a). We will inspect its evolution whenγ 0 Γ 0 . While it is always true that δ 1 < 1, the upper boundary, being given by δ 2 , can enter the upper hemisphere of the cavity reaching a value in the interval 1 < δ 2 < 2. The IHV corresponding to the transition case δ 2 = 1 follows from equation Eq. (14a), where the two lowest roots δ 1 , δ 2 are con-sidered as known. After some manipulation, one can writė where both the boundary values of α c leading to an infinite IHV are obviously not admissible, as could be expected.
Increasing the IHV beyond all limits, the upper limit of the strip approaches a theoretical maximum: It is an asymptotic position, which is monotonously approached as the initial velocity rises to infinity. In this theoretical state, the trajectory becomes planar, having a circular form with the diameter 2R. This plane is inclined, being determined by the horizontal tangent at the SPT and by the cavity center; see Fig. 11a.

Rotation of the osculating plane
In the case whenγ 0 is high but finite, 0 < Γ 0 |γ 0 | < ∞, the root δ 2 is adequately lower: Trajectories maintain their spatial character, but it is worthwhile to define an osculating plane at each point. The osculating plane makes sense from a physical perspective if it enables us to characterize the basic position of the trajectory using simpler elements than in the general case for lower IHV. In the case we are discussing, it can be assumed that the spiral does not differ much from a planar shape, and that its projection into the osculating plane at each point represents a good approximation. In other words, we can define an affine space which meets a sub-manifold at a point in such a way as to have a second order of contact at that point. The osculating plane passes the initial point δ 1 and a point in the upper hemisphere at the upper boundary of the strip δ 2 − , which is slightly below the level of the limit case Eq. (38). Therefore, its inclination is slightly lower than that corresponding to the limit case forγ 0 → |∞|. The osculating plane rotates around the vertical axis of the cavity with rotational speed Ω v , which decreases asγ 0 rises, and vanishes for an infinite IHV (γ 0 ). The velocity Ω v can be estimated based on the position of the contact point of the trajectory and the upper boundary δ 2 evaluated at the conclusion of a single period. Although an explicit formula cannot be expressed, the process can be carried out observing the outline in Fig. 10. A sketch of the osculating plane behavior is graphically demonstrated in Fig. 11  Assuming that velocityγ 0 is high, the second term of the right-hand side in Eq. (36a) can be neglected and also enables E to be reduced. Then, it can be approximately written: where coefficient E γ can be determined due to the equality of values at the lower and upper limits, The length of the period in time is and, therefore, for the rotational speed of the osculation plane around the z axis, we can approximately write These results can be confirmed intuitively by examining both of the graphs in Fig. 8. Indeed, the length of the period approaches the horizontal asymptote on the 2π level more or less exponentially. Because the tangential velocity along the trajectory is approximately constant,γ (t) ≈γ 0 +η, the relation between γ (t) and t can be expressed simply as γ (t) ≈ (γ 0 +η)t. Then, it holds approximately that where κ, [s] is a positive constant. Employing L'Hospital's rule, it is obvious that limγ 0 →|∞| Ω v = 0 holds as before.

Discussion
Let us point out some properties of the trajectory discussed above, such as features of a curve passing through a layer of thickness ; see Eq. (39) and Fig. 10. An analysis using the small parameter approach cannot be done directly because the solution forγ 0 → |∞|, which serves as a zero approximation, does not exist in an explicit form for several reasons (infinite energy, infinite angular momentum, indefinite derivatives with respect to time, etc.). Despite this fact, we have seen that all trajectories are stable, whatever their parameter setting is. Consequently, the existence of the zero approximation can be assumed in an implicit meaning of a certain limit. Therefore, analogously with Eq. we are entitled to writė where α 2 ,γ 0 are relevant to the planar trajectory incident with the inclined plane, see Fig. 11a, andη, ζ are small unidirectional deviations. They enable us to define high values ofγ , α but still finite values of the initial approximations α 2 ,γ 0 . We recall Eq. (11). By the way, let us note that despite the fact that the ISV is not included in this section (ω n = 0), we can see that the proportion of the spin energy at high velocityγ is negligible anyway. Thus, putting approximations Eq. (45) into Eqs. (11a,b), and comparing terms that involve the same powers ofη, ζ , we can writė It is obvious that for high values ofγ 0 , the term ω 2 0 cos α 2 is negligible. The remaining coefficient is always positive. In a ratio to the length of the circular trajectory, it represents a parameter analogous with Eq. (30). It is proportional to the velocity Ω v of the osculating plane rotation around the z axis, which approaches zero forγ 2 → |∞|. This result is identical with the one obtained above in this section.
To demonstrate the character of the upper root δ 3 during the IHV limitation to ±∞, let us assess its value respecting Eq. (39), and the fact that the trajectory is running in a thin layer following the scheme in Fig. 10. The process is limited to within this domain, and, consequently, it can be linearized in the framework of this layer. Making use of these facts and a factored form of the polynomial, we can reformulate the CE, Eq. (14a), because two roots δ c , δ 2 are known: which results in the third root: This formula demonstrates that 2 < δ 3 → ∞ forγ 2 0 → ±∞.
For completeness, let us demonstrate an example of a trajectory based on a high IHV (Fig. 12). The time

Roots of the CF
The CF has a form corresponding to curve (b) in Fig. 4. In general, δ 1 can descend to zero, as we can easily deduce from Eq. (50) or from the original characteristic equation [Eq. (14)], with the vanishing absolute term. This case, together with the neighborhood of this value (0 ≤ δ 1 < ), will be discussed in a separate section (Sect. 6). For these initial settings, the spin-free and the spin-considered cases of initial settings intermingle, and, hence, it is worthwhile to discuss the two of them together. The spherical strip in which the trajectory for IHV γ 0 < Γ 0 emerges is below the SC. The strip is limited by the SC, which forms its upper boundary, δ c = δ 2 , and by the lower boundary δ 1 , which is given by the quadratic term in Eq. (31). The basic analysis is similar to that which has been performed in the beginning of Sect. 4.1. The only difference is that the roots are ordered as follows: 0 < δ 1 ≤ δ c = δ 2 < 1 < δ 3 . The roots δ 1,3 may be symbolically written as in Eq. (32): where the discriminant has the same form as in Eq. (32b) for δ c = δ 2 .
In order to determine velocityγ at the tangent point at the δ 1 boundary, we refer to Eq. (34), where α 1 is to be substituted instead of α 2 : It is obvious thatγ 1T >γ 0 due to a lower potential energy at the δ 1 level. Inspecting Eq. (14), where ω n = 0 is substituted, we can see thatγ is a simple continuous function of δ, which is integrable on any interval δ ∈ (a, b) with 0 < a, b < 2. This implies that no singular points emerge within one period, and the velocity along the trajectory is mildly variable. The trajectory has the shape of a simple curve without any multiple or turnabout points reversing its velocity.

Vertical periodicity of trajectories
As for the length of a single vertical period, a similar deduction can be made like in Sect. 4.1. However, care should be taken when the IHV approaches zero and the system Eqs. (36) becomes unstable or discontinuous in the neighborhood of γ = π/2. For details, see Sect. 6. Nevertheless, in a common case, like in Sect. 4.1, we can assume once again that one period consists of two identical halves symmetrically distributed around the tangent point on the lower boundary δ 1 , and, therefore, it is sufficient to examine only one half of the period. Hence, the differential system Eq. (36) can also be used here, except that the zero initial conditions are formulated for the upper strip boundary, δ(0) = δ c ,δ(0) = 0 and the position of the lower boundary δ 1 should either be evaluated using Eq. (50e) or during solution of system Eq. (36).
Let us follow Fig. 13 demonstrating several trajectories comparable in their initial conditions forγ 0 < Γ 0 . Recalling Eq. (36a): (15), we can see that unlike cases withγ 0 > Γ 0 , the first term of the right-hand side loses dominance, and a significant nonlinear character of the equation emerges. Retaining only the second term, we obtain an equation solvable using elliptic functions. This effect is obvious in Fig. 13, where the curves for decreasing IHV become more and more similar to an elliptic sinus. Returning to properties of the rational or irrational spiral form, see also Sect. 4.1, we can conclude that the vertical period length varies in the interval 2T γ ∈ (π, 2π) throughout all |γ 0 | ∈ (0, ∞), and, therefore, obviously no synchronization of the primary type can occur.

Discussion
We shall point out some properties of the length and shift of periods that are specific for the trajectories below the SC. Characteristics of the vertical period length and the lower strip boundary position δ 1 as a function of the IHV and SPT level are demonstrated in Fig. 14; see the solid orange curves below Γ 0 = 3.63. Both curves are smooth on the whole intervalγ 0 including the SPT. It is typical that the period 2γ T shortens to π forγ 0 → 0, as it also corresponds to the character of an elliptical sinus that characterizes the shape of the trajectory; see [1,21]. The lower boundary of the strip δ 1 evidently tends to zero, i.e., toward the SPC.
With reference to Sect. 6, we can see that the period forγ 0 → 0 does not approach 2π , as could be intuitively expected. However, even cases discussed in this section evince some attributes which are distinctly visible for low IHV values. It is obvious, for instance, that the shape of the relevant curve depends on the SPT level. The phenomenon of period "shift" becomes more prominent as the height of the SPT or the boundary δ 1 level increases. This effect is discussed in more detail for both zero and nonzero ISV for low IHV in Sect. 6.
An example of a trajectory below the SC is plotted in Fig. 15. The time history in this domain appears as a simple periodic curve without any higher harmonics and with a weak multiplicative modulation. The SPC always lies inside individual loops. They run round the SPC and never pass it, whatever the SPT and IHV are. The spatial character of the trajectory is obvious in plot (c) of Fig. 15. As in Sect. 4.1, no visible intervention of higher harmonics is observed, although a very light quasi-periodic character can be noticed. The shape in the top view is helical and close to a hypotrochoid form of a prolate type. As no sharp apexes or loop multiple points are detected, we can conclude that no basic cycloid or curtate hypotrochoid is approached.

Roots of the CE
We again revisit Eq. (14). This time the ISV will be considered nonzero, ω n = 0. The characteristic function Eq. (14a) can be decomposed as in Eq. (31); however, the coefficients L , M will be different: As in Sect. 4, either the lower or upper boundary of the strip is one of the roots δ i , i = 1, 2, which are determined by initial condition δ c , i.e., δ c = δ 1 or δ c = δ 2 for the strip situated above or below the SC, respectively. The remaining two roots δ 1 , δ 3 or δ 2 , δ 3 can be calculated from the quadratic term in Eq. (52a) as in Eqs. (32a) or (50): The expressions above can be reformulated with reference to Eqs. (9, 12a). Provided E, H are specified with respect to initial conditions, cf. Eq. (15), we obtain Discriminant D is positive with the exception of two cases when D = 0: The first option represents two singular cases when ω n = 0 and δ 1,2,3 ∈ {0, 2}, while the latter one describes a particular trajectory which occurs for a given negative spin (assuming positiveγ 0 ) and passes through the NPC. This case, however, generally requires IHVγ 0 = Γ 0 ; only for δ c = 2/3 doeṡ The other particular case occurs when M = 0, i.e., for Then, δ 1 = 0 and δ 3 = 2 +γ 0 (2 − δ c )/ω 2 0 > 2. In other cases, K > 0, M > 0, and thus, 0 < D < L 2 . Hence, both roots are real and positive and fulfill the relation: 0 ≤ δ 1 ≤ δ 2 < 2 < δ 3 . Here, δ c is either δ 1 or δ 2 , depending on what is considered given. For the rest of this section, we consider the IHV to always be equal to positive velocity Γ 0 , as it corresponds to the IHV of the SC.

General considerations
We can see that the trajectories influenced by a negative ISV, ω n < 0, lie within a spherical strip above the SC, with the bottom boundary δ c = δ 1 and the upper boundary δ 2 , which follows from Eq. (53). Note that the third first integral, Eqs. (11c) or (12b), specifies that ω n is constant throughout the whole investigated process.
First, we briefly outline the basic character preliminarily inspecting Fig. 19 and taking into consideration also Figs. 16, 17 and 18. In general, observing the time history of the horizontal and vertical components of the response, it is noticeable that each trajectory consists of two basic components (except for some higher marginal harmonics). They are independent in their frequencies, as the first one is related to the basic spin-free movement and the second follows from the spinning rotation of the ball. The latter one proves to be more or less distinct in its amplitude according to the width of the strip where the particular trajectory is operating. In other words, it is determined by the active area δ ∈ (δ c , δ 2 ). In general, the influence of the spin on the overall shape of the trajectory increases with the value of |ω n | and acquires a significant dominance for ω n above limit ω ns , given by Eq. (60), and in particular for 0 > ω ns ω n or ω n → −∞.
It follows from Eq. (14) that just three types of trajectories can be encountered when ω n = 0. Let us remember that the same equations, presented in Sect. 4.1, conclude that only one type of trajectory can exist within the strip if no spin of the ball is applied. The main reason for this alteration follows from the fact that the right-hand side of Eq. (14b) can have both positive or negative values, when ω n = 0 is considered. The three types of trajectories can be classified with respect to the parameters and initial conditions of the system, namely the ISV. The shapes of the trajectory types differ significantly in the neighborhood of the contact point on the upper boundary of the strip; see Figs. 16 and 17.

Wavy trajectories
Let us now discuss some specific details of the individual trajectory types. The general form of the first type is obvious from Fig. 19 (i). The trajectories reflect the ISV in interval ω n ∈ (0, ω ns ), where ω ns is the ISV of the separating case, Eq. (60). Roughly observed, these trajectories are not too far from those discussed in Sect. 4, although some influence of the response component caused by the spin is discernible. However, the basic form again resembles an irrational spiral with slightly distorted detailed periods. The difference in the basic Fig. 16 Shape of the trajectory above the SC for various initial spin velocities; colors of curves: ω n = 0-black, ω ns < ω n < 0-red, ω n -bold green ("kings crown" shape-separating case), ω n < ω ns -blue. (Color figure online)   Fig. 17 Shapes of trajectories in the neighborhood of the contact point on the upper boundary of the strip (δ 2 ); ω n = 0black, ω ns < ω n < 0-red, ω n -bold green ("kings crown" shape-separating case), ω n < ω ns -blue. The symbol Δγ (horizontal axis) means a local coordinate within one period or an increase/decrease of γ with respect to γ = γ T (position of the tangential point on the δ 2 boundary). (Color figure online) shape is rather quantitative, as the frequency of the spin is more or less related with that generated by the basic movement of the ball, and that is why it is hidden in the primary component influencing its amplitude. The trajectory shape is obvious from Fig. 16, where it is plotted as a function of angle γ (the two red curves relevant to ω ns < ω n < 0).
Details of the trajectory character near the contact point on the upper boundary are demonstrated in Fig. 17 (the two red curves). For the sake of a better visual comparison of individual trajectory behavior in the neighborhood of the contact point on the upper boundary δ 2 , all trajectory graphs have been shifted and concentrated around this point, which serves as the origin of local coordinate Δγ . The starting and finishing points of one period on the lower boundary δ c are denoted γ b1 , γ b2 . The width of the strip increases with descending ω n from zero until a maximum width is reached for ω n = ω ns ; see Figs. 17 and 18a, b. This corresponds to the total energy conservation principle.
The assessment of the length and duration of one period of the trajectory can be done analogously to Sect. 4.1. The reasoning which brought us to the differential system Eq. (36) is more or less the same, being based on the fact that contact points on the strip boundaries represent points with unavoidable singularity. Hence, with reference to Sect. 4.1, we can deduce the following modified system: where The denominator in Eq. (57b) is identical with that in Eq. (36b) and is positive in the considered interval. The system Eq. (57) is solved for initial conditions: δ(0) = δ c ,δ(0) = 0, and, eliminating the time, one obtains δ as the function of γ .
For the selected δ c = 0.3 and the relevant Γ 0 = 3.63, two samples for ω ns < ω n < 0 are plotted in Figs. 16 and 17. Together with the three diagrams in Fig. 18, we can evaluate the character of a single period. Pictures (a) and (b) show the strip width as the difference between relevant angles α and parameters δ, respectively. Picture (c) presents the spatial width of a single period as a function of ω n . It is obvious that a decrease in ISV (ω n < 0) leads to an increase in amplitude δ 2 or α 2 until a maximum is reached for ω n = ω ns . The length of the period simultaneously decreases. The trajectory for ω n ∈ (ω ns , 0) could only have a periodic character if the ratio 2πρ sin α c /2T γ is a rational number, where T γ denotes the half-period in the angular scale.

The limit case-the pointed trajectory
The second type of trajectory, see Fig. 19 (ii), occurs for the spin frequency ω n = ω ns , see the bold green curve in Figs. 16 and 17. It represents a limit case separating groups below and above ω ns . The trajectory shape in the axonometric view resembles a "kings crown." It contains a sharp apex in the midpoint of every period which touches the upper boundary δ 2 . The derivative with respect to the circumferential coordinate γ is discontinuous in this singular point, jumping from ∞ to −∞. All velocity components vanish, and both tangential acceleration components are discontinuous. A spin energy that is proportional to ω 2 n is retained. To determine the ISV needed to achieve this special type of trajectory, we take advantage of the fact that values of the first integrals for E and H , Eqs. (9, 12a), are constant throughout the entire time history. Indeed, at the SPT they can be expressed as follows: while in the sharp apex, at the δ 2 level, it holds that Evaluating the relevant equivalences E δ c = E δ 2 and H δ c = H δ 2 , after some manipulations one obtains Substituting Γ from Eq. (23) only provides real results in an unrealistic situation when δ c > 1. However, when Γ = Γ 0 , Eq. (24), Eq. (60) transforms to Relations Eqs. (59-61) are valid only when the apex occurs at δ 2 . This is characterized by an infinite curvature of the trajectory, i.e., when 0 = (γ 2 +δ 2 ) 3/2 . The necessary condition forγ = 0, given by Eqs. (9), (15), (24), is 0 < δ c ≤ 2/3, cf. also Eq. (55). It holds obviously that δ 2 > δ c , and that the upper boundary of the trajectory can reach into the upper hemisphere of the cavity. This effect occurs when When δ c is increased further, the apex may reach the NPC for δ c = 2/3; this case also nullifies the discriminant Eq. (54b). For even larger δ c , the apex ceases to exist. These remarks are only theoretical, because in the upper hemisphere the contact force becomes negative particularly at the apex point, where all the components of the velocity vector v identically vanish.

Curly trajectories
The third type of trajectories, see Fig. 19 (iii), emerges provided ω n < ω sn and possibly descending as ω n → (a) (b) (c) Fig. 19 Examples of trajectories with negative initial spin velocity ISV (ω n < 0): 1st row: ω ns < ISV < 0; 2nd row: ISV = ω ns ; 3rd row: ISV < ω ns . Column a time history, b top view, c axonometric demonstration −∞. The trajectory has a curly form making a loop every period; see Figs. 16 and 17. As before, one period starts and finishes at the same vertical level, which is given by the root δ c = δ 1 . The highest point, reached in half a period, lies on the upper circle boundary. It is given by the root δ 2 according to Eqs. (53, 54). Let us turn our attention to the monotonously descending width of the strip, Fig. 18a, b for ω n < ω ns depicted either with respect to angle α or to parameter δ.
The sharp apex-encountered in the previous typedelimiting the first and third types of trajectories prefigures a formation of two turnabout points γ v1 , γ v2 with tangents following their relevant meridians and with vanishing velocityγ ; see Figs. 16, 17 (blue curves).
The vertical position of points γ v1 , γ v2 above δ c can be found with respect to conservation of E and H values along the trajectory and due to the fact thatγ = 0 at turnabout points andα = 0 at the SPT. Indeed, using Eqs. (9,12), we can write the following relations for unknown valuesα v , δ v at turnabout points: Simple manipulations result iṅ which indicates that velocityα is positive or negative in γ v1 or γ v2 , respectively, i.e., the trajectory is rising or descending. The level δ v > δ c , but the difference δ v − δ c descends to zero as ω n → −∞; see Sect. 5.4. Both expressions remind us that ω n should exceed a certain limit in order for the formulae to be meaningful; otherwise, a curly form of the trajectory can exist. The horizontal position of points γ v1 , γ v2 can be determined using the same system Eq. (57) as in the case of the first and second trajectory types, where γ b1 = 0 is taken as a reference point. The position of the remaining points γ v2 , γ b2 follows from the symmetry of the second half-period. Evaluation of the period length is already obvious and represents the difference γ b2 − γ b1 . Note that the point γ b1 always precedes the point γ b2 on the advancing coordinate γ like in the case of the first and second types of trajectories, and, consequently, the difference γ b2 − γ b1 is always positive, but monotonously approaches zero for ω n → −∞; see Fig. 18c.
Points γ v1 , γ v2 delimit the boundaries of the upper part of the loop, which is defined in the interval γ ∈ (γ v1 , γ v2 ) in the framework of one trajectory period; see Fig. 18. The velocityγ in the upper part of the loop between the points γ v1 , γ v2 is of the opposite sign than that in the lower part of the loop. The width of interval γ ∈ (γ v1 , γ v2 ) starts from zero for ω n = ω ns , where the curly form arises, and it increases for descending ω n < ω ns . This width reaches the period length and then becomes larger than the distance between the starting and finishing points γ b2 − γ b1 of one period. Although both widths subsequently approach zero, the width of the loop becomes more and more dominant. In addition, the total (curvilinear) length of one loop significantly exceeds the simple distance γ b2 − γ b1 .
It is useful to define a certain average (or effective) velocity Ω av of the ball advancement along the horizontal SC; it is given by the time and length of one period. This ratio drops as ω n increases. Consequently, the average velocity Ω av decreases accordingly. It approaches zero for high values ω n ; see Sect. 5.4. Due to the arrangement of points γ b1 , γ b2 on the γ axis, it holds that Ω av > 0. This positive sign is maintained throughout the whole interval ω n < 0.

General considerations
This section will be devoted to trajectories resulting from the same IHV Γ 0 and a positive ISV. The trajectories lie in a spherical strip below the SC, with the given upper boundary δ c = δ 2 and lower boundary δ 1 , which is determined by means of Eq. (53), i = 1. Similarly as before, we go through the main properties of these trajectories inspecting Figs. 20, 21, 22 and 23 together with three examples of the trajectory with positive ω n , which are presented in Fig. 24.
Trajectories below the SC can also be classified into three types or, in other words, into two groups separated by the special limit case corresponding to a fixed frequency ω n = ω ns like in the previous section. In general, all trajectories apparently have a spiral form of the prolate hypotrochoid type, as we can see in the top view presented in Fig. 24 (pictures in column (b) for the three types). However, the shape of these spirals differs from those above the SC. The difference between particular types for ω n > 0 consists mainly in their relation to point A (SPC). Let us point out that the time history again has the distinct form of a two-component periodic process resulting from the movement of the ball Fig. 22 Shapes of trajectories in the neighborhood of the contact point on the lower boundary of the strip (δ 1 ); ω n = 0-black, 0 < ω n < ω ns -red, ω n = ω ns -bold green (discontinuousseparating case), ω ns < ω n -blue. The symbol Δγ (horizontal axis) means a local coordinate within one period or an increase/decrease of γ with respect to γ = γ T (the position of the tangential point on the δ 1 boundary). (Color figure online) and its rotation around the moving normal specified by the ISV (ω n ); see Fig. 24 column (a).
The shape of the trajectory in the neighborhood of the contact point on the lower boundary of the strip and the development of a curly trajectory are obvious in Figs. 21 and 20. The width of the strip δ c − δ 1 ∈ (0, 1) changes from zero (ω n = 0, δ 1 = δ c ) until reaching δ c (maximal width) when ω n = ω ns and δ 1 = 0. In this latter case, the trajectory passes through the SPC and represents the transition case. A further increase of ω n > ω ns leads to rising of δ 1 backwards to the SC, which is reached for an infinite ISV; the width of the strip vanishes.

Trajectories running around the SPC
Let us point out some important features of particular trajectory types. The first type of trajectory, see Figs. 20a, 22 (red curves) and 24 (i), is related to an ISV in the interval 0 < ω n < ω ns . The trajectory has a spiral shape, where individual loops are prolate and running around the SPC. The influence of the second periodic component is small but still discernible; see Fig. 24(i). This phenomenon becomes more pronounced as the limit case ω ns is approached.
Roughly observed (as in Sect. 5.2), the trajectories do not differ significantly from those discussed in Sect. 4.3. The basic form resembles again an irrational spiral with slightly distorted detailed periods. The trajectory shape is obvious from Fig. 21, where three red curves are plotted relevant to ω ns < ω n < 0 as functions of angle γ .
Details of the trajectory character near the contact point on the lower boundary δ 1 are demonstrated in Fig. 22 (three red curves) within a single period. This depiction in Fig. 22 was used in order to facilitate a comparison of the trajectory behavior throughout all ω n considered in the neighborhood of the contact point on the lower strip boundary δ 1 . The width of the strip increases with ascending 0 < ω n < ω ns from zero until a maximum is reached for ω n = ω ns as it corresponds to the principle of conservation of total energy.
The vertical position of the lower boundary δ 1 was determined using Eqs. (53, 54). To assess the angular length and duration of one period of the trajectory and the horizontal position of other important points, the same procedure as in Sect. 5.2 can be applied.  Fig. 23, we can evaluate the character of one period. The graphs in Fig. 23a, b show the strip width expressed in α or δ variables. Picture (c) presents the spatial width of one period as a function of ω n . It is obvious that increasing the ISV leads to increasing the amplitude δ 1 (or α 1 ) up to a maximum reached for ω n = ω ns . At the same time, the length of the period reaches its maximum for ω n = ω ns . The periodicity of the trajectory in the interval 0 < ω n < ω ns is obvious. The half-period is denoted T γ (angular scale). Doubtlessly, it can be expected that the spiral is irrational like in Sect. 4.1, except for some special cases.

The limit case-the trajectory passing through the SPC
The second trajectory type represents the trajectory that separates the "lower and upper" groups with respect to frequency ω ns , see the bold-green curve in Figs This means that for the fixed elevation δ c of the SPT, and the implicitly defined velocity Γ 0 , the following ISV should be applied: in order to produce a trajectory the loops of which go through the SPC; see also Eq. (56). Coordinate γ becomes discontinuous when approaching the SPC; see the discontinuity in the green curve and the jump of length π in Fig. 22. For the same reason, the derivative of the curve describing the width of the strip in ω ns is non-continuous (Fig. 23a). However, the one-sided derivatives with respect to ω n at point ω ns are finite and equal in absolute value. Therefore, the non-smooth character of the curve is merely a result of maintaining the polar angle α as positive. Similar reasoning regards the discontinuity in the length of the spatial period in Fig. 23c. The 2π jump originates also from the fact that the relevant radius-vector length is always positive.

Trajectories passing by the SPC
The third type of trajectory can be observed when ω n > ω ns and possibly ω n → ∞ as the blue curves in Figs 20c, 21, 22 and also in case (iii) of Fig. 24. These initial conditions result in trajectories that pass by the SPC, which remains outside each loop; see Fig. 20c. The position of the lower boundary δ 1 can be obtained from Eq. (53). It rises monotonously toward δ c for an increasing ω n . Accordingly, as ω n → ∞, we can observe that the width of the spherical strip diminishes; see Fig. 23a, b.
The angular length of one half of the period and the position of the turnabout points γ v1 , γ v2 , see Fig. 22, can be evaluated in a similar way as in the case of the third-type trajectories for the ω n < 0, i.e., to employ system Eq. the coordinate γ as a function of spin frequency ω n > 0; note a jump in the point ω n = ω ns changes direction. Their attributes are similar to those encountered in the previous section, although their geometrical interpretation is slightly different. Also here, we can define a bottom part of the loop within points γ v1 , γ v2 where velocityγ is of the opposite sign. In contrast to Sect. 5.2 (negative ISV), the loops do not include any twofold point, and the period length γ b2 − γ b1 becomes negative for ω ns < ω n ; see the two blue curves with an indication of the movement direction in Fig. 22. Both properties of the loops follow from the ordering of starting and finishing points of the period and are typical for trajectories with a positive ISV. It is obvious that as ω n → ∞, the particular loops approach a zero amplitude; see also Sect. 5.4. The same also applies to the length of period T γ . With reference to points γ v1 , γ v2 and the starting and finishing points of the period γ b1 , γ b2 , we can define the average (or effective) angular velocity Ω av < 0 of the ball movement forward along the SC for ω ns ω n . It is obvious that the clockwise movement of points γ b1 , γ b2 , γ v1 , γ v2 between consecutive loops gets slower for increasing ω n , i.e., as Ω av → 0. The limit case of ω n → ∞ is discussed in detail in Sect. 5.4; see also Fig. 26.

High initial spin values
In this section, we will investigate trajectories starting from the SC, when an extremely high ISV is applied in either the positive or negative sense. We have seen in Sects. 5.2 and 5.3 that the shape of individual trajectories differs significantly if a positive or negative velocity of initial spin is introduced, even though some analogy can be noticed.
Significantly increasing the ISV, the trajectory properties for ω n 0 and ω n 0 become more and more related and finally produce a symmetric image with respect to the SC. Both of them are far from the transition case, in which the curly form trajectories start. For ω n 0, the curly trajectory is located on the upper side of the SC. The upper boundary δ 2 descends from a level above the SC to root δ c = δ 1 as ω n → −∞. Provided ω n 0, the lower boundary (root δ 1 ) moves upwards to the SC, represented by δ c = δ 2 , when ω n → ∞.
Let us outline approximately this process. Either roots δ 2 , δ 3 or δ 1 , δ 3 (depending on which trajectory group is being analyzed) can be determined using the quadratic equation Eq. (52) or its symbolic solution Eq. (53), respectively. For a high |ω n |, the term L becomes dominant due to the higher power of |ω n |, so it holds: L 2 K · M, and, consequently, it can be approximately written for i = 1, 2: It is obvious that δ 3 → ∞ and is therefore of no interest. The other root can be approximated: 24 Examples of trajectories with a positive "Initial Spin Velocity ISV" (ω n > 0): 1st row: 0 < ω n < ω ns ; 2nd row: ω n = ω ns ; 3rd row: ω n > ω ns . Column a time history, b top view, c axonometric demonstration The second and higher terms can be neglected, as they vanish for |ω n | → ∞. Substituting now for K , M from Eq. (54b), and neglecting terms with the zero degree of ω n , we obtain Formula Eq. (68) shows that δ i converges to δ c from above or below depending on the sign of the ISV. The width of the spherical strip decreases for increasing |ω n | and so does the width of the curly trajectory period, as it has been defined in Sects. 5.2 and 5.3; see the plots in Fig. 25. Graph (i) represents the case ω n 0 advancing from left to right, whereas (ii) regards ω n 0 from right to left. The shape of the individual loops for |ω n | 0 approaches a circle, which is followed in a negative or positive angular sense (with respect to the The increment of the distance which the circle performs during one loop decreases rapidly with growing |ω n | because the distance between the starting and finishing points γ b1 , γ b2 of one period (positive or negative) decreases faster than the circle diameter. In general, the effective distance passed along the SC within one period of a curly trajectory is mostly exploited by the approximating circle. Therefore, the effective velocity Ω av of the ball advancement along the SC decreases. Let us remind the reader, referring to Sects 5.2 and 5.3, that the sense of this slow rotation around the z axis is either positive (Ω av > 0, ω n 0)-above the SC, or negative (Ω av < 0, ω n 0), but at any time it holds that Ω a v → 0, when |ω n | → ∞.
This deduction concludes with a phenomenon which seems paradoxical at first glance. As we have seen, as ω n → ±∞, both the diameter and effective horizontal velocity of the approximating circle degenerate to zero. Therefore, the ball does not evince any movement along the SC, the loops reduce to a single point, and the ball is apparently at rest, only spinning around its normal at the SPT with an infinite spin velocity. Its position, which seems to be out of static equilibrium, is fixed by an infinitely strong gyroscopic effect, whatever the IHV is. This type of trajectory is illustrated as a numerical simulation in Fig. 26: (a)  6 Vertical plane-related trajectories

Transformation of the governing system
We examine trajectories which emerge for a very small initial horizontal velocity (0 ≤ |γ 0 | Γ ) and a zero or small initial spin. If both IHV and ISV are zero, the trajectory of the ball defines a vertical plane passing through the initial point, ball center and point A of the cavity (SPC). With an injection of a small IHV or ISP, the trajectory mildly declines from this vertical plane. However, the distance from the vertical plane remains small. To describe the character of this spatial curve, it is worthwhile to consider the difference from the vertical plane as a small parameter. This simplification enables us to get into some special properties of this trajectory family.
With reference to Fig. 3, root δ 1 will either coincide with the origin, i.e., δ 1 = 0, or get a small positive value. The position of root δ 2 = δ c depends on the initial position of the ball.
To conveniently describe the movement of the ball in a narrow strip adjacent to plane yz, where the azimuthal angle γ exhibits a jump of 2π when α passes the zero value, it is advisable to rotate the coordinate system around the y axis.

Approximated governing system and relevant trajectories
The small parameter can be introduced as follows: ξ ≈ π/2 − ε; see Fig. 27. Hence, Eq. (71) gets a modified form: Provided no IHV or ISV is applied, ε vanishes. Equation (72a) is fulfilled identically, and Eq. (72b) degenerates into a nonlinear pendulum equation. This equation can be solved in elliptic functions. Depending on the initial velocityζ (0), the solution is either periodic-for smallζ (0)-or continuously increasing in time with periodically variable velocity. Detailed discussion can be found, for instance, in [29,34].
On the other hand, cases for ζ 0 ≥ π/2 are physically meaningless due to a negative contact force. Note that the case considering ζ = 0 and ε = 0 represents only a perturbation, which can be treated using linearized expressions.
Under the assumption that the amplitudes of ζ are small, the nonlinear terms in Eq. (72) may be approximated as sin ζ ≈ y ρ , cos ζ ≈ 1 andζ 2 ε ≈ 0,εεζ ≈ 0, and the system can be linearized. The following reduced system can be formulated: For a nonzero ISV, equations Eq. (73) are coupled by a pair of gyroscopic forces which cause the spatial character of the trajectory even if the IHV vanishes. For a zero ISV, the system is characterized by two independent linear oscillators with identical eigenfrequencies. Consequently, no rosette form trajectory can occur for a non-trivial IHV, and a simple ellipse-like curve is observed. Indeed, considering initial conditions: IHV: z(0) = z * 0 , and an initial deviation along the meridian in the plane x y: y(0) = y 0 , a solution to linear system Eq. (73) can be expressed as follows: which produces a pair of independent components: w = y 0 cos ω 0 t + iz * 0 /ω 0 sin ωt, if the ISV ω n = 0. Provided ω n = 0, but is rather small, the top view of the trajectory obviously has the form of a strongly prolate hypotrochoid. Therefore, it is worthwhile to define an affine space which meets a sub-manifold at a point in such a way as to have a second order of contact at this point. From a geometrical point of view, this case is similar to the one we discussed in Sect. 4.2. Indeed, we can assume that the shape of the spiral does not deviate much from the vertical plane during one cycle. The osculating plane rotates slowly around the vertical axis with an angular velocity Ω v = −μω n /2. Several plots of trajectories following this linear approach are demonstrated in Fig. 28, row (i), without or with an ISV, pictures (a) or (b), respectively.
When the ζ amplitude is large and the nonlinear character of the system must be respected, the second terms in Eq. (72) remain in force. In a certain sense, they also have the character of gyroscopic forces and lead to a rosette character of the trajectory when a small IHV is applied, see Fig. 28, row (ii), including details near the upper boundary of the strip. Observing picture (b) (a nonlinear approach), we can see that Ω v > 0 when only IHV is considered and no ISV is applied (counterclockwise). If ω n > 0, then Ω v < 0 (clockwise).
As a demonstration, we show a trajectory resulting from a numerical simulation; see Fig. 29. The scheme of this figure is the same as that in Fig. 6. Although an ISV was introduced in this initial setting, we can see in the top view that a slightly counterclockwise rotating spiral emerges. Subsequent simulations for the descending SPT level showed a decreasing value of this rotation velocity and vice versa, just as we have seen above for small initial amplitudes (the level of the SPT) [Eq. (74)], and higher amplitudes [Eq. (72)], where the nonlinear character was respected. Hence, the top view of the trajectory is similar to that of the Foucault pendulum even if no ISV is applied. Note that the individual loops go around point A (SPC) regardless of the initial setting if no ISV is applied, i.e., for ω n = 0.

Conclusions
The dynamic behavior of a non-holonomic system represented by a ball moving inside a spherical cavity has been investigated using an analytical approach based on the Lagrangian governing system. The rolling of the ball is considered to be slipping-less and free of damping at the contact of the ball and cavity. The cavity is assumed to be fixed. Energy is introduced into the system by means of appropriate settings of initial conditions. The reason for this system layout consists in the possibility to differentiate individual groups of trajectory types and to investigate each one separately. This method of investigation indicated a set of trajectory types and limit cases that either separate them or demonstrate a process of limitation of some parameters to certain special values (infinity, zero, etc.). The results were compared to those obtained numerically by the principally different method following from the Appell-Gibbs approach. It can be concluded that the results coincide perfectly. The two sets of results com-(a) (b) Fig. 28 Trajectory at a low IHV: line (i) linear approach-low level SPT, line (ii) nonlinear approach; column a no initial spin, column b initial spin included plement one another well with regard to their strengths and weaknesses. The positions of limit cases match with respect to their relevant parameter settings. Using the analytical approach, the transition zones were also qualitatively and quantitatively examined and interpreted, which was impossible using solely numerical simulations.
For the investigation of particular properties of the system, two types of basic relations were used: (1) a Lagrangian governing system with incorporated Pfafftype non-holonomic constraints and (2) its three first integrals corresponding to the total energy of the ball, the constant angular momentum with respect to fixed vertical axis and the constant angular momentum with respect to the common normal at the contact point of the ball and cavity. The last one showed that spin velocity is constant throughout the whole period the ball is moving along its trajectory.
A cubic algebraic equation makes up the backbone, characterizing the energy flow within the system. This equation always possesses three real roots, where the two lower roots are physically meaningful. They delimit a spherical strip on the spherical surface of the cavity inside of which the relevant trajectory emerges. The properties of these roots enabled us to define individual trajectory types, indicate transition, limit and other cases as well as perform parametric analyses within each group of trajectories.
-The "Separation Circle" (SC), a horizontal circular trajectory at a certain level of the lower hemisphere of the cavity, was revealed to be the main classification element. The selected level on a meridian of the cavity determines the particular critical value of "Initial Horizontal" and/or "Spin Velocity" (IHV or ISV, respectively) necessary to maintain this trajectory. Trajectories within the close neighborhood of the SC were examined, confirming the stability of the SC trajectory with respect to perturbation in initial conditions and to a small cross-impulse imparted to the ball. -Trajectory types above the SC, assuming a positive IHV, are characterized by an IHV that is higher than the critical velocity and/or a nonnegative ISV. As a rule, spin-free trajectories have a form of irrational spirals similar to a prolate hypotrochoid occurring within the spherical strip. Increasing IHV beyond all limits, one approaches the limit state, which is represented by a planar trajectory that is symmetrically distributed with respect to the cavity center. -Settings with a positive ISV can be classified into two sub-groups separated by the case that is represented by a "kings-crown" shaped trajectory. In this particular case, every loop contains a singular point (apex) which lies on the upper boundary of the spherical strip. A closed loop expands from this point for higher spin velocities, imparting a curly form to the trajectory. Lower spin velocities produce spirals without loops analogously to cases without spin, despite the fact that the shape of the spiral itself is different. -Trajectories below the SC have a significantly different character. Depending on the decreasing IHV, the lower limit of spin-free trajectories successively descends down to the "Southern Pole of the Cavity" (SPC), where it degenerates into a point. Regarding cases with negative spin, they can also be classified into two sub-groups. The separating case is represented by a curly trajectory, where every loop passes through the SPC. Higher spin velocities result in trajectories the individual loops of which go round the SPC, while lower spin velocities produce loops that miss this point. Thus, the lower boundary of the strip again rises toward the SC. -An interesting process is encountered when the ISV velocity tends to ±∞. The width of the spherical strip where the trajectory is being traced decreases to zero; the trajectory has a shape approaching a small circle slowly moving along the SC. For infinite spin velocity, the ball is seemingly fixed at the initial point of the trajectory due to infinite energy concentrated in the spin of the ball producing an infinite gyroscopic effect. -The last group includes associated cases characterized by low initial energy of the ball, i.e., the IHV and/or ISP of the ball is small or vanishing when the trajectory starts at a certain height. For zero ISV and low-level initial position, the prob-lem can be linearized. In the top view, a simple ellipse emerges, degenerating into a line segment for zero IHV. If nonlinear terms are respected, then this ellipse-like curve slowly rotates counterclockwise and remotely resembles the Foucault pendulum trajectory. Introducing a slight ISV, the elliptical trajectory shape changes into a significantly prolate hypotrochoid.
Ball-type vibration absorbers are equally effective regardless of the particular direction of excitation, yet they are mostly designed for unidirectional cases, sometimes with the assumption of small angles. It is rare to take into account the auto-parametric effects that stem from the nonlinear nature of the system, and which cause the spatial movement of the sphere in the cavity even under a unidirectional load. Although the movement of the sphere in the presence of external excitation differs from the trajectories described in this work, their character is broadly similar. As shown by numerical simulations of an externally excited case and also by experiments in analogous pendulum absorbers, the spatial motion of the damping mass along paths related to those described above shows an undesirable stability against fluctuations of possible ambient excitation. In the case of pendulum absorbers, such phenomena can be effectively mitigated by applying sufficient damping. However, this can be a problem in the case of balltype absorbers. The spatial resonance movement of the ball in a limit cycle may negatively affect the structure. Thus, it seems necessary to include a detailed autoparametric analysis of the complete system during the design stage.
As regards future investigation in this area, nonconservative systems should be paid attention on the basis of analytical processes. Specific difficulties should be expected related to first integrals, which will need to be generalized significantly in order to overcome variable total energy, angular momentum and spin velocity. The mutual penetration of groups and sub-groups when moving along one particular trajectory will have to be taken into account. Chaotic response processes will doubtlessly emerge, as has already been shown when doing numerical simulations. It will be necessary to reconcile with the fact that some steps used in this paper will become inapplicable. Nevertheless, some alternative methods of analysis are emerging and seem promising.