On Hopf bifurcation in the problem of motion of a heavy particle on a rotating sphere: the viscous friction case

We investigate the Hopf bifurcation of a mass on a rotating sphere under the influence of gravity and viscous friction. After determining the equilibria, we study their stability and calculate the first Lyapunov coefficient to determine the post-critical behavior. It is found that the bifurcating periodic branches are initially stable. For several inclination angles of the sphere’s rotation axis, the periodic solutions are calculated numerically, which shows that for large inclination angles turning points occur, at which the periodic solutions become unstable. We also investigate the limiting case of small friction coefficients, when the mass moves close to the equator of the rotating sphere.


Introduction
The motion of bodies on surfaces is a classical problem of mechanics [1,2] and was investigated in various statements (see, for example, [3][4][5][6]). In the most simple case, a point particle instead of a rigid body could be considered. That kind of problems appears when we study the dynamics of mechanical systems with rotating parts performing different operations such as the mixing, grinding, and drying, of diverse substances. Based on results of computer simulations [7,8], it is possible to investigate the dynamics of systems with a large number of particles. However, the output of such a simulation usually does not represent any analytical results. That is why it is reasonable to consider simple systems, such as a particle moving on a surface under the action of a friction force. Even in these simple cases, some complex dynamic effects can be discovered. If there is a friction force and the surface doesn't move, the system will come to rest. However, if we assume that the surface rotates with a constant angular velocity, steady and periodic motions can appear in the system. This fact also makes it possible to use such problems to identify the friction coefficient. The problem of the motion of a point particle on a rotating surface was studied in [9]. The motion of a particle on a rotating table was investigated analytically and numerically in [10].
The problem of motion of a heavy bead on a circular hoop rotating about its vertical diameter has been studied in [11]. The similar problem for a circular hoop rotating about some other vertical axis has also been investigated [12]. In the present paper, a three-dimensional analogue of this problem is studied under the assumption that there is viscous friction between the point and the sphere.

Definition of the problem and equations of motion
Let P be a heavy particle of mass m which moves on a two-dimensional sphere of radius under the action of a viscous friction force, with the coefficient of friction being c. The sphere rotates with a constant angular velocity ω about a fixed axis. It is assumed that the axis passes through the center of the sphere O.
Let Ox 1 x 2 x 3 be an absolute coordinate system such that the plane Ox 1 x 2 is horizontal and the axis e 3 = Ox 3 is directed along the upward vertical, so the gravity force acting on the particle is −mge 3 . Assume that the rotation axis belongs to the plane Ox 1 x 3 and its angle of inclination is α, 0 ≤ α ≤ π/2. Denote the spherical angles, which specify the position of the moving particle on the sphere, by θ and ϕ, 0 ≤ θ ≤ π, 0 ≤ ϕ < 2π (Fig. 1).
In this case, the motion of the particle P can be described by the following system [15]: where System (1) possesses a symmetry w.r.t. mirror reflection about the (x 1 , x 3 )-plane: Inserting , we obtain the same system of differential equations in the new variables.

Equilibria positions and their stability
Assumingθ = 0,φ = 0,γ = 0,˙ = 0 we obtain the equations that determine the equilibria of the system: The system (2) can be solved for (θ 0 , ϕ 0 ) by introducing the quantities Then (2b) gives which can be inserted in (2a): By squaring both sides in (3) and rearranging, we obtain the quadratic equation for v 2 Since increases monotonically for v > 0, Eq. (4) has one positive solution for v 2 , the second solution of the quadratic equation leads to complex-valued solutions of ϕ 0 and θ 0 .
Due to the mirror reflection symmetry mentioned above, only two of these four solution pairs need to be investigated further, the stability and bifurcation behavior of the other two steady states follows by applying the mirror reflection. The parameter χ is a bifurcation parameter for system (1). However, it is convenient to introduce a parameter λ = 1 χ instead of χ . Let us assume now that whereθ,φ,γ , andˆ are small perturbations of the equilibria. Using (6) to eliminate ϕ 0 , we can write the series expansion of the perturbed system up to the cubic terms as follows: = cos α λ sin θ 0 cos θ 0θ The linear part of the system (7) isẋ where and The characteristic equation of the linearized system is where The Routh-Hurwitz criterion leads to the following stability conditions: Now we will consider the obtained conditions. Taking into account the form of the coefficient b 1 , the first condition is equivalent to The third condition b 3 > 0 is equivalent to Taking into account (11) and (12), the second and the fourth conditions are also satisfied. The last condition R > 0 can be rewritten as follows: This condition holds when Since b 4 = det(F(λ)) > 0 for p = 0 and cos α = 0, there cannot occur a steady-state bifurcation. The boundary of a stability region is given by the condition R = 0. On this boundary the characteristic equation has a pair of purely imaginary roots. The type of the second pair of roots on the boundary depends on the sign of the expression * Thus, on the boundary R = 0 the characteristic equation (9) has the roots As shown in [15], in this case Thus, the system has two purely imaginary roots and two roots with a negative real part. This fact allows us to suggest a branching of a limit cycle when λ = λ 0 . Such a bifurcation was found numerically, and so it is also reasonable to study the bifurcation analytically.
The dependence of the system behavior near the stability boundary R = 0 on the parameter λ is determined by theorems given below. Consider the variation of the parameter λ on a small interval λ 0 − η ≤ λ ≤ λ 0 + η, where λ 0 can be found from R(λ 0 ) = 0. Let a = L 1 (λ 0 ) be the first Lyapunov coefficient given below that determines the stability for λ = λ 0 . Neglecting higher-order terms, the behavior of the system close to the bifurcation point can be described by the one-dimensional differential equation for the radius of the bifurcating periodic solutionṙ where d = − dσ c /dλ ∝ d R/dλ, evaluated at the Hopf bifurcation point.
Thus, with increasing λ, the stable equilibrium becomes unstable, but the image point stays in an εneighborhood of the equilibrium. With decreasing λ, the equilibrium becomes stable again and the image point returns to the equilibrium, as shown in Fig. 2. The behavior of the system is reversible with respect to λ.
Thus, with increasing λ the stable equilibrium becomes unstable, and the point leaves the neighborhood of it. When λ is decreased again below λ 0 , the point will usually not return to the equilibrium. The behavior of the system is therefore irreversible with respect to λ. A typical situation is shown in Fig. 3 for a locally subcritical, but globally supercritical Hopf bifurcation: If λ increases beyond the critical value λ c , the system quickly moves to the stable large amplitude oscillation and remains there, even if λ decreases below λ 0 again. Only if it reaches the limit point cycle ("LP"), it jumps back to the stationary state, exhibiting a hysteretic behavior.
This theorem has also been proved in [13].
Since cos α > 0 and p > 0, this expression is negative. Therefore, the system satisfies the conditions of Theorem 1, and the boundary R = 0 is safe. Periodic modes that appear when λ>λ 0 are stable.

Numerical investigation of the bifurcating solutions
Since by (14) the first Lyapunov equation L 1 (λ 0 ) is negative for all parameter values, the bifurcating periodic solutions exist locally for λ ≥ λ 0 , respectively, for friction coefficients χ ≤ χ 0 = 1/λ 0 . In order to get the global behavior, several branches of periodic solutions were computed by the BVP solver Boundsco [18] and the continuation algorithm Hom [19] for fixed values of p = 3 and inclination angle α for varying values of the distinguished parameter χ . As can be seen in Fig. 4, the periodic solutions are throughout stable for small inclination angles α and extend down to χ = 0. For inclination angles α > α c ≈ 54.8 • , the branches display turning points, where the stability of the periodic solutions changes. For values of χ between the turning points, three different periodic solutions are found, two of these are stable and the medium one is unstable. For α = 75 • and χ = 0.02, the periodic solutions are displayed in Fig. 5. The largest orbit is almost aligned with the equator of the rotating sphere, whereas the smallest solution oscillates close to the bottom of the sphere.
During the path-following along the periodic solutions, we encounter a singularity of the differential equations due to the use of spherical coordinates: Close to the Hopf bifurcation point the solutions are periodic in θ and ϕ; but after the trajectories pass through the south pole of the sphere, the azimuthal angle ϕ increases by 2π during one period. In order to overcome this difficulty, a different coordinate system (Cartesian or spherical coordinates along the rotation axis) has been used close to the crossing of the south pole, when the equations become singular.

Limiting behavior for small friction coefficients
When the parameter χ is set to zero, all gravitational and damping forces vanish and the mass can move freely on the sphere, tracing out arbitrary great circles. For small values of χ , or equivalently, for large rotation speeds ω, the numerical calculations indicate that the periodic solution approaches the equator of the rotating sphere and rotates with the same speed as the sphere. That behavior is also expected by mechanical reasoning: If we neglect the gravitational force, the friction and centrifugal force will cause the mass to move along the equator; a small gravitational force causes a periodic excitation. In order to study this motion, we use spherical angles aligned with the rotation axis of the sphere. In these coordinates, the friction terms become simpler, by setting α = 0 in (1), because in the new coordinate system the sphere rotates about the z-axis. Since now the gravity acts in the direction e g = (sin α, 0, − cos α) T , its potential becomes V = −χ pe g · (sin θ cos ϕ, sin θ sin ϕ, cos θ) T = χ p(cos θ cos α − sin θ cos ϕ sin α).
The equations of motion are therefore given bẏ For p = 0 a solution of (15a)-(15d) is given by the steady rotation θ = π/2, ϕ = t along the equator. Since Eq. (15d) for is dominated by the term χ(1 − ), the azimuthal dynamics will change slightly under the perturbation by p: At leading order, we obtain from (15d) where the integration constant has been chosen such that the average perturbation of vanishes. Setting θ = π/2 + ϑ, 2 ≈ 1 + 2 pχ cos t sin α and expanding Eq. (15c) up to first order in χ , we obtain the linear equation with parametric excitation ϑ + χθ + (1 + 3 pχ sin α cos t)ϑ = pχ cos α. (17) One might think that the leading order expansion yields a good approximation for small χ , but the numerical solutions show a large periodic deviation from this estimate, which occurs, because the parametric excitation frequency in (17) is in 1 : 1 resonance with the eigenfrequency of the unperturbed equation. In order to find the periodic component of ϑ, we set ϑ = ϑ 0 + ϑ 1 and obtain the differential equation In order to find the influence of the parametric excitation term 3 pχ sin α sin t ϑ 1 , we apply two steps of Normal Form reduction to the equation to find the periodic solution which satisfies (19) up to terms of order O(χ 2 ). As demonstrated in Fig. 6, this approximation agrees well with the numerically obtained periodic solution of the nonlinear equation, while the steady approximation is very poor.

Conclusions
The Hopf bifurcation of a moving mass on a rotating sphere has been investigated. By transforming the system to Jordan Normal Form, calculating the Center Manifold and simplifying the system using Normal Form theory we obtained a simple expression for the first Lyapunov coefficient. Since this coefficient is negative, the periodic solutions bifurcate supercritically from the steady state. These bifurcating solutions are also computed numerically, and their limiting behavior for vanishing friction force is studied.
Acknowledgements Open access funding provided by TU Wien (TUW).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A Calculation of the Lyapunov coefficient
For the case under consideration, a pair of eigenvalues ±iω on the imaginary axis and a second complex conjugate pair m ± in with m < 0, the Lyapunov coefficient was already derived in [13]. Since this formula looks extremely complicated, the contributions from the cubic terms in (24) from Normal Form calculations and from Center Manifold Reduction are stated separately.

A.1 Transformation to Jordan Normal Form
The first step to calculate the bifurcating solutions is to introduce a linear change of coordinates, which transforms the linearized equations (8) at the steady state to real Jordan Normal Form: Let v c and v s be the complex valued eigenvectors of the Jacobian F corresponding to the eigenvalues iω and m + in, respectively, then the matrix satisfies transforms system (7) to the new nonlinear systeṁ where F NL (x) contains the nonlinearities of (7). The explicit formulas for the entries of V and f (ξ ) are given below.
In the transformed coordinates, the system takes the following form: where b, n, m are defined by the solutions of the characteristic equation (9).