text:
# 1 Introduction

## 1 Introduction

Suppose that we have a Hamiltonian vector field defined on a Poisson space. A natural problem is to find its equilibrium points and to check whether or not they are stable. Stability is important in mathematical modelling, as it gives an indication which behaviour exhibited by a (mathematical) modelling system is a reliable representative of behaviour in the corresponding (real) modelled system. In the present paper we adopt the following

Definition. Let $`𝐗_H`$ be the Hamiltonian vector field corresponding to the function $`H`$. The equilibrium point $`x`$ of $`𝐗_H`$ is stable if for any neighbourhood $`U`$ of $`x`$, there exists a neighbourhood $`V`$ of $`x`$ such that $`\varphi (0)V`$ and $`\dot{\varphi }(t)=𝐗_H(\varphi (t))`$ implies that $`\varphi (t)Ut0`$; otherwise $`x`$ is an unstable equilibrium of $`𝐗_H`$.

The aim of the present work is to perform a stability analysis for certain integrable Euler equations associated with the group $`SO(n)`$, focusing mainly on the case $`n=4`$, and in addition to examine the heteroclinic orbits. These equations represent a particularly simple special case of the integrable Hamiltonian systems introduced by Manakov in . In the $`n=3`$ case they reduce to the classical Euler equations for the angular momentum of a free rigid body in the moving frame. As explained in several mechanics textbooks, the qualitative behaviour of the solutions of the classical Euler equations is easily visualised in terms of their phase portrait, see for instance the picture on the cover of the book . This system lives on a coadjoint orbit $`S^2`$ of $`SO(3)`$, and it has the interesting feature that the unstable equilibria are connected by heteroclinic orbits that are given by great circles on the sphere $`S^2`$. Let us recall that a heteroclinic orbit in general consists of the points of a nontrivial integral curve of a dynamical system and equilibrium points. Since the great circles on $`S^2`$ are the orbits of the 1-parameter subgroups of $`SO(3)`$, we shall in the $`SO(n)`$ case enquire about the existence of heteroclinic orbits that are orbits of 1-parameter subgroups of $`SO(n)`$; we consider this to give rise to the most interesting results of the paper.

This work is intended as a step towards a stability analysis of the full set of the Manakov systems  which, in addition to the special case studied here, contains for example the $`n`$-dimensional rigid body of , and has many interesting Lie-algebraic generalizations . In fact, most of our results are not difficult to extend. The result for which generalization presents problems is the one described in Section 3.

For convenient reference later, we now recall some standard facts about the linearisation of a Hamiltonian dynamical system and its use in the stability analysis of the original system. Let $`M`$ be a Poisson space. Let $`HC^{\mathrm{}}(M)`$ and suppose that $`x`$ is an equilibrium point of the Hamiltonian vector field $`𝐗_H`$. The linearisation at $`x`$ is a flow in $`T_xM`$ given by

$$\dot{v}=(_V𝐗_H)(x),$$
(1.1)
with $`V`$ any vector field such that $`V(x)=v`$. By choosing any system of local coordinates in a neighbourhood of $`x`$, this becomes a system of the form $`\dot{v}=𝐋v`$, with $`𝐋`$ a square matrix of the same size as the dimension of $`M`$. For a Hamiltonian system the eigenvalues of the linearisation at $`x`$ come in groups of four, in the sense that if $`\lambda `$ is an eigenvalue of $`𝐋`$ then $`\lambda `$ is an eigenvalue of $`𝐋`$ and so is the complex conjugate $`\overline{\lambda }`$. The following statements are well known, see .

1. $`x`$ is unstable if the linearisation at $`x`$ of the system $`\dot{\varphi }=𝐗_H(\varphi )`$ has an eigenvalue with a positive real part. If no eigenvalues of the linearisation have positive real part, then all eigenvalues have to be imaginary; in this case $`x`$ may or may not be stable.

2. $`x`$ is stable if there exists $`fC^{\mathrm{}}(M)`$ for which $`\{f,H\}=0`$ in a neighbourhood of $`x`$ and

$$(i)df(x)=0,(ii)d^2f(x)\text{ is definite.}$$
(1.2)
If the rank of the Poisson bracket is constant in some neighbourhood of $`x`$, then it is sufficient that properties $`(i)`$ and $`(ii)`$ in (1.2) be satisfied with respect to vectors tangent to the symplectic leaf through $`x`$. This is discussed for example in . In this paper we consider such a “regular situation” since the phase spaces of our interest will be generic coadjoint orbits of the Lie group $`SO(n)`$. We will assume all entities appearing in the definition of the systems studied to be generic, since this would be a reasonable assumption in a physical context and it also simplifies the problem.

The organization of the paper and of our results is as follows. The next section contains the definition of the Hamiltonian systems of interest associated with $`SO(n)`$ together with a description of their equilibrium points (Proposition 1). In section 3 we present a complete analysis of the stability of the equilibrium points in the $`n=4`$ case. The outcome of our study is given by Proposition 2. In section 4 we describe a necessary condition (Proposition 3) for the possibility to construct heteroclinic orbits by means of $`1`$-parameter subgroups for Hamiltonian systems living on a coadjoint orbit, and concretely construct such heteroclinic orbits for the systems associated with $`SO(4)`$. In section 5 the main features of the stability analysis are outlined for any $`n`$. In particular, the construction of the heteroclinic orbits is generalized to the $`SO(n)`$ case (see Proposition 4). Section 6 contains a brief summary of the results and some open problems.

## 2 A family of integrable Euler equations for $`SO(n)`$

We define below the Hamiltonian systems to be studied and describe their equilibrium points. As explained at the end of the section, these systems correspond to a special case of the integrable Euler equations introduced in .

Consider the Lie algebra $`so(n)`$ of the real orthogonal group $`SO(n)`$. An element of $`so(n)`$ is an $`n\times n`$ antisymmetric real matrix. The Lie-Poisson bracket of functions on $`so(n)^{}`$ is given by

$$\{\varphi ,\psi \}(\alpha )=<\alpha ,[d_\alpha \varphi ,d_\alpha \psi ]>\alpha so(n)^{},$$
(2.1)
where $`d_\alpha \varphi so(n)`$ is defined by

$$<\beta ,d_\alpha \varphi >=\frac{d}{dt}|_{t=0}\varphi (\alpha +t\beta )\beta so(n)^{},$$
(2.2)
and $`d_\alpha \psi `$ similarly. The symplectic leaves in $`so(n)^{}`$ are the coadjoint orbits of $`SO(n)`$ in $`so(n)^{}`$. It will be convenient to identify $`so(n)^{}`$ with $`so(n)`$ with the aid of a multiple of the standard trace form for $`n\times n`$ matrices, so that $`<\beta ,X>:=\frac{1}{2}\mathrm{tr}(\beta X)`$.

Let us define the Cartan subalgebra $`𝐡`$ in $`so(n)`$ to be the set of all matrices $`x`$ of the form

$$x=\underset{k=1}{\overset{m}{}}x_ke_{kk}\mathrm{i}\sigma _2\text{if }n=2m,\text{or}x=\left(\begin{array}{cc}_{k=1}^mx_ke_{kk}\mathrm{i}\sigma _2& 0\\ 0& 0\end{array}\right)\text{if }n=2m+1,$$
(2.3)
where $`m`$ is any positive integer. Here $`e_{ij}`$ is the $`m\times m`$ matrix having 1 for the term in the $`i`$th row and in the $`j`$th column and all other terms zero. We use the Pauli matrices

$$\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _2=\left(\begin{array}{cc}0& \mathrm{i}\\ \mathrm{i}& 0\end{array}\right),\sigma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\sigma _0:=\mathrm{𝟏}=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(2.4)
An element $`x`$ of $`𝐡`$ is generic if $`x_k0k`$ and $`x_k^2x_l^2`$ if $`kl`$. Using the identification of $`so(n)^{}`$ with $`so(n)`$, a generic symplectic leaf can be written as

$$𝒪_x=\{gxg^1|gSO(n)\}$$
(2.5)
with $`x`$ a generic element in $`𝐡`$. The isotropy subalgebra of $`x`$ in $`so(n)`$ consists of matrices of the same form as $`x`$, i.e. given by the same formula as (2.3) with different values of $`x_i`$. The isotropy subgroup $`SO(n)_x`$ is the exponential of this algebra.

In this paper we are interested in Hamiltonian systems $`(𝒪_x,\{,\},H)`$ on generic coadjoint orbits, where $`H`$ has the form

$$H(\mu ):=\frac{1}{2}\mathrm{tr}(J\mu ^2),\mu 𝒪_x,$$
(2.6)
with some constant matrix $`J=\mathrm{diag}(J_1,\mathrm{},J_n)`$. We assume that $`J_i^2J_j^2`$ if $`ij`$. The generalized Euler equation defined by the Hamiltonian vector field $`𝐗_H`$ can be written as follows:

$$\dot{\mu }=[J,\mu ^2].$$
(2.7)
An equilibrium point on $`𝒪_x`$, for $`x`$ given by (2.3), is a point $`gxg^1`$ such that

$$0=[J,gx^2g^1].$$
(2.8)
Let $`p`$ be an element of the permutation group $`S_n`$ (the Weyl group of $`sl(n)`$), and introduce the permutation matrix $`\overline{p}O(n)`$ by

$$\overline{p}_{ij}=\delta _{i,p(j)}(i,j=1,\mathrm{},n).$$
(2.9)
For any diagonal matrix $`D=\mathrm{d}iag(d_1,\mathrm{},d_n)`$, one has

$$p(D):=\mathrm{d}iag(d_{p^1(1)},\mathrm{},d_{p^1(n)})=\overline{p}D\overline{p}^1,$$
(2.10)
and the parity of $`p`$ satisfies $`\mathrm{s}gn(p)=\mathrm{d}et(\overline{p})`$. Since $`J`$ and $`x^2`$ are diagonal matrices, $`\overline{p}x\overline{p}^1`$ is clearly an equilibrium point whenever it belongs to $`𝒪_x`$. This holds obviously for the even permutations. If $`n=(2m+1)`$ is odd, then $`\overline{p}x\overline{p}^1𝒪_x`$ for any $`pS_n`$, since in this case

$$\overline{p}D\overline{p}^1=\widehat{p}D\widehat{p}^1\text{with}\widehat{p}:=\mathrm{s}gn(p)\overline{p}SO(2m+1).$$
(2.11)
We can prove that the equilibrium points associated in this manner with the permutations exhaust all the equilibria on $`𝒪_x`$.

Proposition 1. The set of equilibrium points on a generic orbit $`𝒪_x`$, for $`x`$ of the form given in (2.3), consists of the matrices $`\overline{p}x\overline{p}^1`$, where $`pS_n`$ is an even permutation if $`n`$ is even, and $`pS_n`$ is an arbitrary permutation if $`n`$ is odd. The equilibrium points associated with different permutations are different.

Proof. Let us consider the set

$$E_x:=\{gxg^1|[J,gx^2g^1]=0,gO(n)\}.$$
(2.12)
Since $`J`$ is a regular diagonal matrix by assumption, $`gx^2g^1`$ must be a diagonal matrix whose entries are obtained by permuting the entries of the diagonal matrix $`x^2`$. We can choose a set of elements of $`S_n`$, say $`\{p_i\}_{i=1}^N`$, for which the matrices $`p_i(x^2)`$ are distinct from each other for $`ij`$ and they contain all matrices that are obtained by permuting the diagonal entries of $`x^2`$. Note that $`N=\frac{n!}{2^m}`$ for $`n=2m`$ or $`n=(2m+1)`$, and the $`p_i`$ are a set of representatives for the coset space $`S_n/S_n^{x^2}`$, where

$$S_n^{x^2}=\{pS_n|p(x^2)=x^2\}.$$
(2.13)
For $`n=2m`$ or $`n=(2m+1)`$, the group $`S_n^{x^2}`$ is generated by the elements

$$\tau ^{1,2},\tau ^{3,4},\mathrm{},\tau ^{2m1,2m},$$
(2.14)
where $`\tau ^{k,l}S_n`$ denotes the transposition that exchanges $`k`$ with $`l`$.

Since any $`gO(n)`$ that appears in (2.12) satisfies $`gx^2g^1=\overline{p}_ix^2\overline{p}_i^1`$ with some $`1iN`$, it follows that the most general such $`g`$ can be written as

$$g=\overline{p}_i\gamma \text{with some}\gamma O(n)_{x^2},$$
(2.15)
$$O(n)_{x^2}:=\{\gamma O(n)|\gamma x^2\gamma ^1=x^2\}.$$
(2.16)
The isotropy group $`O(n)_{x^2}`$ consists of block-diagonal matrices with arbitrary elements of $`O(2)`$ in the $`2\times 2`$ blocks. It is useful to consider also

$$O(n)_x:=\{qO(n)|qxq^1=x\},$$
(2.17)
which consists of block-diagonal matrices with each $`2\times 2`$ block containing an arbitrary element of $`SO(2)`$. The point to notice is that any $`\gamma O(n)_{x^2}`$ can be uniquely written in the form111In fact, $`O(n)_x`$ is a normal subgroup of $`O(n)_{x^2}`$ and $`S_n^{x^2}`$ is the corresponding factor group.

$$\gamma =\overline{\mathrm{\Gamma }}q,\mathrm{\Gamma }S_n^{x^2},qO(n)_x.$$
(2.18)
This follows from the fact that $`O(2)/SO(2)`$ can be identified with the group generated by the transposition matrix $`\sigma _1`$. By using these observations, we see that any element $`gxg^1E_x`$ has the form

$$gxg^1=\overline{p}_i\overline{\mathrm{\Gamma }}x(\overline{p}_i\overline{\mathrm{\Gamma }})^1(1iN,\mathrm{\Gamma }S_n^{x^2}).$$
(2.19)
As $`\overline{p}=\overline{p}_i\overline{\mathrm{\Gamma }}`$ for $`p=p_i\mathrm{\Gamma }`$, this implies that all elements of $`E_x`$ are given by $`\overline{p}x\overline{p}^1`$ with some $`pS_n`$. It is clear from the definitions that any $`pS_n`$ can be decomposed as $`p=p_i\mathrm{\Gamma }`$ with a unique $`p_i`$ and a unique element of $`S_n^{x^2}`$, and one can check directly that different permutations are associated with different points of $`E_x`$.

If $`n=(2m+1)`$, the statement of the proposition follows immediately from the above-established results (the coadjoint orbits of $`O(2m+1)`$ and $`SO(2m+1)`$ coincide). The proof is completed by noting that in the $`n=2m`$ case only those elements $`\overline{p}x\overline{p}^1`$ lie on the coadjoint orbit of $`SO(2m)`$ through $`x`$ for which $`p`$ is an even permutation. This fact can be verified, for example, by performing an analogous analysis as above in the case for which $`g`$ in (2.12) is restricted to $`SO(n)`$ from the beginning.

Remark 1. Suppose that we study the nature of an equilibrium point $`\overline{p}x\overline{p}^1𝒪_x`$. We may then choose a different basis in which this point is represented by the matrix of $`x`$, and $`J`$ is replaced by the matrix $`\overline{p}^1J\overline{p}`$. We may thus assume without loss of generality that the equilibrium point of interest is always represented by the same matrix $`x𝐡`$ in (2.3).

Remark 2. It follows from (2.7) that in our case the linearised system at the point $`x`$ is the flow in $`T_x𝒪_x`$ defined by

$$\dot{v}=[J,vx+xv].$$
(2.20)
Remark 3. In the terminology of generalized rigid bodies  the quantity $`\mu `$ in the Euler equation (2.7) is the angular momentum relative to the body. Correspondingly, the inverse of the moment of inertia operator maps $`\mu `$ to the angular velocity $`\omega `$ relative to the body according to $`\mu \omega =(J\mu +\mu J)`$. Indeed, then (2.7) takes the classical form $`\dot{\mu }=[\mu ,\omega ]`$. This is a special case of the integrable rigid body systems introduced in  by the relation $`\mu _{ij}=\frac{a_ia_j}{b_ib_j}\omega _{ij}`$ with arbitrary constants $`a_i`$, $`b_i`$. The case (2.7) arises by setting $`b_i=a_i^2`$ with $`a_i=J_i`$, while the $`n`$-dimensional rigid body of  is obtained by setting $`a_i=b_i^2`$.

## 3 Stability analysis in the $`n=4`$ case

The Lie algebra $`so(4)`$ is the same as the direct sum $`so(3)so(3)`$. This can be seen by identifying $`so(3)`$ with $`su(2)`$ and then finding two commuting copies of $`su(2)`$ in $`so(4)`$. In terms of the Pauli matrices, we have $`su(2)=span\{\mathrm{i}\sigma _1,\mathrm{i}\sigma _2,\mathrm{i}\sigma _3\}`$, and two commuting $`su(2)`$ subalgebras that together span $`so(4)`$ are

$$su(2)\mathrm{span}\{\sigma _1\mathrm{i}\sigma _2,\mathrm{i}\sigma _2\sigma _0,\sigma _3\mathrm{i}\sigma _2\}\mathrm{span}\{i\sigma _2\sigma _1,\sigma _0\mathrm{i}\sigma _2,\mathrm{i}\sigma _2\sigma _3\}.$$
(3.1)
In the coordinates $`(l,m)`$ on $`so(4)=so(3)so(3)`$ given by

$$\mu =\left[l_1\sigma _1\mathrm{i}\sigma _2+l_2\mathrm{i}\sigma _2\sigma _0+l_3\sigma _3\mathrm{i}\sigma _2\right]\left[m_1\mathrm{i}\sigma _2\sigma _1+m_2\mathrm{i}\sigma _2\sigma _3+m_3\sigma _0\mathrm{i}\sigma _2\right],$$
(3.2)
the Poisson bracket is

$$\{l_i,l_j\}=ϵ_{ijk}l_k,\{l_i,m_j\}=0,\{m_i,m_j\}=ϵ_{ijk}m_k,$$
(3.3)
and $`|l|^2`$ and $`|m|^2`$ are Casimir functions. The rigid body Hamiltonian $`H`$ and an independent commuting integral $`K`$ are now given by

$$H(l,m)=l^T\mathrm{\Lambda }m,K(l,m)=\frac{1}{2}l^T\mathrm{\Lambda }^2l+\frac{1}{2}m^T\mathrm{\Lambda }^2ml^T\mathrm{\Theta }m,$$
(3.4)
where $`\mathrm{\Lambda }`$ and $`\mathrm{\Theta }`$ are constant diagonal matrices

$$\mathrm{\Lambda }=\mathrm{diag}(\mathrm{\Lambda }_1,\mathrm{\Lambda }_2,\mathrm{\Lambda }_3),\mathrm{\Theta }=\mathrm{diag}(\mathrm{\Lambda }_2\mathrm{\Lambda }_3,\mathrm{\Lambda }_1\mathrm{\Lambda }_3,\mathrm{\Lambda }_1\mathrm{\Lambda }_2).$$
(3.5)
The equations of motion corresponding to $`H`$ are

$$\dot{l}=(\mathrm{\Lambda }m)l\dot{m}=(\mathrm{\Lambda }l)m.$$
(3.6)
We apply the usual identification of $`so(3)`$ with $`𝐑^3`$ equipped with the vector-product, denoted by $``$. The formulae in (3.4) may be recovered by a technique due to Manakov : define $`=\lambda J+\mu `$, then the set of coefficients of $`\lambda `$ amongst all traces of powers of $``$ forms a commuting family. Here $`K(\mathrm{t}r^4)|_{\lambda ^2}`$, $`H(\mathrm{t}r^3)|_{\lambda ^1}`$ up to Casimirs and (from (2.6)) $`\mathrm{\Lambda }`$ is related to $`J`$ by $`\mathrm{\Lambda }_1=J_1+J_2+J_3J_4`$, $`\mathrm{\Lambda }_2=J_1+J_2J_3+J_4`$, $`\mathrm{\Lambda }_3=J_1J_2+J_3+J_4`$. From now on we make the genericity assumption that

$$\mathrm{\Lambda }_i{}_{}{}^{2}\mathrm{\Lambda }_j{}_{}{}^{2}\text{if}ij\text{and}\mathrm{\Lambda }_i0i.$$
(3.7)
The first part of these conditions follows from the assumption that $`J_p^2J_q^2`$ for $`pq`$.

Let $`e_k`$ ($`k=1,2,3`$) denote the standard basis of $`𝐑^3`$. The equilibrium points of (3.6) that lie on generic coadjoint orbits are in fact given by

$$(l,m)=b(ae_k,e_k)\text{with }𝐑a,b0,k=1,2,3.$$
(3.8)
We next study the stability of an equilibrium point of the form

$$(l,m)=(ae_3,e_3)\text{ with }e_3=(0,0,1)^T\text{and}𝐑a0,$$
(3.9)
and then the general case (3.8) will be reduced to this one.

The elements of $`T_{(ae_3,e_3)}𝒪_{(ae_3,e_3)}`$ can be parametrized as $`(a\xi e_3,\eta e_3)`$ with $`\xi =\xi _1e_1+\xi _2e_2`$ and $`\eta =\eta _1e_1+\eta _2e_2`$. By putting $`v=(\xi _1,\eta _1,\xi _2,\eta _2)^T`$, the linearised system at $`(ae_3,e_3)`$ is given explicitly by $`\dot{v}=Lv`$ with

$$L=\left(\begin{array}{cccc}0& 0& \mathrm{\Lambda }_3& \mathrm{\Lambda }_1\\ 0& 0& a\mathrm{\Lambda }_1& a\mathrm{\Lambda }_3\\ \mathrm{\Lambda }_3& \mathrm{\Lambda }_2& 0& 0\\ a\mathrm{\Lambda }_2& a\mathrm{\Lambda }_3& 0& 0\end{array}\right).$$
(3.10)
The eigenvalues $`\zeta `$ of $`L`$ satisfy

$$\zeta ^4+\left[(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2\right]\zeta ^2+a^2(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2})(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2})=\mathrm{d}et(\zeta IL)=0.$$
(3.11)
Let $`𝒟`$ be defined by

$`𝒟:=\left[(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2\right]^24a^2(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2})(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2})`$
$`=\mathrm{\Lambda }_3{}_{}{}^{2}[(a^21)^2\mathrm{\Lambda }_3{}_{}{}^{2}+4a^2(\mathrm{\Lambda }_1{}_{}{}^{2}+\mathrm{\Lambda }_2{}_{}{}^{2})+4a(a^2+1)\mathrm{\Lambda }_1\mathrm{\Lambda }_2].`$ (3.12)
Stability of the equilibrium point $`(ae_3,e_3)`$ requires all roots of (3.8) to be imaginary. Hence all three of the following conditions must be fulfilled:

$`(i)(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2})(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2})>0,`$
$`(ii)(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2>0,`$ (3.13)
$`(iii)a)𝒟>0,\text{ or }b)𝒟=0.`$
If any one of the conditions of (3.13) is not satisfied then (3.11) has roots of the form $`\zeta =\pm \alpha \pm i\beta `$, with $`\alpha 0`$, and the equilibrium point (3.9) is unstable.

Suppose that $`𝒟=0`$. Every neighbourhood of $`(ae_3,e_3)`$ contains points of the form $`((a\pm ϵ)e_3,e_3)`$ with $`ϵ>0`$. As $`𝒟<0`$ at one of these two points, it follows that there are unstable equilibrium points arbitrarily close to $`(ae_3,e_3)`$ and hence $`(ae_3,e_3)`$ is unstable. We have instability then if (i), (ii), (iiib) of (3.13) are satisfied despite all the eigenvalues of the linearised system being pure imaginary.

We shall prove stability in the case (i), (ii), (iiia) of (3.13) by exhibiting a constant of motion for which (1.2) holds. As a preparation let us introduce

$$F:=\mathrm{\Lambda }_1\mathrm{\Lambda }_2H+\mathrm{\Lambda }_3K$$
(3.14)
and denote by $`\stackrel{~}{H}`$ and $`\stackrel{~}{F}`$ the restrictions of $`H`$ and $`F`$ to the orbit through the equilibrium point (3.9). One can check that $`d\stackrel{~}{F}=0`$ at $`(ae_3,e_3)`$ and, up to a common constant of proportionality, the Hessians of $`\stackrel{~}{H}`$ and $`\stackrel{~}{F}`$ at this critical point are found to be

$$d^2\stackrel{~}{H}\left(\begin{array}{cc}𝐇_1& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝐇_2\end{array}\right)$$
(3.15)
and

$$d^2\stackrel{~}{F}\left(\begin{array}{cc}𝐅_1& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝐅_2\end{array}\right)$$
(3.16)
where $`𝐇_1,𝐇_2,𝐅_1,𝐅_2`$ are the following $`2\times 2`$ matrices:

$$𝐇_1=\left(\begin{array}{cc}\mathrm{\Lambda }_3\mathrm{\Lambda }_2& 0\\ 0& \mathrm{\Lambda }_3+\mathrm{\Lambda }_2\end{array}\right),𝐇_2=\left(\begin{array}{cc}\mathrm{\Lambda }_3\mathrm{\Lambda }_1& 0\\ 0& \mathrm{\Lambda }_3+\mathrm{\Lambda }_1\end{array}\right),$$
(3.17)
$$𝐅_1=(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2})\left(\begin{array}{cc}a+1& a1\\ a1& a+1\end{array}\right)\left(\begin{array}{cc}\mathrm{\Lambda }_3+\mathrm{\Lambda }_1& 0\\ 0& \mathrm{\Lambda }_3\mathrm{\Lambda }_1\end{array}\right)\left(\begin{array}{cc}a+1& a1\\ a1& a+1\end{array}\right),$$
(3.18)
$$𝐅_2=(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2})\left(\begin{array}{cc}a+1& a1\\ a1& a+1\end{array}\right)\left(\begin{array}{cc}\mathrm{\Lambda }_3+\mathrm{\Lambda }_2& 0\\ 0& \mathrm{\Lambda }_3\mathrm{\Lambda }_2\end{array}\right)\left(\begin{array}{cc}a+1& a1\\ a1& a+1\end{array}\right).$$
(3.19)
Lemma. If (i), (ii), (iiia) of (3.13) are all satisfied then the equilibrium point (3.9) is stable.

Proof. There are two cases to consider.

Case one: $`\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2}>0`$ and $`\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2}>0`$. In this case (i) clearly holds. It is obvious that $`d^2\stackrel{~}{H}`$ is either positive or negative definite at $`(ae_3,e_3)`$ and the same applies to $`d^2\stackrel{~}{F}`$. Of course it can be shown that (ii) and (iiia) also hold.

Case two: $`\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2}<0`$ and $`\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2}<0`$. In this case again (i) clearly holds. Let us additionally suppose that (ii) and (iiia) both hold. We can show that there exists a $`z𝐑`$ such that $`d^2(4z\stackrel{~}{H}+\stackrel{~}{F})`$ is definite at $`(ae_3,e_3)`$.

The details of the proof in case two are as follows. Let us write

$$d^2(4z\stackrel{~}{H}+\stackrel{~}{F})\left(\begin{array}{cc}𝐐_1& 0\\ 0& 𝐐_2\end{array}\right):=𝐐,$$
(3.20)
with the $`2\times 2`$ matrices $`𝐐_i=4z𝐇_i+𝐅_i`$. Now $`𝐐`$ is a positive or negative definite matrix if and only if

$$\mathrm{d}et𝐐_1>0,\mathrm{d}et𝐐_2>0,\text{and }\mathrm{t}r𝐐_1\mathrm{t}r𝐐_2>0.$$
(3.21)
The first and second conditions of (3.21) require

$$z^2+\left[(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2\right]z+a^2(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2})(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2})<0.$$
(3.22)
Notice that (3.22) is similar to (3.11). Now (i), (ii), (iiia) together are equivalent to (3.11) having four distinct, imaginary eigenvalues, and this is obviously equivalent to the solvability of (3.22) for $`z𝐑`$. Let us write a solution $`z`$ in the form

$$z=\frac{1}{2}\left[(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2\right]+\frac{1}{2}\beta .$$
(3.23)
Then (3.22) implies

$$\beta ^2<𝒟$$
(3.24)
and because of (i),

$$\beta <(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2.$$
(3.25)
Using (3.23) we obtain

$$\mathrm{t}r𝐐_1\mathrm{t}r𝐐_2=(8z\mathrm{\Lambda }_3+4(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2})\mathrm{\Lambda }_3(a^2+1))(8z\mathrm{\Lambda }_3+4(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2})\mathrm{\Lambda }_3(a^2+1))=16\mathrm{\Lambda }_3{}_{}{}^{2}XY,$$
(3.26)
with

$`X=(\beta \left[(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2\right]+(a^2+1)(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2}),`$
$`Y=(\beta \left[(a^2+1)\mathrm{\Lambda }_3{}_{}{}^{2}+2a\mathrm{\Lambda }_1\mathrm{\Lambda }_2\right]+(a^2+1)(\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2}).`$ (3.27)
Eq. (3.25) with the assumptions $`\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_1{}_{}{}^{2}<0`$ and $`\mathrm{\Lambda }_3{}_{}{}^{2}\mathrm{\Lambda }_2{}_{}{}^{2}<0`$ imply that $`X<0`$ and $`Y<0`$, and hence $`\mathrm{t}r𝐐_1\mathrm{t}r𝐐_2>0`$. Since all three conditions (3.21) for the definiteness of $`𝐐`$ are satisfied, $`f:=(4z\stackrel{~}{H}+\stackrel{~}{F})`$ satisfies (1.2) at the equilibrium point $`(ae_3,e_3)`$, whereby the proof is complete.

The results proven above imply the following proposition, which provides a characterization of the stability of the equilibrium points of (3.6) on generic coadjoint orbits.

Proposition 2. The equilibrium point $`b(ae_k,e_k)`$ in (3.8) is stable if and only if (i), (ii), (iiia) of (3.13) hold for the constant $`a`$ and the matrix $`\mathrm{\Lambda }`$ replaced by the matrix $`\mathrm{\Lambda }^P:=diag(\mathrm{\Lambda }_{P(1)},\mathrm{\Lambda }_{P(2)},\mathrm{\Lambda }_{P(3)})`$ where $`P`$ is an even permutation of $`(1,2,3)`$ for which $`P(3)=k`$.

The permutation part of the statement follows obviously from (3.6) after checking that the stability of $`b(ae_3,e_3)`$ is equivalent to the stability of $`(ae_3,e_3)`$. In general, equation (2.7) has the property that $`\mu (t)`$ is a solution if and only if $`\mu _b(t):=b\mu (bt)`$ is a solution for any $`b0`$. This implies the required result for $`b>0`$. The $`b=1`$ case is settled by using the facts that the matrix of the linear system (2.20) simply gets multiplied by $`1`$ under such a rescaling of the equilibrium point, while the conserved quantities $`H`$ and $`K`$ in (3.4), and thus also their second variations, remain unchanged.

## 4 Heteroclinic orbits from 1-parameter subgroups

Consider two equilibrium points, $`x_0`$ and $`x_1`$, of a smooth Hamiltonian vector field $`𝐗_H`$ on a coadjoint orbit $`𝒪`$ of a compact Lie group $`G`$ with Lie algebra $`𝐠`$. Let us look for a 1-parameter subgroup of $`G`$ that generates a heteroclinic orbit of $`𝐗_H`$ connecting these equilibria. For $`Y𝐠`$, define

$$\gamma (s)=Ad_{\mathrm{exp}(sY)}^{}x_0.$$
(4.1)
Then our first requirement is that $`\gamma (s_1)=x_1`$ for some $`s_1>0`$. Setting $`s_0:=0`$, our second requirement is that the curve $`\gamma :(s_0,s_1)𝒪`$ yields an integral curve of $`𝐗_H`$ by a suitable reparametrization. In other words, there should exist an increasing diffeomorphism $`T:(s_0,s_1)(\mathrm{},+\mathrm{})`$ for which the curve $`c(t)`$ defined by

$$c(T(s))=\gamma (s)s(s_0,s_1)$$
(4.2)
satisfies $`\dot{c}(t)=𝐗_H(c(t))`$ for any $`t𝐑`$. Denoting the derivative with respect to $`s`$ by prime, it follows that $`s(s_0,s_1)`$ we have

$$\chi (s)\gamma ^{}(s)=𝐗_H(\gamma (s))\text{with}\chi (s)=\frac{1}{T^{}(s)}.$$
(4.3)
Because of the smoothness of the right hand side as a function of $`s𝐑`$, we observe that a unique extension of $`\chi `$ to $`[s_0,s_1]`$ must exist. This extended function must clearly satisfy the conditions

$$\chi (s_0)=\chi (s_1)=0,\chi ^{}(s_0+0)0,\chi ^{}(s_10)0.$$
(4.4)
By using that (4.3) holds on $`[s_0,s_1]`$ and taking the appropriate derivatives of this equality at the endpoints, one arrives at the following statement.

Proposition 3. If $`\gamma (s)`$ in (4.1) yields a heteroclinic orbit in the above-described sense, then the vectors $`ad_Y^{}x_iT_{x_i}𝒪`$ are eigenvectors of the linearisation of $`𝐗_H`$ at $`x_i`$, for $`i=0,1`$, with the respective eigenvalues being $`\chi ^{}(s_0+0)`$ and $`\chi ^{}(s_10)`$.

In particular, notice from the proposition that the existence of a real eigenvalue of the linearisation of $`𝐗_H`$ at $`x_0`$ is a necessary condition for the construction of a heteroclinic orbit through $`x_0`$ by means of a 1-parameter subgroup of $`G`$. For the rigid body systems described in section 2, this is in fact also a sufficient condition. For $`n=3`$ this is a well known result. We verify it below in the $`n=4`$ case by using the explicit analysis of the preceding section.

As before we may assume that the equilibrium point of interest is $`x_0=(ae_3,e_3)`$ in (3.9). Let $`(a\xi e_3,\eta e_3)=[Y,x_0]`$ be an eigenvector of the linearised flow at $`x_0`$ with real eigenvalue $`z>0`$. Note that $`z=0`$ is excluded by (3.11) and that we have $`Y=(\xi ,\eta )`$ by using the identification of the Lie bracket of $`so(3)su(2)`$ with the vector-product. Then we can check that

$$\xi _1^2+\xi _2{}_{}{}^{2}=\eta _1{}_{}{}^{2}+\eta _2{}_{}{}^{2}.$$
(4.5)
This follows from the eigenvector equation $`Lv=zv`$ with $`L`$ in (3.10) and $`v=(\xi _1,\eta _1,\xi _2,\eta _2)^T`$. We set $`\mathrm{\Delta }:=\sqrt{\xi _1^2+\xi _2^2}`$ and consider the curve

$$\gamma (s)=e^{sY}x_0e^{sY}=\mathrm{cos}(s\mathrm{\Delta })(ae_3,e_3)+\mathrm{\Delta }^1\mathrm{sin}(s\mathrm{\Delta })(a\xi e_3,\eta e_3).$$
(4.6)
We can verify that this curve yields a heteroclinic orbit that connects $`x_0`$ with $`x_1:=(ae_3,e_3)`$ for $`s_0=0`$ and $`s_1=\frac{\pi }{\mathrm{\Delta }}`$. Indeed, the functions $`\chi `$ and $`T`$ introduced in (4.3) are found as

$$\chi (s)=\frac{z}{\mathrm{\Delta }}\mathrm{sin}(s\mathrm{\Delta }),$$
(4.7)
$$T(s)=\frac{1}{z}\mathrm{log}\mathrm{tan}\frac{s\mathrm{\Delta }}{2}\text{for}s_0<s<s_1.$$
(4.8)
Note that the adjoint and coadjoint actions are the same for any compact Lie group and any orbit $`𝒪_x=G/G_x`$ carries a canonical $`G`$-invariant Riemannian metric induced by the Cartan-Killing form on $`𝐠`$. It is well known that the geodesics of this metric coincide with the orbits of the 1-parameter subgroups of $`G`$. Thus the heteroclinic orbits considered above are proper generalizations of the heteroclinic orbits of the standard rigid body that are great circles on $`S^2=SO(3)/SO(2)`$.

## 5 On the stability analysis for $`n>4`$

We are able to repeat a large part of the stability analysis performed in the 4-dimensional case. Specifically: we can find the equilibrium points (Proposition 1); we can find the eigenvectors and corresponding eigenvalues of the linearised system at each equilibrium point; we can prove the converse of Proposition 3. However the problem of proving stability (or not) for the equilibrium points having all eigenvalues pure imaginary is more complicated. We present here only an outline of the stability analysis for general $`n`$.

To find the eigenvalues and eigenvectors of the linearised system (2.20) at $`x`$ it is useful to decompose $`so(n)`$ as the vector space direct sum $`so(n)=`$ Ker$`(ad_x)+`$ Im$`(ad_x)`$, whereby we can uniquely parametrize $`vT_x𝒪_x`$ as $`v=[Y,x]`$ with $`Y`$ Im$`(ad_x)`$. The linearised system (2.20) then reads as

$$[\dot{Y},x]=[J,[Y,x^2]]$$
(5.1)
and an eigenvector $`[Y,x]T_x𝒪_x`$ with eigenvalue $`z`$ satisfies

$$[J,[Y,x^2]]=z[Y,x].$$
(5.2)
Let us take $`x`$ to be of the form (2.3) and choose coordinates on $`so(n)`$ according to the natural decomposition into blocks. That is, for $`n=2m`$ write $`YIm(ad_x)`$ as $`Y=AA^T`$ with $`A=_{i<j}e_{ij}\xi _{ij}`$ and $`\xi _{ij}`$ a real $`2\times 2`$ real matrix. If $`n=2m+1`$, then $`Im(ad_x)Y=\left(\begin{array}{cc}AA^T& v\\ v^T& 0\end{array}\right)`$ with $`A`$ as before and $`v^T=(v_1^T,v_2^T,\mathrm{},v_m^T)`$ with $`v_i`$ a real $`2\times 1`$ matrix. Writing (5.2) in these coordinates, we see directly that there are several copies of the eigenvector equation for $`so(4)`$ \- each of which has 4 solutions - and in the odd $`n`$ case also several copies of the eigenvector equation for $`so(3)`$ \- each of which has 2 solutions. In fact we obtain exactly the right number of such decoupled equations to generate all eigenvectors and their eigenvalues. If any eigenvalue is real and nonzero, then we can use either the result described for $`so(4)`$ or a similar one - which has not been explicitly described here, but which is straightforward - for $`so(3)`$, to construct heteroclinic orbits by suitable curves of the form in (4.1). This leads to the following converse of Proposition 3.

Proposition 4. Suppose that $`z`$ is a nonzero, real eigenvalue of the linear system (5.2) at $`x`$. Then there exists a corresponding eigenvector $`[Y,x]`$ for which the curve $`\gamma (s)=e^{sY}xe^{sY}`$ yields a heteroclinic orbit of the rigid body system (2.7).

We now sketch the proof of this proposition in the $`n=2m`$ case. In this case we can write $`J=_{i=1}^me_{ii}D_i`$, where the $`D_i`$ are 2 by 2 diagonal matrices. By putting

$$Y:=\underset{1i<jm}{}Y_{ij}\text{with}Y_{ij}:=e_{ij}\xi _{ij}e_{ji}\xi _{ij}^T,$$
(5.3)
the eigenvector equation (5.2) decouples into separate equations for each pair of indices $`i<j`$,

$$(x_i^2x_j^2)(D_i\xi _{ij}\xi _{ij}D_j)=z(x_j\xi _{ij}Sx_iS\xi _{ij}),S:=\mathrm{i}\sigma _2.$$
(5.4)
For any $`1i<jm`$, consider the $`so(4)`$ subalgebra of $`so(2m)`$ given by

$$so(4)_{ij}:=\mathrm{span}\{e_{ii}S,e_{jj}S,(e_{ij}Qe_{ji}\times Q^T)|Qgl(2,𝐑)\}.$$
(5.5)
The point to notice is that (5.4) coincides with the eigenvalue equation for a rigid body system defined on $`so(4)_{ij}`$ at the corresponding equilibrium point $`x_{ij}:=x_ie_{ii}S+x_je_{jj}S`$. This implies by the $`so(4)`$ result established in section 4 that if $`Y_{ij}`$ is a solution of (5.4) with some real $`z0`$, then the curve

$$\gamma _{ij}(s):=e^{sY_{ij}}x_{ij}e^{sY_{ij}}$$
(5.6)
yields a heteroclinic orbit connecting the unstable equilibria $`\pm x_{ij}`$ of the induced rigid body system on $`so(4)_{ij}`$. Decomposing $`x`$ as $`x=x_{ij}+x_{ij}^{}`$, we can check the relations

$$\gamma (s):=e^{sY_{ij}}xe^{sY_{ij}}=\gamma _{ij}(s)+x_{ij}^{}$$
(5.7)
and

$$[J,\gamma ^2(s)]=[J_{ij},\gamma _{ij}^2(s)],J_{ij}:=e_{ii}D_i+e_{jj}D_j.$$
(5.8)
Eq. (5.8) relates the Hamiltonian vector fields for the rigid body systems on $`so(2m)`$ and on $`so(4)_{ij}`$ along the respective curves $`\gamma (s)`$ and $`\gamma _{ij}(s)`$. By collecting the above remarks, we conclude that $`\gamma (s)`$ in (5.7) yields a heteroclinic orbit that connects the unstable equilibria $`x_{ij}^{}\pm x_{ij}`$.

To illustrate what happens for odd $`n`$, let us look at $`n=5`$. Let us assume that the equilibrium point $`x`$ of interest has the form

$$x=\left(\begin{array}{ccccc}0& a1& 0& 0& 0\\ a+1& 0& 0& 0& 0\\ 0& 0& 0& a1& 0\\ 0& 0& a+1& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),$$
(5.9)
and parametrize $`Y`$ Im$`(ad_x)`$ according to

$$Y=\frac{1}{2}\left(\begin{array}{ccccc}0& 0& \xi _2\eta _2& \xi _1\eta _1& v_1\\ 0& 0& \xi _1\eta _1& \xi _2+\eta _2& v_2\\ \eta _2\xi _2& \eta _1+\xi _1& 0& 0& w_1\\ \eta _1\xi _1& \xi _2\eta _2& 0& 0& w_2\\ v_1& v_2& w_1& w_2& 0\end{array}\right).$$
(5.10)
Then consider the eigenvector equation (5.2) with $`J=\mathrm{diag}(J_1,J_2,J_3,J_4,J_5)`$. By setting $`v_i=0=w_i`$ we reduce to the eigenvector condition for $`so(4)`$; by setting $`\xi _i=\eta _i=w_i=0`$ we reduce to the eigenvector condition for $`so(3)`$ and by setting $`\xi _i=\eta _i=v_i=0`$ we reduce to the eigenvector condition for $`so(3)`$ too. In fact the coordinates have been chosen here so as to agree exactly with those used for the so(4) analysis in section 3. In this way we find all 8 eigenvalues. The problems of checking if the eigenvalues are real, complex or imaginary reduce to those of the $`so(3)`$ and $`so(4)`$ cases. Similarly the construction of heteroclinic orbits as orbits of 1-parameter subgroups reduces to the $`so(3)`$ and $`so(4)`$ cases. To check if all eigenvalues of the linearisation being imaginary is sufficient for stability we could try to prove the convexity at $`x`$ of a function of the form

$$f=\alpha H+\beta H_1+\gamma H_2+\delta H_4$$
(5.11)
with $`H`$, $`H_1`$, $`H_2`$, $`H_3`$ the Hamiltonian together with 3 independent commuting integrals, which can be generated using the Lax matrix of Manakov , where $`\alpha `$, $`\beta `$, $`\gamma `$, $`\delta `$ are expected to depend on the equilibrium point in question. Of course, while even this can be done in principle, there is no strategy telling us how to proceed for general $`n`$.

## 6 Conclusion

In this paper we studied the equilibrium points for the integrable Euler equations in (2.7). In particular, we described the equilibrium points (Proposition 1) and associated heteroclinic orbits with any nonzero, real eigenvalue of the linearised system for any $`n`$ (Proposition 4). We also found a complete characterization of the stability of the equilibrium points for $`n=4`$ (Proposition 2), but our stability analysis is incomplete for $`n>4`$. In this case an open question is to find a criterion for the stability of those equilibrium points for which all eigenvalues of the linearised system are imaginary.

As a final remark, we wish to mention the work of Mishchenko and Fomenko  (for a review, see ) that contains generalizations of the systems of Manakov  to other Lie algebras. Various elements of our results have a general Lie-algebraic nature and thus may be applicable to the systems of . In this respect, it is natural to ask if Proposition 4 is valid only for the special cases (2.7) that we considered here, or can be extended to other systems among those in , too. It would also be interesting to find a general criterion of stability that could be used effectively to analyse these systems.

## Acknowledgments

We wish to thank T.S. Ratiu for posing the problem studied here and for discussions. L.F. was supported in part by the Hungarian Scientific Research Fund (OTKA) under T034170, T030099, T029802 and M036804.


format: text

text:
# ”

## Abstract

Temperature dependence of an electron-nuclear plasma equilibrium density is considered basing on known approaches, which are given in -. It is shown that at a very high temperature, which is characteristic for a star interior, the equilibrium plasma density is almost constant and equals approximately to $`10^{25}`$ particles per $`cm^3`$. At a relatively low temperature, which is characteristic for star surface, the equilibrium plasma density is in several orders lower and depends on temperature as $`T^{3/2}`$.

The temperature dependence of equilibrium plasma density

B.V.Vasiliev

Institute in Physical-Technical Problems, 141980, Dubna, Russia

vasiliev@dubna.ru

PACS: 64.30.+i; 95.30.-k

The electron-nuclear plasma has a high density at high temperature inside stars. It may be supposed that a star exists in a steady state, and the plasma is in a thermodynamic equilibrium. And so it may be expected that an equilibrium plasma should have a density which is determined by its energy minimum.

It will be shown below that an energy minimum of equilibrium plasma exists really. It enables to determinate an equilibrium density of plasma as a function of its temperature. It is important that the density of inner regions of stars is equal to $`a_0^310^{25}`$ particles in $`cm^3`$ approximately ($`a_0`$ is Bohr radius) and it is in several orders lower for surface regions, where it depends on temperature as $`T^{3/2}`$.

Let us consider an ensemble consisting of a permanent number of particles N at constant temperature T. The equilibrium state at a minimum of free energy F will be settled by a ensemble volume V change, i.e. by a change of a particle density n:

$$\left(\frac{F}{n}\right)_{N,T}=0.$$
(1)
The chemical potential describes an assemble equilibrium. The direct interaction between nuclei can be neglected in a dense electrically neutral equilibrium plasma. The plasma energy is determined by electron-electron interaction and electron-nucleus interaction.

The chemical potential of an electron gas at high temperature (in Boltzman approximation) is known :

$$\mu _B=kTln\xi _B$$
(2)
where

$$\xi _=\frac{1}{2}\left[\frac{2\pi h^2}{mkT}\right]^{3/2}n_e.$$
(3)
and the electron gas density

$$n_e=\frac{2^{1/2}m^{3/2}}{\pi ^2\mathrm{}^3}_0^{\mathrm{}}\frac{\epsilon ^{1/2}d\epsilon }{e^{(\epsilon \mu )/kT}+1}.$$
(4)
We will search for the chemical potential of plasma at arbitrary temperature as

$$\mu =kTln\xi $$
(5)
Using the substitution

$$\zeta =\frac{\epsilon }{kT},$$
(6)
from (4) we obtain

$$\xi =\xi _B\frac{\pi ^{1/2}}{2I_1}$$
(7)
where

$$I_1=_0^{\mathrm{}}\frac{\zeta ^{1/2}d\zeta }{e^\zeta +\xi }$$
(8)
Thus

$$\frac{d\xi _B}{d\xi }=\frac{2}{\pi ^{1/2}}I_1\left(1\xi \frac{I_3}{I_1}\right).$$
(9)
where

$$I_3=_0^{\mathrm{}}\frac{\zeta ^{1/2}d\zeta }{(e^\zeta +\xi )^2}.$$
(10)
The free energy of an electron gas is

$$F_e=\frac{N(2m)^{3/2}(kT)^{5/2}\xi }{3\pi ^2\mathrm{}^3n_e}_0^{\mathrm{}}\frac{\zeta ^{3/2}d\zeta }{e^\zeta +\xi }=\frac{2I_2}{3I_1}NkT$$
(11)
where

$$I_2=_0^{\mathrm{}}\frac{\zeta ^{3/2}d\zeta }{e^\zeta +\xi }.$$
(12)
At high temperature (in Boltzman approximation) when $`\xi =\xi _B1`$ and chemical potential $`\mu =\mu _B<0`$ and $`|\mu _B|1`$, the free energy of an electron gas may be expanded in series. If we conserve the two first terms of the series, we obtain

$$F_e=F_{ideal}+N\frac{\pi ^{3/2}a_o^{3/2}e^3}{4(kT)^{1/2}}n_e$$
(13)
Here the first term is the free energy of an ideal gas and the second term is the correction for the identity of electrons. This correction is positive as it takes into account that electrons can not take places which are already occupied by other electrons and it causes an increasing incompressibility of electron gas.

From Eq.(13) we obtain

$`\left({\displaystyle \frac{dF_e}{dn_e}}\right)_{N,T}={\displaystyle \frac{2^{1/2}\pi ^2a_o^{3/2}e^3}{3(kT)^{1/2}}}{\displaystyle \frac{(I_4I_1I_3I_2)}{I_1^3(1\xi I_3/I_1)}}`$ (14)
where $`a_o`$ is the Bohr radius and

$$I_4=_0^{\mathrm{}}\frac{\zeta ^{3/2}d\zeta }{(e^\zeta +\xi )^2}.$$
(15)
In order to take into account a role of nuclei in the plasma energy formation, we must calculate the correlation correction :

$`\delta F_{corr}=N{\displaystyle \frac{2\pi ^{1/2}e^3kT}{3}}\left[Z^2{\displaystyle \frac{n_i}{kT}}+\left({\displaystyle \frac{n_e}{\mu }}\right)_{N,T}\right]^{3/2}`$ (16)
where Z is the nucleus charge, $`n_i=n_e/Z`$ is the nucleus density.

At high temperature (in Boltzman approximation) this correction is

$$\delta F_{corr}=N\frac{2\pi ^{1/2}e^3(Z+1)^{3/2}}{3(kT)^{1/2}}n_e^{1/2}$$
(17)
This correction takes into account that nuclei ”condense” the electron gas in their closeness. On the contrary, the electron gas density decreases away from the nuclei and its compressibility increases. This explains why this correction is negative.

Thus the free energy of plasma at high temperature with account for both corrections

$`F=F_{ideal}+N{\displaystyle \frac{\pi ^{3/2}a_o^{3/2}e^3}{4(kT)^{1/2}}}n_eN{\displaystyle \frac{2\pi ^{1/2}e^3(Z+1)^{3/2}}{3(kT)^{1/2}}}n_e^{1/2}`$ (18)
i.e. the free energy of plasma has two correction added to ideal value. They have different signs and differently depend on particle density, but they equally depend on temperature. This enables to calculate the equilibrium density of a dense plasma at high temperature from balance condition Eq.(1)

$$n_0=\frac{16(Z+1)^3}{9\pi ^2a_o^3}$$
(19)
In this case the Fermi energy of electron gas $`\epsilon _F\frac{Z^2e^2}{a_0}`$ and the Boltzman approximation is applicable when $`kT\epsilon _F`$. In this approximation the equilibrium density does not depend on temperature.

For arbitrary temperature starting from Eq.(16), we obtain

$`\left({\displaystyle \frac{dF_{corr}}{dn_e}}\right)_{N,T}={\displaystyle \frac{\pi ^{1/2}e^3(Z+1)^{3/2}}{3(kTn_e)^{1/2}}}[1{\displaystyle \frac{\xi I_3}{(Z+1)I_1}}]^{3/2}`$
$`\left[1{\displaystyle \frac{3\xi }{(Z+1)I_1}}{\displaystyle \frac{(I_32\xi I_5+\xi I_3^2/I_1)}{(1\xi I_3/I_1)[1\frac{\xi I_3}{(Z+1)I_1}]}}\right]`$ (20)
The equilibrium density of plasma is

$`{\displaystyle \frac{n_e}{n_0}}=\{{\displaystyle \frac{3}{2^{5/2}\pi ^{1/2}}}{\displaystyle \frac{I_1^3(1\xi I_3/I_1)}{(I_4I_1I_2I_3)}}[1{\displaystyle \frac{\xi I_3}{(Z+1)I_1}}]^{3/2}`$
$`[1{\displaystyle \frac{3\xi }{(Z+1)I_1}}{\displaystyle \frac{(I_32\xi I_5+\xi I_3^2/I_1)}{(1\xi I_3/I_1)[1\frac{\xi I_3/I_1}{(Z+1)}]}}]\}^2`$ (21)
where

$$I_5=_0^{\mathrm{}}\frac{\zeta ^{1/2}d\zeta }{(e^\zeta +\xi )^3}.$$
(22)
At the same time

$$T=\frac{2\pi a_0e^2}{k}\left(\frac{n_e}{\pi ^{1/2}\xi I_1}\right)^{2/3}.$$
(23)
The numerically calculated dependence of the equilibrium plasma density on temperature at Z=1 obtained from these equations is shown in Fig.1.

Thus, at a high temperature $`TE_F/k`$ the equilibrium density of plasma approaches to $`n_0`$, and at a relatively low temperature ($`TE_F/k`$)

$$n_eT^{3/2}.$$
(24)
It is a consequence of that the full energy of electron gas in a general case is

$$E_e=N\frac{2^{1/2}m^{3/2}(kT)^{5/2}\xi }{\pi ^2\mathrm{}^3n_e}_0^{\mathrm{}}\frac{\zeta ^{3/2}d\zeta }{e^\zeta +\xi }$$
(25)
or with account for above formulas

$$E_e=NkT\frac{I_2}{I_1}.$$
(26)
For a degenerate electron gas

$$E_e=\frac{3}{10}NE_F.$$
(27)
As in this case the ratio $`I_2/I_11`$, from Eq.(26) we obtain for the equilibrium degenerate plasma

$$E_FkT.$$
(28)
and with account for the definition of $`E_F`$, we obtain Eq.(24).

The obtained equilibrium values of a plasma density enables to explain a series of astrophysical effects: the proportionality of star magnetic moments to their angular momenta (the Blackett effect), the star mass spectrum and others .


text:
# Ground state and elementary excitations of single and binary Bose-Einstein condensates of trapped dipolar gases

## I Introduction

The nature and stability of a Bose-Einstein condensate (BEC)  are strongly influenced by the interparticle interactions, which are described by the $`s`$-wave scattering length $`a`$. If $`a>0`$ the interactions are repulsive, and condensates with an arbitrary large number of particles are stable. On the contrary, spatially homogeneous condensates with $`a<0`$ are absolutely unstable with regard to local collapses . The presence of a trapping potential changes the situation drastically, as revealed in successful experiments with magnetically trapped atomic 7Li ($`a=14`$ Å) . As found in theoretical studies , there will be a metastable BEC if the number of condensed particles is sufficiently small, such that the spacing between the trap levels exceeds the mean field interparticle interaction $`n_0|g|`$ (where $`n_0`$ is the condensate density, $`g=4\pi \mathrm{}^2a/m`$, and $`m`$ is the atom mass). In other words, the BEC is stabilized if the negative pressure caused by the interparticle attraction is compensated by the quantum pressure imposed by the trapping potential.

The effects of the interparticle interactions on the condensate properties have been mainly discussed for the case of van der Waals (short-range) interactions. However, the BEC in the presence of dipole-dipole interactions has recently raised a considerable interest . Novel physics is expected for dipolar BEC, since the dipole-dipole interactions are long-range, anisotropic and partially attractive. The non-trivial task of achieving and controlling dipolar BECs is thus particularly challenging.

The interest on dipolar gases has been partially motivated by the recent success in creating ultra-cold molecular clouds . This success opens fascinating prospects to achieve quantum degeneracy in trapped gases of heteronuclear molecules, which could interact via electric dipole-dipole forces after being oriented in a sufficiently high electric field. On the other hand, the ultra-cold gases of atoms with large magnetic dipole moments, as chromium  or europium , have also been subject of growing interest. In this case, the dipole-dipole interactions are not expected to be dominant, although for a relatively small $`a`$ the BEC may reflect the interplay between short-range and dipole-dipole interactions . Interestingly, these effects can be amplified by reducing $`a`$ via Feshbach resonances . Similar effects have been discussed in Refs.  for ground-state atoms with electric dipole moments induced by a high dc electric fields (of the order of 106 V/cm). It has been also suggested that laser-induced dipole-dipole interactions could be achieved by exciting atoms to Rydberg states . In this case applications to quantum information processing have been discussed (see e.g. ).

The stability of the condensate is significantly modified by the presence of dipole-dipole interactions . In particular, a BEC of particles dominantly interacting via dipole forces is, similarly to condensates with $`a<0`$, unstable in a spatially homogeneous case and can be stabilized by confinement in a trap. It has been shown  that the sign and the value of the dipole-dipole interaction energy is strongly influenced by the trapping geometry and, hence, the stability diagram depends crucially on the trap anisotropy. This opens new interesting possibilities for controlling and engineering macroscopic quantum states. In particular, for dipoles oriented along the axis of a cylindrical trap there exists a critical value $`l_{}=0.4`$ for the square root of the ratio of the radial to axial frequency $`l=(\omega _\rho /\omega _z)^{1/2}`$. Pancake traps with $`l<l_{}`$ provide mostly a repulsive mean field of the dipole-dipole interaction, and thus the dipolar condensate in such traps will be stable at any number of particles $`N`$. For $`l>l_{}`$ the stability requires $`N<N_c`$, where the critical value $`N_c`$ at which the collapse occurs is determined by the condition that (on average) the mean field interaction energy is attractive and close in absolute value to $`\mathrm{}\omega _\rho `$.

The study of the condensate properties would be incomplete without the analysis of the excitation spectrum, which determines the dynamical behavior of the system in the regime of weak perturbations. Zero-temperature excitation frequencies have been extensively studied in the case of condensates in dilute alkali gases, both experimentally  and theoretically . In this paper we discuss the low-lying collective excitation frequencies of trapped dipolar BEC, complementing the analysis of Refs. . We first consider the case of dominant dipole-dipole interactions, and later on we discuss the situation where the short-range interaction is also relevant. In particular, we discuss in detail the nature of instability and demonstrate that one of the excitations frequencies tends to zero at the criticality as a power of $`(N_cN)^\beta `$. We discuss the qualitative character of the low-lying modes, and show that the exponent $`\beta `$ undergoes a crossover from $`1/4`$ for $`ll_{}`$ to $`2`$ at $`ll_{}`$. To illustrate the case of mixed short-range and dipole-dipole interactions, we present predictions for the excitation frequencies for the particular case of a chromium gas.

The second part of the paper is devoted to the analysis of binary dipolar BECs. In recent years the development of trapping techniques has allowed for creation of multicomponent condensates, formed by atoms in different internal (electronic) states . The physics of multicomponent BEC, far from being a trivial extension of the single-component one, presents novel and fundamentally different scenarios for its ground-state wave function  and excitations . In particular, it has been experimentally observed that a BEC can reach an equilibrium state characterized by phase separation of the species in different domains . The analysis of multicomponent BEC has been so far mostly limited to the case of short-range interparticle interactions (a model long-range interaction has been considered in Ref. ). One of the main goals of this paper is to analyze the properties of multicomponent dipolar BEC. Such a mixture can be achieved in different physical contexts. In particular, it would be in general the case for experiments in ultracold polar molecules , and in chromium  in which different magnetic-moment species are simultaneously trapped in an optical dipole trap. It would also be the case of atomic electric dipoles created by laser-induced pumping to two different Rydberg states. Finally, the same situation would appear in condensates of heteronuclear Hund A diatomic molecules, for which the direction of the magnetic moment is correlated (parallel or antiparallel) with the direction of the molecular axis. Thus, if the magnetic moments are oriented in a magnetic field, the electric moments can acquire two possible directions. We show below that a binary dipolar BEC of two antiparallel dipole components differs qualitatively from the case of a short-range interacting binary BEC.

Our paper is organized as follows. In Sec. II we briefly review the ground-state properties of a single component BEC of trapped dipolar gases . Sec. III is devoted to the analysis of the excitation spectrum of single component trapped dipolar condensed gases. Sec. IV briefly discusses the ballistic expansion of a dipolar BEC. In Sec. V the ground state of a BEC of two different dipolar species is considered. The excitation spectrum for this case is discussed in Sec. VI. We conclude in Sec. VII.

## II Ground-state properties of a single-component dipolar BEC

### A Description of the system

In this section we briefly review the results of Ref. . We consider a condensate of dipolar particles in a cylindrical harmonic trap. In the following we consider the case of electric dipoles, although the same physics is expected for magnetic ones . All dipoles are assumed to be oriented along the trap axis by a sufficiently strong external field. Accordingly, the dipole-dipole interaction potential between two dipoles is given by $`V_d(𝐑)=(d^2/R^3)(13\mathrm{cos}^2\theta )`$, where $`d`$ characterizes the dipole moment, $`𝐑`$ the vector between the dipoles ($`R=|𝐑|`$ being its length), and $`\theta `$ the angle between $`𝐑`$ and the dipole orientation. Similarly as in Refs. , we describe the dynamics of the condensate wave function $`\psi (𝐫,t)`$ by using the time-dependent Gross-Pitaevskii equation (GPE):

$`i\mathrm{}{\displaystyle \frac{}{t}}\psi (𝐫,t)=\{{\displaystyle \frac{\mathrm{}^2}{2m}}^2+{\displaystyle \frac{m}{2}}(\omega _\rho ^2\rho ^2+\omega _z^2z^2)+\mathrm{\Gamma }0`$ (1)
$`\mathrm{\Gamma }0+g|\psi (𝐫,t)|^2+d^2{\displaystyle }d𝐫^{}{\displaystyle \frac{13\mathrm{cos}^2\theta }{|𝐫𝐫^{}|^3}}|\psi (𝐫^{},t)|^2\}\psi (𝐫,t).`$ (2)
where $`\psi (𝐫,t)`$ is normalized to the total number of condensate particles $`N`$. The third term in the rhs of Eq. (2) corresponds to the mean field of short-range (van der Waals) forces and the last term to the mean field of the dipole-dipole interaction. In this section, we omit the term $`g|\psi (𝐫,t)|^2\psi (𝐫,t)`$, assuming that the interparticle interaction is dominated by dipole-dipole forces ($`d^2|g|=4\pi \mathrm{}^2|a|/m`$) and the system is away from shape resonances of $`V_d(𝐑)`$. The effects of the short-range interactions on the excitation spectrum are discussed in detail in Sec. III.

The wave function of the relative motion of a pair of dipoles is influenced by the dipole-dipole interaction at interparticle distances $`|𝐫𝐫^{}|r_{}=2md^2/\mathrm{}^2`$. This influence is ignored in the dipole-dipole term of Eq.(2), as the main contribution to the integral comes from distances $`|𝐫𝐫^{}|`$ of order the spatial size of the condensate, which we assume to be much larger than $`r_{}`$.

Away from the shape resonances the dipolar condensate is unstable in the spatially homogeneous case. For all dipoles parallel to each other, by using the Bogoliubov method, one finds an anisotropic dispersion law for elementary excitations: $`\epsilon (𝐤)=[E_k^2+8\pi E_kn_0d^2(\mathrm{cos}^2\theta _k1/3)]^{1/2}`$, where $`E_k=\mathrm{}^2k^2/2m`$, $`n_0`$ is the condensate density, and $`\theta _k`$ is the angle between the excitation momentum $`𝐤`$ and the direction of the dipoles. The instability is clearly seen from the fact that at small $`k`$ and $`\mathrm{cos}^2\theta _k<1/3`$ imaginary excitation energies $`\epsilon `$ emerge.

### B Numerical results

Equation (2), contrary to the usually employed GPE with short-range interactions, is an integro-differential equation. The evaluation of the integral term deserves special attention, since the integrand diverges at relative interparticle distances tending to zero. Fortunately, the calculation of the integral term can be simplified by means of the convolution theorem:

$$d^2𝑑𝐫^{}\frac{13\mathrm{cos}^2\theta }{|𝐫𝐫^{}|^3}|\psi (𝐫^{})|^2=^1\left\{[V](𝐪)[|\psi |^2](𝐪)\right\},$$
(3)
where $``$ and $`^1`$ indicate the Fourier transform and the inverse Fourier transform, respectively. The Fourier transform of the dipole-dipole potential reads:

$$[V](𝐪)=4\pi d^2(13\mathrm{cos}^2\alpha )[\frac{\mathrm{cos}(bq)}{(bq)^2}\frac{\mathrm{sin}(bq)}{(bq)^3}],$$
(4)
where $`\alpha `$ is the angle between the momentum $`𝐪`$ and the dipole direction, and $`b`$ is a cutoff distance corresponding to the atomic radius (a few Bohr radii). Since $`b`$ is much smaller than any significant length scale of the system, one can safely perform the limit

$$\underset{b0}{lim}(V(𝐫))=\frac{4\pi }{3}d^2(3\mathrm{cos}^2\alpha 1).$$
(5)
In order to evaluate the Fourier transform of Eq. (3) $`(|\psi |^2)`$ is numerically evaluated by means of a standard Fast Fourier Transform (FFT) algorithm and multiplied by $`(V(𝐫))`$.

The ground-state of the system is obtained by employing a standard split-operator technique in imaginary time. The split-operator technique is also based in an FFT algorithm and, consequently, for each time step four FFT’s are needed: two for the calculation of the integral term and two for the evolution. We would like to stress that this procedure constitutes a non trivial computational task. Additionally, the FFT algorithm must be evaluated in Cartesian coordinates and, as a consequence, computationally demanding fully three-dimensional calculations are required.

To understand the influence of the trapping potential on the dipolar condensate, we have simulated Eq.(2) for various values of the number of particles $`N`$, dipole moment $`d`$, and the trap aspect ratio $`l`$. We have found the conditions under which the condensate is stabilized by the trapping field and investigated static properties of this Bose-condensed state.

For a stationary condensate the wave function $`\psi (𝐫,t)=\psi _0(𝐫)\mathrm{exp}(i\mu t/\mathrm{})`$, where $`\mu `$ is the chemical potential, and the lhs of Eq. (2) becomes $`\mu \psi _0(𝐫)`$. The important energy scales of the problem are the trap frequencies $`\omega _z`$, $`\omega _\rho `$ and the dipole-dipole interaction energy per particle defined as $`V=(1/N)V_d(𝐫𝐫^{})\psi _0^2(𝐫)\psi _0^2(𝐫^{})𝑑𝐫𝑑𝐫^{}`$. Accordingly, the trap frequencies, and the (renormalized) number of particles $`\sigma =Nr_{}/a_{\mathrm{max}}`$ (with $`a_{\mathrm{max}}=(\mathrm{}/2m\omega _{\mathrm{min}})^{1/2}`$ being the maximal oscillator length of the trap) form the necessary set of parameters allowing us to determine the chemical potential and give a full description of the ground state of a trapped dipolar condensate.

We have found that the dipolar condensate is stable either at $`V>0`$ or at $`V<0`$ with $`|V|<\mathrm{}\omega _\rho `$. This requires $`N<N_c`$, where the critical number $`N_c`$ depends on the trap aspect ratio $`l=(\omega _\rho /\omega _z)^{1/2}`$. The calculated dependence $`N_c(l)`$ clearly indicates the presence of a critical point $`l_{}0.43`$ . In pancake traps with $`l<l_{}`$ the condensate is stable at any $`N`$, because $`V`$ always remains positive. For small $`N`$ the shape of the cloud is Gaussian in all directions. With increasing $`N`$, the quantity $`V`$ increases and the cloud first becomes Thomas-Fermi in the radial direction and then, for a very large $`N`$, also axially. The ratio of the axial to radial size of the cloud, $`L=L_z/L_\rho `$, continuously decreases with increasing number of particles and reaches a limiting value at $`N\mathrm{}`$ (see Fig. 3 of Ref.). In this respect, for a very large $`N`$ we have a pancake Thomas-Fermi condensate.

For $`l1`$ the mean field dipole-dipole interaction is always attractive. The quantity $`|V|`$ increases with $`N`$ and the shape of the cloud changes. In spherical traps the cloud becomes more elongated in the axial direction and near $`N=N_c`$ the shape of the cloud is close to Gaussian with the aspect ratio $`L=2.1`$. In cigar-shaped traps ($`l1`$) especially interesting is the regime where $`\mathrm{}\omega _z|V|\mathrm{}\omega _\rho `$. In this case the radial shape of the cloud remains the same Gaussian as in a non-interacting gas, but the axial behavior of the condensate will be governed by the dipole-dipole interaction which acquires a quasi 1-dimensional (1D) character. Thus, one has a (quasi) 1D gas with attractive interparticle interactions and is dealing with a stable (bright) soliton-like condensate where attractive forces are compensated by the kinetic energy . With increasing $`N`$, $`L_z`$ decreases. Near $`N=N_c`$, where $`|V|`$ is close to $`\mathrm{}\omega _\rho `$, the axial shape of the cloud also becomes Gaussian and the aspect ratio takes the value $`L3.0`$. For $`l_{}ł<1`$ the dipole-dipole interaction energy is positive for a small number of particles and increases with $`N`$. The quantity $`V`$ reaches its maximum and further increase in $`N`$ reduces $`V`$ and makes the cloud less pancake-shaped. At the critical point $`N=N_c`$ the shape of the cloud is close to Gaussian and the aspect ratio $`L<3.0`$.

In the previous analysis, the case of dominant dipole-dipole interactions was considered. However , in the general case the effects of the short-range interactions are comparable or even larger than those related to the dipole-dipole interactions. In such situations the short-range term must be maintained in Eq. (2). Provided that $`a`$ is sufficiently small and positive, the system can become unstable and undergo a collapse in a similar way to what was observed in experiments with 7Li  and 85Rb . For the case of negative $`a`$ the dipolar gas is expected to be highly unstable, but the dipole-dipole interaction could be employed to stabilize the gas by the trap geometry in a way analogous to the one presented in Ref. .

## III Excitations of a single-component dipolar BEC

In this section we analyze the collective excitations of a dipolar BEC. Since this is a potentially unstable system, there is a fundamental question about the qualitative and quantitative nature of this instability. Another important question concerns quantitative aspects, in particular how relevant are the effects of the dipole-dipole interaction in the excitation spectrum and to what extent they can observed in the experiments.

In the subsection III A, we briefly summarize the variational approach introduced by Pérez-García et al.  and later used by Yi and You  to describe the low-lying excitation spectrum of a dipolar condensate (for a different variational method using the self-similarity assumption see ). We employ the notation introduced in the Refs. . Subsection III B is devoted to the analysis of the behavior of the excitations close to the instability. In subsection III C, we calculate numerically the response of the system to small perturbations and compare the results with those yielded by the method of subsection III B. Finally, in subsection III C, we also discuss how the effects of the dipole-dipole interactions can manifest themselves in the excitation spectrum for various kinds of dipolar gases and trapping geometries.

### A Variational approach

The problem of solving the time dependent GPE Eq. (2) can be restated as a variational problem corresponding to the stationary point of the action $`S=𝑑t𝑑𝐫`$ related to the Lagrangian density :

$``$ $`=`$ $`{\displaystyle \frac{i}{2}}\mathrm{}\left[\psi (𝐫){\displaystyle \frac{\psi ^{}(𝐫)}{t}}\psi ^{}(𝐫){\displaystyle \frac{\psi (𝐫)}{t}}\right]`$ (6)
$`+`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}|\psi (𝐫)|^2+V_\mathrm{t}(𝐫)|\psi (𝐫)|^2`$ (7)
$`+`$ $`{\displaystyle \frac{g}{2}}|\psi (𝐫)|^4+{\displaystyle \frac{u_2}{2}}|\psi (𝐫)|^2{\displaystyle 𝑑𝐫^{}\frac{Y_{20}(\theta )}{|𝐫𝐫^{}|^3}|\psi (𝐫^{})|^2},`$ (8)
where $`u_2=4\sqrt{\pi /5}d^2`$ and $`V_\mathrm{t}(𝐫)`$ describes the trapping potential. We consider the following Gaussian ansatz for the condensate wave function:

$`\psi (x,y,z,t)=A(t){\displaystyle \underset{\eta =x,y,z}{}}e^{\eta ^2/2w_\eta ^2}e^{i\eta ^2\beta _\eta (t)},`$ (9)
where in contrast to Ref.  we omit the possibility of the sloshing motion of the condensate. By inserting the ansatz (9) into Eq. (8) and integrating over the spatial coordinates, one obtains the effective Lagrangian. From this effective Lagrangian one finds the equations of motion of the variational parameters, i.e. the corresponding Euler-Lagrange equations. In particular, if $`w_\eta =v_\eta \sqrt{\mathrm{}/m\omega }`$, $`\omega =\omega _{x,y}`$, $`\omega _z=\lambda \omega _x`$, $`v_{x,y}=v`$, $`P=\sqrt{2/\pi }Na/\sqrt{\mathrm{}/m\omega }`$, and $`\tau =\omega t`$, the Euler-Lagrange equations take the following form:

$`{\displaystyle \frac{d^2}{d\tau ^2}}v_\eta +\lambda _\eta ^2v_\eta ={\displaystyle \frac{1}{v_\eta ^3}}`$ (10)
$`P{\displaystyle \frac{}{v_\eta }}\left[{\displaystyle \frac{1}{v_xv_yv_z}}\left(1+{\displaystyle \frac{u_2}{g}}{\displaystyle 𝑑𝐫\mathrm{exp}\left(\underset{\eta }{}\frac{\eta ^2}{2v_\eta ^2}\right)\frac{Y_{20}(\theta )}{r^3}}\right)\right],`$ (11)
where $`\lambda _{x,y}=1`$, $`\lambda _z=\lambda `$. This equation describes the motion of a particle with coordinates $`(v_x,v_y,v_z)`$ in an effective potential

$`U(v_x,v_y,v_z)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\lambda _x^2v_x^2+\lambda _y^2v_y^2+\lambda _z^2v_z^2\right)`$ (12)
$`+`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{v_x^2}}+{\displaystyle \frac{1}{v_y^2}}+{\displaystyle \frac{1}{v_z^2}}\right)+{\displaystyle \frac{P}{v_xv_yv_z}}`$ (13)
$`\times `$ $`\left[1+{\displaystyle \frac{u_2}{g}}{\displaystyle 𝑑𝐫\mathrm{exp}\left(\underset{\eta }{}\frac{\eta ^2}{2v_\eta ^2}\right)\frac{Y_{20}(\theta )}{r^3}}\right].`$ (14)
Small amplitude oscillations around the stationary solution can be found by studying second derivatives of $`U`$. This procedure provides the excitation frequency for the first three compressional excitation modes of the system . The geometry of these modes is depicted in Fig. 1. For cylindrical traps with the axis along the dipole direction, the projection of the angular momentum, $`m`$, on the $`z`$ axis is a good quantum number. For modes $`2`$ and $`3`$ we have $`m=0`$, whereas $`m=1`$ for mode $`1`$. In the following we call mode $`2`$ ($`3`$) as the breathing (quadrupole) mode.

### B The nature of the instability

To study the nature of the instability one needs to determine which mode becomes unstable when the number of atoms reaches the critical value and to describe the behavior of the frequency of the corresponding mode close to the criticality. These two issues have been first discussed for the case of short-range interacting Bose gases with $`a<0`$. Bergeman  observed numerically that, as the ratio $`\gamma `$ of the nonlinear interaction energy to the trap frequency approaches a given critical value $`\gamma _c`$, the frequency of the breathing mode tends to zero and merges with the frequency of the Goldstone mode corresponding to the overall phase of the condensate. Above the criticality, the breathing mode becomes unstable and attains complex frequency. Singh and Rokhsar  analyzed this instability using self-similar solutions describing the modes (equivalently one may employ the variational approach of the previous subsection). They have shown that close to the criticality the frequency of the breathing mode $`2`$ vanishes as $`|\gamma \gamma _c|^{1/4}`$.

In the case of a dipolar gas with dominant dipole interactions the situation is completely different and much more complex. Only for aspect ratios far above the criticality, $`ll^{}`$ ($`l>1.29`$) the situation resembles that of a gas with $`a<0`$. The mode corresponding to the lowest frequency is the breathing mode. This mode becomes unstable when the parameter $`\sigma \sigma _c`$. The scaling behavior of the frequency of this mode can be analyzed employing the variational approach of the previous section . We find that $`\omega _2`$ goes to zero as $`(\sigma _c\sigma )^\beta `$, with $`\beta 1/4`$. For intermediate values of $`l>l`$ ($`0.75<l<1.29`$) the exponent $`\beta `$ is still $`1/4`$, but the geometry of the mode which drives the instability depends on $`\sigma `$. For $`\sigma `$ far below $`\sigma _c`$ the mode corresponding to the lowest frequency has a breathing symmetry, whereas as one approaches the critical value of $`\sigma `$ the lowest mode becomes quadrupole-like. For $`l`$ close to $`l^{}`$ ($`l<0.75`$) the situation changes and the mode corresponding to the lowest frequency is quadrupolar-like. For $`l`$ not too close to $`l^{}`$ the exponent $`\beta `$ is still close to $`1/4`$. However, as $`l`$ approaches $`l^{}`$, the exponent $`\beta `$ departs from $`1/4`$ towards a greater value, which becomes $`2`$ at $`l=l_{}`$. The latter prediction must be understood as a qualitative argument, since for $`l`$ very close to $`l_{}`$, the system enters a Thomas-Fermi regime where the Gaussian ansatz is no more valid.

### C Numerical results

In this subsection we study the low-lying collective excitations of a trapped dipolar BEC by analyzing numerically the response of a BEC after applying an external perturbation. We compare our numerical results with those obtained from the variational approach described above. In the following, we employ the method introduced in Ref.  for the case of a short-range interacting BEC. After generating the BEC ground-state wave function using imaginary time evolution of the GPE, we first slightly perturb the trapping potential in a periodic way:

$$V(𝐫,t)=\frac{1}{2}m\underset{\eta =x,y,z}{}[1+A_\eta \mathrm{sin}(\omega _{mod}t+\alpha _\eta )]^2\omega _\eta ^2\eta ^2,$$
(15)
where $`\omega _{mod}`$ is the modulation frequency, $`A_\eta `$ are the amplitudes, and $`\alpha _\eta `$ are the initial phases. A perturbation of a chosen symmetry can be accomplished by properly selecting these parameters. In particular, sufficiently small amplitudes are necessary near the instabilities, since the system can be easily driven into collapse. On the other hand, by setting large amplitudes one may probe various nonlinear effects . The response of the system is enhanced by selecting $`\omega _{mod}`$ in the vicinity of expected excitation frequencies.

In a second stage, the condensate evolves in an unperturbed trap (i.e. $`A_\eta =0`$). The condensate widths are monitored and subsequently Fourier-transformed to reveal the excitation spectrum. Our aim is to determine the (potentially small) deviations of the excitation frequencies with respect to those expected for purely short-range interacting gases, as precisely as they can be experimentally measured (typically $`\mathrm{\Delta }\omega /\omega _{trap}0.01`$ ). Long integration times are necessary to accomplish the desired spectral resolution. This introduces serious technical difficulties, since a large number of integration steps is needed to guarantee the energy conservation during the whole evolution. Note that additional complications arise from the evaluation of the interaction integral in Eq. (2) at each time step.

#### 1 Dominant dipole-dipole interactions

For the case of dominant dipole-dipole interactions , the short-range part of Eq. (2) can be safely omitted. We consider two different trap geometries, that were employed at JILA in Ref. , with a trap aspect ratio $`l=8^{1/4}`$ (referred to as pancake JILA trap), and a trap with $`l=8^{1/4}`$ (referred to as cigar JILA trap). The later trap has been recently employed for experiments in 85Rb .

In the following, we study the excitation frequencies as a function of the dipolar parameter

$$\zeta =N\frac{m}{\mathrm{}^2}\sqrt{\frac{m\omega _\rho }{\mathrm{}}}d^2.$$
(16)
Note the relation $`\sigma =2\sqrt{2\omega _\rho /\omega _{min}}\zeta `$. Since in the pancake JILA trap the breathing mode has a relatively high frequency and it is not excited by the perturbation of the GPE, we concentrate in this case on the modes $`1`$ and $`3`$ of Fig. 1. The results for the pancake JILA trap are presented in Fig. 2. For $`\zeta =0`$ the ideal-gas result $`\omega /\omega _\rho =2`$ is retrieved for both modes $`1`$ and $`3`$. As $`\zeta `$ increases mode $`1`$ is shifted upwards, whereas the frequency of mode $`3`$ (quadrupole) goes down, vanishing for a critical $`\zeta _{cr}`$ at which the system becomes unstable. Probing the system in this regime is very difficult, since even a slight disturbance drives the collapse. For $`\zeta >\zeta _{cr}3.68`$ it is not possible to obtain stable ground state solutions of the time-independent GPE . The variational analysis reproduces the numerical results for relatively weak dipolar interactions. However, it does not describe well the numerical results close to the instability.

Fig. 3 depicts a similar dependence for the cigar JILA trap. In this case the breathing mode is the lowest one. The frequency of mode $`1`$ displays an upward shift whereas that of the quadrupole mode stays essentially untouched. Again, the instability threshold is not correctly predicted by the variational analysis.

#### 2 Interplay between short-range and dipole-dipole interactions

Let us now consider the case in which the short-range interactions and the dipole-dipole ones have a comparable strength . In particular, we have performed our simulations for the particular case of 52Cr, which has drawn some experimental interest , since it possesses a large magnetic dipole moment of $`6\mu _B`$ (Bohr magnetons). The value of $`a`$ for this element is at the moment unknown, and consequently we have explored different values of $`a`$ in the calculations below. In the following, the trap parameters correspond to those of an ongoing experiment in Stuttgart , namely $`\omega _z=2\pi \times 40`$ Hz and $`\omega _\rho =2\pi \times 485`$ Hz.

Fig. 4 shows the three lowest modes for different values of $`a`$, for the case of $`N=10000`$ atoms. When $`a`$ is smaller than a critical $`a_{crit}`$ the system becomes unstable due to an unbalanced attractive component of dipole-dipole interactions . When approaching $`a^{crit}/a_{\mathrm{Na}}0.23`$ from above, the frequency of the lowest mode decreases to $`0`$ and the system collapses. For $`0.23<a/a_{Na}<0.26a_{Na}`$ the lowest mode (triangles) has a breathing geometry, whereas the highest one (circles) is quadrupole-like. The opposite is true for $`a/a_{Na}>0.26a_{Na}`$. The highest mode shows virtually no dependence on the scattering length whereas the frequency of mode $`1`$ is shifted upwards. Again, close to the instability, the results obtained from the variational calculation differ from those obtained numerically. Moreover, the behavior of the lowest mode is not reproduced correctly even qualitatively. Fig. 4 also presents the excitation frequencies for the same system in the absence of dipole-dipole interactions. Apart from the obvious presence of the instability in the dipolar case (and the corresponding behavior of the lowest mode) we observe large frequency differences for mode $`1`$ between the cases with and without dipole-dipole interactions.

The above result suggests that the effects of dipole-dipole interactions may be detectable. Having discussed the role of $`a`$, it is thus important to discuss possible values of $`N`$ and $`\omega _\rho /\omega _z`$ which could maximize the dipole-induced frequency shifts.

Fig. 5 shows the deviation between the case with and without dipole-dipole interactions, for the mode $`1`$, $`\omega _z/\omega _\rho =40/485`$, and $`a=a_{Na}`$. For this calculation we have employed the variational method as it is very exact away from the instability. We observe that for number of atoms $`N`$ greater than $`5000`$ the frequency shift is observable (i.e. greater than $`0.01`$), and reaches its maximum ($`3.5\%`$) for $`N10000`$ which should be the number available in ongoing experiments . In this case, shifts for modes $`2`$ and $`3`$ are substantially lower ($`0.001`$).

Finally, we have analyzed (using the variational method) the dependence of the spectrum on the trap aspect ratio $`l`$ for $`a=a_{Na}`$ and $`N=10000`$. The results are presented in Fig. 6. We observe that in pancake traps with $`l<0.7`$, the quadrupole mode experiences a shift of up to $`10\%`$ with respect to the nondipolar case. We have confirmed this conclusion by performing an exact numerical simulation of the perturbed GPE for $`l=0.58`$, $`a=a_{Na}`$, and $`N=10000`$. For these parameters the excitation frequency of mode $`3`$ is $`1.76\omega _\rho `$ ($`1.86\omega _\rho `$) with (without) dipole-dipole interactions. The remaining two modes also display large shifts over a wide range of aspect ratios.

## IV Ballistic expansion of a dipolar BEC

In typical BEC experiments, the measurements are performed after removing the trapping potential and allowing for a ballistic expansion of the condensate. It is the aim of this section to briefly discuss such an expansion, by means of the previously introduced variational approach. After evaluating the minimum of the potential provided by Eq. (14), i.e. the ground state of the trapped gas, we set $`\lambda _{x,y,z}=0`$ in Eqs. (11) and evolve these equations using the ground state as an initial condition. In this section we restrict ourselves to the case of dominant dipole-dipole interactions.

The ballistic expansion of a short-range interacting BEC is characterized by an inversion of the aspect ratio of the condensate cloud, i.e. cigar-shaped condensates become pancake-like after expansion, and vice versa. This is not necessarily the case in dipolar condensates. Fig. 7 shows the condensate aspect ratio for a spherical trap, before releasing the trap (dashed), and once the aspect ratio reaches a stationary value after the expansion (solid). The aspect ratio of the cloud decreases during the expansion, but never becomes smaller than $`1`$, i.e. the condensate keeps its original cigar-shaped character during the expansion. A similar behavior is observed for $`1<l<1.12`$, where there exits a range of $`\zeta `$ values for which the BEC remains cigar-like. For $`l_{}<l<1`$, for $`\zeta `$ close to $`\zeta _{cr}`$ the BEC also remains cigar-shaped during the time of flight. On the contrary, for $`l>1.12`$ the cloud becomes pancake-like after expanding for all $`\zeta `$, and for $`l<l_{}`$ the initial pancake cloud always becomes cigar-shaped.

## V Ground state of a binary dipolar BEC

In the previous sections we have shown that various novel phenomena are expected in dipolar BEC. One can therefore expect that even a richer behavior can be displayed by a multicomponent dipolar BEC. In the case of two-component dipolar gases, the system properties depend on a large number of control parameters, including the number of particles in each component, the strengths and orientations of the dipoles, and the trap geometries. In the following, for simplicity, we focus on the case of a BEC of two dipolar identical components polarized in opposite directions. As discussed in Sec. I this should be the case for Hund A molecules in a magnetic field, where the magnetic moment is oriented by the field, but the (possibly large) electric moment can be parallel or antiparallel to the direction of the applied magnetic field. To simplify the analysis even further, we consider only a spherical trap. Thus, the system is determined by $`\zeta `$ and $`\eta =N_1/N`$, where as before $`N`$ is the total number of particles and $`N_1`$ is the number of particles in the component $`1`$. The system can be described by a system of two coupled GPEs:

$`i\dot{\psi }_1=`$ (18)
$`\{{\displaystyle \frac{^2}{2}}+{\displaystyle \frac{r^2}{2}}++\zeta (\eta V_1(𝐫)(1\eta )V_2(𝐫))\}\psi _1,`$ (19)
$`i\dot{\psi }_2=`$ (20)
$`\{{\displaystyle \frac{^2}{2}}+{\displaystyle \frac{r^2}{2}}++\zeta (\eta V_1(𝐫)+(1\eta )V_2(𝐫))\}\psi _2,`$ (21)
where $`V_i(𝐫)=𝑑𝐫^{}(13\mathrm{cos}^2\theta )|\psi _i(𝐫^{},t)|^2/|𝐫𝐫^{}|^3`$. In the above expression we have employed harmonic oscillator units of length $`d=\sqrt{\mathrm{}/m\omega }`$. For equal population of the components ($`\eta =0.5`$) both ground-state wave functions are the same. As a consequence, the nonlocal interaction terms vanish and the system behaves like an ideal gas, the ground-state wave functions being that of the corresponding harmonic oscillator. However, as the interaction strength $`\zeta `$ increases, the system becomes unstable. The reason of this instability is rooted in the excitation spectrum and is discussed below.

For each value of $`\eta `$ it is possible to compute a critical $`\zeta _{cr}`$ above which no stable solution of Eq. (V) exists. The corresponding stability diagram, resulting from a numerical solution of Eqs. (V), is presented in Fig. 8. We observe that the more symmetric the mixture is the more stable it is.

For $`\eta 0.5`$ the ground state has a more complicated structure. In this respect it resembles solutions for short-range interacting binary condensates . In particular, the component whose self interaction is stronger partially wraps around the other. However, contrary to short-range interacting binary BEC, due to the anisotropy of the dipole-dipole interactions, the density dip in one component appears only along the radial direction, whereas along $`z`$ both components possess Gaussian profiles. An example of such a structure is shown in Fig. 9.

In the analysis presented in this section we have neglected the influence of short-range interactions. Since the two components are identical except for the opposite dipole moments, the short-range interactions (both intercomponent and the intracomponent ones) are the same. Therefore, in the absence of dipole-dipole interactions a miscible mixture is expected, which differs from the ideal-gas Gaussian case considered above (eventually acquiring a Thomas-Fermi profile). However, in the same way as above, the dipole-dipole mean fields will be exactly canceled for $`\eta =0.5`$. Therefore, the cloud is expected to remain in an unperturbed Thomas-Fermi profile until the $`\zeta `$ coefficient reaches a critical value $`\zeta _{cr}`$ (only quantitatively different from that obtained in absence of short-range interactions). For $`\zeta >\zeta _{cr}`$ the system becomes unstable.

## VI Excitations of a binary dipolar BEC

In this section we analyze the excitation spectrum of the previously discussed gas of two oppositely oriented dipolar components. We first discuss the homogeneous case, and then study the case of harmonically trapped binary dipolar BEC.

### A Homogeneous space

Let us consider two species of dipoles $`𝐝_1=d_1𝐮_z`$ and $`𝐝_2=d_2𝐮_z`$ with respective densities $`n_1`$ and $`n_2`$ in a homogeneous box of volume $`V`$. The Hamiltonian in second quantization is:

$`\widehat{H}={\displaystyle \underset{𝐤\alpha =1,2}{}}{\displaystyle \frac{\mathrm{}^2𝐤^2}{2m}}a_{𝐤\alpha }^{}a_{𝐤\alpha }`$ (22)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1,2}{}}{\displaystyle \underset{𝐤_1,𝐤_2,𝐤_1^{},𝐤_2^{}}{}}g(𝐤_1^{},𝐤_2^{};𝐤_1,𝐤_2)a_{𝐤_1^{}\alpha }^{}a_{𝐤_2^{}\alpha }^{}a_{𝐤_2\alpha }a_{𝐤_1\alpha }`$ (23)
$`{\displaystyle \underset{𝐤_1,𝐤_2,𝐤_1^{},𝐤_2^{}}{}}g(𝐤_1^{},𝐤_2^{};𝐤_1,𝐤_2)a_{𝐤_1^{}1}^{}a_{𝐤_2^{}2}^{}a_{𝐤_22}a_{𝐤_11},`$ (24)
with

$`g(𝐤_1^{},𝐤_2^{};𝐤_1,𝐤_2)=`$ (25)
$`{\displaystyle 𝑑𝐫𝑑𝐫^{}\psi _{𝐤_1^{}}^{}(𝐫)\psi _{𝐤_2^{}}^{}(𝐫^{})V(𝐫𝐫^{})\psi _{𝐤_2}(𝐫^{})\psi _{𝐤_1}(𝐫)},`$ (26)
and $`a_{𝐤i}`$ the annihilation operator of a particle of the $`i`$-th component with momentum $`𝐤`$. Performing the Bogoliubov approximation ($`a_{01}^{}=a_{01}=\sqrt{N_1}`$ and $`a_{02}^{}=a_{02}=\sqrt{N_2}`$), one obtains

$`\widehat{H}={\displaystyle \underset{𝐤0}{}}\{[ϵ_𝐤+2B(\widehat{𝐤})]a_{𝐤1}^{}a_{𝐤1}+[ϵ_𝐤+2\lambda ^2B(\widehat{𝐤})]a_{𝐤2}^{}a_{𝐤2}`$ (27)
$`+B(\widehat{𝐤})[(a_{𝐤1}^{}a_{𝐤1}^{}+a_{𝐤1}a_{𝐤1})+\lambda ^2(a_{𝐤2}^{}a_{𝐤2}^{}+a_{𝐤2}a_{𝐤2})`$ (28)
$`2\lambda (a_{𝐤1}a_{𝐤2}+a_{𝐤1}^{}a_{𝐤2}+a_{𝐤1}^{}a_{𝐤2}+a_{𝐤2}^{}a_{𝐤1})]\},`$ (29)
with $`ϵ_𝐤=\mathrm{}^2𝐤^2/2m`$, $`\sqrt{n_2d_2^2}=\lambda \sqrt{n_1d_1^2}`$ and $`B(\widehat{𝐤})=\frac{8\pi }{3}\sqrt{\frac{\pi }{5}}(n_1d_1^2)Y_{20}(\widehat{𝐤})`$. The diagonalization of the Hamiltonian provides two different branches of quasiparticle excitations:

$`E_1=ϵ_𝐤,`$ (31)
$`E_2=\sqrt{ϵ_𝐤^2+{\displaystyle \frac{32\pi }{3}}\sqrt{{\displaystyle \frac{\pi }{5}}}(1+\lambda ^2)(n_1d_1^2)Y_{20}(\widehat{𝐤})ϵ_𝐤}.`$ (32)
For the case of $`\lambda 0`$ the spectrum of a homogeneous single-component dipolar BEC is discussed in Sec. II. As observed, the spectra (31) is that of a soft mode, whereas the second branch (32) is even more unstable than that of a single component BEC.

### B Trapped case

For simplicity we limit ourselves to the case of a spherical trap with frequency $`\omega `$, but similar arguments can be employed for a cylindrical trap. We also consider $`|𝐝_1|=|𝐝_2|`$ and $`\eta =0.5`$ (i.e. $`N_1=N_2=N/2`$). As we have already noted, in this case the ground state of the system is the ground state of a harmonic oscillator. It is therefore more convenient to rewrite the Hamiltonian in terms of the creation and annihilation operators of the harmonic oscillator modes for the dipole components $`1`$ and $`2`$:

$`\widehat{H}={\displaystyle \underset{𝐧}{}}\mathrm{}\omega _𝐧(a_{𝐧1}^{}a_{𝐧1}+a_{𝐧2}^{}a_{𝐧2})+{\displaystyle \frac{N_1d^2}{2}}{\displaystyle \underset{𝐧_1,𝐧_20}{}}`$ (33)
$`\{g(𝐧_1,𝐧_2;0,0)a_{𝐧_11}^{}a_{𝐧_21}^{}+g(0,0;𝐧_1,𝐧_2)a_{𝐧_11}a_{𝐧_21}`$ (34)
$`+g(𝐧_1,𝐧_2;0,0)a_{𝐧_12}^{}a_{𝐧_22}^{}+g(0,0;𝐧_1,𝐧_2)a_{𝐧_12}a_{𝐧_22}`$ (35)
$`+2g(𝐧_1,0;𝐧_2,0)a_{𝐧_11}^{}a_{𝐧_21}+2g(𝐧_1,0;𝐧_2,0)a_{𝐧_12}^{}a_{𝐧_22}`$ (36)
$`2g(0,0;𝐧_1,𝐧_2)a_{𝐧_11}a_{𝐧_22}2g(𝐧_1,𝐧_2;0,0)a_{𝐧_11}^{}a_{𝐧_22}^{}`$ (37)
$`2g(𝐧_1,0;0,𝐧_2)a_{𝐧_11}^{}a_{𝐧_22}2g(𝐧_1,0;0,𝐧_2)a_{𝐧_12}^{}a_{𝐧_21}\},`$ (38)
where $`𝐧(n,l,m)`$ denotes the set of spherical quantum numbers for the corresponding eigenstate of the harmonic oscillator, $`\mathrm{}\omega _𝐧`$ is the energy of the eigenstate $`𝐧`$, and

$`g(𝐧_1,𝐧_2;𝐧_3,𝐧_4)=`$ (39)
$`{\displaystyle 𝑑r𝑑r^{}\psi _{𝐧_1}^{}(r)\psi _{𝐧_2}^{}(r^{})V(rr^{})\psi _{𝐧_3}(r^{})\psi _{𝐧_4}(r)}.`$ (40)
Let us perform the Bogoliubov transformation:

$$b_\nu =\underset{𝐧0}{}(u_{\nu 𝐧}a_{𝐧1}+v_{𝐧\nu }a_{𝐧1}^{}+\stackrel{~}{u}_{\nu 𝐧}a_{𝐧2}+\stackrel{~}{v}_{𝐧\nu }a_{𝐧2}^{}),$$
(41)
which satisfies $`ϵ_\nu b_\nu =[b_\nu ,\widehat{H}]`$, leading to the corresponding Bogoliubov-de Gennes (BdG) equations. Such equations can be simplified by taking the more appropriate variables $`A_{\nu 𝐧}=u_{\nu 𝐧}+\stackrel{~}{u}_{\nu 𝐧}`$, $`B_{\nu 𝐧}=v_{\nu 𝐧}+\stackrel{~}{v}_{\nu 𝐧}`$, $`C_{\nu 𝐧}=u_{\nu 𝐧}\stackrel{~}{u}_{\nu 𝐧}`$, $`D_{\nu 𝐧}=v_{\nu 𝐧}\stackrel{~}{v}_{\nu 𝐧}`$. Then, the BdG equations take the form:

$`ϵ_\nu A_{\nu 𝐧}=\mathrm{}\omega _𝐧A_{\nu 𝐧},`$ (43)
$`ϵ_\nu B_{\nu 𝐧}=\mathrm{}\omega _𝐧B_{\nu 𝐧},`$ (44)
$`ϵ_\nu C_{\nu 𝐧}=\mathrm{}\omega _𝐧C_{\nu 𝐧}`$ (45)
$`+2N_1d^2{\displaystyle \underset{𝐧^{}0}{}}\{g(𝐧^{},0;𝐧,0)C_{\nu 𝐧^{}}g(0,0;𝐧,𝐧^{})D_{\nu 𝐧^{}}\},`$ (46)
$`ϵ_\nu D_{\nu 𝐧}=\mathrm{}\omega _𝐧D_{\nu 𝐧}`$ (47)
$`+2N_1d^2{\displaystyle \underset{𝐧^{}0}{}}\{g(𝐧,𝐧^{};0,0)C_{\nu 𝐧}g(𝐧,0;𝐧^{},0)D_{\nu 𝐧^{}}\}.`$ (48)
In the following, we omit for simplicity $`\nu `$ in the notation for the $`A`$, $`B`$, $`C`$ and $`D`$ coefficients. We observe that, like in the homogeneous case, we have two different branches of quasi-particles, one of them being a soft mode. The last two equations are the BdG equations for the non-soft mode. From the properties of spherical harmonics and the $`3j`$-symbols, it is possible to obtain that $`g(0,0;nlm,n^{}l^{}m^{})=g(nlm,n^{}l^{}m^{};0,0)=g(0,0;nlm,n^{}l^{}m)\delta _{m^{},m}`$ $`g(nlm,0;n^{}l^{}m^{},0)=g(n^{}l^{}m^{},0;nlm,0)=g(0,0;nlm,n^{}l^{}m)\delta _{m^{},m}`$. Also, $`l^{}=l,l\pm 2`$ must be fulfilled. Employing these properties, and changing to the variables $`F_{nlm}=C_{nlm}(1)^mD_{nlm}`$, $`G_{nlm}=C_{nlm}+(1)^mD_{nlm}`$, one obtains that $`G_{nlm}=\frac{ϵ}{\mathrm{}\omega _n}F_{nlm}`$, and

$`ϵ^2F_{nlm}=(\mathrm{}\omega _n)^2F_{nlm}`$ (49)
$`+4N_1d^2\mathrm{}\omega _n{\displaystyle \underset{n0}{}}g(0,0;nlm,n^{}l^{}m)F_{n^{}l^{}m}.`$ (50)
One can write:

$`g(0,0;nlm,n^{}l^{}m^{})=`$ (51)
$`(1)^{(l^{}l)/2}{\displaystyle \frac{1}{2\sqrt{2}a_0^3}}h(n,l;n^{},l^{})l^{}m|l,m|20,`$ (52)
where $`a_0=\sqrt{\mathrm{}/2m\omega }`$. In Eq. (52) the angular contribution is given by

$`l^{},m|lm|20=`$ (53)
$`(1)^m\left({\displaystyle \genfrac{}{}{0pt}{}{l^{}l\mathrm{\hspace{0.17em}2}}{mm0}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{l^{}l\mathrm{\hspace{0.17em}2}}{000}}\right)\left[{\displaystyle \frac{5(2l+1)(2l^{}+1)}{4\pi }}\right]^{1/2},`$ (54)
whereas the coefficients $`h(n,l;n^{},l^{})`$ constitute the radial contribution and are of the form:

$`h(n,l;n^{},l^{})={\displaystyle \frac{8}{3\sqrt{5}}}{\displaystyle \frac{1}{2^{(n+n^{})/2}}}`$ (55)
$`\times {\displaystyle \frac{\mathrm{\Gamma }\left(\frac{n+n^{}+3}{2}\right)}{\sqrt{\mathrm{\Gamma }\left(\frac{n+l+3}{2}\right)\mathrm{\Gamma }\left(\frac{nl+2}{2}\right)\mathrm{\Gamma }\left(\frac{n^{}+l^{}+3}{2}\right)\mathrm{\Gamma }\left(\frac{n^{}l^{}+2}{2}\right)}}}.`$ (56)
Then the problem reduces to finding the eigenvalues of the following equation:

$`ϵ^2F_{nlm}=n^2F_{nlm}+\sqrt{2}\zeta n`$ (57)
$`{\displaystyle \underset{l^{}=l,l\pm 2}{}}{\displaystyle \underset{n^{}l^{}}{}}(1)^{(ll^{})/2}h(nl;n^{}l^{})l^{}m|lm|20F_{n^{}l^{}m},`$ (58)
where $`\zeta `$ is defined in Eq. (16), and $`ϵ`$ is in units of $`(\mathrm{}\omega )^2`$ . From the analysis of the spectrum one can observe that for $`\zeta >\zeta _{cr}=6.7`$ the energy of the lowest state of the subspace $`n=2`$, $`m=0`$ becomes imaginary. Therefore, although the ground-state of the mixture of antiparallel dipolar components is that of an ideal gas, the system will eventually become unstable for $`\eta >\eta _{cr}`$ due to the appearance of imaginary excitation frequencies. The value of $`\zeta _{cr}`$ compares very well with the numerically found value of $`\zeta _{cr}7.3`$ for the onset of the instability.

Finally, it is interesting to analyze the situation in which $`\zeta 1`$, since in this situation the states with different $`n`$ do not mix, and one can analytically diagonalize in each $`n`$ subspace. The first excited states are (in units of $`(\mathrm{}\omega )^2`$):

$`ϵ^2(n=1,l=1,m=0)=1+4\sqrt{2/\pi }\zeta /15,`$ (59)
$`ϵ^2(n=1,l=1,m=\pm 1)=12\sqrt{2/\pi }\zeta /15,`$ (60)
$`ϵ^2(n=2,l=\{2,0\},m=0)=4+(2\pm \sqrt{102})\sqrt{2/\pi }\zeta /21,`$ (61)
$`ϵ^2(n=2,l=2,m=\pm 1)=4+2\sqrt{2/\pi }\zeta /21.`$ (62)
From these expressions it becomes clear that, for example, in the subspace of states with $`n=2`$, $`m=0`$ (which is an admixture of the noninteracting quadrupole and monopole modes) even for $`\zeta 0.1`$ deviations of more than $`1\%`$ are expected. One can also observe that vortices with vorticity in the dipole direction are less energetic than those with vorticity in the plane perpendicular to the dipole direction.

## VII conclusions

In this paper we have analyzed the ground state and the excitation spectrum of single- and two-component dipolar condensates. For the case of single-component BECs we have shown that their stability properties are determined by the trapping geometry. In particular, for sufficiently pancake-like traps the condensate is always stable independent of the number of particles. We have then analyzed the excitation spectrum of the single-component condensate by analyzing the response of the condensate to small perturbations and comparing the results with analytical calculations based on a variational approach. We have discussed in detail the nature of the instability and associated it with vanishing of frequency of one of the excitation modes. The scaling behavior of this frequency was also analyzed. We have discussed different possible scenarios in which the dipole-dipole effects could have observable effects even in the presence of a dominant short-range interaction. This analysis could be of special interest for ongoing experiments on atoms with large magnetic moments, such as chromium . We have provided guidelines for experimental parameters corresponding to largest discrepancies between the cases with and without dipole-dipole interactions. In the second part of the paper we have analyzed the properties of a two-component BEC of dipolar particles. In particular, we have studied the stability of the ground state as a function of the relative density of both species as well as the appearance of phase separation in the binary condensate. Finally, we have obtained the excitation spectrum for this particular physical system and discussed the nature of its instability.

We acknowledge support from the Alexander Von Humboldt Foundation, the ZIP Program of the German Government, the Deutscher Akademischer Austauschdienst (DAAD), the Deutsche Forschungsgemeinschaft, the RTN Cold Quantum Gases, ESF Program BEC2000+, and the subsidy of the Foundation for Polish Science. K.G. acknowledges support by the Polish KBN grant 5 P03B 102 20. We thank M. Gajda, M. Lewenstein, G. Meijer, T. Pfau, K. Rza̧żewski, G. V. Shlyapnikov, E. Tiemann and L. You for fruitful discussions. Part of the results was obtained using computers at the Interdisciplinary Center for Mathematical and Computational Modeling at Warsaw University. K. G. and L. S. are grateful for the kind hospitality extended to them in Hannover and Warsaw, respectively.