# Multi-component states for trapped spin-1 Bose-Einstein Condensates in the
presence of magnetic field

Projjwal K. Kanjilal<EMAIL_ADDRESS>Department of
Physics, Indian Institute of Science Education and Research, Dr. Homi Bhabha
Road, Pune 411 008, India  A. Bhattacharyay<EMAIL_ADDRESS>Department of Physics, Indian Institute of Science Education and Research, Dr.
Homi Bhabha Road, Pune 411 008, India

(August 30, 2024)

###### Abstract

In presence of a magnetic field, the multi-component ground states appear for
trapped spin-1 Bose-Einstein condensates for both ferromagnetic and anti-
ferromagnetic types of spin-spin interaction. We aim to produce an accurate
analytical description of the multi-components states which is of fundamental
importance. First, we show the analytical results provided by the Thomas-Fermi
approximation, which is one of the most trusted and widely used analytical
method. Despite being in the regime of Thomas-Fermi approximation, it provides
wrong physical conclusions. We generalize the variational method introduced in
the article Kanjilal, Projjwal Kanti and Bhattacharyay, A. (2022) and
demonstrate that the variational method is crucial in giving a proper
analytical description of the multi-component states which is also backed by
numerical simulation.

##  I Introduction

The successful experimental realization of Bose-Einstein condensates (BEC)
with alkali atoms Davis  _et al._ (1995); Anderson  _et al._ (1995); Bradley
_et al._ (1995) inside a magnetic trap spurred renewed interest Giorgini  _et
al._ (1996); Dalfovo and Stringari (1996); Stringari (1996) in ultracold
atomic physics. Soon it attracted a lot of attention from both atomic- and
condensed-matter physics communities as it provided an ideal test bed as a
quantum simulator Dupont  _et al._ (2021); Senaratne  _et al._ (2018);
Yamamoto and Takahashi (2016); Bravo  _et al._ (2015); Cirac and Zoller
(2012); Georgescu  _et al._ (2014) and precision measurements Becker  _et al._
(2018); Szigeti  _et al._ (2020); Gupta  _et al._ (2002); Jamison  _et al._
(2011); Hanneke  _et al._ (2008); Dutta  _et al._ (2022) for its unprecedented
experimental control.

The early experiments on BEC Ensher  _et al._ (1996a, b); Ketterle  _et al._
(1996) were done in the magnetic trap which only captures atoms with weak-
field-seeking hyperfine states, thus, the magnetic degrees of freedom were
frozen in the resulting BEC. Later, with the optical trapping technique this
limitation was overcome, and spinor BECs were created with all the hyperfine
states of the constituent spin-f atoms (f is an integer) Stamper-Kurn  _et
al._ (1998); Chang  _et al._ (2004); Stamper-Kurn and Ketterle . The order
parameter of such a system constitutes $(2f+1)$ components. Due to the
interplay of magnetic field and interatomic interaction, the spinor BEC shows
a rich variety of phenomena including spin-textures Kawaguchi  _et al._
(2006); Leanhardt  _et al._ (2003), domain structures Timmermans (1998);
Świsłocki and Matuszewski (2012); Xi  _et al._ (2011); Bandyopadhyay  _et al._
(2017); Isoshima  _et al._ (1999); Matuszewski  _et al._ (2008, 2009); Ho and
Shenoy (1996); Bhuvaneswari  _et al._ (2018); Wen  _et al._ (2012); Liu
(2009); Zhu and Li (2017); Tojo  _et al._ (2010); Li  _et al._ (2018); Lee
_et al._ (2016); Sabbatini  _et al._ (2011); Gautam and Adhikari (2014);
Gautam and Angom (2011); Saito and Ueda (2005); Kanjilal and Bhattacharyay
(2020); *Kanjilal_corri, and topological phases Ueda (2014). Spinor BEC also
attracted a lot of attention due to its complex soliton structures Nistazakis
_et al._ (2008); Bersano  _et al._ (2018); Katsimiga  _et al._ (2021); Hegde
_et al._ (2022), interesting few-body physics in low dimensions Mistakidis
_et al._ (2022).

The role of the multi-component ground state where more than one spin-
projections is populated is crucial for the plethora of interesting phenomena
that occurs in spinor BEC. As a result, there has been a lot of seminal
studies on multi-component ground states Ho (1998); Yi  _et al._ (2002); Zhang
_et al._ (2003); Gautam and Adhikari (2015). But most of the analytical
studies focus on the single-mode approximation which is only followed in
absence of a magnetic field. In this article, we look at the ground state
structures involving the multi-modal multi-component stationary states of
confined spin-1 BEC in the presence of the magnetic field without any finite
temperature Phuc  _et al._ (2011); Kawaguchi  _et al._ (2012); Serrano-
EnsÃ¡stiga and Mireles (2022) or inter-particle correlation effects Mittal
_et al._ (2020). The contribution of the magnetic field comes into play in
terms of linear and quadratic Zeeman terms ($p$, $q$ respectively) that are
present in the Gross-Pitaevskii (GP) equation which governs the dynamics of
the three component order parameter (the mean fields $\psi_{1}$, $\psi_{0}$
and $\psi_{-1}$). In absence of trapping, solving the GP equation one can get
to the phase diagram in $p$, $q$ parameter space Stenger  _et al._ (1998);
Kawaguchi and Ueda (2012); Stamper-Kurn and Ketterle . In a realistic
situation, the potential energy due to trapping and interaction along with the
kinetic energy play important role in determining the ground state properties.
To get an analytical description for such a situation one takes resorts to the
Thomas-Fermi (T-F) approximation which is a well-known method for its
simplicity and accuracy at large particle number limit. It is well-accepted
that the T-F approximation produces ground state structures with reasonable
accuracy if the density is high enough to neglect the kinetic energy term in
comparison to interaction terms. However, there exist reasons for the
breakdown of T-F approximation for competing ground state candidates of
comparable free energy even for large condensates as well as small ones.

In this article, we present such a situation where the kinetic term is
actually negligible in comparison to the interaction terms and hence in the
range of applicability of the T-F approximation. But, the T-F provides the
wrong physical conclusion. Such a situation arises in a multi-modal multi-
component ground state where the density of the sub-components varies quite a
lot. T-F approximation in such cases wrongly predicts a domain formation of
multi-component states along with a peripheral single-component state.
However, a detailed investigation beyond T-F approximation would reveal the
actual scenario which is shown in this paper. We have extended and generalized
the formalism of the multi-modal variational method, introduced in the article
Kanjilal, Projjwal Kanti and Bhattacharyay, A. (2022) for a simplistic
situation in absence of the magnetic field. We will show that this generalized
formalism provides quite an accurate description in comparison with the
numerically obtained number density profiles.

One of the key findings of this article is that the T-F approximation requires
significant correction for multi-component stationary states that also become
the ground state for spin-1 BEC under trapping. It also comes out from the
present analysis that one cannot avoid a multi-modal analysis while dealing
with such multi-component ground states to correct for limitations of T-F
approximation. We consider quasi-one-dimensional ${}^{87}Rb$ and ${}^{23}Na$
systems under harmonic trapping. For the ${}^{87}Rb$, T-F approximation
predicts a domain-like structure between the phase-matched state near the
center of the trap followed by the polar state outside. Similarly, for a
specific choice of Zeeman terms, the T-F approximation predicts the ground
state to be a domain structure between the anti-ferromagnetic and
ferromagnetic states. Using the generalized multi-modal variational method, we
show that there is no domain structure in the ground states of the above-
mentioned cases. The multi-component stationary states at the core of the trap
are actually the ground states. The results of the variational method are
validated through a comparison with the numerical simulation.

The article is organized as follows: In section II, we discuss the mean-field
theory of the spin-1 trapped BEC in presence of the magnetic field. In section
III we focus on the analytical description of the phase-matched state and the
anti-ferromagnetic state that becomes the ground state for ${}^{87}Rb$ and
${}^{23}Na$. respectively. We start with the T-F results and then we provide a
description based on the variational method compare it with the numerical
simulation. In this section, we also discuss why T-F approximation cannot be
trusted and the variational method must be used for multi-component states.
This follows with a short discussion about the generality of the treatment
presented here and possible directions to explore.

##  II Mean field theory: GP equation

The dynamics of the mean fields for a trapped spin-1 BEC in presence of the
magnetic field is captured in the Gross-Pitaevskii (GP) equation Ho (1998);
Ohmi and Machida (1998); Kanjilal and Bhattacharyay (2020); Kawaguchi and Ueda
(2012),

$\begin{split}i\hbar\dfrac{\partial\psi_{m}}{\partial
t}=\Bigg{(}-\dfrac{\hbar^{2}\nabla^{2}}{2M}+U(\vec{r})&-pm+qm^{2}+c_{0}n\Bigg{)}\psi_{m}\\\
&+c_{1}\sum_{m^{\prime}=-1}^{1}\vec{F}.\vec{f}_{mm^{\prime}}\psi_{m^{\prime}},\end{split}$
(1)

where the first term on the r.h.s of Eq.1 is the kinetic energy contribution
for particles of mass M. The second term is due to the confining potential and
the presence of the magnetic field is captured by the linear and quadratic
Zeeman terms $p$ and $q$ respectively. We assume a two-body contact
interaction that can be decomposed into spin-independent (the term involving
$c_{0}$) and inter-spin interactions (the term with $c_{1}$). As this is a
spin-1 system, the suffix $m$ and $m^{\prime}$ take the values $1,\ 0$ and,
$-1$.

The total density n is defined as,

$n(r)=\sum^{1}_{m=-1}|\psi_{m}|^{2},$ (2)

and the local spin density $\vec{F}$ is

$\vec{F}=\sum^{1}_{i,j=-1}\psi^{*}_{i}\vec{f}_{ij}\psi_{j},$ (3)

where the $\vec{f}$ is defined via the spin-1 Pauli matrices Kawaguchi and
Ueda (2012). So, Eq.1 is the GP equation which is a set of three coupled non-
linear partial differential equations which yields the order parameter
$(\psi_{1},\psi_{0},\psi_{-1})$ as its solution.

The mean-field approximated total energy, Kawaguchi and Ueda (2012)

$\begin{split}E=\int
d\vec{r}\sum^{1}_{m=-1}\psi_{m}^{*}&\left(-\dfrac{\hbar^{2}\nabla^{2}}{2M}+U_{trap}(\vec{r})-pm+qm^{2}\right)\psi_{m}\\\
&\qquad\qquad\qquad+\dfrac{c_{0}}{2}n^{2}+\dfrac{c_{1}}{2}|\vec{F}|^{2},\end{split}$
(4)

can be compared for different stationary states to get to the ground state of
the system. In this paper, we represent the stationary states by notation
$(\mathbbm{n}_{1},\mathbbm{n}_{0},\mathbbm{n}_{-1})$, where, $\mathbbm{n}_{m}$
is the placeholder for the binary notation, 0 or 1. If the sub-components are
populated we represent it as 1 and if unpopulated we represent it as 0. In
this notation, for example, the ferromagnetic state is represented as
$(1,0,0)/(0,0,1)$ where the sub-component corresponding to $m=1$ or $m=-1$ is
populated.

To simplify the GP equations further, one can write the mean fields in terms
of the density and corresponding phase,

$\psi_{m}(\vec{r},t)=\sqrt{n_{m}(\vec{r})}exp(-\dfrac{i\mu
t}{\hbar})exp(-i\theta_{m}),$ (5)

where the parameter $\mu$ stands for chemical potential. One can get the
number and phase dynamics separately Kanjilal and Bhattacharyay (2020) by
using the ansatz Eq.5 in the GP equation Eq.1,

$\dot{n}_{0}(\vec{r})=-\dfrac{4c_{1}n_{0}\sqrt{n_{1}n_{-1}}\sin\theta_{r}}{\hbar},$
(6) $\dot{n}_{\pm
1}(\vec{r})=\dfrac{2c_{1}n_{0}\sqrt{n_{1}n_{-1}}\sin\theta_{r}}{\hbar},$ (7)
$\begin{split}\hbar\dot{\theta}_{0}=\dfrac{1}{\sqrt{n_{0}(\vec{r})}}&\left(\mathcal{H}-\mu\right)\sqrt{n_{0}(\vec{r})}\\\
&+c_{1}\left(n_{1}+n_{-1}+2\sqrt{n_{-1}n_{1}}\cos\theta_{r}\right),\end{split}$
(8) $\begin{split}\hbar\dot{\theta}_{\pm 1}=\dfrac{1}{\sqrt{n_{\pm
1}(\vec{r})}}&\left(\mathcal{H}-\mu\right)\sqrt{n_{\pm 1}(\vec{r})}\pm
c_{1}\left(n_{1}-n_{-1}\right)+q\\\ &\quad\mp
p+c_{1}n_{0}\left(1+\sqrt{\dfrac{n_{\mp 1}(\vec{r})}{n_{\pm
1}(\vec{r})}}\cos\theta_{r}\right),\end{split}$ (9)

where, $\mathcal{H}=-\dfrac{\hbar^{2}\nabla^{2}}{2M}+U(\vec{r})+c_{0}n$. The
relative phase $\theta_{r}$ in Eq.8-9 is defined as,
$\theta_{r}=\theta_{1}+\theta_{-1}-2\theta_{0}$ Kawaguchi and Ueda (2012). The
same ansatz Eq.5 makes the energy (Eq.4) a function of the sub-component
number density and the relative phase,

$\begin{split}E&=\int d\vec{r}e(\vec{r})\\\ &=\int
d\vec{r}\Bigg{(}-\sum^{1}_{m=-1}\sqrt{n_{m}(\vec{r})}\dfrac{\hbar^{2}\nabla^{2}}{2M}\sqrt{n_{m}(\vec{r})}\\\
&\qquad+U(\vec{r})n(\vec{r})-p(n_{1}-n_{-1})+q(n_{1}+n_{-1})\\\
&\qquad+\dfrac{c_{0}}{2}n^{2}(\vec{r})+\dfrac{c_{1}}{2}\left(n_{1}-n_{-1}\right)^{2}\\\
&\qquad+c_{1}n_{0}\left[n_{1}+n_{-1}+2\sqrt{n_{1}n_{-1}}cos\theta_{r}\right]\Bigg{)},\end{split}$
(10)

where $e(\vec{r})$ is the energy density. Note that we are not interested in
vortex solutions, so in Eq.6-10 we have neglected spatial variation of the
sub-component phases assuming phases to be either constant or varying slowly.

For the stationary states, there is no temporal variation of the sub-component
number density and the sub-component phases i.e., the l.h.s of the Eq.6-9 can
be equated to zero. From the number density variation equations (Eq.6-7), one
can conclude that at least one of the sub-components should be empty to
satisfy the equations. Otherwise, if all the sub-components are populated,
then the relative phase can either be 0 or $\pi$. When all the sub-components
are populated, the stationary state where $\theta_{r}=0$ is satisfied is also
known as the phase-matched (PM) state whereas for $\theta_{r}=\pi$ one gets
the anti-phase-matched (APM) state. The sub-component phase equations (Eq.8-9)
for a particular stationary state can be solved to get the sub-component
number densities. Before going into that, let us get these equations in a non-
dimensional form.

The system considered is in quasi-one-dimensional harmonic confinement, i.e.,
the condensate is elongated along the x-axis. This means the trapping
frequency along the x-direction is much less than the geometric mean of the
trapping frequency along the other two directions i.e.,
$\omega_{x}<<\sqrt{\omega_{xy}}$, where
$\omega_{yz}=\sqrt{\omega_{y}\omega_{z}}$. One can scale the number density
and the interaction parameters as Kanjilal, Projjwal Kanti and Bhattacharyay,
A. (2022),

$c_{0}=2\pi l^{2}_{yz}l_{x}\lambda^{1D}_{0}\hbar\omega_{x},\quad c_{1}=2\pi
l^{2}_{yz}l_{x}\lambda^{1D}_{1}\hbar\omega_{x},$ (11) $u_{m}=2\pi
l^{2}_{yz}l_{x}n_{m},\quad r=l_{x}\zeta$ (12)

where, $l_{x}^{2}=\hbar/(m\omega_{x})$, $l_{yz}^{2}=\hbar/(m\omega_{yz})$, and
$N$ is the total number of particles in the condensate. As a result, the
parameters $\lambda^{1D}_{0}$, $\lambda^{1D}_{1}$, $\zeta$ and $u_{m}$ become
all dimensionless.

States |  Variation of density |  Energy density |  Restriction
---|---|---|---
(1,0,0) $F1$ |  $(\lambda_{0}+\lambda_{1})u(\zeta)=\mu^{\prime}+p^{\prime}-q^{\prime}-\zeta^{2}/2$ |  $\frac{\left[\zeta^{2}/2-p^{\prime}+q^{\prime}\right]\left[\mu^{\prime}+p^{\prime}-q^{\prime}-\zeta^{2}/2\right]}{(\lambda_{0}+\lambda_{1})}+\frac{\left[\mu^{\prime}+p^{\prime}-q^{\prime}-\zeta^{2}/2\right]^{2}}{2(\lambda_{0}+\lambda_{1})}$ |  $none$
(0,1,0) $P$ |  $\lambda_{0}u(\zeta)=\mu^{\prime}-\zeta^{2}/2$ |  $\frac{\zeta^{2}/2\left[\mu^{\prime}-\zeta^{2}/2\right]}{\lambda_{0}}+\frac{\left[\mu^{\prime}-\zeta^{2}/2\right]^{2}}{2\lambda_{0}}$ |  $none$
(0,0,1) $F2$ |  $(\lambda_{0}+\lambda_{1})u(\zeta)=\mu^{\prime}-p^{\prime}-q^{\prime}-\zeta^{2}/2$ |  $\frac{\left[\zeta^{2}/2+p^{\prime}+q^{\prime}\right]\left[\mu^{\prime}-p^{\prime}-q^{\prime}-\zeta^{2}/2)\right]}{(\lambda_{0}+\lambda_{1})}+\frac{\left[\mu^{\prime}-p^{\prime}-q^{\prime}-\zeta^{2}/2\right]^{2}}{2(\lambda_{0}+\lambda_{1})}$ |  $none$
(1,0,1) $AF$ |  $\lambda_{0}u(\zeta)=\mu^{\prime}-q^{\prime}-\zeta^{2}/2$ and $(u_{1}-u_{-1})\equiv F_{z}=\frac{p^{\prime}}{\lambda_{1}}$ |  $\frac{\left[\zeta^{2}/2+q^{\prime}\right]\left[\mu^{\prime}-q^{\prime}-\zeta^{2}/2\right]}{\lambda_{0}}+\frac{\left[\mu^{\prime}-q^{\prime}-\zeta^{2}/2\right]^{2}}{2\lambda_{0}}-\frac{p^{\prime 2}}{2\lambda_{1}}$ |  $none$
(1,1,1) $(A)PM$ |  $(\lambda_{0}+\lambda_{1})u(\zeta)=k_{1}-\zeta^{2}/2$ where, $k_{1}=\mu^{\prime}+\frac{(p^{\prime 2}-q^{\prime 2})}{2q^{\prime}}$ |  $\frac{\zeta^{2}/2\left[k_{1}-\zeta^{2}/2)\right]}{\lambda_{0}+\lambda_{1}}+\frac{\lambda_{1}}{2}\left[\frac{k_{1}-\zeta^{2}/2}{\lambda_{0}+\lambda_{1}}-\frac{p^{\prime 2}-q^{\prime 2}}{2q^{\prime}\lambda_{1}}\right]^{2}+\frac{\lambda_{0}}{2}\left[\frac{k_{1}-\zeta^{2}/2}{\lambda_{0}+\lambda_{1}}\right]^{2}$ |  $PM(|p^{\prime}|<|q^{\prime}|)$ $APM\quad(|p^{\prime}|>|q^{\prime}|)$

Table 1: The density and the energy density expressions corresponding to
different stationary states at $\lambda_{1}\neq 0$ obtained via T-F
approximation are shown here Kanjilal and Bhattacharyay (2020);
*Kanjilal_corri. All the parameters in this table are in dimensionless form.
One can use Eq.11-12 to convert expressions into dimensional forms. The energy
expressions and the density expressions for PM and APM states are identical.
However, PM and APM states are restricted in space where APM state exists if
the absolute value of the linear Zeeman term is higher than that of the
quadratic Zeeman term and PM state exists otherwise.

Three phase-stationary equations Eq.8-9 can now be written as,

$\begin{split}\bigg{\\{}&-\dfrac{1}{2}\dfrac{d^{2}}{d\zeta^{2}}+\dfrac{1}{2}\zeta^{2}+\lambda_{0}u-\mu^{\prime}\\\
&+\lambda_{1}\left(u_{1}+u_{-1}+2\sqrt{u_{-1}u_{1}}\cos\theta_{r}\right)\bigg{\\}}\sqrt{u_{0}}=0,\end{split}$
(13)
$\begin{split}\bigg{\\{}&-\dfrac{1}{2}\dfrac{d^{2}}{d\zeta^{2}}+\dfrac{1}{2}\zeta^{2}+\lambda_{0}u-\mu^{\prime}\pm\lambda_{1}\left(u_{1}-u_{-1}\right)\\\
&\pm p^{\prime}+q^{\prime}\bigg{\\}}\sqrt{u_{\pm
1}}+\lambda_{1}u_{0}\left(\sqrt{u_{\pm 1}}+\sqrt{u_{\mp
1}}\cos\theta_{r}\right)=0.\end{split}$ (14)

where, $\mu^{\prime}$, $p^{\prime}$, and $q^{\prime}$ are the scaled
parameters for chemical potential, the linear and quadratic Zeeman terms
respectively. The scaling is done by dividing these parameters with the factor
$(\hbar\omega_{x})$. The sub-component densities add up to provide the total
density of the system, i.e. $u=u_{1}+u_{0}+u_{-1}$.

If the kinetic term is negligible in comparison to interaction terms, then one
can use the T-F approximation and solve these three equations to get the sub-
component number densities for the stationary state. Thus, by substituting the
sub-component densities in Eq.10 one can get the energy or the energy density
of the stationary states. The T-F predicted ground states of our present
interest are detailed in Table 1.

##  III Multi-component stationary states:

Consider a quasi-one dimensional cigar-shaped harmonic confinement of trapping
frequency along the elongated direction $\omega_{x}=2\pi\times 50\hskip
2.84544ptHz$ and the geometric mean of the trapping frequencies along the
transverse direction $\omega_{yz}=2\pi\times 1261\hskip 2.84544ptHz$. In the
article Kanjilal, Projjwal Kanti and Bhattacharyay, A. (2022), it was shown
that for 1-D trapping geometry with the same trapping frequency, the T-F
approximation gives reasonably good results in predicting the sub-component
densities for single component stationary states for $N\geq 500$ in the
absence of magnetic field. There is reason to believe that T-F approximation
might produce fairly accurate results for multi-component stationary states in
presence of a small magnetic field if $N\geq 500$. For the case studies in
this paper, the number of condensate particles is fixed at $N=5000$ for which
T-F should produce even better results. However, in what follows we will show
that even at such a high particle number T-F results falter and corrections
are needed.

Our present focus is on the PM state which is a multi-component state that
appears as a ground state for a range of $p$ and $q$ values in condensates
with the ferromagnetic spin-spin interaction e.g., ${}^{87}Rb$. The anti-
ferromagnetic state that becomes the ground state for ${}^{23}Na$ where the
spin-spin interaction is of anti-ferromagnetic type Kawaguchi and Ueda (2012)
is also discussed in this paper. Note that, the quasi-one-dimensional
confinement is taken for convenience in numerical simulation. However, the
analytic formalism developed and validated is general and can be extended to
higher dimensional cases where numerical analysis could be problematic. In the
following, we show the results of the T-F study for these two cases, in order
to be able to draw a detailed comparison with the beyond T-F results later.

###  III.1 PM state: T-F study

For ${}^{87}Rb$ condensate with ferromagnetic type interaction, the numerical
values of the parameters corresponding to the trap geometry are
$l_{x}=1.53\hskip 2.84544pt\micro m$, $l_{yz}=0.30\hskip 2.84544pt\micro m$,
$\lambda_{0}=17.66\times 10^{-3}$ and $\lambda_{1}=-6.22\times 10^{-4}$. In
this section, to compare the T-F result with the numerical profile as a case
study the linear and quadratic Zeeman terms are fixed at $p^{\prime}=0.01$ and
$q^{\prime}=0.3$. The stationary state that is energetically favorable to be
ground state at these parameter values is the PM state. Note that, the
conclusions, however, are not dependent on a specific choice of parameter
values but rather will remain valid for a range of $p^{\prime}$ and
$q^{\prime}$ values for which the PM state is favorable as the ground state.

The sub-component number densities of the PM state are Kanjilal and
Bhattacharyay (2020); *Kanjilal_corri,

$u_{1}^{TF}=\dfrac{(p^{\prime}+q^{\prime})^{2}}{4q^{\prime
2}}\left[\dfrac{\mu^{\prime}+\dfrac{(p^{\prime 2}-q^{\prime
2})}{2q^{\prime}}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}+\lambda_{1}}+\dfrac{q^{\prime
2}-p^{\prime 2}}{2\lambda_{1}q^{\prime}}\right],$ (15)
$u_{-1}^{TF}=\dfrac{(p^{\prime}-q^{\prime})^{2}}{4q^{\prime
2}}\left[\dfrac{\mu^{\prime}+\dfrac{(p^{\prime 2}-q^{\prime
2})}{2q^{\prime}}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}+\lambda_{1}}+\dfrac{q^{\prime
2}-p^{\prime 2}}{2\lambda_{1}q^{\prime}}\right],$ (16)
$u_{0}^{TF}=\dfrac{(q^{\prime 2}-p^{\prime 2})}{2q^{\prime
2}}\left[\dfrac{\mu^{\prime}+\dfrac{(p^{\prime 2}-q^{\prime
2})}{2q^{\prime}}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}+\lambda_{1}}-\dfrac{q^{\prime
2}+p^{\prime 2}}{2\lambda_{1}q^{\prime}}\right],$ (17)

which can be compared with the numerically simulated profiles. For the purpose
of numerical simulation, imaginary time propagation in the split-step Fourier
method Kaur  _et al._ (2021) is used. The number density profiles as predicted
by the T-F are used as initial seeds.

(a) T-F profile of PM state and numerically obtained sub-component number
densities (b) T-F predicted domain-like structure and numerically obtained
profile of $u_{0}$ component

Figure 1: Subfig-(a): The T-F sub-component number densities for PM state is
compared with the numerical densities for the parameter value
$p^{\prime}=0.01$ and $q^{\prime}=0.3$. The T-F approximated $u_{0}$
expression corresponding to PM state (shown in dash-dot line) starts to
disagree with that of the numerical simulation beyond $|\zeta|>10.3$, which is
the T-F radius of the $u_{\pm 1}$ component. Subfig-(b): The T-F prediction of
the domain-like situation is plotted where, the $u_{0}$ component (solid line)
follows the T-F expression of the PM state when all the sub-components are
populated ($|\zeta|<10.3$ for this case), followed by a polar state-like
behavior. Note that discontinuity appears in the $u_{0}$ (inset) while the
numerical result is smooth (bubble markers).

Comparison of the sub-component densities of T-F profiles with the numerical
ones reveals very good agreement for the sub-components $u_{1}$ and $u_{-1}$.
The profile for the $u_{0}$ component also agrees in the high-density region
of the trap Fig.1LABEL:sub@subfig-1:_sub-component-TF, but starts to deviate
when the other components $u_{1}$ and $u_{-1}$ vanish which for this parameter
values is at around $\zeta^{TF}_{\pm 1}=\pm 10.3$. Note that, $\zeta^{TF}_{\pm
1}$ is the T-F radius of the $u_{\pm 1}$ components. Naturally, one might
expect that, beyond the point, $\zeta^{TF}_{\pm 1}$, the summation of the T-F
approximated sub-component number densities (i.e.,
$u_{tot}=u_{1}+u_{0}+u_{-1}$) will not agree with the numerical profile as
there is a significant mismatch in the $u_{0}$ component.

This is quite natural as the PM state is only defined as long as all the sub-
components are populated. However, PM state exists only until the $u_{1}$ and
$u_{-1}$ vanish (according to T-F calculations) at $\zeta=10.3$. So, within
the T-F picture, it can be said that near the center of the harmonic trap, the
PM state is occupied but as the $u_{\pm 1}$ goes to zero, the PM state ceases
to exist. For $|\zeta|>10.3$, only the $u_{0}$ component is present, which
signifies that it is the polar state that occupies the low-density region
Fig.1LABEL:sub@subfig-2:_PM-pol-density.

The sub-component number densities for such a construct can be given as,

$u_{1}^{TF}=\dfrac{(p^{\prime}+q^{\prime})^{2}}{4q^{\prime
2}}\left[\dfrac{\mu^{\prime}+\dfrac{(p^{\prime 2}-q^{\prime
2})}{2q^{\prime}}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}+\lambda_{1}}+\dfrac{q^{\prime
2}-p^{\prime 2}}{2\lambda_{1}q^{\prime}}\right],$ (18)
$u_{-1}^{TF}=\dfrac{(p^{\prime}-q^{\prime})^{2}}{4q^{\prime
2}}\left[\dfrac{\mu^{\prime}+\dfrac{(p^{\prime 2}-q^{\prime
2})}{2q^{\prime}}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}+\lambda_{1}}+\dfrac{q^{\prime
2}-p^{\prime 2}}{2\lambda_{1}q^{\prime}}\right],$ (19)
$u_{0}^{TF}=\begin{dcases}\dfrac{(q^{\prime 2}-p^{\prime 2})}{2q^{\prime
2}}\Bigg{[}\dfrac{\mu^{\prime}+\dfrac{(p^{\prime 2}-q^{\prime
2})}{2q^{\prime}}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}+\lambda_{1}}-\dfrac{q^{\prime
2}+p^{\prime 2}}{2\lambda_{1}q^{\prime}}\Bigg{]},&\\\ \hskip
142.26378pt\text{if }|\zeta|\leq\zeta_{\pm 1}^{TF}\\\
\dfrac{\mu^{\prime}_{polar}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}},&\text{otherwise},\end{dcases}$

where, $\zeta_{\pm 1}^{TF}$ is the Thomas-Fermi radius of the $u_{\pm 1}$
component for this 1-D geometry.

As can be seen in Fig.1LABEL:sub@subfig-2:_PM-pol-density this domain-like
explanation works really well as we compare the sub-component density with the
numerical $u_{0}$, but there will be a discontinuity at $|\zeta|=\zeta_{\pm
1}^{TF}$. The slope of analytical $u_{0}$ also changes drastically around this
point resulting in a lot of kinetic energy cost.

(a)  Figure 2: The total density obtained from numerical simulation matches
quite well with the expression Eq.20, which is the T-F expression
corresponding to the total density of the PM state. Surprisingly, this
expression of total number density is only valid as long as all the sub-
components are populated. (inset) The same expression matches with the
numerical total density even beyond the T-F radius of the $u_{\pm 1}$
component which is roughly at $|\zeta|\simeq 10.3$.

If we add the sub-component densities of the PM state given by expressions
Eq.15-17 one arrives at the total number density,

$u_{tot}^{PM}=\dfrac{\mu^{\prime}+\dfrac{(p^{2}-q^{2})}{2q}-\dfrac{1}{2}\zeta^{2}}{\lambda_{0}+\lambda_{1}},$
(20)

as long as all the sub-components are populated, a necessary condition for the
validity of PM state. As the number density cannot be negative, at
$|\zeta|>|\zeta_{\pm 1}^{TF}|$ where the $u_{\pm 1}^{TF}$ goes to zero, the
above expression is not valid beyond the TF radius of the $u_{\pm 1}$
components. Also, the PM state does not exist beyond this point.

In Fig.2, we notice that the total number density obtained from the numerical
simulation can be aptly described by the analytical expression Eq.20. The
element of surprise is that even beyond $|\zeta_{\pm 1}^{TF}|$, where
according to T-F the PM state ceases to exist, the numerical profile of total
number density agrees with Eq.20.

The above conclusion about the discontinuity in the $u_{0}$ component and
hence the limitations of T-F approximation can be further strengthened if we
compare numerical simulation and T-F prediction for a different spin-spin
interaction, $c_{1}\xrightarrow{}5c_{1}^{Rb}$ while keeping all other
parameter values the same. Experimentally, this is possible via Feshbach
resonance Papoular  _et al._ (2010).

(a) T-F predicted domain-like structure and numerically obtained profile of
$u_{0}$ component for a different spin-spin interction (b) T-F and numerically
obtained total number density

Figure 3: Subfig-(a): The T-F sub-component number density $u_{0}$ for PM
state for the same the parameter values $p^{\prime}=0.01$ and $q^{\prime}=0.3$
and $\lambda_{0}$ but for a different spin-dependent interaction coefficient
$c_{1}$ (see Eq.4) which is now 5 times of that of ${}^{87}Rb$ Kawaguchi and
Ueda (2012). The PM-like followed by polar-like $u_{0}$ has a distinct
discontinuity at $|\zeta_{\pm 1}^{TF}|=13.13$ (see inset). Subfig-(b): The
numerical total number density follows Eq. 20. (See inset) The analytical and
the numerical total number density agree across $|\zeta_{\pm 1}^{TF}|$ but a
visible deviation comes near the tail of the condensate where kinetic energy
contribution is the most significant.

For this choice of parameter values, the domain-like structure predicted by
T-F has a more visible discontinuity (Fig.3LABEL:sub@subfig-1:_sub-component-
TF-5c1) for the sub-component $u_{0}$ at $|\zeta_{\pm 1}^{TF}|=13.13$. The
total density for this domain-like structure will also have that same
discontinuity. Whereas, in Fig.3LABEL:sub@subfig-2:_TF-vs-num-density-5c1 we
verify that the numerical density and the Eq.20 are in agreement near the T-F
radius of $u_{\pm 1}$ and only deviates near the T-F radius of the $u_{0}$
component itself (inset of Fig.3LABEL:sub@subfig-2:_TF-vs-num-density-5c1).

To explain this, one has to remember that T-F approximation tells that the
$u_{\pm 1}$ sharply goes to zero at $|\zeta|=\zeta_{\pm 1}^{TF}$ but in
reality that cannot be the case. We need to include the kinetic energy term in
the analysis. As a result, near this point the Laplacian terms in Eq.14 cannot
be neglected. However, one can ignore the Laplacian term in Eq.13 as the
$u_{0}$ component is still in the high-density region near this point. As the
phase dynamics GP equations are coupled, one has to solve the Eq.13-14 by
keeping the Laplacian terms corresponding to $u_{\pm 1}$ components. This
shows that the T-F approximation is bound to produce an inaccurate description
of multi-component stationary states even in the so-called "T-F regime". This
warrants a need for an accurate analytical description of the multi-component
states appearing for the spin-1 BEC.

###  III.2 PM state: The variational method

In the article Kanjilal, Projjwal Kanti and Bhattacharyay, A. (2022), a multi-
modal variational method was introduced which incorporates the kinetic
contribution and produces the sub-component number-density profile of a
stationary state with great accuracy. The variational method (VM) was
developed for a spin-1 BEC in absence of a magnetic field which works really
well even for condensates with a small number of particles where the T-F
approximation was shown to be no longer valid. The PM state was also analyzed
using the multi-modal VM, but for $p=0$ and $q=0$ the PM state sub-component
densities are of the same spatial mode, therefore, reducing significant
complexity.

We extend the same procedure here in presence of a magnetic field. The number
densities of the PM state are multi-modal in spatial variation, thus, the
situation is much more complex but the extended VM tackles it with ease. In
this section we will briefly discuss the procedure; for a more elaborate
description see appendix A.

For the variational method to work, we need to estimate the sub-component
densities in the high-density region from the G-P equations Eq.13-14 by
neglecting the kinetic energy contribution. This is followed by an assumption
of the Gaussian-like tail in the low-density region. So the sub-component
densities are defined as,

$u^{in}_{\pm 1,0}=g_{\pm 1,0}(\mu^{\prime},\zeta),\hskip
8.5359pt\text{for}\quad|\zeta|\leq\zeta_{\pm 1,0}^{mat}$ (21)
$\begin{split}u^{out}_{\pm 1,0}=(a_{\pm 1,0}+c_{\pm 1,0}|\zeta|+d_{\pm
1,0}&\zeta^{2})exp\left(-\dfrac{\zeta^{2}}{b_{\pm 1,0}}\right)\\\
&\text{for}\quad|\zeta|\geq\zeta_{\pm 1,0}^{mat};\end{split}$ (22)

where $g_{\pm 1,0}(\mu^{\prime},\zeta)$ is the functional form of the sub-
component density $u_{\pm 1,0}$ near the center of the trap. In the low-
density region, we assume the number density (or the wave function) taking
into account the first few oscillator states.

To determine the unknown coefficients $a_{m}$, $b_{m}$, $c_{m}$, $d_{m}$ with
$m=\pm 1,0$, we impose a smooth matching condition at $|\zeta|=\zeta^{mat}$
where the $\sqrt{u_{\pm 1,0}}$ and their first three derivatives become equal
for low- and high-density expressions. The third derivative matching ensures a
smooth kinetic energy profile as well. From these four conditions, one can
determine all four coefficients appearing in Eq.22, which now become functions
of $\mu^{\prime}$ and $\zeta^{mat}_{\pm 1,0}$. The parameter $\mu^{\prime}$
can be expressed in terms of $\zeta^{mat}_{\pm 1,0}$ from the number
conservation equation (see appendix A). Following this step, the densities
(and as a result), the total energy can be written only as a function of the
matching points $\zeta^{mat}_{\pm 1,0}$ which are the only variables left.
These matching points can be found from the minimization of the total energy.
Once the matching points are obtained we can write the analytical expressions
for all the sub-component densities.

(a) The VM approximated sub-component number densities and the numerical
profile (b) The VM approximated and numerically obtained total number density

Figure 4: Subfig-(a): The sub-component number density of the PM state
obtained via the numerical and variational methods are plotted against the
distance from the trap center $\zeta$. The VM produces a very good analytical
profile that describes the numerical data quite well even in the tail part of
the condensate. The VM rules out any domain-like possibility and analytically
estimates the $u_{0}$ component (shown in the dash-dot line) that is quite
accurate in comparison to the numerical (circles) profile (see inset).
Subfig-(b): The variationally obtained total number density provides an
excellent analytical profile to match the numerical result near the core of
the trap as well as in the low-density region, where it gives an analytic
estimate of the condensate density which asymptotically goes to zero with the
increase of the distance $\zeta$.

Following the numerical evidence that the total number density follows Eq.20
(see Fig.2), it is sensible to apply the multi-modal VM for total number
density and the $u_{\pm 1}$. The $u_{0}$ expression can be later found out by
subtracting the sum of variational profiles of $n_{\pm 1}$ from the total
number density profile, i.e.,
$u_{0}^{VM}=u_{tot}^{VM}-u_{1}^{VM}-u_{-1}^{VM}$. So for simplicity, instead
of the $u_{0}$ component, we will focus on total number density which will
also provide a matching point $\zeta^{mat}_{tot}$. Note that, for the PM
state, the spatial mode for the $u_{1}$ and $u_{-1}$ are equivalent, so the
matching point for those two sub-components will be the same i.e.,
$\zeta^{mat}_{1}=\zeta^{mat}_{-1}$. Thus, the minimization of the total energy
in the two-dimensional parameter space of $\zeta^{mat}_{1}$ and
$\zeta^{mat}_{tot}$ will determine the energy itself as well as approximate
values of these parameters.

From the minimization of total energy for this particular case, we find that
the matching points are, $\zeta^{mat}_{1}=\zeta^{mat}_{-1}=7.98$ and
$\zeta^{mat}_{tot}=13.41$ and the corresponding value of the parameter,
$\mu^{\prime}=95.46$. Thus, the analytical number densities can be written as,

$u_{1}^{var}=\begin{dcases}1.5173(53.0977-0.5\zeta^{2}),\quad\text{if
}|\zeta|\leq 7.978\\\ \begin{split}1.5173&(9308.6737-2838.3002|\zeta|\\\
&+229.4058\zeta^{2})exp(-0.0642\zeta^{2})\end{split},&\\\ \hskip
128.0374pt\text{otherwise},\end{dcases}$ (23)
$u_{-1}^{var}=\begin{dcases}1.3278(53.0977-0.5\zeta^{2}),\quad\text{if
}|\zeta|\leq 7.978\\\ \begin{split}1.3278&(9308.6737-2838.3002|\zeta|\\\
&+229.4058\zeta^{2})exp(-0.0642\zeta^{2})\end{split},&\\\ \hskip
128.0374pt\text{otherwise},\end{dcases}$ (24)
$u_{tot}^{var}=\begin{dcases}5.6839(95.4569-0.5\zeta^{2}),\quad\text{if
}|\zeta|\leq 13.412\\\ \begin{split}5.6839&(8.895\times 10^{24}-1.3539\times
10^{24}|\zeta|\\\ &+5.1533\times
10^{22}\zeta^{2})exp(-0.2692\zeta^{2})\end{split},&\\\ \hskip
128.0374pt\text{otherwise},\end{dcases}$ (25)

where the numbers are rounded up to four decimal places as it may cause a
large rounding error in the density expressions. Note that, in Eq.25 the
coefficients of the total number density might look very large, but for
$|\zeta|\geq 13.412$ where the expression is valid, the contribution from the
exponential part is so small that the combined contribution asymptotically
goes to zero at large distances and matches quite accurately with the
numerical profile (see the inset of Fig.4LABEL:sub@subfig-2:_var-vs-num-
density)

We have made a case study of the PM state which is a multi-component
stationary state that becomes the ground state for a range of linear and
quadratic Zeeman strengths. For the purpose of comparison with numerical
simulation, we have chosen 1-D harmonic trapping and particular values of
${}^{\prime}p$, $q^{\prime}$, and the number of condensate particles $N$.

Note that the VM is an approximation scheme that works really well in
estimating the sub-component number densities (also the mean fields) which
produces a very good estimation of the vector order parameter of the spin-1
system. Like other approximate methods, it has some limitations as well. For
example, at a large distance from the center of the trap (very large $\zeta$),
where the total density $u_{tot}^{VM}$ and $u_{\pm 1}^{VM}$ are very close to
zero and can be considered negligible, we find that the total density is
slightly lesser than the combined contribution of the $\pm 1$ sub-components
hence, making $u^{VM}_{0}$ slightly negative which is not physical. For this
reason, we have taken the contribution up to a large $\zeta$ after which we
assume that $u_{0}$ goes to zero. Thus the kinetic energy contribution is
included and considered up to a large distance without discontinuity.

In the next section, we will do a brief case study on the anti-ferromagnetic
state which is the other possible multi-component stationary state that
becomes the ground state for ${}^{23}Na$.

###  III.3 Anti-ferromagnetic state

For the ${}^{23}Na$-condensate which has an anti-ferromagnetic type of spin-
spin interaction ($c_{1}>0$), we set the same trapping frequencies
corresponding to 1-D confinement as mentioned earlier for the ferromagnetic
type condensate. The oscillator length in elongated direction
$l_{x}=2.97\hskip 2.84544pt\micro m$ and in the transverse direction
$l_{yz}=0.59\hskip 2.84544pt\micro m$. Note that, although we consider the
same trapping geometry, the oscillator length scale for ${}^{23}Na$ and
${}^{87}Rb$ condensates are different due to the different masses of the
species. The non-dimensional spin-independent and spin-spin interaction
parameters are $\lambda_{0}=46.16\times 10^{-3}$ and $\lambda_{1}=7.43\times
10^{-4}$ corresponding to the values given in Kawaguchi and Ueda (2012). For a
range of linear and quadratic Zeeman terms, the anti-ferromagnetic (AF) state
is found to be favorable to be the ground state. For the purpose of numerical
study, we focus on the case these terms are fixed at $p^{\prime}=0.2$ and
$q^{\prime}=-0.5$.

As long as the $u_{1}$ and $u_{-1}$ sub-components are non-zero, the T-F
approximation gives an estimation of the total as well as sub-component number
densities (see Table 1),

$u^{TF}_{1}=\dfrac{\mu^{\prime}-q^{\prime}-\zeta^{2}/2}{2\lambda_{0}}+\dfrac{p^{\prime}}{2\lambda_{1}},$
(26)
$u^{TF}_{-1}=\dfrac{\mu^{\prime}-q^{\prime}-\zeta^{2}/2}{2\lambda_{0}}-\dfrac{p^{\prime}}{2\lambda_{1}}.$
(27)

As the value of $p^{\prime}$ is chosen as positive (and $\lambda_{1}$ is
positive as well) the $u_{-1}$ component goes to zero much faster than the
other one. So beyond the T-F radius of the $u_{-1}$ component the AF state
ceases to exist. As the $u_{1}$ component is still non-zero, so according to
T-F, the situation is domain-like with AF state at the center of the trap
followed by the ferromagnetic state outside ( for $|\zeta|>\zeta^{TF}_{-1}$).

Just like the PM state that we have seen earlier, the numerical simulation
does not vindicate the domain-like prediction. Rather the AF state is found to
be present for all values of $\zeta$. The VM can be applied to get an
analytical explanation of the situation for this multi-component state.

(a)  Figure 5: Sub-component density expressions obtained via VM are compared
with the numerical profile for the anti-ferromagnetic state with the Zeeman
terms fixed at $p^{\prime}=0.2$ and $q^{\prime}=-0.5$. The density expressions
obtained from VM for $u_{\pm 1}$ are shown in dashed and solid lines
respectively, which agrees quite well with the numerical results both in high-
and low-density regions.

Following the same procedure as described for PM state one can apply the VM
for the total density and the sub-component density $u_{-1}$. The $u_{0}$
component is empty and there will be two matching points $\zeta_{tot}^{mat}$
and $\zeta_{-1}^{mat}$ that are to be found out from the minimization of the
total energy. For the previously mentioned $p^{\prime}$ and $q^{\prime}$
values, we get these two matching points to be, $\zeta_{tot}^{mat}=8.36$ and
$\zeta_{-1}^{mat}=5.79$. These also produce the analytical formulae of the
total density and the $u_{-1}$ components. The density expression for $u_{1}$
component can be obtained by subtracting the other sub-component density from
the total density.

The VM shows that the domain-like situation that the T-F predicts is incorrect
and justifies the fact that the kinetic energy terms cannot be neglected at
least near the T-F radius for the sub-component which is of smaller density,
in this case, the $u_{-1}$ component. The VM also produces a low-density
expression of the $u_{-1}$ component which has a small but non-zero presence
beyond $u_{-1}^{TF}$. Thus, it is only the AF state that is present for all
regions of space. Moreover, the analytic number density expressions obtained
from VM corresponding to each sub-component are in fair agreement with the
numerically obtained profiles (see Fig.5).

##  IV Discussion

In this article, we have discussed the inaccuracy of the T-F as the most
widely used analytical method to describe the multi-component states of a
trapped spin-1 BEC in presence of the magnetic field. The T-F approximation
leads to a wrong conclusion about the existence of the domain-like structure
while explaining the ground state structure involving multi-component
stationary states even in the so-called "T-F regime", where the density of the
condensate is high enough to neglect the kinetic energy contribution in
comparison to the interaction terms. We have extended the multi-modal
variational method introduced in Kanjilal, Projjwal Kanti and Bhattacharyay,
A. (2022) to the presence of the magnetic field. We have shown that the
variational method successfully describes the density distribution of the
multi-component states which is also backed by numerical simulation. The VM
properly captures the multi-component states because it can estimate the tail
part of the components with smaller densities where T-F cuts it off. These
tail parts of the densities reside in a region where the kinetic energy
contributions are most significant and cannot be neglected even if it is
small.

For this study, we focus on multi-component states with fixed values of linear
and quadratic Zeeman terms. We also considered a fixed value of the total
number of condensate particles, $N$. However, this procedure is applicable for
different values of these parameters. The VM can be easily implemented for 3-D
harmonic trapping to get the details of all the stationary states including
all the multi-component states. Whereas, doing numerical simulation for a 3-D
system is well-known to be computationally expensive. Note that for single-
component stationary states (like ferromagnetic and polar states), the T-F
gives fairly accurate results in comparison to the numerical simulation in the
T-F regime. Still, the variational method finds its importance which
incorporates the kinetic energy contribution and predicts the condensate
profile even in the low-density region with fair accuracy. Specifically for
condensates with a smaller number of particles where the kinetic energy term
is comparable to the interaction energy, T-F is in general, not applicable and
the VM is indispensable as an analytical method Kanjilal, Projjwal Kanti and
Bhattacharyay, A. (2022). This is true even for non-zero contributions of $p$
and $q$.

The VM is obviously an approximation scheme but it provides a good analytical
profile of the condensate. Note that, many studies require a numerically exact
ground state profile of a spin-1 condensate under trapping in presence of a
magnetic field. This analytical method provided here would also find its
relevance in such cases if one uses the VM approximated density profile as an
initial condition for quick convergence to the numerical profile.

The ground state structure and the phase diagram of a trapped spin-1
condensate in the $p$, $q$ parameter space are of fundamental relevance. The
total energy corresponding to different stationary states can be compared
which reveals the stationary state that becomes the ground state for different
values of $p$ and $q$. Thus, an accurate analytical description of all the
stationary states, specifically the multi-component states is required.
Therefore, the method involving the variational procedure shown here is of
fundamental importance and becomes a prerequisite to such a tedious study.

##  V Acknowledgement

PKK would like to thank the Council of Scientific and Industrial Research
(CSIR), India for providing funding during this research. The support and the
resources provided by PARAM Brahma Facility under the National Supercomputing
Mission, Government of India at the Indian Institute of Science Education and
Research; Pune are gratefully acknowledged.

Data availability statement: The datasets generated and analyzed during the
current study are available from the corresponding author at a reasonable
request.

## References

  * Kanjilal, Projjwal Kanti and Bhattacharyay, A. (2022) Kanjilal, Projjwal Kanti and Bhattacharyay, A., Eur. Phys. J. Plus 137, 547 (2022).
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995).
  * Anderson  _et al._ (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995).
  * Bradley  _et al._ (1995) C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995).
  * Giorgini  _et al._ (1996) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 54, R4633 (1996).
  * Dalfovo and Stringari (1996) F. Dalfovo and S. Stringari, Phys. Rev. A 53, 2477 (1996).
  * Stringari (1996) S. Stringari, Phys. Rev. Lett. 77, 2360 (1996).
  * Dupont  _et al._ (2021) N. Dupont, G. Chatelain, L. Gabardos, M. Arnal, J. Billy, B. Peaudecerf, D. Sugny, and D. Guéry-Odelin, PRX Quantum 2, 040303 (2021).
  * Senaratne  _et al._ (2018) R. Senaratne, S. V. Rajagopal, T. Shimasaki, P. E. Dotti, K. M. Fujiwara, K. Singh, Z. A. Geiger, and D. M. Weld, Nat Commun 9, 2065 (2018).
  * Yamamoto and Takahashi (2016) Y. Yamamoto and Y. Takahashi, “Bose-einstein condensation: A platform for quantum simulation experiments,” in  _Principles and Methods of Quantum Information Technologies_, edited by Y. Yamamoto and K. Semba (Springer Japan, Tokyo, 2016) pp. 265–307.
  * Bravo  _et al._ (2015) T. Bravo, C. Sabín, and I. Fuentes, EPJ Quantum Technology 2, 3 (2015).
  * Cirac and Zoller (2012) J. I. Cirac and P. Zoller, Nature Physics 8, 264 (2012).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * Becker  _et al._ (2018) D. Becker, M. D. Lachmann, S. T. Seidel, H. Ahlers, A. N. Dinkelaker, J. Grosse, O. Hellmig, H. Müntinga, V. Schkolnik, T. Wendrich, A. Wenzlawski, B. Weps, R. Corgier, T. Franz, N. Gaaloul, W. Herr, D. Lüdtke, M. Popp, S. Amri, H. Duncker, M. Erbe, A. Kohfeldt, A. Kubelka-Lange, C. Braxmaier, E. Charron, W. Ertmer, M. Krutzik, C. Lämmerzahl, A. Peters, W. P. Schleich, K. Sengstock, R. Walser, A. Wicht, P. Windpassinger, and E. M. Rasel, Nature 562, 391 (2018).
  * Szigeti  _et al._ (2020) S. S. Szigeti, S. P. Nolan, J. D. Close, and S. A. Haine, Phys. Rev. Lett. 125, 100402 (2020).
  * Gupta  _et al._ (2002) S. Gupta, K. Dieckmann, Z. Hadzibabic, and D. E. Pritchard, Phys. Rev. Lett. 89, 140401 (2002).
  * Jamison  _et al._ (2011) A. O. Jamison, J. N. Kutz, and S. Gupta, Phys. Rev. A 84, 043643 (2011).
  * Hanneke  _et al._ (2008) D. Hanneke, S. Fogwell, and G. Gabrielse, Phys. Rev. Lett. 100, 120801 (2008).
  * Dutta  _et al._ (2022) P. Dutta, S. S. Maurya, K. Patel, K. Biswas, J. Mangaonkar, S. Sarkar, and U. D. Rapol, Journal of the Indian Institute of Science (2022), 10.1007/s41745-022-00335-8.
  * Ensher  _et al._ (1996a) J. R. Ensher, D. S. Jin, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 77, 4984 (1996a).
  * Ensher  _et al._ (1996b) J. R. Ensher, D. S. Jin, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 77, 4984 (1996b).
  * Ketterle  _et al._ (1996) W. Ketterle, M. R. Andrews, K. B. Davis, D. S. Durfee, D. M. Kurn, M. O. Mewes, and N. J. van Druten, Physica Scripta 1996, 31 (1996).
  * Stamper-Kurn  _et al._ (1998) D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle, Phys. Rev. Lett. 80, 2027 (1998).
  * Chang  _et al._ (2004) M.-S. Chang, C. D. Hamley, M. D. Barrett, J. A. Sauer, K. M. Fortier, W. Zhang, L. You, and M. S. Chapman, Phys. Rev. Lett. 92, 140403 (2004).
  * (25) D. M. Stamper-Kurn and W. Ketterle, in  _Les Houches - Ecole d’Ete de Physique Theorique_ (Springer Berlin Heidelberg) pp. 139–217.
  * Kawaguchi  _et al._ (2006) Y. Kawaguchi, H. Saito, and M. Ueda, Phys. Rev. Lett. 97, 130404 (2006).
  * Leanhardt  _et al._ (2003) A. E. Leanhardt, Y. Shin, D. Kielpinski, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 90, 140403 (2003).
  * Timmermans (1998) E. Timmermans, Phys. Rev. Lett. 81, 5718 (1998).
  * Świsłocki and Matuszewski (2012) T. Świsłocki and M. Matuszewski, Phys. Rev. A 85, 023601 (2012).
  * Xi  _et al._ (2011) K.-T. Xi, J. Li, and D.-N. Shi, Phys. Rev. A 84, 013619 (2011).
  * Bandyopadhyay  _et al._ (2017) S. Bandyopadhyay, A. Roy, and D. Angom, Phys. Rev. A 96, 043603 (2017).
  * Isoshima  _et al._ (1999) T. Isoshima, K. Machida, and T. Ohmi, Phys. Rev. A 60, 4857 (1999).
  * Matuszewski  _et al._ (2008) M. Matuszewski, T. J. Alexander, and Y. S. Kivshar, Phys. Rev. A 78, 023632 (2008).
  * Matuszewski  _et al._ (2009) M. Matuszewski, T. J. Alexander, and Y. S. Kivshar, Phys. Rev. A 80, 023602 (2009).
  * Ho and Shenoy (1996) T.-L. Ho and V. B. Shenoy, Phys. Rev. Lett. 77, 3276 (1996).
  * Bhuvaneswari  _et al._ (2018) S. Bhuvaneswari, K. Nithyanandan, and P. Muruganandam, Journal of Physics Communications 2, 025008 (2018).
  * Wen  _et al._ (2012) L. Wen, W. M. Liu, Y. Cai, J. M. Zhang, and J. Hu, Phys. Rev. A 85, 043602 (2012).
  * Liu (2009) Z. Liu, Journal of Mathematical Physics 50, 102104 (2009), https://doi.org/10.1063/1.3243875 .
  * Zhu and Li (2017) L. Zhu and J. Li, Modern Physics Letters B 31, 1750215 (2017), https://doi.org/10.1142/S0217984917502153 .
  * Tojo  _et al._ (2010) S. Tojo, Y. Taguchi, Y. Masuyama, T. Hayashi, H. Saito, and T. Hirano, Phys. Rev. A 82, 033609 (2010).
  * Li  _et al._ (2018) J. Li, Y.-M. Yu, K.-J. Jiang, and W.-M. Liu, arXiv preprint arXiv:1802.00138 (2018).
  * Lee  _et al._ (2016) K. L. Lee, N. B. Jørgensen, I.-K. Liu, L. Wacker, J. J. Arlt, and N. P. Proukakis, Phys. Rev. A 94, 013602 (2016).
  * Sabbatini  _et al._ (2011) J. Sabbatini, W. H. Zurek, and M. J. Davis, Phys. Rev. Lett. 107, 230402 (2011).
  * Gautam and Adhikari (2014) S. Gautam and S. K. Adhikari, Phys. Rev. A 90, 043619 (2014).
  * Gautam and Angom (2011) S. Gautam and D. Angom, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 025302 (2011).
  * Saito and Ueda (2005) H. Saito and M. Ueda, Phys. Rev. A 72, 023610 (2005).
  * Kanjilal and Bhattacharyay (2020) P. K. Kanjilal and A. Bhattacharyay, Physica Scripta 95, 045702 (2020).
  * Kanjilal and Bhattacharyay (2022) P. K. Kanjilal and A. Bhattacharyay, Physica Scripta 97, 129501 (2022).
  * Ueda (2014) M. Ueda, Reports on Progress in Physics 77, 122401 (2014).
  * Nistazakis  _et al._ (2008) H. E. Nistazakis, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and R. Carretero-González, Phys. Rev. A 77, 033612 (2008).
  * Bersano  _et al._ (2018) T. M. Bersano, V. Gokhroo, M. A. Khamehchi, J. D’Ambroise, D. J. Frantzeskakis, P. Engels, and P. G. Kevrekidis, Phys. Rev. Lett. 120, 063202 (2018).
  * Katsimiga  _et al._ (2021) G. C. Katsimiga, S. I. Mistakidis, P. Schmelcher, and P. G. Kevrekidis, New Journal of Physics 23, 013015 (2021).
  * Hegde  _et al._ (2022) G. Hegde, S. M. Jose, and R. Nath, Phys. Rev. A 106, 043307 (2022).
  * Mistakidis  _et al._ (2022) S. I. Mistakidis, A. G. Volosniev, R. E. Barfknecht, T. Fogarty, T. Busch, A. Foerster, P. Schmelcher, and N. T. Zinner, “Cold atoms in low dimensions – a laboratory for quantum dynamics,” (2022).
  * Ho (1998) T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998).
  * Yi  _et al._ (2002) S. Yi, O. E. Müstecaplıoğlu, C. P. Sun, and L. You, Phys. Rev. A 66, 011601 (2002).
  * Zhang  _et al._ (2003) W. Zhang, S. Yi, and L. You, New Journal of Physics 5, 77 (2003).
  * Gautam and Adhikari (2015) S. Gautam and S. K. Adhikari, Phys. Rev. A 92, 023616 (2015).
  * Phuc  _et al._ (2011) N. T. Phuc, Y. Kawaguchi, and M. Ueda, Phys. Rev. A 84, 043645 (2011).
  * Kawaguchi  _et al._ (2012) Y. Kawaguchi, N. T. Phuc, and P. B. Blakie, Phys. Rev. A 85, 053611 (2012).
  * Serrano-EnsÃ¡stiga and Mireles (2022) E. Serrano-EnsÃ¡stiga and F. Mireles, “Spinor bose-einstein condensates: self-consistent symmetries and characterization,” (2022).
  * Mittal  _et al._ (2020) K. M. Mittal, S. I. Mistakidis, P. G. Kevrekidis, and P. Schmelcher, Phys. Rev. A 102, 013302 (2020).
  * Stenger  _et al._ (1998) J. Stenger, S. Inouye, D. Stamper-Kurn, H.-J. Miesner, A. Chikkatur, and W. Ketterle, Nature 396, 345 (1998).
  * Kawaguchi and Ueda (2012) Y. Kawaguchi and M. Ueda, Physics Reports 520, 253 (2012), spinor Bose–Einstein condensates.
  * Ohmi and Machida (1998) T. Ohmi and K. Machida, Journal of the Physical Society of Japan 67, 1822 (1998).
  * Kaur  _et al._ (2021) P. Kaur, A. Roy, and S. Gautam, Computer Physics Communications 259, 107671 (2021).
  * Papoular  _et al._ (2010) D. J. Papoular, G. V. Shlyapnikov, and J. Dalibard, Phys. Rev. A 81, 041603 (2010).

##  Appendix A Variational Method

To work with the non-zero contribution of $p$ and $q$, we extend the
variational method introduced in the article Kanjilal, Projjwal Kanti and
Bhattacharyay, A. (2022). First, we will present the general method in a brief
manner and then we will implement that for the PM state and the AF state which
are the multi-component states of interest to this article.

To implement the variational method, in presence of quasi-one-dimensional
harmonic trapping, one needs to solve the GP equations Eq.13-14 by getting rid
of the kinetic energy terms which will yield the sub-component densities in
the high-density region near the center of the trap, where the kinetic energy
can be neglected in comparison to the interaction terms. The sub-component
densities will get a functional form specific to different stationary states.
Next, we assume that near the low-density region where the kinetic energy is
of relevance and the interaction terms are very small, the mean fields can be
described in terms of the first few oscillator states. Thus the density
distributions can be written as,

$u^{in}_{\pm 1,0}=g_{\pm 1,0}(\mu^{\prime},\zeta),\hskip
8.5359pt\text{for}\quad|\zeta|\leq\zeta_{\pm 1,0}^{mat}$ (28)
$\begin{split}u^{out}_{\pm 1,0}=(a_{\pm 1,0}+c_{\pm 1,0}|\zeta|+d_{\pm
1,0}&\zeta^{2})exp\left(-\dfrac{\zeta^{2}}{b_{\pm 1,0}}\right)\\\
&\text{for}\quad|\zeta|\geq\zeta_{\pm 1,0}^{mat},\end{split}$ (29)

where, $u^{in(out)}_{\pm 1,0}$ is the sub-component densities in the
high-(low-)density region. Now we impose the condition that for each sub-
component the low-density $\sqrt{u^{out}_{\pm 1,0}}$ and the high-density
$\sqrt{u^{in}_{\pm 1,0}}$ expressions match at a point $\zeta^{mat}$. Not only
do they match but their first three derivatives also match. These four
constraints provide the four unknowns $a$, $b$, $c$, and $d$ for each sub-
components in terms of the matching points and the parameter $\mu^{\prime}$.
Note that, imposing the matching condition up to three derivatives also gives
a smooth profile of the corresponding kinetic energy.

Once all the coefficients in Eq.29 are known, the sub-component density
profile only depends on the parameter $\mu^{\prime}$ and the matching points.
The parameter $\mu^{\prime}$ can be obtained as a function of the matching
points from,

$\begin{split}\sum_{m=-1}^{1}\Bigg{[}\int_{0}^{\zeta_{m}^{mat}}&u_{m}^{in}(\mu^{\prime},\zeta)d\zeta\\\
&+\int_{\zeta_{m}^{mat}}^{\infty}u_{m}^{out}(\mu^{\prime},\zeta,\zeta_{m}^{mat})d\zeta\Bigg{]}=N,\end{split}$
(30)

where N is the total number of condensate particles. Note that, one might
expect that the right side should be $N/2$ as the integration is running in
only one direction from the center of the trap, but it is N in the right side
due to the Eq.11-12 which we used to write the GP equation in non-dimensional
form. From this Eq.30 the parameter $\mu^{\prime}$ can be written as a
function of the matching points for a particular N.

Thus the sub-component number densities and hence the total energy of a
stationary state (Eq.10) also becomes the function of the matching points
only. Now minimization of the total energy in the parameter space of the
matching points one can get the matching points as well as the total energy
itself.

###  A.1 Variational method for the PM state:

For the PM state, all the sub-components are populated followed by the phase
matching condition, i.e., the relative phase being $\theta_{r}=0$. One can
solve the phase stationary equations (Eq.13-14) by ignoring the kinetic part
to get the sub-component densities in the high-density region as,

$u^{in}_{m}=k_{m}\Bigg{[}\dfrac{\mu^{\prime}_{m}-\zeta^{2}/2}{\lambda_{0}+\lambda_{1}}\Bigg{]},$
(31)

where,

$\begin{split}&k_{1}=\dfrac{(p^{\prime}+q^{\prime})^{2}}{4q^{\prime 2}}\\\
&k_{0}=\dfrac{q^{\prime 2}-p^{\prime 2}}{2q^{\prime 2}}\\\
&k_{-1}=\dfrac{(p^{\prime}-q^{\prime})^{2}}{4q^{\prime 2}}\end{split}$ (32)

and,

$\begin{split}&\mu^{\prime}_{\pm
1}=\mu^{\prime}_{eff}+(\lambda_{0}+\lambda_{1})\dfrac{q^{\prime 2}-p^{\prime
2}}{2\lambda_{1}q^{\prime}},\\\
&\mu^{\prime}_{0}=\mu^{\prime}_{eff}-(\lambda_{0}+\lambda_{1})\dfrac{q^{\prime
2}+p^{\prime 2}}{2\lambda_{1}q^{\prime}},\\\
&\mu^{\prime}_{eff}=\mu^{\prime}+\dfrac{p^{\prime 2}-q^{\prime
2}}{2q^{\prime}}.\end{split}$ (33)

Once, we write the sub-component densities in the high-density region in such
a compact form like Eq.31 one can apply the four matching conditions mentioned
earlier-. This would allow us to get the four unknown coefficients in the low-
density expression for each sub-component,

$a_{m}=\dfrac{1}{-8\mu^{\prime}_{m}+4\zeta_{m}^{2}}\Bigg{(}\mu^{\prime}_{m}\Big{(}-56\mu^{\prime}_{m}+70\zeta_{m}^{2}+4\kappa_{m}\Big{)}-3\zeta_{m}^{2}\Big{(}14\zeta_{m}^{2}+\kappa_{m}-6\mu^{\prime}_{m}\Big{)}\Bigg{)}exp\Bigg{(}\dfrac{12\zeta_{m}^{2}}{\kappa_{m}}\Bigg{)},$
(34)

$b_{m}=\dfrac{\kappa_{m}}{12},$ (35)

$c_{m}=\dfrac{48\zeta_{m}^{3}\Big{(}-12\mu^{\prime}_{m}+6\zeta_{m}^{2}+\kappa_{m}\Big{)}}{\kappa_{m}^{2}}exp\Bigg{(}\dfrac{12\zeta_{m}^{2}}{\kappa_{m}}\Bigg{)},$
(36)
$\begin{split}d_{m}=\dfrac{1}{2\zeta_{m}^{2}\Big{(}-2\mu^{\prime}_{m}+\zeta_{m}^{2}\Big{)}}\Bigg{(}-6(\mu^{\prime}_{m})^{2}-\zeta_{m}^{2}&\Big{(}\kappa_{m}+13\zeta_{m}^{2}-6\mu^{\prime}_{m}\Big{)}\\\
&+\mu^{\prime}_{m}\Big{(}14\zeta_{m}^{2}+\kappa_{m}-6\mu^{\prime}_{m}\Big{)}\Bigg{)}exp\Bigg{(}\dfrac{12\zeta_{m}^{2}}{\kappa_{m}}\Bigg{)},\end{split}$
(37)

where, $\zeta_{m}$ is an abbreviation for the matching point $\zeta_{m}^{mat}$
and,

$\begin{split}\kappa_{m}=6\mu^{\prime}_{m}&-9\zeta_{m}^{2}\\\
&+\sqrt{36\mu_{m}^{2}-12\mu_{m}\zeta_{m}^{2}+33\zeta_{m}^{4}},\end{split}$
(38)

given one writes the sub-component densities in the low-density region as,

$u^{out}_{m}=\dfrac{k_{m}}{\lambda_{0}+\lambda_{1}}\Big{(}a_{m}+c_{m}|\zeta|+d_{m}\zeta^{2}\Big{)}exp\Bigg{(}-\dfrac{\zeta^{2}}{b_{m}}\Bigg{)}.$
(39)

Now applying Eq.30 one can find the parameter $\mu^{\prime}$ for different
values of $\zeta_{m}$. Note that, as $\mu^{\prime}_{1}=\mu^{\prime}_{-1}$ (see
Eq.33) the matching points are the same for these two components, i.e.,
$\zeta_{1}^{mat}=\zeta_{-1}^{mat}$. Thus, the total energy (Eq.10) for the PM
state becomes only a function of the matching points.

Note that, from a physical perspective, the high-density expressions of
$u_{0}$ and $u_{\pm 1}$ were found by getting rid of the kinetic terms in the
Eq.13-14. So the high-density expressions given in Eq.31 are true as long as
all the sub-components are in the high-density region. If we focus on the case
presented in the article, for $p^{\prime}=0.01$ and $q^{\prime}=0.3$, near T-F
radius the $u_{\pm 1}$ sub-component does not follow the high-density
expressions. As a result, one can neglect the derivative term in Eq.13 but
cannot neglect the same in Eq.14. This precisely makes the multi-component
stationary states beyond the reach of T-F approximation as long as the sub-
components are multi-modal in nature.

Precisely because of the reason stated above one has to shift the focus toward
the total density. So, instead of using the high-density expression
$u^{in}_{0}$ we will use the total density expression,

$u_{tot}^{in}=k_{tot}\Bigg{[}\dfrac{\mu^{\prime}_{eff}-\zeta^{2}/2}{\lambda_{0}+\lambda_{1}}\Bigg{]},$
(40)

written in the same fashion as Eq.31, where $k_{tot}=1$ and the $\mu_{eff}$
can be found in Eq.33. Now instead of putting $m=0$, one can put $m=tot$ in
Eq.34-37 which would provide the total density expression in the low-density
region (putting $m=tot$ in Eq.39). The $u_{0}$ component can be found by
subtracting the other two component densities from the total density.

Now the total energy of the PM state only becomes a function of the two
matching points, $\zeta_{tot}^{mat}$ and $\zeta_{\pm 1}^{mat}$. Minimizing the
total energy with respect to these two parameter variations one can get the
matching points and the total energy itself. And as the matching points are
found, the analytical density expressions are also obtained.

###  A.2 Variational method for the AF state:

When the $u_{1}$ and the $u_{-1}$ components are populated (even though they
are unequally populated for non-zero values of $p^{\prime}$ (see Table 1)) the
stationary state is referred to as the anti-ferromagnetic or in short, AF
state. When both the sub-components are in the high-density regions, one can
write the densities by neglecting the derivative terms in the phase equations
as,

$u^{in}_{\pm 1}=k_{\pm 1}\Bigg{[}\dfrac{\mu^{\prime}_{\pm
1}-\zeta^{2}/2}{\lambda_{0}}\Bigg{]},$ (41)

where, $k_{\pm 1}=1/2$ and,

$\mu^{\prime}_{\pm
1}=\mu^{\prime}-q^{\prime}\pm\dfrac{\lambda_{0}}{\lambda_{1}}p^{\prime}.$ (42)

The sub-component density $u_{0}$ is zero throughout.

Note that, the high energy expressions are valid as long as both the sub-
component density is high enough so that the derivative terms can be safely
ignored. But for $p^{\prime}=0.2$ and $q^{\prime}=-0.5$ (the case we
discussed), the $u_{-1}$ component has a lesser T-F radius than the other
component. As a result of it near the low-density region of the $u_{-1}$
component the high-density expression of the $u_{1}$ component would be
invalid, for the reasons stated earlier. So, we will take the total density
and the $u_{-1}$ component to implement the variational method.

In the high-density region, the total density can be written as,

$u^{in}_{tot}=k_{tot}\Bigg{[}\dfrac{\mu^{\prime}_{tot}-\zeta^{2}/2}{\lambda_{0}}\Bigg{]},$
(43)

where $k_{tot}=1$. Now if one writes the low-density expressions as,

$u^{out}_{m}=\dfrac{k_{m}}{\lambda_{0}}\Big{(}a_{m}+c_{m}|\zeta|+d_{m}\zeta^{2}\Big{)}exp\Bigg{(}-\dfrac{\zeta^{2}}{b_{m}}\Bigg{)}$
(44)

where the subscript can be $m=-1,tot$. Now the coefficients will have the same
expressions as Eq.34-38.

Following the same method as explained earlier one can minimize the total
energy corresponding to this stationary state in the parameter space of the
matching points $\zeta_{tot}$ and $\zeta_{-1}$ which provides the analytical
form of the full profile of the condensate in terms of the total density and
the $u_{-1}$ component. By subtracting the $u_{-1}$ component from the total
density profile one can get to the $u_{1}$ component.