\documentclass[reqno]{amsart}
\usepackage{amssymb,amsmath}
\usepackage{hyperref}

\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2008(2008), No. 79, pp. 1--15.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2008 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2008/79\hfil Three species Gross--Pitaevskii systems]
{Spatial patterns for the three species Gross--Pitaevskii system in the plane}

\author[M. Caliari, M. Squassina\hfil EJDE-2008/79\hfilneg]
{Marco Caliari, Marco Squassina}  % in alphabetical order

\address{Marco Caliari \newline
Dipartimento di Informatica \\
Universit\`a di Verona \newline
C\`a Vignal 2, Strada Le Grazie 15, 37134 Verona, Italy}
\email{marco.caliari@univr.it}

\address{Marco Squassina \newline
Dipartimento di Informatica \\
Universit\`a di Verona \newline
C\`a Vignal 2, Strada Le Grazie 15, 37134 Verona, Italy}
\email{marco.squassina@univr.it}


\thanks{Submitted March 14, 2008. Published May 28, 2008.}
\subjclass[2000]{35B40, 35Q55, 81V05, 81V45}

\keywords{Mixtures of Bose-Einstein condensates; spatial segregation;
\hfill\break\indent Gross-Pitaevskii vector equation; 
ground state solutions, Thomas-Fermi approximation}

\begin{abstract}
 In this paper we highlight some particular spatial patterns of ground 
 state solutions for the three species Gross--Pitaevskii system in the 
 plane having physical coefficients with particular attention to the 
 cases where the inter-species  coefficients become large. 
 The solutions models least energy stationary states of a mixture 
 of three Bose--Einstein condensates.
\end{abstract}


\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks


\section{Introduction}

Although Bose--Einstein condensates were predicted by Einstein \cite{einstein} 
around 1925, their successful experimental realization for atomic gases
was firstly achieved in 1995, see \cite{BecR}. Next, in 1997, the condensation for a mixture 
of two interacting species with the same mass
was realized, see \cite{BecSR}. Finally, around 2003, triplet species states were observed in \cite{triplet}.
In two recent papers \cite{CORT08,CalSqu} we investigated the numerical approximation (via spectral methods) and 
the large interaction patterns (via variational arguments) of ground state solutions for a class of vector Gross--Pitaevskii 
equations in $\mathbb{R}^2$ modelling a binary mixture of Bose--Einstein condensates \cite{Bigpap,PiTaBook}. As known, 
depending upon the anisotropy of the trapping potentials, there are various situations where the full 
physical model in $\mathbb{R}^3$ can be reduced, with a good approximation, to the planar case (see \cite[Section 2.2]{Baobab}),
which, therefore, is physically meaningful.
The aim of this paper is to complete the work initiated in \cite{CalSqu}, by showing some particular spatial 
patterns, in the large interaction regime, for the ground state solutions of
the three species 2D Gross--Pitaevskii system 
\begin{gather*}
\hbar i\partial_t\psi_1=-\frac{\hbar^2}{2m_1}\Delta \psi_1
+V_1(x_1,x_2)\psi_1+\theta_{11}\hbar^2|\psi_1|^2\psi_1
 +\sum_{j\neq 1}^3\theta_{1j}\hbar^2|\psi_j|^2\psi_1, \\
\hbar i\partial_t\psi_2=-\frac{\hbar^2}{2m_2}\Delta \psi_2
+V_2(x_1,x_2)\psi_2+\theta_{22}\hbar^2|\psi_2|^2\psi_2
 +\sum_{j\neq 2}^3\theta_{2j}\hbar^2|\psi_j|^2\psi_2, \\
\hbar i\partial_t\psi_3=-\frac{\hbar^2}{2m_3}\Delta \psi_3
+V_3(x_1,x_2)\psi_3+\theta_{33}\hbar^2|\psi_3|^2\psi_3
 +\sum_{j\neq 3}^3\theta_{3j}\hbar^2|\psi_j|^2\psi_3, 
\end{gather*}
for the unknown $\psi_i:\mathbb{R}^2\to\mathbb{C}$, where $\hbar$ denotes the reduced Planck constant 
and $m_i$ are the masses of the atomic species composing the Bose--Einstein triple mixture.
The coefficients of the coupling matrix $(\theta_{ij})$, which is symmetric so as to give
the system a variational structure, are positive and play the 
role of intra-species $(\theta_{ii})$ and inter-species $(\theta_{i\neq j})$ coefficients
respectively and can be represented as
$$
\theta_{ij}=2\pi\frac{\sigma_{ij}}{m_{ij}},\quad
\frac{1}{m_{ij}}=\frac{1}{m_i}+\frac{1}{m_j},\quad 
\sigma_{ij}=\sigma_{ji}, \quad i,j=1,2,3,
$$
where the constants $\sigma_{ij}$ are related to the scattering lengths for the $i$-$j$ species,
depending on the interaction potential between atoms. We point out that, due to Feshbach resonance,
the interspecies scattering lengths can be made large, by applying a suitable external magnetic
field (see \cite{feshbach}). Concerning the potential functions, we consider the general harmonic off-centered case, that is 
there are three centers $(x_1^i,x_2^i)$ and six positive constants $\omega_{ix},\omega_{iy}$, $i=1,2,3$, such that
$$
V_i(x_1,x_2)={\tfrac{m_i}{2}}\big(\omega_{ix}^2(x_1-x_1^i)^2
+\omega_{iy}^2(x_2-x_2^i)^2\big).
$$
The potential $V_i$s are often taken with the same centers, typically, without loss of generality, the origin.
On the other hand, there are some relevant physical situations which lead to consider the off-centered 
case (see e.g.\ \cite{ribolimod}). 


We will prove that, when one of the inter-species coefficients, say $\theta_{i_0j_0}$, becomes very large, 
then phase separation behaviour between the wave densities $\psi_{i_0}$ and $\psi_{j_0}$ tends to appear.
We shall highlight analytically (see Proposition \ref{Gsspatseg}) and numerically (see Figures \ref{as1}, \ref{as2}, \ref{as3}, \ref{as4}, \ref{as5}
within Section \ref{NumSim}) the spatial segregation of the solution components $\phi_i$
of the ground state. In general, this phenomenon can appear by two possibly coexisting causes, 
that is the separation of the trapping potential centers on one side (see Section \ref{TFapp} and Figure~\ref{as2})
and the large interaction regime on the other (see Section \ref{SIr} and Figure~\ref{as3}), the second one persisting also in absence
of external potentials.

\section{Functional setting}

Let $\mathcal{H}$ be the Hilbert subspace of $H^1(\mathbb{R}^2,\mathbb{C}^3)$ defined by
$$
\mathcal{H}=\Big\{(\psi_1,\psi_2,\psi_3)\in H^1(\mathbb{R}^2,\mathbb{C}^3):
\int_{\mathbb{R}^2} V_i(x_1,x_2)|\psi_i|^2<\infty,\; i=1,2,3\Big\},
$$
which is the natural framework for bound state solutions, endowed with the norm
$$
\|(\psi_1,\psi_2,\psi_3)\|_{\mathcal{H}}^2=\sum_{i=1}^3\int_{\mathbb{R}^2} \frac{\hbar^2}{2m_i}|\nabla 
\psi_i|^2+V_i(x_1,x_2)|\psi_i|^2,
$$
and consider the total energy associated with the system, given by the
Hamiltonian $E=E_\infty+J$, where $E_\infty,J:\mathcal{H}\to\mathbb{R}$ are defined by
\begin{gather*}
E_\infty(\psi_1,\psi_2,\psi_3) =\sum_{i=1}^3 E_\infty^i(\psi_i),   \\
J(\psi_1,\psi_2,\psi_3)  =\sum_{i\neq j}^3 J^{ij}(\psi_i,\psi_j),
\end{gather*}
being, for any $i,j=1,2,3$,
\begin{gather*}
E_\infty^i(\psi_i) =\int_{\mathbb{R}^2} \frac{\hbar^2}{2m_i}|\nabla\psi_i|^2+
V_i(x_1,x_2)|\psi_i|^2+\frac{\theta_{ii}\hbar^2}{2}|\psi_i|^4, \\
J^{ij}(\psi_i,\psi_j) =\theta_{ij}\hbar^2\int_{\mathbb{R}^2} |\psi_i|^2|\psi_j|^2.
\end{gather*}
We assume that $\theta_{ij}>0$ for any $i,j$. It is standard to see that, along a solution, the energy map
$$
\{t\mapsto E(\psi_1(\cdot,t),\psi_2(\cdot,t),\psi_3(\cdot,t))\},\quad t\geq 0
$$ 
is a constant and that the total particle numbers are time independent, 
\begin{equation}
\label{totmasses}
\int_{\mathbb{R}^2} |\psi_i(\cdot,t)|^2=N_i,\quad t\geq 0,\;  i=1,2,3.
\end{equation}
The \emph{ground state}  solution (also called \emph{least energy}  solution) of the Gross--Pitaevskii system 
is a solution $(\psi_1,\psi_2,\psi_3)\in \mathcal{H}$ with ansatz 
\begin{equation}
\label{ansatz}
\psi_i(x_1,x_2,t)=e^{-{\rm i}\frac{\mu_it}{\hbar}}\phi_i(x_1,x_2),\quad 
(x_1,x_2)\in\mathbb{R}^2,\; t\geq 0,\;  i=1,2,3
\end{equation}
where $(\phi_1,\phi_2,\phi_3)\in \mathcal{H}$ is real valued and minimizes the functional $E$
constrained to the normalization conditions \eqref{totmasses} (with $\phi_i$ in place of $\psi_i$).
Consequently, the functions $\phi_i$s solve the nonlinear eigenvalue problem
\begin{equation}
\label{GPSv}
\begin{gathered}
-\frac{\hbar^2}{2m_1}\Delta \phi_1
+V_1(x_1,x_2)\phi_1+\theta_{11}\hbar^2|\phi_1|^2\phi_1+\sum_{j\neq 1}^3
\theta_{1j}\hbar^2|\phi_j|^2\phi_1=\mu_1\phi_1, \\
-\frac{\hbar^2}{2m_2}\Delta \phi_2
+V_2(x_1,x_2)\phi_2+\theta_{22}\hbar^2|\phi_2|^2\phi_2+\sum_{j\neq 2}^3\theta_{2j}
\hbar^2|\phi_j|^2\phi_2=\mu_2\phi_2, \\
-\frac{\hbar^2}{2m_3}\Delta \phi_3
+V_3(x_1,x_2)\phi_3+\theta_{33}\hbar^2|\phi_3|^2\phi_3+\sum_{j\neq 3}^3\theta_{3j}
\hbar^2|\phi_j|^2\phi_3=\mu_3\phi_3, \\
{\int_{\mathbb{R}^2} \phi_i^2=N_i},\quad  i=1,2,3. 
\end{gathered}
\end{equation}
A direct computation yields the representation formula for the eigenvalues
\begin{equation}
\label{formulautov}
N_i\mu_i=E^i_\infty(\phi_i)+\frac{\theta_{ii}\hbar^2}{2}\int_{\mathbb{R}^2}|\phi_i|^4+\sum_{j\neq i}^3 J^{ij}(\phi_i,\phi_j),
\end{equation}
for any $i=1,2,3$. The existence of nontrivial solutions of the nonlinear system \eqref{GPSv} is straightforward
as all the coupling constants are positive, which
makes the Hamiltonian $E$ coercive and weakly lower semicontinuous on the $L^2\times L^2$ 
sphere \eqref{totmasses} in $\mathcal{H}$.  
In addition, by the standard gradient inequality $\int_{\mathbb{R}^2}|\nabla |\phi||^2\leq \int_{\mathbb{R}^2}|\nabla\phi |^2$, it follows that the ground 
state solutions can be sought among nonnegative functions, so the $\phi_i$ are positive.

\section{Large inter-species parameters}
\label{SIr}

Let $\mathcal{H}\subset H^1(\mathbb{R}^2)\times H^1(\mathbb{R}^2)$ be the realization of the Hilbert subspace
given in the introduction. For any index pair $i\neq j$, we set 
\begin{gather*}
{\mathcal S}  =\Big\{(\phi_1,\phi_2,\phi_3)\in \mathcal{H}:\,\int_{\mathbb{R}^2}\phi_i^2=N_i,\;\text{$\forall i=1,2,3$}\Big\}, \\
{\mathcal S}_{ij}^\infty  =\Big\{(\phi_1,\phi_2,\phi_3)\in {\mathcal S}:\, \int_{\mathbb{R}^2}\phi_i^2\phi_j^2=0\Big\}, \\
{\mathcal S}^\infty  =\bigcap_{i\neq j,\; i,j=1}^3{\mathcal S}_{ij}^\infty.
\end{gather*}
Assume now that one of the inter-species parameter, say 
$\theta_{i_0j_0}$ with $i_0\neq j_0$, gets very large, say
$\theta_{i_0j_0}=\kappa\to\infty$ while the other remain bounded, say $\theta_{lm}\in(0,1]$ for any other $l,m=1,2,3$ with $l\neq m$.
The least energy level of the ground state solutions is then defined and denoted as follows
$$
c_{\kappa}^{i_0j_0}=\inf_{(\phi_1,\phi_2,\phi_3)\in {\mathcal S}} 
\big[E_\infty(\phi_1,\phi_2,\phi_3)+J_\kappa(\phi_1,\phi_2,\phi_3)\big],
$$
where the Hamiltonian is $E_\infty+J_\kappa=E_\kappa=E:\mathcal{H}\to\mathbb{R}$, with 
$$
J_\kappa(\phi_1,\phi_2,\phi_3)=\kappa\hbar^2\int_{\mathbb{R}^2} |\phi_{i_0}|^2|\phi_{j_0}|^2
+\sum_{n\neq m,\,n\neq i_0,\,m\neq j_0}^3 J^{nm}(\phi_n,\phi_m).
$$
We also define the candidate for the limiting (as $\kappa\to\infty$) energy $c_\infty^{i_0j_0}$,
\begin{equation}
\label{specialenerlevel}
c_\infty^{i_0j_0}=\!\!\!\inf_{(\phi_1,\phi_2,\phi_3)\in {\mathcal S}_{i_0j_0}^\infty} \Big[E_\infty(\phi_1,\phi_2,\phi_3)
+\sum_{n\neq m,\,n\neq i_0,\,m\neq j_0}^3 J^{nm}(\phi_n,\phi_m)\Big] .
\end{equation}
With obvious modifications one can define the energy levels corresponding to the case
where more than one parameter diverges.
In the case where $\theta_{ij}\to\infty$ for \emph{all} $i\neq j$ then the limiting energy is $c_\infty$,
\begin{equation}
\label{specialenerlevelTT}
c_\infty=\!\!\!\inf_{(\phi_1,\phi_2,\phi_3)\in{\mathcal S}^\infty}E_\infty(\phi_1,\phi_2,\phi_3).
\end{equation}
As ${\mathcal S}^\infty\subset {\mathcal S}_{i_0j_0}^\infty\subset {\mathcal S}$, taking
into account the definition of $c_{\kappa}^{i_0j_0}$, $c_\infty^{i_0j_0}$ and $c_\infty$ it holds
\begin{equation}
\label{controlenergg}
c_{\kappa}^{i_0j_0}\leq c_\infty^{i_0j_0}\leq c_\infty,
\end{equation}
for any $\kappa>0$. In this setting the following result holds.

\begin{proposition}
\label{Gsspatseg}
As $\kappa$ goes to infinity, the sequence of ground state solutions $(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)\subset {\mathcal S}$ 
converges in $\mathcal{H}$ to a function $(\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)\in {\mathcal S}_{i_0j_0}^\infty$ 
at energy level $c_\infty^{i_0j_0}$. Moreover, there exist $\mu_i^\infty>0$ such that the variational inequalities hold
\begin{equation}
\label{varineqqq}
-\frac{\hbar^2}{2m_i}\Delta \phi_i^\infty+V_i(x_1,x_2)\phi_i^\infty
+\theta_{ii}\hbar^2|\phi_i^\infty|^2\phi_i^\infty
\leq \mu_i^\infty\phi_i^\infty \quad \text{in }\mathbb{R}^2,
\end{equation}
for all $i=1,2,3$.
\end{proposition}

\begin{remark} \label{rmk3.2} \rm
It is natural to wonder if the limit function $\phi_i^\infty$ solves the equation
\begin{equation*}
-\frac{\hbar^2}{2m_i}\Delta \phi_i^\infty+V_i(x_1,x_2)\phi_i^\infty+\theta_{ii}\hbar^2(\phi_i^\infty)^3
=\mu_i^\infty\phi_i^\infty \quad \text{in $\Omega_i=\{\phi_i^\infty>0\}$}
\end{equation*}
when $\Omega_i\subset\mathbb{R}^2$ is an open set.
In other words, taken any positive and compactly supported function $\varphi$ with support in $\Omega_i$,
do we have $\kappa\int_{\mathbb{R}^2} |\phi_{j}^\kappa|^2\phi_{i}^\kappa\varphi\to 0$ when $\kappa\to\infty$? 
We believe this is true.
\end{remark}

\begin{proof}[Proof of Proposition~\ref{Gsspatseg}]
In light of the first inequality of \eqref{controlenergg}, if 
$(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)\in \mathcal{H}$, $\phi_i^\kappa\not\equiv 0$ for any $i$ is the ground state solution,
we have $E_\kappa(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)=c_\kappa^{i_0j_0}$ and
\begin{align*}
\label{BddkappaInt}
\kappa\hbar^2\int_{\mathbb{R}^2} |\phi_{i_0}^\kappa|^2|\phi_{j_0}^\kappa|^2 
&\leq \sum_{n\neq m}^3\theta_{nm}\hbar^2\int_{\mathbb{R}^2} |\phi_n|^2|\phi_m|^2 \\
& \leq (J_\kappa+E_\infty)(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)\\
&=c_\kappa^{i_0j_0}\leq c_\infty^{i_0j_0},
\end{align*}
for every $\kappa>0$. As a consequence, we obtain $\int_{\mathbb{R}^2} |\phi_{i_0}^\kappa|^2|\phi_{j_0}^\kappa|^2\to 0$
as $\kappa\to\infty$. In addition, we have $\|(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)\|_\mathcal{H}^2\leq  
E_\kappa(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)\leq c_\infty^{i_0j_0}$
for any $\kappa>0$. Hence,
the sequences $(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)$ is bounded in $\mathcal{H}$, with respect to $\kappa$.
In particular there exist $(\phi_1^\infty,\phi_2^\infty,\phi_3^\kappa)$ in $\mathcal{H}$ 
such that, up to a subsequence, $(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)\rightharpoonup (\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)$ in $\mathcal{H}$ 
and $\phi_i^\kappa(x_1,x_2)\to\phi_i^\infty(x_1,x_2)$ for a.e.\ $(x_1,x_2)\in\mathbb{R}^2$ as $\kappa\to\infty$. Hence,
by Fatou's Lemma, we get $\int_{\mathbb{R}^2} (\phi_{i_0}^\infty)^2(\phi_{j_0}^\infty)^2=0$, namely
$\phi_{i_0}^\infty\phi_{j_0}^\infty=0$ a.e.\ in $\mathbb{R}^2$. 
Since by definition $\int_{\mathbb{R}^2}|\phi_i^\kappa|^2=N_i$ for any $\kappa>0$ and $\mathcal{H}$
in compactly embedded into the space $L^r(\mathbb{R}^2)\times L^r(\mathbb{R}^2)$ for any $r\geq 2$ (combining an inequality like \eqref{embedkeyineq} 
with the Gagliardo--Nirenberg inequalities), up to a further subsequence, we have
$\int_{\mathbb{R}^2}|\phi_i^\infty|^2=N_i$ for $i=1,2,3$. Whence
\begin{equation}
\label{appartenenzaclass}
(\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)\in {\mathcal S}_{i_0j_0}^\infty.
\end{equation}
Observe also that $E^i_\infty(\phi_i^\kappa)\leq c_\infty^{i_0j_0}$ for any $i=1,2,3$ and
\begin{equation*}
\sup_{\kappa\geq 1}\mu_i^\kappa=\frac{1}{N_i}\sup_{\kappa\geq 1}\Big\{E^i_\infty(\phi_i^\kappa)+\frac{\theta_{ii}\hbar^2}{2}
\int_{\mathbb{R}^2} |\phi_i^\kappa|^4+
\sum_{m\neq i}^3\theta_{im}\hbar^2\int_{\mathbb{R}^2} |\phi_m^\kappa|^2|\phi_i^\kappa|^2 
\Big\}<\infty,
\end{equation*}
denoting $\mu_i^\kappa$ the eigenvalues corresponding to  $\phi_i^\kappa$.
Hence, up to a subsequence, $\mu_i^\kappa\to \mu_i^\infty$ as $\kappa\to\infty$. By testing the equations of 
the system by an arbitrary compactly supported positive function $\eta$, we get
\begin{equation*}
\frac{\hbar^2}{2m_i}\int_{\mathbb{R}^2}\nabla \phi_i^\kappa\cdot\nabla\eta+\int_{\mathbb{R}^2} V_i(x_1,x_2)\phi_i^\kappa\eta+\theta_{ii}\int_{\mathbb{R}^2}|\phi_i^\kappa|^2
\phi_i^\kappa\eta \leq \mu_i^\kappa\int_{\mathbb{R}^2} \phi_i^\kappa\eta,
\end{equation*}
for all $\kappa>0$. Hence, letting $\kappa\to\infty$, it turns out that $\phi_i^\infty$
satisfies the variational inequalities \eqref{varineqqq}.
Notice that, by Fatou's lemma and the first inequality of \eqref{controlenergg}, 
we have
\begin{align*}
 & \sum_{i=1}^3\frac{\hbar^2}{2m_i}\int_{\mathbb{R}^2} |\nabla \phi_i^\infty|^2+ \sum_{i=1}^3\int_{\mathbb{R}^2} V_i|\phi_i^\infty|^2+ \sum_{i=1}^3\frac{\theta_{ii}\hbar^2}{2}  
 \int_{\mathbb{R}^2} |\phi_i^\infty|^4  \\
 &\quad+\lim_{\kappa\to\infty}\kappa\hbar^2\int_{\mathbb{R}^2} |\phi_{i_0}^\kappa|^2|\phi_{j_0}^\kappa|^2 
 +\!\!\!\sum_{n\neq m,\,n\neq {i_0},\,m\neq {j_0}}^3 
 J^{nm}(\phi_n^\infty,\phi_m^\infty)\\
& \leq  \sum_{i=1}^3\frac{\hbar^2}{2m_i}\liminf_{\kappa\to\infty}\int_{\mathbb{R}^2} |\nabla \phi_i^\kappa|^2+  
\sum_{i=1}^3\liminf_{\kappa\to\infty}\int_{\mathbb{R}^2} V_i|\phi_i^\kappa|^2+ 
 \sum_{i=1}^3 \frac{\theta_{ii}\hbar^2}{2}\liminf_{\kappa\to\infty}\int_{\mathbb{R}^2} 
 |\phi_i^\kappa|^4  \\
&\quad +\lim_{\kappa\to\infty}\kappa\hbar^2\int_{\mathbb{R}^2} |\phi_{i_0}^\kappa|^2|\phi_{j_0}^\kappa|^2 
  +\!\!\!\sum_{n\neq m,\,n\neq {i_0},\,m\neq {j_0}}^3 \liminf_{\kappa\to\infty}J^{nm}(\phi_n^\kappa,\phi_m^\kappa)\\
&  \leq \liminf_{\kappa\to\infty}
E_\kappa(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)
=\liminf_{\kappa\to\infty} c_\kappa^{i_0j_0}
  \leq c_\infty^{i_0j_0}.
\end{align*}  
Recalling formula \eqref{appartenenzaclass}, by the definition of $c_\infty^{i_0j_0}$ the above inequalities rewrite as
\begin{align*}
 & E_\infty(\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)
+\lim_{\kappa\to\infty}\kappa\hbar^2\int_{\mathbb{R}^2} |\phi_{i_0}^\kappa|^2|\phi_{j_0}^\kappa|^2 
 +\!\!\!\sum_{n\neq m,\,n\neq {i_0},\,m\neq {j_0}}^3 J^{nm}(\phi_n^\infty,\phi_m^\infty)\\
 & \leq c_\infty^{i_0j_0} \leq E_\infty(\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)  
+\!\!\!\sum_{n\neq m,\,n\neq {i_0},\,m\neq {j_0}}^3 J^{nm}(\phi_n^\infty,\phi_m^\infty)
\end{align*}
which yields
\begin{equation}
\label{strongerconcl}
\lim_{\kappa\to\infty}\kappa\int_{\mathbb{R}^2} |\phi_{i_0}^\kappa|^2 |\phi_{j_0}^\kappa|^2=0.
\end{equation}
Therefore, the convergence of $\phi_i^\kappa$ to $\phi_i^\infty$ in $\mathcal{H}$
is strong, otherwise, assuming by contradiction that this is not the case,
the previous inequalities would become strict, yielding immediately
a contradiction with \eqref{strongerconcl}. Finally, as a further
consequence, $c_\infty^{i_0j_0}=E_\infty(\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)$, 
concluding the proof. 
\end{proof}

\begin{remark}\rm
In the assumptions and notations of Proposition~\ref{Gsspatseg}, not all the components become mutually phase segregated.
If not, then, $\phi_i^\infty\phi_j^\infty=0$ a.e.\
for {\em any} index $i\neq j$, and so, by definition we would have $(\phi_1^\infty,\phi_2^\infty,
\phi_3^\infty)\in{\mathcal S}^\infty$. In turn, by \eqref{specialenerlevelTT}, \eqref{strongerconcl} and Fatou's Lemma, we obtain 
$$
c_\infty\leq E_\infty(\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)
\leq\liminf_{\kappa\to\infty}E_\kappa(\phi_1^\kappa,\phi_2^\kappa,\phi_3^\kappa)\leq c_\infty^{i_0j_0}<c_\infty,
$$
which is a contradiction. This is confirmed by numerical experiments.
\end{remark}

\begin{remark} \label{e3.3} \rm
The strong convergence of $\phi^\kappa_i$ to $\phi_i^\infty$ in $\mathcal{H}$, \eqref{formulautov} and \eqref{appartenenzaclass} yield
\begin{gather*}
N_i\mu_i =E^i_\infty(\phi_i^\infty)+\frac{\theta_{ii}\hbar^2}{2}\int_{\mathbb{R}^2}
|\phi_i^\infty|^4+\sum_{m\neq i}^3 J^{im}(\phi_i^\infty,\phi_m^\infty),\quad 
i\neq i_0, \\
N_{i_0}\mu_{i_0} =E^{i_0}_\infty(\phi_{i_0}^\infty)+\frac{\theta_{i_0i_0}\hbar^2}{2}
\int_{\mathbb{R}^2}|\phi_{i_0}^\infty|^4+\sum_{m\neq i_0,j_0}^3 J^{i_0m}(\phi_{i_0}^\infty,\phi_m^\infty).
\end{gather*}
\end{remark}


\begin{remark} \label{rmk3.4} \rm
Assume that one of the parameters $\omega_i$ in the trapping potentials gets 
very large, say
$\omega_{i_0x}=\Lambda\to+\infty$ while the other remain bounded, say 
$\omega_{ix},\omega_{iy}\in(0,1]$ 
for any $i=1,2,3$ with $i\neq i_0$. Then, numerical simulations show that 
the corresponding component
$\phi_{i_0}^\Lambda$ of the ground state tends to assume a cigar-like shape along one direction, and the bigger $\Lambda$ is,
the thinner is the profile of $\phi_{i_0}$ (see e.g.\ Figure~\ref{as1}). We show that, in the asymptotic process $\lambda\to+\infty$, 
contrary to what happens in the strong interaction limit $\kappa\to\infty$, the energies of the ground state solutions
cannot remain bounded from above. More precisely, set
$$
V_{i_0}^\Lambda(x_1,x_2)={\tfrac{m_{i_0}}{2}}\big(\Lambda^2(x_1-x_1^{i_0})^2+\omega_{i_0y}^2(x_2-x_2^{i_0})^2\big).
$$
Hence, we denote the least energy of the ground state solution as follows
$$
c_{\Lambda}=\inf_{(\phi_1,\phi_2,\phi_3)\in {\mathcal S}}\Big( E_\Lambda^{i_0}(\phi_{i_0})+\sum_{i\neq i_0}^3E_\infty^{i}(\phi_i)+J(\phi_1,\phi_2,\phi_3)\Big),
$$
where $E^{i_0}_\Lambda=E^{i_0}_\infty$ with $V_{i_0}=V^\Lambda_{i_0}$. Then we have
$$
\Gamma:=\sup_{\Lambda>0}c_{\Lambda}=+\infty.
$$
Indeed, assume by contradiction that this is not the case, namely $\Gamma<\infty$. Hence, the sequence of ground 
state solutions $(\phi_1^\Lambda,\phi_2^\Lambda,\phi_3^\Lambda)$ is bounded in $H^1(\mathbb{R}^2,\mathbb{R}^3)$. In particular,
up to a subsequence, it converges weakly in $H^1(\mathbb{R}^2,\mathbb{R}^3)$ and pointwise a.e.\ to a 
function $(\phi_1^\infty,\phi_2^\infty,\phi_3^\infty)$. Moreover, since it holds
\begin{equation}
\label{embedkeyineq}
\sup_{\Lambda\geq 1}\sup_{\varrho>0} \Big[\varrho^2\int_{\mathbb{R}^2\setminus B_\varrho(x_1^{i},x_2^{i})}\!\!\!|\phi_{i}^\Lambda|^2\Big]\leq\Gamma,
\end{equation}
for $i=1,2,3$, it follows that $(\phi_1^\Lambda,\phi_2^\Lambda,\phi_3^\Lambda)$ also converges, strongly, in $L^2(\mathbb{R}^2,\mathbb{R}^3)$.
Since, for any $\Lambda>0$ and $i=1,2,3$,
$$
\int_{\mathbb{R}^2} |\phi_{i}^\Lambda|^2=N_i,
$$
taking the limit as $\Lambda\to+\infty$ entails $\int_{\mathbb{R}^2} |\phi_{i}^\infty|^2=N_i$. Whence $\phi_{i}^\infty\neq 0$ in $H^1(\mathbb{R}^2)$
for every $i=1,2,3$. On the other hand, we also have
$$
\int_{\mathbb{R}^2} V_{i_0}^\Lambda(x_1,x_2)|\phi_{i_0}^\Lambda|^2\leq \Gamma,
$$
for all $\Lambda>0$, yielding in particular
$$
\int_{\mathbb{R}^2} \frac{m_{i_0}}{2} (x_1-x_1^{i_0})^2|\phi_{i_0}^\Lambda |^2\leq \frac{\Gamma}{\Lambda^2}.
$$
By Fatou's lemma this entails  $|x_1-x_1^{i_0}||\phi_{i_0}^\infty(x_1,x_2) |=0$ for a.e.\ $(x_1,x_2)\in\mathbb{R}^2$,
namely $\phi_{i_0}^\infty=0$ in $H^1(\mathbb{R}^2)$, which produces a contradiction.
\end{remark}
\medskip

\section{Location of components}
\label{TFapp}

In the so called Thomas--Fermi regime, a very good approximation of the ground state solutions of~\eqref{GPSv}
which holds for sufficiently large values of the coupling constants $\theta_{ij}$, can be obtained by simply
dropping the diffusion terms $-\Delta\phi_i$, the kinetic contributions, namely assuming the wave functions
to be slowly varying (cf.\ \cite{fermi,thomas,lieb}). In turn, system \eqref{GPSv} reduces to 
the algebraic system ($\hbar=m_i=1$)
\begin{equation}
\label{systemGPGenTF}
\begin{gathered}
\theta_{11}|\phi_1|^2+\theta_{12}|\phi_2|^2+\theta_{13}|\phi_3|^2=\mu_1-V_1(x_1,x_2), \\
\noalign{\vskip3pt}
\theta_{21}|\phi_1|^2+\theta_{22}|\phi_2|^2+\theta_{23}|\phi_3|^2=\mu_2-V_2(x_1,x_2), \\
\noalign{\vskip3pt}
\theta_{31}|\phi_1|^2+\theta_{32}|\phi_2|^2+\theta_{33}|\phi_3|^2=\mu_3-V_3(x_1,x_2).
\end{gathered}
\end{equation}
Let us denote by $\Theta=(\theta_{ij})$ the symmetric coupling matrix and set
$\eta_i=\phi^2_i$ and $\chi_i(x_1,x_2)=\mu_i-V_i(x_1,x_2)$,
where the eigenvalues $\mu_i$s should be computed through the normalization conditions \eqref{totmasses}.
Moreover, assume that $|\Theta|>0$ (positive determinant). Then, we obtain
\begin{gather}
\label{kramers}
|\Theta|\,\eta_1(x_1,x_2)=
\begin{vmatrix}
\chi_1(x_1,x_2) & \theta_{12} & \theta_{13} \\
\chi_2(x_1,x_2) & \theta_{22} & \theta_{23} \\
\chi_3(x_1,x_2) & \theta_{32} &\theta_{33}
\end{vmatrix},
\\
|\Theta|\,\eta_2(x_1,x_2)=
\begin{vmatrix}
\theta_{11} & \chi_1(x_1,x_2) & \theta_{13} \\
\theta_{21} & \chi_2(x_1,x_2) & \theta_{23} \\
\theta_{31} & \chi_3(x_1,x_2) &\theta_{33}
\end{vmatrix},\notag
\\
|\Theta|\,\eta_3(x_1,x_2)=
\begin{vmatrix}
\theta_{11} & \theta_{12} & \chi_1(x_1,x_2) \\
\theta_{21} & \theta_{22} & \chi_2(x_1,x_2) \\
\theta_{31} & \theta_{32} & \chi_3(x_1,x_2)  
\end{vmatrix}.\notag
\end{gather}
As the coupling coefficients are positive, if we set $r_i=\sqrt{2\mu_i}$ for $i=1,2,3$, it is evident that system
\eqref{systemGPGenTF} makes sense only if the right hand sides of each equation in it is positive, that is in the set
$$
{\mathcal D}=\cap_{i=1}^3{\mathcal D}_i,\quad
{\mathcal D}_i =\big\{(x_1,x_2)\in\mathbb{R}^2:\;\omega_{ix}^2(x_1-x_{i1})^2+\omega_{iy}^2(x_2-x_{i2})^2\leq r_i^2\big\}.
$$
Furthermore, taking into account that, for any $i$, the $\eta_i$s are positive and are a combination
of quadratic polynomials (due to the structure of $\chi_i$), there exist positive constants  
$\Omega_{ix}$, $\Omega_{iy}$ and $R_i$ and centers $(y_{i1},y_{i2})$ 
which allow to define the (possibly empty) overlap region of the components of the wave functions
$$
{\mathcal O}=\cap_{i=1}^3{\mathcal O}_i,\quad
{\mathcal O}_i =\big\{(x_1,x_2)\in{\mathcal D}:\;\Omega_{ix}^2(x_1-y_{i1})^2+\Omega_{iy}^2(x_2-y_{i2})^2\leq R_i^2\big\}.
$$
Then, for ${\mathcal O}\neq\emptyset$, there is $\alpha_i>0$ such that a non-smooth approximation of the 
$i$-th component of the ground state is given by
\begin{equation}
\label{TFrappppree}
\phi_i^2(x_1,x_2):=
\begin{cases}
\alpha_i(R_i^2-\Omega_{ix}^2(x_1-y_{i1})^2-\Omega_{iy}^2(x_2-y_{i2})^2),
 & \text{in ${\mathcal O}$}, \\[3pt]
\frac{r_i^2-\omega_{ix}^2(x_1-x_{i1})^2-\omega_{iy}^2(x_2-x_{i2})^2}
{2\theta_{ii}}, & \text{in ${\mathcal D}_i\setminus{\mathcal O}$}, \\[3pt]
0, & \text{in $\mathbb{R}^2\setminus {\mathcal D}_i$}.
\end{cases}
\end{equation}
Thinking for instance to the case where
the diagonal coefficients $\theta_{ii}$ are much larger than the $\theta_{ij}$s, i.e.\ $\theta_{ii}\gg\theta_{ij}\gg 1$, we have 
from \eqref{kramers}, e.g.\ for the component $\phi_1$,
$$
\theta_{11}\theta_{22}\theta_{33}\,\phi_1^2(x_1,x_2)\approx
\begin{vmatrix}
\chi_1(x_1,x_2) & \theta_{12} & \theta_{13} \\
\chi_2(x_1,x_2) & \theta_{22} & \theta_{23} \\
\chi_3(x_1,x_2) & \theta_{32} & \theta_{33}
\end{vmatrix}
\approx\chi_1(x_1,x_2)\theta_{22}\theta_{33}.
$$
This clarifies why it makes sense to extend to the set
${\mathcal D}_1\setminus{\mathcal O}$ the Thomas--Fermi approximation defined in ${\mathcal O}$
according to the second line of formula \eqref{TFrappppree} (for $i=1$).

\medskip

\section{Numerical computation of solutions}
\label{NumSim}

As done in \cite{CalSqu}, for the sake of self-containedness, we briefly describe the 
numerical algorithm used for the computation of the ground 
states. For further details, we refer the interested reader to \cite{CORT08}. It is sufficient
to consider the single one-dimensional Gross--Pitaevskii equation. In fact
the extension to the case with any number of equations is straightforward. 
Moreover, without loss of generality, we reduce to the case $\hbar=m=1$. 
The main idea is to directly minimize the energy $E(\phi)$ associated to a 
wave function $\psi(x)=e^{-i\mu t}\phi(x)$, discretized by 
Hermite functions, with a 
normalization constraint for the wave function. 
As it is known, the Hermite functions $(\mathcal{H}^\beta_l)_{l\in\mathbb{N}}$ are
defined by
$$
\mathcal{H}^\beta_l(x)=H^\beta_{l}(x)e^{-\frac12\beta^2x^2},\quad l\in\mathbb{N},
$$
where $(H_{l}^\beta)_{l\in\mathbb{N}}$ are the \emph{Hermite polynomials} 
\cite{boyd}, orthonormal in $L^2$ with respect to the weight $e^{-\beta^2x^2}$.
The Hermite functions are the solutions (ground state, for $l=0$, 
and excited states, if else) to the eigenvalue problem for the linear 
Schr\"odinger equation with \emph{standard} harmonic potential
$$
\frac{1}{2}\Big(-\frac{d^2}{dx^2}+(\beta^2x)^2\Big)\mathcal{H}_l=
\lambda_l\mathcal{H}_l,\quad
\lambda_l=\beta^2\big(l+\frac{1}{2}\big).
$$
If we set
$$
\phi=\sum_{l\in\mathbb{N}}\phi_l\mathcal{H}_l,
$$
where
$$
\phi_l=(\phi,\mathcal{H}_l)_{L^2}=\int_\mathbb{R} \phi\mathcal{H}_l,
$$
the energy functional rewrites as
\begin{equation*}
E(\phi)=\sum_{l\in\mathbb{N}}\lambda_l\phi_l^2+
\int_\mathbb{R} \Big(V(x)-\frac{(\beta^2x)^2}{2}\Big)
\Big(\sum_{l\in\mathbb{N}}\phi_l\mathcal{H}_l \Big)^2+
\frac12\theta\int_\mathbb{R}\Big(\sum_{l\in\mathbb{N}}\phi_l\mathcal{H}_l \Big)^4,
\end{equation*}
and the chemical potential turns into
\begin{equation}
\label{frapreMu}
N\mu=E(\phi)+
\frac12\theta\int_\mathbb{R}\Big(\sum_{l\in\mathbb{N}}\phi_l\mathcal{H}_l \Big)^4   
\end{equation}
By minimizing $E$, under the constraint $\|\phi\|_{L^2}^2=N$, 
we look for local minima of
\[
E(\phi;\lambda)=E(\phi)+\lambda\Big(N-\sum_{l\in\mathbb{N}}\phi_l^2\Big)
\]
which solve the system, with $k\in\mathbb{N}$, 
\begin{gather*}
(\lambda_\kappa-\lambda)\phi_\kappa+
\int_\mathbb{R} \Big(V(x)-\frac{(\beta^2x)^2}{2}\Big)\mathcal{H}_k
\Big(\sum_{l\in\mathbb{N}}\phi_l\mathcal{H}_l\Big)+
\theta\int_\mathbb{R} \mathcal{H}_k\Big(\sum_{l\in\mathbb{N}}\phi_l\mathcal{H}_l \Big)^3
=0,\\
\sum_{l\in\mathbb{N}}\phi_l^2=N.
\end{gather*}
We notice that, if $\phi$ is a solution of the above system, then it is 
immediately seen, by multiplying times $\phi_k$, summing up over $k$ and using \eqref{frapreMu},
that the Lagrange multiplier $\lambda$ equals the chemical potential $\mu$.
Next, we truncate the Hermite series 
to degree $L-1$ and introduce an additional parameter
$\rho=1$ in front of the first integral (its usage will be clear later), to
obtain a corresponding truncated energy functional 
$E_L(\phi;\lambda;\rho)$, whose local minima 
solve the system, with $0\le k\le L-1$,
\begin{gather*}
(\lambda_\kappa-\lambda)\phi_\kappa+
\rho\int_\mathbb{R} \Big(V(x)-\frac{(\beta^2x)^2}{2}\Big)\mathcal{H}_k
\Big(\sum_{l=0}^{L-1}\phi_l\mathcal{H}_l\Big)+
\theta\int_\mathbb{R} \mathcal{H}_k\Big(\sum_{l=0}^{L-1}\phi_l\mathcal{H}_l\Big)^3
=0,\\
\sum_{l=0}^{L-1}\phi_l^2=N.
\end{gather*}
To approximate the integrals, we used a Gauss--Hermite quadrature 
formula with $2L-1$ nodes relative to the weight $e^{-2\beta^2x^2}$.
The system is solved by a modified Newton method
with backtracking line-search, which guarantees global
convergence to the ground states. We refer to \cite{Baoal,CalThal} and, 
in particular, to \cite{CORT08} for the details. 
Here we just mention that 
the initial guess for the Newton iteration is obtained by a continuation 
technique over $\rho$ and $\theta$, starting from the ground state of
the Schr\"odinger equation with the standard harmonic potential,
which corresponds to $\rho=\theta=0$. 
Using the tensor basis of the Hermite functions, 
the extension to the two-dimensional case is straightforward. 

In the following figures we show some typical spatial patters of the 
ground states solution triplet with respect some relevant features as:
\begin{enumerate}
\item the anisotropy of the trapping potentials (Figure~\ref{as1});
\item the phase separation via potential off-centering (Figure~\ref{as2});
\item the phase separation via large inter-atomic interactions (Figure~\ref{as3}); 
\item the shape of supports with respect to the number of atoms $N_i$ (Figure~\ref{as4}); 
\item the shape of supports with respect to the size of the masses $m_i$ (Figure~\ref{as5}).
\end{enumerate}


\begin{figure}[!ht]
\includegraphics[scale=0.44]{figures/phi1AS1_1}\hfill
\includegraphics[scale=0.44]{figures/phi2AS1_1}\hfill
\includegraphics[scale=0.44]{figures/phi3AS1_1}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS1_2}\hfill
\includegraphics[scale=0.44]{figures/phi2AS1_2}\hfill
\includegraphics[scale=0.44]{figures/phi3AS1_2}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS1_3}\hfill
\includegraphics[scale=0.44]{figures/phi2AS1_3}\hfill
\includegraphics[scale=0.44]{figures/phi3AS1_3}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS1_4}\hfill
\includegraphics[scale=0.44]{figures/phi2AS1_4}\hfill
\includegraphics[scale=0.44]{figures/phi3AS1_4}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS1_5}\hfill
\includegraphics[scale=0.44]{figures/phi2AS1_5}\hfill
\includegraphics[scale=0.44]{figures/phi3AS1_5}
\caption{(anisotropy) ground state $(\phi_1,\phi_2,\phi_3)$ (left to right);
$\omega_{y1}$ and $\omega_{x2}$ assume values $\pi,1.1\pi,1.5\pi,2\pi,10\pi$ (top to bottom),
other $\omega_{yi}=\pi$, $\omega_{xi}=\pi$, $m_i = 1.44\cdot 10^{-25}$, $N_i=10^7$, $\sigma_{11}=\sigma_{22}=\sigma_{33}=10^{-6}$
and $\sigma_{12}=\sigma_{23}=\sigma_{13}=10\sigma_{ii}$.}
\label{as1}
\end{figure}
\clearpage

\begin{figure}[!ht]
\includegraphics[scale=0.42]{figures/phi1AS2_1}\hfill
\includegraphics[scale=0.42]{figures/phi2AS2_1}\hfill
\includegraphics[scale=0.42]{figures/phi3AS2_1}\vskip5pt
\includegraphics[scale=0.42]{figures/phi1AS2_2}\hfill
\includegraphics[scale=0.42]{figures/phi2AS2_2}\hfill
\includegraphics[scale=0.42]{figures/phi3AS2_2}\vskip5pt
\includegraphics[scale=0.42]{figures/phi1AS2_3}\hfill
\includegraphics[scale=0.42]{figures/phi2AS2_3}\hfill
\includegraphics[scale=0.42]{figures/phi3AS2_3}\vskip5pt
\includegraphics[scale=0.42]{figures/phi1AS2_4}\hfill
\includegraphics[scale=0.42]{figures/phi2AS2_4}\hfill
\includegraphics[scale=0.42]{figures/phi3AS2_4}\vskip5pt
\includegraphics[scale=0.42]{figures/phi1AS2_5}\hfill
\includegraphics[scale=0.42]{figures/phi2AS2_5}\hfill
\includegraphics[scale=0.42]{figures/phi3AS2_5}
\caption{(off-centering) ground state $(\phi_1,\phi_2,\phi_3)$ (left to right); $V_2$ with center
$(0,0)$; $V_1$ with centers $(-4,4)$, $(-3,3)$, $(-2,2)$,
$(-1,1)$, $(-0.4,0.4)$ and $V_3$ with centers $(4,4)$, $(3,3)$, $(2,2)$,
$(1,1)$, $(0.4,0.4)$ (up to $10^{-5}$, top to bottom);
$\omega_{yi}=\omega_{xi}=\pi$, $m_i = 1.44\cdot 10^{-25}$, $N_i=10^7$, $\sigma_{11}=\sigma_{33}=2\cdot10^{-7}$, $\sigma_{22}=100\sigma_{11}$
and $\sigma_{ij}=50\sigma_{11}$.}
\label{as2}
\end{figure}
\clearpage


\begin{figure}[!ht]
\includegraphics[scale=0.44]{figures/phi1AS3_1}\hfill
\includegraphics[scale=0.44]{figures/phi2AS3_1}\hfill
\includegraphics[scale=0.44]{figures/phi3AS3_1}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS3_2}\hfill
\includegraphics[scale=0.44]{figures/phi2AS3_2}\hfill
\includegraphics[scale=0.44]{figures/phi3AS3_2}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS3_3}\hfill
\includegraphics[scale=0.44]{figures/phi2AS3_3}\hfill
\includegraphics[scale=0.44]{figures/phi3AS3_3}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS3_4}\hfill
\includegraphics[scale=0.44]{figures/phi2AS3_4}\hfill
\includegraphics[scale=0.44]{figures/phi3AS3_4}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS3_5}\hfill
\includegraphics[scale=0.44]{figures/phi2AS3_5}\hfill
\includegraphics[scale=0.44]{figures/phi3AS3_5}
\caption{(phase segregation) ground state $(\phi_1,\phi_2,\phi_3)$ (left to right); 
$\omega_{yi}=\omega_{xi}=\pi$, $m_i = 1.44\cdot 10^{-25}$, $N_i=10^7$, $\sigma_{11}=\sigma_{22}=\sigma_{33}=10^{-6}$
and $\sigma_{12}=0,0.3,0.8,1.4,2\cdot 10^{-6}$, $\sigma_{23}=0,0.5,1,1.8,5\cdot 10^{-6}$,
$\sigma_{13}=0,0.7,1.8,5,50\cdot 10^{-6}$ (top to bottom).}
\label{as3}
\end{figure}
\clearpage

\begin{figure}[!ht]
\includegraphics[scale=0.44]{figures/phi1AS4_1}\hfill
\includegraphics[scale=0.44]{figures/phi2AS4_1}\hfill
\includegraphics[scale=0.44]{figures/phi3AS4_1}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS4_2}\hfill
\includegraphics[scale=0.44]{figures/phi2AS4_2}\hfill
\includegraphics[scale=0.44]{figures/phi3AS4_2}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS4_3}\hfill
\includegraphics[scale=0.44]{figures/phi2AS4_3}\hfill
\includegraphics[scale=0.44]{figures/phi3AS4_3}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS4_4}\hfill
\includegraphics[scale=0.44]{figures/phi2AS4_4}\hfill
\includegraphics[scale=0.44]{figures/phi3AS4_4}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS4_5}\hfill
\includegraphics[scale=0.44]{figures/phi2AS4_5}\hfill
\includegraphics[scale=0.44]{figures/phi3AS4_5}
\caption{(numbers of atoms) ground state $(\phi_1,\phi_2,\phi_3)$ (left to right); 
$\omega_{yi}=\omega_{xi}=\pi$, $m_i = 1.44\cdot 10^{-25}$, $\sigma_{ii}=10^{-6}$,
$\sigma_{ij}=4\sigma_{ii}$, $N_1=1,0.8,0.4,0.3,0.1\cdot 10 ^7$ and 
$N_3=1,1.3,2,4,6\cdot 10^7$ (top to bottom).}
\label{as4}
\end{figure}
\clearpage

\begin{figure}[!ht]
\includegraphics[scale=0.44]{figures/phi1AS5_1}\hfill
\includegraphics[scale=0.44]{figures/phi2AS5_1}\hfill
\includegraphics[scale=0.44]{figures/phi3AS5_1}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS5_2}\hfill
\includegraphics[scale=0.44]{figures/phi2AS5_2}\hfill
\includegraphics[scale=0.44]{figures/phi3AS5_2}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS5_3}\hfill
\includegraphics[scale=0.44]{figures/phi2AS5_3}\hfill
\includegraphics[scale=0.44]{figures/phi3AS5_3}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS5_4}\hfill
\includegraphics[scale=0.44]{figures/phi2AS5_4}\hfill
\includegraphics[scale=0.44]{figures/phi3AS5_4}\vskip5pt
\includegraphics[scale=0.44]{figures/phi1AS5_5}\hfill
\includegraphics[scale=0.44]{figures/phi2AS5_5}\hfill
\includegraphics[scale=0.44]{figures/phi3AS5_5}
\caption{(atomic masses)  ground state $(\phi_1,\phi_2,\phi_3)$ (left to right);
$\omega_{yi}=\omega_{xi}=\pi$, $N_i = 10^{7}$, $\sigma_{ii}=10^{-6}$,
$\sigma_{ij}=4\sigma_{ii}$, $m_1=1.44,1,0.8,0.5,0.3\cdot 10^{-25}$, $m_2=1.44\cdot 10^{-25}$ and 
$m_3=1.44,1.8,1.9,2,2.1\cdot 10^{-25}$ (top to bottom).}
\label{as5}
\end{figure}
\clearpage

\begin{thebibliography}{00}

\bibitem{BecR} M. H.\ Anderson, J. R.\ Ensher, M. R.\ Matthews,
 C. E.\ Wieman, E. A.\ Cornell;
Observation of Bose--Einstein condensation in a diluite atomic vapor, 
\emph{Science} {\bf 269} (1995), 198--201.

\bibitem{Baobab} W.\ Bao;
\textit{Ground states and dynamics of multicomponent Bose--Einstein condensates},
Multiscale Model.\ Simul.\ {\bf 2} (2004), 210--236.

\bibitem{Baoal} W.\ Bao, W.\ Tang;
Ground-state solution of Bose--Einstein condensate by directly minimizing the energy functional,
\emph{J.\ Comp.\ Phys.} {\bf 187} (2003), 230--254.  

\bibitem{boyd} J. P. Boyd; 
Chebyshev and Fourier spectral methods, Dover, New York, 2001.

\bibitem{CORT08} M.~Caliari, A.~Ostermann, S.~Rainer, M.~Thalhammer;
A minimisation approach for computing the ground state of
{G}ross--{P}itaevskii systems, (2008). Submitted.

\bibitem{CalSqu} M. Caliari, M. Squassina;
Location and phase segregation of ground
and excited states for 2D Gross--Pitaevskii systems, 
\emph{Dyn.\ Partial Diff.\ Equat.} (2008). To appear.

\bibitem{CalThal} M. Caliari, M. Thalhammer;
High-order time-splitting Fourier-Hermite spectral methods for the Gross--Pitaevskii equation,
preprint, (2007).

\bibitem{Bigpap} F.\ Dalfovo, S.\ Giorgini, L. P.\ Pitaevskii, S.\ Stringari;
Theory of trapped Bose--condensed gases,
\emph{Rev.\ Mod.\ Phys.} {\bf 71} (1999), 463--512.

\bibitem{einstein} A.\ Einstein;
 \emph{Sitz.\ Ber.\ Preuss.\ Akad. Wiss.} (Berlin) 1, 3 (1925).

\bibitem{fermi} E.\ Fermi; Un metodo statistico per la determinazione
di alcune propriet\`a dell'atomo, 
\emph{Rend.\ Lincei}  {\bf 6}, (1927), 602--607. 

\bibitem{gross} E. P.\ Gross;
Structure of a quantized vortex in boson systems, 
\emph{Nuovo Cimento} {\bf 20} (1961), 454--477.

\bibitem{feshbach} S.\ Gupta, Z.\ Hadzibabic, M. W.\ Zwierlein, C. A.\ Stan, 
K.\ Dieckmann, C. H.\ Schunck, E. G. M.\ van Kempen, B. J.\ Verhaar, 
W.\ Ketterle; Radio-frequency spectroscopy of ultracold fermions,
\emph{Science } {\bf 300}, 1723--1726.

\bibitem{lieb} E. H.\ Lieb;
Thomas--Fermi and related theories of atoms and molecules, 
\emph{Rev. Mod. Phys.} {\bf 53} (1981), 603--641.
 
\bibitem{BecSR} C. J.\ Myatt, E. A.\ Burt, R. W.\ Ghrist, E. A.\ Cornell, 
C. E.\ Wieman; 
Production of two overlapping Bose--Einstein condensates by sympathetic cooling,
\emph{Phys.\ Rev.\ Lett.} {\bf 78} (1997), 586--589.
 
\bibitem{PiTaBook} L. P.\ Pitaevskii,  S.\ Stringari,
Bose--Einstein condensation,
Oxford Science Publications,
Int.\ Series of Monographs on Physics, 2003.

\bibitem{ribolimod} F.\ Riboli, M.\ Modugno;
Topology of the ground state of two interacting Bose--Einstein condensates
\emph{Phys.\ Rev.\ A} {\bf 65} 063614 (2002).

\bibitem{triplet} Ch.\ R\"uegg,  N.\ Cavadini,  A.\ Furrer,  H-U.\ G\"udel,  
K. Kr\"amer,  H.\ Mutka,  A.\ Wildes,  K.\ Habicht,  P.\ Vorderwisch;
Bose--Einstein condensation of the triplet states in the magnetic 
insulator TlCuCl3,
\emph{Nature} {\bf 423} (6935):62-5, (2003).

\bibitem{thomas} L.H.\ Thomas; 
The calculation of atomic fields, 
\emph{Proc.\ Cambridge Phil.\ Soc.\ } {\bf 23} (1927), 542--548.

\bibitem{timmermans} E.\ Timmermans, 
Phase separation of Bose--Einstein condensates, 
\emph{Phys.\ Rev.\ Lett.} {\bf 81} (1998), 5718--5721.


\end{thebibliography}

\end{document}



