# Bifurcation analysis of the problem of a “rubber” ellipsoid of revolution
rolling on a plane

A. A.Kilin, E. N.Pivovarova
Ural Mathematical Center,
Udmurt State University,
ul. Universitetskaya 1, Izhevsk, 426034 Russia

###### Abstract

This paper is concerned with the problem of an ellipsoid of revolution rolling
on a horizontal plane under the assumption that there is no slipping at the
point of contact and no spinning about the vertical. A reduction of the
equations of motion to a fixed level set of first integrals is performed.
Permanent rotations corresponding to the rolling of an ellipsoid in a circle
or in a straight line are found. A linear stability analysis of permanent
rotations is carried out. A complete classification of possible trajectories
of the reduced system is performed using a bifurcation analysis. A
classification of the trajectories of the center of mass of the ellipsoid
depending on parameter values and initial conditions is performed.

###### Contents

  1. 1 Introduction
  2. 2 Equations of motion for a body of revolution on a plane
    1. 2.1 Configuration space and kinematic relations
    2. 2.2 Constraint equations and dynamical equations
    3. 2.3 Invariants and symmetries
    4. 2.4 Reduction to a system with one degree of freedom
  3. 3 Explicit form of the equations of motion for an ellipsoid of revolution
  4. 4 Bifurcation analysis
    1. 4.1 General definitions
    2. 4.2 Critical sets
    3. 4.3 Example of a bifurcation diagram
    4. 4.4 Bifurcations of critical sets and types of bifurcation diagrams
  5. 5 Analysis of the trajectories of an ellipsoid

##  1 Introduction

Analysis of the motion of rigid bodies on a plane is one of the classical
problems in mechanics. To describe the motion of a rigid body in the first
approximation, use is often made of conservative models: the model of an
absolutely smooth plane (the point of contact can slip) and the model of an
absolutely rough plane (the classical nonholonomic model in which there is no
slipping at the point of contact). Another conservative model of contact is
that of a rolling “rubber” body — it assumes that there is no slipping at the
point of contact and no spinning about the vertical. Early studies of the
rolling motion of bodies without spinning go back to the classical works on
modern mechanics: this model was proposed by Hadamard [1] and developed by
Beghin [2]. In [3], the term “rubber rolling” was introduced to describe this
motion model. A systematization of results on the rolling motion of rigid
bodies using the above-mentioned model was presented in [4, 5], where existing
tensor invariants of the system depending on the shape of the body and the
mass distribution in it were found.

The main feature of the problem of a body of revolution moving on a plane is
its integrability regardless of the shape of the body. The same holds true for
an absolutely smooth plane [6, 7], for an absolutely rough plane [8], and for
the rubber body rolling model [4]. Moreover, if the system is subject to one
nonholonomic constraint which prohibits sliding in the direction of the
projection of the axis of rotation, then such a problem is also integrable.
This was shown in [9], where the authors addressed the problem of the rolling
motion of a sphere with axisymmetric mass distribution. It was shown that this
system admits a redundant set of first integrals and an invariant measure, and
that one of the found integrals is similar to the additional integral found in
[10] in the problem of a circular disk rolling on smooth ice. This allowed a
reduction to a system with one degree of freedom and made it possible to show
that all nonsingular trajectories are periodic functions of time. It is
interesting that, if the system is subject to a constraint prohibiting sliding
in the direction perpendicular to the axis of rotation, then the problem
becomes nonintegrable [11]. For the case of the sphere with axisymmetric mass
distribution it is shown that, in the general case, this system is
nonintegrable (one additional integral of motion and an invariant measure are
needed for the system to be integrable). It is numerically shown that the
system exhibits chaotic dynamics. In addition, the problem of a body of
revolution rolling on a plane remains integrable after adding an additional
potential that is invariant under rotation about the vertical [12, 13, 14].

There is a related problem, that of a homogeneous sphere rolling without
slipping on a surface of revolution, which similar to the problem of a body of
revolution rolling on a plane. This problem is also integrable regardless of
the form of the surface. This is shown in [15, 16] using the same methods of
analysis (in particular, reduction by symmetries). In [17, 18] the rolling
motion of a sphere on an arbitrary rotating surface of revolution is
investigated, and the authors of [19, 20] examine integrable cases of this
problem: the rolling motions of a sphere on the rotating surface of a cone and
a circular cylinder. In [21, 22] a new interesting problem concerning the
rolling of spheres between two spherical surfaces is addressed and cases of
its integrability are found, and in [23] the integrability of the related
problem of the motion of a rigid body in a spherical suspension is shown [24].

Compared to the study of body motion using the classical nonholonomic model
prohibiting only slipping at the point of contact, much less attention has
been given in the literature to the dynamics of bodies using the rubber body
model. In particular, analysis of rigid body dynamics using this model was
made in [25, 26, 27, 3]. In [28], the rubber body model was employed to
examine the dynamics of an unbalanced disk. The rolling motion of spheres with
different mass distributions was studied in [25, 26, 27, 3, 29, 30, 31, 32,
33], and the motion of a triaxial ellipsoid was examined in [34].

However, the “rubber” motion model is applied not only in handling some model
problems concerning, e.g., disks or spheres, but also in simulating the motion
of mobile robots. In particular, the problem of using various motion models
for simulating the dynamics of mobile robots is discussed in [35], where it is
noted that the rolling motion of a sphere without slipping or spinning is the
classical ball-plate problem [36, 37]. In [35], the authors also note that in
many studies the vertical component of the angular velocity of the spherical
robot is assumed to be zero since the possibility of omnidirectional motion
allows motion control without using vertical rotations. Such a motion model is
used in analyzing the dynamics of a spherical robot with a pendulum drive [38,
39], a spherical robot with an internal platform [40, 41], a spherical robot
with internal omniwheels [42, 43], and a spherical robot actuated by changing
the gyrostatic momentum [44]. In addition, some authors [41, 42, 45, 46], who
use the “rubber” motion model consider only kinematics without taking into
account the system dynamics.

In this paper, we carry out a qualitative analysis of the motion of an
ellipsoid of revolution on a horizontal plane using the model of rubber body.
Although we restrict our attention to a special case of body shape, many
conclusions of this paper will be valid for a body of an arbitrary shape as
well. In addition, the analysis of the dynamics in specific problems remains
of interest to researchers (see, e.g., the recent paper on qualitative
analysis of the rolling motion of a torus on a plane [47]).

The paper is organized as follows. Section 2.1 presents equations of motion,
invariants, and a reduction to a system with one degree of freedom for a body
of revolution of arbitrary form. In Section 3, an explicit form of the
equations of motion for an ellipsoid of revolution is presented. Section 4
provides a qualitative analysis of the dynamics of the ellipsoid and a
complete classification of motions depending on parameters and values of first
integrals. In Section 5, an analysis of the trajectories of the center of mass
of the ellipsoid in absolute space is made and conclusions concerning the
boundedness and unboundedness of the trajectories are drawn.

##  2 Equations of motion for a body of revolution on a plane

###  2.1 Configuration space and kinematic relations

Consider the rolling motion of a dynamically symmetric rigid body of
revolution of mass $m$ on an absolutely rough horizontal plane (Fig. 1). We
will assume that the body rolls without slipping (the velocity of the point of
contact is zero) and without spinning (the projection of the angular velocity
of the body onto the vertical is zero), and has only one point of contact,
$P$, with the supporting plane.

Figure 1: A body of revolution on a plane.

We introduce two coordinate systems:

  1. –

a fixed (inertial) coordinate system $Oxyz$ with unit vectors
$\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}$, where
$\boldsymbol{\gamma}$ is a vertical unit vector;

  2. –

a moving coordinate system $Cx_{1}x_{2}x_{3}$ with vectors
$\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}$ rigidly attached to
the body, and with origin at its center of mass, $C$, such that the axis
$Cx_{3}$ is directed along the symmetry axis.

We will specify the position of the body of revolution by the coordinates of
its center of mass $\boldsymbol{r}_{c}=(x_{c},y_{c},z_{c})\in\mathbb{R}^{3}$
in the coordinate system $Oxyz$, and its orientation in space, by the
orthogonal matrix ${\bf Q}$, whose columns are the projections of the unit
vectors $\boldsymbol{\alpha},\,\boldsymbol{\beta},\,\boldsymbol{\gamma}$ onto
the axes of the moving coordinate system $Cx_{1}x_{2}x_{3}$

${\bf Q}=\left(\begin{aligned} &\alpha_{1}&\beta_{1}&&\gamma_{1}\\\
&\alpha_{2}&\beta_{2}&&\gamma_{2}\\\
&\alpha_{3}&\beta_{3}&&\gamma_{3}\end{aligned}\right)\in SO(3).$

Hence, the configuration space of the problem of the free motion of the body
of revolution is the product $\mathcal{N}=\\{(\boldsymbol{r}_{c},{\bf
Q})\\}=\mathbb{R}^{3}\times SO(3)$.

To describe the dynamics of the body on a plane, we introduce the vector of
the angular velocity $\boldsymbol{\omega}=(\omega_{1},\omega_{2},\omega_{3})$
and the vector of the velocity of the center of mass
$\boldsymbol{v}=(v_{1},v_{2},v_{3})$ and refer them to the moving coordinate
system $Cx_{1}x_{2}x_{3}$111Here and below, unless otherwise specified, all
vectors will be referred to the moving coordinate system $Cx_{1}x_{2}x_{3}$..
Using these velocities (quasi-velocities), the evolution of the vector
$\boldsymbol{r}_{c}$ and the matrix ${\bf Q}$ can be written as the kinematic
relations

$\begin{gathered}\dot{\boldsymbol{\alpha}}=\boldsymbol{\alpha}\times\boldsymbol{\omega},\quad\dot{\boldsymbol{\beta}}=\boldsymbol{\beta}\times\boldsymbol{\omega},\quad\dot{\boldsymbol{\gamma}}=\boldsymbol{\gamma}\times\boldsymbol{\omega},\\\
\dot{\boldsymbol{r}}_{c}={\bf Q}^{\top}\boldsymbol{v}.\end{gathered}$ (1)

###  2.2 Constraint equations and dynamical equations

We assume that the body rolls without loss of contact with the plane. This
assumption implies that the following holonomic constraint is imposed on the
system:

$z_{c}+(\boldsymbol{r},\boldsymbol{\gamma})=0,$ (2)

where $\boldsymbol{r}$ is the radius vector of the point of contact (see Fig.
1), which can be expressed in terms of the components of the vector
$\boldsymbol{\gamma}$. Here and in what follows, by
$(\boldsymbol{\cdot},\boldsymbol{\cdot})$ we denote the scalar product. In the
case of an arbitrary axisymmetric surface of the body, the dependence
$\boldsymbol{r}(\boldsymbol{\gamma})$ can be represented as

$\boldsymbol{r}(\boldsymbol{\gamma})=(\chi_{1}(\gamma_{3})\gamma_{1},\chi_{1}(\gamma_{3})\gamma_{2},\chi_{2}(\gamma_{3})),$

where $\chi_{1}(\gamma_{3})$ and $\chi_{2}(\gamma_{3})$ are arbitrary
functions related to each other by [14]

$\frac{d\chi_{2}}{d\gamma_{3}}=\chi_{1}-\frac{1-\gamma_{3}^{2}}{\gamma_{3}}\frac{d\chi_{1}}{d\gamma_{3}}.$

This relation is obtained from the condition
$(\dot{\boldsymbol{r}},\boldsymbol{\gamma})=0$.

The no slip condition at the point of contact is described by a constraint of
the form

$\boldsymbol{f}=\boldsymbol{v}+\boldsymbol{\omega}\times\boldsymbol{r}={\boldsymbol{0}}.$
(3)

######  Remark.

Two projections of equation (3), $(\boldsymbol{f},\boldsymbol{\alpha})$ and
$(\boldsymbol{f},\boldsymbol{\beta})$, are nonholonomic constraints, and the
third, $(\boldsymbol{f},\boldsymbol{\gamma})$, is the time derivative of
equation (2).

The condition that there be no spinning about the vertical is described by
another nonholonomic constraint:

$g=(\boldsymbol{\omega},\boldsymbol{\gamma})=0.$ (4)

The kinetic energy $T$ and the potential energy $U$ of the system have the
form

$T=\frac{1}{2}(\boldsymbol{\omega},{\bf
I}\boldsymbol{\omega})+\frac{1}{2}m(\boldsymbol{v},\boldsymbol{v}),\quad
U=-m{\rm g}(\boldsymbol{r},\boldsymbol{\gamma}),$

where ${\bf I}=\operatorname{diag}(i_{1},i_{1},i_{3})$ is the principal
central tensor of inertia of the body and ${\rm g}$ is the free-fall
acceleration.

######  Remark.

Note that, due to the axial symmetry, the potential energy $U$ of the system
depends only on the inclination of the symmetry axis relative to the vertical,
that is, the component $\gamma_{3}$.

We write the equations of motion of the system in the absence of external
forces as Lagrange equations in quasi-velocities with constraints [48, 49]

$\displaystyle\dfrac{\rm d}{{\rm
d}t}\left(\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{v}}\right)+\boldsymbol{\omega}\times\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{v}}=\left(\dfrac{\partial\boldsymbol{f}}{\partial\boldsymbol{v}}\right)^{\\!\\!\top}\boldsymbol{\lambda},$
(5) $\displaystyle\dfrac{\rm d}{{\rm
d}t}\\!\left(\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{\omega}}\right)\\!+\boldsymbol{\omega}\times\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{\omega}}+\boldsymbol{v}\times\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{v}}+\boldsymbol{\gamma}\times\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{\gamma}}=\left(\dfrac{\partial\boldsymbol{f}}{\partial\boldsymbol{\omega}}\right)^{\\!\\!\top}\boldsymbol{\lambda}+\left(\dfrac{\partial
g}{\partial\boldsymbol{\omega}}\right)\mu,$

where $\mathcal{L}=T-U$ is the Lagrangian function of the system, and
$\boldsymbol{\lambda}$, $\mu$ are the undetermined multipliers corresponding
to the reactions of the constraints (3) and (4), respectively. The
undetermined multipliers can be found from the common solution of equations
(5) and the time derivatives of the constraints (3) and (4):

$\begin{gathered}\boldsymbol{\lambda}=-m((\boldsymbol{\omega}\times\boldsymbol{r})^{\boldsymbol{\cdot}}+\boldsymbol{\omega}\times(\boldsymbol{\omega}\times\boldsymbol{r})),\\\
\mu=\frac{(\boldsymbol{\omega}\times{\bf
J}\boldsymbol{\omega}+m\boldsymbol{r}\times(\boldsymbol{\omega}\times\dot{\boldsymbol{r}})+m{\rm
g}\boldsymbol{\gamma}\times\boldsymbol{r},{\bf
J}^{-1}\boldsymbol{\gamma})}{(\boldsymbol{\gamma},{\bf
J}^{-1}\boldsymbol{\gamma})},\end{gathered}$

where ${\bf J}={\bf I}+m\boldsymbol{r}^{2}{\bf
E}_{3}-m\boldsymbol{r}\otimes\boldsymbol{r}$ is the tensor of inertia of the
body relative to the point of contact. Here ${\bf E}_{3}$ is the $3\times 3$
identity matrix and the symbol $\otimes$ stands for a tensor multiplication of
the vectors which associates to the pair of vectors
$\boldsymbol{a},\,\boldsymbol{b}$ the matrix ${\bf A}$ with components
$A_{i,j}=a_{i}b_{j}$.

After finding the undetermined multipliers and eliminating the velocity
$\boldsymbol{v}$ (using the constraint (3)), the equations of motion (5)
reduce to one vector equation

${\bf J}\dot{\boldsymbol{\omega}}+\boldsymbol{\omega}\times{\bf
J}\boldsymbol{\omega}+m\boldsymbol{r}\times(\boldsymbol{\omega}\times\dot{\boldsymbol{r}})+m{\rm
g}\boldsymbol{\gamma}\times\boldsymbol{r}-\mu\boldsymbol{\gamma}=0.$ (6)

Adding the kinematic equation of evolution of the vector $\boldsymbol{\gamma}$
to equation (6), we obtain the following closed reduced system of equations:

$\displaystyle{\bf
J}\dot{\boldsymbol{\omega}}+\boldsymbol{\omega}\\!\times\\!{\bf
J}\boldsymbol{\omega}+m\boldsymbol{r}\\!\times\\!(\boldsymbol{\omega}\\!\times\\!\dot{\boldsymbol{r}})+m{\rm
g}\boldsymbol{\gamma}\\!\times\\!\boldsymbol{r}-\mu\boldsymbol{\gamma}=0,$ (7)
$\displaystyle\dot{\boldsymbol{\gamma}}=\boldsymbol{\gamma}\times\boldsymbol{\omega}.$

By “construction”, this system admits a first integral of motion defined by
the constraint (4)
$(\boldsymbol{\omega},\boldsymbol{\gamma})=\operatorname{const}$. The flow
defined by equations (7) and restricted to the zero level set of this integral
describes the evolution of the variables $\boldsymbol{\omega}$ and
$\boldsymbol{\gamma}$ in the problem of a body of revolution rolling without
slipping or spinning on a horizontal plane.

To reconstruct the velocity of the center of mass, use should be made of the
constraint equation (3). The coordinate $z_{c}$ is defined from equation (2),
and the remaining variables, from the quadratures (1).

###  2.3 Invariants and symmetries

The equations of motion (7) admit the following integrals of motion:

  1. –

the energy integral

$\mathcal{E}=\frac{1}{2}({\bf J}\boldsymbol{\omega},\boldsymbol{\omega})-m{\rm
g}(\boldsymbol{r},\boldsymbol{\gamma});$

  2. –

the geometric integral

$\mathcal{F}_{0}=(\boldsymbol{\gamma},\boldsymbol{\gamma})=1;$

  3. –

the no-spin constraint (4)

$\mathcal{F}_{1}=(\boldsymbol{\omega},\boldsymbol{\gamma})=0;$

  4. –

the additional integral linear in the angular velocity

$\mathcal{F}_{2}=J(\gamma_{3})\omega_{3},\,\,J(\gamma_{3})=\sqrt{i_{1}\gamma_{3}^{2}+i_{3}(1-\gamma_{3}^{2})+m(\boldsymbol{r},\boldsymbol{\gamma})^{2}}.$

######  Remark.

The system (7) admits the energy integral $\mathcal{E}$ and the additional
integral $\mathcal{F}_{2}$ only on the zero level set of the integral
$\mathcal{F}_{1}$.

######  Remark.

We note that in the presence of the constraint (4) the additional integral is
expressed in terms of elementary functions for an arbitrary body of
revolution, whereas in the classical nonholonomic case the additional integral
is written in elementary functions only for the spherical surface of the body.

In addition, equations (7) also admit the invariant measure $\rho
d\boldsymbol{\omega}d\boldsymbol{\gamma}$ with density [14]

$\rho(\gamma_{3})=(i_{1}+m\boldsymbol{r}^{2})J(\gamma_{3})$

and possess the symmetry field

$\boldsymbol{\varsigma}=\omega_{1}\frac{\partial}{\partial\omega_{2}}-\omega_{2}\frac{\partial}{\partial\omega_{1}}+\gamma_{1}\frac{\partial}{\partial\gamma_{2}}-\gamma_{2}\frac{\partial}{\partial\gamma_{1}},$
(8)

which corresponds to invariance of the system under rotations about the axis
of dynamical symmetry of the body.

Thus, the problem of a body of revolution rolling on a plane without slipping
or spinning about the vertical is integrable by the Euler – Jacobi theorem and
reduces to quadratures [4].

###  2.4 Reduction to a system with one degree of freedom

The existence of four integrals of motion allows us to reduce the system (7)
to a system with one degree of freedom.

To perform a reduction, we parameterize the vector $\boldsymbol{\gamma}$ by
the nutation angle $\vartheta$ and the angle of proper rotation $\varphi$ in
the form

$\gamma_{1}=\sin\vartheta\sin\varphi,\quad\gamma_{2}=\sin\vartheta\cos\varphi,\quad\gamma_{3}=\cos\vartheta.$
(9)

It is easy to show that on the fixed level set of the first integrals

$\mathcal{F}_{0}=1,\quad\mathcal{F}_{1}=0,\quad\mathcal{F}_{2}=k$

the system of equations (7) reduces to a system with one degree of freedom

$\displaystyle\dot{\vartheta}=p_{\vartheta},\qquad\dot{p}_{\vartheta}=\frac{G(\vartheta)}{B(\vartheta)},$
(10) $\displaystyle
G(\vartheta)=\left(\frac{k^{2}\cot\vartheta}{\sin^{2}\\!\vartheta}-\frac{B^{\prime}(\vartheta)}{2}p_{\vartheta}^{2}-\frac{\partial
U(\vartheta)}{\partial\vartheta}\right),$

where $B(\vartheta)=i_{1}+m\boldsymbol{r}^{2}$, $U(\vartheta)$ is the
potential energy of the system. The system (10) admits the energy integral

$\mathcal{E}=\frac{B(\vartheta)p_{\vartheta}^{2}}{2}+\frac{k^{2}}{2\sin^{2}\\!\vartheta}+U(\vartheta).$

In this case, the quadrature for angle $\varphi$ has the form

$\dot{\varphi}=\frac{\omega_{3}}{\sin^{2}\\!\vartheta},$

where

$\omega_{3}=\frac{k}{J(\gamma_{3})}=\frac{k}{J(\vartheta)}.$

The values of the variables $\omega_{1}$ and $\omega_{2}$ can be reconstructed
from the relations

$\displaystyle\omega_{1}=\dot{\vartheta}\cos\varphi-\omega_{3}\cot\vartheta\sin\varphi,$
(11)
$\displaystyle\omega_{2}=-\dot{\vartheta}\sin\varphi-\omega_{3}\cot\vartheta\cos\varphi.$

Next, we consider the case of an ellipsoid of revolution and perform a
complete classification of possible motions of this system.

##  3 Explicit form of the equations of motion
for an ellipsoid of revolution

We now turn our attention to a specific body of revolution, an ellipsoid of
revolution with semiaxes $b_{1}$ and $b_{3}$ and with a center of mass
displaced along the symmetry axis by distance $a$ (see Fig. 2).

Figure 2:  (a) Diagrammatic representation of an ellipsoid on a plane. The
plane $\rm\Sigma$ is a vertical plane passing through the symmetry axis of the
ellipsoid. The red elliptic curve is the intersection of the ellipsoid and the
plane $\rm\Sigma$. (b) Representation of the ellipsoid on the plane in a
meridional section.

The radius vector of the point of contact of the ellipsoid can be represented
as

$\boldsymbol{r}=-\frac{{\bf
B}\boldsymbol{\gamma}}{\sqrt{(\boldsymbol{\gamma},{\bf
B}\boldsymbol{\gamma})}}-a\boldsymbol{e}_{3},\qquad{\bf
B}=\operatorname{diag}(b_{1}^{2},b_{1}^{2},b_{3}^{2}).$

######  Remark.

In this case, the explicit expressions for the functions
$\chi_{1}(\gamma_{3})$ and $\chi_{2}(\gamma_{3})$ for the case of an ellipsoid
of revolution have the form

$\displaystyle\chi_{1}(\gamma_{3})=-\frac{b_{1}^{2}}{\sqrt{b_{1}^{2}(1-\gamma_{3}^{2})+b_{3}^{2}\gamma_{3}^{2}}},$
$\displaystyle\chi_{2}(\gamma_{3})=-\frac{b_{3}^{2}\gamma_{3}}{\sqrt{b_{1}^{2}(1-\gamma_{3}^{2})+b_{3}^{2}\gamma_{3}^{2}}}-a.$

The equations of motion (7) written in the new variables depend on many
parameters. In order to decrease their number, we introduce dimensionless
parameters: by choosing the length of semiaxis $b_{3}$ as a unit of length,
the mass of the ellipsoid $m$ as a unit of mass, and the quantity
$\sqrt{b_{3}/{\rm g}}$ as a unit of time, we obtain four dimensionless
parameters

$\displaystyle\alpha=\frac{a}{b_{3}}\in[0,1],\quad\beta=\frac{b_{1}}{b_{3}}\in(0,\infty),$
$\displaystyle\nu=\frac{i_{3}}{i_{1}}\in(0,2],\quad\eta=\frac{mb_{3}^{2}}{i_{1}}\in(0,\infty)$

and two dimensionless constants of the integrals

$\varepsilon=\frac{\mathcal{E}}{m{\rm
g}b_{3}},\quad\kappa^{2}=\frac{k^{2}}{m{\rm g}b_{3}}.$

The values $\beta\in(0,1)$ correspond to an “prolate” ellipsoid, and
$\beta\in(1,\infty)$, to an “oblate” one. In the limiting case $\beta=0$ we
obtain a thin rod, and for $\beta\rightarrow\infty$ we obtain a thin disk.
However, these cases correspond to the motion of a body with a sharp edge,
which we do not consider in this paper.

After transformation to the new variables and nondimensionalization for the
case of an ellipsoid of revolution, the potential energy $U(\vartheta)$ takes
the form

$\displaystyle U(\vartheta)=\alpha\cos\vartheta+Z(\vartheta),$ $\displaystyle
Z(\vartheta)=\sqrt{\beta^{2}\sin^{2}\\!\vartheta+\cos^{2}\\!\vartheta}>0\quad\forall\,\vartheta\in[0,\pi],$

and the equations of motion (10) can be written as

$\displaystyle\dot{\vartheta}=p_{\vartheta},$
$\displaystyle\dot{p}_{\vartheta}=\frac{1}{B(\vartheta)}\left(\frac{\kappa^{2}\cot\vartheta}{\sin^{2}\\!\vartheta}-\frac{B^{\prime}(\vartheta)}{2}p_{\vartheta}^{2}+\alpha\sin\vartheta+\frac{(1-\beta^{2})\sin\vartheta\cos\vartheta}{Z(\vartheta)}\right),$

where

$B(\vartheta)=\frac{1}{\eta}+\frac{1}{Z(\vartheta)^{2}}\left(\beta^{4}\sin^{2}\\!\vartheta+(\alpha
Z(\vartheta)-\cos\vartheta)^{2}\right).$

We note here that the physical meaning of $Z(\vartheta)$ is the height of the
geometric center of the ellipsoid above the supporting plane.

Next, we perform a complete classification of possible motions of the
resulting system.

##  4 Bifurcation analysis

###  4.1 General definitions

To classify the motions of the system, we use the topological approach
described in [50, 51]. For this we consider the system (7) in the redundant
coordinates $\boldsymbol{z}=(\boldsymbol{\omega},\boldsymbol{\gamma})$, which
defines the phase flow in the space $\mathcal{M}_{6}=\\{\boldsymbol{z}\\}$.
This system possesses four integrals of motion $\mathcal{F}_{0}$,
$\mathcal{F}_{1}$, $\mathcal{F}_{2}$, and $\mathcal{E}$. We restrict the
system (7) to the common level set of the integrals $\mathcal{F}_{0}$ and
$\mathcal{F}_{1}$

$\mathcal{M}_{4}=\\{\boldsymbol{z}\,|\,\mathcal{F}_{0}=1,\,\mathcal{F}_{1}=0\\}\subset\mathcal{M}_{6}.$

Let us define the integral map of this system
${\rm\Phi}\,:\,\mathcal{M}_{4}\rightarrow\mathbb{R}^{2}$

$\boldsymbol{z}\rightarrow{\rm\Phi}(\boldsymbol{z})=(\mathcal{F}_{2}(\boldsymbol{z}),\mathcal{E}(\boldsymbol{z}))\in\mathbb{R}^{2},\quad\boldsymbol{z}\in\mathcal{M}_{4}.$

######  Remark.

As local coordinates on $\mathcal{M}_{4}$ we can choose, for example, the
quantities $(\vartheta,\,\varphi,\,\omega_{1},\,\omega_{2})$. However, some
solutions of interest lie on the singularities of the chosen coordinates
($\vartheta=0,\pi$). Therefore, it is more convenient to perform calculations
in the redundant coordinates.

The complete image of the phase space ${\rm\Phi}(\mathcal{M}_{4})$ defines the
region of possible motions $($RPM$)$ of the system on the plane of first
integrals $(\kappa,\varepsilon)$. Each point on the plane of first integrals
$(\kappa,\varepsilon)$ corresponds to the integral manifold

$\mathcal{M}_{\kappa,\varepsilon}\\!=\\!\\{\boldsymbol{z}\,|\,\mathcal{F}_{0}=1,\mathcal{F}_{1}=0,\mathcal{F}_{2}(\boldsymbol{z})=\kappa,\mathcal{E}(\boldsymbol{z})=\varepsilon\\}\\!\subset\\!\mathcal{M}_{6},$

which can contain several connectedness components (generally two-dimensional
tori).

Changes in the values of the integrals $\kappa$ and $\varepsilon$ lead to
bifurcations of the integral manifolds. These bifurcations correspond to a
critical set of the integral map which is defined as follows:

$S=\\{\boldsymbol{z}\in\mathcal{M}_{4}\,|\,\operatorname{rank}{\mathrm{d}}{\rm\Phi}<2\\}.$

The image of the critical set $S$ in the space of first integrals
${\rm\Phi}(S)$ and the region of possible motions ${\rm\Phi}(\mathcal{M}_{4})$
form a bifurcation diagram of the integrable system.

The set $S$ is a union of two subsets:

$\displaystyle
S_{0}=\\{\boldsymbol{z}\in\mathcal{M}_{4}\,|\,\operatorname{rank}{\mathrm{d}}{\rm\Phi}=0\\},$
$\displaystyle
S_{1}=\\{\boldsymbol{z}\in\mathcal{M}_{4}\,|\,\operatorname{rank}{\mathrm{d}}{\rm\Phi}=1\\}.$

In the case at hand the subset $S_{0}$ consists of separate points, and the
subset $S_{1}$ is a one-parameter family of closed curves. The points
belonging to $S_{0}$ are the fixed points of the system (7), and the set
$S_{1}$ is filled by periodic solutions of the system (7), which correspond to
fixed points of the reduced system (10).

###  4.2 Critical sets

To find the sets $S_{1}$ or $S_{0}$ in redundant coordinates, it is necessary
to solve the system of equations given by the equations

$\left.\operatorname{rank}\frac{\partial(\mathcal{F}_{0},\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{E})}{\partial\boldsymbol{z}}\right|_{\mathcal{M}_{4}}=3$
(12)

or

$\left.\operatorname{rank}\frac{\partial(\mathcal{F}_{0},\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{E})}{\partial\boldsymbol{z}}\right|_{\mathcal{M}_{4}}=2,$
(13)

respectively. This can be done, for example, using the method of undetermined
multipliers (see, e.g., [52]). We will not present here this procedure in
explicit form since it is standard and the intermediate expressions are
cumbersome.

Next, we consider each of the critical subsets separately.

####  4.2.1 The critical subset $S_{1}$

Analysis of equations (12) shows that they can have different solutions,
depending on whether or not the parameter $\alpha$ is equal to zero. We first
consider the general case $\alpha\neq 0$.

Solving the system of equations given by (12), we can prove the following
proposition.

######  Proposition.

In the case $\alpha\neq 0$ the critical subset $S_{1}$ of the integral map
$\Phi$ has the form

$S_{1}=S_{1}^{\alpha}=\\{\boldsymbol{z}\,|\,\boldsymbol{\omega}=\omega_{0}(\vartheta)\boldsymbol{\gamma}\times(\boldsymbol{e}_{3}\times\boldsymbol{\gamma}),\,\boldsymbol{\gamma}=\boldsymbol{\gamma}(\vartheta,\varphi)\\},$

where the dependence $\boldsymbol{\gamma}(\vartheta,\varphi)$ is defined by
the expressions (9), $\omega_{0}(\vartheta)$ has the form

$\omega_{0}^{2}=\frac{-\eta\sin^{2}\\!\vartheta\left(\cos\vartheta(1-\beta^{2})+\alpha
Z(\vartheta)\right)}{\cos\vartheta
Z(\vartheta)\left(\eta(Z(\vartheta)+\alpha\cos\vartheta)^{2}+\cos^{2}\\!\vartheta+\nu\sin^{2}\\!\vartheta\right)},$
(14)

$\varphi\in[0,2\pi)$, and the range of angle $\vartheta$ is determined from
the condition that the expression (14) be nonnegative for $\omega_{0}^{2}$.

The critical set $S_{1}^{\alpha}$ is filled by periodic trajectories for which
the dependences $\vartheta(t)$ and $\varphi(t)$ have the form

$\vartheta=\vartheta_{0}=\operatorname{const},\quad\varphi=\varphi_{0}+\frac{\omega_{0}}{\sin\vartheta_{0}}t,$
(15)

where $\varphi_{0}\in[0,2\pi)$ defines the initial angle of proper rotation.
Due to the axial symmetry of the system, without loss of generality, one can
set $\varphi_{0}=0$. As is customary, we will call the solutions (15)
permanent rotations. In the reduced system (10) they correspond to a one-
parameter family of fixed points

$\vartheta=\vartheta_{0}=\operatorname{const},\quad p_{\vartheta}=0.$ (16)

In the problem we consider here, under permanent rotations, the trajectory of
the ellipsoid in the fixed coordinate system is closed and the following
proposition holds.

######  Proposition.

The fixed points (16) correspond to the rolling of an ellipsoid in a circle
with a constant inclination angle of the symmetry axis, $\vartheta_{0}$. The
point of contact traces out a circle of radius
$\rho_{p}=\frac{\beta^{2}}{Z(\vartheta_{0})}\tan\vartheta_{0}$, and the center
of mass traces out a circle of radius
$\rho_{c}=Z(\vartheta_{0})\tan\vartheta_{0}+\alpha\sin\vartheta_{0}$.

###### Proof.

Let us parameterize the rotation matrix ${\bf Q}$ by Euler’s angles: the
nutation angle $\vartheta$, the precession angle $\psi$ and the angle of
proper rotation $\varphi$ (we have already introduced two of the three angles
above for parameterization of the vector $\boldsymbol{\gamma}$) [49]

${\bf
Q}\\!=\\!{\begin{pmatrix}\cos\psi\cos\varphi-\cos\vartheta\sin\psi\sin\varphi\\!\\!\\!&\\!\\!\\!\sin\psi\cos\varphi+\cos\vartheta\cos\psi\sin\varphi\\!\\!\\!&\\!\\!\\!\sin\varphi\sin\vartheta\\\
-\cos\psi\sin\varphi-\cos\vartheta\sin\psi\cos\varphi\quad\\!\\!\\!&\\!\\!\\!-\sin\psi\sin\varphi+\cos\vartheta\cos\psi\cos\varphi\quad\\!\\!\\!&\\!\\!\\!\cos\varphi\sin\vartheta\\\
\sin\vartheta\sin\psi\\!\\!&\\!\\!-\sin\vartheta\cos\psi\\!\\!&\\!\\!\cos\vartheta\end{pmatrix}}.$

The quadrature for angle $\psi$, with (1) and (11) taken into account, has the
form

$\dot{\psi}=-\frac{\kappa\cos\vartheta}{J(\vartheta)\sin^{2}\\!\vartheta}.$
(17)

Thus, for the solution (15) the dependence $\psi(t)$ has the form

$\psi=\psi_{0}+\omega_{\psi}t,\quad\omega_{\psi}=-\frac{\kappa\cos\vartheta_{0}}{J(\vartheta_{0})\sin^{2}\\!\vartheta_{0}}=\operatorname{const},$

where $\psi_{0}$ is the initial value of angle $\psi$, which without loss of
generality can be set equal to zero $\psi_{0}=0$.

In accordance with (15), angle $\varphi$ also changes linearly in time

$\varphi=\varphi_{0}+\omega_{\varphi}t,\quad\omega_{\varphi}=\frac{\kappa}{J(\vartheta_{0})\sin^{2}\\!\vartheta_{0}}=\operatorname{const},$

where the initial value of angle $\varphi_{0}$ will also be assumed to be
zero.

Substituting the resulting expressions into the quadratures (1) and
integrating them over time, we obtain

$x_{c}=\rho_{c}\sin\psi,\quad y_{c}=-\rho_{c}\cos\psi,\quad
z_{c}=\alpha\cos\vartheta_{0}+Z(\vartheta_{0}),$

where $\rho_{c}=Z(\vartheta_{0})\tan\vartheta_{0}+\alpha\sin\vartheta_{0}$.

Knowing the trajectory of the center of mass $\boldsymbol{r}_{c}(t)$, we find
the trajectory of the point of contact from the relation

$\boldsymbol{r}_{p}=\boldsymbol{r}_{c}+{\bf Q}\boldsymbol{r},$

where $\boldsymbol{r}_{p}=(x_{p},y_{p},0)$ is the radius vector of the point
of contact referred to the fixed coordinate system. This yields

$x_{p}=x_{c}+(\boldsymbol{r},\boldsymbol{\alpha})=\rho_{p}\sin\psi,\quad
y_{p}=y_{c}+(\boldsymbol{r},\boldsymbol{\beta})=-\rho_{p}\cos\psi,$

where $\rho_{p}=\frac{\beta^{2}}{Z(\vartheta_{0})}\tan\vartheta_{0}$.

∎

From the resulting expressions for the radius of the circle of the trajectory
of the contact point and the center of mass it can be seen that, as
$\vartheta_{0}\rightarrow\pi/2$, $\rho_{c}\rightarrow\infty$ and
$\rho_{p}\rightarrow\infty$. That is, as the symmetry axis approaches the
horizontal position, the radius of the circle increases without bound. Figure
3 shows a schematic representation of permanent rotations of an ellipsoid of
revolution on a plane.

Figure 3: Example of permanent rotations of an ellipsoid

Note that, for fixed values of the parameters $\alpha$ and $\beta$, permanent
rotations are not possible for all inclination angles of the symmetry axis of
the ellipsoid $\vartheta_{0}$. Possible inclination angles are defined by the
inequality $\omega_{0}^{2}(\vartheta)\geqslant 0$. Figure 4 shows the regions
on the plane $(\beta^{2},\vartheta)$ that are given by this inequality for a
fixed value of $\alpha$.

Figure 4: Regions of positive definiteness of the expression (14) on the
parameter plane $(\beta^{2},\vartheta)$ for a fixed value of $\alpha$.

######  Remark.

In [53] conditions are found for the parameters of an ellipsoid under which
permanent rotations may occur at any point of the ellipsoid’s surface, and in
[54] the stability of such rotations is analyzed.

Analyzing the expression (14) and Fig. 4, we can formulate the following
statements.

######  Statement 1.

Permanent rotations under which the center of mass of the ellipsoid lies above
its geometric center exist only for a fairly oblate ellipsoid
$(\beta^{2}>1+\alpha)$.

######  Statement 2.

Permanent rotations under which the center of mass of the ellipsoid lies below
its geometric center exist regardless of the form of the ellipsoid. However,
for a fairly prolate ellipsoid $(\beta^{2}<1-\alpha)$ such rotations are
possible only for a fairly large inclination angle of its symmetry axis.

In addition, it follows from the analysis of the expression (14) and Fig. 4
that there exist inclined equilibrium positions of the ellipsoid. These
equilibrium positions are defined from the condition $\omega_{0}^{2}=0$. The
inclination angle for such equilibrium positions depends on the parameters of
the ellipsoid as follows:

$\vartheta_{*}=\arccos\frac{\alpha\beta}{\sqrt{(1-\beta^{2})(1+\alpha^{2}-\beta^{2})}}.$
(18)

In Fig. 4, these equilibrium positions correspond to the boundary of the
region $\omega_{0}^{2}\geqslant 0$. These positions exist only under some
restrictions on the mass-geometric parameters of the ellipsoid. These
restrictions have a simple mechanical interpretation.

######  Statement 3.

The inclined equilibrium positions of an ellipsoid of revolution with
axisymmetric mass distribution exist only when its center of mass lies between
the centers of curvature of its vertices.

In the case $\alpha=0$ (without displacement of the center of mass),
additional symmetry leads to an increase in the critical set $S_{1}$.

######  Proposition.

In the case $\alpha=0$ the critical subset $S_{1}$ of the integral map
${\rm\Phi}$ is a union of two subsets

$\displaystyle S_{1}=S_{1}^{\alpha}\cup S_{1}^{0},$ $\displaystyle
S_{1}^{0}=\\{(\boldsymbol{\omega},\boldsymbol{\gamma}(\vartheta,\varphi))\,|\,\boldsymbol{\omega}=(0,0,\omega_{3}),\,\vartheta=\frac{\pi}{2},$
$\displaystyle\hskip
150.79968pt\varphi\in[0,2\pi),\,\omega_{3}\in\mathbb{R}\\}.$

The critical subset $S_{1}^{0}$ is filled by periodic trajectories

$\omega_{3}=\omega_{0}=\operatorname{const},\quad\vartheta=\frac{\pi}{2},\quad\varphi=\varphi_{0}+\omega_{0}t.$
(19)

In the reduced system (10) all solutions (19) correspond to one fixed point

$\vartheta=\frac{\pi}{2},\quad p_{\vartheta}=0.$

In the fixed coordinate system, the solutions (19) correspond to the rolling
along the ellipsoid’s equator in a straight line with an arbitrary velocity.

####  4.2.2 The critical subset $S_{0}$

We now consider the critical subset $S_{0}$. Solving the system of equations
given by equation (13), we can prove the following proposition.

######  Proposition.

The critical subset $S_{0}$ of the integral map $\Phi$ consists of two points

$S_{0}=S_{0}^{\pm}=\\{(\boldsymbol{\omega},\boldsymbol{\gamma})\,|\,\boldsymbol{\omega}=0,\,\boldsymbol{\gamma}=(0,0,\pm
1)\\}.$

The critical subset $S_{0}$ is filled by fixed points of the system (7). These
points correspond to the equilibrium positions of the ellipsoid in which it
rests on one of its vertices. In Fig. 4, the fixed points $S_{0}$ correspond
to the horizontal straight lines $\vartheta=0$ and $\vartheta=\pi$.

######  Remark.

The fixed points of the complete system (7) $S_{0}$ correspond to fixed points
of the reduced system (10). This is due to the fact that such equilibrium
positions lie on singular zero-dimensional orbits of the action of the
rotation group (8).

###  4.3 Example of a bifurcation diagram

Consider an example of constructing a bifurcation diagram which is an image of
the critical set $S$. We present the explicit form of bifurcation curves and
points corresponding to the image ${\rm\Phi}(S)$:

  1. –

the image of the subset $S_{1}^{\alpha}$ consists of curves on the plane of
first integrals $(\kappa,\varepsilon)$, which can be represented as

$\sigma_{\vartheta}:\left\\{\begin{aligned}
&\varepsilon=\frac{3Z(\vartheta_{0})^{2}-1}{2Z(\vartheta_{0})}+\alpha\frac{3\cos^{2}\\!\vartheta_{0}-1}{2\cos\vartheta_{0}},\\\
&\kappa=\pm\sin^{2}\\!\vartheta_{0}\sqrt{\frac{\beta^{2}-1}{Z(\vartheta_{0})}-\frac{\alpha}{\cos\vartheta_{0}}},\end{aligned}\right.$
(20)

where

$\displaystyle\vartheta_{0}\in(\pi/2,\vartheta_{*})\quad\text{for}\quad\beta^{2}<1-\alpha,$
$\displaystyle\vartheta_{0}\in(\pi/2,\pi)\quad\text{for}\quad
1-\alpha<\beta^{2}<1+\alpha,$
$\displaystyle\vartheta_{0}\in(0,\vartheta_{*}]\cup(\pi/2,\pi)\quad\text{for}\quad\beta^{2}>1+\alpha;$

  2. –

the image of the subset $S_{1}^{0}$ is a parabola

$\sigma_{\pi/2}:\,\varepsilon=\frac{\kappa^{2}}{2}+\beta;$ (21)

  3. –

the image of the subset $S_{0}$ on the plane of first integrals
$(\kappa,\varepsilon)$ consists of a pair of point

$\sigma_{0}\,:\,(0,1+\alpha),\qquad\sigma_{\pi}\,:\,(0,1-\alpha).$ (22)

Depending on the values of the parameters $\alpha$ and $\beta$, these points
can be isolated. It can be seen from Fig. 4 that the point $\sigma_{0}$ is
isolated when $\beta^{2}<1+\alpha$, and the point $\sigma_{\pi}$, when
$\beta^{2}<1-\alpha$. We note that in the case $\alpha=0$ these points lie on
the same level set of the integral $\kappa=0$, $\varepsilon=1$, but do not
coincide in $\mathcal{M}_{6}$.

Figure 5 shows an example of constructing a bifurcation diagram on the plane
of first integrals $(\kappa,\varepsilon)$ for the general case $\alpha\neq 0$.
This example of a bifurcation diagram is given for the case of an “oblate”
ellipsoid $\beta^{2}>1$ with a relatively small displacement of the center of
mass. In this section, we explain the notation used in the diagram, and in the
following sections we will use them without repeating our explanations.

Figure 5: Example of a bifurcation diagram. This type of diagram is
characteristic of an “oblate” ellipsoid with a fairly small displacement of
the center of mass $\alpha$.

In the diagram, gray denotes the region of possible motions of the system. In
the light gray region, the integral manifold is two-dimensional torus
$\mathbb{T}^{2}$. Dark gray denotes the region which corresponds to two
connectedness components of the integral manifold (two two-dimensional tori
$2\mathbb{T}^{2}$). These two connectedness components are glued together to
each other along the curve $\sigma_{u}$. In the white region, no motions of
the system are possible ($\varnothing$).

The heavy solid and dotted lines in the diagram show parts of the bifurcation
curves (20) (and (21) for the case $\alpha=0$). The solid curves correspond to
stable fixed points, and the dotted curves indicate unstable ones. The curve
(20) is partitioned by critical values of angle $\vartheta_{0}$ into separate
segments with the same type of stability, which are denoted in the diagram by
different symbols. The curve $\sigma_{s\pi}$ corresponds to
$\vartheta_{0}\in(\pi/2,\pi)$, the curve $\sigma_{s0}$ corresponds to
$\vartheta_{0}\in(0,\vartheta_{c})$, and the curve $\sigma_{u}$, to
$\vartheta_{0}\in(\vartheta_{c},\vartheta_{*}]$. Here $\vartheta_{c}$ is the
critical value of angle $\vartheta_{0}$. At this value,
$\vartheta_{0}=\vartheta_{c}$ (which corresponds to the point
$(\kappa_{c},\varepsilon_{c})$ in the diagram), one can observe a change in
the type of stability of fixed points corresponding to the curve
$\sigma_{\vartheta}$. The fixed points (22) corresponding to equilibrium
positions are shown as colored black points (in the case of their stability)
or punctured points (in the case of their instability).

As the parameters of the system $\alpha$, $\beta$ change, the critical subsets
$S_{1}$ and $S_{0}$ can change, resulting in a qualitative change in the type
of the diagram. The bifurcations of critical points described above occur
simultaneously with this change. In order to determine possible types of
bifurcation diagrams, we investigate the bifurcations of critical points in
what follows.

###  4.4 Bifurcations of critical sets and types of bifurcation diagrams

####  4.4.1 Bifurcations of the critical subset $S_{0}$

Consider the bifurcations of the critical subset $S_{0}$.

When $\beta^{2}<1-\alpha$ (see Fig. 4), the subset $S_{0}$ consists of two
isolated points. When $\beta^{2}=1-\alpha$, a bifurcation occurs, resulting in
the component $S_{0}^{-}$ merging with the component $S_{1}^{\alpha}$.

Next, when $1-\alpha<\beta^{2}<1+\alpha$, one of the points of the subset
$S_{0}$ is isolated ($S_{0}^{+}$) and the other is not isolated ($S_{0}^{-}$).

When $\beta^{2}=1+\alpha$, a second bifurcation occurs during which the
component $S_{1}^{\alpha}$ is born from $S_{0}^{+}$, and both points of the
subset $S_{0}$ are not isolated now.

We note that, simultaneously with the bifurcations of critical sets described
above, changes in the stability of the fixed points occur. To illustrate this
fact, we analyze the linear stability of the fixed points

$\vartheta=0,\quad p_{\vartheta}=0\qquad\text{and}\qquad\vartheta=\pi,\quad
p_{\vartheta}=0$ (23)

of the reduced system (10).

To do so, we represent the system (10) as

$\dot{\boldsymbol{q}}={\bf f}(\boldsymbol{q}),$

where $\boldsymbol{q}=(\vartheta,p_{\vartheta})$, ${\bf f}(\boldsymbol{q})$ is
the vector function whose elements depend on the components of
$\boldsymbol{q}$. Introducing the perturbation vector
$\boldsymbol{\tilde{q}}=\boldsymbol{q}-\boldsymbol{q}_{0}$, where
$\boldsymbol{q}_{0}=(\vartheta_{0},0)$ is the equilibrium position, and
representing the equation for $\boldsymbol{\tilde{q}}$ as a series expansion
near the equilibrium position, we obtain

$\dot{\tilde{\boldsymbol{q}}}={\bf\Lambda}\boldsymbol{\tilde{q}},\quad{\bf\Lambda}=\left.\frac{\partial{\bf
f}(\boldsymbol{q})}{\partial\boldsymbol{q}}\right|_{\boldsymbol{q}=\boldsymbol{q}_{0}}.$

The matrix $\bf\Lambda$ has the form

${\bf\Lambda}=\begin{pmatrix}0&1\\\
\frac{G_{0}^{\prime}(\vartheta_{0})}{B(\vartheta_{0})}&0\end{pmatrix},$

where the notation
$G_{0}(\vartheta_{0})=\left.G(\vartheta_{0})\right|_{p_{\vartheta}=0}$ is
used, i.e., $G_{0}^{\prime}(\vartheta_{0})$ denotes the derivative of the
function $G(\vartheta)$ with respect to $\vartheta$, calculated at the point
$\vartheta_{0}$ for $p_{\vartheta}=0$.

We write the characteristic equation

$\det({\bf\Lambda-\lambda{\bf E}_{2}})=0,$

where $\lambda$ are the eigenvalues of the matrix $\bf\Lambda$ and ${\bf
E}_{2}$ is the $2\times 2$ identity matrix. This equation reduces to the
simple quadratic equation

$\lambda^{2}-\frac{G_{0}^{\prime}(\vartheta_{0})}{B(\vartheta_{0})}=0.$

Since $B(\vartheta_{0})>0$ for any values of $\vartheta_{0}$, the eigenvalues
of $\lambda$ are defined by the value of $G_{0}^{\prime}(\vartheta_{0})$. If
$G_{0}^{\prime}(\vartheta_{0})<0$, the corresponding fixed point is a stable
center-type point. If $G_{0}^{\prime}(\vartheta_{0})>0$, then the
corresponding fixed point is an unstable saddle point.

For the lower equilibrium position ($\vartheta_{0}=\pi$, $\kappa=0$) the value
of $G_{0}^{\prime}(\vartheta_{0})$ is negative for $\beta^{2}>1-\alpha$, and
for the upper equilibrium position ($\vartheta_{0}=0$, $\kappa=0$) it is
negative for $\beta^{2}>1+\alpha$. Thus, depending on the type of stability of
the fixed points (23), the parameter plane $(\alpha,\beta^{2})$ can be divided
into three regions:

  1. (a)

$\beta^{2}<1-\alpha$. In this region, both fixed points $\sigma_{0}$,
$\sigma_{\pi}$ are unstable and isolated;

  2. (b)

$1-\alpha<\beta^{2}<1+\alpha$. In this region, the lower position of the
ellipsoid $\sigma_{\pi}$ is stable and is nonisolated, and the upper one,
$\sigma_{0}$, is unstable and isolated;

  3. (c)

$\beta^{2}>1-\alpha$. In this region, both fixed points are stable
$\sigma_{0}$, $\sigma_{\pi}$ and nonisolated.

######  Remark.

We note that the resulting boundaries of the stability of the solutions
$\sigma_{0}$ and $\sigma_{\pi}$ are defined only by the geometry of the body
and the position of the center of mass. For the stability of the vertical
position of the ellipsoid it is necessary that its center of mass
($1\pm\alpha$) be below the center of curvature ($\beta^{2}$) of the
ellipsoid’s surface at the point of contact. The center of curvature in this
case is correctly determined because the body is axisymmetric and the
principal radii of curvature at the vertices of the ellipsoid are the same.
This statement is in agreement with the results obtained in the book [55] for
an arbitrary body of revolution and in our recent paper [13] on the dynamics
of an ellipsoid on a vibrating plane.

####  4.4.2 Bifurcations of critical subsets $S_{1}$

We now turn our attention to bifurcations of the critical subset $S_{1}$,
which corresponds to the fixed points (16) of the reduced system (10). The
topological type of the subset $S_{1}$ is uniquely related to the type of the
bifurcation curve (20). The type of bifurcation curve is determined, in its
turn, by its critical points, which satisfy the system of equations

$G_{0}(\vartheta)=0,\quad\dfrac{\partial
G_{0}(\vartheta)}{\partial\vartheta}=0.$ (24)

In the bifurcation diagram, these critical points correspond to the cusp
points and tangency points of the bifurcation curves; also, the bifurcation
curves can lose smoothness at these points.

Let us analyze the bifurcations of the critical points depending on the values
of the parameters $\alpha$ and $\beta$.

The emergence of new critical points is defined by the additional condition

$\dfrac{\partial^{2}G_{0}(\vartheta)}{\partial\vartheta^{2}}=0.$ (25)

The common solution of the system of equations (24)–(25) defines the curve on
the plane $(\alpha,\beta^{2})$ which separates regions with different numbers
of critical points and hence with different types of bifurcation diagrams.

In the case $\alpha\neq 0$, the system of equations (24)–(25) gives two
solutions:

  1. 1.

$\beta^{2}=1+\alpha$. This straight line on the plane $(\alpha,\beta^{2})$
separates the regions with one cusp point (with $\beta^{2}>1+\alpha$) and
without cusp points (with $\beta^{2}<1+\alpha$). In addition, as shown in the
preceding section, when this straight line is intersected, a bifurcation curve
$\sigma_{\vartheta}$ (branches $\sigma_{s0}$ and $\sigma_{u}$) is born from
point $\sigma_{0}$, and $\sigma_{0}$ ceases to be isolated.

  2. 2.

$\beta^{2}=1-\alpha$. This straight line on the plane $(\alpha,\beta^{2})$
corresponds to the loss of smoothness of the bifurcation curve $\sigma_{s\pi}$
with $\kappa=0$, in which case the point $\sigma_{\pi}$ descends to the
bifurcation curve $\sigma_{s\pi}$ and ceases to be isolated.

As we see, the bifurcations of $S_{1}$ listed above occur at the same time as
those of $S_{0}$. This system has no other bifurcations of $S_{1}$ which do
not affect $S_{0}$.

As was shown above, when $\alpha=0$ the subset $S_{1}$ bifurcates. Next, we
consider the question of bifurcations of $S_{1}$ for constant $\alpha=0$ and
for the varying parameter $\beta$. In the case $\alpha=0$ the system of
equations (24)–(25) has a unique solution $\beta^{2}=1$. For $\beta^{2}<1$
there exists only one curve $\sigma_{\pi/2}$ on the plane of first integrals,
and the fixed points $\sigma_{0}$ and $\sigma_{\pi}$ are isolated. When
$\beta^{2}=1$, bifurcation curves $\sigma_{s0}$ and $\sigma_{s\pi}$ are born
from the points $\sigma_{0}$ and $\sigma_{\pi}$, respectively. The points
$\sigma_{0}$ and $\sigma_{\pi}$ cease to be isolated, and on the curve
$\sigma_{\pi/2}$ tangency points are born at which the curves $\sigma_{s0}$
and $\sigma_{s\pi}$ are tangent to the curve $\sigma_{\pi/2}$.

Simultaneously, the stability of the fixed point corresponding to the curve
$\sigma_{\pi/2}$ changes. Indeed, following the above-mentioned algorithm of
the linear stability analysis of fixed points, the condition for stability is
defined by the sign $G_{0}^{\prime}(\pi/2)$, whence we find that the point
(19) is stable for

$|\kappa|<\kappa_{c}=\sqrt{\frac{\beta^{2}-1}{\beta}}.$ (26)

From this, we obtain the critical value of the parameter

$\beta^{2}=1.$

Thus, when $\beta^{2}<1$ the fixed point

$\vartheta=\frac{\pi}{2},\quad p_{\vartheta}=0$

is stable for any values of $\kappa$, and for $\beta^{2}>1$ it is stable only
when condition (26) is satisfied.

####  4.4.3 Types of bifurcation diagrams

Summarizing, we note that the entire parameter space $(\alpha,\beta^{2})$ can
be divided into five regions corresponding to different types of bifurcation
diagrams (see Fig. 6). Figure 7 shows all possible types of bifurcation
diagrams for the problem under consideration. Diagrams (a)–(c) correspond to
the general case $\alpha>0$ and have been constructed for the fixed value
$\alpha=0.5$ and different values of $\beta$, and diagrams (c) and (d), for
the case $\alpha=0$.

Figure 6: Regions with different types of bifurcation diagrams on the
parameter plane $(\alpha,\beta^{2})$. Figure 7: Possible types of bifurcation
diagrams for (a)–(c) $\alpha=0.5$ and (d), (e) $\alpha=0$.

Next, we consider in more detail each of the types of diagrams. In the names
of the types of diagrams (a)–(e) we refer to Fig. 7. In addition, for each
type of diagrams we construct all possible types of phase portraits of the
reduced system (10) on different level sets of the integral $\mathcal{F}_{2}$.

####  4.4.4 Bifurcation diagram (a)

Figure 8 shows a bifurcation diagram and the corresponding phase portraits for
the case where the bifurcation diagram has no cusp points for $\alpha\neq 0$.
As can be seen from this figure, both points $\sigma_{0}$ and $\sigma_{\pi}$
do not lie on the bifurcation curve $\sigma_{\vartheta}$ and are isolated.

Figure 8: Bifurcation diagram (b) and examples of phase portraits (a), (c)
corresponding to different values of the integral $\kappa$ for $\alpha=0.5$,
$\beta=0.5$.

For this case, two types of qualitatively different phase portraits are
possible: for $\kappa=0$ and for $|\kappa|>0$. In Fig. 8, the points $s_{0}$
and $s_{1}$ denote fixed points of the family $\sigma_{\vartheta}$ which lie
on the chosen level sets of the first integral $\kappa=0$ and $\kappa=1$,
respectively.

As can be seen from Fig. 8a, for the case $\kappa=0$ the phase trajectories
pass through the poles $\vartheta=0$ and $\vartheta=\pi$. For a correct
passage through the pole, it is necessary to glue together the trajectories in
the phase portrait: $\dot{\vartheta}\rightarrow-\dot{\vartheta}$. Thus, in
this case each level set of the first integrals corresponds to one
connectedness component, which is in complete agreement with the bifurcation
diagram.

####  4.4.5 Bifurcation diagram (b)

As the parameter $\beta$ increases, the curve $\sigma_{s\pi}$ approaches the
point $\sigma_{\pi}$. When $\beta^{2}=1-\alpha$, the point $\sigma_{\pi}$
descends to the bifurcation curve $\sigma_{s\pi}$. Also, $\sigma_{\pi}$ ceases
to be isolated, and the branch $\sigma_{s\pi}$ loses smoothness at this point
($\kappa=0$). The corresponding type of diagram and the phase portraits are
presented in Fig. 9.

######  Remark.

A similar bifurcation, not when the parameter, but the level set of the
integrals is changed, can be observed in the classical problem of the Lagrange
top. The corresponding conditions for stability of a sleeping top are called
the Maievsky condition [49, 55].

Figure 9: Bifurcation diagram (b) and examples of phase portraits (a), (c)
corresponding to different values of the integral $\kappa$ for $\alpha=0.5$,
$\beta=1.1$.

For this case also, only two types of qualitatively different phase portraits
similar to diagram (a) are possible. The only difference is that the point
$\sigma_{\pi}$ is stable for $\kappa=0$.

####  4.4.6 Bifurcation diagram (c)

Further increase in the parameter $\beta$, with $\beta^{2}=1+\alpha$, leads to
a new bifurcation. The second isolated point $\sigma_{0}$ gives rise to
bifurcation branches $\sigma_{s0}$ and $\sigma_{u}$. The corresponding diagram
and the phase portraits are presented in Fig. 10.

Figure 10: Bifurcation diagram (b) and examples of phase portraits (a),
(c)–(e) corresponding to different values of the integral $\kappa$ for
$\alpha=0.5$, $\beta=3$.

As can be seen from the figure, in this case there exist four different types
of phase portraits. When $\kappa=0$ and $0<|\kappa|<\kappa_{c}$ (Figs. 10a,
10d) the phase portrait has three fixed points: two of them are stable
($\sigma_{0}$, $\sigma_{\pi}$ for $\kappa=0$ and $s_{11}$, $s_{12}$ for
$\kappa=\kappa_{1}$) and one is unstable ($u_{0}$ and $u_{1}$, respectively).
When $\kappa=\kappa_{c}$ (Fig. 10e) the curves $\sigma_{u}$ and $\sigma_{s0}$
merge at the point $u_{c}$. In this case, the phase portrait exhibits a
saddle-node bifurcation — the stable point and the unstable point merge and
disappear. When $|\kappa|>\kappa_{c}$ (Fig. 10c), only one fixed point remains
in the phase portrait.

####  4.4.7 Bifurcation diagram (d)

In the case $\alpha=0$, $\beta<1$ the bifurcation diagram consists of the
curve $\sigma_{\pi/2}$, which corresponds to rolling along the ellipsoid’s
equator, and of two isolated fixed points $\sigma_{0}$ and $\sigma_{\pi}$,
which lie on the same level set of the integrals.

Figure 11: Bifurcation diagram (b) and examples of phase portraits (a), (c)
corresponding to different values of the integral $\kappa$ for $\alpha=0$,
$\beta=0.5$.

The phase portrait has, for any value of the integral $\kappa$, a fixed point
for $\vartheta=\pi/2$, and is symmetric about this point. For $\beta<1$, two
qualitatively different types of phase portrait are possible (see Fig. 11).
When $\kappa=0$, the phase plane has three fixed points: two unstable ones,
$\sigma_{0}$ and $\sigma_{\pi}$, and a stable fixed point, $s_{0}$ (which
belongs to the family $\sigma_{\pi/2}$). When $\kappa>0$, only one stable
fixed point, $s_{1}$, which belongs to the family $\sigma_{\pi/2}$, remains in
the phase portrait.

####  4.4.8 Bifurcation diagram (e)

When $\beta=1$ (with constant $\alpha=0$), the fixed points $\sigma_{0}$ and
$\sigma_{\pi}$ corresponding to the vertical positions of the ellipsoid
descend to the bifurcation curve $\sigma_{\pi/2}$ and give rise to bifurcation
branches $\sigma_{s0}$ and $\sigma_{s\pi}$, respectively. In this case, they
become stable and cease to be isolated (see Fig. 12). unstable for $\kappa=0$.

Figure 12: Bifurcation diagram (b) and examples of phase portraits (a), (c),
(d) corresponding to different values of the integral $\kappa$ for $\alpha=0$,
$\beta=1.5$.

When $\kappa=0$, the phase portrait has three fixed points: two stable ones,
$\sigma_{0}$ and $\sigma_{\pi}$, and an unstable one, $u_{0}$, at
$\vartheta=\pi/2$.

When $0<|\kappa|<\kappa_{c}$, the phase portrait has three fixed points: two
stable ones and an unstable one. Next, when $|\kappa|=\kappa_{c}$, a pitch-
fork bifurcation occurs, such that the three fixed points merge into a stable
fixed point222We do not present here the phase portrait at the moment of
bifurcation for $\kappa=\kappa_{c}$ since visually it does not differ from the
case $\kappa>\kappa_{c}$.. As a result, when $|\kappa|>\kappa_{c}$, only one
stable point ($\vartheta=\pi/2$) remains in the phase portrait.

##  5 Analysis of the trajectories of an ellipsoid

In this section, we analyze how an ellipsoid behaves during motion on a plane,
depending on the initial conditions and parameters of the problem, that is, to
what motions of the ellipsoid the periodic trajectories of the reduced system
(10) correspond.

We have already ascertained that permanent rotations correspond to the rolling
motion of the ellipsoid (both its center of mass and the point of contact) in
a circle. In the limiting cases when $\vartheta=0$ or $\vartheta=\pi$, the
circle degenerates to a point, and as $\vartheta\rightarrow\pi/2$, it
degenerates to a straight line. Let us find characteristic trajectories of the
ellipsoid for other solutions.

In Section 2.4 we have obtained the system with one degree of freedom (10),
which admits an energy integral. In the general position case, the solution of
this system is periodic motion with frequency $\omega_{\vartheta}$ along
closed trajectories. From the expressions for the quadrature of the angle
$\psi$ (17) it follows that the function $\dot{\psi}$ is also periodic with
the same frequency $\omega_{\vartheta}$. Following our previous work [56, 57],
we represent the right-hand side of the quadrature for angle $\psi$ (17) as a
Fourier series expansion

$\dot{\psi}=\omega_{\psi}+\sum\limits_{n\in\mathbb{Z}\backslash\\{0\\}}{\rm\Psi}_{n}{\rm
e}^{\mathrm{i}n\omega_{\vartheta}t},$ (27)

where ${\rm\Psi}_{n}$ are the coefficients of Fourier series expansion,
$\mathrm{i}$ is an imaginary unit, and $\omega_{\psi}$ is the zeroth term of
the expansion which is calculated by the formula

$\omega_{\psi}=\frac{1}{T_{\vartheta}}\int\limits_{0}^{T_{\vartheta}}\dot{\psi}(t)dt,\quad
T_{\vartheta}=\frac{2\pi}{\omega_{\vartheta}}.$

Integrating the expression (27) over time, we can represent the expression for
$\psi(t)$ as

$\psi(t)=\omega_{\psi}t+\psi_{\vartheta}(t),\quad\psi_{\vartheta}(t)=\sum\limits_{n\in\mathbb{Z}\backslash\\{0\\}}\frac{{\rm\Psi}_{n}}{\mathrm{i}n\omega_{\vartheta}}{\rm
e}^{\mathrm{i}n\omega_{\vartheta}t}.$ (28)

Here $\psi_{\vartheta}(t)$ is also a periodic function with frequency
$\omega_{\vartheta}$.

Next, we analyze the trajectory of the center of mass of the ellipsoid. We
introduce the complex variable $\zeta=x_{c}+\mathrm{i}y_{c}$, where $x_{c}$
and $y_{c}$ are the coordinates of the center of mass of the ellipsoid. Using
(1), we can write the equation for changing the variable $\zeta$ in the form

$\displaystyle\dot{\zeta}=\frac{\kappa}{J(\vartheta)\sin\vartheta}(Z(\vartheta)+\alpha\cos\vartheta){\rm
e}^{\mathrm{i}\psi(t)}.$

In view of (28) and the Fourier series expansion of $\dot{\zeta}$, the
expression for $\zeta$ has the form

$\zeta=\zeta(0)+\sum\limits_{n\in\mathbb{Z}}\frac{v_{n}}{\mathrm{i}(\omega_{\psi}+n\omega_{\vartheta})}{\rm
e}^{\mathrm{i}(\omega_{\psi}+n\omega_{\vartheta})t},$

where $v_{n}$ are the coefficients of the Fourier series expansion of
$\dot{\zeta}$. Thus, in the general case the trajectory of the center of mass
of the ellipsoid is a bounded curve except for the cases where the resonance
relation $\omega_{\psi}+n\omega_{\vartheta}=0$ is satisfied.

Next, we introduce the concept of rotation number $N$ as a ratio between
frequencies $\omega_{\psi}$ and $\omega_{\vartheta}$

$N=-\frac{\omega_{\psi}}{\omega_{\vartheta}}=-\frac{1}{2\pi}\int\limits_{0}^{T_{\vartheta}}\dot{\psi}(t)dt.$

The value of the rotation number determines the type of trajectory of the
ellipsoid in absolute space. To each trajectory of the system (10) there
corresponds its own rotation number depending on the values of frequencies
$\omega_{\vartheta}$ and $\omega_{\psi}$. The frequencies $\omega_{\psi}$ and
$\omega_{\vartheta}$, in their turn, are uniquely defined by the values of the
first integrals $\kappa$, $\varepsilon$ and by the connectedness component of
the integral manifold $\mathcal{M}_{\kappa,\varepsilon}$, on which the
trajectory lies. Thus, the rotation number $N$ can be defined as a function
given in the (many-sheeted) bifurcation diagram. Hence, the type of trajectory
of the ellipsoid in absolute space is uniquely (up to the chosen sheet)
defined by the values of the first integrals $\kappa$ and $\varepsilon$.

We now turn to the question of possible types of trajectories of the ellipsoid
depending on $\kappa$ and $\varepsilon$.

We note that it follows from (17) that, when $\kappa=0$, one has $N=0$ for any
parameter values and initial conditions. As a result, the integral manifolds
$\mathcal{M}_{0,\varepsilon}$ are completely filled by periodic solutions,
i.e., they are families of resonant tori. The motions on these tori can be of
two types. For small energies $\varepsilon<\varepsilon_{\rm min}$ they are
oscillatory motions near an equilibrium position. For large energies
$\varepsilon>\varepsilon_{\rm min}$ the ellipsoid rolls to infinity along the
meridian crossing both of its vertices. The value of energy $\varepsilon_{\rm
min}$ is determined from the condition

$\varepsilon_{\rm min}=\left\\{\\!\begin{aligned}
&\left.\varepsilon\right|_{\sigma_{0}}=1+\alpha&\text{for
}\,\beta^{2}<1+\alpha,\\\
&\left.\varepsilon\right|_{\sigma_{u},\kappa=0}=\sqrt{\frac{1+\alpha^{2}-\beta^{2}}{1-\beta^{2}}}\beta&\text{for
}\,\beta^{2}>1+\alpha.\end{aligned}\right.$

The critical invariant manifolds corresponding to the intersection of
$\kappa=0$ with bifurcation curves are filled by inclined equilibrium
positions (18) and, in the case of instability, by trajectories asymptotic to
them.

Analyzing the resulting dependences, we can formulate the following
proposition as a conclusion from the above reasoning (see also [58, 59]).

######  Proposition.

The trajectory of the center of mass of the ellipsoid of revolution on a plane
can be of one of the following types:

  1. 1.

a point in the case of equilibrium positions $\sigma_{0}$, $\sigma_{\pi}$ and
an inclined equilibrium position for $\kappa=0$,
$\varepsilon=U(\vartheta_{*})$.

  2. 2.

a circle in the case of permanent rotations (15) with a constant inclination
angle $\vartheta_{0}=\operatorname{const}$ $($see Fig. 3$)$.

  3. 3.

an unbounded straight line in the case of permanent rotations $\sigma_{\pi/2}$
for $\alpha=0$ or in the case of rolling along the meridian for $\kappa=0$ and
$\varepsilon>\varepsilon_{\rm min}$.

  4. 4.

a straight-line segment in the case of pendulum motions near an equilibrium
position along the meridian for $\kappa=0$ and $\varepsilon<\varepsilon_{\rm
min}$ and in the case of trajectories asymptotic to unstable equilibrium
positions.

  5. 5.

a trajectory going to infinity without bound for integer rotation numbers
$N\in\mathbb{Z}$ $($see Fig. 13$)$.

  6. 6.

a closed periodic trajectory for fractionally rational rotation numbers
$N\in\mathbb{Q}\backslash\mathbb{Z}$ $($see Fig. 14a$)$.

  7. 7.

a bounded quasi-periodic curve for all other trajectories of the system (10),
except for singular trajectories asymptotic to unstable fixed points
$($separatrices$)$ $($see Fig. 14b$)$.

  8. 8.

a bounded curve asymptotically tending in forward and backward time to two
circles; these trajectories of the ellipsoid correspond to the trajectories of
the reduced system (10) which are asymptotic to unstable permanent rotations
$\sigma_{u}$.

  9. 9.

an unbounded curve asymptotically tending in forward and backward time to two
straight lines; these trajectories of the ellipsoid correspond to the
trajectories of the reduced system (10) which are asymptotic to unstable
equilibrium positions $\sigma_{0}$.

Since the rotation number $N$ depends on the values of the first integrals
$\kappa$ and $\varepsilon$, the conditions of resonance define some resonance
curves on the plane of first integrals $(\kappa,\varepsilon)$. Examples of
such curves for a symmetrically truncated sphere with a displaced center of
mass were constructed in [57].

Figure 13 shows an enlarged fragment of the bifurcation diagram for
$\alpha=0.5$ and $\beta=3$ (see Fig. 7c). The light blue curve corresponding
to the resonance $N=0$ lies near the curve $\sigma_{u}$ and “goes round the
curve” as it were. But since the lower part of the curve runs very close to
$\sigma_{u}$, it seems — when one looks at the figure — as if they coincide.
Higher-order resonances $N>0$ lie still closer to $\sigma_{u}$, and so we do
not represent them here. Also, Fig. 13 shows the trajectory of the center of
mass of the ellipsoid (the dark red solid curve in the left part of Fig. 13)
and of its point of contact (dotted dark blue curve) for the initial
conditions, which correspond to the resonance $N=0$. It can be seen that, for
these initial conditions, the ellipsoid does not, on average, drift along the
axis $Oy$, but goes to infinity along the axis $Ox$.

Figure 13: Right: an enlarged fragment of the bifurcation diagram (Fig. 7c,
$\alpha=0.5$, $\beta=3$) with the curve corresponding to resonance $N=0$.
Left: the trajectory of the center of mass (solid curve) and of the point of
contact (dotted curve) of the ellipsoid for the case where the resonance
relation $N=0$ is satisfied and for $\kappa=0.5$.

The numerical analysis has enabled us to find integer-valued resonances only
near the curve corresponding to unstable fixed points ($\sigma_{u}$), i.e.,
only for the case of diagram (c) in Fig. 7. Based on the results obtained, we
can formulate the following hypothesis.

######  Hypothesis.

The trajectories of an ellipsoid of revolution on a plane are bounded for

  1. $\bullet$

$\beta^{2}<1+\alpha$ except for the case $\kappa=0$ at
$\varepsilon>\varepsilon_{\rm min}$;

  2. $\bullet$

$\beta^{2}>1+\alpha$, $\alpha>0$ and $\kappa>\kappa_{\rm max}$, where
$\kappa_{\rm max}$ is a solution of the system

$N(\kappa,\varepsilon)=0,\quad\frac{\partial N}{\partial\varepsilon}=0.$

Figure 14 shows the trajectories of the center of mass and of the point of
contact of the ellipsoid at resonance $N=-1/7$ (Fig. 14a) and in the case of a
quasi-periodic trajectory (Fig. 14b). The values of the first integrals for
which these trajectories are constructed are shown in Fig. 13 as points (a)
and (b), respectively. Both of these cases correspond to bounded motion of the
ellipsoid on a plane.

Figure 14: (a) Example of a periodic trajectory (the case $N=-1/7$) for
$\kappa=0.8$ and $\vartheta(0)=0.4678$. (b) Example of a quasi-periodic
trajectory for $\kappa=1$ and $\vartheta(0)=2.1$. The solid dark red curve
denotes the trajectories of the center of mass, and the dotted dark blue curve
indicates the trajectories of the point of contact. Both parts of Fig. 14 are
constructed for $\alpha=0.5$, $\beta=3$, $\nu=0.5$, $\eta=0.5$ and
$\dot{\vartheta}(0)=0$.

## Summary

In this work, we have carried out a complete qualitative analysis of the
dynamics of an ellipsoid of revolution with a displaced center of mass on a
plane with no-slip and no-spin conditions. As is well known [4, 5], the
problem of a body of revolution rolling on a plane is integrable for the
above-mentioned constraints and, moreover, all integrals are represented in
elementary functions. For an ellipsoid of revolution, we have found particular
solutions corresponding to the rolling motion in a circle. Using these
solutions, we have constructed bifurcation diagrams and performed a complete
classification of these diagrams depending on the parameters of the problem.

In addition, we have shown that in the general case the motion of the
ellipsoid of revolution occurs along a quasi-periodic bounded curve. However,
depending on the values of the first integrals, there can also exist both
unbounded and closed periodic trajectories.

A possible avenue for future research is to simulate the motion of a mobile
robot of elliptic form by adding internal propulsion devices, for example,
gyrostats or a pendulum/a platform or their combination. In addition, to
simulate the motion of such a robot under conditions close to experimental
ones, it would be interesting to add to the model irregularity of the
supporting surface and/or forces and torques of rolling resistance.

## Acknowledgments

The authors thank Ivan Mamaev for useful comments and discussions.

## Funding

The work of A. Kilin (Sections 1–3) was carried out within the framework of
the state assignment of the Ministry of Science and Higher Education of Russia
(FEWS-2020-0009). The work of E. Pivovarova (Sections 4 and 5) is supported by
the Russian Science Foundation under grant 23-71-01045.

## References

  * [1] Hadamard, J., Sur les mouvements de roulement, _Mémoires de la Société des sciences physiques et naturelles de Bordeaux, 4 sér._ , 1895, vol. 5, pp. 397–417.
  * [2] Beghin, H., Sur les conditions d’application des équations de Lagrange à un système non holomorphe, _Bulletin de la Société Mathématique de France_ , 1929, vol. 57, pp. 118–124.
  * [3] Ehlers, K. and Koiller, J., Rubber rolling: Geometry and dynamics of 2-3-5 distributions, in _IUTAM Symposium on Hamiltonian Dynamics, Vortex Structures, Turbulence: Proceedings of the IUTAM Symposium held in Moscow, 25–30 August, 2006_ , Springer, pp. 469–479.
  * [4] Borisov, A. V. and Mamaev, I. S., Conservation laws, hierarchy of dynamics and explicit integration of nonholonomic systems, _Regular and Chaotic Dynamics_ , 2008, vol. 13, no. 5, pp. 443–490.
  * [5] Borisov, A. V., Mamaev, I. S., and Bizyaev, I. A., The hierarchy of dynamics of a rigid body rolling without slipping and spinning on a plane and a sphere, _Regular and Chaotic Dynamics_ , 2013, vol. 18, no. 3, pp. 277–328.
  * [6] Appell, P., Sur L’intégration des équations du mouvement dun corps pesant de révolution roulant par une arête circulaire sur un plan horizontal; cas particulier du cerceau, _Rendiconti del Circolo Matematico di Palermo (1884-1940)_ , 1900, vol. 14, no. 1, pp. 1–6.
  * [7] Appell, P., _Traité de mécanique rationnelle_ , vol. 2, Gauthier-Villars, 1904.
  * [8] Chaplygin, S. A., On a motion of a heavy body of revolution on a horizontal plane, _Regular and Chaotic Dynamics_ , 2002, vol. 7, no. 2, pp. 119–130.
  * [9] Kilin, A. A. and Ivanova, T. B., The Integrable Problem of the Rolling Motion of a Dynamically Symmetric Spherical Top with One Nonholonomic Constraint, _Russian Journal of Nonlinear Dynamics_ , 2023, vol. 19, no. 1, pp. 3–17.
  * [10] Kozlov, V. V. and Kolesnikov, N. N., On theorems of dynamics, _Journal of Applied Mathematics and Mechanics_ , 1978, vol. 42, no. 1, pp. 26–31.
  * [11] Kilin, A. A. and Ivanova, T. B., The Problem of the Rolling Motion of a Dynamically Symmetric Spherical Top with One Nonholonomic Constraint, _Russian Journal of Nonlinear Dynamics_ , 2023, vol. 19, no. 4, pp. 533–543.
  * [12] Borisov, A. V. and Ivanov, A. P., A Top on a Vibrating Base: New Integrable Problem of Nonholonomic Mechanics, _Regular and Chaotic Dynamics_ , 2022, vol. 27, no. 1, pp. 2–10.
  * [13] Kilin, A. A. and Pivovarova, E. N., Stability of Vertical Rotations of an Axisymmetric Ellipsoid on a Vibrating Plane, _Mathematics_ , 2023, vol. 11, no. 18, p. 3948.
  * [14] Borisov, A. and Mamaev, I., The Rolling Motion of a Rigid Body on a Plane and a Sphere. Hierarchy of Dynamics, _Regular and Chaotic Dynamics_ , 2002, vol. 7, no. 2, pp. 177–200.
  * [15] Routh, E. J., _The Advanced Part of a Treatise on the Dynamics of a System of Rigid Bodies_ , London: Macmillan and Co., 1892.
  * [16] Borisov, A. V., Mamaev, I. S., and Kilin, A. A., The rolling motion of a ball on a surface. New integrals and hierarchy of dynamics, _Regular and Chaotic Dynamics_ , 2002, vol. 7, no. 2, pp. 201–219.
  * [17] Dalla Via, M., Fassò, F., and Sansonetto, N., On the dynamics of a heavy symmetric ball that rolls without sliding on a uniformly rotating surface of revolution, _Journal of Nonlinear Science_ , 2022, vol. 32, no. 6, p. 84.
  * [18] Fassò, F. and Sansonetto, N., On Some Aspects of the Dynamics of a Ball in a Rotating Surface of Revolution and of the Kasamawashi Art, _Regular and Chaotic Dynamics_ , 2022, vol. 27, no. 4, pp. 409–423.
  * [19] Borisov, A. V., Ivanova, T. B., Kilin, A. A., and Mamaev, I. S., Nonholonomic rolling of a ball on the surface of a rotating cone, _Nonlinear Dynamics_ , 2019, vol. 97, no. 2, pp. 1635–1648.
  * [20] Ivanova, T. B., Non-holonomic rolling of a ball on the surface of a rotating cylinder, _ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik_ , 2020, vol. 100, no. 12, p. e202000067.
  * [21] Dragović, V., Gajić, B., and Jovanović, B., Spherical and Planar Ball BearingsNonholonomic Systems with Invariant Measures, _Regular and Chaotic Dynamics_ , 2022, vol. 27, no. 4, pp. 424–442.
  * [22] Dragović, V., Gajić, B., and Jovanović, B., Spherical and Planar Ball Bearingsa Study of Integrable Cases, _Regular and Chaotic Dynamics_ , 2023, vol. 28, no. 1, pp. 62–77.
  * [23] Borisov, A. V., Kilin, A. A., and Mamaev, I. S., Generalized Chaplygins transformation and explicit integration of a system with a spherical support, _Regular and Chaotic Dynamics_ , 2012, vol. 17, pp. 170–190.
  * [24] Fedorov, Y. N., Motion of a rigid body in a spherical suspension, _Vestnik Moskovskogo Universiteta. Seriya 1. Matematika. Mekhanika_ , 1988, , no. 5, pp. 91–93.
  * [25] Cendra, H. and Etchechoury, M., Rolling of a symmetric sphere on a horizontal plane without sliding or spinning, _Reports on Mathematical Physics_ , 2006, vol. 57, no. 3, pp. 367–374.
  * [26] Cendra, H., Etchechoury, M., and Ferraro, S. J., Impulsive control of a symmetric ball rolling without sliding or spinning, _Journal of Geometric Mechanics_ , 2010, vol. 2, no. 4, pp. 321–342.
  * [27] Koiller, J. and Ehlers, K., Rubber rolling over a sphere, _Regular and Chaotic Dynamics_ , 2007, vol. 12, pp. 127–152.
  * [28] Bizyaev, I. A., Borisov, A. V., and Mamaev, I. S., An Invariant Measure and the Probability of a Fall in the Problem of an Inhomogeneous Disk Rolling on a Plane, _Regular and Chaotic Dynamics_ , 2018, vol. 23, no. 6, pp. 665–684.
  * [29] Bolsinov, A. V., Borisov, A. V., and Mamaev, I. S., Rolling of a Ball without Spinning on a Plane: The Absence of an Invariant Measure in a System with a Complete Set of Integrals, _Regular and Chaotic Dynamics_ , 2012, vol. 17, pp. 571–579.
  * [30] Borisov, A. V., Mamaev, I. S., and Treschev, D. V., Rolling of a rigid body without slipping and spinning: Kinematics and dynamics, _J. Appl. Nonlinear Dyn_ , 2013, vol. 2, no. 2, pp. 161–173.
  * [31] Kazakov, A. O., Strange attractors and mixed dynamics in the problem of an unbalanced rubber ball rolling on a plane, _Regular and Chaotic Dynamics_ , 2013, vol. 18, no. 5, pp. 508–520.
  * [32] Borisov, A. V., Kazakov, A. O., and Pivovarova, E. N., Regular and chaotic dynamics in the rubber model of a Chaplygin top, _Regular and Chaotic Dynamics_ , 2016, vol. 21, pp. 885–901.
  * [33] Mamaev, I. S. and Vetchanin, E. V., Dynamics of rubber Chaplygin sphere under periodic control, _Regular and Chaotic Dynamics_ , 2020, vol. 25, pp. 215–236.
  * [34] Bizyaev, I. A. and Kazakov, A. O., Integrability and stochastic behavior in some nonholonomic dynamics problems, _Russian Journal of Nonlinear Dynamics_ , 2013, vol. 9, no. 2, pp. 257–265.
  * [35] Ylikorpi, T., Forsman, P., Halme, A. et al., Gyroscopic precession in motion modelling of ball-shaped robots, in _Proceedings of the 28th European Conference on Modelling and Simulation_ , ECMS.
  * [36] Jurdjevic, V., The geometry of the plate-ball problem, _Archive for Rational Mechanics and Analysis_ , 1993, vol. 124, pp. 305–328.
  * [37] Bicchi, A., Prattichizzo, D., and Sastry, S. S., Planning motions of rolling surfaces, in _Proceedings of 1995 34th IEEE conference on decision and control_ , vol. 3, IEEE, pp. 2812–2817.
  * [38] Zheng, M., Qiang, Z., Jinkun, L., and Yao, C., Control of a spherical robot: Path following based on nonholonomic kinematics and dynamics, _Chinese Journal of Aeronautics_ , 2011, vol. 24, no. 3, pp. 337–345.
  * [39] Cai, Y., Zhan, Q., and Xi, X., Path tracking control of a spherical mobile robot, _Mechanism and Machine Theory_ , 2012, vol. 51, pp. 58–73.
  * [40] Alves, J. and Dias, J., Design and control of a spherical mobile robot, _Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering_ , 2003, vol. 217, no. 6, pp. 457–467.
  * [41] Zhan, Q., Cai, Y., and Yan, C., Design, analysis and experiments of an omni-directional spherical robot, in _2011 IEEE International Conference on Robotics and Automation_ , IEEE, pp. 4921–4926.
  * [42] Chen, W.-H., Chen, C.-P., Yu, W.-S., Lin, C.-H., and Lin, P.-C., Design and implementation of an omnidirectional spherical robot Omnicron, in _2012 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM)_ , IEEE, pp. 719–724.
  * [43] Karavaev, Y. L., Spherical robots: An up-to-date overview of designs and features, _Russian Journal of Nonlinear Dynamics_ , 2022, vol. 18, no. 4, pp. 709–750.
  * [44] Artemova, E. M., Karavaev, Y. L., Mamaev, I. S., and Vetchanin, E. V., Dynamics of a spherical robot with variable moments of inertia and a displaced center of mass, _Regular and Chaotic Dynamics_ , 2020, vol. 25, pp. 689–706.
  * [45] Marigo, A. and Bicchi, A., Rolling bodies with regular surface: Controllability theory and applications, _IEEE Transactions on Automatic Control_ , 2000, vol. 45, no. 9, pp. 1586–1599.
  * [46] Mukherjee, R., Minor, M. A., and Pukrushpan, J. T., Motion planning for a spherical mobile robot: Revisiting the classical ball-plate problem, _J. Dyn. Sys., Meas., Control_ , 2002, vol. 124, no. 4, pp. 502–511.
  * [47] García-Agúndez Blanco, A., García Vallejo, D., and Freire Macías, E., Analytical and numerical stability analysis of a toroidal wheel with nonholonomic constraints, _Nonlinear Dynamics_ , 2024, vol. 112, no. 4, pp. 2453–2476.
  * [48] Poincaré, H., Sur le forme nouvelle des equations de la mecanique, _C. R. Acad. Sci. Paris_ , 1901, vol. 132, pp. 369–371.
  * [49] Borisov, A. and Mamaev, I. S., _Rigid body dynamics_ , _De Gruyter Studies in Mathematical Physics_ , vol. 52, Walter de Gruyter GmbH & Co KG, 2018\.
  * [50] Bolsinov, A. V., Borisov, A. V., and Mamaev, I. S., Topology and stability of integrable systems, _Russian Mathematical Surveys_ , 2010, vol. 65, no. 2, p. 259.
  * [51] Bizyaev, I. A., Bolsinov, A., Borisov, A., and Mamaev, I., Topology and bifurcations in nonholonomic mechanics, _International Journal of Bifurcation and Chaos_ , 2015, vol. 25, no. 10, p. 1530028.
  * [52] Borisov, A. V. and Mamaev, I. S., Topological Analysis of an Integrable System Related to the Rolling of a Ball on a Sphere, _Regular and Chaotic Dynamics_ , 2013, vol. 18, pp. 356–371.
  * [53] Karapetyan, A., Families of permanent rotations of triaxial ellipsoid on rough horizontal plane and their branchings, _Actual Problems of Classical and Celestial Mechanics_ , 1998, pp. 46–51.
  * [54] Bizyaev, I. A. and Mamaev, I. S., Permanent Rotations in Nonholonomic Mechanics. Omnirotational Ellipsoid, _Regular and Chaotic Dynamics_ , 2022, vol. 27, no. 6, pp. 587–612.
  * [55] Markeev, A. P., The dynamics of a body contacting a rigid surface, _M.: Nauka_ , 1992, p. 336.
  * [56] Kilin, A. A. and Pivovarova, E. N., The rolling motion of a truncated ball without slipping and spinning on a plane, _Regular and Chaotic Dynamics_ , 2017, vol. 22, no. 3, pp. 298–317.
  * [57] Kilin, A. A. and Pivovarova, E. N., Qualitative Analysis of the Nonholonomic Rolling of a Rubber Wheel with Sharp Edges, _Regular and Chaotic Dynamics_ , 2019, vol. 24, pp. 212–233.
  * [58] Borisov, A. V., Mamaev, I. S., and Kilin, A. A., Dynamics of rolling disk, _Regular and Chaotic Dynamics_ , 2003, vol. 8, no. 2, pp. 201–212.
  * [59] Borisov, A. V. and Mamaev, I. S., An inhomogeneous Chaplygin sleigh, _Regular and Chaotic Dynamics_ , 2017, vol. 22, pp. 435–447.