\section{Introduction \label{sec:introduction}}

Since Bose-Einstein condensates (BECs) can be routinely prepared in laboratories,
ultracold gases became a very important testbed for many predictions of condensed
matter physics \cite{Bloch2008}. The experimental examination of binary condensates
started nearly the same time as for single condensates by using two different
quantum states of the same species, such as $^{87}$Rb \cite{Myatt1997} or $^{23}$Na
\cite{Stamper-Kurn1998}. With the development of sympathetic cooling ultracold
mixtures have been assembled from two different alkalies, $^{41}$K--$^{87}$Rb
\cite{Modugno2001, Thalhammer2008}, $^{7}$Li--$^{133}$Cs \cite{Mudrich2002} and
$^{87}$Rb--$^{133}$Cs \cite{Pilch2009} or for different isotopes of the same alkali
atom $^{85}$Rb--$^{87}$Rb \cite{Papp2008}. The tunability of the inter- and intra-species
scattering lengths via driving the mixture through a Feshbach-resonance
has also been experimentally demonstrated \cite{Thalhammer2008, Pilch2009}.

The ability to create Bose-Einstein condensate(s), a highly coherent form of matter,
also facilitated the convergence of two fields of physics: condensed matter physics
and quantum optics, and therefore BECs became favourable candidates for examining
the effects of nonlinearity in matter waves, where this nonlinearity originates from
the mean-field representation of the interatomic interaction. The similarity between
electromagnetic waves in nonlinear medium and coherent matter waves is also expressed
in the equations of motion which is the nonlinear Schr{\"o}dinger equation (NLS) for
the former and the Gross-Pitaevskii (GP) equation for the latter. Although the physical
interpretation of these equations is different, their structures are the same apart
from the external potential term. Furthermore, in some cases this extra term can even
be removed~\cite{Schumayer2002} and the GP equation is transformed into a form
coinciding with the NLS equation. Consequently all results for the NLS equation
known in nonlinear optics can be readily adapted to Bose-Einstein condensates.

One of the surprising phenomena of nonlinear optics is the existence of particle-like
wave-forms, the so-called solitons~\cite{Hasegawa2002}. Such excitations have already
been experimentally observed in single- or two-component Bose-Einstein condensates:
dark solitons \cite{Burger1999, Denschlag2000}, bright solitons \cite{Khaykovich2002,
Strecker2002}, their two-component coupled analogues, the dark-dark \cite{Ohberg2001},
bright-bright \cite{Liu2002} or even dark-bright \cite{Busch2001, Becker2008}
multi-component solitary waves \cite{Ginsberg2005, Becker2008}.

However, the question of existence of solitons needs more attention than simply
recognising the similarity between the two governing equations. The existence of
solitons is strongly related to the integrability of the given physical model.
An usual test to determine whether or not an equation is integrable is the
Painlev{\'e}-test (P-test) \cite{Steeb1984}. It was shown that without any potential
term or inhomogeneity the one-dimensional NLS equation, $i u_{t} + u_{xx} \pm 2
\abs{u}^{2} u = 0$, is completely integrable both in its one-component
\cite{Zakharov1972, Steeb1984} or multi-component \cite{Zakharov1974, Zakharov1985}
form. However, the inclusion of a potential term, $v(x,t) u$, in the one-component
NLS equation or different coupling strengths in the multi-component NLS/GP equations
fundamentally changes their integrability \cite{Sahadevan1986, Clarkson1988}. It
was shown that integrability is preserved provided the external potential, $v(x,t)$,
has a specific form \cite{Clarkson1988}. The authors examined \cite{Schumayer2001}
the integrability of the two-component coupled Gross-Pitaevskii (CGP) equations and
lead to similar conclusion: the scattering lengths and the external potentials cannot
be arbitrary if the integrability of CGP is to be preserved. The inter- and intraspecies
scattering lengths must satisfy the following equation
\begin{equation} \label{eq:ConditionForIntegrability}
   \frac{%
         2 \xi_{1} \xi_{2}
         - \kappa_{1} \xi_{1}
         - \kappa_{2} \xi_{2}
        }%
        {
         \xi_{1} \xi_{2}
         - \kappa_{1} \kappa_{2}
        }
   =
   \!\frac{(2n+1)^{2} + 7}{16},
\end{equation}
where $\xi_{1}=a_{11}/a_{21}$, $\xi_{2}=a_{22}/a_{12}$, $\kappa_{1}=\mu_{11}/ \mu_{21}$,
$\kappa_{2}=\mu_{22}/ \mu_{12}$ and $\mu_{ij}$ denotes the reduced mass of pair of
particles composed by an atom from the $i$th and $j$th species. On the right hand side
of Eq.~(\ref{eq:ConditionForIntegrability}) $n$ is a non-negative integer number. One
may call $n$ a {\emph{classification number}}, because it determines the form of the
external potentials for which CGP equations remain integrable. For example, for $n=2$,
the external potential, apart from the quadratic trapping potential, may even contain
an imaginary time-dependent term \cite{Schumayer2001}. This term can mimic the loss or
gain in the number of particles of the given species. We note here that usually
dissipation works against long-living coherent matter waves, however, the importance of
this imaginary potential term has been analysed in \cite{Brazhnyi2009, Rajendran2009}
and shown to permit exact soliton solution \cite{Nakkeeran2001}. In the context of BECs
at finite temperature, the interaction of the condensate with the thermal cloud could
also be taken into account as an imaginary term in the governing GP equation. This
interaction, due to its stochastic nature, can influence the dynamics of the solitons,
via density-fluctuation.

In this paper we carry out an analysis on how the intra- and interspecies interaction
influences the dynamics of a binary mixture of Bose-Einstein condensates.  We select
out of the many possible systems of two component BECs the pairs $^{87}$Rb--$^{87}$Rb
(prepared in two distinct hyperfine states), $^{23}$Na--$^{87}$Rb, and $^7$Li--$^{39}$K.
In the former two systems there is a possibility to study stability of bright-dark and
dark-dark soliton pairs, while the last pair is capable to sustain bright-bright
and bright-dark excitations.

The organization of the paper is as follows. Section~\ref{sec:TheModel} defines a quasi
one-dimensional model derived from the general three-dimensional, coupled Gross-Pitaevskii
equations assuming cigar-like harmonic oscillator trap potential. In the first part of
section \ref{sec:StabilityTest} we shall perform a stability analysis based on coupled
soliton excitations, and on Eq.~(\ref{eq:ConditionForIntegrability}) of the P-test. In
the second part of section \ref{sec:StabilityTest}, possible new interesting modes
exhibited by the bright-bright solitons will be shown. Section \ref{sec:summary} is
devoted to a conclusion and the summary.

\section{The model \label{sec:TheModel}}

In ultracold gases the interaction between two particles can usually be well
described by a scalar parameter, the scattering length. For a two-component
Bose-Einstein condensate one has to introduce three, possibly different,
scattering lengths, characterising the intra-species interactions ($a_{11}$,
$a_{22}$) and the inter-species coupling ($a_{12} = a_{21}$).

In the mean-field approximation a two-component BEC is described by the
coupled Gross-Pitaevskii equations \cite{Gross1961, Pitaevskii1961} in
the form:
\begin{equation} \label{eq:coupled_GP_math_notation}
   i \hbar \, \frac{\partial}{\partial t} \Psi_{i}
   =
   \left \lbrack
      - \frac{\hbar^{2}}{2m_{i}} \, \Delta
      + \sum_{j=1}^{2}{\Omega_{ij} \abs{\Psi_{j}}^{2}}
      + V_{i}
   \right \rbrack \! \Psi_{i},
\end{equation}
where $m_{i}$ denotes the individual mass of the $i$th atomic species, $\Omega_{ij}
= 2 \pi \hbar^{2} a_{ij} /\mu_{ij}$ with $a_{ij}$ being the 3D scattering length,
$\mu_{ij}=m_{i} m_{j}/(m_{i} + m_{j})$ is the reduced mass, and $V_{i}$ denotes
the external trapping potential. In the following, indices $i$ and $j$ label the
components, therefore, take only two values, 1 and 2. In the case of real trap
potentials the normalisation of the wave functions reads as $N_{i} = \int{\!\!
\abs{\Psi_{i}}^{2} \! dV}$ with $N_{i}$ denoting the number of atoms in the $i$th
component. We exclude those cases from our analysis where the species can transform
into each other, therefore the number of atoms in each component hereafter is
conserved.

\subsection{Transformed equations \label{subsec:TransformedEquations}}

If the three dimensional quadratic trapping potential is weak in one direction, i.e.
\begin{equation} \label{eq:3D-HO}
   V_{i}
   =
   \frac{1}{2} m_{i}
   \Bigl \lbrack
      \omega_{i, x}^{2}      x^{2} +
      \omega_{i, \perp}^{2} \left ( y^{2} + z^{2} \right )
   \Bigr \rbrack,
\end{equation}
where $\omega_{ix} \ll \omega_{i\bot}$, one may replace the three-dimensional equations
(\ref{eq:coupled_GP_math_notation}) with a coupled system of quasi one-dimensional GP
equations. Although Eq.~(\ref{eq:coupled_GP_math_notation}) is a nonlinear equation,
physically we may assume that the weak $x$-direction decouples from the strong [$yz$]--plane,
therefore the macroscopic wave functions can be written as
\begin{equation} \label{eq:1D-wf}
   \Psi_{i} ({\bf r},t)
   =
   \sqrt{N_{1}}\,
   \psi_{i} (x,t)\, \chi_{i, \perp} (y,z,t)
\end{equation}
 with
$\chi_{i, \perp}$ represents the ground-state solution of the corresponding two-dimensional
Schr{\"o}dinger equation in the [$yz$]--plane. The external potential introduces suitable
units of length and time as $a_{\perp} = \sqrt{ {\hbar}/{m_{1} \omega_{1, \perp}}}$ and
$\tau = 1 / \omega_{1, \perp}$, respectively. By rescaling the spatial and temporal variable
with $a_{\perp}$ and $\tau$ one obtains two quasi one-dimensional GP equations
\begin{subequations}
   \begin{eqnarray}
      \label{eq:coupled_GP1}
      i {\psi}_{1, t}
      &=&
      \!\!
      \left \lbrack
         - \frac{1}{2} \, \partial_{xx}
         +  \frac{\lambda_1^2}{2} \,  x^2
         + {b_{11} \abs{{\psi}_{1}}^{2}}
         + {b_{12} \abs{{\psi}_{2}}^{2}}
      \right \rbrack \! {\psi}_{1},
      \hspace*{7mm}
      \\
   \label{eq:coupled_GP2}
      i {\psi}_{2,t}
      &=&
      \!\!
      \left \lbrack
         - \frac{\kappa}{2} \, \partial_{xx}
         +  \frac{\lambda_2^2}{2\kappa} \,  x^2
         + {b_{21} \abs{{\psi}_{1}}^{2}}
         + {b_{22} \abs{{\psi}_{2}}^{2}}
      \right \rbrack \! {\psi}_{2},
   \end{eqnarray}
\end{subequations}
where $b_{11} = 2a_{11} N_{1}$, $b_{22} = 2 a_{22} N_{1} \kappa/
\gamma$, $b_{12} = b_{21} = 2 a_{12} N_{1} (1 + \kappa)/(1 +
\gamma)$, $\gamma = \omega_{2, \perp}/ \omega_{1, \perp}$, $\kappa
= m_{1}/m_{2}$, $\lambda_{1} = \omega_{1,x} / \omega_{1, \perp}$,
$\lambda_{2} = \omega_{2,x} / \omega_{1, \perp}$. The
normalisation is such that $\int{\abs{\psi_{1}}^{2}dx} = 1$ and
$\int{\abs{\psi_{2}}^{2}dx} = N_{2}/N_{1}$. Moreover, the relation
$\gamma^{2} = \kappa$ must hold if both species experience the
same harmonic potential. (Note the slight departure from
Ref.~\cite{Liu2009} in the definition of $b_{22}$ which, however,
may result in large difference of values of $b_{22}$ if a
two-component condensate contains species with different masses
$m_{1} \ne m_{2}$.)

\subsection{Thomas-Fermi background \label{sec:TFBackground}}

If the kinetic energy term is negligible compared to the potential
energy terms in (\ref{eq:coupled_GP1}-b), then one may apply the
Thomas-Fermi approximation to determine the density distribution
of the ground state. Following~\cite{Schumayer2004} we write the
corresponding wave-functions as
\begin{equation} \label{eq:TF-1}
   \psi_i(x,t)
   \approx
   \Phi^{\mathrm{TF}}_{i} (x)
   \exp{(-i E_{i}^{\mathrm{TF}} t/ \hbar)},
\end{equation}
resulting in TF densities
\begin{equation} \label{eq:TF-2}
   \abs{\Phi^{\mathrm{TF}}_{i}}^{2}
   =
   \frac{A_{i}}{\Delta}
   \left ( x_{i}^{2} - x^{2} \right ),
   \quad
   (\abs{x} < x_{i})
\end{equation}
and TF energies
\begin{subequations}
   \begin{eqnarray}
      \label{eq:TF-E1}
      E_{1}^{\mathrm{TF}}
      &=&
      \left ( b_{11} A_{1} x_{1}^{2} + b_{12} A_{2} x_{2}^{2} \right ) / \Delta,
      \\
      \label{eq:TF-E2}
      E_{2}^{\mathrm{TF}}
      &=&
      \left ( b_{12} A_{1} x_{1}^{2} + b_{22} A_{2} x_{2}^{2} \right ) / \Delta,
   \end{eqnarray}
\end{subequations}
where $\Delta = b_{11} b_{22}-b_{12}^{2}$. The parameters $A_i$ and $x_{i}$ represent
the amplitude of the density and the extension of the condensates, respectively. All
these quantities are determined by the system parameters $b_{ij}$, $N_{i}$, $\lambda_{i}$
according to the following relations:
\begin{equation} \label{eq:A-1}
   A_{1} = \frac{b_{22}}{2} \lambda_{1}^{2} - \frac{b_{12}}{2\kappa} \lambda_{2}^{2},
   \quad \mbox{and} \quad
   A_{2} = \frac{b_{11}}{2 \kappa} \lambda_{2}^{2} - \frac{b_{12}}{2} \lambda_1^2,
\end{equation}
while the extensions are
\begin{equation} \label{eq:x-1}
   x_{1} = \left ( \frac{3}{4} \frac{\Delta}{A_{1}} \right )^{1/3}
   \quad \mbox{and} \quad
   x_{2} = \left ( \frac{3}{4} \frac{\Delta}{A_{2}} \frac{N_{2}}{N_{1}} \right )^{1/3}.
\end{equation}
Although the Thomas-Fermi density distribution is not physical at $x=x_{i}$,
it still provides a good starting point for analytical calculations. In our numerical
treatment we will not use this approximation, rather start our simulations from the
appropriate ground state solution of the one-dimensional GP equation.

\subsection{Coupled soliton excitations \label{sec:CoupledSolitons}}

Now we are seeking solutions of Eqs.~(\ref{eq:coupled_GP1}-b) which support
soliton excitations. A static soliton excitation can be written as
\cite{Schumayer2004}
\begin{equation} \label{sol-excit}
   {\widetilde{\psi}}_{i} (x,t)
   =
   \Phi^{\mathrm{TF}}_i(0)
   \varphi_{i} (x)
   \exp{(-i {\widetilde{E}}_{i} t)}.
\end{equation}
By inserting this ansatz into the GP Eqs.~(\ref{eq:coupled_GP1}-b) and
neglecting the small potential contributions one obtains the coupled soliton
equations as follows
\begin{subequations}
   \begin{eqnarray}
      \label{eq:coupled_SOL1}
      \widetilde{E}_{1} {\varphi}_{1}
      &=&
      \left \lbrack
         - \frac{1}{2} \, \partial_{xx}
         + \widetilde{b}_{11} \abs{{\varphi}_{1}}^{2}
         + \widetilde{b}_{12} \abs{{\varphi}_{2}}^{2}
      \right \rbrack \! {\varphi}_{1}
      \\
      \label{eq:coupled_SOL2}
      \widetilde{E}_2 {\varphi}_{2}
      &=&
      \left \lbrack
       - \frac{\kappa}{2} \, \partial_{xx}
       + \widetilde{b}_{21} \abs{{\varphi}_{1}}^{2}
       + \widetilde{b}_{22} \abs{{\varphi}_{2}}^{2}
      \right \rbrack \! {\varphi}_{2}
   \end{eqnarray}
\end{subequations}
with $\widetilde{b}_{ij} = b_{ij} A_{j} x_{j}^{2} / \Delta$. The normalisation
of the soliton solutions reads as follows
\begin{subequations}
   \begin{eqnarray}
      \label{eq:norm_SOL1}
      \int_{-L_{1}}^{L_{1}}{\abs{\varphi_{1}}^{2} dx}
      &=&
      \frac{\Delta}{A_{1} x_{1}^{2}},
      \\
      \label{eq:norm_SOL2}
      \int_{-L_{2}}^{L_{2}}{\abs{\varphi_{2}}^{2} dx}
      &=&
      \frac{\Delta}{A_{2} x_{2}^{2}}
      \frac{N_{2}}{N_{1}},
   \end{eqnarray}
\end{subequations}
where the integrations, in both cases, are over the spatial extension,
$L_{i}$, of the solitons.

The above coupled equations admit generic moving soliton solutions
of the types: bright-bright (BB), bright-dark (BD), dark-bright (DB)
and dark-dark (DD). We shall investigate here a simple static bright-dark
soliton pair solution by taking the first component to be a static bright
soliton
\begin{subequations}
   \begin{equation} \label{eq:_BSOL1}
      \varphi_{1}^{\mathrm{BD}} (x) = q_{1} \sech{(k_{1} x)},
      \quad
      \varphi_{1}^{\mathrm{BD}}(x \rightarrow \pm \infty) = 0,
      \quad
   \end{equation}
and the second component to be the static dark soliton
   \begin{equation} \label{eq:_BSOL2}
      \varphi_{2}^{\mathrm{BD}} (x) = q_{2} \tanh{(k_{2} x)},
      \quad
      \varphi_{2}^{\mathrm{BD}}(x \rightarrow \pm \infty) = \pm q_{2}
   \end{equation}
\end{subequations}
with yet unknown complex amplitudes $q_{i}$ and wave-vector, $k_{i}$. The
latter one is related to the width of the soliton $k_{i} \sim 1/L_{i}$.
By inserting the above ansatz into the Eqs.~(\ref{eq:coupled_SOL1}-b)
and equating the coefficients of the constant and $x$-dependent terms,
respectively, one may conclude that the wave-vectors of the dark and
bright solitons must be equal, $k_{1} = k_{2} \equiv k$. The amplitudes
are expressed by the system parameters as
\begin{subequations}
   \begin{eqnarray}
      \label{eq:_q1}
      \abs{q_{1}}^{2} = \frac{k^{2}}{A_{1} x_{1}^{2}} (\kappa b_{12} - b_{22}),
      \\
      \label{eq:_q2}
      \abs{q_{2}}^{2} = \frac{k^{2}}{A_{2} x_{2}^{2}} (\kappa b_{11} - b_{12}).
   \end{eqnarray}
\end{subequations}
The energy of these excitations read as
\begin{subequations}
   \begin{eqnarray}
      \label{eq:_tildeE1BD}
      \widetilde{E}_{1}^{\mathrm{BD}}
      &=&
      k^{2} \left (
               \frac{\kappa b_{11} - b_{12}}{\Delta} b_{12} - \frac{1}{2}
            \right ),
      \\
      \label{eq:_tildeE2BD}
      \widetilde{E}_{2}^{\mathrm{BD}}
      &=&
      k^{2} \, \frac{\kappa b_{11} - b_{12}}{\Delta}b_{22}.
   \end{eqnarray}
\end{subequations}

Suppose now that the width parameter $k$ is a real number. The
modulus of the amplitudes of the bright-dark soliton pair superimposed
on the Thomas-Fermi background must be real numbers, thus one obtains
a the following set of conditions for the existence of this bright-dark
soliton excitations
\begin{eqnarray}
      \label{eq:_BD1cond}
      C_{1} \equiv \frac{\kappa b_{12} - b_{22}}{\Delta} \ge 0
      \quad \mbox{and} \quad
      C_{2} \equiv \frac{\kappa b_{11} - b_{12}}{\Delta} \ge 0.
      \qquad
   \end{eqnarray}

One may apply the same method to generate static bright-bright or
dark-dark soliton excitations. In the bright-bright case one obtains
the following solutions
\begin{subequations}
   \begin{eqnarray}
      \label{eq:_bb1}
      \varphi_{1}^{\mathrm{BB}} (x)
      &=&
      \frac{k}{\sqrt{A_{1}} x_{1}} \sqrt{\kappa b_{12} - b_{22}}\, \sech{(kx)},
      \\
      \label{eq:_bb2}
      \varphi_{2}^{\mathrm{BB}} (x)
      &=&
      \frac{k}{\sqrt{A_{2}} x_{2}} \sqrt{b_{12} - \kappa b_{11}}\, \sech{(kx)},
   \end{eqnarray}
\end{subequations}
with the conditions
\begin{equation}
   \label{eq:_BBcond}
   C_1\ge 0 \quad \textrm{and} \quad C_2\le 0,
\end{equation}
while the energies are $\widetilde{E}_{1}^{\mathrm{BB}} =
\widetilde{E}_{2}^{\mathrm{BB}}/\kappa  = -k^{2}/2$. The
dark-dark coupled soliton solutions read as follows
\begin{subequations}
   \begin{eqnarray}
      \varphi_{1}^{\mathrm{DD}} (x)
      &=&
      \frac{k}{\sqrt{A_{1}} \, x_{1}} \sqrt{b_{22} - \kappa b_{12}} \, \tanh{(kx)},
      \\
      \varphi_{2}^{\mathrm{DD}} (x)
      &=&
      \frac{k}{\sqrt{A_{2}} \, x_{2}} \sqrt{\kappa b_{11} - b_{12}} \, \tanh{(kx)},
   \end{eqnarray}
\end{subequations}
with the conditions
\begin{equation}
   \label{eq:_DDcond}
   C_{1} \le 0
   \quad \mathrm{and} \quad
   C_{2} \ge 0,
\end{equation}
and energies $\widetilde{E}_{1}^{\mathrm{DD}} =
\widetilde{E}_{2}^{\mathrm{DD}}/\kappa = k^{2}$.
Our interesting result shows that in the bright-bright and dark-dark
cases the energies are uniquely determined by the wave-vector and the
mass ratio. Note that the existence conditions (\ref{eq:_BD1cond}),
(\ref{eq:_BBcond}), and (\ref{eq:_DDcond}) are just the same as obtained
in Ref. \cite{Liu2009} for the existence of moving soliton pairs, while
the constraints for static excitations were published in \cite{Schumayer2004}.

\section{Stability tests by simulation \label{sec:StabilityTest}}

Below we are going to numerically investigate the stability of soliton
pairs. To solve the time dependent coupled Gross-Pitaevskii equations
(\ref{eq:coupled_GP1}-b), a third-order accurate split-step Fourier
transform method is used as described in Ref. \cite{Javanainen2006} for
a single-component condensate. Here we  solve the time independent
coupled Gross-Pitaevskii equations for their numerically exact ground
states using imaginary time method \cite{Lehtovaara2007} combined with
the split-step operator technique. Choosing initial distributions is
necessary to this method, and the Thomas-Fermi approximate solution
proved to be an effective initial guess for this purpose.

The procedure explained in the previous section can be generalised
for solitons moving with velocity $v$. Such a solution is given by
\begin{widetext}
\begin{subequations}
   \begin{eqnarray}
      \widetilde{\psi}_{1}^{\mathrm{BD}}
      &=&
      \Phi^{\mathrm{TF}}_1 (0)
      \varphi_{1}^{\mathrm{BD}} (x-vt)\,
      \exp{\left(-i {\widetilde{E}}_{1}^{\mathrm{BD}} t\right)}
      \exp{\!\left \{
                    -i \left \lbrack
                          v^{2} t \left ( \frac{b_{12} C_{2}}{\kappa^{2}} - \frac{1}{2} \right ) -
                          v (x - vt)
                       \right \rbrack
          \right \},
          }
      \hspace*{5mm}
      \\
      \widetilde{\psi}_{2}^{\mathrm{BD}}
      &=&
      \left \lbrack
         i \frac{\sqrt{C_2}}{\kappa}\,v + \Phi^{\mathrm{TF}}_2 (0) \varphi_{2}^{\mathrm{BD}} (x-vt)
      \right \rbrack \,
      \exp{\left(-i {\widetilde{E}}_{2}^{\mathrm{BD}} t\right)}
      \exp{\!\left (
                    -i \,\frac{b_{22} C_{2}}{\kappa^{2}} v^{2} t
          \right )
          }.
   \end{eqnarray}
\end{subequations}
\end{widetext}
At $v=0$ we obtain  the static soliton excitation solution as given by
\eqref{sol-excit} in the preceding section.

The existence conditions (\ref{eq:_BD1cond}) prescribe various relations
between domains of the inter- and intra coupling strengths, $b_{ij}$, where
it is possible to create bright-dark soliton pairs. These domains are listed
in Table \ref{tab:ExistenceOfBDSolitons}. Although the creation of a bright
soliton is generally associated with attractive interaction ($b_{ii} < 0$)
among the particles, Table \ref{tab:ExistenceOfBDSolitons} clearly shows
that, due to the appropriate other couplings, it may also be possible to
create a bright soliton in case of repulsive interaction $(b_{ii}>0)$. Such
a situation occurs in case 4a of Tab.~\ref{tab:ExistenceOfBDSolitons} which
we are going to analyse below.

\subsection{Moving bright-dark soliton pairs \label{sec:MovingBDSolitonPairs} }

Let us now investigate the stability of a bright-dark soliton pair for
the case 4a of Table \ref{tab:ExistenceOfBDSolitons} by considering two
experimentally accessible two-component BECs. The first is composed of
the two hyper-fine states of $^{87}$Rb atoms, the second is obtained
from $^{23}$Na and $^{87}$Rb atoms. In both cases the scattering lengths
are well known and can be tuned over a broad limit by using the Feshbach
resonance method. Our aim here is to explore the sensitivity of the temporal
evolution of soliton pair when the inter-atomic coupling strength, $b_{12}$,
is varied around the value prescribed by the ratio formula
(\ref{eq:ConditionForIntegrability}) obtained by performing
a Painlev{\'e}-analysis of the coupled GP equations \cite{Schumayer2001}.
In this respect we fix the intra-atomic strengths, $b_{ii}$, to values
easily accessible to the experiments and vary the inter-atomic interaction
values, $b_{12}$, within a small range allowed by the case 4a listed in
Table~\ref{tab:ExistenceOfBDSolitons}. In order to represent a more realistic
situation, a small velocity of $v=0.04$ is given the solitons, a weak harmonic
trapping potential is added along the longitudinal direction and the  solitons
are superimposed on the ground state density distribution. This procedure
spatially confines the solitons without affecting their essential stability
properties \cite{Liu2009}.

\begin{table}[h]
   \begin{tabular}{r<{\hspace*{1mm}}|p{5mm}p{5mm}cc}
      Case & $b_{11}$ & $b_{22}$ & Constraint on $b_{12}$ \\
      \hline\hline
       1   & $+$ & $-$ & no bright-dark soliton pair\\
       2   & $-$ & $+$ & $b_{11}\kappa<b_{12}<b_{22}/\kappa$\\
       3a  & $-$ & $-$ & $-\sqrt{b_{11}b_{22}}<b_{12}<b_{11}\kappa$ & if $\kappa^2 b_{11}>b_{22} $\\
       3b  &      &      & $b_{11}\kappa<b_{12}<-\sqrt{b_{11}b_{22}}$ & if $\kappa^2 b_{11}<b_{22} $\\
       4a  & $+$ & $+$ & $\sqrt{b_{11}b_{22}}<b_{12}<b_{22}/\kappa$ & if $\kappa^2 b_{11}<b_{22} $\\
       4b  &      &      & $b_{22}/\kappa<b_{12}<\sqrt{b_{11}b_{22}}$ & if $\kappa^2 b_{11}>b_{22} $\\
   \end{tabular}
   \caption{\label{tab:ExistenceOfBDSolitons}%
            Various domains of inter- and intra-atomic interaction strengths,
            $b_{ij}$, permitting the existence of a bright-dark soliton pair.}
\end{table}

\begin{figure}[hbt!]
   \includegraphics[width=\textwidth/2-3mm]{./FIGURES/BD/Figure1a.eps}
   \includegraphics[width=\textwidth/2-3mm]{./FIGURES/BD/Figure1b.eps}
   \caption{\label{fig:RbRb}%
            (Color online) Evolution of the bright-dark solitons of two-component
            BEC composed of two hyperfine states of the $^{87}$Rb atom. Intraspecies
            scattering length $a_{12}=5.5$\,nm ($b_{12}=13.6$) satisfying ratio
            formula (\ref{eq:ConditionForIntegrability}) (top two panels), $a_{12}=5.4$\,nm
            ($b_{12}=13.35$) disobeying the ratio formula (\ref{eq:ConditionForIntegrability})
            (bottom two panels). Other parameters $a_{11}=0.8 \times 5.5$\,nm ($b_{11}=10.88$),
            $a_{22}=1.2 \times 5.5$\,nm ($b_{22}=16.32$), $m_{1} = m_{2} = 87$\,au,
            $\omega_{1, \perp}=2\pi\times 710$ Hz, $\lambda_{1} = \lambda_{2}=0.2$,
            $v = 0.04$, $N_{1} = 500$, $N_{2}=6600$. The snapshots are depicting the
            soliton pairs at $t=20$ and $t=200$, respectively.}
\end{figure}

\begin{figure}[hbt!]
   \includegraphics[width=\textwidth/2-4mm]{./FIGURES/BD/Figure2a.eps}
   \includegraphics[width=\textwidth/2-4mm]{./FIGURES/BD/Figure2b.eps}
   \caption{\label{fig:NaRb}%
            (Color online) Evolution of a bright-dark soliton pair in a two-component
            BEC composed of $^{23}$Na and $^{87}$Rb atoms. Intraspecies scattering
            length $a_{12}= 6.831$\,nm ($b_{12}=7.25$) satisfies formula
            (\ref{eq:ConditionForIntegrability}) (top two panels), while
            $a_{12}= 3.769$\,nm ($b_{12}=4$) deviates from Eq.~(\ref{eq:ConditionForIntegrability})
            (bottom two panels). Other parameters are: $a_{11}=2.7$\,nm $(b_{11}=3.4)$,
            $a_{22}= 5.5$\,nm $(b_{22}=3.6)$, $m_{1}=23$\,au, $m_{2}=87$\,au,
            $\omega_{1, \perp}=2 \pi \times 710$\,Hz, $\lambda_{1} = \lambda_{2}=0.25$,
            $v=0.15$, $N_{1}=500$, $N_{2}=58000$. The snapshots are taken at $t=20$,
            and $t=2$, respectively.}
\end{figure}

In Fig.~\ref{fig:RbRb}, we plotted a bright-dark soliton pair composed of
two hyperfine states of $^{87}$Rb atoms. The top two panels show the case
when the inter-species interaction is chosen according to
Eq.~(\ref{eq:ConditionForIntegrability}). Both solitons oscillate in the
harmonic trap in a stable manner with an angular frequency slightly less
than $\omega_{x}/\sqrt{2}$ what the one-dimensional Thomas-Fermi model
predicts for a single dark soliton \cite{Busch2000}. The difference is
probably caused by the presence of the bright soliton, since the bright
component fills the dip of the dark soliton, therefore the dark soliton
has to drag this extra mass as well. This effect has recently been observed
\cite{Becker2008} with a $^{87}$Rb-$^{87}$Rb condensate, prepared in the
$\vert F=2, m_{F} = 0 \rangle$ (bright soliton) and  $\vert F=1, m_{F} = 0
\rangle$ (dark soliton) hyperfine states.

However, if $a_{12}$ (or equivalently $b_{12}$) is tuned away from this
particular value, the stability is lost, and the initial forms of the
solitons are destroyed by the destructive interference of the constantly
emitted and re-captured sound waves. It is worthwhile to mention, although
only as a qualitative statement, that the appearance of sound waves made
the solitons' oscillation faster (see the different range of time in the
top two and bottom two panels of Fig.~\ref{fig:RbRb}). The sound waves
travel faster than the solitons, and after being reflected back from the
edge of the condensate, they collide with the solitons. The subsequent
collisions possibly speed up the oscillation and turn it into an irregular
sloshing. Despite the irregular movement of the solitons, the dark component
still captures the bright soliton during the motion.

In Fig.~\ref{fig:NaRb} the evolution of a bright-dark soliton pair in binary
BEC composed of $^{23}$Na and $^{87}$Rb atoms is plotted with corresponding
snapshots of the density distribution. Upper two panels exhibit the bright
($^{23}$Na) and dark ($^{87}$Rb) excitations, respectively, which are stable
for a long period due to the fine tuning of scattering lengths $a_{ij}\,
(i,j=1,2)$ which satisfy the P-test formula \eqref{eq:ConditionForIntegrability}.
On the lower two panels of Fig.~\ref{fig:NaRb} we see, however, that the
initial soliton excitations do not remain stable but are rapidly dissolved
due to detuning the interspecies scattering length $a_{12}$ from the value
obeying the integrability condition expressed by Eq. \eqref{eq:ConditionForIntegrability}.

Summarising, both examples exhibited in Figs.~\ref{fig:RbRb} and \ref{fig:NaRb}
show that a long lived stability of the bright-dark soliton pairs can be achieved
only in case if the interatomic interaction parameters $a_{12}$ (or $b_{12}$)
is chosen according to Eq.~\eqref{eq:ConditionForIntegrability}. This observation
emphasises that there may be situations in the creation of two-component BECs
when fine tuning of scattering lengths according to the P-test formula
(\ref{eq:ConditionForIntegrability}) may prove useful.

At the end of this section we briefly mention an early investigation of the
stability of a heterogeneous two-component Bose-Einstein condensate by Law
et al. \cite{Law1997}. This work carried out a linear stability analysis by
calculating the lowest eigenvalue of an excitation. In this approach the
appearance of a negative eigenvalue signals instability. It was shown for
a sodium-rubidium system, just as above, that stability occurs only in finite
range of the inter-species scattering length. The direct quantitative comparison
with our work, however, is less straightforward because Law et al.'s analysis
assumes a spherically symmetric condensate, while we analyse a quasi one-dimensional
model.

\subsection{Static bright-bright soliton pairs\label{sec:StaticBBSolitonPairs}}

Let us now investigate the temporal stability of a static ($v=0$) bright-bright
soliton pair obtained as the exact solution of Eqs.~(\ref{eq:coupled_GP1}-b) in
the absence of a trapping potential ($\lambda_{i}=0$). The solutions have the
form as in Eq.~(\ref{sol-excit}) where $\varphi_{i}(x)$ are chosen to have the
functional forms given in Eqs.~(\ref{eq:_bb1}-b).

The existence conditions (\ref{eq:_BBcond}) prescribe the relations between
domains of the inter- and intra-species coupling strengths, $b_{ij}$, which
are listed in Table \ref{tab:ExistenceOfBBSolitons}. Moreover, the common
wave-vector is $k = 1/ (2 C_{1})$.

\begin{table}[h]
   \begin{tabular}{r<{\hspace*{1mm}}|p{5mm}p{5mm}cc}
      Case & $b_{11}$ & $b_{22}$ & Constraint on $b_{12}$ \\
      \hline\hline
      1    & $+$ & $-$ & $b_{12}< b_{22}/\kappa$ \\
      2    & $-$ & $+$ & $b_{12}< b_{11}/\kappa$ \\
      3a   & $-$ & $-$ & $b_{22}/\kappa < b_{12}<\sqrt{b_{11}b_{22}}$ & if $\kappa^2 b_{11}<b_{22} $\\
      3b   &      &      & $b_{12}<b_{11}\kappa$ & if $\kappa^2 b_{11}<b_{22} $ \\
      3c   &      &      & $b_{11}\kappa < b_{12}<\sqrt{b_{11}b_{22}}$& if $\kappa^2 b_{11}>b_{22} $\\
      3d   &      &      & $b_{12}<b_{22}/\kappa$ & if $\kappa^2 b_{11}>b_{22} $ \\
      4    & $+$ & $+$ & $b_{12}<-\sqrt{b_{11}b_{22}}$\\
   \end{tabular}
   \caption{\label{tab:ExistenceOfBBSolitons}
            Various domains of inter- and intra-atomic interaction strengths,
            $b_{ij}$, permitting the existence of a bright-bright soliton pair.
           }
\end{table}

\begin{figure}[htb]
   \includegraphics[width=\textwidth/2-3.7mm]{./FIGURES/BB/Figure3a.eps}
   \includegraphics[width=\textwidth/2-3.7mm]{./FIGURES/BB/Figure3b.eps}
   \includegraphics[width=\textwidth/2-3.7mm]{./FIGURES/BB/Figure3c.eps}
   \caption{\label{fig:BrightBright}%
            (Color online) Evolution of a bright-bright soliton pair in a binary BEC
            composed of $^{7}$Li and $^{39}$K atoms. Intraspecies scattering lengths
            and atom numbers are $a_{12}=0.2$\,nm ($b_{12}=0.46$), $N_{2}=2029$ (first
            two panels), $a_{12}= 0.275$\,nm ($b_{12}=0.64$) $N_{2}=2270$ (second pair
            of the panels), $a_{12}= 0.3$\,nm ($b_{12}=0.7$) $N_2=2347$ (third pair of
            the panels). Other parameters are: $a_{11}=-1.4$\,nm $(b_{11}=-3.93)$,
            $a_{22}= -0.9$\,nm $(b_{22}=-1.07)$, $m_{1}=7$\,au, $m_{2}=39$\,au,
            $\omega_{1, \perp}=2\pi\times 710$\,Hz, $\lambda_{1} = \lambda_{2}=0$,
            $v=0$, $N_{1}=2000$. The snapshots are taken at $t=35$, $t=35$, and
            $t=20$, respectively.}
\end{figure}
\clearpage

Scenario 4, for example, describes two condensates for which the intra-species
interactions are repulsive. This situation, using Hartree-Fock calculation, has
been theoretically analysed \cite{Esry1997} soon after the observation of overlapping
condensates prepared from two hyperfine states of $^{87}$Rb \cite{Myatt1997}.
For this case, i.e. $m_{1}=m_{2}$, it was established that the two condensates
cannot co-exist if $\abs{a_{12}} < \sqrt{a_{11}a_{22}}$. Our approach reproduces
and extends this result for the case of different species. This surprisingly simple
relation can be understood using energetic arguments; if $b_{12}$ overcomes the
geometric mean of the intra-species interaction strengths, the repulsion between
the two condensates will separate the two condensate completely and they will
not overlap any more. However, if the two species attract each other enough,
i.e. $b_{12}< -\sqrt{b_{11}b_{22}}$, the attraction will dominate and can counteract
the individual repulsion present in each component.

Another interesting scenario here is the one listed under the case 3a in
Table \ref{tab:ExistenceOfBBSolitons} showing that it is possible to create
a bright-bright pair, within the range $0 < b_{12} < \sqrt{b_{11} b_{22}}$,
in spite of the repulsive inter-atomic interaction. In order to investigate
the stability of the bright-bright soliton pair in this domain we simulate
the temporal evolution of the BEC system composed of $^7$Li and $^{39}$K
atoms accessible for experiments. As Fig.~\ref{fig:BrightBright} shows,
depending on the value of the interspecies interaction, two types of
instability may occur. At values of $a_{12}$ less than a critical value of
about $0.275$\,nm, the two standing solitons begin to repel each other
(first two panels) and depart from each other as they were particles,
preserving the total zero momentum. Above the critical $a_{12}$ value
the soliton with constituents of the smaller mass splits into two equal
parts and the heavier component becomes a breather keeping its original
place (third pair of panels). At the critical value of $a_{12}$ both
types of instability can be observed (third and fourth panel)
simultaneously.

Similar instability effects of a bright-bright and other soliton pairs
have also been observed in Ref.~\cite{Beitia2010}. The instability is
explained by Kevrekidis at al.~\cite{Kevrekidis2005} as the appearance
of a negative eigenvalue pair obtained by a linear stability analysis.
It has also been shown that higher order bright-bright soliton pairs
would exhibit instability irrespectively of the system parameters.
Interestingly, it is possible for one of the higher-order bright solitons,
i.e., for the one which has the stronger self-attraction, to recover its
stability by collapsing into one bright soliton and expel the other
component from its original position.

Finally, it is important to note here, that the phases of the divided
solitons are equal, therefore these particle-like wave-packets remain
coherent with each other. Similarly to an optical beam splitter where
light is divided into two coherent beams, one could divide these matter
waves and use one of them as a probe and the other one as a control
packet. The probe packet could undergo transformations, while the control
packet is left to evolve freely, thereby a phase-difference could build
up between the two packets. If the two packets are brought together again,
the phase-difference could cause interference pattern which allows one to
quantify coherence. This may potentially be helpful for calibrational
purposes as well.

\subsection{Dark-dark soliton pairs \label{sec:DDSolitonPairs}}

At the end we are examining the third possible combination of soliton pairs;
a dark soliton is excited in each condensate. Interestingly this pairing can
be stable even if the interaction inside each condensate is attractive ($a_{11}$,
$a_{22} <0$). The dark-dark solitons are described by the following formulas
\begin{widetext}
\begin{subequations}
   \begin{eqnarray}
      \widetilde{\psi}_{1}^{\mathrm{DD}}
      &=&
      \left \lbrack
         i \sqrt{-C_{1}} \,v +
         \Phi^{\mathrm{TF}}_{1} (0) \varphi_{1}^{\mathrm{DD}} (x-vt)
      \right \rbrack \,
      \exp{ \left ( -i {\widetilde{E}}_{1}^{\mathrm{DD}} t \right ) }
      \exp{ \left \lbrack -i \left( \frac{b_{12} C_{2}}{\kappa^2} -
                          b_{11} C_{1} \right ) v^{2} t
            \right \rbrack},
      \\
      \widetilde{\psi}_{2}^{\mathrm{DD}}
      &=&
      \left [ i \frac{\sqrt{C_2}}{\kappa}\,v +
      \Phi^{\mathrm{TF}}_{2} (0)
      \varphi_{2}^{\mathrm{DD}}(x-vt)\right]\,
      \exp{\left(-i {\widetilde{E}}_{2}^{\mathrm{DD}} t\right)}
      \exp{\left[-i\left( \frac{b_{22} C_{2}}{\kappa^2}-b_{12} C_{1}
                   \right ) v^{2} t \right]}.
\end{eqnarray}
\end{subequations}
\end{widetext}
For our numerical investigations a $^{87}$Rb--$^{87}$Rb system has been
chosen, with repulsive intra-species interactions, $a_{11}=5.335$\,nm
and $a_{22} = 5.665$\,nm taken from \cite{Hall1998}. This choice of
coupling corresponds to scenario 4 of Table~\ref{tab:ExistenceOfDDSolitons}.
Tuning $a_{12}$ into the positive regime results in a stable dark-dark
soliton pairs, see top two panels of Fig.~\ref{fig:DarkDark}. In the
bright-bright case we have found that if the intra-species interactions
are attractive and the inter-species interactions are repulsive, the
solutions are unstable. We have tested numerically that this observation
is valid in the dark-dark case as well, if the sign of the interactions
is inverted. Furthermore, our numerical calculation provided an interesting
phenomenon in the dynamics of these solitons. Preparing the dark solitons
in the same way as before (i.e., by superimposing onto the ground state)
but reversing the sign of $a_{12}$, resulted not only in losing the
long-term stability, but also in the decay of a dark soliton in one
of the components, and emerging as a secondary dark soliton in the
other component (bottom two panels of Fig.~\ref{fig:DarkDark}).

\begin{figure}[hbt!]
   \includegraphics[width=\textwidth/2-3mm]{./FIGURES/DD/Figure4a.eps}
   \includegraphics[width=\textwidth/2-3mm]{./FIGURES/DD/Figure4b.eps}
   \caption{\label{fig:DarkDark}%
            (Color online) Evolution of a dark-dark soliton pair in a
            binary BEC composed of $^{87}$Rb and $^{87}$Rb atoms.
            Intraspecies scattering lengths are $a_{12}=1.5$\,nm
            ($b_{12}=3.62$) (top panel), $a_{12}= -1.5$\,nm ($b_{12}=
            -3.62$) (bottom panel). Other parameters are: $a_{11}=0.97
            \times 5.5$\,nm $(b_{11}=13.04)$, $a_{22}=1.03\times 5.5$\,nm
            $(b_{22}=13.52)$, $m_{1}=87$\,au, $m_{2}=87$\,au,
            $\omega_{1, \perp}=2\pi\times 710$\,Hz, $\lambda_{1} =
            \lambda_{2}=0.2$, $v=0.2$, $N_{1}=500$, $N_{2}=1600$.
            The snapshots are taken at $t=50$.
           }
\end{figure}

\begin{table}[h]
   \begin{tabular}{r<{\hspace*{1mm}}|p{5mm}p{5mm}cc}
      Case & $b_{11}$ & $b_{22}$ & Constraint on $b_{12}$ \\
      \hline\hline
      1    & $+$ & $-$ &                        $b_{12} > b_{11} \kappa$        \\
      2    & $-$ & $+$ &                        $b_{12} > b_{22}/\kappa$        \\
      3    & $-$ & $-$ &                        $b_{12} > \sqrt{b_{11}b_{22}} $ \\
      4a   & $+$ & $+$ & $-\sqrt{b_{11}b_{22}} < b_{12} < b_{22}/\kappa$ & if $\kappa^2 b_{11}>b_{22}$\\
      4b   &     &     &                        $b_{12} > b_{11} \kappa$ & if $\kappa^2 b_{11}>b_{22}$\\
      4c   &     &     & $-\sqrt{b_{11}b_{22}} < b_{12} < b_{11}\kappa$  & if $\kappa^2 b_{11}<b_{22}$\\
      4d   &     &     &                        $b_{12} > b_{22}/\kappa$ & if $\kappa^2 b_{11}<b_{22}$\\
\end{tabular}
   \caption{\label{tab:ExistenceOfDDSolitons}
            Various domains of inter- and intra-atomic interaction strengths,
            $b_{ij}$, permitting the existence of a dark-dark soliton pair.
           }
\end{table}

\section{Summary \label{sec:summary}}

We have considered the existence and stability of soliton excitations in a
two-component Bose-Einstein condensate both analytically and numerically.
Our model allows the components to represent different elements ($m_{1}
\ne m_{2}$), but we also included those cases when two hyperfine states of
the same element constitute the condensates. We excluded the possibility
that these components can transmute into each other, i.e. the hyperfine
states cannot be driven into each other. The dynamics of these condensates,
within the mean-field zero-temperature approximation, are governed by the
coupled Gross-Pitaevskii equations. We chose our presentation to be suitable
for combining analytical results with earlier numerical investigations,
e.g. \cite{Schumayer2001, Xun-Xu2010, Liu2009}.

The occurrence of particle-like excitations together with conserved quantities
are associated with the integrability of a nonlinear evolution field-equation,
such as the coupled Gross-Pitaevskii equations. Using the results of a recent
Painlev{\'e} analysis of the coupled Gross-Pitaevskii equation
\cite{Schumayer2001}, we showed how the system parameters determine the
integrability of this system. However, for the CGP equations the studies
so far restricted themselves for either equal coupling coefficients ($b_{11}=
b_{12}=b_{22}$) and/or equal masses. We note here that excitations with
long lifetime may exist for a nonlinear evolution equation even when the
integrability conditions are violated \cite{Novoa2008}, but these cases
are possibly exceptional and do not represent the generic behaviour.

We examined those ranges of system parameters which permit coupled soliton
solutions of the governing equations (\ref{eq:coupled_GP1}-b). Here, we
utilised the two well-known one-soliton solutions of the one-component
nonlinear Schr{\"o}dinger equation; the bright and dark solitons. These
excitations then spatially modify the Thomas-Fermi ground state in our
theoretical description or by the appropriate ground state in our numerical
simulation. We found analytically that each possible pair of these
solitons have a range of system parameters where they are stable, and
presented these ranges in tabulated form in Tabs.
\ref{tab:ExistenceOfBDSolitons}, \ref{tab:ExistenceOfBBSolitons},
and \ref{tab:ExistenceOfDDSolitons}. One can present these
findings differently, if the intra-species coupling coefficients
are held relatively fixed, and only the inter-species one is
tuned. For example, if we assume $b_{11}$ and $b_{22}$ to be
positive (effective repulsive interaction) then the three tables
\ref{tab:ExistenceOfBDSolitons}, \ref{tab:ExistenceOfBBSolitons},
and \ref{tab:ExistenceOfDDSolitons} can be combined into one
graph (see Fig.~\ref{fig:DifferentRepresentation}).

\begin{figure}[hbt!]
   \includegraphics[width=\textwidth/2-4mm]{./FIGURES/Figure5.eps}
   \caption{\label{fig:DifferentRepresentation}%
            Ranges of the inter-species coupling coefficient, $b_{12}$, in
            which different types of soliton pairs can exist and are stable.
            The top graph is valid if $\kappa^{2} b_{11} < b_{22}$, while
            the lower graph is for the complementary case, i.e. $\kappa^{2}
            b_{11} > b_{22}$. The dotted line shows the result of the
            Painlev{\'e}-test (\ref{eq:ConditionForIntegrability}) in
            the case of harmonic trapping potentials.
           }
\end{figure}

Moreover, we also examined how the stability of these soliton pairs is
changing if the inter-species coupling coefficient is detuned from the
value predicted by the Painle{\'e}-analysis via
Eq.~(\ref{eq:ConditionForIntegrability}). It was shown, irrespectively
of the pairing, that the stability is lost, although the pairs were not
equally sensitive to the detuning, e.g. the motion of the bright-dark
pair became erratic and preserved its periodicity only qualitatively
after changing $a_{12}$ from 5.5\,nm to 5.4\,nm.

Two types of instability of static bright-bright soliton pairs composed
of species with unequal masses have been observed. When the interspecies
interaction lies below a critical value ($0 < b_{12} < b_{12}^{\mathrm{cr}}$),
the static bright-bright pair evolves into a repulsive, momentum conserving,
moving soliton pair. When the value of $b_{12} > b_{12}^{cr}$ then one of
the bright soliton (the constituent with smaller mass) splits into two equal
portion of same phase while the other bright soliton becomes a breather.

Well below the critical temperature of the Bose-Einstein condensate our
description is expected to be adequate and the results could help
experimentalist to modify the scattering lengths via Feshbach resonance
into a range where stable soliton pairs exist. As the temperature increases,
however, the interaction with the thermal cloud becomes more and more
important, and could not be neglected any more. We have not yet examined
how the interplay between the condensates and the thermal clouds (for each
component) would modify the dynamics. This needs further research beyond
the mean field description.

{\emph{Acknowledgment}} ---  D. Schumayer acknowledges financial
support from NERF-UOOX0703 (NZ) and also by the University of
Otago. G. Csire and B. Apagyi thank DFG for subsidising their
stay at Institute of Theoretical Physics, University of Giessen.


\section{Introduction}
Cooling of atoms/ions is an important precondition for quantum logic and computation \cite{Wineland2005,sasura01}, wherein an initial ground state preparation of the vibrational state is required, for high resolution spectroscopy beyond the Doppler limitations and for the analysis of low temperature phenomena (e.g. BEC condensates). Cooling by light (laser cooling) for atoms and ions was originally proposed in \cite{haensch74,wineland75}.  For trapped ions, sideband cooling \cite{stenholm86,wineland78}  is the standard cooling method, achieved by tuning the laser frequency to the red sideband of the bare atomic frequency, which leads to a resonant cooling contribution whereas heating  is only induced by off-resonant carrier and blue-sideband excitations. To obtain well-resolved sidebands requires the linewidth of the atomic transition to be much smaller than the trapping frequency $\Gamma\ll\nu$ (strong confinement). Furthermore, taking into account that off-resonant heating excitations increase with increasing laser intensity   the coupling Rabi-frequency has to be smaller than the trapping frequency $\Omega\ll\nu$ (weak coupling) \cite{sasura01}. Cooling to lower temperatures was reached using the concept of Raman sideband cooling \cite{lindberg86}, applicable for atoms \cite{kerman00} and ions \cite{lee96,monroe95}, which allows the realisation of very small cooling linewidths by Raman two-photon transitions. Further improvements on the way to reach the ground state are achieved with dark state concepts based on eliminating off-resonant excitations via quantum interference.  In EIT cooling \cite{morigi00,roos00} the concept of electromagnetically induced transparency \cite{fleischhauer05} is used to eliminate the carrier transition, which gives the main heating contribution in the sideband schemes. It uses a three level system coupled by two lasers, wherein one of the lasers (the stronger one) creates the absorption spectrum,  a Fano like profile with a zero at the carrier transition and a peak at the red sideband transition. Therefore stronger couplings can be realised, no strong confinement is required due to the three level nature of the system and a large cooling bandwidth (i.e. frequency range for cooling) allows simultaneous cooling of several modes. An alternative approach, based on the Stark shift gate \cite{jonathan00}, is the Stark shift cooling \cite{retzker07}, cancelling the carrier transition for a correctly tuned Rabi frequency, which is of the order of the trapping frequency. This is achieved in the dressed state picture, i.e. taking the Stark shift of the laser interaction (originating from the zeroth order in the Lamb-Dicke parameter contributions of the resonant coupling laser) into account. The robust laser cooling scheme \cite{cerrillo10} combines the methods of EIT and Stark shift cooling resulting in a dark state for the combination of EIT and Stark shift Hamiltonian contribution and leading to an elimination of both the carrier transition and the blue sideband heating transition by interference.

The scheme presented here also eliminates the carrier and blue sideband transitions in lowest order in the Lamb-Dicke parameter, using an atomic three level system, a combinations of four laser beams and  applying a magnetic field gradient. The application of an inhomogeneous magnetic field to trapped ions resulting in a position dependent resonance shift  and mediating the coupling between internal and external degrees of freedom (e.g. allowing the coupling by microwaves with a non-vanishing effective Lamb-Dicke parameter) was discussed in \cite{mintert01,wunderlich03}.  For our method only the proper choice of the gradient strength and the phase relations of the laser beams   form the conditions for a pure red sideband interaction up to first order in the Lamb-Dicke parameter. \par
The paper is organized as follows: In section\,\ref{one} the cooling scheme is presented and the influence of the magnetic field gradient is analysed. It turns out, that in a Schrieffer-Wolff transformed picture, both the carrier as well as the blue sideband transition are eliminated by interference and thus cooling to the vibrational ground state in zeroth order is achieved. In section\,\ref{three} the resonance condition for the cooling process is derived. Moreover it is shown that the maximal cooling amplitude is limited by the trapping frequency, a fact that arises from the amplitudes of the dressed states, which are tuned to resonance. Section\,\ref{four} presents an analytic formula for the cooling rate obtained by the method of adiabatic elimination and valid in the weak coupling case. Indeed a comparison to numerical results shows a good agreement in the range of small Rabi frequencies. In section\,\ref{five} the dependence of the cooling rate on the spontaneous decay rate as well as on the magnitude of the Lamb-Dicke parameter is analysed. Furthermore the dependence of the cooling rate and final population on phase fluctuations of the coupling lasers is presented resulting in the conclusion that the scheme, albeit becoming slower, still provides fast cooling and low final populations even for relatively large deviations from the optimal phase relation. Section\,\ref{six} deals with the cooling of multiple vibrational modes, i.e. with multiple trapped ions. Only for the center of mass mode a complete elimination of heating terms (carrier and blue sideband) is possible; however the strength of the heating terms for other modes is smaller than in standard cooling schemes.

\section{Cooling scheme and Hamiltonian}\label{one}
The scheme presented here is based on a complete elimination of the blue sideband heating terms in a Schrieffer-Wolff transformed picture (see eq.\,(\ref{sw1}) and (\ref{sw2}) below) in lowest order in the Lamb-Dicke parameter, which is achieved by a magnetic gradient field together with a four laser coupling scheme. It requires an atomic three level system consisting of a (fast decaying) excited state $\ket{e}$ and a Zeeman split ground state $\ket{+1}$ and $\ket{-1}$ induced by the external position dependent magnetic field. Each of the Zeeman split levels is coupled by two lasers with Rabi frequency $\Omega$  and opposite Lamb-Dicke parameters $\eta$ ($\eta=k\,\sqrt{\hbar/(2\,m\,\nu)}\,\cos\theta$, with $\theta$ the angle of the laser beam to the trap axis and $k$ the absolute value of the wavevector), obtained by irradiating the beams from opposite directions, to the excited level $|e\rangle$ having a phase difference $\varphi=\pi/2$ (see fig.\ref{fsetup}).\\
\begin{figure*}
\includegraphics[scale=0.6]{bothschemes.png}
\caption{\label{fsetup} Cooling scheme: a) Coupling scheme in the standard basis. Each of the Zeeman split ground state levels is coupled to the excited level $\ket{e}$ by two lasers with Rabi frequency $\Omega$, Lamb Dicke parameter $\eta$, phase $\varphi$ and detuning $\delta$.  Note that the external magnetic field gives, beside the position independent $\omega_0$ splitting, an additional position dependent energy shift via the magnetic gradient contribution. b) Scheme in dressed state basis after Schrieffer-Wolff transformation. See also eq.\,(\ref{ham4}).  }
\end{figure*}
The ion is assumed to be trapped in an harmonic trapping potential with frequency $\nu$ and for simplicity only the motion along the trap axis is taken into account. Moreover, the expansion up to first order in terms of the Lamb-Dicke parameter requires the atom/ion to be localized on the scale provided by the laser wavelength (Lamb-Dicke limit, $\eta \ll 1$).
The position dependent magnetic field $B(x)$ has the form
\begin{equation} B(x) \approx B_0+\left.\dfrac{\mathrm{d}B}{\mathrm{d}x}\right|_{x=0}\,x    \end{equation}
and higher terms are not present or can be neglected in the following expansion in terms of the effective Lamb-Dicke parameter respectively.
The corresponding energy shift due to the magnetic Zeeman splitting is (with the position operator $x=\sqrt{\frac{\hbar}{2\,m\,\nu}}(b+b^\dagger)$, wherein $b^\dagger$ and $b$ are the phonon creator and annihilator operators respectively) given by
\begin{equation}\label{en_split} \Delta E=-\mu\,B=\hbar\,\omega_0\,m_s+\eta_{eff}\,\hbar\,\nu\,\left(b+b^\dagger\right)\,m_s   \end{equation}
with $\mu$ the magnetic moment, $m_s$  the magnetic quantum number and $\eta_{eff}=\frac{\sqrt{\hbar}}{\sqrt{2\,m\,\nu}}\frac{\omega_0}{\nu B_0}\,\frac{\partial B}{\partial x}$, providing a coupling between the spatial position and the internal energy, i.e. phonons and internal atomic states. The effective Lamb Dicke parameter has a physical meaning as the energy difference between the left and right bound of the ground state wavefunction extent in units of $\hbar\,\nu$. The energy shift of the homogeneous part is given by $\hbar\,\omega_0=-g\,\mu_B\,B_0$, wherein $g$ denotes the g-factor and $\mu_B$ the Bohr magneton.
Thus the Hamiltonian of the system (cp. to fig.\,\ref{fsetup}) is given by (setting $\hbar=1$)
\begin{equation}\label{ham1}\begin{split}
H=&\nu\,b^\dagger\,b+\omega_e\,\ket{e}\bra{e}+\omega_0\,\left(\ket{+1}\bra{+1}-\ket{-1}\bra{-1}  \right)
		+\eta_{eff}\,\nu\,\left(\ket{+1}\bra{+1} -\ket{-1}\bra{-1} \right)\,\left(b+b^\dagger  \right)\\
		&+\Omega\,\left(\ket{e}\bra{+1}+\ket{+1}\bra{e}  \right)\cdot\left[\cos\left( k\,x-\omega_L\,t \right)+\cos\left(-k\,x-\omega_L\,t+\pi/2  \right)\right]\\
		&+\Omega\,\left(\ket{e}\bra{-1}+\ket{-1}\bra{e}  \right)\cdot\left[\cos\left(-k\,x-\omega^*_L\,t  \right)+\cos\left(k\,x-\omega^*_L\,t+\pi/2  \right)  \right]
\end{split}
\end{equation}
with $\omega_e$ being the energy of the upper level $\ket{e}$.
In order to end up with a time independent Hamiltonian we assume for the detuning: $\delta=\omega_e-\omega_0-\omega_L=\omega_e+\omega_0-\omega_L^*$, wherein $\omega_L$ is the frequency of the first laser pair (coupling $\ket{+1}$ to $\ket{e}$) and $\omega_L^*$ the frequency of the second pair respectively (coupling $\ket{-1}$ to $\ket{e}$).
In an interaction picture, expanding the Hamiltonian up to the first order in the Lamb-Dicke parameter $\eta$ and applying the rotating wave approximation, the Hamiltonian takes the form
\begin{equation}
\begin{split}
\label{ham2}
H=&\nu\,b^\dagger b+\delta\,\ket{e}\bra{e}+\eta_{eff}\,\nu \left(\ket{+1}\bra{+1}-\ket{-1}\bra{-1}\right)\,\left(b+b^\dagger \right)\\
  &+\left[\dfrac{\Omega}{2}\,(1+i)\left(\ket{e}\bra{+1}+\ket{e}\bra{-1}\right)\right.
  \left .+\eta\,\dfrac{\Omega}{2}\,(1+i)\,\left(b+b^\dagger  \right)\left(\ket{e}\bra{+1}-\ket{e}\bra{-1}\right)+\text{h.c.}\right].
\end{split}
\end{equation}
In order to cancel  the blue sideband heating contributions completely, it is important to adjust the magnetic gradient field in a way so that
\begin{equation}\label{ass1} \eta=\eta_{eff}     \end{equation}
as will be explained later. Values for the necessary magnetic gradient to fulfill condition (\ref{ass1}) for Ytterbium-172 and Calcium-40 ions are given in table\,\ref{tabelle1}. \\
\begin{table}
\caption{\label{tabelle1} Magnetic gradients needed to fulfill condition\,(\ref{ass1}) for different trapping frequencies. Values are given for the $S_{1/2}$ to $P_{1/2}$ transitions of Ytterbium-172 ions (369\,nm) and Calcium-40 ions (397\,nm). In this case the magnetic gradient is given by $\frac{\partial B}{\partial x}=\frac{2\,\nu\,\hbar}{(-g_e)\,\mu_B}\,k\cos\theta$, wherein $\theta$ denotes the angle between the laser beams and the trap axis ($\eta=\sqrt{\hbar\,k^2/(2\,m\,\nu)}\,\cos\theta$), which in the following is assumed to be zero and $g_e$ denotes the electron g-factor.}
\renewcommand{\arraystretch}{1.5}
\begin{tabular}{|l||c|c||c|c|} \hline
 &\multicolumn{2}{|c||}{$\boldsymbol{^{172}Yb^+}$}         &\multicolumn{2}{c|}{$\boldsymbol{^{40}Ca^+}$}\\\hline\hline
\textbf{Trapping frequency} $\boldsymbol{[\nu/(2\,\pi)]/Hz}$ & $\qquad 5\cdot10^5 \qquad $& $\qquad 1\cdot10^6 \qquad$ & $\qquad 5\cdot 10^5 \qquad$ & $\qquad 1\cdot 10^6 \qquad$\\\hline
\textbf{Magnetic gradient} $\boldsymbol{\frac{\partial B}{\partial x}}$\textbf{/(T/m)} & 607.6 & 1215 & 564.7 & 1130 \\\hline
\textbf{Lamb-Dicke parameter} $\boldsymbol{\eta}$  & 0.13 & 0.09 & 0.25 & 0.18  \\\hline
\end{tabular}
\end{table}
The basic idea of the cancellation mechanism is that $\sigma_x (b+b^+) + i \sigma_y (b-b^+)$ generates the red sideband and cancels the blue one.
Since the phonon coupling introduces a term that is proportional to $b+b^+$ and the magnetic gradient introduces a term that is proportional   to $i(b-b^+),$ we can expect that if the phases are adjusted appropriately  we can get a complete interference of the blue sideband.

The effect of the magnetic field is best seen by applying a unitary Schrieffer-Wolff transformation $U$\,\cite{schrieffer66} to the Hamiltonian (\ref{ham2})
\begin{equation}\label{sw1}  U=\mathrm{exp}\left\{-\eta\,\left(\ket{+1}\bra{+1}-\ket{-1}\bra{-1}\right)\,\left(b-b^\dagger\right)\right\}      \end{equation}
resulting in the following transformations
\begin{equation} \label{sw2}
\begin{split}
b \qquad &\rightarrow \quad b-\eta\,\left( \ket{+1}\bra{+1}-\ket{-1}\bra{-1}  \right)\\
\ket{+1}  \qquad &\rightarrow \quad \ket{+1}\,\mathrm{exp}\left[-\eta\left(b-b^\dagger\right)\right]    \\
\ket{-1} \qquad &\rightarrow \quad \ket{-1}\,\mathrm{exp}\left[\phantom{-}\eta\left(b-b^\dagger\right)\right]  \,\, .
\end{split}
\end{equation}

Application of the Schrieffer-Wolff transformation to the Hamiltonian (\ref{ham2}) and keeping only terms up to first order in $\eta$ we get (fig.\,\ref{fsetup}b, note that $\eta$ is equal to $\eta_{eff}$):
\begin{equation}\label{ham4}
H'=U\,H\,U^{\dagger}
  =\nu\,b^\dagger\,b+\delta\,\ket{e}\bra{e}
   +\left[\dfrac{\Omega}{\sqrt{2}}\,(1+i)\,\ket{e}\bra{B}\right.
   +\left.\eta\,\Omega\,\sqrt{2}\,(1+i)\,b\,\ket{e}\bra{D}+\text{h.c.}\right].
\end{equation}
Herein we have defined:
\begin{equation}
\begin{split}
\ket{B}&=\dfrac{1}{\sqrt{2}}\,\left(\ket{+1}+\ket{-1}\right),\\
\ket{D}&=\dfrac{1}{\sqrt{2}}\,\left(\ket{+1}-\ket{-1}\right).
\end{split}
\end{equation}

Note that the first order expansion terms of the Schrieffer-Wolff transformation ($\eta_{eff}$) of the zeroth order term in $\eta$ (Lamb-Dicke-Parameter) Hamiltonian contributions (see \ref{ham2}) lead to a cancellation of the first order in $\eta$ blue sideband laser interaction contributions under the assumption (\ref{ass1}).

The complete time evolution in the Schrieffer-Wolff transformed picture is given by
\begin{equation}\label{p_evol1} \dfrac{\partial}{\partial t}\rho=-i\left[ H',\rho \right]+\mathcal{L}^d \,\rho.  \end{equation}
The Liouvillian Superoperator $\mathcal{L}^d$ describes the spontaneous decay from the excited level $\ket{e}$ to the levels $\ket{+1}$ and $\ket{-1}$. It exhibits the following form:
\begin{equation} \mathcal{L}^d=\mathcal{L}_0^d+\mathcal{L}_1^d \,.    \end{equation}
$\mathcal{L}_0^d\,\rho$ describes the general spontaneous decay terms in zeroth order in the Lamb-Dicke Parameter:
\begin{equation}\label{spdec1} \mathcal{L}_0^d\,\rho=\sum_{i=\{(+1),(-1)\}}\dfrac{\gamma_i}{2}\,\Bigl[2\, \ket{i}\bra{e}\rho \ket{e}\bra{i}-\ket{e}\bra{e}\,\rho-\rho\,\ket{e}\bra{e}   \Bigr]\,.    \end{equation}
$\mathcal{L}_1^d\,\rho$ is the Schrieffer-Wolff transformation contribution of the spontaneous decay term (in first order in $\eta$):
\begin{equation}\label{spdec2} \mathcal{L}_1^d\,\rho = \sum_{i=\{(+1),(-1)\}} \gamma_i\,(-1)^i\,\eta\,\Bigl[\left(b-b^\dagger\right) \bm{,} \ket{i}\bra{e}\rho\ket{e}\bra{i} \Bigr]\,. \end{equation}
Note that recoil only contributes in order $\eta^4$, leading to a very small correction to the cooling rate (decreasing it) due to the dark state nature of the cooling scheme and therefore can be neglected for small $\eta$ (which was tested by numerical simulations).  This is due to the fact, that the probability for being in the excited state is given by $\eta^2$ and the spontaneous decay including recoil  from this state is a process in second order in $\eta$ and therefore recoil only contributes with a strength $\propto \eta^4 $.

If the final phonon population $\ket{n=0}$ is reached, the system will be decoupled from the laser interaction and will end up in the state $\ket{D,n=0}$ (`` dark state''). Note that due to the presence of only red sideband terms the phonon state $\ket{n=0}$ will always be reached (unless $\Omega=0$).  In the original (i.e. non Schrieffer-Wolff transformed) picture, this final state has the form:
\begin{equation}\label{state_origianl} U^\dagger\,\ket{D,n=0}=\ket{D,n=0}-\eta\,\ket{B,n=1}+\mathcal{O}\left(\eta^2\right)\,.    \end{equation}
Therefore the presented scheme cools to the vibrational ground state (n=0) in zeroth order in $\eta$. Furthermore the Schrieffer-Wolff transformation preserves pure states and thus does not increase the entropy (which remains zero).  A simple (projective) measurement of the internal state $\ket{D}$ as well as a unitary rotation could further improve the achieved final population, that is they would result in a state with $n=0$ up to second order in $\eta$.\par

\section{Conditions for optimal cooling rates and limitations}\label{three}
At first sight it might seem possible to increase the cooling rate by an arbitrary amount just by increasing the laser Rabi frequency $\Omega$. However a detailed analysis of the $\Omega\,\ket{e}\bra{B}$-coupling terms (cp. Hamiltonian (\ref{ham4})) shows that this is not the case. These terms are changing the $\ket{e}$-state  (shifting it and mixing it with the $\ket{B}$-state) and therefore crucially influence the cooling term  ($\sim \eta\Omega b \ket{e}\bra{D}$) . Optimal rates are achieved for tuning the cooling interaction in resonance with one of the dressed states of the $\Omega\,\ket{e}\bra{B}$ interaction. The Hamiltonian (\ref{ham4}) in the dressed state basis takes the form
\begin{equation}
\begin{split}
H'=&\nu\,b^\dagger\,b+\omega_u\,\ket{u}\bra{u}+\omega_d\,\ket{d}\bra{d}+\eta\,\Omega\,\sqrt{2}\,(1+i)\,b\,\left\{\left(a^u_e\ket{u}+a^d_e\,\ket{d}\right)\bra{D}+\text{h.c.}      \right\}
\end{split}
\end{equation}
wherein $\ket{u}$ and $\ket{d}$ is the dressed state basis (each having contributions of $\ket{e}$ and $\ket{B}$ states), $\omega_u$ and $\omega_d$ are the energies of those states (i.e. eigenvalues of the diagonalization procedure)
\begin{equation}\label{ds1}  \omega_{u/d}=\dfrac{1}{2}\,\left(\delta\pm\sqrt{\delta^2+4\,\Omega^2}\right)   \end{equation}
and $a_e^u$ and $a_e^d$ the amplitudes of $\ket{u}$ and $\ket{d}$ in $\ket{e}$ respectively (see also eq.\,(\ref{eamp1})), i.e. the excited level $\ket{e}$ expressed in the dressed state basis is given by
\begin{equation}\label{eamp1} \ket{e}=a^u_e\,\ket{u}+a^d_e\,\ket{d}\,.    \end{equation}
In the interaction picture with respect to the first three Hamiltonian contributions, the resonance condition (for tuning the $\ket{u}$ state to resonance) is easily seen to be fulfilled by setting $\omega_u=\nu $.
 This procedure leads to the following resonance condition for optimal cooling rates:
\begin{equation}\label{resonance_cond1}  \delta=\dfrac{\nu^2-\Omega^2}{\nu}\,.    \end{equation}
Consequently the detuning is - for a given trapping and laser Rabi frequency - uniquely determined in the optimal cooling case. \\
To see how the magnitude of the cooling term increases with increasing $\Omega$ it is useful to define an ``effective Rabi frequency`` $\Omega_{eff}$.
 If $\ket{u}$ corresponds to the resonant state (the one that has been tuned to resonance), then $\Omega_{eff}$ is defined as
\begin{equation}\label{oeffdef} \Omega_{eff}=\Omega\,a^u_e     \end{equation}
  and is a direct measure of the  strength of the resonant (and therefore relevant) cooling contribution. In case of fulfilling condition (\ref{resonance_cond1}) it follows that
\begin{equation}\label{oeff1}  |\Omega_{eff}|=\dfrac{\nu\,|\Omega|}{\sqrt{\nu^2+\Omega^2}}    \end{equation}
(see figure \ref{foeff}a).
\begin{figure}
\includegraphics[scale=0.16]{limit.png}
\caption{\label{foeff} a) $|\Omega_{eff}|/\nu$ is plotted over the Rabi frequency $\Omega$ in units of the trapping frequency $\nu$. It turns out that the maximal achievable $|\Omega_{eff}|$  is limited to the trapping frequency $\nu$ (see also eq.\,(\ref{oeff1})).
b) Absolute value of $\ket{e}$ amplitudes in the dressed states; the red curve represents the resonant state. The minimal amplitude determines the dynamics for a given $\Omega_{eff}$. }
\end{figure}
From equation (\ref{oeff1}) it is clearly seen, that it is impossible to increase the cooling term contribution above a limiting upper value, or more precisely, $\Omega_{eff}$ is limited to a maximal value, which corresponds to the trapping frequency $\nu$ ($|\Omega_{eff}|_{max}=\nu$). Therefore the trapping frequency is responsible for an upper limit to the maximal achievable cooling rate.

A second limitation arises from the decay rate $\gamma$ and the contributions of the excited state $\ket{e}$ in the new dressed state basis. Note that spontaneous decay to both dressed states and the $\ket{D}$ state takes place. The cooling process relies on the decay to the $\ket{D}$ state, and thus - to keep the process going - it is necessary that the other states decay fast enough, i.e. have a sufficient large $\ket{e}$ contribution. Otherwise the system will be trapped (at least for a certain amount of time) in one of the dressed states and thus the cooling rate is reduces significantly. For equal decay rates $\gamma$ in (\ref{spdec1}),(\ref{spdec2}), equal $\ket{e}$ contributions in the dressed states would lead to the optimal case when taking only the spontaneous decay influence into account. The  $\ket{e}$-amplitudes in the dressed states for the resonance case (\ref{resonance_cond1}) are plotted in figure \ref{foeff}b. Note that the minimal amplitude determines the dynamics. After a point of equal $\ket{e}$-amplitudes in the dressed states, the minimal $\ket{e}$-amplitude decreases again for further increasing $\Omega_{eff}$, which sets the optimal rate to a value $\Omega_{eff}<\nu$, depending on the relation between $\gamma$ and $\nu$.

\section{Rate equation formula and comparison to the numerical results}\label{sec_formula}\label{four}
The procedure described in \cite{cirac92,lindberg86} allows to calculate an analytic formula for the cooling rate by adiabatically eliminating the internal degrees of freedom under the assumption that the cooling dynamics is much slower than the internal dynamics. Moreover it is assumed that the trapping frequency is large compared to the cooling dynamics (i.e. $\eta\,\Omega \ll\nu$), which allows a projection on the phonon $\ket{n}\bra{n}$ subspace; i.e., a derivation of a rate equation. This second condition is only  fulfilled for small values of $\Omega$ in our case (otherwise the cooling term is of the order of the trapping frequency or for further increasing $\Omega$ even bigger).
The rate $W$ can be written in the form $W=A_--A_+$ and the final steady state population is given by $\langle n \rangle_{ss}=A_+/(A_--A_+)$, wherein $A_+$ /$ A_-$- are the transition rates for increasing or lowering the phonon number respectively. Those quantities consist of two time correlation functions of internal operators, which can be calculated using the quantum regression theorem and the optical Bloch equations.  It follows that, due to the fact that there exist no heating terms in first order in the Lamb-Dicke Parameter, $A_+$ is equal to zero and thus $W=A_-$ and the final phonon population is given by $\langle n  \rangle_{ss}=0$.\\
Following the procedure given in \cite{cirac92}, we get for the rate $W$ (see appendix for the calculation)
\begin{equation}\label{rate_formula}  W= \dfrac{8\,\eta^2\,\Omega^2\,\gamma\,\nu^2}{\gamma^2\,\nu^2+\left[(\delta_L-\nu)\,\nu+\Omega^2) \right]^2}  \,  ,  \end{equation}
where the rate $W$ is defined via:
\begin{equation}\label{rate_def}  \dfrac{\mathrm{d}}{\mathrm{d}t} \langle n\rangle =-W\, \langle n\rangle\, .   \end{equation}
This rate (\ref{rate_formula}) gets optimal if the term in brackets is equal to zero and the resulting condition is equivalent to the resonance condition given in eq.\,(\ref{resonance_cond1}). In that case, the rate takes the following form:
\begin{equation}\label{rate2}  W=\dfrac{8\,\eta^2\,\Omega^2}{\gamma} \,.   \end{equation}
The fact that the rate is proportional to $1/\gamma$ is a direct consequence of the adiabatic elimination in high $\gamma$. More precisely, the rate $W$ is proportional  to the spontaneous decay rate $\gamma$ of the excited level times the probability $p_u$ for achieving this level. The latter quantity is proportional to $1/\gamma^2$, i.e. the probability for reaching the excited dressed state (here: $\ket{u}$) decreases with increasing dissipation.  Therefore the rate is given by: $W\propto (1/\gamma^2)\,\gamma=1/\gamma$.
However it does not mean that the optimal achievable rate is getting slower for increasing $\gamma$, because it is possible to reach higher final values for $\Omega$ due to the fact, that the (internal) dynamics is less limited by the decreasing $\ket{e}$-amplitude of the resonant state, which occurs for increasing $\Omega$ beyond a point of equal dressed state amplitudes.
Figure\,\ref{frcomp} shows the rate obtained by numerical simulations and that one calculated using eq.\,(\ref{rate2}) for different Rabi frequencies. As expected, the rate formula describes the rate quite well in the limit of small $\Omega$ values whereas it fails completely in the opposite limit (when the cooling contribution is getting stronger and the conditions used for the formula derivation are no longer fulfilled).
\begin{figure}
\includegraphics[scale=0.9]{rcomp.png}
\caption{\label{frcomp} Cooling rate $W$ as a function of the Rabi frequency $\Omega$ both in units of the trapping frequency $\nu$ for $\eta=0.1$ and $\gamma=5\,\nu$ (left plot) / $\gamma=0.5\,\nu$ (right plot). The detuning $\delta$ for a given $\Omega$ is chosen to fulfill the resonance condition eq.\,(\ref{resonance_cond1}). The blue line is the rate obtained by numerical simulations, the red dashed line the one calculated using the rate equation formula (\ref{rate2}).}
\end{figure}

\section{Simulation and results}\label{five}
\begin{figure}
\includegraphics[scale=0.2]{sim16.png}
\caption{\label{fsim} Average population $\langle n \rangle$ vs. time in units of the trapping frequency $\nu$ for optimally chosen cooling parameters. The left plot describes the optimal case for $\gamma=0.5\,\nu$ ($\Omega=0.85\,\nu, \delta=0.28\,\nu$, $\eta=0.1$), the right plot for  $\gamma=5\,\nu$ ($\Omega=2.8\,\nu$, $\delta=-6.84\,\nu$, $\eta=0.1$). The insets show the final cooling stage of the corresponding figures using a logarithmic scale for the population, which clearly reveals the final population of the order of $\eta^2$. }
\end{figure}
Figure \ref{fsim} shows optimal cooling cycles for different spontaneous decay rates ($\gamma=0.5\,\nu<\nu$ and $\gamma=5\,\nu>\nu$). In the first case ($\gamma<\nu$, the rate is of the magnitude of $\gamma$) the optimal cooling rate is achieved for $\Omega_{eff}\approx\nu$ (cp. to eq.\,(\ref{resonance_cond1})) , i.e. the case of equal $\ket{e}$-amplitudes in the dressed states. That is what one would expect, because with further increasing $\Omega_{eff}$ to its maximal value $\nu$ there is one dressed state with a small $\ket{e}$ amplitude and therefore the slower spontaneous decay from this state leads to a decreasing rate despite the increasing strength of the cooling term contribution. In the $\gamma=5\,\nu, \gamma>\nu$, $W\ll\gamma$ (cooling rate) case, the spontaneous decay rate is much faster than the cooling dynamics. Therefore, as expected, it is possible to increase $\Omega_{eff}$ beyond the ''equal amplitude`` situation and thus a slightly higher effective Rabi frequency can be realized. The spontaneous decay time of the dressed states is much less limiting in this case. Note that this higher $\Omega$ is needed to compensate the influence of the quantum Zeno effect, which reduces the rate compared to the $\gamma=0.5\,\nu$ case for a given $\Omega$.
The optimal achieved rates $W$ (see also eq.\,(\ref{rate_def})) are $0.08\,\nu$ ($\Omega=0.85\,\nu, \Omega_{eff}=0.65\,\nu$) for the $\gamma=0.5\,\nu$ case and $0.09\,\nu$ ($\Omega=2.8\,\nu, \Omega_{eff}=0.94\,\nu$) for the $\gamma=5\,\nu$ case respectively.
Rates for different values of $\Omega$ are presented in figure \ref{frcomp}.

The dependence of the cooling rate $W$ on the spontaneous decay rate $\gamma$ is plotted in fig.\,\ref{fdiffgamma}. As discussed before, for higher spontaneous decay rates higher values of $\Omega$ are possible because of less restrictions due to the spontaneous decay and therefore the rate increases with increasing $\gamma$. For small values of $\gamma$, when the cooling dynamics is limited by the spontaneous decay rate and not by the $\Omega_{eff}$ coupling, the rate increases significantly with increasing $\gamma$. For high values of $\gamma$ the rate converges to a limiting value for the rate. In this case the cooling dynamics is clearly limited by the fact, that $\Omega_{eff}$ reaches its maximal value $\nu$ (the trapping frequency (\ref{oeff1}), see right inset of fig.\,\ref{fdiffgamma} for the corresponding $\Omega_{eff}$ values). Thus further increasing $\Omega$ in the region $\Omega_{eff} \lesssim\nu$ for increasing $\gamma$ (which is indeed necessary to compensate for the higher quantum Zeno effect) leads to a very small increasing of $\Omega_{eff}$ and therefore leads to a saturation of the rate with $\Omega_{eff}$ converging to its maximal value.
\begin{figure}
\includegraphics[scale=0.6]{diffgamma4.png}
\caption{\label{fdiffgamma} Cooling rate $W$ in units of the trapping frequency $\nu$ for different values of the spontaneous decay rate $\gamma$ (in units of the trapping frequency $\nu$). The inset shows the optimal value of $\Omega$ (left) and $\Omega_{eff}$ (right) for a given $\gamma$. While the rate increases fast in regions where the cooling is limited by the spontaneous decay rate (small values of $\gamma$) it saturates for high spontaneous decay rates, where the rate is limited by the limitation of $\Omega_{eff}$ (see also fig.\,\ref{foeff}). Values have been calculated for $\eta=0.1$, optimal values for the Rabi frequency and the detuning chosen according to the resonance condition eq.\,(\ref{resonance_cond1}).}
\end{figure}

The dependence on $\eta$ (for otherwise fixed parameters) is shown in fig.\,\ref{fdiffeta}. For $\eta\ll 1$, i.e. in the regime where heating terms of higher orders in $\eta$ can be neglected, the rate scales approximately linearly with  increasing Lamb-Dicke parameter $\eta$. Note that the rate equation (\ref{rate2}) includes a quadratic contribution of $\eta$, but also the optimal Rabi frequency depends on $\Omega$, scaling approximately $\propto 1/\sqrt{\eta}$ and thus leading to an approximately linear scaling\footnote{More exactly, numerical simulations show that the approximate scaling is given by $W\propto \eta^{1.28}$ and $\Omega\propto\eta^{-0.45}$.} of the rate dependence on $\eta$. Thus decreasing $\eta$ leads to increasing  $\Omega$ for the optimal condition, which can be explained as compensating for the quantum Zeno effect although it deteriorates the spontaneous decay amplitude behaviour. By increasing $\eta$ to higher values it reaches some optimal point and decreases again afterwards due to the fact that the heating terms in second order (and by further increasing also in higher orders) in $\eta$ become more and more important.
\begin{figure}
\includegraphics[scale=0.35]{diffeta4.png}
\caption{\label{fdiffeta} Cooling rate $W$ in units of the trapping frequency $\nu$ for different values of the Lamb-Dicke parameter $\eta$. The corresponding $\Omega$ dependence (for optimal rates) is shown in the inset. For small values of $\eta$ (when the influence of higher order heating terms is small) the rate scales approximately linear in $\eta$. All values have been calculated for a fixed spontaneous decay rate $\gamma=5\,\nu$ and optimal detunings according to the resonance condition\,(\ref{resonance_cond1}), including terms up to second order in $\eta$.}
\end{figure}

The cooling scheme concept is designed for the Lamb-Dicke regime ($\eta\ll1$), because only in this regime higher order contributions can be neglected. Numerical simulations for higher Lamb-Dicke parameters show that the (minimal) final population always increases with increasing Lamb-Dicke parameter. While the final population for $\eta=0.1$ is given by $\langle n\rangle = 0.01$ corresponding to $\eta^2$ as expected according to\,(\ref{state_origianl}), it is given by $\ex{n}\simeq0.06$ for $\eta=0.2$, $\ex{n}\simeq0.64$ for $\eta=0.3$ and $\ex{n}\simeq4.5$ for $\eta=0.6$ ($\Omega=13\,\nu$, $\gamma=5\,\nu$). Therefore only in the Lamb-Dicke regime the goal to achieve the vibrational ground state can be achieved. With respect to the new final population the cooling exhibits for higher Lamb-Dicke parameters in average again an exponential decay, leading e.g. for the above described $\eta=0.6$ case to a cooling rate of $W\simeq 0.006\,\nu$. Note however that the cooling rate and the final temperature are not optimized simultaneously for the same parameters, which especially for high Lamb-Dicke parameters leads to significant differences in the final populations and rates. The fact that both the final population and rate are worse for increasing Lamb-Dicke parameters underlines the restriction to the Lamb-Dicke regime.\vspace{1ex}

{Robustness to phase fluctuations --- } The presented scheme relies on destructive interference and thus may be fragile to dephasing, i.e., fluctuations in the phase difference between the laser pairs around the optimal $\pi/2$ value. The numerical analysis of the phase fluctuation is shown in fig. \ref{phase}.  It can be seen that even at fluctuations which are larger than $20\%$ the population is still less than $10^{-2}$ and the rates decrease only by a factor of $2.$
\begin{figure}
\includegraphics[scale=0.6]{phase.png}
\caption{\label{phase}
The left plot shows the final population $\langle n \rangle_{ss}$, the right plot the rate $W$ in units of $\nu$ for different phase differences $\varphi$ between the two coupling laser pairs. Note that $\varphi=\pi/2$ corresponds to the optimal case allowing the complete elimination of heating terms in first order in $\eta$. The simulation is performed in the Schrieffer-Wolff transformed picture and an Hamiltonian expansion up to second order in the Lamb-Dicke parameter with the parameters  $\Omega=2.8\,\nu, \delta=-6.84\,\nu$ (according to the resonance condition eq.\,(\ref{resonance_cond1})), $\eta=0.1$ and $\gamma=5\,\nu$ is used. Note that the final population corresponds to the population in the transformed picture.
}
\end{figure}

The proposed cooling scheme is faster  than previously proposed continuous cooling schemes. The fastest continuous scheme that was implemented is the EIT cooling scheme \cite{roos00} which cooled at a rate of $10^{-4} \nu$. The robust cooling scheme \cite{cerrillo10} can in principle reach a rate of $10^{-2} \nu.$ In contrary the proposed cooling scheme can reach rates which are one order of magnitude less than the trap frequency.

\section{Multimode cooling} \label{six}
The presented cooling scheme was tested for the case of several trapped ions. Furthermore the cooling scheme, i.e. the laser coupling of internal and external degrees of freedom as presented in section\,\ref{one}, is applied to one of the trapped ions (named ion\,1 in the following). It is advantageous to express the Hamiltonian in terms of normal modes, wherein the creator and annihilator of the n-th normal mode is denoted by $a_n^\dagger$ and $a_n$ respectively. In contrast, the ones of the i-th local mode are denoted by $b_i$ and $b_i^\dagger$. Expanding equation\,(\ref{ham2}) to the multi particle case leads to the Hamiltonian
\begin{equation}\label{sevm1}\begin{split} H=&\delta\,\ket{e}\bra{e}+\sum_{n} \nu_n\,a_n^\dagger\,a_n+ \eta\,\nu\,\left[\ket{+1}\bra{+1}-\ket{-1}\bra{-1} \right]\,\left(b_1+b_1^\dagger  \right) \\
					  &+\left[\dfrac{\Omega}{\sqrt{2}}\,(1+i)\,\ket{e}\bra{B}
					  +\eta\,\dfrac{\Omega}{\sqrt{2}}\,(1+i)\,\ket{e}\bra{D}\left(b_1+b_1^\dagger \right)+\text{h.c.}\right]+\mathcal{O}(\eta^2)\, .
 \end{split} \end{equation}
Herein $\nu_n$ describes the frequencies of the corresponding normal modes, the summation over $n$ equals a summation over all normal modes and all internal states correspond to those of ion\,1. With the transformation $x_i=\sum_{n}M_i^n\,Q_n$ between the local coordinates $\{x_i\}$ and the normal coordinates $\{Q_n\}$ the local operators in (\ref{sevm1}) can be expressed in terms of normal ones using the relation ($x_i=\sqrt{\hbar/(2m\nu)}\,(b_i+b_i^\dagger)$)
\begin{equation}\label{sevm3}  \left(b_i+b_i^\dagger \right)=\sum_{n}\,M_i^n\,\sqrt{\frac{\nu}{\nu_n}}\,\left(a_n+a_n^\dagger \right)     \end{equation}
and thus the Hamiltonian\,(\ref{sevm1}) reads
\begin{equation}\label{sevm4}\begin{split} H=&\delta\,\ket{e}\bra{e}+\sum_{n}\left[\nu_n\,a_n^\dagger\,a_n +\eta_n'\,\nu\,\left[\ket{+1}\bra{+1}-\ket{-1}\bra{-1}\right]\, \left(a_n+a_n^\dagger  \right) \right]\\ &+\left[\dfrac{\Omega}{\sqrt{2}}\,(1+i)\,\ket{e}\bra{B}
+\sum_{n}\eta_n' \,\dfrac{\Omega}{\sqrt{2}}\,(1+i)\,\ket{e}\bra{D}\,(a_n+a_n^\dagger)+\text{h.c.}\right]+\mathcal{O}(\eta^2)\, .
 \end{split}\end{equation}
wherein $\eta_n'=\eta\,M_1^n\,\sqrt{\nu/\nu_n}$.\\
Again, the effect of the magnetic gradient terms is best seen by moving to a Schrieffer-Wolff transformed picture in which the magnetic gradient displacements of the individual normal mode oscillators are removed. This transformation is given by
\begin{equation}\label{sevm8}\begin{split}  U'&=\text{exp}\left[-\sum_{n}\tilde{\eta}_n\,[\ket{+1}\bra{+1}-\ket{-1}\bra{-1}]\,(a_n-a_n^\dagger)\right] \\
		&=\prod_n\,\text{exp}\left[-\tilde{\eta}_n\,[\ket{+1}\bra{+1}-\ket{-1}\bra{-1}]\,\left(a_n-a_n^\dagger \right)\right]	 \end{split}\end{equation}
with $\tilde{\eta}_{n}=\eta_n'\,\nu/\nu_n=\eta\,M_1^n\,(\nu/\nu_n)^{3/2}$. Note that this transformation differs from the transformation used in sect.\,\ref{one} what can easily be checked  by inserting (\ref{sevm3}) into (\ref{sw1}).
Applying this transformation (cp. to eq.\,(\ref{sw2})) to the Hamiltonian\,(\ref{sevm1}) leads to $H'=U'\,H\,U'^\dagger$
\begin{equation}\label{sevm10}\begin{split} H'
  =&\delta \,\ket{e}\bra{e}+\sum_{n}\nu_n\,a_n^\dagger\,a_n+\biggl[\dfrac{\Omega}{\sqrt{2}}\,(1+i)\,\ket{e}\bra{B}\\
		      &+\eta\,\dfrac{\Omega}{\sqrt{2}}\,(1+i)\,\ket{e}\bra{D}\,\sum_{n}\,M_1^n\,\sqrt{\dfrac{\nu}{\nu_n}}\left(\left[1+\dfrac{\nu}{\nu_n}  \right]\,a_n+\left[1-\dfrac{\nu}{\nu_n} \right]\,a_n^\dagger  \right)+\text{h.c.}\biggr]+\mathcal{O}(\eta^2)\,.
 \end{split}\end{equation}
Consequently, the amplitudes of the blue and red sideband terms for a specific mode $n$ depend on the relation of the trapping frequency $\nu$ to the normal mode frequency $\nu_n$. A complete cancellation of blue sideband heating terms is only possible for the center of mass mode (COM), because only in this case $\nu_{COM}=\nu$. For higher modes heating terms are always present, although they are off-resonant by tuning the laser frequency to the red sideband. In general, for increasing $\nu_n$ the heating amplitude increases and the cooling amplitude decreases, which leads to equal amplitudes for $\nu_n\to\infty$. The latter case corresponds to the situation in cooling schemes that do not cancel blue sideband terms, e.g. EIT cooling and sideband cooling. Therefore the presented cooling scheme results in an improved cooling behaviour with respect to those schemes also in the multiparticle case albeit not being able to cool other modes than the COM-mode without heating terms. \par
\begin{figure}
\includegraphics[scale=0.2]{severalmodes.png}
\caption{\label{bsevmod3} Cooling of  three trapped ions. The line colours blue, red and green correspond to the single modes with increasing frequency in this order. The simulation was performed in the original (non transformed) picture using a second order in $\eta$ Hamiltonian expansion and the values $\Omega=2.8\,\nu, \gamma=5\,\nu$ and $\eta=0.1$. The insets show the corresponding cooling processes with a logarithmic scaling of the population. The left plot shows the cooling of the three modes all starting with an initial population $\langle n \rangle=2$. The detuning $\delta$ was chosen in a way to set the red sideband terms of the first and second mode to resonance respectively, corresponding to continuous and dashed lines. Continuous lines describe the situation when the laser is tuned to resonance with the red sideband of the center of mass mode and dashed lines the analogous situation with the breathing mode tuned to resonance (see resonance condition eq.\,(\ref{resonance_cond1})).  The right plot shows the heating of the excited modes during the cooling process of the center of mass mode, which is set to resonance. The final populations are: $\langle n_1 \rangle=0.005$, $\langle n_2 \rangle=0.031$ and $\langle n_3 \rangle=0.14$. }
\end{figure}
The cooling of the individual modes (COM mode (blue), first excited mode (breathing mode, red) and second excited mode (egyptian mode, green)) for the case of three ions is presented in fig.\,\ref{bsevmod3}. For three ions \cite{sasura01} $\vec{M}^1=(1/\sqrt{3},-1/\sqrt{2},1/\sqrt{6})$, $\nu_1=\nu$, $\nu_2=\sqrt{3}\,\nu$, $\nu_3=\sqrt{29/5}\,\nu$. The left plot shows the cooling starting with an initial population $\langle n \rangle=2$ for each mode.  Due to the fact that in the higher excited modes heating terms are always present while in the COM mode those contributions are eliminated completely in first order in $\eta$, the COM mode is always cooled with a faster rate also in the case when the cooling terms of any of the excited modes is set to resonance, i.e. the COM cooling transitions are off-resonant. The right plot shows the heating of the non-resonant excited modes originally in $\langle n \rangle =0$ during the resonant cooling process of the center of mass mode starting from $\langle n_{1} \rangle=2$. The modes are cooled/heated to final populations of $\langle n_1 \rangle=0.005$, $\langle n_2 \rangle=0.031$ and $\langle n_3 \rangle=0.14$. That is the heating turns out to be very small and increases with $\nu_n$ as expected from eq.\,(\ref{sevm10}).

\section{Conclusion and comments}
We presented a cooling scheme in the Lamb-Dicke regime, which completely eliminates the influence of the blue sideband heating transitions as well as the carrier transitions in lowest order in the Lamb-Dicke parameter,  resulting in a good cooling behaviour and providing cooling rates of the order of one magnitude below the trapping frequency. Albeit the Lamb-Dicke regime required for the Hamiltonian expansion and the resonance condition limit the cooling rate well below the trapping frequency, the achieved rates are faster than in comparable existing schemes working in the same regime (e.g. EIT cooling\,\cite{morigi00}, Robust cooling scheme\,\cite{cerrillo10}).  This is achieved by a proper choice of the magnetic gradient strength whose influence was analysed in a Schrieffer-Wolff transformed picture.  Optimal cooling conditions were derived together with a formula valid in the weak laser coupling case. The dependence on the Lamb-Dicke  parameter as well as on the spontaneous decay rate was presented. Moreover the cooling behaviour under deviations around the optimal phase differences between the coupling laser pairs was analysed showing that deviations of about 20 percent only reduce the rates by a factor of two and the final population is still of the order of $\eta^2$, both very satisfying values compared to other cooling schemes working in the Lamb-Dicke regime. Furthermore an analysis of the multiparticle case revealed good cooling behaviour also for several modes albeit not being able to eliminate the carrier and blue sideband terms completely up to first order in $\eta$ in other modes than the center of mass mode.  \par

\begin{acknowledgments}
Support by the European Commission (STREPs HIP and PICC), the Alexander von Humboldt foundation, the Deutsche Forschungsgemeinschaft, and secunet AG is acknowledged.
\vspace{5ex}
\end{acknowledgments}