\section{Introduction}

The purpose of this paper is to study the geometry and the dynamics of the free rigid body on the Lie algebra $\mathfrak{so}(5)$. The general free rigid body on $\mathfrak{so}(n)$ has been studied as a completely integrable system in the classical works \cite{Mishchenko70}, \cite{Manakov76}, \cite{Ratiu80}. The general case and more extensively the even case $n=2m$ have been analyzed by Feh\'{e}r and Marshall \cite{FeMa03}, Spiegler
\cite{Spiegler04}, where a method for determining a certain class of equilibria is given. They also study the stability of these equilibria using the Energy-Casimir method. As mentioned in \cite{Spiegler04}, the odd case $n=2m+1$ is significantly different and only some indications on stability results have been given.

In \cite{noi} it has been proved that this class of equilibria comes as the intersection of coordinate type Cartan subalgebras with regular adjoint orbits. Consequently these equilibria are unions of Weyl orbits and they correspond to long-short axis type of equilibria known from the dynamics on $\mathfrak{so}(3)$.

The stability of these equilibria, for the case of $\mathfrak{so}(4)$, has been studied in \cite{noi}, using Williamson normal form \cite{BoFo04}.
For the case of $\mathfrak{so}(4)$ there have been discovered two new continuous families of equilibria on every regular adjoint orbit. Also, in
\cite{noi} it has been shown that these are nonlinearly stable as a family, that is, if a solution of the $\mathfrak{so}(4)$-free rigid body equation starts near an equilibrium on such a curve, at any later time it will stay close to this curve but in the direction of the curve itself it may drift.

First of all, in this paper we will give a list of integrals of motion for the general case of the $\mathfrak{so}(n)$ free rigid body, which proves to be a set of generating functions for the Mishchenko's quadratic integrals of motion and which will play a crucial role in the approach of the stability problem using energy methods.

For the case of $\mathfrak{so}(5)$ rigid body we found fifteen families of coordinate type Cartan subalgebras, which intersected with a regular adjoint orbit give 120=15$\times$8 equilibria corresponding to long-short axis type of equilibria (8 being the cardinal of the Weyl group of $\mathfrak{so}(5)$). We also found ten continuous families of equilibria on every regular adjoint orbit. As in the case of $\mathfrak{so}(4)$, these two types of equilibria give all the equilibria of $\mathfrak{so}(5)$ free rigid body.

The stability problem is studied using the linearization method and Arnold's method, which is equivalent to Energy-Casimir method \cite{BiPu07}. In order to do this, we explicitly compute the integrals of motion using Mishchenko's and Manakov's methods and we point out the need for the additional constants of motion introduced in Section 2 of this paper. Extended indications on how Arnold's method applies for each studied equilibria are given.

\section{A new family of integrals of motion for the free rigid body on $\mathfrak{so}(n)$}

The equations of the rigid body on $\mathfrak{so}(n)$ are given by
\begin{equation}\label{rigideq}
\dot M=[M,\Omega], \end{equation} where $\Omega\in \mathfrak{so}(n)$,
$M=\Omega J+J\Omega\in \mathfrak{so}(n)$ with $J=\operatorname{diag}(\lambda_i)$, a real constant diagonal matrix satisfying $ \lambda _i + \lambda _j \geq 0 $, for all $i, j=1, \ldots, n$, $i \neq j $ (see, for example, \cite{Ratiu80}).
Note that $M=[m_{ij}]$ and $\Omega=[ \omega_{ij}]$ determine each other if and only if $\lambda_i + \lambda _j >0$ since $m_{ij} = ( \lambda_i+ \lambda_j) \omega_{ij}$ which physically means that the rigid body is not concentrated on a lower dimensional subspace of $\mathbb{R}^n$.

It is well known and easy to verify that equations \eqref{rigideq} are Hamiltonian relative to the minus Lie-Poisson bracket
\begin{equation}
\label{LP_n}
\{F,G\}(M) : = \frac{1}{2}\operatorname{Trace}(M [\nabla F(M), \nabla G(M)]),
\end{equation}
and the Hamiltonian function
\begin{equation}
\label{ham_n}
H(M) := - \frac{1}{4} \operatorname{Trace}(M \Omega).
\end{equation}
Here $F,G, H \in C ^{\infty}(\mathfrak{so}(n))$ and the gradient is taken relative to the Ad-invariant inner product
\begin{equation}
\label{ip_n}
\left\langle X, Y \right\rangle : = - \frac{1}{2}\operatorname{Trace}(XY),
\qquad X, Y \in \mathfrak{so}(n)
\end{equation}
which identifies $(\mathfrak{so}(n))^*$ with $\mathfrak{so}(n)$. This means that $\dot{F} = \{F, H\}$ for all $F \in C^{\infty}(\mathfrak{so}(n))$, where $\{\cdot , \cdot \}$ is given by \eqref{LP_n} and $H $ by \eqref{ham_n}, if and only if \eqref{rigideq} holds.
Note that the linear isomorphism $X \in \mathfrak{so}(n) \mapsto XJ + JX
\in \mathfrak{so}(n)$ is self-adjoint relative to the inner product \eqref{ip_n} and thus $\nabla H(M) =  \Omega$.

In what follows, we assume that $\lambda_i$ are all distinct. As stated in \cite{Mishchenko70}, the system \eqref{rigideq} admits a sequence $m_r$ ($r=1,2,3,\dots$) of integrals of motion, all of them depending quadratically on the angular momentum $M$. It is easy to compute an explicit form for these integrals \cite{MoPi} as follows:
$$m_r(M)=\sum\limits_{\stackrel{i,k=1}{i<k}}^n\frac{\lambda_i^r-\lambda_k^r}{\lambda_i^2-\lambda_k^2}m_{ik}^2.$$
Next we will find an additional family of $n$ integrals of motion for the system \eqref{rigideq}, which have simple and elegant expressions and which prove to generate Mishchenko's integrals of motion. The complete result is contained in the following
\begin{thm}
\label{generators}
The functions
$$F_i(M)=\sum\limits_{\stackrel{k=1}{k\not= i}}^{n}\frac{1}{\lambda_i^2-\lambda_k^2}m_{ik}^2,~~~i=1,2,\dots,n$$
are integrals of motion for the system \eqref{rigideq}; moreover, the following identities hold for each $r$:
$$m_r(M)=\sum\limits_{i=1}^n\lambda_i^rF_i(M).$$
\end{thm}

\noindent {\bf Proof.} It is easy to see that
\begin{align*}
m_r(M)&=\sum\limits_{\stackrel{i,k=1}{i<k}}^n\frac{\lambda_i^r}{\lambda_i^2-\lambda_k^2}m_{ik}^2-\sum\limits_{\stackrel{i,k=1}{i<k}}^n\frac{\lambda_k^r}{\lambda_i^2-\lambda_k^2}m_{ik}^2\\
&=\sum\limits_{\stackrel{i,k=1}{i<k}}^n\frac{\lambda_i^r}{\lambda_i^2-\lambda_k^2}m_{ik}^2+\sum\limits_{\stackrel{i,k=1}{k<i}}^n\frac{\lambda_i^r}{\lambda_i^2-\lambda_k^2}m_{ki}^2\\
&=\sum_{i=1}^n\lambda_i^r\left(\sum\limits_{\stackrel{k=1}{k>i}}^n\frac{1}{\lambda_i^2-\lambda_k^2}m_{ik}^2+\sum\limits_{\stackrel{k=1}{k<i}}^n\frac{1}{\lambda_i^2-\lambda_k^2}m_{ik}^2\right)\\
&=\sum_{i=1}^n\lambda_i^r\cdot\sum\limits_{\stackrel{k=1}{k\not= i}}^n\frac{1}{\lambda_i^2-\lambda_k^2}m_{ik}^2\\
&=\sum_{i=1}^n\lambda_i^r F_i(M).
\end{align*}
Now, writing the above identities for $r=1,2,\dots,n$ we obtain a Vandermonde-type linear system, which solved leads to the conclusion that the functions $F_i$ are linear combinations of Mishchenko's integrals of motion $m_r$ ($r=1,2,\dots,n$) and thus are integrals of motion for our system \eqref{rigideq}. \rule{0.5em}{0.5em}

\section{The free rigid body on the Lie algebra $\mathfrak{so}(5)$}

We begin by making some considerations on $\mathfrak{so}(5)$, which is the Lie algebra of the subgroup
$\operatorname{SO}(5)=\{A\in \mathfrak{gl}(5, \mathbb{R})
\mid A^tA=I_5,\det (A)=1\}$ of the special linear Lie group
$\operatorname{SL}(5,\mathbb{R})$.

We choose as basis of $\mathfrak{so}(5)$ the matrices
$$E_1=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&-1&0&0\\
0&1&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0
\end{array}\right];~E_2=\left[\begin{array}{ccccc}
0&0&1&0&0\\
0&0&0&0&0\\
-1&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\end{array}\right];$$
$$E_3=\left[\begin{array}{ccccc}
0&-1&0&0&0\\
1&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\end{array}\right];
E_4=\left[\begin{array}{ccccc}
0&0&0&1&0\\
0&0&0&0&0\\
0&0&0&0&0\\
-1&0&0&0&0\\
0&0&0&0&0\end{array}\right];$$
$$E_5=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&0&1&0\\
0&0&0&0&0\\
0&-1&0&0&0\\
0&0&0&0&0\end{array}\right];E_6=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&1&0\\
0&0&-1&0&0\\
0&0&0&0&0\end{array}\right];$$
$$
E_7=\left[\begin{array}{ccccc}
0&0&0&0&1\\
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
-1&0&0&0&0\end{array}\right];E_8=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&0&0&1\\
0&0&0&0&0\\
0&0&0&0&0\\
0&-1&0&0&0\end{array}\right];$$
$$E_9=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&1\\
0&0&0&0&0\\
0&0&-1&0&0\end{array}\right];
E_{10}=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&1\\
0&0&0&-1&0\end{array}\right]$$
and hence we represent $\mathfrak{so}(5)$ as
\begin{equation}
\label{so4_representation}
\mathfrak{so}(5)=\left\{\, M=\left|\, \left[\begin{array}{ccccc}
0&-x_3&x_2&y_1&z_1\\
x_3&0&-x_1&y_2&z_2\\
-x_2&x_1&0&y_3&z_3\\
-y_1&-y_2&-y_3&0&z_4\\
-z_1&-z_2&-z_3&-z_4&0
\end{array}\right] \, \right|\, x_1,x_2,x_3,y_1,y_2,y_3,z_1,z_2,z_3,z_4\in \mathbb{R}\right\}.
\end{equation}
Relative to this basis, the Lie algebra structure of $\mathfrak{so}(5)$ is given by the following table
\medskip
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
\hline
$[\cdot,\cdot]$&$E_1$&$E_2$&$E_3$&$E_4$&$E_5$&$E_6$&$E_7$&$E_8$&$E_9$&$E_{10}$\\
\hline
$E_1$&0&$E_3$&$-E_2$&0&$E_6$&$-E_5$&0&$E_9$&$-E_8$&0\\
\hline
$E_2$&$-E_3$&0&$E_1$&$-E_6$&0&$E_4$&$-E_9$&0&$E_7$&0\\
\hline
$E_3$&$E_2$&$-E_1$&0&$E_5$&$-E_4$&0&$E_8$&$-E_7$&0&0\\
\hline
$E_4$&0&$E_6$&$-E_5$&0&$E_3$&$-E_2$&$-E_{10}$&0&0&$E_7$\\
\hline
$E_5$&$-E_6$&0&$E_4$&$-E_3$&0&$E_1$&0&$-E_{10}$&0&$E_8$\\
\hline
$E_6$&$E_5$&$-E_4$&0&$E_2$&$-E_1$&0&0&0&$-E_{10}$&$E_9$\\
\hline
$E_7$&0&$E_9$&$-E_8$&$E_{10}$&0&0&0&$E_3$&$-E_2$&$-E_4$\\
\hline
$E_8$&$-E_9$&0&$E_7$&0&$E_{10}$&0&$-E_3$&0&$E_1$&$-E_5$\\
\hline
$E_9$&$E_8$&$-E_7$&0&0&0&$E_{10}$&$E_2$&$-E_1$&0&$-E_6$\\
\hline
$E_{10}$&0&0&0&$E_7$&$-E_8$&$-E_9$&$E_4$&$E_5$&$E_6$&0\\
\hline
\end{tabular}
\end{center}
\medskip
In the basis $\{E_1, \ldots, E_{10}\}$, the matrix of the Lie-Poisson
structure \eqref{LP_n} is
\begin{equation}
\label{LP}
\Gamma_-=\left[\begin{array}{cccccccccc}
0&-x_3&x_2&0&-y_3&y_2&0&-z_3&z_2&0\\
x_3&0&-x_1&y_3&0&-y_1&z_3&0&-z_1&0\\
-x_2&x_1&0&-y_2&y_1&0&-z_2&z_1&0&0\\
0&-y_3&y_2&0&-x_3&x_2&z_4&0&0&-z_1\\
y_3&0&-y_1&x_3&0&-x_1&0&z_4&0&-z_2\\
-y_2&y_1&0&-x_2&x_1&0&0&0&z_4&-z_3\\
0&-z_3&z_2&-z_4&0&0&0&-x_3&x_2&y_1\\
z_3&0&-z_1&0&-z_4&0&x_3&0&-x_1&y_2\\
-z_2&z_1&0&0&0&-z_4&-x_2&x_1&0&y_3\\
0&0&0&z_1&z_2&z_3&-y_1&-y_2&-y_3&0
\end{array}\right].
\end{equation}
Since $\hbox{rank}\,\mathfrak{so}(5) = 2$, there are two functionally independent Casimir functions which are given respectively by
\[
C_1(M) := -\frac{1}{4}\operatorname{Trace}(M^2)
= \frac{1}{2}\left(\sum\limits\limits_{i=1}^3 x_i^2+\sum\limits\limits_{i=1}^3 y_i^2+\sum\limits\limits_{i=1}^4 z_i^2 \right)
\]
and
\[
C_2(M) := \frac{1}{8}\operatorname{Trace}(M^4)=
\]
$$=\frac{1}{8}\left[(x_2^2+x_3^2+y_1^2+z_1^2)^2+(x_1^2+x_3^2+y_2^2+z_2^2)^2+(x_1^2+x_2^2+y_3^2+z_3^2)^2+\right.$$
$$+(y_1^2+y_2^2+y_3^2+z_4^2)^2+(z_1^2+z_2^2+z_3^2+z_4^2)^2+2(y_1z_1+y_2z_2+y_3z_3)^2+$$
$$+2(x_2y_3-x_3y_2-z_1z_4)^2+2(x_3y_1-x_1y_3-z_2z_4)^2+2(x_1y_2-x_2y_1-z_3z_4)^2+$$
$$+2(x_3z_2-x_2z_3-y_1z_4)^2+2(x_1z_3-x_3z_1-y_2z_4)^2+2(x_2z_1-x_1z_2-y_3z_4)^2+$$
$$\left. +2(x_1x_2-y_1y_2-z_1z_2)^2+2(x_1x_3-y_1y_3-z_1z_3)^2+2(x_2x_3-y_2y_3-z_2z_3)^2 \right].$$
Thus the generic adjoint orbits are the level sets
$$
\operatorname{Orb}_{c_1c_2}(M)=(C_1\times C_2)^{-1}(c_1,c_2),
\qquad (c_1,c_2) \in \mathbb{R}^2.
$$

In all that follows we will denote by $\operatorname{Orb}_{c_1;c_2}$ the regular adjoint orbit $\operatorname{Orb}_{c_1c_2}$,
where $c_1>0,c_2>0$ and $2c_2>c_1^2>c_2$.

Using the Lie bracket table in the chosen basis given above, it is immediately seen that the coordinate type Cartan subalgebras of $\mathfrak{so}(5)$ are
$\mathfrak{t}_k,~~1\leq k\leq 15$, where
$$\mathfrak{t}_1:=\operatorname{span}(E_3,E_6)=\left\{\left. M_{a,b}^1:=\left[\begin{array}{ccccc}
0&a&0&0&0\\
-a&0&0&0&0\\
0&0&0&b&0\\
0&0&-b&0&0\\
0&0&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_2:=\operatorname{span}(E_6,E_8)=\left\{\left. M_{a,b}^2:=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&0&0&a\\
0&0&0&-b&0\\
0&0&b&0&0\\
0&-a&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_3:=\operatorname{span}(E_6,E_7)=\left\{\left. M_{a,b}^3:=\left[\begin{array}{ccccc}
0&0&0&0&a\\
0&0&0&0&0\\
0&0&0&-b&0\\
0&0&b&0&0\\
-a&0&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_4:=\operatorname{span}(E_5,E_7)=\left\{\left. M_{a,b}^4:=\left[\begin{array}{ccccc}
0&0&0&0&a\\
0&0&0&b&0\\
0&0&0&0&0\\
0&-b&0&0&0\\
-a&0&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_5:=\operatorname{span}(E_1,E_7)=\left\{\left. M_{a,b}^:=\left[\begin{array}{ccccc}
0&0&0&0&a\\
0&0&-b&0&0\\
0&b&0&0&0\\
0&0&0&0&0\\
-a&0&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_6:=\operatorname{span}(E_2,E_5)=\left\{\left. M_{a,b}^6:=\left[\begin{array}{ccccc}
0&0&a&0&0\\
0&0&0&-b&0\\
-a&0&0&0&0\\
0&b&0&0&0\\
0&0&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_7:=\operatorname{span}(E_5,E_9)=\left\{\left. M_{a,b}^7:=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&0&a&0\\
0&0&0&0&-b\\
0&-a&0&0&0\\
0&0&b&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_8:=\operatorname{span}(E_1,E_{10})=\left\{\left. M_{a,b}^8:=\left[\begin{array}{ccccc}
0&0&0&0&0\\
0&0&a&0&0\\
0&-a&0&0&0\\
0&0&0&0&-b\\
0&0&0&b&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_9:=\operatorname{span}(E_1,E_4)=\left\{\left. M_{a,b}^9:=\left[\begin{array}{ccccc}
0&0&0&a&0\\
0&0&-b&0&0\\
0&b&0&0&0\\
-a&0&0&0&0\\
0&0&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_{10}:=\operatorname{span}(E_2,E_8)=\left\{\left. M_{a,b}^{10}:=\left[\begin{array}{ccccc}
0&0&a&0&0\\
0&0&0&0&-b\\
-a&0&0&0&0\\
0&0&0&0&0\\
0&b&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_{11}:=\operatorname{span}(E_4,E_8)=\left\{\left. M_{a,b}^{11}:=\left[\begin{array}{ccccc}
0&0&0&a&0\\
0&0&0&0&-b\\
0&0&0&0&0\\
-a&0&0&0&0\\
0&b&0&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_{12}:=\operatorname{span}(E_3,E_{10})=\left\{\left. M_{a,b}^{12}:=\left[\begin{array}{ccccc}
0&a&0&0&0\\
-a&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&b\\
0&0&0&-b&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_{13}:=\operatorname{span}(E_3,E_9)=\left\{\left. M_{a,b}^{13}:=\left[\begin{array}{ccccc}
0&a&0&0&0\\
-a&0&0&0&0\\
0&0&0&0&-b\\
0&0&0&0&0\\
0&0&b&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_{14}:=\operatorname{span}(E_4,E_9)=\left\{\left. M_{a,b}^{14}:=\left[\begin{array}{ccccc}
0&0&0&a&0\\
0&0&0&0&0\\
0&0&0&0&-b\\
-a&0&0&0&0\\
0&0&b&0&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\},$$
$$\mathfrak{t}_{15}:=\operatorname{span}(E_2,E_{10})=\left\{\left. M_{a,b}^{15}:=\left[\begin{array}{ccccc}
0&0&a&0&0\\
0&0&0&0&0\\
-a&0&0&0&0\\
0&0&0&0&b\\
0&0&0&-b&0\end{array}\right]\, \right| \,a,b\in \mathbb{R}
 \right\}.$$

The intersection of a regular adjoint orbit and a coordinate
type Cartan subalgebra has eight elements which represents a Weyl group orbit. Specifically, we have the following result:

\begin{thm} \label{Weil_group_orbit}
The following equalities hold:
$$\mathfrak{t}_k\cap \operatorname{Orb}_{c_1;c_2}=\left\{
M_{a,b}^k,M_{-a,-b}^k,M_{b,a}^k,M_{-b,-a}^k,M_{-a,b}^k,M_{a,-b}^k,M_{b,-a}^k,M_{-b,a}^k
\right\}$$
for all $k$, $1\leq k\leq 15$,
where
\begin{equation}
\label{values_a_b}
\left\{
\begin{aligned}
a&=\sqrt{c_1+\sqrt{2c_2-c_1^2}}\\
b&=\sqrt{c_1-\sqrt{2c_2-c_1^2}}.
\end{aligned}
\right.
\end{equation}
\end{thm}

\noindent {\bf Proof.}
Let $M_{\alpha,\beta}^1\in \mathfrak{t}_1\cap \operatorname{Orb}_{c_1;c_2}$.
Then $M_{\alpha,\beta}^1\in\mathfrak{t}_1$, $c_1=C_1(M_{\alpha,\beta}^1)=\frac{1}{2}\left(\alpha^2+\beta^2\right)$ and $c_2=C_2(M_{\alpha,\beta}^1)=\frac{1}{4}\left(\alpha^4+\beta^4\right)$.
We obtain the system
$$\left\{\begin{array}{l}
\alpha^2+\beta^2=2c_1\\
\alpha^4+\beta^4=4c_2\end{array}\right.$$
which leads immediately to the result. For $\mathfrak{t}_k\cap \operatorname{Orb}_{c_1;c_2}$, $2\leq k\leq 15$, we proceed in a similar manner.
\rule{0.5em}{0.5em}
\medskip

\section{Equilibria for the $\mathfrak{so}(5)$-rigid body}

We will work from now on with a generic $\mathfrak{so}(5)$-rigid body, that is, $\lambda_i+ \lambda_j>0$ for $i \neq j$ and all $\lambda_i $ are distinct.
The relationship between $\Omega=[\omega_{ij}]\in \mathfrak{so}(5)$ and the matrix $M\in \mathfrak{so}(5)$ in the representation \eqref{so4_representation} is hence given by
$$(\lambda_3+\lambda_2)\omega_{32}=x_1; \quad
(\lambda_1+\lambda_3)\omega_{13}=x_2; \quad
(\lambda_2+\lambda_1)\omega_{21}=x_3;$$
$$(\lambda_1+\lambda_4)\omega_{14}=y_1; \quad
(\lambda_2+\lambda_4)\omega_{24}=y_2; \quad
(\lambda_3+\lambda_4)\omega_{34}=y_3;$$
$$(\lambda_1+\lambda_5)\omega_{15}=z_1; \quad
(\lambda_2+\lambda_5)\omega_{25}=z_2; \quad
(\lambda_3+\lambda_5)\omega_{35}=z_3; \quad
(\lambda_4+\lambda_5)\omega_{45}=z_4.$$
and thus the equations of motion \eqref{rigideq} are equivalent to the system
\begin{equation}\label{soM}
\left\{\begin{array}{l}
\vspace{.1cm}
\dot x_1=(\lambda_2-\lambda_3)\left[\frac{y_2y_3}{(\lambda_2+\lambda_4)(\lambda_3+\lambda_4)}-\frac{x_2x_3}{(\lambda_1+\lambda_2)(\lambda_1+\lambda_3)}+\frac{z_2z_3}{(\lambda_2+\lambda_5)(\lambda_3+\lambda_5)}\right]\\\\
\dot x_2=(\lambda_3-\lambda_1)\left[\frac{y_1y_3}{(\lambda_1+\lambda_4)(\lambda_3+\lambda_4)}-\frac{x_1x_3}{(\lambda_1+\lambda_2)(\lambda_2+\lambda_3)}+\frac{z_1z_3}{(\lambda_1+\lambda_5)(\lambda_3+\lambda_5)}\right]\\\\
\dot x_3=(\lambda_1-\lambda_2)\left[\frac{y_1y_2}{(\lambda_1+\lambda_4)(\lambda_2+\lambda_4)}-\frac{x_1x_2}{(\lambda_2+\lambda_3)(\lambda_1+\lambda_3)}+\frac{z_1z_2}{(\lambda_1+\lambda_5)(\lambda_2+\lambda_5)}\right]\\\\
\dot y_1=(\lambda_1-\lambda_4)\left[\frac{x_2y_3}{(\lambda_1+\lambda_3)(\lambda_3+\lambda_4)}-\frac{x_3y_2}{(\lambda_1+\lambda_2)(\lambda_2+\lambda_4)}-\frac{z_1z_4}{(\lambda_1+\lambda_5)(\lambda_4+\lambda_5)}\right]\\\\
\dot y_2=(\lambda_2-\lambda_4)\left[\frac{x_3y_1}{(\lambda_1+\lambda_2)(\lambda_1+\lambda_4)}-\frac{x_1y_3}{(\lambda_2+\lambda_3)(\lambda_3+\lambda_4)}-\frac{z_2z_4}{(\lambda_2+\lambda_5)(\lambda_4+\lambda_5)}\right]\\\\
\dot y_3=(\lambda_3-\lambda_4)\left[\frac{x_1y_2}{(\lambda_2+\lambda_3)(\lambda_2+\lambda_4)}-\frac{x_2y_1}{(\lambda_1+\lambda_3)(\lambda_1+\lambda_4)}-\frac{z_3z_4}{(\lambda_3+\lambda_5)(\lambda_4+\lambda_5)}\right]\\\\
\dot z_1=(\lambda_1-\lambda_5)\left[\frac{x_2z_3}{(\lambda_1+\lambda_3)(\lambda_3+\lambda_5)}-\frac{x_3z_2}{(\lambda_1+\lambda_2)(\lambda_2+\lambda_5)}+\frac{y_1z_4}{(\lambda_1+\lambda_4)(\lambda_4+\lambda_5)}\right]\\\\
\dot z_2=(\lambda_2-\lambda_5)\left[
\frac{x_3z_1}
{(\lambda_1+\lambda_2)(\lambda_1+\lambda_5)}-\frac{x_1z_3}{(\lambda_2+\lambda_3)(\lambda_3+\lambda_5)}+\frac{y_2z_4}{(\lambda_2+\lambda_4)(\lambda_4+\lambda_5)}\right]\\\\
\dot z_3=(\lambda_3-\lambda_5)\left[\frac{x_1z_2}{(\lambda_2+\lambda_3)(\lambda_2+\lambda_5)}-\frac{x_2z_1}{(\lambda_1+\lambda_3)(\lambda_1+\lambda_5)}+\frac{y_3z_4}{(\lambda_3+\lambda_4)(\lambda_4+\lambda_5)}\right]\\\\
\dot z_4=(\lambda_4-\lambda_5)\left[-\frac{y_1z_1}{(\lambda_1+\lambda_4)(\lambda_1+\lambda_5)}-\frac{y_2z_2}{(\lambda_2+\lambda_4)(\lambda_2+\lambda_5)}-\frac{y_3z_3}{(\lambda_3+\lambda_4)(\lambda_3+\lambda_5)}\right].
\end{array}\right.
\end{equation}
The Hamiltonian \eqref{ham_n} has in this case the expression
\begin{align*}
H(M)&=-\frac{1}{4}\hbox{Trace}(M\Omega)\\
&=\frac{1}{2}\left(\frac{1}{\lambda_2+\lambda_3}x_1^2+\frac{1}{\lambda_1+\lambda_3}x_2^2+\frac{1}{\lambda_1+\lambda_2}x_3^2+\frac{1}{\lambda_1+\lambda_4}y_1^2+\frac{1}{\lambda_2+\lambda_4}y_2^2+\frac{1}{\lambda_3+\lambda_4}y_3^2+\right.\\
&+\left.\frac{1}{\lambda_1+\lambda_5}z_1^2+\frac{1}{\lambda_2+\lambda_5}z_2^2+\frac{1}{\lambda_3+\lambda_5}z_3^2+\frac{1}{\lambda_4+\lambda_5}z_4^2\right).
\end{align*}
The Hamiltonian nature of system \eqref{soM} can be checked in this case directly, writing
$$(\dot{x}_1, \dot{x}_2, \dot{x}_3, \dot{y}_1, \dot{y}_2, \dot{y}_3,\dot{z}_1,\dot{z}_2,\dot{z}_3,\dot{z}_4)^{\sf T} = \Gamma_- (\nabla H)^{\sf T},$$ where the Poisson structure $\Gamma_- $ is given by \eqref{LP}.

\begin{thm} \label{echilibre}
If $\mathcal{E}$ denotes the set of the equilibrium points
of \eqref{soM}, then
${\cal E}=\left(\bigcup\limits_{k=1}^{15}\mathfrak{t}_k\right)\cup \left(\bigcup\limits_{l=1}^{10}\mathfrak{s}_l\right)$,
where $\mathfrak{s}_l,~1\leq l\leq 10$, are the three dimensional vector subspaces given by\\

\medskip

$\mbox{\fontsize{10}{10}\selectfont $
\mathfrak{s}_{1,2}:=\operatorname{span}_{\mathbb{R}}\left\{\left(\frac{1}{\lambda_1+\lambda_4}E_1 \pm \frac{1}{\lambda_2+\lambda_3}E_4\right),\,\left(\frac{1}{\lambda_2+\lambda_4}E_2\pm \frac{1}{\lambda_1+\lambda_3}E_5\right), \,\left(\frac{1}{\lambda_3+\lambda_4}E_3\pm \frac{1}{\lambda_1+\lambda_2}E_6\right) \right\};$}$

\medskip

$\mbox{\fontsize{10}{10}\selectfont $
\mathfrak{s}_{3,4}:=\operatorname{span}_{\mathbb{R}}\left\{\left(\frac{1}{\lambda_4+\lambda_5}E_1 \pm \frac{1}{\lambda_2+\lambda_3}E_{10}\right),\,\left(\frac{1}{\lambda_3+\lambda_5}E_5\pm \frac{1}{\lambda_2+\lambda_4}E_9\right), \,\left(\frac{1}{\lambda_2+\lambda_5}E_6\mp \frac{1}{\lambda_3+\lambda_4}E_8\right) \right\};$}$

\medskip

$\mbox{\fontsize{10}{10}\selectfont $
\mathfrak{s}_{5,6}:=\operatorname{span}_{\mathbb{R}}\left\{\left(\frac{1}{\lambda_4+\lambda_5}E_2 \pm \frac{1}{\lambda_1+\lambda_3}E_{10}\right),\,\left(\frac{1}{\lambda_1+\lambda_5}E_6\pm \frac{1}{\lambda_3+\lambda_4}E_7\right), \,\left(\frac{1}{\lambda_3+\lambda_5}E_4\mp \frac{1}{\lambda_1+\lambda_4}E_9\right) \right\};$}$

\medskip

$\mbox{\fontsize{10}{10}\selectfont $
\mathfrak{s}_{7,8}:=\operatorname{span}_{\mathbb{R}}\left\{\left(\frac{1}{\lambda_1+\lambda_5}E_1 \pm \frac{1}{\lambda_2+\lambda_3}E_7\right),\,\left(\frac{1}{\lambda_2+\lambda_5}E_2\pm \frac{1}{\lambda_1+\lambda_3}E_8\right), \,\left(\frac{1}{\lambda_3+\lambda_5}E_3\pm \frac{1}{\lambda_1+\lambda_2}E_9\right) \right\};$}$

\medskip

$\mbox{\fontsize{10}{10}\selectfont $
\mathfrak{s}_{9,10}:=\operatorname{span}_{\mathbb{R}}\left\{\left(\frac{1}{\lambda_4+\lambda_5}E_3 \pm \frac{1}{\lambda_1+\lambda_2}E_{10}\right),\,\left(\frac{1}{\lambda_2+\lambda_5}E_4\pm \frac{1}{\lambda_1+\lambda_4}E_8\right), \,\left(\frac{1}{\lambda_1+\lambda_5}E_5\mp \frac{1}{\lambda_2+\lambda_4}E_7\right) \right\}.$}$
\end{thm}

\noindent {\bf Proof.} The proof follows after a long, but straightforward computation. \rule{0.5em}{0.5em}

\section{Constants of motion and nonlinear stability}
In this section we study the nonlinear stability of the
equilibrium states ${\cal E}\cap \operatorname{Orb}_{c_1;c_2}$ for the dynamics \eqref{soM} on a generic adjoint orbit.

Using Mishchenko's method \cite{Mishchenko70}, \cite{Ratiu80}, \cite{MoPi}, we obtain the following additional independent constants of the motion for Eqs. \eqref{soM}:
\begin{align*}
K_1(M)&=-\frac{1}{4}\operatorname{Trace}\sum\limits_{p=0}^3 J^{p}MJ^{3-p}\Omega\\
&=\frac{1}{2}[(\lambda_2^2+\lambda_3^2)x_1^2+(\lambda_1^2+\lambda_3^2)x_2^2+(\lambda_1^2+\lambda_2^2)x_3^2+(\lambda_1^2+\lambda_4^2)y_1^2+(\lambda_2^2+\lambda_4^2)y_2^2+(\lambda_3^2+\lambda_4^2)y_3^2\\
&+(\lambda_1^2+\lambda_5^2)z_1^2+(\lambda_2^2+\lambda_5^2)z_2^2+(\lambda_3^2+\lambda_5^2)z_3^2+(\lambda_4^2+\lambda_5^2)z_4^2]
\end{align*}
and respectively
\begin{align*}
K_2(M)&=-\frac{1}{4}\operatorname{Trace}\sum\limits_{p=0}^5 J^{p}MJ^{5-p}\Omega\\
&=\frac{1}{2}[(\lambda_2^4+\lambda_3^4+\lambda_2^2\lambda_3^2)x_1^2+(\lambda_1^4+\lambda_3^4+\lambda_1^2\lambda_3^2)x_2^2+(\lambda_1^4+\lambda_2^4+\lambda_1^2\lambda_2^2)x_3^2+\\
&+(\lambda_1^4+\lambda_4^4+\lambda_1^2\lambda_4^2)y_1^2+(\lambda_2^4+\lambda_4^4+\lambda_2^2\lambda_4^2)y_2^2+(\lambda_3^4+\lambda_4^4+\lambda_3^2\lambda_4^2)y_3^2+\\
&+(\lambda_1^4+\lambda_5^4+\lambda_1^2\lambda_5^2)z_1^2+(\lambda_2^4+\lambda_5^4+\lambda_2^2\lambda_5^2)z_2^2+(\lambda_3^4+\lambda_5^4+\lambda_3^2\lambda_5^2)z_3^2+(\lambda_4^4+\lambda_5^4+\lambda_4^2\lambda_5^2)z_4^2],
\end{align*}
which are, more precisely, the Mishchenko's integrals of order 4 and respectively 6.\\
Thus Mishchenko's method provides two constants of motion (the above constants $K_1$ and $K_2$), which add to the Hamiltonian $H$. As mentioned in \cite{Ratiu80}, the number of independent constants of motion generated by this method is half of the dimension of the adjoint orbit only for $\mathfrak{so}(4)$. It follows that using Mishchenko's method we cannot obtain other independent constants of motion for the rigid body on $\mathfrak{so}(5)$ except the already mentioned ones.\\
In what follows, we will analyze the results of Manakov's method for finding constants of motion (\cite{Manakov76}, \cite{Ratiu80}, \cite{MoPi}), which are to be found as coefficients of the powers of $\gamma$ in the expansion of
$$\frac{1}{2r}\operatorname{Trace}(M+\gamma J^2)^r$$
for $r=2,3,4,5$.\\
For $r=2$ we obtain the first Casimir, $C_1$. For $r=3$ we obtain the constant of motion $K_1$. For $r=4$ we obtain the constant of motion $K_2$ and the second Casimir, $C_2$. Finally, for $r=5$ we obtain a constant of motion which can also be obtained by Mishchenko's method for order 8, but is dependent of the previously obtained ones, and the new functionally independent constant of motion
$$K_3(M)=\frac{1}{10}\left(\sum\limits_{j=1}^5T_j\lambda_j^2\right),$$
where
\begin{align*}
T_1&=(x_3z_2-x_2z_3-y_1z_4)^2+(x_2y_3-x_3y_2-z_1z_4)^2+(x_1x_3-y_1y_3-z_1z_3)^2+\\
&+(x_1x_2-y_1y_2-z_1z_2)^2+(x_2^2+x_3^2+y_1^2+z_1^2)^2;
\end{align*}
\begin{align*}
T_2&=(x_1z_3-x_3z_1-y_2z_4)^2+(x_3y_1-x_1y_3-z_2z_4)^2+(x_2x_3-y_2y_3-z_2z_3)^2+\\
&+(x_1x_2-y_1y_2-z_1z_2)^2+(x_1^2+x_3^2+y_2^2+z_2^2)^2;
\end{align*}
\begin{align*}
T_3&=(x_2z_1-x_1z_2-y_3z_4)^2+(x_1y_2-x_2y_1-z_3z_4)^2+(x_1x_3-y_1y_3-z_1z_3)^2+\\
&+(x_2x_3-y_2y_3-z_2z_3)^2+(x_1^2+x_2^2+y_3^2+z_3^2)^2;
\end{align*}
\begin{align*}
T_4&=(y_1z_1+y_2z_2+y_3z_3)^2+(x_3y_1-x_1y_3-z_2z_4)^2+(x_2y_3-x_3y_2-z_1z_4)^2+\\
&+(x_1y_2-x_2y_1-z_3z_4)^2+(y_1^2+y_2^2+y_3^2+z_4^2)^2;
\end{align*}
\begin{align*}
T_5&=(x_3z_2-x_2z_3-y_1z_4)^2+(x_1z_3-x_3z_1-y_2z_4)^2+(x_2z_1-x_1z_2-y_3z_4)^2+\\
&+(y_1z_1+y_2z_2+y_3z_3)^2+(z_1^2+z_2^2+z_3^2+z_4^2)^2.
\end{align*}
We find, in the particular case of $\mathfrak{so}(5)$, the result from \cite{Uga}, \cite{MoPi}, that the Mishchenko integrals of even order ($m_2=C_1$, $m_4=K_1$, $m_6=K_2$) are a subsystem of the Manakov integrals.

We obtained the following list of functionally independent and Poisson commuting constants of motion for \eqref{soM}: $H,C_1,C_2,K_1,K_2,K_3$, which proves that our system is completely integrable \cite{Manakov76}.

Additionally, as proved in Section 2, the system \eqref{soM} admits the following "generators" integrals of motion:
$$F_1(M)=\frac{x_2^2}{\lambda_1^2-\lambda_3^2}+\frac{x_3^2}{\lambda_1^2-\lambda_2^2}+\frac{y_1^2}{\lambda_1^2-\lambda_4^2}+\frac{z_1^2}{\lambda_1^2-\lambda_5^2};$$
$$F_2(M)=\frac{x_1^2}{\lambda_2^2-\lambda_3^2}+\frac{x_3^2}{\lambda_2^2-\lambda_1^2}+\frac{y_2^2}{\lambda_2^2-\lambda_4^2}+\frac{z_2^2}{\lambda_2^2-\lambda_5^2};$$
$$F_3(M)=\frac{x_1^2}{\lambda_3^2-\lambda_2^2}+\frac{x_2^2}{\lambda_3^2-\lambda_1^2}+\frac{y_3^2}{\lambda_3^2-\lambda_4^2}+\frac{z_3^2}{\lambda_3^2-\lambda_5^2};$$
$$F_4(M)=\frac{y_1^2}{\lambda_4^2-\lambda_1^2}+\frac{y_2^2}{\lambda_4^2-\lambda_2^2}+\frac{y_3^2}{\lambda_4^2-\lambda_3^2}+\frac{z_4^2}{\lambda_4^2-\lambda_5^2};$$
$$F_5(M)=\frac{z_1^2}{\lambda_5^2-\lambda_1^2}+\frac{z_2^2}{\lambda_5^2-\lambda_2^2}+\frac{z_3^2}{\lambda_5^2-\lambda_3^2}+\frac{z_4^2}{\lambda_5^2-\lambda_4^2}.$$

We continue by presenting the list of factorizations of the characteristic polynomials for the linearized equations on the tangent space to the orbit, corresponding to the equilibria from
$\mathfrak{t}_k\cap \operatorname{Orb}_{c_1;c_2}.$\\
For the equilibrium $M_{a,b}^1$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_1t^2+U_1'=0;~~V_1t^2+V_1'=0;~~W_1t^4+W_1't^2+W_1''=0;\\
U_1=(\lambda_3+\lambda_4)^2(\lambda_3+\lambda_5)(\lambda_4+\lambda_5);~~U_1'=b^2(\lambda_3-\lambda_5)(\lambda_4-\lambda_5);\\
V_1=(\lambda_1+\lambda_2)^2(\lambda_1+\lambda_5)(\lambda_2+\lambda_5);~~V_1'=a^2(\lambda_1-\lambda_5)(\lambda_2-\lambda_5);\\
W_1=(\lambda_1+\lambda_2)^4(\lambda_3+\lambda_4)^4(\lambda_1+\lambda_3)(\lambda_1+\lambda_4)(\lambda_2+\lambda_3)(\lambda_2+\lambda_4);\\
W_1''=(\lambda_1-\lambda_3)(\lambda_1-\lambda_4)(\lambda_2-\lambda_3)(\lambda_2-\lambda_4)[a^2(\lambda_3+\lambda_4)^2-b^2(\lambda_1+\lambda_2)^2]^2.
\end{array}\right.$$
For $M_{a,b}^2$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_2t^2+U_2'=0;~~V_2t^2+V_2'=0;~~W_2t^4+W_2't^2+W_2''=0;\\
U_2=(\lambda_2+\lambda_5)^2(\lambda_1+\lambda_2)(\lambda_1+\lambda_5);~~U_2'=b^2(\lambda_1-\lambda_2)(\lambda_1-\lambda_5);\\
V_2=(\lambda_3+\lambda_4)^2(\lambda_1+\lambda_3)(\lambda_1+\lambda_4);~~V_2'=a^2(\lambda_1-\lambda_3)(\lambda_1-\lambda_4);\\
W_2=(\lambda_2+\lambda_5)^4(\lambda_3+\lambda_4)^4(\lambda_2+\lambda_3)(\lambda_2+\lambda_4)(\lambda_3+\lambda_5)(\lambda_4+\lambda_5);\\
W_2''=(\lambda_2-\lambda_3)(\lambda_2-\lambda_4)(\lambda_3-\lambda_5)(\lambda_4-\lambda_5)[a^2(\lambda_2+\lambda_5)^2-b^2(\lambda_3+\lambda_4)^2]^2.
\end{array}\right.$$
For $M_{a,b}^3$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_3t^2+U_3'=0;~~V_3t^2+V_3'=0;~~W_3t^4+W_3't^2+W_3''=0;\\
U_3=(\lambda_1+\lambda_5)^2(\lambda_1+\lambda_2)(\lambda_2+\lambda_5);~~U_3'=-b^2(\lambda_1-\lambda_2)(\lambda_2-\lambda_5);\\
V_3=(\lambda_3+\lambda_4)^2(\lambda_2+\lambda_3)(\lambda_2+\lambda_4);~~V_3'=a^2(\lambda_2-\lambda_3)(\lambda_2-\lambda_4);\\
W_3=(\lambda_1+\lambda_5)^4(\lambda_3+\lambda_4)^4(\lambda_1+\lambda_3)(\lambda_1+\lambda_4)(\lambda_3+\lambda_5)(\lambda_4+\lambda_5);\\
W_3''=(\lambda_1-\lambda_3)(\lambda_1-\lambda_4)(\lambda_3-\lambda_5)(\lambda_4-\lambda_5)[a^2(\lambda_1+\lambda_5)^2-b^2(\lambda_3+\lambda_4)^2]^2.
\end{array}\right.$$
For $M_{a,b}^4$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_4t^2+U_4'=0;~~V_4t^2+V_4'=0;~~W_4t^4+W_4't^2+W_4''=0;\\
U_4=(\lambda_1+\lambda_5)^2(\lambda_1+\lambda_3)(\lambda_3+\lambda_5);~~U_4'=-b^2(\lambda_1-\lambda_3)(\lambda_3-\lambda_5);\\
V_4=(\lambda_2+\lambda_4)^2(\lambda_2+\lambda_3)(\lambda_3+\lambda_4);~~V_4'=-a^2(\lambda_2-\lambda_3)(\lambda_3-\lambda_4);\\
W_4=(\lambda_1+\lambda_5)^4(\lambda_2+\lambda_4)^4(\lambda_1+\lambda_2)(\lambda_1+\lambda_4)(\lambda_2+\lambda_5)(\lambda_4+\lambda_5);\\
W_4''=(\lambda_1-\lambda_2)(\lambda_1-\lambda_4)(\lambda_2-\lambda_5)(\lambda_4-\lambda_5)[a^2(\lambda_1+\lambda_5)^2-b^2(\lambda_2+\lambda_4)^2]^2.
\end{array}\right.$$
For $M_{a,b}^5$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_5t^2+U_5'=0;~~V_5t^2+V_5'=0;~~W_5t^4+W_5't^2+W_5''=0;\\
U_5=(\lambda_1+\lambda_5)^2(\lambda_1+\lambda_4)(\lambda_4+\lambda_5);~~U_5'=-b^2(\lambda_1-\lambda_4)(\lambda_4-\lambda_5);\\
V_5=(\lambda_2+\lambda_3)^2(\lambda_2+\lambda_4)(\lambda_3+\lambda_4);~~V_5'=a^2(\lambda_2-\lambda_4)(\lambda_3-\lambda_4);\\
W_5=(\lambda_1+\lambda_5)^4(\lambda_2+\lambda_3)^4(\lambda_1+\lambda_2)(\lambda_1+\lambda_3)(\lambda_2+\lambda_5)(\lambda_3+\lambda_5);\\
W_5''=(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)(\lambda_2-\lambda_5)(\lambda_3-\lambda_5)[a^2(\lambda_1+\lambda_5)^2-b^2(\lambda_2+\lambda_3)^2]^2.
\end{array}\right.$$
For $M_{a,b}^6$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_6t^2+U_1'=0;~~V_6t^2+V_6'=0;~~W_6t^4+W_6't^2+W_6''=0;\\
U_6=(\lambda_2+\lambda_4)^2(\lambda_2+\lambda_5)(\lambda_4+\lambda_5);~~U_6'=b^2(\lambda_2-\lambda_5)(\lambda_4-\lambda_5);\\
V_6=(\lambda_1+\lambda_3)^2(\lambda_1+\lambda_5)(\lambda_3+\lambda_5);~~V_6'=a^2(\lambda_1-\lambda_5)(\lambda_3-\lambda_5);\\
W_6=(\lambda_1+\lambda_3)^4(\lambda_2+\lambda_4)^4(\lambda_1+\lambda_2)(\lambda_1+\lambda_4)(\lambda_2+\lambda_3)(\lambda_3+\lambda_4);\\
W_6''=-(\lambda_1-\lambda_2)(\lambda_1-\lambda_4)(\lambda_2-\lambda_3)(\lambda_3-\lambda_4)[a^2(\lambda_2+\lambda_4)^2-b^2(\lambda_1+\lambda_3)^2]^2.
\end{array}\right.$$
For $M_{a,b}^7$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_7t^2+U_7'=0;~~V_7t^2+V_7'=0;~~W_7t^4+W_7't^2+W_7''=0;\\
U_7=(\lambda_3+\lambda_5)^2(\lambda_1+\lambda_3)(\lambda_1+\lambda_5);~~U_7'=b^2(\lambda_1-\lambda_3)(\lambda_1-\lambda_5);\\
V_7=(\lambda_2+\lambda_4)^2(\lambda_1+\lambda_4)(\lambda_1+\lambda_2);~~V_7'=a^2(\lambda_1-\lambda_2)(\lambda_1-\lambda_4);\\
W_7=(\lambda_2+\lambda_4)^4(\lambda_3+\lambda_5)^4(\lambda_2+\lambda_3)(\lambda_2+\lambda_5)(\lambda_3+\lambda_4)(\lambda_4+\lambda_5);\\
W_7''=-(\lambda_2-\lambda_3)(\lambda_2-\lambda_5)(\lambda_3-\lambda_4)(\lambda_4-\lambda_5)[a^2(\lambda_3+\lambda_5)^2-b^2(\lambda_2+\lambda_4)^2]^2.
\end{array}\right.$$
For $M_{a,b}^8$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_8t^2+U_1'=0;~~V_8t^2+V_8'=0;~~W_8t^4+W_8't^2+W_8''=0;\\
U_8=(\lambda_4+\lambda_5)^2(\lambda_1+\lambda_4)(\lambda_1+\lambda_5);~~U_8'=b^2(\lambda_1-\lambda_4)(\lambda_1-\lambda_5);\\
V_8=(\lambda_2+\lambda_3)^2(\lambda_1+\lambda_2)(\lambda_1+\lambda_3);~~V_8'=a^2(\lambda_1-\lambda_2)(\lambda_1-\lambda_3);\\
W_8=(\lambda_2+\lambda_3)^4(\lambda_4+\lambda_5)^4(\lambda_2+\lambda_4)(\lambda_2+\lambda_5)(\lambda_3+\lambda_4)(\lambda_3+\lambda_5);\\
W_8''=(\lambda_2-\lambda_4)(\lambda_2-\lambda_5)(\lambda_3-\lambda_4)(\lambda_3-\lambda_5)[a^2(\lambda_4+\lambda_5)^2-b^2(\lambda_2+\lambda_3)^2]^2.
\end{array}\right.$$
For $M_{a,b}^9$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_9t^2+U_1'=0;~~V_9t^2+V_9'=0;~~W_9t^4+W_9't^2+W_9''=0;\\
U_9=(\lambda_1+\lambda_4)^2(\lambda_1+\lambda_5)(\lambda_4+\lambda_5);~~U_9'=b^2(\lambda_1-\lambda_5)(\lambda_4-\lambda_5);\\
V_9=(\lambda_2+\lambda_3)^2(\lambda_2+\lambda_5)(\lambda_3+\lambda_5);~~V_9'=a^2(\lambda_2-\lambda_5)(\lambda_3-\lambda_5);\\
W_9=(\lambda_1+\lambda_4)^4(\lambda_2+\lambda_3)^4(\lambda_1+\lambda_2)(\lambda_1+\lambda_3)(\lambda_2+\lambda_4)(\lambda_3+\lambda_4);\\
W_9''=(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)(\lambda_2-\lambda_4)(\lambda_3-\lambda_4)[a^2(\lambda_1+\lambda_4)^2-b^2(\lambda_2+\lambda_3)^2]^2.
\end{array}\right.$$
For $M_{a,b}^{10}$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_{10}t^2+U_1'=0;~~V_{10}t^2+V_{10}'=0;~~W_{10}t^4+W_{10}'t^2+W_{10}''=0;\\
U_{10}=(\lambda_2+\lambda_5)^2(\lambda_2+\lambda_4)(\lambda_4+\lambda_5);~~U_{10}'=-b^2(\lambda_2-\lambda_4)(\lambda_4-\lambda_5);\\
V_{10}=(\lambda_1+\lambda_3)^2(\lambda_1+\lambda_4)(\lambda_3+\lambda_4);~~V_{10}'=a^2(\lambda_1-\lambda_4)(\lambda_3-\lambda_4);\\
W_{10}=(\lambda_1+\lambda_3)^4(\lambda_2+\lambda_5)^4(\lambda_1+\lambda_2)(\lambda_1+\lambda_5)(\lambda_2+\lambda_3)(\lambda_3+\lambda_5);\\
W_{10}''=-(\lambda_1-\lambda_2)(\lambda_1-\lambda_5)(\lambda_2-\lambda_3)(\lambda_3-\lambda_5)[a^2(\lambda_2+\lambda_5)^2-b^2(\lambda_1+\lambda_3)^2]^2.
\end{array}\right.$$
For $M_{a,b}^{11}$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_{11}t^2+U_1'=0;~~V_{11}t^2+V_{11}'=0;~~W_{11}t^4+W_{11}'t^2+W_{11}''=0;\\
U_{11}=(\lambda_2+\lambda_5)^2(\lambda_2+\lambda_3)(\lambda_3+\lambda_5);~~U_{11}'=-b^2(\lambda_2-\lambda_3)(\lambda_3-\lambda_5);\\
V_{11}=(\lambda_1+\lambda_4)^2(\lambda_1+\lambda_3)(\lambda_3+\lambda_4);~~V_{11}'=-a^2(\lambda_1-\lambda_3)(\lambda_3-\lambda_4);\\
W_{11}=(\lambda_1+\lambda_4)^4(\lambda_2+\lambda_5)^4(\lambda_1+\lambda_2)(\lambda_1+\lambda_5)(\lambda_2+\lambda_4)(\lambda_4+\lambda_5);\\
W_{11}''=-(\lambda_1-\lambda_2)(\lambda_1-\lambda_5)(\lambda_2-\lambda_4)(\lambda_4-\lambda_5)[a^2(\lambda_2+\lambda_5)^2-b^2(\lambda_1+\lambda_4)^2]^2.
\end{array}\right.$$
For $M_{a,b}^{12}$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_{12}t^2+U_1'=0;~~V_{12}t^2+V_{12}'=0;~~W_{12}t^4+W_{12}'t^2+W_{12}''=0;\\
U_{12}=(\lambda_4+\lambda_5)^2(\lambda_3+\lambda_4)(\lambda_3+\lambda_5);~~U_{12}'=b^2(\lambda_3-\lambda_4)(\lambda_3-\lambda_5);\\
V_{12}=(\lambda_1+\lambda_2)^2(\lambda_1+\lambda_3)(\lambda_2+\lambda_3);~~V_{12}'=a^2(\lambda_1-\lambda_3)(\lambda_2-\lambda_3);\\
W_{12}=(\lambda_1+\lambda_2)^4(\lambda_4+\lambda_5)^4(\lambda_1+\lambda_4)(\lambda_1+\lambda_5)(\lambda_2+\lambda_4)(\lambda_2+\lambda_5);\\
W_{12}''=(\lambda_1-\lambda_4)(\lambda_1-\lambda_5)(\lambda_2-\lambda_4)(\lambda_2-\lambda_5)[a^2(\lambda_4+\lambda_5)^2-b^2(\lambda_1+\lambda_2)^2]^2.
\end{array}\right.$$
For $M_{a,b}^{13}$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_{13}t^2+U_1'=0;~~V_{13}t^2+V_{13}'=0;~~W_{13}t^4+W_{13}'t^2+W_{13}''=0;\\
U_{13}=(\lambda_3+\lambda_5)^2(\lambda_3+\lambda_4)(\lambda_4+\lambda_5);~~U_{13}'=-b^2(\lambda_3-\lambda_4)(\lambda_4-\lambda_5);\\
V_{13}=(\lambda_1+\lambda_2)^2(\lambda_1+\lambda_4)(\lambda_2+\lambda_4);~~V_{13}'=a^2(\lambda_1-\lambda_4)(\lambda_2-\lambda_4);\\
W_{13}=(\lambda_1+\lambda_2)^4(\lambda_3+\lambda_5)^4(\lambda_1+\lambda_3)(\lambda_1+\lambda_5)(\lambda_2+\lambda_3)(\lambda_2+\lambda_5);\\
W_{13}''=(\lambda_1-\lambda_3)(\lambda_1-\lambda_5)(\lambda_2-\lambda_3)(\lambda_2-\lambda_5)[a^2(\lambda_3+\lambda_5)^2-b^2(\lambda_1+\lambda_2)^2]^2.
\end{array}\right.$$
For $M_{a,b}^{14}$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_{14}t^2+U_1'=0;~~V_{14}t^2+V_{14}'=0;~~W_{14}t^4+W_{14}'t^2+W_{14}''=0;\\
U_{14}=(\lambda_3+\lambda_5)^2(\lambda_2+\lambda_3)(\lambda_2+\lambda_5);~~U_{14}'=b^2(\lambda_2-\lambda_3)(\lambda_2-\lambda_5);\\
V_{14}=(\lambda_1+\lambda_4)^2(\lambda_1+\lambda_2)(\lambda_2+\lambda_4);~~V_{14}'=-a^2(\lambda_1-\lambda_2)(\lambda_2-\lambda_4);\\
W_{14}=(\lambda_1+\lambda_4)^4(\lambda_3+\lambda_5)^4(\lambda_1+\lambda_3)(\lambda_1+\lambda_5)(\lambda_3+\lambda_4)(\lambda_4+\lambda_5);\\
W_{14}''=-(\lambda_1-\lambda_3)(\lambda_1-\lambda_5)(\lambda_3-\lambda_4)(\lambda_4-\lambda_5)[a^2(\lambda_3+\lambda_5)^2-b^2(\lambda_1+\lambda_4)^2]^2.
\end{array}\right.$$
For $M_{a,b}^{15}$ the eigenvalues are the roots of equations of the following form:
$$\left.\begin{array}{l}
U_{15}t^2+U_1'=0;~~V_{15}t^2+V_{15}'=0;~~W_{15}t^4+W_{15}'t^2+W_{15}''=0;\\
U_{15}=(\lambda_4+\lambda_5)^2(\lambda_2+\lambda_5)(\lambda_2+\lambda_4);~~U_{15}'=b^2(\lambda_2-\lambda_4)(\lambda_2-\lambda_5);\\
V_{15}=(\lambda_1+\lambda_3)^2(\lambda_1+\lambda_2)(\lambda_2+\lambda_3);~~V_{15}'=-a^2(\lambda_1-\lambda_2)(\lambda_2-\lambda_3);\\
W_{15}=(\lambda_1+\lambda_3)^4(\lambda_4+\lambda_5)^4(\lambda_1+\lambda_4)(\lambda_1+\lambda_5)(\lambda_3+\lambda_4)(\lambda_3+\lambda_5);\\
W_{15}''=(\lambda_1-\lambda_4)(\lambda_1-\lambda_5)(\lambda_3-\lambda_4)(\lambda_3-\lambda_5)[a^2(\lambda_4+\lambda_5)^2-b^2(\lambda_1+\lambda_3)^2]^2.
\end{array}\right.$$
Without loss of generality, we can choose an ordering for $\lambda_i$'s, namely
$$\lambda_1>\lambda_2>\lambda_3>\lambda_4>\lambda_5.$$

\medskip

Since for the equilibria $\mathfrak{t}_3 \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_4 \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_5 \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{10} \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{11} \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{13} \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{14} \cap \operatorname{Orb}_{c_1;c_2}$ and $\mathfrak{t}_{15} \cap \operatorname{Orb}_{c_1;c_2}$ at least one of the eigenvalues is real and strictly positive we have nonlinear {\it instability} for these equilibria.

\medskip

It is immediate to prove that if
$$a^2(\lambda_2+\lambda_4)^2\not= b^2(\lambda_1+\lambda_3)^2$$
is equivalent from Eqs. \ref{values_a_b} with
\begin{equation}
\label{special1}
c_1^2\left[(\lambda_2+\lambda_4)^4+(\lambda_1+\lambda_3)^4\right]\not= c_2\left[(\lambda_2+\lambda_4)^2+(\lambda_1+\lambda_3)^2\right]^2,
\end{equation}
and in this case we have $W_6''<0$.
Since $W_6>0$ and $W_6''<0$, it follows that if the condition \eqref{special1} holds, then the equilibrium $M_{a,b}^6$ from $\mathfrak{t}_6 \cap \operatorname{Orb}_{c_1;c_2}$ (together with the equilibria $M_{-a,b}^6$, $M_{a,-b}^6$, $M_{-a,-b}^6$) are also nonlinear {\it unstable}.

On the other hand, since
$$b^2(\lambda_2+\lambda_4)^2< a^2(\lambda_1+\lambda_3)^2,$$
(because $a>b$ and $\lambda_1+\lambda_3>\lambda_2+\lambda_4$)
it follows that the equilibrium $M_{b,a}^6$ from $\mathfrak{t}_6 \cap \operatorname{Orb}_{c_1;c_2}$ (together with the equilibria $M_{b,-a}^6$, $M_{-b,a}^6$, $M_{-b,-a}^6$) are nonlinear {\it unstable}, regardless if condition \eqref{special1} holds or not.

\medskip

A similar discussion holds for the equilibria from $\mathfrak{t}_7\cap \operatorname{Orb}_{c_1;c_2}$.

\medskip

We continue with the study of the remaining families of equilibria, which is $\mathfrak{t}_1\cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_2\cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_8\cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_9\cap \operatorname{Orb}_{c_1;c_2}$ and $\mathfrak{t}_{12}\cap \operatorname{Orb}_{c_1;c_2}$.

For start, let us notice that energy methods using the "canonical" integrals of motion (the Hamiltonian and/or the constants of motion $K_1,K_2,K_3$) work only for a few of the above equilibria. But choosing particular linear combinations of some "generators" integrals of motion introduced in Section 2 proves to be very effective in solving the stability problem with energy methods.

To begin with, we note that the method involving linearization used above is inconclusive for the stability of the $M_{a,b}^1$ equilibrium.
Therefore, we will use Arnold's method \cite{Arnold65}, which is equivalent with the other energy methods \cite{BiPu07}.

Consider the smooth function $G_{mn}\in
C^{\infty}(\mathfrak{so}(5),\mathbb{R})$, where $m,n$ are real numbers,
$$G_{mn}(M)=F_1(M)+F_5(M)+m C_1(M)+n C_2(M).$$
Choosing $m,n$ such that $\mathbf{d}G_{mn}(M_{a,b}^1)=0$, namely
$$m=\frac{2b^2}{(\lambda_1^2-\lambda_2^2)(a^2-b^2)},~~~n=-\frac{2}{(\lambda_1^2-\lambda_2^2)(a^2-b^2)}$$
and taking into account that
$$W:=\ker dC_1(M_{a,b}^1)\cap \ker dC_2(M_{a,b}^1)=\operatorname{span} (E_1,E_2,E_4,E_5,E_7,E_8,E_9,E_{10})$$
we obtain the determinants associated to all upper-left submatrices of the
Hessian $\mathbf{d}^2G_{mn}(M_{a,b}^1)|_{W\times W}$ as follows:
\begin{align*}
D_1&=-\frac{2a^2}{(\lambda_1^2-\lambda_2^2)(a^2-b^2)}<0;\\
D_2&=\frac{4a^2[a^2(\lambda_2^2-\lambda_3^2)+b^2(\lambda_1^2-\lambda_2^2)]}{(\lambda_1^2-\lambda_2^2)^2(\lambda_1^2-\lambda_3^2)(a^2-b^2)^2}>0;\\
D_3&=-\frac{8a^2(\lambda_2^2-\lambda_4^2)[a^2(\lambda_2^2-\lambda_3^2)+b^2(\lambda_1^2-\lambda_2^2)]}{(\lambda_1^2-\lambda_2^2)^3(\lambda_1^2-\lambda_3^2)(\lambda_1^2-\lambda_4^2)(a^2-b^2)^2}<0;\\
D_4&=\frac{16a^4(\lambda_2^2-\lambda_3^2)(\lambda_2^2-\lambda_4^2)}{(\lambda_1^2-\lambda_2^2)^4(\lambda_1^2-\lambda_3^2)(\lambda_1^2-\lambda_4^2)(a^2-b^2)^2}>0;\\
D_5&=-\frac{32a^4(\lambda_2^2-\lambda_4^2)(\lambda_2^2-\lambda_3^2)}{(\lambda_1^2-\lambda_2^2)^5(\lambda_1^2-\lambda_3^2)(\lambda_1^2-\lambda_4^2)(a^2-b^2)^2}<0;\\
D_6&=\frac{64a^4(\lambda_1^2-\lambda_5^2)(\lambda_2^2-\lambda_3^2)(\lambda_2^2-\lambda_4^2)}{(\lambda_1^2-\lambda_2^2)^6(\lambda_1^2-\lambda_3^2)(\lambda_1^2-\lambda_4^2)(\lambda_2^2-\lambda_5^2)(a^2-b^2)^2}>0;\\
D_7&=-\frac{128a^4(\lambda_1^2-\lambda_5^2)(\lambda_2^2-\lambda_3^2)(\lambda_2^2-\lambda_4^2)}{(\lambda_1^2-\lambda_2^2)^6(\lambda_1^2-\lambda_3^2)(\lambda_1^2-\lambda_4^2)(\lambda_2^2-\lambda_5^2)(\lambda_3^2-\lambda_5^2)(a^2-b^2)^2}<0;\\
D_8&=\frac{256a^4(\lambda_1^2-\lambda_5^2)(\lambda_2^2-\lambda_3^2)(\lambda_2^2-\lambda_4^2)}{(\lambda_1^2-\lambda_2^2)^6(\lambda_1^2-\lambda_3^2)(\lambda_1^2-\lambda_4^2)(\lambda_2^2-\lambda_5^2)(\lambda_3^2-\lambda_5^2)(\lambda_4^2-\lambda_5^2)(a^2-b^2)^2}>0.
\end{align*}
Consequently, $\mathbf{d}^2G_{mn}(M_{a,b}^1)|_{W\times W}$ is negative definite, which implies nonlinear {\it stability} for the equilibrium $M_{a,b}^1$. The same computations lead to nonlinear stability for the equilibria $M_{a,-b}^1$, $M_{-a,b}^1$ and $M_{-a,-b}^1$ from $\mathfrak{t}_1\cap \operatorname{Orb}_{c_1;c_2}$.

For the remaining equilibria $M_{b,a}^1$, $M_{-b,a}^1$, $M_{b,-a}^1$ and $M_{-b,-a}^1$ from $\mathfrak{t}_1\cap \operatorname{Orb}_{c_1;c_2}$ a convenient energy function for applying Arnold's method is $G_{mn}(M)=F_4(M)+F_5(M)+m C_1(M)+n C_2(M)$. Thus, all eight equilibria in $\mathfrak{t}_1\cap \operatorname{Orb}_{c_1;c_2}$ are nonlinear {\it stable}.

\medskip

For the equilibria $M_{a,b}^2$, $M_{-a,b}^2$, $M_{a,-b}^2$ and $M_{-a,-b}^2$ from $\mathfrak{t}_2\cap \operatorname{Orb}_{c_1;c_2}$ a convenient energy function for applying Arnold's method is $G_{mn}(M)=F_1(M)-F_4(M)+m C_1(M)+n C_2(M)$. Thus, the equilibria $M_{a,b}^2$, $M_{-a,b}^2$, $M_{a,-b}^2$ and $M_{-a,-b}^2$ are nonlinear {\it stable}.

\medskip

For the equilibria $M_{a,b}^8$, $M_{-a,b}^8$, $M_{a,-b}^8$ and $M_{-a,-b}^8$ from $\mathfrak{t}_8\cap \operatorname{Orb}_{c_1;c_2}$ a convenient energy function for applying Arnold's method is $G_{mn}(M)=F_1(M)+F_5(M)+m C_1(M)+n C_2(M)$. For the equilibria $M_{b,a}^8$, $M_{-b,a}^8$, $M_{b,-a}^8$ and $M_{-b,-a}^8$ from $\mathfrak{t}_8\cap \operatorname{Orb}_{c_1;c_2}$ a convenient energy function for applying Arnold's method is $G_{mn}(M)=F_1(M)+F_2(M)+m C_1(M)+n C_2(M)$. Thus, all eight equilibria in $\mathfrak{t}_8\cap \operatorname{Orb}_{c_1;c_2}$ are nonlinear {\it stable}.

\medskip

For the equilibria $M_{b,a}^9$, $M_{-b,a}^9$, $M_{b,-a}^9$ and $M_{-b,-a}^9$ from $\mathfrak{t}_9\cap \operatorname{Orb}_{c_1;c_2}$ a convenient energy function for applying Arnold's method is $G_{mn}(M)=F_4(M)-F_5(M)+m C_1(M)+n C_2(M)$. Thus, the equilibria $M_{b,a}^9$, $M_{-b,a}^9$, $M_{b,-a}^9$ and $M_{-b,-a}^9$ from $\mathfrak{t}_9\cap \operatorname{Orb}_{c_1;c_2}$ are nonlinear {\it stable}.

\medskip

For the equilibria $M_{a,b}^{12}$, $M_{-a,b}^{12}$, $M_{a,-b}^{12}$ and $M_{-a,-b}^{12}$ from $\mathfrak{t}_{12}\cap \operatorname{Orb}_{c_1;c_2}$ a convenient energy function for applying Arnold's method is $G_{mn}(M)=F_1(M)-F_3(M)+m C_1(M)+n C_2(M)$. For the equilibria $M_{b,a}^{12}$, $M_{-b,a}^{12}$, $M_{b,-a}^{12}$ and $M_{-b,-a}^{12}$ from $\mathfrak{t}_{12}\cap \operatorname{Orb}_{c_1;c_2}$ a convenient energy function for applying Arnold's method is $G_{mn}(M)=F_3(M)-F_5(M)+m C_1(M)+n C_2(M)$. Thus, all eight equilibria in $\mathfrak{t}_{12}\cap \operatorname{Orb}_{c_1;c_2}$ are nonlinear {\it stable}.

\medskip

For the equilibria $M_{b,a}^2$, $M_{-b,a}^2$, $M_{b,-a}^2$ and $M_{-b,-a}^2$ from $\mathfrak{t}_2\cap \operatorname{Orb}_{c_1;c_2}$ and for the equilibria $M_{a,b}^9$, $M_{-a,b}^9$, $M_{a,-b}^9$ and $M_{-a,-b}^9$ from $\mathfrak{t}_9\cap \operatorname{Orb}_{c_1;c_2}$ the stability problem remains open and it is likely that a bifurcation phenomenon occurs. Such a phenomenon appears for some equilibria in the case of
the $\mathfrak{so}(4)$ free rigid body and is extensively studied in \cite{noi}.

\medskip

We proved the following result:
\begin{thm}
\begin{itemize}
\item[{\rm (i)}] The equilibria from $\mathfrak{t}_3 \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_4 \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_5 \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{10} \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{11} \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{13} \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_{14} \cap \operatorname{Orb}_{c_1;c_2}$ and $\mathfrak{t}_{15} \cap \operatorname{Orb}_{c_1;c_2}$ are unstable.
\item[{\rm (ii)}] The equilibria $M_{b,a}^6$, $M_{b,-a}^6$, $M_{-b,a}^6$ and $M_{-b,-a}^6$ from $\mathfrak{t}_6 \cap \operatorname{Orb}_{c_1;c_2}$ and respectively the equilibria $M_{b,a}^7$, $M_{b,-a}^7$, $M_{-b,a}^7$ and $M_{-b,-a}^7$ from $\mathfrak{t}_7 \cap \operatorname{Orb}_{c_1;c_2}$ are unstable.
\item[{\rm (iii)}] a) If the condition
$$c_1^2\left[(\lambda_1+\lambda_3)^4+(\lambda_2+\lambda_4)^4\right]\not= c_2\left[(\lambda_1+\lambda_3)^2+(\lambda_2+\lambda_4)^2\right]^2$$
holds, then the equilibria $M_{a,b}^6$, $M_{-a,b}^6$, $M_{a,-b}^6$ and $M_{-a,-b}^6$ from $\mathfrak{t}_6 \cap \operatorname{Orb}_{c_1;c_2}$ are unstable.\\
b) If the condition
$$c_1^2\left[(\lambda_2+\lambda_4)^4+(\lambda_3+\lambda_5)^4\right]\not= c_2\left[(\lambda_2+\lambda_4)^2+(\lambda_3+\lambda_5)^2\right]^2$$
holds, then the equilibria $M_{a,b}^7$, $M_{-a,b}^7$, $M_{a,-b}^7$ and $M_{-a,-b}^7$ from $\mathfrak{t}_7 \cap \operatorname{Orb}_{c_1;c_2}$ are unstable.
\item[{\rm (iv)}] The equilibria from $\mathfrak{t}_1 \cap \operatorname{Orb}_{c_1;c_2}$, $\mathfrak{t}_8 \cap \operatorname{Orb}_{c_1;c_2}$ and $\mathfrak{t}_{12} \cap \operatorname{Orb}_{c_1;c_2}$ are nonlinear stable.
\item[{\rm (v)}] The equilibria $M_{a,b}^2$, $M_{-a,b}^2$, $M_{a,-b}^2$ and $M_{-a,-b}^2$ from $\mathfrak{t}_2 \cap \operatorname{Orb}_{c_1;c_2}$ and respectively the equilibria $M_{b,a}^9$, $M_{-b,a}^9$, $M_{b,-a}^9$ and $M_{-b,-a}^9$ from $\mathfrak{t}_9\cap \operatorname{Orb}_{c_1;c_2}$ are nonlinear {\it stable}.
\end{itemize}
\end{thm}

\medskip

\noindent {\bf Acknowledgements.} The author acknowledges the very helpful discussions with Professors Tudor Ra\c tiu and Petre Birtea.


\subsection*{Abstract}
\textbf{
Anyons -- particles carrying fractional statistics that interpolate between bosons and fermions -- have been conjectured to exist in low dimensional systems. In the context of the fractional quantum Hall effect (FQHE), quasi-particles made of electrons take the role of anyons whose statistical exchange phase is fixed by the filling factor.
Here we propose an experimental setup to create anyons in one-dimensional lattices with fully tuneable exchange statistics. In our setup, anyons are created by bosons with occupation-dependent hopping amplitudes, which can be realized by assisted Raman tunneling. The statistical angle can thus be controlled in situ by modifying the relative phase of external driving fields. This opens the fascinating possibility of smoothly transmuting bosons via anyons into fermions and of inducing a phase transition by the mere control of the particle statistics as a free parameter.
In particular, we demonstrate how to induce a quantum phase transition from a superfluid into an exotic Mott-like state where the particle distribution exhibits plateaus at fractional densities.
}

\subsection*{Introduction}

Usually, every particle in quantum theory is classified as either a boson -- a particle joining any number of identical particles in a single quantum state -- or a fermion, characterized by the sole occupancy of its state.
The exchange of two fermions leads due to the Pauli principle to a phase factor -1 in the total wavefunction, while the wavefunction of two bosons remains invariant under particle exchange.
More than 30 years ago, researchers proposed a third fundamental category of particles living in two-dimensional (2D) systems, ``anyons'' \cite{anyon1,anyon2,anyon3,Tsui, Laughlin}.  For two anyons, the wavefunction acquires a fractional phase $e^{i \theta}$ under particle exchange, giving rise to \emph{fractional statistics}, with $0<\theta<\pi$.
For a few years the physics of anyons remained restricted to the
2D world\cite{Girvin}, until Haldane presented the concept of fractional
statistics in arbitrary dimensions\cite{Haldane}.

\begin{figure} []
  \centering
  \includegraphics[width=1\linewidth]{Fig1.pdf}
  \vspace{-0.5cm}
  \caption{
  \textbf{Anyon-Boson Mapping and schematic of the proposed experiment.} (a) Anyons in 1D can be mapped onto bosons featuring occupation-dependent hopping amplitudes.  (b) Assisted Raman tunneling can selectively address hopping processes connecting different occupational states and induce a relative phase, realizing a fully tuneable particle exchange statistics angle $\theta$. Energies are not in scale.
}
  \label{fig:scheme}
\end{figure}

Anyons in one dimension (1D) are still unexplored to a wide extent.
Recently, it has been put forward to create fractional statistics in a 1D
Hubbard model of fermions with correlated hopping processes\cite{bond-charge}. Anyons are realized in this case as low-energy elementary excitations.

In this Article, we introduce an exact mapping between anyons and bosons in 1D. We show that a Hubbard model of anyons is equivalent to a variant of the Bose-Hubbard model\cite{BHM} where the bosonic hopping amplitudes are state-dependent. This conditional-hopping phase factor breaks reflection parity in the system, which is an important ingredient to realize fractional statistics\cite{Wilczek}. We propose to realize bosons with conditional-hopping amplitudes using assisted Raman tunneling in an optical lattice\cite{Ignacio,Greiner} (OL). We discuss how the direct control of the statistical phase can induce a quantum phase transition from a bosonic superfluid into a Mott-like state, exhibiting exotic Mott shells at fractional densities. The ``statistical ramp'' transmutes bosons smoothly into ``pseudofermions'', with anyons as intermediate steps.

Anyons in 1D are defined\cite{Batchelor, Kundu} by the generalized
commutation relations \begin{equation}
a_{j}a_{k}^{\dagger}-e^{-i\theta\mathrm{sgn}(j-k)}a_{k}^{\dagger}a_{j}=\delta_{jk},\hspace{1em}a_{j}a_{k}=e^{i\theta\mathrm{sgn}(j-k)}a_{k}a_{j},\label{CCR}\end{equation}
where the operators $a_{j}^{\dagger},a_{j}$ create or annihilate
an anyon on site $j$, respectively. The sign function is such that
$\mathrm{sgn}(j-k)=0$ for $j=k$, hence two particles on the
same site behave as ordinary bosons. In consequence, anyons with statistics
$\theta=\pi$ are pseudofermions:  while being bosons on-site, they are fermions off-site.

\subsection*{Results}

\subsubsection*{Mapping between Anyons and Bosons}

We introduce an exact mapping between anyons and bosons in 1D.
Let us define the \emph{fractional version} of a Jordan-Wigner transformation,
\begin{equation}
a_{j}=b_{j}\exp\left(i\theta\sum_{i=1}^{j-1}n_{i}\right)\label{mapping}\end{equation}
with $n_{i}=a_{i}^{\dagger}a_{i}=b_{i}^{\dagger}b_{i}$ the number
operator for both particle types. Provided that the particles of type
$b$ are bosons, $[b_{j},b_{i}^{\dagger}]=\delta_{ji}$ and $[b_{j},b_{i}]=0$,
we can show that the mapped operators $a$ indeed obey the anyonic
commutation relations as introduced in equation (\ref{CCR}). For a proof see Methods. This mapping
elucidates that anyons in 1D are indeed non-local quasi-particles, made of bosons with an attached string operator.

Our ultimate goal is to propose a realistic setup for demonstrating
an interacting gas of anyons in 1D OLs. We therefore introduce the Anyon-Hubbard model
\begin{equation}
H^{a}=-J\sum_{j}^L (a_{j}^{\dagger}a_{j+1}+\mathrm{h.c.})+\frac{U}{2}\sum_{j}^L n_{j}(n_{j}-1),\label{AHM}\end{equation}
where $J$ is the tunneling amplitude connection two neighbouring sites
and $U$ is the on-site interaction energy. By inserting the Anyon-Boson
mapping, equation (\ref{mapping}), the Hamiltonian $H^{a}$ can be rewritten
in terms of bosonic operators: \begin{equation}
H^{b}=-J\sum_{j}^L (b_{j}^{\dagger}b_{j+1}e^{i\theta n_{j}}+\mathrm{h.c.})+\frac{U}{2}\sum_{j}^L n_{j}(n_{j}-1).\label{bosons}\end{equation}

The mapped, bosonic Hamiltonian thus describes bosons with a occupation-dependent amplitude $Je^{i\theta n_{j}}$ for hopping processes from right to left ($j+1 \rightarrow j$). If the target site $j$ is unoccupied,
the hopping amplitude is simply $J$. If it is occupied by one boson,
the amplitude reads $Je^{i\theta}$, and so on. The conditional hopping scheme is depicted in Fig. 1a.
We emphasize that the non-local mapping between anyons and bosons, equation (\ref{mapping}), leads luckily to a purely local, and thus viable Hamiltonian.

As expected from anyons, the reflection parity symmetry is broken\cite{Wilczek} at the level of the commutation relations, equation (\ref{CCR}).
The fractional Jordan-Wigner transformation, equation (\ref{mapping}), transfers this asymmetry also to the bosonic case:
The resulting Hamiltonian, equation (\ref{bosons}), features a phase factor acting only on the left site $j$ and thus violates parity.
Even in the absence of the on-site interaction, $U=0$, the exponential operator in equation (\ref{bosons}) gives rise to many-body interactions, as expected for anyons\cite{Wilczek}.

In the limit $U/J\to\infty$, bosons are impenetrable and each site
contains at most one particle. In this case, the phase factor disappears
and the bosonic Hamiltonian $H^{b}$ reduces to an ordinary Tonks-Girardeau
gas\cite{Tonks1,Tonks2,Tonks3,Tonks4}. However, we consider local occupation
numbers beyond the hard-core limit, $n_{j}>1$. Thus anyons can exchange
positions, changing the phase of the total wavefunction, and show
non-trivial features.

\subsubsection*{Experimental realization}

In our experimental proposal, the key point for realizing anyonic statistics is
to induce a hopping term with a phase shift which depends on the
occupation of the left-hand site $j$, Fig.~1b displays the basic concept. In order to distinguish between different local occupational states, we require a non-zero on-site interaction $U$.
For simplicity, let us consider lattice site occupations that are restricted to $n_j=0,1,2$ (higher local truncations are also possible, see Methods). The occupation-dependent tunneling and thus the conditional-hopping model equation (\ref{bosons}) can be implemented in OLs with present experimental techniques.
We propose to employ an assisted tunneling scheme, based on ideas by Jaksch and Zoller\cite{jaksch03} and Juzeliunas \emph{et al.}\cite{fleischhauer}.
The OL is tilted, with an energy offset $\Delta$ between neighbouring sites, this additional field gradient breaks reflection parity.
Two different occupational states (note that the occupational state $n_j=0$ is not relevant) in either of the two sites form in total a 4-dimensional atomic ground state manifold, which we propose to couple to an excited state $|e\rangle$ via four external driving fields (labeled 1, 2, 3 and 4 in Fig. 1b). According to this notation, singly- and doubly-occupied states are
 coupled by fields 2 and 1 in the left site and by 3 and 4 in the right site, respectively.
The excited state can be experimentally realized in at least two alternative ways.\\
First, two spin-dependent lattices\cite{Mandel, Campbell, McKay} can be employed. In the case of Rb$^{87}$, one lattice for example traps atoms in the $F=1, m_F=-1$ hyperfine state, assigned to the ground state manifold. The excited state $|e\rangle$ can then be engineered as a vibrational state of a second lattice, trapping atoms in the $F=1, m_F=0$ hyperfine state. Atoms in the excited state would then be localized between the left and right wells of the $F=1, m_F=-1$ lattice, but not necessarily in their center.\\
This implementation would offer the advantage of external driving fields in the radio-frequency regime. Such frequencies could then still resolve\cite{Campbell} the typical energies $U$ and $\Delta$ (both of the order of a few kHz), which is necessary to selectively couple to the four different states in the ground state manifold.\\
Second, one can use two optical lattices, and trap ground state manifold atoms in the red-detuned lattice, while the excited state would live in the blue-detuned one. The driving fields necessary in this case would be then, however, typically in the THz frequency range, making a precise resolution of $U$ and $\Delta$ more challenging for the experimentalist: in principle a laser with a linewidth $\delta_{\textrm{linewidth}}  \ll U, \Delta$ would be needed.\\
The effective tunneling rates $J_{ab}$ ($a \in \{1,2\}, b \in \{3,4\}$) are calculated for four $\Lambda$-systems, laser frequencies $\omega_i$, and Rabi frequencies $\Omega_i$ ($i \in \{1,2,3,4\}$) via adiabatic elimination, see Methods. We emphasize that the tilt energy $\Delta$ disappears in the effective Hamiltonian after rotating out time-dependent phase factors: indeed this energy offset is absorbed (or released) by the external radiation field, yielding a total Hamiltonian without a tilt term (see also Jaksch and Zoller\cite{jaksch03} on this issue).  \\
The following conditions on the effective tunneling rates $J_{ab}$ have to be satisfied in order to realize our model in equation (\ref{bosons}):
\begin{eqnarray}
J_{23} &=& J_{24} \equiv J \, ,\label{condition1}  \\
J_{13} &=& J_{14} \equiv J e^{i \theta} \label{condition2} \, ,
\end{eqnarray}
where $\theta$ is the anyonic exchange phase.
For a more detailed consideration of realistic energy scales and appropriate parameter regimes, see Methods.

\begin{figure} []
  \centering
  \includegraphics[width=1.05\linewidth]{Fig2.pdf}
  \vspace{-0.5cm}
  \caption{
  \textbf{Density distribution in quasi-momentum space \braket{n_k} as a function of particle statistics $\theta$.}  The shift of the density peaks with increasing $\theta$ displays a characteristic quadratic behaviour. The  fit to the trace of density maxima is depicted in white and yields fitting parameters $k_0= 0.9828 \pi$, $\alpha= -1.0511/ \pi$, $\beta= 0.9982 \pi$. The inset displays the decrease of the peak occupations with $\theta$ (yellow circles) and indicates the \emph{statistically induced} phase transition from a superfluid to a Mott-like state. Parameters: $L=30, N=31, U/J=0.2$.
}
  \label{fig:nk}
\end{figure}

\subsubsection*{Density in quasi-momentum space}

We have computed the ground state wave function for the conditional-hopping Bose-Hubbard model, equation (\ref{bosons}), using the Density Matrix Renormalization Group (DMRG)\cite{White,DMRG}.
In Fig. 2 we plot the quasi-momentum distribution
\begin{equation}
\langle n_k \rangle = \frac{1}{L} \sum_{ij} e^{ik(x_i-x_j)} \langle b_{i}^\dagger b_{j}\rangle
\end{equation}
as a function of the statistical phase angle $\theta$. The case $\theta = 0$ corresponds to ordinary bosons, which for the parameters chosen quasi-condense. The density distribution in quasi-momentum space is thus peaked at $k=0$. Increasing $\theta$ to non-zero values, we find that the position of the peaks $\theta_{max}(k)$ is shifted as a non-linear function of $k$. Indeed, for fillings $N/L>1$ one finds a quadratic behaviour $\theta_{max}(k) = \alpha (k-k_0)^2 + \beta$. For fillings close to $N/L=1$, $\beta\rightarrow \pi, k_0 \rightarrow \pi$ and $\alpha \rightarrow -1/\pi$. For higher fillings $N/L\rightarrow 2$, these fitting parameters are altered, however, the characteristic quadratic dependence is conserved.\\
In this analysis, we find two important characteristics of conditional hopping bosons and thus anyons in 1D. The quadratic dependence of $\theta_{max}$ on $k$ contrasts ordinary magnetic fields (with a constant phase factor $e^{i \theta}$ in the kinetic Hamiltonian). In this case, the shift of the peaks depends linearly on the phase angle $\theta$.
In the anyonic case, however, the growth of correlations with increasing $\theta$ may induce the characteristic quadratic trace, which could be directly observed in OL experiments using standard time-of-flight imaging.

An even more important observation is that the contrast of the peaks (and the phase coherence) decays with increasing $\theta$. The peak values $\braket{n_k(\theta=\theta_{max})}$ are plotted in the inset of Fig. 2, in the pseudofermionic limit ($\theta \rightarrow \pi$) the peak is almost washed out. This suggests that an increase of $\theta$ transforms the initial quasi-condensate into a quantum state where all phase coherence is lost. It will become evident in the subsequent paragraphs that this quantum state will turn out to be a Mott-like state, with Mott plateaus emerging at fractional densities. We emphasize that this quantum phase transition is only driven by the statistical angle $\theta$, all other parameters such as $J/U$ are fixed. The loss of coherence can be understood as follows. With increasing $\theta$ the occupation-dependent phase factor in equation (\ref{bosons}) becomes more and more important: Hopping processes connecting sites with different occupations will contribute different phases and will cancel out in the kinetic Hamiltonian due to incoherent superpositions. This destructive interference effect is amplified by an increasing $\theta$ and induces the localization of particles, yielding an insulating phase.
We emphasize that the present analysis of the density distribution in momentum space refers to the bosonic particles only, $\langle n_k \rangle=\langle b_k^{\dagger} b_k \rangle$. Namely, Fig. 2 represents what would be observed in the experiment. However, while the mapping in equation (\ref{mapping}) establishes a 1-1 correspondence between the number operators in real space, $\langle a_j^{\dagger}a_j\rangle=\langle b_j^{\dagger}b_j\rangle$, the density distributions in momentum space may differ significantly. In this sense, a study of $\tilde{n}_k=\langle a_k^{\dagger}a_k\rangle$ and the superfluid order parameter in the original anyonic model (\ref{AHM}) would be very interesting, but it is beyond the scope of this paper. A study of density distributions in momentum and real space was recently presented for the particular case of hard-core anyons\cite{Hao}.

\subsubsection*{Phase diagram}

Next, we present the phase diagram for conditional-hopping bosons in the $(\mu, J/U)$-plane.  Just as in the case of ordinary bosons ($\theta=0$) -- the celebrated phase diagram of the Bose-Hubbard model\cite{BHM} -- the anyonic version will display insulating and superfluid regions. However, we find that the size of the insulating regions (Mott lobes) grows with increasing statistical angle $\theta$ -- a fact that will be central in our proposal to design phase transitions by tuning the particle statistics. \\
The phase diagram is calculated as follows. We start with a unit filling ground state $\ket{N=L}$ and compare the energies with the ground states of $\ket{N=L\pm1}$. This yields the gap energies $\Delta E^{\pm} (L)$, corresponding to the energy required to add or subtract one particle, respectively. These gap energies were calculated using DMRG for system sizes $L=15, 30, 60, 90$ and $120$. From the finite-size scaling we extrapolate the infinite-size values $\Delta E^{\pm}$, which are plotted in Fig. 3 in the $(\mu, J/U)$-plane. Note that a Mott insulator is associated with a non-zero gap $\epsilon = \Delta E^+ - \Delta E^-$, while for the superfluid phase $\epsilon =0$ in the thermodynamic limit. For the bosonic case $\theta=0$ we recover the first Mott lobe as in the work by K\"uhner and Monien\cite{BHM1D}. Turning on the statistical angle and transmuting the bosons into anyons, we observe an expansion of the Mott-like insulating phase in both dimensions $\mu$ and $J/U$. The phase transition points $(U/J)_{crit}$ (defined by the cusps of the Mott lobes) are plotted in the inset of Fig. 3 as a function of $\theta$.
In the pseudofermion limit  $\theta \rightarrow \pi$, the Mott-like phase seems to extend to very large values of $J/U$. In contrast, ordinary bosons would form a superfluid state in this parameter regime.
The expansion of the Mott lobes with $\theta$ is also observed in mean-field calculations for our model, equation (\ref{bosons}). The mean field solution produces very interesting patterns
for the transition lines, as shown in Fig. \ref{fig:lobes MF}. For details of this calculation, see Methods. \\
We note that the phase diagrams for conditional-hopping bosons and for anyons are the same. Due to the unitarity of the mapping (\ref{mapping}), the two models (\ref{AHM}) and (\ref{bosons}) are isospectral, and thus feature the same energy gaps and phase diagrams.

\begin{figure} []
  \centering
  \includegraphics[width=1.08\linewidth]{Fig3.pdf}
  \vspace{-0.5cm}
  \caption{
  \textbf{Phase diagram of the anyonized gas.} The Mott lobe (corresponding to $\braket{n}=1$) expands with increasing statistical angle $\theta$ in both directions in the $(\mu, J/U)$-plane. This demonstrates the novel possibility to induce a quantum phase transition from the superfluid into the insulating, Mott-like phase, by simply changing the particle statistics. The energies $\Delta E^{\pm}$ are symbolized by circles (crosses), respectively. The critical phase transition points $(U/J)_{crit}$ are shown on the inset, as a function of $\theta$. In the extreme limit $\theta \rightarrow \pi$, $(U/J)_{crit}$ tends to zero, i.e. pseudofermions seem to be always in the insulating phase.
}
  \label{fig:phases}
\end{figure}

\subsection*{Discussion}

We envision the following procedure to demonstrate the first \emph{statistically induced phase transition}:
We fix the parameters at $J/U = 0.5$, $N/L = 1$ and $\mu \simeq 0$. We start with a phase detuning $\theta=0$ between the external driving fields (see Methods), and thus realize a superfluid bosonic gas as the ground state. The phase angle $\theta$ is now continuously increased, leading to anyonization of the gas and growth of the Mott-like phase. At a critical value $\theta_c$ the phase border will surpass the fixed point in parameter space, which will be then located inside the Mott phase. The critical angle can be estimated from the phase diagram to be in the range $\theta_c \in [\pi/2,3\pi/4]$. For laser detunings beyond this critical angle, $\theta > \theta_c$, the gas will be in an insulating, Mott-like state. In this way, the variation of the particle statistics in our proposal directly realizes a novel superfluid-to-Mott quantum phase transition.

An intriguing aspect of the Mott-like state emerges when a harmonic trapping potential is added to the system, $H^b_{tr} = H^b + V \sum_i ((L+1)/2-i)^2 n_i$. This simulates the experimental conditions also to a more realistic extent. In a harmonic trap, ordinary Mott insulators form real-space density distributions $\braket{n_i}$ that resemble ``wedding cakes'' \cite{BHM}. Due to the broken translational invariance (induced by the trap), the chemical potential in the local density approximation now is a function of the lattice sites. The distribution $\braket{n_i}$ thus exhibits plateaus at integer densities, for ordinary bosons in a Mott-insulating state.

\begin{figure} []
  \centering
  \includegraphics[width=1.02\linewidth]{Fig4.pdf}
  \vspace{-0.5cm}
  \caption{
  \textbf{Fractional Mott plateaus.} The density distributions $\braket{n_i}$ are plotted for increasing statistical angle $\theta$.  In each panel the density distribution for ordinary bosons ($\theta=0$, blue line) is drawn for reference. The distribution at high statistical angles ($\theta=\frac{3\pi}{4}, \pi$, marked black and green) exhibit a Mott-like plateau at the fractional value $n \simeq \frac{1}{2}$ (marked by arrows). Parameters: $N=L=30$, $J/U=0.5$, $V/U=0.01$.
}
  \label{fig:mott}
\end{figure}

We have computed the density distribution for a system with additional harmonic confinement, at the same fixed parameter point as in the procedure outlined above. The result is plotted in Fig. 4. For low statistical phase angles $\theta$ (and thus a superfluid state), the density distribution displays a smooth, quadratic profile centered around the minimum of the trap.
Beyond the critical phase angle, $\theta > \theta_c$, Mott-like plateaus appear in $\braket{n_i}$. Surprisingly, in addition to the integer density Mott plateau at $n=1$, a new plateau emerges at the fractional value $n \simeq \frac{1}{2}$. Fractional Mott plateaus persist also in other parameter regimes, when $J/U$ and the trapping potential $V$ are varied, and seem to form a universal feature of 1D anyonic gases at large exchange phases $\theta$, which will be subject of a future study.\\
Recent progress on the experimental side have made direct imaging of density distributions possible, using ``quantum gas microscopy''\cite{Microscope1, Microscope2}. A few weeks ago, the plateau structure of the Mott insulator was directly observed for the first time\cite{Motti} at the single atom level.
This new technique opens up the possibility to directly demonstrate fractional Mott plateaus and to image statistically induced phase transitions \emph{in situ}.

In summary, we have shown how to realize fractional statistics in 1D optical lattices, using bosons in a realistic experimental environment.
The experiment we propose features the full control and tuneability of the particles' exchange statistics -- paving the way to the first \emph{statistically induced} quantum phase transition.

\subsection*{Methods}

\subsubsection*{Realizing conditional-hopping bosons}

In this section we discuss how to realize the conditional-hopping Hamiltonian equation (\ref{bosons}), using four different and independent $\Lambda$-transitions\cite{Scully}.
In general, we assume a deep optical lattice potential, giving rise to a negligible bare kinetic tunneling amplitude  $J_{kin}$.
Let us focus for a moment on one $\Lambda$-scheme, where
two ground state levels $|a\rangle$ and $|b\rangle$ are coupled through
a Raman process via an excited state $|e\rangle$.
In our case, $|a\rangle$ and $|b\rangle$ correspond to the wavefunctions
of atoms at distance $d$ localized in neighbouring sites of
the tilted lattice $V(x)$ ($d$ being the lattice constant), while $|e\rangle$ experiences another potential
$V'(x)$ (for a brief discussion of two realistic experimental possibilites to realize $V'(x)$, see the main text).
The levels have energies $E_{i}=\hbar\omega_{i}$,
$i=a,b,e$ and the transition between $a (b)$ and $e$ is driven
by an external radiation field with frequency $\omega_{e}-\omega_{a(b)}-\delta$,
with detuning $\delta$. Note that in our scheme the energies $E_{i}$ depend on $\Delta$ and $U$. The Hamiltonian for the three-level $\Lambda$-system
reads \begin{equation}
H=\sum_{i=a,b,e}\hbar\omega_{i}|i\rangle\langle i|+\frac{\hbar}{2}\left(\gamma_{a}|e\rangle\langle a|+\gamma_{b}|e\rangle\langle b|+\mathrm{h.c.}\right)\label{eq:Ham lasers}\end{equation}
where $\gamma_{a(b)}=\Omega_{a(b)}^{e}W_{a(b)}^{e}e^{-i(\omega_{e}-\omega_{a(b)}-\delta)t}$.
Here, the quantity $\Omega_{a(b)}^{e}$ is the Rabi frequency for the transition
$a(b)\to e$ with the atom centered in the same position. However,
since ground and excited states feel different lattices, the $x$
components of the Wannier functions $w(x)$ are slightly displaced.
Thus, the off-diagonal elements in equation (\ref{eq:Ham lasers}) contain
the integrals\cite{jaksch03} \begin{eqnarray}
W_{a}^{e} & = & e^{ik_{a}x_a}\int w_{e}^{*}(x+x_e)e^{ik_{a}x}w_{a}(x) dx \,,\nonumber \\
W_{b}^{e} & = & e^{ik_{b} (x_a +d)}\int w_{e}^{*}(x+x_e)e^{ik_{b}x}w_{b}(x+d)dx \,,\label{eq:Ws}\end{eqnarray}
where $k_{a(b)}$ is the $x$-component of the momentum of the driving
radiation field, with modulus $|\boldsymbol{\mathrm{k}}_{a(b)}|=(\omega_{e}-\omega_{a(b)}-\delta)/c$. $x_a$ is the atom position in the left well, while $x_e$ refers to the position in the excited state.
Note that the integrals in (\ref{eq:Ws}) are in general non-zero,
as the two Wannier functions belong to different lattices.
The quantities $\gamma_{a(b)}$ are complex numbers, whose modulus
and phase can be freely tuned by choosing the appropriate intensity, polarization
and direction of the driving fields. For sufficiently large detunings $\delta > |\gamma_{a(b)}|$, the
level $|e\rangle$ can be adiabatically eliminated and in the
rotating wave approximation the effective Hamiltonian in the subspace
$\{|a\rangle,|b\rangle\}$ reads
\begin{equation}
H_{\textrm{eff}}=-\frac{\hbar}{4\delta}\left(\begin{array}{cc}
|\tilde{\gamma}_{a}|^{2} & \tilde{\gamma}_{a}^{*}\tilde{\gamma}_{b}\\
\tilde{\gamma}_{b}^{*}\tilde{\gamma}_{a} & |\tilde{\gamma}_{b}|^{2}\end{array}\right) \,,
\end{equation}
where $\tilde{\gamma}$ is the non-rotating version of $\gamma$, i.e. without the time-dependent phase factors.
In order to realize equation (\ref{bosons}), we suggest to employ four driving fields (see labels in Fig. 1b) with different frequencies in order to avoid interferences. This situation would correspond to a maximum of two atoms per site.
If the local density truncation were set to a higher number, more external driving fields would become necessary. As all fields can be tuned independently from each other, this poses no problems besides potential budgetary considerations.\\
The couplings $J_{13}$, $J_{14}$, $J_{23}$ and $J_{24}$, between
the four different levels are then obtained in terms of the effective Rabi
frequencies, $J_{ab}=\tilde{\gamma}_{a}^{*}\tilde{\gamma}_{b}/2\delta$. \\
Our aim is to satisfy the conditions set in equations (\ref{condition1}-\ref{condition2}) in order to engineer a state-dependent phase factor. This implies $\tilde{\gamma}_{2} \equiv \tilde{\gamma}_{1}e^{i\theta}$, which can be achieved by the free tuneability of each driving fields's frequency, intensity, polarization
and direction. Furthermore this choice of parameters implies $|\tilde{\gamma}_{a}| = |\tilde{\gamma}_{b}|$, i.e. the diagonal elements of the effective Hamiltonian are now equivalent. Thus, the tilt energy $\Delta$ has vanished via adiabatic elimination, and consequently also does not appear in equation (\ref{bosons}). Indeed this offset energy between neighbouring sites is compensated by the external radiation field. As the assisted tunneling proposed in this Article is the only mechanism for hopping in the lattice, unwanted effects like Bloch oscillations do not appear in our system (the bare kinetic tunneling amplitude $J_{kin}$ is assumed to be negligible compared to all energy scales discussed here).\\
In summary, the parameters discussed here have to satisfy the following conditions:\\
First, $\delta_{\textrm{linewidth}} \ll \Delta, U$, so that the external driving fields resolve the different levels of the ground state manifold. Second, large detunings $\delta > |\gamma_{a(b)}|$ are required for a short-lived excited state and the validity of the adiabatic elimination. Third, $\Delta$ and $U$ can be in the same frequency regime (a few kHz), but should not be identical (their difference should be $\gg \delta_{\textrm{linewidth}}$). As an example, $\Delta \simeq 2$ kHz, $U \simeq$ 3 kHz, $|J_{ab}| = J \simeq$ 5 kHz  and $|\gamma_{a(b)}| \simeq 20$ kHz would be sufficient if the linewidth of the radiation field were $\delta_{\textrm{linewidth}} \simeq 50$ Hz, which is a realistic assumption for typical radio-frequency driving fields (see e.g. the works by Campbell \emph{et al.}\cite{Campbell} and McKay and DeMarco\cite{McKay}).

\subsubsection*{Fractional Jordan-Wigner mapping}
In the following we prove that anyons are isomorphic to bosons in 1D. \\
In particular we prove that the operators $a$, as defined in the non-local mapping, equation (\ref{mapping}), indeed obey the anyonic commutation relations of equation (\ref{CCR}), provided that the particles of type $b$ are bosons.

For the case $i<j$ we wish to rewrite products of anyonic operators in terms of the bosonic ones:
\begin{eqnarray}
a_i a_j^\dagger &=& b_i e^{-i \theta \sum_{i \leq k < j}n_k} b^\dagger_j \nonumber \\
&=&  e^{-i \theta \sum_{i < k < j}n_k} b_i  b^\dagger_j e^{-i \theta n_i}   \nonumber \, , \\
f(\theta) a_j^\dagger a_i  &=& e^{-i \theta \sum_{i < k < j}n_k} e^{-i \theta n_i} f(\theta) b^\dagger_j b_i \nonumber \\
&=& e^{-i \theta \sum_{i < k < j}n_k} e^{-i \theta (n_i+1)} b^\dagger_j b_i  \, .
\end{eqnarray}
Here we have defined $f(\theta) \equiv e^{i \theta \mathrm{sgn}(i-j)}$ and used that $f(\theta) = e^{-i \theta}$ since $i<j$ was assumed.
We can now evaluate the LHS of equation (\ref{CCR}):
\begin{eqnarray}
a_i a_j^\dagger - f(\theta) a_j^\dagger a_i &=& e^{-i \theta \sum_{i < k < j}n_k} (b_i  b^\dagger_j e^{-i \theta n_i} - e^{-i \theta (n_i+1)} b^\dagger_j b_i)  \nonumber \\
 &=& e^{-i \theta \sum_{i < k < j}n_k} e^{-i \theta (n_i+1)} [b_i,b^\dagger_j] \nonumber \\
 &=& 0
\end{eqnarray}
Thus the anyonic commutation relations have been proven for the case $i<j$. The proof for the case $i>j$ is very similar.
For the case $i=j$ one just has to note that $a_i^\dagger a_i = b_i^\dagger b_i$ and $f(\theta)=1$.
Note that the resulting conditional-hopping bosonic Hamiltonian, equation (\ref{bosons}), resembles the exactly solvable model of q-bosons\cite{qbosons}. However, this model is not equivalent to our model.

\begin{figure} []
  \centering
  \includegraphics[width=1.0\linewidth]{Fig5.pdf}
  \vspace{-0.5cm}
  \caption{
  \textbf{Mean-field solution for the Mott-superfluid transition.}  The transition lines are displayed for several values of $\theta$. For comparison
we display also the data for the bosonic case $\theta=0$ (blue curve).
}
  \label{fig:lobes MF}
\end{figure}

\subsubsection*{Mean-field calculation}

The conditional-hopping bosonic Hamiltonian is given by equation (\ref{bosons}). In addition, we now include also a chemical potential term.
Expressing energies in units of $U$ (we fix $U\equiv1$), we have

\begin{eqnarray}
H & = & \sum_{j}\left[\frac{1}{2}n_{j}(n_{j}-1)-\mu n_{j}-J(c_{j}^{\dagger}b_{j+1}+b_{j+1}^{\dagger}c_{j})\right]\label{eq:anyon Hubbard}\end{eqnarray}
where for convenience we have defined $c_{j}=e^{-i\theta n_{j}}b_{j}$.

In absence of hopping $J=0$, all the sites are independent and the
ground state is of the Gutzwiller type\[
|\Psi_{0}\rangle=|\psi\rangle^{\otimes L},\qquad|\psi\rangle=\sum_{\nu=0}^{\infty}c_{\nu}\frac{(b^{\dagger})^{\nu}}{\sqrt{\nu!}}|0\rangle,\]
where $\nu=N/L$ is the filling factor. The local energies are $\epsilon(\nu)=\frac{1}{2}\nu(\nu-1)-\mu\nu$,
while the local gaps for adding or subtracting one boson are, respectively
\[
\epsilon(\nu+1)-\epsilon(\nu)=\nu-\mu,\qquad\epsilon(\nu-1)-\epsilon(\nu)=-(\nu-1)+\mu.\]
So, the ground state in every site has $\nu$ particles in the interval $\mu_{-}^{(\nu)}<\mu<\mu_{+}^{(\nu)}$
, with $\mu_{-}^{(\nu)}=\nu-1$ and $\mu_{+}^{(\nu)}=\nu$. The gap
in the whole system at a given number of particles is given by $\Delta=1$,
obtained by removing an atom at some site an putting it in another
one already occupied.

Here the mean field (MF) is obtained by decoupling the hopping term
as $c_{j}^{\dagger}b_{j+1}\approx-\alpha_{2}^{*}\alpha_{1}+\alpha_{2}^{*}b_{j+1}+\alpha_{1}c_{j}^{\dagger}$,
where the order parameters are $\alpha_{1}=\langle b_{j}\rangle$
and $\alpha_{2}=\langle c_{j}\rangle$. Accordingly, the Hamiltonian
(\ref{eq:anyon Hubbard}) in MF becomes

\[
H=\sum_{j}H_{j}+LJ(\alpha_{2}^{*}\alpha_{1}+\alpha_{1}^{*}\alpha_{2})\]
with\[
H_{j}=\frac{1}{2}n_{j}(n_{j}-1)-\mu n_{j}-J\sum_{j}(\alpha_{2}b_{j}^{\dagger}+\alpha_{2}^{*}b_{j}+\alpha_{1}c_{j}^{\dagger}+\alpha_{1}^{*}c_{j}).\]
The two parameters $\alpha_{1}$ and $\alpha_{2}$ are not completely
independent as as they are both vanishing or non vanishing at the
same time. Of course there is a trivial solution corresponding to
$\alpha_{l}=0$, $l=1,2$, which corresponds to the Mott insulating phase. The
occurrence of $\alpha_{l}\neq0$ signals the instability towards superfluid
correlations. On inhomogeneous lattices a further situation can in
principle occur, where $\alpha_{l}\neq0$ only on a fraction of the
lattice sites. The self-consistent relation defines a map $\alpha_{l}=\Lambda_{ll'}\alpha_{l'}$,
obtained by linearizing about the solution $\alpha_{l}=0$. The instability
of the trivial solution sets in when the maximal eigenvalue of $\Lambda$
is greater than 1. Close to the trivial point, it holds $|\alpha_{l}|\ll1$,
hence the kinetic term can be treated perturbatively. Up to first
perturbative order, the wavefunction can be written as $|\psi\rangle=|\psi^{(0)}\rangle+|\psi^{(1)}\rangle$,
where $|\psi^{(0)}\rangle=|\nu\rangle$ and \begin{eqnarray*}
|\psi^{(1)}\rangle & = & -J\sum_{\nu'}\frac{\langle\nu'|\alpha_{2}b_{j}^{\dagger}+\alpha_{2}^{*}b_{j}+\alpha_{1}c_{j}^{\dagger}+\alpha_{1}^{*}c_{j}|\nu\rangle}{\epsilon(\nu)-\epsilon(\nu^{\prime})}|\nu^{\prime}\rangle\\
 & = & J\frac{\sqrt{\nu}(\alpha_{2}^{*}+\alpha_{1}^{*}e^{-i\theta(\nu-1)})}{\mu-\nu+1}|\nu-1\rangle\\
 &  & +J\frac{\sqrt{\nu+1}(\alpha_{2}+\alpha_{1}e^{i\theta\nu})}{\nu-\mu}|\nu+1\rangle\end{eqnarray*}

Hence using the self-consistency relations $\alpha_{1}=\langle\psi|b_{j}|\psi\rangle$
and $\alpha_{2}=\langle\psi|c_{j}|\psi\rangle$, \begin{eqnarray*}
\alpha_{1}/J & = & \frac{\nu(\alpha_{2}+\alpha_{1}e^{i\theta(\nu-1)})}{\mu-\nu+1}+\frac{(\nu+1)(\alpha_{2}+\alpha_{1}e^{i\theta\nu})}{\nu-\mu}\\
\alpha_{2}/J & = & \frac{\nu(\alpha_{2}e^{-i\theta(\nu-1)}+\alpha_{1})}{\mu-\nu+1}+\frac{(\nu+1)(\alpha_{2}e^{-i\theta\nu}+\alpha_{1})}{\nu-\mu}\end{eqnarray*}
The matrix $\Lambda$ is then \[
\Lambda=J\left(\begin{array}{cc}
f(\theta) & A\\
A & f(-\theta)\end{array}\right),\qquad f(\theta)=e^{i\theta\nu}\left[A+(e^{-i\theta}-1)B\right]\]
with \[
A=\frac{\mu+1}{(\mu-[\mu])([\mu]-\mu+1)},\qquad B=\frac{[\mu]+1}{\mu-[\mu]},\]
since every lobe is labelled by $\nu=[\mu]+1$. The eigenvalues of
$\Lambda$ are given by

\[
\lambda_{\pm}=\frac{J}{2}\left[f(\theta)+f(-\theta)\pm\sqrt{4A^{2}+(f(\theta)-f(-\theta))^{2}}\right].\]
The condition $\max\{\lambda_{+},\lambda_{-}\}=1$ signals the onset
of instability of the trivial solution and hence establishes the critical coupling
$J_{crit}$ along the Mott-superfluid transition line. In Fig. \ref{fig:lobes MF} the phase diagrams for several values of $\theta$ are shown. The Mott lobes expand for non-zero statistical angles $\theta$, a fact central for designing statistically induced phase transitions.