# Riemannian metrics on 2D-manifolds related to the Euler-Poinsot rigid body motion

Bernard Bonnard, Olivier Cots, Jean-Baptiste Pomet, Nataliya Shcherbakova

► To cite this version:

Bernard Bonnard, Olivier Cots, Jean-Baptiste Pomet, Nataliya Shcherbakova. Riemannian metrics on 2D-manifolds related to the Euler-Poinsot rigid body motion. ESAIM: Control, Optimisation and Calculus of Variations, 2014, 20 (3), pp.864-893. 10.1051/cocv/2013087. hal-00918587v2

HAL Id: hal-00918587

https://inria.hal.science/hal-00918587v2

Submitted on 14 Jan 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
RIEMANNIAN METRICS ON 2D-MANIFOLDS RELATED TO THE
EULER-POINSOT RIGID BODY MOTION *

BERNARD BONNARD<fnref n="1" />, OLIVIER COTS<fnref n="2" />, JEAN-BAPTISTE POMET<fnref n="2" /> AND NATALIYA SHCHERBAKOVA<fnref n="3" />

**Abstract.** The Euler-Poinsot rigid body motion is a standard mechanical system and it is a model for left-invariant Riemannian metrics on SO(3). In this article using the Serret-Andoyer variables we parameterize the solutions and compute the Jacobi fields in relation with the conjugate locus evaluation. Moreover, the metric can be restricted to a 2D-surface, and the conjugate points of this metric are evaluated using recent works on surfaces of revolution. Another related 2D-metric on S<fnref n="2" /> associated to the dynamics of spin particles with Ising coupling is analysed using both geometric techniques and numerical simulations.

**AMS Subject Classification.** 49K15, 53C20, 70Q05, 81Q93.

January 14, 2014.

# 1. INTRODUCTION

The Euler-Poinsot rigid body motion describing the inertial revolution of a rigid body is a standard mechanical system associated to the attitude control in space mechanics. The model is the following: let {$e_1, e_2, e_3$} be a fixed orthonormal frame in $\mathbb{R}^3$, and denote by {$E_1, E_2, E_3$} an orthonormal frame attached to the body. The attitude of the body is represented by the matrix of directional cosines $R(t) = (E_1(t), E_2(t), E_3(t))$, which transforms {$e_i$}$_{i=1}^3$ into {$E_i$}$_{i=1}^3$. If the moving frame {$E_1, E_2, E_3$} coincides with the principal axes of inertia of the body, its trajectories are solutions to the following optimal control system on SO(3):

$$ \frac{dR(t)}{dt} = R(t) \begin{pmatrix} 0 & -u_3 & u_2 \\ u_3 & 0 & -u_1 \\ -u_2 & u_1 & 0 \end{pmatrix} = \sum_{i=1}^{3} u_i \vec{A}_i(R(t)), \\ \frac{1}{2} \int_{0}^{T} \sum_{i=1}^{3} I_i u_i^2(t) dt \rightarrow \min_{u(\cdot)} $$

<footnote>Keywords and phrases: Euler-Poinsot rigid body motion, conjugate locus on surfaces of revolution, Serret-Andoyer metric, spins dynamics</footnote>

<footnote>* Work supported by ANR Geometric Control Methods (project No. NT09-504490).</footnote>

<footnote n="1">¹ Institut de Mathématiques de Bourgogne, 9 avenue Savary, 21078 Dijon, France; e-mail: bernard.bonnard@u-bourgogne.fr</footnote>

<footnote n="2">² INRIA Sophia-Antipolis Méditerranée, 2004, route des Lucioles, 06902 Sophia Antipolis, France; e-mail: olivier.cots@inria.fr</footnote>

<footnote n="3">³ Université de Toulouse, INPT, UPS, Laboratoire de Génie Chimique, 4 allée Emile Monso, 31432 Toulouse, France; e-mail: nataliya.shcherbakova@inp-toulouse.fr</footnote>
where $I_i$'s are the principal momenta of inertia of the body. Geometrically it amounts to define a left-invariant metric on $SO(3)$, where the basis $\vec{A}_i(R) = RA_i$ forms an orthonormal frame, and conversely, every invariant metric on $SO(3)$ can be set in this normal form, the ratio of the principal momenta of inertia $I_1/I_3$ and $I_2/I_3$ being the invariants of the metric.

Numerous articles are devoted to the analysis of the extremal curves of the described problem: geometric properties of the extremals, integrability properties using the Poinsot representation, and explicit computations of the solution by means of the Jacobi or Weierstrass elliptic functions [2, 17, 20].

According to the Pontryagin maximum principle ([23]), the optimal solutions of the problem are projections of the extremal curves, i.e., the integral curves of the Hamiltonian vector field $\vec{H}$ associated to the function $H = \frac{1}{2}(M_1^2 I_1^{-1} + M_2^2 I_2^{-1} + M_3^2 I_3^{-1})$, where the vector $M = (M_1, M_2, M_3)$, which represents the angular momentum of the body measured in the moving frame, is related to the angular velocity vector $\Omega = (u_1, u_2, u_3)$: $M_i = I_i u_i$. In this representation the limit case $I_2 \to +\infty$ leads to the sub-Riemannian problem associated to $H = \frac{1}{2}(M_1^2 I_1^{-1} + M_2^2 I_2^{-1})$. Conversely, every left invariant sub-Riemannian metric on $SO(3)$ can be set in this normal form where $k^2 = I_1/I_3$ is the only invariant of the problem [26].

Beyond the standard integrability analysis of the extremals, an important geometric control problem is to compute the conjugate and cut loci of such a metric. Recall that given a complete Riemannian manifold $(M, g)$, the conjugate locus $C(q_0)$ of a point $q_0 \in M$ is the set of points, where the geodesic curves emanating from $q_0$ loose their optimality for the $C^1$ topology on the set of curves, while the cut locus $C_{\text{cut}}(q_0)$ is the set of points where they loose optimality globally. Computing these sets is equivalent to solving the Hamilton-Jacobi-Bellman equation. In general, it is a very complicated problem and little is known in the literature about the explicit construction of conjugate and cut loci ([22]). Very recently a detailed analysis was done for geodesic flows of Riemannian metrics on 2-spheres of revolution, which are reflectionally symmetric with respect to the equator ([6,25]) with an application to the case of an ellipsoid of revolution. It was extended to a general ellipsoid [16], such computations were also done in the geometric control context [8,11] and in [3] for the axisymmetric ellipsoid of inertia.

This article is a first step to perform such a computation for a left-invariant Riemannian or sub-Riemannian metric on $SO(3)$ by considering the reduction on two-dimensional manifolds. The first case to analyze is the Riemannian metric on a two-dimensional surface related to the Serret-Andoyer variables [15], which define a two-dimensional projection of the full problem. Indeed, assuming $I_1 > I_2 > I_3$, on can introduce suitable symplectic coordinates $(x, y, w, p_x, p_y, p_w)$, such that the Hamiltonian takes the form $H = \frac{1}{4}[z(y)p_x^2 + (2I_3^{-1} - z(y))p_y^2]$, where $z = 2(I_1^{-1} \sin^2 y + I_2^{-1} \cos^2 y)$. In such a representation the Hamiltonian $H$, which does not depend explicitly on $w$ and $p_w$, describes a Riemannian metric in the $(x, y)$-variables. In particular, this reduction is crucial for the integration of the system and for the analysis of the Jacobi equation associated to the left-invariant metric on $SO(3)$.

The second case relies on the sub-Riemannian situation, which can be realized by making $I_2 \to +\infty$. It defines a reduced metric $g = (dr_1^2 + k^2 dr_3^2)r_2^{-2}$ on the two sphere $S^2$, where $k^2 = I_1 I_3^{-1}$, and $r_i$'s are the coordinates on the sphere such that $r_3^2 = 1 - r_1^2 - r_2^2$. Such a metric has a singularity at $r_2 = 0$, and it is related to the dynamics of a three spin system with Ising coupling [28]. Introducing the spherical coordinates $r_2 = \cos \varphi$, $r_1 = \sin \varphi \cos \theta$, $r_3 = \sin \varphi \sin \theta$, one gets

$$g = (\cos^2 \theta + k^2 \sin^2 \theta)d\varphi^2 + 2(k^2 - 1) \tan \varphi \sin \theta \cos \theta d\varphi d\theta + \tan^2 \varphi(\sin^2 \theta + k^2 \cos^2 \theta)d\theta^2,$$

which in the case $k=1$ corresponds to the standard Grushin metric on $S^2$ introduced in [10], whose conjugate and cut loci were described in [6], [7].

In the two examples above: the Serret-Andoyer problem and in the limit case $k=1$ of the spin problem, the reduction leads to a metric of the form $g = d\varphi^2 + \mu(\varphi)d\theta^2$ on a surface of revolution. As it was shown in the previous articles in case of a sphere, under some suitable assumptions the description of conjugate and cut loci can be made explicitly in terms of the period mapping $p_\theta \to T(p_\theta)$ (see Section 3.3.1 for the formal definition) of the $\varphi$ variable.
Besides the geometric analysis, an important role in our study is played by the Hampath code [12], which allows the explicit numerical computation based on the evaluation of conjugate points and determination of the cut locus using the continuation methods. This will be illustrated by the two examples described above.

The organization of this article is the following. In Section 2 we present the preliminary results to conduct our analysis, based on [17], and introduce the Serret-Andoyer formalism in Euler-Poinsot problem [15]. Section 3 contains the main contribution of the article: the study of the two-dimensional Riemannian metric $g = \frac{2dx^2}{z(y)} + \frac{2dy^2}{2C-z(y)}$, which is performed in two steps. First, applying the methodology proposed in [8], we parameterize the extremals and estimate the conjugate times. Then we compute the Darboux normal form $d\varphi^2 + \mu(\varphi)d\theta^2$ of the Serret-Andoyer metric in order to derive the global geometric properties of the conjugate locus in relation with the Gauss curvature. Our analysis extends the previous study of metrics on two dimensional spheres of revolution. In Section 4 we introduce the optimal control problem of the dynamics of spin systems [27, 28] in connection with the left-invariant sub-Riemannian metric on SO(3) and present the preliminary geometric results. In Section 5 we illustrate our results by numerical computations using the Hampath code.

## 2. LEFT-INVARIANT METRICS ON SO(3) AND THE SERRET-ANDOYER FORMALISM IN THE EULER-POINSOT RIGID BODY MOTION

### 2.1. Left-invariant metrics on SO(3)

For the presentation and notations concerning the left-invariant metric on SO(3) and its relation with geo-metric optimal control we use notations of V. Jurdjevic [17]. Denote {$e_1, e_2, e_3$} the (inertial) Euclidean frame in $\mathbb{R}^3$. The rigid body position can be represented by the matrix $R(t) = (E_1(t), E_2(t), E_3(t))$ describing the transformation from $e_i$'s to $E_i$'s, where $E_1, E_2, E_3$ form an orthonormal frame attached to the body. In addition, denote by $A_i$ the standard basis of SO(3)

$$A_1 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{pmatrix}, \quad A_3 = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix},$$

and by $[\cdot, \cdot]$ the Lie bracket computed with the convention $[A, B] = AB - BA$. The vector fields $\vec{A}_i(R) = RA_i$ form a basis of left-invariant vector fields on SO(3) and one has

$$[A_1, A_2] = -A_3, \quad [A_1, A_3] = A_2, \quad [A_2, A_3] = -A_1.$$

Consider now the following optimal control problem on SO(3):

$$\frac{dR}{dt} = \sum_{i=1}^{3} u_i \vec{A}_i(R), \quad \frac{1}{2} \int_{0}^{T} \sum_{i=1}^{3} I_i u_i^2 dt \rightarrow \min_{u(.)}, \quad (1)$$

where the $I_i$'s are the principal momenta of inertia of the body. The Euler-Poinsot dynamics can be derived using the Pontryagin maximum principle [23] and the appropriate (Poincaré) coordinates in the following way. Let $\lambda$ be an element of $T_R^*SO(3)$ and denote $H_i = \lambda(\vec{A}_i(R)), i = 1, 2, 3$, the symplectic lifts of the vector fields $\vec{A}_i$. The pseudo-Hamiltonian associated to the problem (1) takes the form

$$H = \sum_{i=1}^{3} u_i H_i - \frac{1}{2} \sum_{i=1}^{3} I_i u_i^2.$$
The maximization condition of the Pontryagin maximum principle implies that along the extremal solutions $\frac{\partial H}{\partial u} = 0$, and hence $u_i = H_i I_i^{-1}$, $i = 1, 2, 3$. Plugging such $u_i$'s into $H$ we get the true Hamiltonian

$$H_n = \frac{1}{2} \left( \frac{H_1^2}{I_1} + \frac{H_2^2}{I_2} + \frac{H_3^2}{I_3} \right).$$

From the point of view of Mechanics, the vector $H = (H_1, H_2, H_3)$ describes the angular momentum of the body measured in the moving frame and related to the angular velocity in the following way: $H_i = I_i u_i$. Then the well-known Euler equation describing the rigid body evolution takes the form

$$\frac{dH_i}{dt} = dH_i(\vec{H}_n) = \{H_i, H_n\},$$

where $\{\cdot, \cdot\}$ denotes the Poisson bracket. Using the standard relation between Poisson and Lie brackets: $\{H_i, H_j\} = \lambda([A_i, A_j])$, we immediately get the following differential equations

$$\begin{aligned}
\frac{dH_1}{dt} &= \frac{1}{I_2}\{H_1, H_2\}H_2 + \frac{1}{I_3}\{H_1, H_3\}H_3 = H_2H_3\left(\frac{1}{I_3} - \frac{1}{I_2}\right), \\
\frac{dH_2}{dt} &= H_1H_3\left(\frac{1}{I_1} - \frac{1}{I_3}\right), \quad \frac{dH_3}{dt} = H_1H_2\left(\frac{1}{I_2} - \frac{1}{I_1}\right).
\end{aligned}
\qquad (2)$$

**Proposition 2.1.** Euler equations (2) are integrable by quadratures using the two first integrals: the Hamiltonian $H_n$ and the Casimir $|H|^2 = H_1^2 + H_2^2 + H_3^2$. The solutions to Euler's equations can be seen as the curves on the energy ellipsoid formed by its intersection with the sphere of the constant angular momentum.

Such curves are called polhodes in Classical Mechanics.

**Remark 2.2.** As it was pointed in [17], the above property holds for every left-invariant Hamiltonian of the form $H_n = Q(H_1, H_2, H_3)$ on $SO(3)$. In particular, it is true in the SR-case, which is established next as a limit case of the Riemannian one.

## 2.2. The sub-Riemannian case

Setting $u_2 = \epsilon v_2$, with $\epsilon \to 0$, one gets a control system with two control inputs only. Since $u_i = H_i/I_i$, this is equivalent to the assumption $I_2 \to +\infty$. Then

$$H_n = \frac{1}{2} \left( \frac{H_1^2}{I_1} + \frac{H_3^2}{I_3} \right)$$

and the Euler equations take the form

$$\frac{dH_1}{dt} = \frac{H_2H_3}{I_3}, \quad \frac{dH_2}{dt} = H_1H_3\left(\frac{1}{I_1} - \frac{1}{I_3}\right), \quad \frac{dH_3}{dt} = -\frac{H_1H_2}{I_1}. \qquad (3)$$

The parameter $k^2 = I_1/I_3$ corresponds to the invariant classifying the SR-metrics on $SO(3)$ described in [27].

## 2.3. The Euler angles

Once the Euler equation is integrated, the next step is to find the solution of the full system. It relies on the following general property [17]:

**Proposition 2.3.** For each invariant Hamiltonian $H_n$ on $SO(3)$ the full system is integrable by quadratures using the four first integrals: the Hamiltonian $H_n$ and the Hamiltonian lifts of the right invariant vector fields $A_i R$.
Note that the remaining integrals are simply deduced from the Noether theorem in the Hamiltonian form. An explicit description of the quadratures can be easily obtained using the chart on SO(3) associated to the Euler angles $\Phi_1, \Phi_2, \Phi_3$, which can be defined by the following decomposition of the rotation matrix:

$$R = \exp(\Phi_1 A_3) \circ \exp(\Phi_2 A_2) \circ \exp(\Phi_3 A_3).$$

As it is shown in [17], the angles $\Phi_2$ and $\Phi_3$ can be found from the relations

$$H_1 = -|H| \sin \Phi_2 \cos \Phi_3, \quad H_2 = |H| \sin \Phi_2 \sin \Phi_3, \quad H_3 = |H| \cos \Phi_2, \quad (4)$$

while $\Phi_1$ is computed by integrating the differential equation

$$\frac{d\Phi_1}{dt} = \frac{|H|}{H_2 \sin \Phi_3 - H_1 \cos \Phi_3} \left( \sin \Phi_3 \frac{\partial H_n}{\partial H_2} - \cos \Phi_3 \frac{\partial H_n}{\partial H_1} \right). \quad (5)$$

In terms of the Euler angles, the Hamiltonian takes the form

$$H_n = \frac{1}{2I_1} \left( p_2 \sin \Phi_3 - \frac{\cos \Phi_3}{\sin \Phi_2} (p_1 - p_3 \cos \Phi_2) \right)^2 + \frac{1}{2I_2} \left( p_2 \cos \Phi_3 + \frac{\sin \Phi_3}{\sin \Phi_2} (p_1 - p_3 \cos \Phi_2) \right)^2 + \frac{1}{2I_3} p_3^2,$$

where $p_i$ denotes the canonical impulse associated to $\Phi_i$, $i = 1, 2, 3$. The corresponding Riemannian metric reads

$$g = \left( (I_1 \cos^2 \Phi_3 + I_2 \sin^2 \Phi_3) \sin \Phi_2^2 + I_3 \cos^2 \Phi_2 \right) d\Phi_1^2 + \left( I_2 \cos^2 \Phi_3 + I_1 \sin^2 \Phi_3 \right) d\Phi_2^2 + I_3 d\Phi_3^2 + \frac{(I_2 - I_1)}{2} \sin 2\Phi_3 \sin \Phi_2 d\Phi_1 d\Phi_2 + I_3 \cos \Phi_2 d\Phi_1 d\Phi_3.$$

## 2.4. The Serret-Andoyer variables

The complete integration of the Euler-Poinsot problem can be performed by choosing symplectic coordinates adapted to the problem by applying the Serret-Andoyer transformation. Indeed, it is well known from the Euler-Poinsot analysis that, except the solutions of the Euler equation corresponding to separating polhodes, every trajectory evolves on a two dimensional torus, and the motion is defined by two angular variables. This leads to the Liouville action-angle representation in a degenerated form, where one frequency is zero.

A complete description of the Serret-Andoyer transformation and its relation to the Euler angles is given, for instance, in [15]. Symplectic coordinates $(g, k, l, G, K, L)$, where $G, K, L$ denote the canonical impulses associated to the variables $g, k, l$, are defined by

$$H_1 = \sqrt{G^2 - L^2} \sin l, \quad H_2 = \sqrt{G^2 - L^2} \cos l, \quad H_3 = L.$$

So, $G = |H|$ and the Hamiltonian $H_n$ takes the form

$$H_n(g, k, l, G, K, L) = \frac{1}{2} \left( \frac{\sin^2 l}{I_1} + \frac{\cos^2 l}{I_2} \right) (G^2 - L^2) + \frac{L^2}{2I_3}. \quad (6)$$

This yields the Hamiltonian given in the Introduction:

$$H_a = \frac{1}{2} \left( (p_x^2 - p_y^2)(A \sin^2 y + B \cos^2 y) + C p_y^2 \right),$$

where $y = l = \arctan(\frac{H_1}{H_2})$, $p_x = G$, $p_y = L$, $A = I_1^{-1}$, $B = I_2^{-1}$, $C = I_3^{-1}$.
Note that by construction $G \ge |L|$. The inverse transformation towards the Euler angles can be found as follows:

$$ 
\begin{aligned}
p_1 &= K, \\
p_3 &= L, \\
p_2 \sin \Phi_2 &= -\frac{1}{G} \sqrt{G^2 - L^2} \sqrt{G^2 - K^2} \sin g, \\
\cos \Phi_2 &= \frac{1}{G^2} \left( KL - \sqrt{G^2 - L^2} \sqrt{G^2 - K^2} \cos g \right),
\end{aligned}
$$

which finally yields

$$ \Phi_1 = k + \arctan \left( \frac{L - K \cos \Phi_2}{p_2 \sin \Phi_2} \right), \quad \Phi_3 = l + \arctan \left( \frac{K - L \cos \Phi_2}{p_2 \sin \Phi_2} \right). $$

An additional canonical transformation, based on the Hamilton-Jacobi method, described in a recent work [19], leads to a standard action-angle representation of the Euler-Poinsot problem. Instead, in this paper we will use the Serret-Andoyer formalism to transform the Euler-Poinsot problem into a Riemannian problem on a two-dimensional surface. Denoting $z(y) = 2(A \sin^2 y + B \cos^2 y)$, we rewrite the Serret-Andoyer Hamiltonian as follows:

$$ H_a = \frac{1}{2} \left( \frac{z(y)}{2} p_x^2 + \left( C - \frac{z(y)}{2} \right) p_y^2 \right). \qquad (7) $$

Observe that $H_a$ is positive. Moreover, it defines a Riemannian metric

$$ g_a = \frac{2}{z(y)} dx^2 + \frac{2}{2C - z(y)} dy^2, \qquad (8) $$

provided the momenta of inertia of the body are ordered as $A < B < C$. Indeed, in this case $z(y) \in [2A, 2B]$ and $2C - z(y) > 0$. In what follows we will call $g_a$ the Serret-Andoyer metric.

### 3. THE SERRET-ANDOYER RIEMANNIAN METRIC

#### 3.1. The pendulum representation

We start by recalling some well known facts concerning the Euler-Poinsot rigid body dynamics. From now on we assume<fnref n="1" /> $I_1 > I_2 > I_3$, or equivalently, $A < B < C$. According to Proposition 2.1, the physical motion occurs if $2h \in [Ap_x^2, Cp_x^2]$. Every polhode is periodic, except the pair of separating polhodes formed by two circles contained in the planes of intersection of the energy ellipsoid $AH_1^2 + BH_2^2 + CH_3^2 = 2h$ with the sphere of angular momentum $H_1^2 + H_2^2 + H_3^2 = p_x^2 = 2h/B$. These two planes are defined by the equation $H_3 = \pm \sqrt{\frac{B-A}{C-B}}H_1$. If $Ap_x^2 = 2h$ or $Cp_x^2 = 2h$, then the corresponding periodic polhodes degenerate into the pairs of points on the energy ellipsoid, which describe the steady rotations of the body around the axes of its maximal and minimal momenta of inertia.

In the Serret-Andoyer representation, the Hamiltonian

$$ H_a(x, y, w, p_x, p_y, p_w) = \frac{1}{2} \left[ (p_x^2 - p_y^2) (A \sin^2 y + B \cos^2 y) + Cp_y^2 \right] $$

does not depend on $w$ and $p_w$. Hence

$$ \frac{dw}{dt} = \frac{dp_w}{dw} = 0, $$

<footnote n="1">¹This assumption is not restrictive due to the spherical symmetry of the Hamiltonian $H_n = \frac{1}{2} (AH_1^2 + BH_2^2 + CH_3^2)$.</footnote>
and the reduced dynamics on the $(x, y)$ plane is described by the system of equations

$$ 
\begin{aligned}
\frac{dx}{dt} &= p_x(A \sin^2 y + B \cos^2 y), & \frac{dp_x}{dt} &= 0, \\
\frac{dy}{dt} &= p_y(C - A \sin^2 y - B \cos^2 y), & \frac{dp_y}{dt} &= (B - A)(p_x^2 - p_y^2) \sin y \cos y.
\end{aligned}
\qquad (9) $$

The canonical impulse $p_x$ is a first integral of this system corresponding to the cyclic variable $x$, while on the $(y, p_y)$ plane we obtain a 2D pendulum-type phase portrait. Indeed, the Hamiltonian function $H_a$ is $\pi$-periodic with respect to $y$-variable, it verifies the following symmetry relations:

$$ H_a(y, p_y) = H_a(y, -p_y), \qquad H_a(y, p_y) = H_a(-y, p_y), \qquad (10) $$

and $p_y = 0$, $y = k\pi/2$, $k = 0, 1$ are the equilibrium points. The standard computation shows that in the neighborhood of the point $y = 0$, $p_y = 0$ the eigenvalues of the linearized system are solutions to the equation $\lambda^2 = (B - A)(C - B)p_x^2$, and hence both eigenvalues are real since $(B - A)(C - B) > 0$. Similarly, in the neighborhood of the point $y = \pi/2$, $p_y = 0$ one gets the equation $\lambda^2 = -(B - A)(C - A)p_x^2$, whose roots are purely imaginary.

To obtain a more detailed picture, consider the fixed energy level set $H_a = h$ with $h \in [Ap_x^2/2, Cp_x^2/2]$. Since $\dot{y} = p_y(C - z/2)$, the evolution of the variable $y$ satisfies the natural mechanical system of the form

$$ \dot{y}^2 + V(y) = 0, \qquad V(y) = \left(C - \frac{z(y)}{2}\right) \left(\frac{p_x^2 z(y)}{2} - 2h\right). \qquad (11) $$

In other words, for all initial conditions, $y(\cdot)$ is a trajectory of the natural mechanical system $\dot{y}^2 + V(y) = \mathcal{H}$ lying on the zero energy level $\mathcal{H} = 0$. Observe that $V(y) = 0$ implies $\cos^2 y = \frac{2h - Ap_x^2}{p_x^2(B-A)} = \xi_1$, where $\xi_1 \in [0, 1]$ if $h \in [Ap_x^2/2, Cp_x^2/2]$, and $\xi_1 \notin [0, 1]$ otherwise. Combining these facts with the analysis of the equilibria we obtain the following classification of the trajectories of (11) for $y \in [0, \pi]$:

- $h \in [Ap_x^2/2, Cp_x^2/2]$. System (11) admits three types of trajectories: a stable equilibrium $y = \pi/2$ if $h = Ap_x^2/2$, an unstable equilibrium and a separatrice $y = 0$ mod $\pi$ if $h = Cp_x^2/2$, which describes the motion along the separating polhode, and oscillating periodic trajectory with $y(t) \in [\arccos \sqrt{\xi_1}, \pi - \arccos \sqrt{\xi_1}]$ describing the polhodes around the major semi-axis of the energy ellipsoid.

- $h \in [Cp_x^2/2, Ap_x^2/2]$. The system (11) admits rotational trajectories for $h \in [Cp_x^2/2, Ap_x^2/2]$, the limit case $h = Cp_x^2/2$ describing the steady rotation around the minor axis of the energy ellipsoid.

We summarize our analysis in the following

**Proposition 3.1.** The pendulum motion in the $(y, p_y)$ plane can be interpreted on the cylinder $y \in [0, \pi]$ mod $\pi$ with a stable equilibrium at $y = \pi/2$ and an unstable one at $y = 0$ mod $\pi$. There are two types of periodic trajectories on the cylinder: oscillating trajectories homotopic to zero, and rotating trajectories, while non-periodic trajectories are separatrices joining 0 to $\pi$, which correspond to separating polhodes. Moreover, all these trajectories share the reflectional symmetry with respect to the axes $y = 0$ and $p_y = 0$.

In particular, it follows that in order to parameterize all phase trajectories on the $(y, p_y)$ plane it would be sufficient to integrate (9) with the initial condition $y(0) = \pi/2$.

On what follows we will extend the terms oscillating (or rotating) also to the trajectories on the $(x, y)$ plane according to the behaviour of the $y$ component.

## 3.2. Parameterization of the extremal curves and the conjugate locus

The parameterization of the extremal curves of the Euler-Poinsot motion using elliptic functions is standard (see, for instance, in [20]). The aim of this section is to use a different method introduced in [8] in connection with the computation of the conjugate locus.
### 3.2.1. Parameterization of the periodic trajectories

Consider first the trajectories corresponding to the ascending branch of (11): $\dot{y} = \sqrt{-V(y)}$. Denoting $\xi = \cos^2 y$, we get

$$ \dot{\xi} = -\sqrt{P(\xi)}, \qquad (12) $$

where

$$ \begin{aligned} P(\xi) &= 4\xi(1-\xi)(C - (B-A)\xi - A)(2h - p_x^2(B-A)\xi - p_x^2A) \\ &= 4p_x^2(B-A)^2\xi(1-\xi)(\xi-\xi_1)(\xi-\xi_3), \end{aligned} \qquad (13) $$

and $\xi_3 = \frac{C-A}{B-A} > 1$. It is easy to show that $0 \le \xi_1 < \xi_3$, and $\xi(t) \in [0, \min\{\xi_1, 1\}]$, which yields oscillating trajectories if $\xi_1 < 1$, separatrice is $\xi_1 = 1$, and rotations otherwise. Changing the sign in the right-hand part of (12) (or, equivalently, taking $t < 0$) we obtain the parameterization for the descending branch $\dot{y} = -\sqrt{-V(y)}$.

**Case** $\xi_1 < 1$. In order to reduce the integration to a known elliptic integral we will apply the method described by H. Davis in [14]. Set

$$ \eta = \sqrt{\frac{\xi(1 - \xi_1)}{\xi_1(1 - \xi)}}. \qquad (14) $$

Observe that $\eta$ is a monotone increasing function of $\xi$ and $\eta|_{\xi=0} = 0$ and $\eta|_{\xi=\xi_1} = 1$. The direct computation yields

$$ \frac{d\eta}{\sqrt{(1 - \eta^2)(m' + m\eta^2)}} = -Mdt, $$

where

$$ \begin{aligned} m &= \frac{\xi_1(\xi_3 - 1)}{\xi_3 - \xi_1} = \frac{(C - B)(2h - Ap_x^2)}{(B - A)(Cp_x^2 - 2h)}, \\ m' &= 1 - m, \\ M^2 &= (B - A)^2 p_x^2 (\xi_3 - \xi_1) = (B - A)(Cp_x^2 - 2h). \end{aligned} $$

Integrating, we obtain

$$ \eta(t) = \operatorname{cn}(Mt + \psi_0|m), \quad \psi_0 = \operatorname{cn}^{-1}(\eta(0)|m), $$

and by inverting (14), we finally get

$$ \xi(t) = \frac{\xi_1 \operatorname{cn}^2(Mt + \psi_0|m)}{1 - \xi_1 \operatorname{sn}^2(Mt + \psi_0|m)}. $$

Recall now that $\dot{x} = p_x((B-A)\cos^2 y + A) = p_x((B-A)\xi + A)$. The direct integration of this equation allows to express $x$ in terms of an elliptic integral of the third kind<fnref n="2" />:

$$ x(t) - x(0) = p_x \left[ Bt - \frac{(B-A)(1-\xi_1)}{M} \left( \Pi(\xi_1|\operatorname{am}(Mt + \psi_0|m)|m) - \Pi(\xi_1|\operatorname{am}(\psi_0|m)|m) \right) \right]. $$

$$ E(x|m) = \int_{0}^{x} \sqrt{1 - m \sin^2 u} du, \quad \Pi(\xi|x|m) = \int_{0}^{x} \frac{du}{(1 - \xi \sin^2 u)\sqrt{1 - m \sin^2 u}}. $$

<footnote n="2">²Throughout this paper we use the following definitions of the elliptic integrals of the second and of the third kind:</footnote>
**Case** $\xi_1 > 1$. In this case the parameterization formulae can be derived exactly in the same way as before using the substitution

$$ \eta = \sqrt{\frac{\xi(\xi_1 - 1)}{(\xi_1 - \xi)}} $$

which leads simply to the permutation of $\xi_1$ and 1 in all expressions. We omit the details of this computation and sum the results in the following

**Proposition 3.2.** Denote $\xi_1 = \frac{2h-Ap_x^2}{p_x^2(B-A)}$, $\xi_3 = \frac{C-A}{B-A}$. Periodic trajectories of (9) on the $(x, y)$-plane starting at the point $(x_0, y_0)$ can be parameterized as follows:

i). oscillating trajectories:

$$ \cos^2 y(t) = \frac{\xi_1 \operatorname{cn}^2(Mt + \psi_0|m)}{1 - \xi_1 \operatorname{sn}^2(Mt + \psi_0|m)}, \\ x(t) - x_0 = p_x \left[ Bt - \frac{(B-A)(1-\xi_1)}{M} \left( \Pi(\xi_1|\operatorname{am}(Mt + \psi_0|m)|m) - \Pi(\xi_1|\operatorname{am}(\psi_0|m)|m) \right) \right], $$

where

$$ m = \frac{\xi_1(\xi_3 - 1)}{\xi_3 - \xi_1}, \quad M = (B-A)p_x \sqrt{\xi_3 - \xi_1}, \quad \operatorname{cn}(\psi_0|m) = \frac{(1 - \xi_1) \cos^2 y_0}{\xi_1 \sin^2 y_0}; $$

ii). rotating trajectories:

$$ \cos^2 y(t) = \frac{\xi_1 \operatorname{cn}^2(Mt + \psi_0|m)}{\xi_1 - \operatorname{sn}^2(Mt + \psi_0|m)}, \\ x(t) - x_0 = p_x \left[ (A + (B - A)\xi_1)t - \frac{(B - A)(\xi_1 - 1)}{M} \left( \Pi(\xi_1^{-1}|\operatorname{am}(Mt + \psi_0|m)|m) - \Pi(\xi_1^{-1}|\operatorname{am}(\psi_0|m)|m) \right) \right], $$

where

$$ m = \frac{\xi_3 - \xi_1}{\xi_1(\xi_3 - 1)}, \quad M = (B - A)p_x \sqrt{\xi_1(\xi_3 - 1)}, \quad \operatorname{cn}(\psi_0|m) = \frac{(\xi_1 - 1) \cos^2 y_0}{\xi_1 - \cos^2 y_0}. $$

**Corollary 3.3.** i). Solutions of (9) are quasi-periodic trajectories on the plane $(x, y)$ for all $h \in [Ap_x^2/2, Cp_x^2/2] \setminus \{Bp_x^2/2\}$, the y-variable being $\frac{4K(m)}{M}$-periodic along oscillating trajectories and $\frac{2K(m)}{M}$-periodic along rotations, where $K(m)$ denotes the complete elliptic integral of the first kind;

ii). Two oscillating trajectories of (9) issued from the same initial point $(x_0, y_0)$ with initial conditions $(p_x, p_y(0))$ and $(p_x, -p_y(0))$ intersect after a half period $\frac{2K(m)}{M}$ at $y(\frac{2K(m)}{M}) = \pi - y_0$.

In order to shorten the expressions, in what follows we will write $\operatorname{cn}$ for $\operatorname{cn}^2(Mt + \psi_0|m)$, $\operatorname{cn}_0$ for $\operatorname{cn}^2(\psi_0|m)$, and similarly for other elliptic functions.

### 3.2.2. Conjugate times along oscillating trajectories

According to the definition, the time $t_*$ is conjugate to $t_0 = 0$ if the differential of the end-point mapping

$$ e_{(x_0, y_0)}^t : (p_x(0), p_y(0)) \mapsto (x(t), y(t)), $$

degenerates at $t = t_*$, where the extremals of the Hamiltonian system (9) are parameterized by the initial values of $p(0) = (p_x, p_y(0))$ and $p_x$ is the first integral. In what follows we denote by $t_*^1 = \min|t_*|$ the first conjugate (to 0) time.

In order to simplify the further calculation, we first make a change of variables in the space of the canonical impulses.
**Lemma 3.4.** Assume that the initial point $(x_0, y_0)$ is such that $\cos^2(y_0) \neq \xi_1$ and $h \neq Cp_x^2/2$. Then the map

$$ \Phi : (p_x, p_y(0)) \mapsto (M, \xi_1) $$

is non-degenerate.

*Proof.* Indeed, $\Phi$ can be written as a composition map $\Phi = \Phi_2 \circ \Phi_1$, where

$$ (p_x, p_y(0)) \xrightarrow{\Phi_1} (p_x, h) \xrightarrow{\Phi_2} (M, \xi_1). $$

Then

$$ \det D_{(p_x, p_y(0))} \Phi_1 = p_y(0) \left( C - \frac{z(0)}{2} \right). $$

The expression of the Hamiltonian function $H_a$ implies $p_y^2(0) = \frac{p_x^2(\xi_1 - \cos^2(y_0))}{(\xi_3 - \cos^2(y_0))}$. Hence $\Phi_1$ is non-degenerate if $p_y(0) \neq 0$, i.e., if $\cos^2(y_0) \neq \xi_1$, which is also equivalent to the condition $\psi_0 \neq 0 \mod 2K(m)$. Further, we have

$$ D_{(p_x, h)} \Phi_2 = \frac{2}{p_x^3} \sqrt{\frac{Cp_x^2 - 2h}{(B - A)}} \neq 0 $$

if $h \neq Cp_x^2/2$. The statement follows now from the chain rule for composition maps. $\square$

**Lemma 3.5.** Under the assumptions of Lemma 3.4, the conjugate times along oscillating trajectories starting at $(x_0, y_0)$ are solution to the equation

$$ \Lambda Mt = I(Mt), \qquad (15) $$

where

$$ I(T) = \int_{\psi_0}^{T+\psi_0} \frac{du}{\mathrm{sn}^2(u|m)}, \qquad \Lambda = \frac{B\xi_1}{B\xi_1 + A(1 - \xi_1)}. $$

*Proof.* Indeed, we have

$$ D_{(M, \xi_1)} e_{(x_0, y_0)}^t = \frac{\partial y}{\partial \xi} \bigg|_{y=y(t)} \Delta, $$

where

$$ \Delta = \begin{pmatrix} \frac{\partial x(t)}{\partial M} & \frac{\partial x(t)}{\partial \xi_1} \\ \frac{\partial \xi(t)}{\partial M} & \frac{\partial \xi(t)}{\partial \xi_1} \end{pmatrix}. $$

Since $\xi = \cos^2 y$,

$$ \frac{\partial y}{\partial \xi} = \frac{1}{2\sqrt{\xi(1-\xi)}} = \frac{1}{2\sqrt{\xi_1(1-\xi_1)}} \frac{1-\xi_1\mathrm{sn}^2}{\mathrm{cn}}. \qquad (16) $$

On the other hand, a straightforward but rather tedious computation<fnref n="3" /> yields

$$ \det \Delta = - \frac{t \operatorname{cn} \operatorname{dn} \operatorname{sn}}{(B-A)\sqrt{x_3 - x_1(1-\xi_1 \operatorname{sn}^2)^2}} \Delta_{Mt}, $$

where

$$ \begin{gathered} \Delta_T = B\xi_1(E_{T+\psi_0} - E_{\psi_0} + \frac{\mathrm{cn} \, \mathrm{dn}}{\mathrm{sn}} - \frac{\mathrm{cn}_0 \, \mathrm{dn}_0}{\mathrm{sn}_0}) + A(1-\xi_1)(E_{T+\psi_0} - E_{\psi_0} + \frac{\mathrm{cn} \, \mathrm{dn}}{\mathrm{sn}} - \frac{\mathrm{cn}_0 \, \mathrm{dn}_0}{\mathrm{sn}_0} - T) \\ T = Mt, \end{gathered} $$

<footnote n="3">³The authors used the Mathematica package</footnote>
where we denote $E_{\psi} = E(\mathrm{am}(\psi|m)|m)$, $E(\cdot|m)$ being the standard elliptic integral of the second kind of modulus $m$. Since

$$E_{T+\psi_0} - E_{\psi_0} + \frac{\mathrm{cn}\,\mathrm{dn}}{\mathrm{sn}} - \frac{\mathrm{cn}_0\,\mathrm{dn}_0}{\mathrm{sn}_0} = T - \int_{\psi_0}^{T+\psi_0} \frac{du}{\mathrm{sn}^2(u|m)}, \quad (17)$$

and taking into account (16), we finally found that for all initial points $(x_0, y_0)$ such that $\cos^2(y_0) \neq \xi_1$ and $h \neq Cp_x^2/2$ the conjugate times of our problem are solutions to the implicit equation

$$\mathrm{sn} \cdot \left( B\xi_1 T - (B\xi_1 + A(1-\xi_1)) \int_{\psi_0}^{T+\psi_0} \frac{du}{\mathrm{sn}^2(u|m)} \right) = 0, \quad T = Mt. \quad (18)$$

From (17) it is clear that solutions to (18) are zeros of the bracket term only, which leads to equation (15). $\square$

**Lemma 3.6.** Equation $I(\tau) = \Lambda\tau$ has exactly one root on each interval of the form $[2iK(m), 2(i+1)K(m)-\psi_0]$ (for $\tau > 0$) or $[-2iK(m), -2(i+1)K(m)-\psi_0]$ (for $\tau < 0$), $i = 1, 2, \dots$.

*Proof.* It is easy to see that $I(\cdot)$ is a monotone increasing function, $I(0) = 0$, $I(\tau) \to +\infty$ as $T \to 2iK(m) - \psi_0$ from the left and $I(\tau) \to -\infty$ as $T \to 2iK(m) - \psi_0$ from the right for $i = 0, \pm 1, \dots$. In addition, on the interval $(0, 2K(m) - \psi_0]$ the function $I(\tau)$ is positive and $I'(\tau) = \mathrm{sn}^{-2}(\tau + \psi_0|m)) \ge 1$, while $\Lambda \le 1$, the equality taking place only for $2h = Bp_x^2$, and hence this interval contains no root. Inverting the signs, we see that the interval $[-\psi_0, 0)$ contains no root either. Therefore, the positive roots of the equation $I(\tau) = \Lambda\tau$ belong to the intervals $[2iK(m) - \psi_0, 2(i+1)K(m) - \psi_0]$, $i = 1, 2, \dots$, and the negative ones to $[-2iK(m) - \psi_0, -2(i+1)K(m) - \psi_0]$, $i = 0, 1, 2, \dots$.

Let us now show that there are no roots on the intervals of the form $(2iK(m) - \psi_0, 2iK(m)]$ and of the form $[-2iK(m), -2iK(m) - \psi_0)$. Consider first the interval $(2K(m) - \psi_0, 2K(m)]$. The function $I(\cdot)$ is monotone increasing, moreover, $I(2K(m)) = 2(K(m) - E(m)) > 0$ for all $m > 0$, $E(m)$ being the complete elliptic integral of the second kind of modulus $m$. Further, observe that $\Lambda|_{\xi_1=0} = 0$, $\Lambda|_{\xi_1=1} = 1$, and moreover, $\Lambda(\xi_1)$ is a monotone increasing and strictly concave function of $\xi_1$. On the other hand, since $m = \frac{\xi_1(\xi_3-1)}{(\xi_3-\xi_1)}$, $m(\xi_1)$ is monotone increasing and strictly convex. Therefore for all $\xi_1 \in (0, 1)$ we have $\Lambda > m$. Since for all $m \in [0, 1]$ one has $mK(m) \ge K(m) - E(m)$ (the equality taking place iff $m = 0$), we deduce that the graph of $I$ lies below the line $mT$ for all $T \in (2K(m) - \psi_0, 2K(m)]$, and hence there is no root on this interval. Since $I(-2K(m)) = -2(K(m) - E(m)) < 0$, the same argument shows the absence of roots on the interval $[-2K(m), -2K(m) - \psi_0]$. Due to the $\frac{2K(m)}{M}$-periodicity of equation (15), the proof is valid for all intervals of the form $(2iK(m) - \psi_0, 2iK(m)]$ and $[-2iK(m), -2iK(m) - \psi_0)$. $\square$

Since by construction $\psi_0 \in [0, K(m)]$, we immediately get the following upper bound:

**Corollary 3.7.** The first conjugate time verifies $t_*^1 < \frac{3K(m)}{M}$.

**Lemma 3.8.** Conjugate times of the trajectories issued from $(x_0, y_0)$ such that $\cos^2(y_0) = \xi_1$ are of the form $t_*^i = \frac{2iK(m)}{M}$, $i = \pm 1, \dots$

*Proof.* The condition $\cos^2(y_0) = \xi_1$ implies $\eta(0) = 0$, and hence $\psi_0 = 0 \mod 2K(m)$. We limit ourselves to the case $\psi_0 = 0$, the other case being the simple inversion of sign. We have $\mathrm{cn}_0 = \mathrm{dn}_0 = 1$, $\mathrm{sn}_0 = 0$, and according to our previous calculation,

$$\det D_{(p_x, p_y(0))} e_{(x_0, y_0)}^t = - \sqrt{\frac{p_x^2(\xi_1 - \cos^2(y_0))}{\xi_3 - \cos^2(y_0)}} (C - \frac{z(0)}{2}) \times \frac{2}{p_x^3} \sqrt{\frac{Cp_x^2 - 2h}{(B-A)}} \times \frac{1}{2\sqrt{\xi_1(1-\xi_1)}} \frac{1-\xi_1\mathrm{sn}^2}{\mathrm{cn}} \times \frac{t \mathrm{cn} \, \mathrm{dn} \, \mathrm{sn}}{(B-A)\sqrt{x_3 - x_1(1-\xi_1\mathrm{sn}^2)^2}} \Delta_{Mt}, \quad (19)$$
where

$$ \Delta_{Mt} = \frac{1}{\mathrm{sn}_0} \left[ B \mathrm{sn}_0 \xi_1 \left( E(T + \psi_0|m) - E(\psi_0|m) + \frac{\mathrm{cn} \, \mathrm{dn}}{\mathrm{sn}} \right) - B \xi_1 \mathrm{cn}_0 \mathrm{dn}_0 + A(1 - \xi_1) \mathrm{sn}_0 \left( E(T + \psi_0|m) - E(\psi_0|m) + \frac{\mathrm{cn} \, \mathrm{dn}}{\mathrm{sn}} - T \right) - A(1 - \xi_1) \mathrm{cn}_0 \mathrm{dn}_0 \right]. $$

On the other hand, since

$$ \mathrm{sn}_0^2 = \frac{\xi_1 - \cos^2(y_0)}{\xi_1(1 - \cos^2(y_0))}, \quad \sqrt{\frac{\xi_1 - \cos^2(y_0)}{\xi_3 - \cos^2(y_0)}} = \mathrm{sn}_0 \sqrt{\frac{\xi_1(1 - \xi_1)}{\xi_3 - \xi_1}}. $$

Substituting these expressions into (19), we see that $\det D_{p_x, p_y(0)} e_{(x_0, y_0)}^t = 0$ if and only if

$$ t(B\xi_1 + A(1 - \xi_1)) \mathrm{sn}(Mt|m) = 0. $$

Thus the conjugate times for the solutions verifying $\cos^2(y_0) = \xi_1$ are of the form $\frac{2iK(m)}{M}$, $i = 1, 2, \dots$ $\square$

Summing up, we get the following

**Theorem 3.9.** 1. The first conjugate time (to 0) along any oscillating solution of the Hamiltonian system (9) satisfies $\frac{2K(m)}{M} \le t_*^1 < \frac{3K(m)}{M}$, and moreover, $t_*^1 = \frac{2K(m)}{M}$ along trajectories starting with $p_y(0) = 0$;

2. If $p_x \in [\sqrt{2hA^{-1}}, \sqrt{2hB^{-1}}]$, then

$$ \min_{p_x} t_*^1 = \frac{\pi \sqrt{A}}{\sqrt{2h(B-A)(C-A)}}, \quad \lim_{p_x \to \sqrt{\frac{2h}{B}}} t_*^1 = +\infty. $$

*Proof.* The first point summarizes the results of Lemmata 3.4-3.8. For the second point observe that $\frac{d\xi_1}{dp_x} = \frac{-4h}{p_x^3}$. On the other hand, $\xi_1 \ge m$ and

$$ \frac{d}{d\xi_1} \left( \frac{2K(m(\xi_1))}{M(\xi_1)} \right) = \frac{E(m) - K(m) + \xi_1 K(m)}{(1 - \xi_1)\xi_1(B - A)\sqrt{p_x^2(\xi_3 - \xi_1)}} \ge \frac{E(m) - K(m) + mK(m)}{(1 - \xi_1)\xi_1(B - A)\sqrt{p_x^2(\xi_3 - \xi_1)}} \ge 0. $$

Hence $\frac{2K(m)}{M}$ is a decreasing function of $p_x$ if $p_x > 0$, and the statement follows. $\square$

### 3.2.3. Conjugate times along rotating trajectories

The computation of the conjugate times along rotating trajectories is essentially the same as in the oscillating case but using the second part of Proposition 3.2. We omit technical details of this computation and present directly the result:

**Lemma 3.10.** The conjugate times along rotating trajectories are solutions to the equation

$$ BT\xi_1 - (B\xi_1 - A(\xi_1 - 1)) \int_{\psi_0}^{\psi_0+T} \mathrm{sn}^2(u|m)du = 0, \quad T = Mt. \qquad (20) $$

**Corollary 3.11.** Rotating trajectories have no conjugate points on the plane $(x, y)$.
*Proof.* Indeed, rewriting the left hand side of (20) we get

$$B\xi_1 \int_{\psi_0}^{\psi_0+T} \mathrm{cn}^2(u|m)du + A(\xi_1 - 1) \int_{\psi_0}^{\psi_0+T} \mathrm{sn}^2(u|m)du = 0,$$

which clearly has no non-trivial solution.

<figure><img src="image_1.png" /></figure>

### 3.3. Geometric analysis using the (Darboux) polar representation of the metric

In this section we reduce our analysis to the study of the Riemannian metric on a surface of revolution, see [24] for the general framework and [6, 25] for the case of two spheres of revolution.

#### 3.3.1. Riemannian metric on a two-dimensional surface of revolution: generalities

Taking a local chart in an open domain $U$, the Riemannian metric on a two-dimensional surface of revolution can be written in polar coordinates in the form

$$ds^2 = d\varphi^2 + \mu(\varphi)d\theta^2,$$

where $\mu(\varphi) > 0$. We will use the Hamiltonian formalism to describe the geodesics of such a metric. Denote by $q = (\varphi, \theta) \in U$ the state variables and by $p = (p_\varphi, p_\theta)$ the corresponding coordinates on the fibers. In addition, at each point $z = (p, q) \in T^*U$ we denote by $\partial_\varphi$ and $\partial_\theta$ the vertical and the horizontal parts of the linear space $T_z(T^*U)$. Then the Liouville form reads $\omega = pdq$. The Hamiltonian associated to the metric is given by

$$H = \frac{1}{2} \left( p_{\varphi}^{2} + \frac{p_{\theta}^{2}}{\mu(\varphi)} \right),$$

and the geodesics are projections on $U$ of the extremal curves solutions to the Hamiltonian system

$$\dot{\varphi} = p_{\varphi}, \quad \dot{p}_{\varphi} = \frac{1}{2} \frac{\mu'(\varphi)}{\mu^2(\varphi)}, \quad \dot{\theta} = \frac{p_{\theta}}{\mu(\varphi)}, \quad \dot{p}_{\theta} = 0. \qquad (21)$$

Fixing $H = \frac{1}{2}$ one obtains the geodesics parameterized by the arc-length. The condition $p_\theta = \text{const}$ defines the Clairaut relation on a surface of revolution. There are two particular types of solutions of (21):

- meridian curves: $p_\theta = 0$, $\theta(t) = \theta_0$;

- parallel solution: $\dot{\varphi}(0) = p_\varphi(0) = 0$ and $\varphi(t) = \varphi_0$, which, taking into account the Hamiltonian equations, is equivalent to $\mu'(\varphi) = 0$.

Solutions of (21) can be seen as trajectories of the one-parametric families of mechanical systems $\dot{\varphi}^2 + V(p_\theta, \varphi) = 1$, where $V(p_\theta, \varphi) = p_\theta^2/\mu(\varphi)$ is a potential. Note that the parallel solutions correspond to local extrema of $V(p_\theta, \cdot)$ for each fixed value of the Clairaut constant $p_\theta$.

In what follows we assume that

(A1) $\mu'(0) = 0$, i.e., $\varphi = 0$ is a parallel solution (called the equator);

(A2) The metric is reflectionally symmetric with respect to the equator: $\mu(\varphi) = \mu(-\varphi)$, and $\mu''(0) < 0$.

Then in a neighborhood of $\varphi = 0$ on has a family of periodic trajectories solution to the equation

$$\dot{\varphi}^2 = 1 - \frac{p_\theta^2}{\mu(\varphi)},$$

which describes the evolution of the $\varphi$ variable along the geodesics parameterized by arc-length. Denote

$$g(\varphi, p_\theta) = \sqrt{\frac{\mu(\varphi)}{\mu(\varphi) - p_\theta^2}}.$$
Then $\dot{\varphi} = \pm 1/g(\varphi, p_\theta)$, where the sign "+" describes the ascending branch of the trajectory $(\varphi(t), \theta(t))$.

In what follows, taking into account the symmetry of the problem, we assume $\theta_0 = 0$ and $p_\theta \in (0, \mu(\varphi_0))$. Let us consider the ascending branch $\dot{\varphi} > 0$, and moreover, assume $\varphi_0 = 0$, the generalization being straightforward. Since the arc-length parameterized geodesics verify the equations

$$ \dot{\varphi} = \frac{1}{g(\varphi, \theta)}, \qquad \dot{\theta} = \frac{p_{\theta}}{\mu(\varphi)}, $$

we can parameterize $\theta$ by $\varphi$, which yields

$$ \frac{d\theta}{d\varphi} = \frac{g(\varphi, p_{\theta})p_{\theta}}{\mu(\varphi)} = \frac{p_{\theta}}{\sqrt{\mu(\varphi)\sqrt{\mu(\varphi) - p_{\theta}^2}}} = f(\varphi, p_{\theta}). \qquad (22) $$

Let us now recall the main known results concerning the Jacobi equation and the conjugate locus for the Darboux-type metrics. Denote by $\vec{H}$ the Hamiltonian vector field associated to $H$, and by $e^{t\vec{H}}$ the corresponding Hamiltonian flow in $T^*U$. The function $H$ defines a quadratic form on the cotangent bundle. Computing its variation along the solutions parameterized by the arc-length we get $\langle p, \delta p \rangle = 0$.

**Definition 3.12.** Let $t \to \gamma(t)$, $t \in [0, T]$, be a solution of (21). The variational system associated to (21) along $\gamma(\cdot)$ is called the Jacobi equation. The Jacobi field $J: [0, T] \to T_{\gamma(\cdot)}(T^*U)$ is a non-trivial solution of the Jacobi equation $J(t) = e_*^{t\vec{H}}J(0)$.

Assume that $J(t) = (\delta q(t), \delta p(t))$ is a Jacobi field vertical at $t=0$, i.e., $\delta q(0) = 0$, and tangent to the level set of $H = H^{-1}(1/2) : \langle p(0), \delta p(0) \rangle = 0$. Note that since $e^{t\vec{H}}$ preserves the restriction of $\omega$ on $H$ and $\omega(J(0)) = 0$, we have $\omega(J(t)) = 0$ for all $t$. A complementary Jacobi field $P(t)$, transversal to $H$ and called the Poincaré field, can be obtained by considering the vertical curve $\lambda \to (q_0, \exp \lambda p_0)$ and its infinitesimal generator $V = p \partial_p$. Since $[\vec{H}, V] = -\vec{H}$, we have $P(t) = e_*^{t\vec{H}}V = V - t[\vec{H}, V] = V + t\vec{H}$. Hence $\omega(P(t)) = tp \frac{\partial H}{\partial p} = 2tH(p, q) = t \neq 0$ for $t > 0$. Thus, denoting $\pi: (q, p) \to q$ the standard projection, and by $\pi_*$ its differential, we get the following result:

**Lemma 3.13.** Conjugate points are given by $\pi_*(J(t_*)) = 0$, where $\pi_*(J(t)) = (\frac{\partial \varphi(t, p_\theta)}{\partial p_\theta}, \frac{\partial \theta(t, p_\theta)}{\partial p_\theta})$. In addition, for all $t$ we have the collinearity relation

$$ p_{\varphi} \frac{\partial \varphi}{\partial p_{\theta}} + p_{\theta} \frac{\partial \theta}{\partial p_{\theta}} = 0. \qquad (23) $$

In order to describe the conjugate locus to a point on the equator, we first recall the following result (see [24], Corollary 7.2.1):

**Lemma 3.14.** Let $p_\theta \in [0, \sqrt{\mu(\varphi(0))})$ be such that $\dot{\varphi} > 0$ on $(0, t)$. Then this interval contains no conjugate point.

**Definition 3.15.** Let $I = (a, \sqrt{\mu(\varphi(0))})$, $a \ge 0$, be an interval such that for any $p_\theta \in I$ the ascending branch starting at $\varphi_0 = 0$ one has $\dot{\varphi}(\frac{T}{4}) = 0$, and let $\varphi_+ = \varphi(\frac{T}{4})$. Then the trajectory $t \to \varphi(t)$ is periodic with period $T$ given by

$$ \frac{T}{4} = \int_{0}^{\varphi_{+}} g(p_{\theta}, \varphi) d\varphi, $$

its first return to the equator occurs at the time $\frac{T}{2}$, while $\theta$ changes by the quantity

$$ \Delta\theta(p_{\theta}) = \theta(\frac{T}{2}) - \theta_0 = 2 \int_{0}^{\varphi_{+}} f(\varphi, p_{\theta}) d\varphi. $$
The mapping $p_{\theta} \in I \to T(p_{\theta})$ is called the period mapping and $R: p_{\theta} \to \Delta\theta(p_{\theta})$ is called the first return mapping.

**Definition 3.16.** The extremal flow is called tame on $I$ if the first return mapping verifies $R'(p_{\theta}) < 0$.

**Proposition 3.17.** An extremal curve of a tame flow corresponding to $p_{\theta} \in I$ contains no conjugate point on $(0, T(p_{\theta})/2)$.

*Proof.* As it was shown in [6], if $R' < 0$, then the extremal curves issued from $\varphi_0 = 0$ with $p_{\theta} \in I$ do not intersect before coming back to the equator. As conjugate points are limits of the intersecting geodesics, conjugate points cannot occur before the return to the equator. $\square$

The minimal distance to the conjugate locus can be estimated in terms of the Gauss curvature of the metric, which in polar coordinates can be computed via the formula

$$G(\varphi) = -\frac{1}{\sqrt{\mu(\varphi)}} \frac{\partial^2 \sqrt{\mu(\varphi)}}{\partial \varphi^2}.$$

The following fact is standard and follows from the classical result by Poincaré [22]:

**Lemma 3.18.** Assume that $G(0) = \max_{\varphi} G(\varphi) > 0$. Then the first conjugate point along the equator occurs at $t_*^1 = \pi/\sqrt{G(0)}$. Moreover, $(\varphi(t_*), \theta(t_*))$ realizes the minimum distance from $(0,0)$ to the conjugate locus $C_{(0,0)}$, which has a cusp at this point.

Assume now $p_{\theta} \in I$ and consider $t \in (T(p_{\theta})/2, 3T(p_{\theta})/4)$. According to our assumptions, on this interval we have $\varphi(t) < 0$ and $\dot{\varphi}(t) < 0$. Thus

$$\theta(t, p_{\theta}) = \Delta\theta(p_{\theta}) + \int_{\frac{T}{2}}^{t} \frac{p_{\theta} d\tau}{\mu(\varphi(\tau))} = \Delta\theta(p_{\theta}) + \int_{\varphi(t, p_{\theta})}^{0} f(\zeta, p_{\theta}) d\zeta.$$

The collinearity relation (23) implies that

**Lemma 3.19.** For $p_{\theta} \in I$ the conjugate time in the open interval $(\frac{T}{2}, \frac{T}{2} + \frac{T}{4})$ is solution to the equation

$$\frac{\partial \theta}{\partial p_{\theta}}(\varphi, p_{\theta}) = 0, \qquad (24)$$

where

$$\theta(\varphi, p_{\theta}) = \Delta\theta(p_{\theta}) + \int_{\varphi}^{0} f(\zeta, p_{\theta}) d\zeta.$$

To conclude we recall the following relation between the period and the first return map, which we will use in our further computations:

**Lemma 3.20.**

$$R'(p_{\theta}) = \frac{T'(p_{\theta})}{2p_{\theta}}.$$
### 3.3.2. The Darboux normal form of the Serret-Andoyer metric

We now apply the described technique to the analysis of the Serret-Andoyer metric. As before, we assume that $A < B < C$. In order to put metric (8) into the Darboux normal form $g = d\varphi^2 + \mu(\varphi)d\theta^2$, we have to integrate the following equations

$$ \frac{dy}{\sqrt{C - (A \sin^2 y + B \cos^2 y)}} = d\varphi, \quad d\theta = dx. $$

The second equation implies that $\theta = x$ modulo a rotation by a fixed angle. As for the $\varphi$ variable, according to the argument of Section 2, we choose the initial condition at $y_0 = \pi/2$.

Set $X = \sin y$. Thus $X(0) = 1$ and

$$ \int_{\varphi_0}^{\varphi} d\varphi = \int_{1}^{X} \frac{dX}{\sqrt{(1-X^2)((C-B)+(B-A)X^2)}} = -\int_{X}^{1} \frac{dX}{\sqrt{(C-A)(1-X^2)(\frac{C-B}{C-A}+\frac{B-A}{C-A}X^2)}} \\ = -\int_{X}^{1} \frac{dX}{\sqrt{(C-A)(1-X^2)(k'+kX^2)}}, $$

where

$$ k = \frac{(B-A)}{(C-A)} = \frac{1}{\xi_3}, \quad k' = 1-k = \frac{(C-B)}{(C-A)}. $$

Denote $\alpha = \sqrt{C-A}$. Choosing the initial condition $\varphi_0 = 0$ we finally obtain a standard elliptic integral

$$ -\alpha\varphi = \int_{X}^{1} \frac{dX}{\sqrt{(1-X^2)(k'+kX^2)}} = \text{cn}^{-1}(X|k), $$

and hence

$$ \sin y = X = \text{cn}(-\alpha\varphi|k) = \text{cn}(\alpha\varphi|k), $$

which implies

$$ z = 2(A \sin^2 y + B \cos^2 y) = 2(A \operatorname{cn}^2(\alpha\varphi|k) + B \operatorname{sn}^2(\alpha\varphi|k)), \\ \mu(\varphi) = \frac{2}{z} = \frac{1}{A \operatorname{cn}^2(\alpha\varphi|k) + B \operatorname{sn}^2(\alpha\varphi|k)}. $$

So, we get the following

**Proposition 3.21.** i). The Andoyer-Serret metric $g_a$ can be put into the Darboux normal form $d\varphi^2 + \mu(\varphi)d\theta^2$ with $\mu(\varphi) = (A\operatorname{cn}^2(\alpha\varphi|k) + B\operatorname{sn}^2(\alpha\varphi|k))^{-1} \in [B^{-1}, A^{-1}]$, $\alpha = \sqrt{C-A}$.

ii). The Gauss curvature of the Andoyer-Serret metric is given by

$$ G(\varphi) = \frac{(A + B + C)(z(\varphi) - z_{-})(z(\varphi) - z_{+})}{z(\varphi)^2}, $$

where $z(\varphi) = 2(A\operatorname{cn}^2(\alpha\varphi|k) + B\operatorname{sn}^2(\alpha\varphi|k))$ and

$$ z_{\pm} = \frac{2(AB + AC + BC \pm \sqrt{(AB + AC + BC)^2 - 3ABC(A + B + C)})}{A + B + C}. $$

The second statement of the proposition can be verified by a straightforward computation using the Gauss curvature formula mentioned above. In particular, it follows that $G$ is $\frac{2K(k)}{\alpha}$-periodic and reflectionally symmetric: $G(\varphi) = G(-\varphi)$. It also easy to show that $z_{-} \in [2A, 2B]$, while $z_{+} > 2B$, thus on $[-\frac{K(k)}{\alpha}, \frac{K(k)}{\alpha}]$ the
Gauss curvature changes its sign twice at $\pm\varphi_1$, $\varphi_1 = \alpha^{-1}\text{sn}^{-1}(\sqrt{\frac{z_- - 2A}{2(B-A)}})$. In addition, it has local maxima at $\varphi = \pm\frac{K(k)}{\alpha}$ and $\varphi = 0$, and a minimum at the point $\varphi_*$, which is the solution to the equation $z(\varphi_*) = \frac{2z_-z_+}{z_- + z_+}$. More precisely,

$$ \max_{[0, \frac{2}{\alpha} K(k)]} G(\varphi) = G(0) = \frac{(B-A)(C-A)}{A}, \qquad \min_{[0, \frac{2}{\alpha} K(k)]} G(\varphi) = -\frac{(A+B+C)(z_+z_-)^2}{4z_-z_+}. $$

The graph of $G$ is represented by a continuous curve in the left subplot of Fig. 1. The right subplot of Fig. 1 shows the extremal trajectories of the Serret-Andoyer metric on the extended interval $p_\theta \in [0, A^{-1/2}]$. The thick dashed curves describe the permanent rotations around the minor axis of inertia, while the thick continuous curves correspond to the separating polhodes. Note that physically realizable solutions of the Euler-Poisot problem concern $p_\theta \in [\sqrt{C^{-1}}, \sqrt{A^{-1}}]$ (trajectories in the horizontal sector bounded by the thick dashed curves). The Gauss curvature is positive in the gray stripe along the horizontal axis. It changes sign along rotational trajectories, while it is positive along oscillating trajectories which remain sufficiently close to the horizontal axis. Such trajectories correspond to the polhodes around the major axis of the energy ellipsoid.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. The function $\mu(\varphi)$ (dashed curve) and Gauss curvature (continuous curve) of the Serret-Andoyer metric (left); Extremal curves of the Serret-Andoyer metric (right).</figcaption></figure>

### 3.3.3. Conjugate locus of the Serret-Andoyer metric

We limit our analysis to the physically interesting case $p_\theta \in [\sqrt{C^{-1}}, \sqrt{A^{-1}}]$, and consider the geodesics parameterized by the arc-length ($h = \frac{1}{2}$). According to Corollary 3.11, conjugate points belong to oscillating trajectories only. In order to put the analysis of the Serret-Andoyer case in the same framework as in Section 3.3.1, we set $\theta \equiv x$ and $p_\theta \equiv p_x$. The general results in Section 3.3.1 can be extended to the Serret-Andoyer case, provided the symmetry assumptions $(A_1)$ and $(A_2)$ are verified.

Using the explicit form of the function $\mu(\varphi)$ given in Proposition 3.21 (dashed curve in left subplot of Fig. 1), we immediately derive the following symmetry properties of $g_a$:

**Proposition 3.22.** The Serret-Andoyer metric $g_a$ is invariant with respect to rotation by angle $\theta$ ($p_\theta = \text{const}$), it possesses the reflectional symmetry induced by $\mu(-\varphi) = \mu(\varphi)$, and $\mu'(\varphi) \neq 0$ for $\varphi \in (0, K(k)\alpha^{-1})$.

In particular, assumptions $(A_1)$, $(A_2)$ are true. The period mapping of oscillating trajectories is $T = \frac{4K(m)}{M}$, where $m$ and $M$ are given in the first part of Proposition 3.2. Using the explicit parameterization formulae, the first return mapping along oscillating trajectories of the Serret-Andoyer metric can be easily computed:

$$ \Delta\theta(p_\theta) = \theta\left(t + \frac{2K(m)}{M}\right) - \theta(t) = \frac{2p_\theta}{M}[K(m)B - (B-A)(1 - \xi_1)\Pi(\xi_1|m)], $$
Π(·|m) being a complete elliptic integral if the third kind of modulus m.

**Proposition 3.23.** For the non-negative values of the Clairaut constant $p_θ$ such that $1/B < p_θ^2 < 1/A$ the first return mapping of the Serret-Andoyer metric is a strictly convex monotone decreasing function. In particular, the corresponding geodesic flow is tame.

*Proof.* Denoting ′ the derivatives with respect to $p_θ$ and differentiating the period map $T = \frac{4K(m)}{M}$ we get

$$ \frac{1}{4} \frac{dT}{dp_\theta} = \frac{m'}{M} \frac{dK(m)}{dm} - \frac{K(m)M'}{M^2}, $$

$$ \frac{1}{4} \frac{d^2 T}{dp_\theta^2} = \frac{K(m)}{M^3} (2M'^2 - M M'') + \frac{1}{M^2} \frac{dK(m)}{dm} (Mm'' - 2m'M') + \frac{d^2 K(m)}{dm^2} \frac{m'^2}{M}. $$

By definition, $K(m) = \int_0^{\frac{\pi}{2}} \frac{dx}{\sqrt{1-m \sin^2 x}}$, hence for all $m \in [0, 1] \frac{dK(m)}{dm} > 0$ and $\frac{d^2 K(m)}{dm^2} > 0$. Further, if $I_2 < p_\theta^2 \le I_1$, then

$$ m' = - \frac{2(C - A)(C - B)p_\theta}{(B - A)(Cp_\theta^2 - 1)^2} < 0, \quad M' = \frac{Cp_\theta\sqrt{B - A}}{\sqrt{Cp_\theta^2 - 1}} > 0, $$

and hence $T$ is a monotone decreasing function if $p_\theta^2 \in (1/B, 1/A)$. In addition, we have

$$ Mm'' - 2m'M' = \frac{2(C - A)(C - B)(1 + 5Cp_\theta^2)}{\sqrt{B - A(Cp_\theta^2 - 1)^5/2}} > 0, $$

$$ 2M'^2 - MM'' = \frac{(B - A)C(1 + 2Cp_\theta^2)}{Cp_\theta^2 - 1} > 0, \quad m'^2 > 0. $$

Hence $T$ is a strictly convex function of $p_θ$. The statement now follows from Lemma 3.20. □

As usual, the dynamics of $\varphi$-variable is obtained by integrating the natural mechanical system $\dot{\varphi}^2 + \frac{z(\varphi)}{2} p_\theta^2 = 1$, whose potential is given by $V(\varphi, p_\theta) = p_\theta^2 z(\varphi)/2$, and hence $V(\cdot, p_\theta)$ is an increasing function of $p_\theta$. Varying $p_\theta$, we obtain a one parametric family of potentials represented in Fig. 2. The gray zone corresponds to $p_\theta \in [\sqrt{B^{-1}}, \sqrt{A^{-1}}]$, where the limit value $p_\theta = \sqrt{B^{-1}}$ describes the separatrice. Oscillating trajectories starting at $\varphi_0$ are generated by the values $p_\theta \in (\sqrt{B^{-1}}, \sqrt{2z(\varphi_0)^{-1}}]$ (dark gray stripe).

<figure><img src="image_3.png" /><figcaption>FIGURE 2</figcaption></figure>
**Theorem 3.24.** The first conjugate locus to a point $(0, \varphi_0)$ of the Serret-Andoyer metric $g_a$ consist of two components, symmetric with respect to the vertical line $\theta = 0$. Each component is composed by two smooth branches, which asymptotically tend to the horizontal lines $\varphi = \pm \frac{K(k)}{\alpha}$ and form a unique horizontal cusp on the line $\varphi = -\varphi_0$.

*Proof.* The vertical symmetry of the conjugate locus can be easily derived from (10), and it is easy to see that the left component of the locus corresponds to the negative values of $p_\theta$. Let us focus on the right component with $p_\theta > 0$.

Consider the ascending branch $(\dot{\varphi}(0) > 0)$ of the trajectory starting at $\varphi_0 \in [-K(k)\alpha^{-1}, 0]$, $\theta(0) = 0$. According to Theorem 3.9, the first conjugate time along any trajectory $t_1^* \in [\frac{T}{2}, \frac{3T}{4})$. Taking into account that $\dot{\varphi} = 0$ at $\varphi = \pm \frac{K(k)}{\alpha}$, on the interval $[0, \frac{3T}{4})$ we can parameterize $\theta$ by $\varphi$ as follows:

$$ \theta(\varphi, p_{\theta}) = \int_{\varphi_0}^{0} f(\zeta, p_{\theta})d\zeta + \Delta\theta(p_{\theta}) - \int_{0}^{\varphi(p_{\theta})} f(\zeta, p_{\theta})d\zeta. $$

The conjugacy condition $\frac{\partial\theta(\varphi, p_\theta)}{\partial p_\theta} = 0$ implies

$$ - \int_{0}^{\varphi_*(p_{\theta})} f'(\zeta, p_{\theta})d\zeta = - \int_{\varphi_0}^{0} f'(\zeta, p_{\theta})d\zeta - \Delta\theta'(p_{\theta}), $$

where $\varphi_*(p_\theta) = \varphi(t_1^*(p_\theta))$. Computing, we have

$$ - \int_{\varphi_0}^{0} f'(\zeta, p_{\theta})d\zeta - \Delta\theta'(p_{\theta}) = \int_{-\frac{K(k)}{\alpha}}^{\varphi_0} f'(\zeta, p_{\theta})d\zeta - \frac{3}{2}\Delta\theta'(p_{\theta}) > 0 $$

since $f'(\zeta, p_\theta) > 0$ and $\Delta\theta$ is monotone decreasing. In particular, it follows that $\varphi_*(p_\theta) < 0$. Differentiating again we get

$$ -f'(\varphi_*(p_\theta), p_\theta)\varphi'_*(p_\theta) = -\Delta\theta''(p_\theta) - \int_{\varphi_0}^{0} f''(\zeta, p_\theta)d\zeta + \int_{0}^{\varphi_*(p_\theta)} f''(\zeta, p_\theta)d\zeta < 0, $$

since $\Delta\theta(p_\theta)$ is strictly convex, $f''(\zeta, p_\theta) > 0$ and $\varphi_*(p_\theta) < 0$. Therefore $\varphi'_*(p_\theta) > 0$, which implies that the conjugate locus can be smoothly parameterized by $p_\theta$.

Denote $\gamma = (\theta(\varphi_*(p_\theta), p_\theta), \varphi_*(p_\theta))$ such a parameterization of the conjugate locus. Then its tangent vector is given by $\varphi'_*(\frac{\partial\theta}{\partial\varphi}\partial\theta + \partial_\varphi)$, where $p_\theta \in (\sqrt{B^{-1}}, \sqrt{2/z(\varphi_0)})$. The global structure of the conjugate locus can be understood from the asymptotic behaviour of the curve $\gamma$ as $p_\theta$ tends to its limits.

a). Case $p_\theta \to \frac{1}{\sqrt{B}}$. Then

$$ \lim_{p_{\theta} \to \frac{1}{\sqrt{B}} +} f(\varphi_*(p_{\theta}), p_{\theta}) = \frac{A \cosh^2(\alpha \varphi_* |k|) + B \sinh^2(\alpha \varphi_* |k|)}{\sqrt{(B-A)\cosh^2(\alpha \varphi_* |k|)}}. $$

By definition, $\varphi_* = \varphi(t_1^*)$. Since when $p_\theta \to \frac{1}{\sqrt{B}}$, $\varphi(t) \to \varphi_s(t)$, where $\varphi_s$ describes the dynamics of the $\varphi$ variable along the separatrice, we have $\lim_{p_\theta \to \frac{1}{\sqrt{B}} +} \varphi_* = \varphi_s(t_1^*) = \pm \frac{K(k)}{\alpha}$, because, according to Theorem 3.9,
$t_1^* \to +\infty$ as $p_\theta \to \frac{1}{\sqrt{B}} +$. Therefore, since $\frac{\partial \theta}{\partial \varphi} = \pm f(\varphi_*(p_\theta), p_\theta) \to \pm \infty$ as $p_\theta \to \frac{1}{\sqrt{B}} +$, the extremities of the conjugate locus asymptotically tend to the horizontal lines $\varphi = \pm \frac{K(k)}{\alpha}$.

b). Case $p_\theta \to \sqrt{\frac{2}{z(\varphi_0)}}$. Though the existence of the horizontal cusp in this case follows from [22], here we derive it directly from Proposition 3.2. We have

$$ \lim_{p_\theta \to \sqrt{\frac{2}{z(\varphi_0)}}} f(\varphi_*(p_\theta), p_\theta) = \frac{A \operatorname{cn}^2(\alpha \varphi_* |k) + B \operatorname{sn}^2(\alpha \varphi_* |k)}{\sqrt{A(\operatorname{cn}^2(\alpha \varphi_0 |k) - \operatorname{cn}^2(\alpha \varphi_* |k)) + B(\operatorname{sn}^2(\alpha \varphi_0 |k) - \operatorname{sn}^2(\alpha \varphi_* |k))}} $$

Observe that $p_\theta = \sqrt{\frac{2}{z(\varphi_0)}}$ implies $\dot{\varphi}_0 = 0$, and hence $\psi_0 = 0 \mod \frac{2K(m)}{M}$, where $\psi_0, m$ and $M$ are described in the first part of Proposition 3.2. Then, according to Theorem 3.9, $t_*^1 = \frac{2K(m)}{M}$ and $\varphi_*(p_\theta) = -\varphi_0$. Applying the formulae of the first part of Proposition 3.2, we obtain

$$ \lim_{p_\theta \to \sqrt{\frac{2}{z(\varphi_0)}}} \operatorname{cn}^2(\alpha \varphi_* |k) = \sin^2(y(t)) \Big|_{t=\frac{2K(m)}{M}} = 1 - \xi(t) \Big|_{t=\frac{2K(m)}{M}} = 1 - \xi_1. $$

On the other hand, since $\eta(0) = \operatorname{cn}^{-1}(0|m) = 1$, $\xi(0) = \xi_1$ and hence

$$ \lim_{p_\theta \to \sqrt{\frac{2}{z(\varphi_0)}}} \operatorname{cn}^2(\alpha \varphi_0 |k) = \sin^2(y_0) = 1 - \xi(0) = 1 - \xi_1. $$

Therefore $\frac{\partial \theta}{\partial \varphi} = \pm f(\varphi_*(p_\theta), p_\theta) \to \pm \infty$ as $\varphi \to -\varphi_0$, forming a horizontal cusp. $\square$

**Remark 3.25.** Let $\theta_*(\varphi_0) = \theta(t_*^1)$. We know that $\bar{t}_* = \min_{\varphi_0} t_*^1 = \frac{\pi \sqrt{A}}{\sqrt{(B-A)(C-A)}}$ (Theorem 3.9 or Lemma 3.18). Moreover, Lemma 3.18 implies

$$ \min_{\varphi_0} \theta_*(\varphi_0) = \theta_*(0) = \theta(\bar{t}_*) = \frac{\pi A \sqrt{A}}{\sqrt{(B-A)(C-A)}}. $$

## 4. OPTIMAL CONTROL OF A THREE SPIN SYSTEM WITH INEQUAL ISING COUPLING AND LEFT-INVARIANT SUB-RIEMANNIAN METRICS ON SO(3)

In this section we analyse the optimal control of three spin systems with Ising coupling introduced in [27,28] in relation with invariant SR-metrics on SO(3) [18].

### 4.1. The physical problem

Spin dynamics with application in NMR are described in [21] and the general mathematical model is described in [13]. Here we consider a system of three spins with Ising coupling analyzed in [27,28]. In a chosen rotating frame, the Hamiltonian takes the form

$$ H_0 = 2\pi(J_{12}I_{1z}I_{2z} + J_{23}I_{2z}I_{3z}), $$

where $J_{12}, J_{23}$ are the coupling Ising constants and $I_{k\alpha}$, $\alpha = x, y, z$ are the standard tensor products of the respective Pauli matrix $I_\alpha$:

$$ I_{k\alpha} = 1 \otimes \dots \otimes I_\alpha \otimes 1 \cdots \otimes 1, $$

where $I_\alpha$ appears at the k-th position and 1 is the identity matrix.
Following [28], we consider the transfer from $I_{1x}$ to $4I_{1y}I_{2y}I_{3z}$ in minimum time, which realizes an intermediate step of the fastest transfer of $I_{1x}$ to $I_{3x}$. After reduction to the transfer from $x_1$ to $x_4$ identified below, the control system can be written in the form

$$ \frac{dx(t)}{dt} = \begin{pmatrix} 0 & -1 & 0 & 0 \\ 1 & 0 & -u(t) & 0 \\ 0 & u(t) & 0 & -k \\ 0 & 0 & k & 0 \end{pmatrix} x(t), \quad k = \frac{J_{23}}{J_{12}} > 0, \qquad (25) $$

where $x = (x_1, x_2, x_3, x_4)$, $x_1 = \langle I_{1x} \rangle$, $x_2 = \langle 2I_{1y}I_{2z} \rangle$, $x_3 = \langle 2I_{1y}I_{2x} \rangle$, $x_4 = \langle 4I_{1y}I_{2y}I_{3z} \rangle$ and $\langle \mathcal{O} \rangle$ denotes the expectation value of the operator $\mathcal{O}$. The original optimal control problem then becomes the problem of transfer $(1, 0, 0, 0)$ to $(0, 0, 0, 1)$ in minimum time. Introducing the coordinates $r_1 = x_1$, $r_2 = \sqrt{x_2^2 + x_3^2}$, $r_3 = x_4$ and $\tan \alpha = x_3/x_2$, we get the system

$$ \frac{dr(t)}{dt} = \begin{pmatrix} 0 & -\cos \alpha & 0 \\ \cos \alpha & 0 & -k \sin \alpha \\ 0 & k \sin \alpha & 0 \end{pmatrix} r(t) \qquad (26) $$

where $r = (r_1, r_2, r_3) \in S^2$ and $\alpha(t)$ plays the role of control. In these coordinates the minimum time problem for (25) becomes the problem of the fastest transfer of the point $(1, 0, 0)$ to the point $(0, 0, 1)$ on the sphere.

## 4.2. The related almost-Riemannian and sub-Riemannian problems

The original optimal control problem leads to analyze the problem of transfer $r(0) = (1, 0, 0)$ to $r(T) = (0, 0, 1)$ for the system

$$ \dot{r}_1 = u_3 r_2, \quad \dot{r}_2 = -u_3 r_1 + u_1 r_3, \quad \dot{r}_3 = -u_1 r_2 \qquad (27) $$

by minimizing the functional

$$ \int_{0}^{T} (I_{1}u_{1}^{2}(t) + I_{3}u_{3}^{2}(t))dt. \qquad (28) $$

Since $r \in S^2$, this problem defines a singular metric on $S^2$ called almost-Riemannian in the literature ([1]):

$$ g = \frac{dr_1^2 + k^2 dr_3^2}{r_2^2}, \quad k^2 = \frac{I_1}{I_3} > 0. $$

We have the following

**Lemma 4.1.**

1) In the spherical coordinates $r_2 = \cos \varphi$, $r_1 = \sin \varphi \cos \theta$, $r_3 = \sin \varphi \sin \theta$ the metric $g$ takes the form

$$ g = (\cos^2 \theta + k^2 \sin^2 \theta)d\varphi^2 + 2(k^2 - 1) \tan \varphi \sin \theta \cos \theta d\varphi d\theta + \tan^2 \varphi(\sin^2 \theta + k^2 \cos^2 \theta)d\theta^2, $$

with associated Hamiltonian

$$ H = \frac{1}{4k^2} (p_{\varphi}^2(\sin^2 \theta + k^2 \cos^2 \theta) + p_{\theta}^2 \cot^2 \varphi(\cos^2 \theta + k^2 \sin^2 \theta) - 2(k^2 - 1)p_{\varphi}p_{\theta} \cot \varphi \sin \theta \cos \theta). $$

2) If $k=1$, the Hamiltonian takes the form $H = 1/4 (p_{\varphi}^2 + p_{\theta}^2 \cot^2 \varphi)$ and describes the standard Grushin metric on $S^2$.

Moreover, we have
**Lemma 4.2.** The family of metrics g depending upon the parameter k have a fixed singularity on the equator and a discrete symmetry group defined by the two reflexions: $H(\varphi, p_{\varphi}) = H(\pi - \varphi, -p_{\varphi})$ and $H(\theta, p_{\theta}) = H(-\theta, -p_{\theta})$.

Optimal control problem (4.2), (28) admits the following interpretation. Introduce the following matrix: $R(t) = (r_{ij}(t)) \in SO(3)$, and denote $r_1 = r_{11}$, $r_2 = r_{12}$, $r_3 = r_{13}$ the components of its first line. Consider the solutions of the right-invariant control system

$$ \frac{dR^t(t)}{dt} = \begin{pmatrix} 0 & u_3(t) & 0 \\ -u_3(t) & 0 & u_1(t) \\ 0 & -u_1(t) & 0 \end{pmatrix} R^t(t), \qquad (29) $$

minimizing the cost:

$$ \int_{0}^{T} (I_{1}u_{1}^{2}(t) + I_{3}u_{3}^{2}(t))dt. \qquad (30) $$

Define the two sub-manifolds of $SO(3)$:

$$ M_0 = \{R \in SO(3); R^t = (r(0), \cdot, \cdot), r(0) = (1, 0, 0)\}, $$

$$ M_1 = \{R \in SO(3); R^t = (r(T), \cdot, \cdot), r(T) = (0, 0, 1)\}. $$

Denoting by $M^\perp$ the symplectic lift of a sub-manifold $M \in SO(3)$, and using the analysis of Section 2, one has:

**Proposition 4.3.** The extremals of the almost-Riemannian metric g on $S^2$ verifying the boundary conditions $r(0) = (1, 0, 0)$, $r(T) = (0, 0, 1)$ are extremals of the sub-Riemannian problem defined by (29), (30) with $k^2 = I_1/I_3$, which satisfy the boundary conditions $(R(0), \lambda(0)) \in M_0^\perp$, $(R(T), \lambda(T)) \in M_1^\perp$, $\lambda$ being the adjoint vector.

Introducing the isothermal coordinates $x = r_1$, $y = kr_3$ of the metric g, so that $g = \lambda(x, y)(dx^2 + dy^2)$ with

$$ \lambda(x, y) = \left(1 - x^2 - \frac{y^2}{k^2}\right)^{-1} > 0, $$

we can easily compute the Gaussian curvature $G(x, y)$:

$$ G = -\frac{1}{2\lambda(x, y)}\Delta \ln(\lambda(x, y)) = -1 - \frac{1}{k^2} - 2\lambda(x, y) \left(x^2 + \frac{y^2}{k^4}\right), $$

where $\Delta = \partial_{xx} + \partial_{yy}$. This leads to the following

**Lemma 4.4.** The cut locus of the equatorial point is the equator minus zero.

*Proof.* The Gaussian curvature G is strictly negative in each hemisphere of the sphere, where the almost-Riemannian metric is Riemannian. Hence a geodesic starting from the equator cannot have a conjugate point before returning to the equator. In view of the reflectional symmetry with respect to the equator, the cut point along each geodesic starting from the equatorial point is on the equator. $\square$

In conclusion we observe, that according to [6], the conjugate and cut loci near the equator can be easily computed and the whole conjugate locus is determined numerically using a continuation method. We also remark that when this paper was already accepted for publication, we came to know about work [9], where, though in a different physical context, the sub-Riemannian problem (29),(30) was treated in detail, and the cut locus was computed using a different argument. Our results below provide additional details on the structure of the conjugate locus of the problem.
### 4.3. Computation of the extremal curves

The identification provided by Proposition 4.3 allows to parameterize the extremal curves of the Riemannian metric $g$ using the extremals of the left-invariant problem with the Hamiltonian

$$H_n = \frac{1}{2} \left( \frac{H_1^2}{I_1} + \frac{H_3^2}{I_3} \right).$$

According to Section 2.2, Euler's equations take the form

$$\frac{dH_1}{dt} = \frac{H_2 H_3}{I_3}, \quad \frac{dH_2}{dt} = H_1 H_3 \left( \frac{1}{I_1} - \frac{1}{I_3} \right), \quad \frac{dH_3}{dt} = -\frac{H_1 H_2}{I_1}.$$

Fixing $H_n = 1/2$ and introducing $\alpha$ such that $\cos \alpha = H_1 I_1^{-1/2}$, $\sin \alpha = H_3 I_3^{-1/2}$, we get the pendulum equation

$$\frac{d^2 \alpha}{dt^2} = \frac{k^2 - 1}{2I_1} \sin 2\alpha, \quad k^2 = \frac{I_1}{I_3}, \qquad (31)$$

which can be easily integrated. Then, using (4), one can compute the Euler angles $\Phi_2$ and $\Phi_3$, and finally $R(t)$ can be found by integrating the remaining equation for $\Phi_1$:

$$\frac{d\Phi_1}{dt} = \frac{G}{H_2 \sin \Phi_3 - H_1 \cos \Phi_3} \left( \sin \Phi_3 \frac{\partial H_n}{\partial H_2} - \cos \Phi_3 \frac{\partial H_n}{\partial H_1} \right) = \frac{GH_1^2}{I_1(H_1^2 + H_2^2)}, \qquad (32)$$

where $G = |H|$ denotes the angular momentum constant,

The integration of the pendulum equation (31) is a standard exercise of Classical Mechanics (see for instance in [20]). In our case the final result depends on the sign of $k^2 - 1$. i.e., on the sign of $I_1 - I_3$. In what follows we assume $k \neq 1$ (or equivalently $I_1 \neq I_3$). Then (31) has non-trivial periodic solutions, and in particular, if $I_1 < I_3$, we get the well-known pendulum phase portrait with stable equilibria at $0 \mod \pi$ and the instable ones at $\alpha = \pm \pi/2 \mod \pi$. In the case $I_3 < I_1$ this picture is shifted by $\pi/2$, so that the stable equilibria are at $\alpha = \pi/2 \mod \pi$.

Setting $\nu = 2\alpha$, we can rewrite (31) as a natural mechanical system

$$\frac{1}{2}\dot{\nu}^2 + \frac{I_1 - I_3}{I_1 I_3} \cos \nu = c, \qquad (33)$$

which has solutions if and only if the energy constant verifies $c \geq -\frac{|I_1 - I_3|}{I_1 I_3}$. The value of $c$ is completely defined by $I_1, I_3$ and $G$. Indeed, Euler's equations imply $H_2 = -\sqrt{I_1 I_3} \dot{\alpha}$. Recalling that $G^2 = H_1^2 + H_2^2 + H_3^2$ and using the definition of $\alpha$, we get

$$c = \frac{2G^2 - I_1 - I_3}{I_1 I_3}.$$

The standard analysis of the pendulum phase portrait now yields

**Proposition 4.5.** *Pendulum equation (31) admits oscillating solutions if $I_1 > G^2 > I_3$ or $I_3 > G^2 > I_1$, and rotating solutions if $G^2 > I_1 > I_3$ or $G^2 > I_3 > I_1$.*

We skip the technical details of our final computations since they are rather standard, and present directly the parametrization formulae for the extremals of (29). As before, here we extend the terms "oscillating/rotating" to the extremals in the full space taking into account the behavior of the variable $\alpha$.

**Case $k < 1$:** we have $I_1 < I_3$, and denoting $\omega = \sqrt{\frac{I_3 - I_1}{I_1 I_3}}$ we get the following parametrization:
- oscillating extremals:

$$ \sin \alpha(t) = \sqrt{m} \operatorname{sn}(\omega t + \psi_0 | m), \quad \psi_0 = \operatorname{sn}^{-1}(\sin \alpha(0) m^{-1/2} | m), \quad m = \frac{G^2 - I_1}{I_3 - I_1}, $$

$$ \Phi_1(t) - \Phi_1(0) = \frac{G}{I_3} \left[ t + \frac{I_3 - G^2}{\omega G^2} \left( \Pi\left(\frac{I_3 m}{G^2} \middle| \operatorname{am}(\omega t + \psi_0 | m) \middle| m\right) - \Pi\left(\frac{I_3 m}{G^2} \middle| \operatorname{am}(\psi_0 | m) \middle| m\right) \right) \right], $$ $$ \cos \Phi_2 = \frac{\sqrt{I_3 m}}{G} \operatorname{sn}(\omega t + \psi_0 | m), \quad \tan \Phi_3 = \sqrt{I_3 m \omega} \frac{\operatorname{cn}(\omega t + \psi_0 | m)}{\operatorname{dn}(\omega t + \psi_0 | m)}. $$

- rotating extremals:

$$ \sin \alpha(t) = \operatorname{sn}(\bar{\omega} t + \psi_0 | m), \quad \psi_0 = \operatorname{sn}^{-1}(\sin \alpha(0) | m), \quad m = \frac{I_3 - I_1}{G^2 - I_1}, \quad \bar{\omega} = \frac{\omega}{\sqrt{m}}, $$

$$ \Phi_1(t) - \Phi_1(0) = \frac{G}{I_3} \left[ t + \frac{I_3 - G^2}{\bar{\omega} G^2} \left( \Pi\left(\frac{I_3}{G^2} \middle| \operatorname{am}(\bar{\omega} t + \psi_0 | m) \middle| m\right) - \Pi\left(\frac{I_3}{G^2} \middle| \operatorname{am}(\psi_0 | m) \middle| m\right) \right) \right], $$ $$ \cos \Phi_2 = \frac{\sqrt{I_3}}{G} \operatorname{sn}(\bar{\omega} t + \psi_0 | m), \quad \tan \Phi_3 = \sqrt{I_3 \bar{\omega}} \frac{\operatorname{dn}(\bar{\omega} t + \psi_0 | m)}{\operatorname{sn}(\bar{\omega} t + \psi_0 | m)}. $$

Case $k > 1$: we have $I_1 > I_3$, and setting $\omega = \sqrt{\frac{I_1 - I_3}{I_1 I_3}}$ we obtain the following formulae:

- oscillating extremals:

$$ \cos \alpha(t) = -\sqrt{m} \operatorname{sn}(\omega t + \psi_0 | m), \quad \psi_0 = -\operatorname{sn}^{-1}(\cos \alpha(0) m^{-1/2} | m), \quad m = \frac{G^2 - I_3}{I_1 - I_3}, $$

$$ \Phi_1(t) - \Phi_1(0) = \frac{G}{I_3} \left[ t - \frac{1}{\omega} \left( \Pi\left(\frac{I_3}{I_3 - I_1} \middle| \operatorname{am}(\omega t + \psi_0 | m) \middle| m\right) - \Pi\left(\frac{I_3}{I_3 - I_1} \middle| \operatorname{am}(\psi_0 | m) \middle| m\right) \right) \right], $$ $$ \cos \Phi_2 = \frac{\sqrt{I_3}}{G} \operatorname{dn}(\omega t + \psi_0 | m), \quad \tan \Phi_3 = -\sqrt{I_3 \omega} \frac{\operatorname{cn}(\omega t + \psi_0 | m)}{\operatorname{sn}(\omega t + \psi_0 | m)}. $$

- rotating extremals:

$$ \sin \alpha(t) = -\operatorname{cn}(\bar{\omega} t + \psi_0 | m), \quad \psi_0 = \operatorname{cn}^{-1}(-\sin \alpha(0) | m), \quad m^2 = \frac{I_1 - I_3}{G^2 - I_3}, \quad \bar{\omega} = \frac{\omega}{\sqrt{m}}, $$

$$ \Phi_1(t) - \Phi_1(0) = \frac{G}{I_3} \left[ t - \frac{1}{\bar{\omega}} \left( \Pi\left(\frac{I_3}{I_3 - G^2} \middle| \operatorname{am}(\bar{\omega} t + \psi_0 | m) \middle| m\right) - \Pi\left(\frac{I_3}{I_3 - G^2} \middle| \operatorname{am}(\psi_0 | m) \middle| m\right) \right) \right], $$ $$ \cos \Phi_2 = -\frac{\sqrt{I_3}}{G} \operatorname{cn}(\bar{\omega} t + \psi_0 | m), \quad \tan \Phi_3 = \sqrt{I_3 \bar{\omega}} \frac{\operatorname{dn}(\bar{\omega} t + \psi_0 | m)}{\operatorname{sn}(\bar{\omega} t + \psi_0 | m)}. $$

### 4.4. Integrability

We can write the Hamiltonian of the problem in the form $H = H_0 + (1 - k^2)H'$, where $H_0$ is the Hamiltonian of the Grushin metric and

$$ H' = \frac{1}{4k^2} (p_{\varphi} \sin \theta + p_{\theta} \cot \varphi \cos \theta)^2. $$

Since $G^2 = H_1^2 + H_2^2 + H_3^2$ is the Casimir function on $SO(3)$, we get the following
**Proposition 4.6.** *Denote by*

$$F = p_{\varphi}^{2} + \frac{p_{\theta}^{2}}{\sin^{2} \varphi}$$

the Hamiltonian of the round metric on $S^2$ written in spherical coordinates. Then $\{F, H_0\} = \{F, H'\} = 0$ and hence $\{F, H\} = 0$ for every $k \neq 0$.

A direct integration of the geodesics equation of the metric $g$ can be done using the following method ([4], [5]). The existence of the first integral quadratic on fibers allows to find normal coordinates $(u, v)$ such that the metric $g$ takes the Liouville normal form

$$g(u, v) = (f(u) + h(v))(du^2 + dv^2),$$

whose integration is standard. To this end one needs first to write the first integral $F$ in the isothermal coordinates $x, y$ of the metric $g$:

$$F(x, y) = b_1(x, y)p_x^2 + 2b_2(x, y)p_x p_y + b_3(x, y)p_y^2.$$

Then the coordinates $(u, v)$ can be found by solving the equation $\Phi'(w) = \sqrt{R(\Phi(w))}$ where $w = u + iv$, $\Phi: w \to z = x + iy$ and

$$R(z) = (b_1(x, y) - b_3(x, y)) + 2i b_2(x, y).$$

# 5. NUMERICAL ANALYSIS

An important tool in the analysis of the conjugate and cut loci is the Hampath code [12]. A series of commented numerical results is presented next. In all tests below we took $A = 1.5$, $B = 2$ and $C = 2.8$.

## 5.1. Serret-Andoyer Riemannian metric

In Fig. 3 we show the results of the numerical calculation of the conjugate loci (red curves) for the Serret-Andoyer problem using the Hampath code. The thick dashed line represents the separatrix. The cut locus is formed by the horizontal ray starting at the extremity of the horizontal cusp with $x \ge \theta_*(\varphi_0)$. It is formed by the self-intersections of the wave front, which is represented by the isocost curves (green dotted curves).

## 5.2. The spin case

Next the conjugate and cut loci are computed for the fixed initial conditions: $\varphi(0) = \pi/2$, $\theta(0) = 0$, and are represented via the deformation of the parameter $k$ starting from $k = 1$ using the method described in [8]. There are two different cases to be analyzed: $k > 1$ and $k < 1$. Similar computations for the cut locus of this problem can be found in [9]. Starting from the axis of symmetry, the Hamiltonian reduces to $H(\theta(0), \varphi(0), p_\theta(0), p_\varphi(0)) = p_\varphi^2(0)/4$, and restricting the extremals to $H = 1$, we can parameterize the geodesics by $p_\varphi(0) = \pm 2$, $p_\theta(0) \in \mathbb{R}$. By symmetry we can fix $p_\varphi(0) = -2$ and consider $p_\theta(0) \ge 0$. For any $k$, the conjugate locus has a contact of order two at the initial point, as $p_\theta(0) \to \infty$.

We study the deformation of the conjugate locus for $k \ge 1$ in Figs. 4–6. The key point is: when $k > 1$, $\theta$ is not monotonous for all the trajectories. This is true even for small $k$, like $k = 1.01$, taking $p_\theta(0) = 0.1$ and $t_f > 14$.

### 5.2.1. $k \ge 1$

We denote $t_1(p_\theta, k)$ the first conjugate time and $q_1(p_\theta, k) = (\theta, \varphi)|_{t=t_1(p_\theta, k)}$ the associated conjugate point. In Fig. 4, we represent the map $k \in [1, 1.5] \mapsto q_1(k)$ for $p_\theta$ fixed to $10^{-4}$. The value 1.5 is heuristically chosen to simplify the analysis. We can notice that $\theta(t_1(k))$ only takes approximately the values 0 and $\pi$ and so it is on the same meridian as the initial point. It switches three times at $1 < k_1 < k_2 < k_3 < 1.5$, with $k_2 - k_1 \neq k_3 - k_2$. We then restrict the study of the conjugate locus to $k \le k_3$ to simplify. We can see in Fig. 5, three subplots
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Conjugate locus of Serret-Andoyer metric.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. The first conjugate point with respect to $k$, for $p_{\theta}(0)$ fixed to $10^{-4}$. In red is plotted $\theta(t_1(p_{\theta}, k))$ while we have in blue $\varphi(t_1(p_{\theta}, k))$. The $\theta$-variable takes the values 0 and $\pi$. The values $k_1, k_2, k_3$ are approximately and respectively 1.061, 1.250, 1.429.</figcaption></figure>

which represent the deformation of one branch ($p_{\varphi}(0) = -2$ and $p_{\theta}(0) \ge 0$) of the conjugate locus resp. for $k$ in $[1, k_1]$, $[k_1, k_2]$ and $[k_2, k_3]$. For any $k \in [1, k_3]$, the branch is located in the half-plane $\theta \ge 0$. If we denote $k_1 < \bar{k} < k_2$, the parameter value such that $\varphi(t_1(\bar{k})) = \pi/2$, then the branch form a loop for $\bar{k} \le k \le k_3$.

The deformation of the conjugate locus can be explained analysing the behaviors of the trajectories. We describe four types of trajectories in $(\theta, \varphi)$-coordinates (see Fig. 6), limiting the study to $k \le k_3$ to simplify and $p_{\theta}(0) \ge 0$ by symmetry. These trajectories clarify the evolution of the conjugate locus.

i) The first type occuring for any $k$ such that $1 \le k \le k_3$, is represented in the top left subplot of Fig. 6. Its characteristic is that the $\theta$-variable is monotonous non-decreasing on $[0, t_1]$.

The three others trajectories do not have a monotonous $\theta$-variable on $[0, t_1]$. We denote $\bar{t}$ the first time when the trajectory leaves the domain $0 \le \theta \le \pi$.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. The deformation of one branch ($p_{\varphi}(0) = -2$ and $p_{\theta}(0) \ge 0$) of the conjugate locus with respect to the parameter $k \in [1, k_3]$. (top) $k = 1.0, 1.05$. (left) $k = 1.1, 1.2$. (right) $k = 1.3, 1.4$.</figcaption></figure>

ii) The second type (top right) existing for $k_1 \le k \le k_3$ has no self-intersection on $[0, \bar{t}]$ and is such that $\theta(\bar{t}) = 0$.

The last types of extremals have a self-intersection in the state-space in $[0, \bar{t}]$.

iii) The third kind of trajectories (bottom left) is such that $\theta(\bar{t}) = 0$ and occurs for $\bar{k} \le k \le k_3$.

iv) The last one (bottom right) exists only for $k_2 \le k \le k_3$ and has $\theta(\bar{t}) = \pi$.

### 5.2.2. $k \le 1$

The deformation of the conjugate locus in the case $k < 1$ is easier to analyze. We give on Fig. 7 the conjugate locus for $k \in \{0.8, 0.5, 0.2, 0.1\}$ with 15 chosen trajectories. The key point is the non-monotony of the $\theta$-variable for $k < 1$.

### 5.2.3. Deformation of the conjugate locus on the sphere

The deformation of the conjugate locus on the sphere is given Fig. 8. Only the half: $p_{\varphi}(0) = -2, p_{\theta}(0) \in \mathbf{R}$ is plotted to clarify the figures. The deformation is clear: the cusp moves along the meridian with respect to the parameter $k$. It does not cross the equator for $k < 1$ while for $k > 1$ it first crosses the North pole ($k = k_1$), then the equator ($k = \bar{k}$). For $k \ge \bar{k}$, the conjugate locus has self-intersections. Then, it crosses poles again for $k = k_2$ and $k_3$. This is repeated for greater values of $k$ making the loops smaller and smaller.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. The four types of trajectories which clarify the evolution of the conjugate locus.</figcaption></figure>

## 5.3. Riemannian metric on SO(3)

A preliminary computation of conjugate points is shown in Fig. 9. We consider two test trajectories defined by the initial conditions:

case 1: $\Phi_1 = 0, \Phi_2 = 2.8506, \Phi_3 = -0.32175, p_1 = p_2 = 0, p_3 = 0.3155946$;

case 2: $\Phi_1 = 0, \Phi_2 = 2.327, \Phi_3 = -0.32175, p_1 = 0.55, p_2 = 0, p_3 = 0.3155946$;

In Fig. 9 we show the evolution of the Euler angles along these trajectories. The cross indicates the starting point, while the thick red points mark the conjugate points.

By inverting the Serret-Andoyer transformation, one can show that both trajectories considered above project on the same curve defined by

$$p_x = 1.1, p_y(0) = 0.3155946, x_0 = 0, y_0 = \arccos \sqrt{0.1},$$

which lies on the energy level $h = 1$. The first conjugate time of this trajectory is $t_1^* \approx 5.4177$, but the comparison of the conjugate times $t$ and $\tau$ is not really meaningful, since the Serret-Andoyer transformation mixes the state and the co-state variables.
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Conjugate locus with 15 trajectories for $k = 0.8, 0.5, 0.2, 0.1$ from top left-hand to bottom right-hand.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8. Half of the conjugate locus on the sphere. (left) For $k = 1.0$ in magenta and $k = 0.8, 1.15$ in red. (right) For $k = 1.0$ in magenta and $k = 1.18$ in red.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Trajectories of the Euler-Poinsot problem in $SO(3)$. The first conjugate times are: $\tau = 2.42185846$ (case 1), $\tau = 2.4235292$ (case 2).</figcaption></figure>

REFERENCES

[1] A.A. Agrachev, U. Boscain, M. Sigalotti, A Gauss-Bonnet-like Formula on Two-Dimensional Almost-Riemannian Manifolds. Discrete and Continuous Dynamical Systems-A, vol. 20 (2008), 801–822.

[2] V. I. Arnold, *Mathematical Methods of Classical Mechanics*, vol. 60 of Graduate Texts in Mathematics. Translated from the Russian by K. Vogtmann and A. Weinstein, 2^nd^ edition, Springer-Verlag, New York, 1989, xvi+508.

[3] L. Bates, F. Fassò, *The conjugate locus for the Euler top. I. The axisymmetric case*, Int. Math. Forum, **2** (2007), no. 41–44, 2109–2139.

[4] G. D. Birkhoff, *Dynamical Systems*, AMS Colloquium Publications, vol. IX, 1927.

[5] A. V. Bolsinov, A. T. Fomenko, *Integrable Hamiltonian Systems. Geometry, Topology, Classification.*, Translated from the 1999 Russian original. Chapman & Hall/CRC, Boca Raton, FL, 2004, 730 pp.

[6] B. Bonnard, J.-B. Caillau, R. Sinclair and M. Tanaka, *Conjugate and cut loci of a two-sphere of revolution with application to optimal control*, Ann. Inst. H. Poincaré Anal. Non Linéaire, **26** (2009), no. 4, 1081–1098.

[7] B. Bonnard, J.-B. Caillau, G. Janin, *Conjugate-cut loci and injectivity domains on two-spheres of revolution*, ESAIM COCV, **19** (2013), no. 2, 533–554.

[8] B. Bonnard, O. Cots, N. Shcherbakova and D. Sugny, *The energy minimization problem for two-level dissipative quantum systems*, J. Math. Phys., **51** (2010), no. 9 092705, 44.

[9] U. Boscain, T. Chambrion, G. Charlot, *Nonisotropic 3-level Quantum Systems: Complete Solutions for Minimum Time and Minimum Energy*. Discrete and Continuous Dynamical Systems-B, **5** (2005), no. 4, 957–990.

[10] U. Boscain, G. Charlot, J.-P. Gauthier, S. Guérin, H.-R. Jauslin, *Optimal Control in laser-induced population transfer for two and three-level quantum systems*. J. Math. Phys., **43** (2002), 2107–2132.

[11] U. Boscain and F. Rossi, *Invariant Carnot-Caratheodory metrics on $S^3$, $SO(3)$, $SL(2)$, and lens spaces*. SIAM J. Control Optim., **47** (2008), no. 4, 1851–1878.

[12] J.-B. Caillau, O. Cots and J. Gergaud, *Differential continuation for regular optimal control problems*, Optimization Methods and Software, **27** (2011), no. 2, 177–196.

[13] D. D'Alessandro, *Introduction to quantum control and dynamics*, Chapman & Hall/CRC applied and nonlinear sciences series, 2008

[14] H. T. Davis, *Introduction to nonlinear differential and integral equations*. Dover Publications, Inc., New York, 1962, xv+566pp.
[15] P. Gurfil, A. Elipe, W. Tangren and M. Efroimsky, *The Serret-Andoyer formalism in rigid-body dynamics. I. Symmetries and perturbations*, Regul. Chaotic Dyn., **12** (2007), no. 4, 389–425.

[16] J. Itoh and K. Kiyohara, *The cut loci and the conjugate loci on ellipsoids*, Manuscripta Math., **114** (2004), no. 2, 247–264.

[17] V. Jurdjevic, “Geometric Control Theory”, vol. 52 of *Cambridge Studies in Advanced Mathematics*, Cambridge University Press, Cambridge, 1997.

[18] N. Khaneja, R. Brockett and S. J. Glaser, *Sub-Riemannian geometry and time optimal control of three spin systems: Quantum gates and coherence transfer*, Physical Review A, vol. 65, 032301.

[19] M. Lara, S. Ferrer, *Closed form Integration of the Hitzl-Breakwell problem in action-angle variables*. IAA-AAS-DyCoSS1-01-02 (AAS 12-302), 27–39.

[20] D. Lawden, “Elliptic Functions and Applications”, vol. 80 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1989, xiv+334.

[21] M. H. Levitt, “Spin dynamics, basis of Nuclear Magnetic Resonance”, John Wiley and sons, 2007 (2nd edition).

[22] H. Poincaré, *Sur les lignes géodésiques des surfaces convexes*. Trans. Amer. Math. Soc., **6** (1905), no. 3, 237–274.

[23] L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze, E.F. Mishchenko, The mathematical theory of optimal processes. Inter-science Publishers John Wiley & Sons, Inc. New York-London 1962.

[24] K. Shiohama, T. Shioya and M. Tanaka, “The Geometry of Total Curvature on Complete Open Surfaces”, vol. 159 of *Cambridge Tracts in Mathematics*, Cambridge University Press, Cambridge, 2003, x+284.

[25] R. Sinclair and M. Tanaka, *The cut locus of a two-sphere of revolution and Toponogov's comparison theorem*. Tohoku Math. J. (2), **59** (2007), no. 3, 379–399.

[26] A.M. Vershik and V.Ya. Gershkovich, *Nonholonomic Dynamical Systems, Geometry of Distributions and Variational Problems*. in Dynamical Systems VII, Encyclopedia of Mathematical Sciences vol. 16, Springer Verlag, 10–81.

[27] H. Yuan, *Geometry, optimal control and quantum computing*, Phd Thesis, Harvard, (2006).

[28] H. Yuan, R. Zeier, N. Khaneja. *Elliptic functions and efficient control of Ising spin chains with unequal coupling*, Physical Review A77, 032340, 2008.