††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Breakup of rotating asymmetric quartic-quadratic trapped condensates

Leonardo Brito1 Alex Andriati1 Lauro Tomio2 Arnaldo Gammal1 1Instituto de
Física, Universidade de São Paulo, 05508-090 São Paulo, Brazil.
2Instituto de Física Teórica, Universidade Estadual Paulista, 01156-970 São
Paulo, SP, Brazil.

###### Abstract

The threshold conditions for a rotating pancakelike asymmetric quartic-
quadratic confined condensate to break in two localized fragments, as well as
to produce giant vortex at the center within the vortex-pattern distributions,
are investigated within the Thomas-Fermi (TF) approximation and full-numerical
solution of the corresponding Gross-Pitaevskii (GP) formalism. By comparing
the TF predictions with the GP solutions, in our investigation with two
different quartic-quadratic trap geometries, of particular relevance is to
observe that the TF approach is not only very useful to display the averaged
density distribution, but also quite realistic in establishing the critical
rotational conditions for the breakup occurrence and possible giant-vortex
formation. It provides almost exact results to define the contour of the
condensate distribution, even for high rotating system, after the system split
in two (still confined) clouds. The applicability of the Feynman rule to the
vortex distribution (full-numerical GP solutions) is also being confirmed for
these non-homogeneous asymmetric trap configurations. This study is expected
to be relevant for manipulating the rotation and trap parameters in addition
to Feshbach resonance techniques. It can also be helpful to define initial
conditions for any further studies on dynamical evolution of vortex pattern
distributions.

##  I Introduction

The fragmentation phenomenon of objects due to their fast rotation is
ubiquitous in physical systems, ranging from classical astronomical objects,
such as galaxies and supermassive rotating stars Maeder2009 ; 2017Uchida ;
2017Chiappini , to the microscopic world of quantum phenomena, in nuclear
physics 1979Bohr ; 1998Clark . In atomic-molecular physics, this possibility
could also be visualized by considering the experimental realization in Ref.
2020Guo  of superfluid flow of quantum gas rotating in an anharmonic
potential. Theoretically, the occurrence of breakup in a rotating dipolar
condensate was suggested in Ref. Kumar2016 , when studying three-dimensional
(3D) vortex structures in a rotating dipolar Bose-Einstein condensate (BEC).
For coupled condensed systems, the breakup can be easily understood as related
to the miscibility properties of the mixture, such that the system split in
two spatially separated condensates when going to the immiscible regime, in
which the rotation can help to reduce the confinement of the two species
together 1998Barankov ; 2019Wilson ; 2012Web ; 2014Pattinson ; 2014Filatrella
; 2017Kumar ; 2019Kumar . With single species atomic condensates, there are
plenty of studies with different kind of rotating trap confinements, which
allows the generation of multiple quantized vortices, when considering enough
strong power-law and anharmonic potentials Fetter2001 ; 2002Lundh ;
2003KavoulakisNJP ; Fetter2005 ; 2009Fetter . Among them, of particular
interest are the investigations with the more experimentally feasible quartic-
quadratic trap combinations 2002Lundh ; 2003KavoulakisNJP , such that one can
identify formations of circular super-flow and stable giant vortex, when
considering fast and overcritical rotations 2001Recati ; 2002Kasamatsu ;
2003Engels ; 2003Fischer ; 2004Bretin ; 2004Aftalion ; 2004Simula ; 2004Baym ;
2011Correggi ; 2013Correggi ; 2019Adhikari . However, in these cases, without
interspecies interaction to play, we are not aware of specific approaches
concerning the possibilities to fragment a condensate by introducing high
rotation frequencies in a controlled way, with the fragments still kept
confined. In this direction, we can foresee the experimental possibilities,
which are accessible, as one can follow from the recent studies performed with
high rotation superfluids and condensates, generating long-lived dynamical
rings 2020Guo ; 2019Pandey . In order to generate long-lived stable vortices,
the experiments are usually performed with repulsive two-body interactions, as
being already well-known that vortices are generally unstable in rotating BEC
with attractive interactions 2004Lundh ; 2006Carr ; 2008Sakaguchi ;
2015Sudharsan , which are closely related to the condensate requirements for
stability 1999Dalfovo ; 2009Fetter .

In the present work, we consider a BEC with repulsive two-body interactions,
as our purpose is to study more closely the possibilities to break a stable
condensate by increasing the rotation. For that an originally simple
pancakelike two-dimensional (2D) trap confinement is assumed, which is
obtained by a much stronger confinement in one of the directions, which can be
easily implemented experimentally. Our aim with such study is to fill a gap in
the basic literature concerning the properties of rotating BECs 1999Dalfovo ;
2006Aftalion ; 2008Pethick ; 2009Fetter ; 2016Pitaevskii , which can be
relevant to control and calibrate possible experimental realizations when
considering more involved studies on properties of ultra-cold condensed
systems. By introducing such a breakup mechanism in an experimental apparatus,
with the fragments of the condensate still being kept confined, one could have
access to another mechanism to control the condensed cloud, in addition to the
Feshbach resonance techniques 1998inouye ; 2010chin . In view of the above
motivations for considering a condensate breakup, we introduce a theoretical
demonstration that such procedure can be easily implemented in a single non-
dipolar condensate, induced by the rotation in an asymmetric trap. For that,
we rely in the Thomas-Fermi (TF) approximation, for a 2D confined system,
supported by the corresponding full-numerical results obtained via the usual
Gross-Pitaevskii (GP) formalism. The numerical results are performed by
imaginary time propagation approach 2006-Brtka  and, complementary, refined by
fixed point methods 2009-Yang .

In order to simplify our analysis and the system be experimentally feasible,
we are assuming a strongly confined quasi-2D trap, considering two distinct
quartic-quadratic non-symmetric shapes in the pancakelike perpendicular plane.
Both trap non-symmetric geometries are equally simple to be implemented in
actual cold-atom experiments, which are slightly different from the more
symmetric quartic-quadratic confinements being so far investigated
2003KavoulakisNJP ; 2009Fetter ; Fetter2005 ; 2005Cozzini . In both the cases
that we are studying, a general double-well structure is assumed for the
asymmetric confinement, such that just by manipulating the rotation frequency
parameter $\Omega$, one can deform and split the condensed cloud in two
fragments, which can still be kept strongly confined within fixed two-well
positions, which are affected by the $\Omega$ value. Of particular relevance
for possible experimental proposals, as it will be further discussed, is the
outcome of our study on the reliability of the TF approximation, which is
quite robust and useful to define the limiting rotational conditions, even for
the asymmetric trap structures under analysis, before and after the condensate
splits in two fragments.

By full-numerical solution of the GP formalism, we are also investigating the
vortex structure patterns, which occur in the two kind of asymmetric quartic-
quadratic trap configurations being considered. The applicability of the well-
known Feynman rule (FR) 1955-Feynman  in these quite asymmetric geometries is
well established, by verifying the vortex density numbers inside the confined
condensed cloud, even when applying the rule to each fragment of the
condensate, once it is appropriately defined regions with enough vorticity.
The role played by hidden vortices (located in low-density regions) is also
being characterized in our analysis.

In the next Sect. II, we present some details on the BEC formalism, with
quartic-quadratic confinement, in which we are paying particular attention to
the related TF description. All the numerical results for the GP formalism, as
well as for the TF approach, are presented and discussed along this section,
following the formalism. In the Sect. III, we present our analysis on the
vortex-pattern distributions verified in the densities, in which we discuss
the validity of the well-known Feynman rule for asymmetric shaped BEC
confinement. Finally, in Sec. IV we resume the work with our conclusions and
perspectives for further related studies.

##  II Quartic-quadratic trapped BEC

In the next, we present the usual basic quasi-2D GP formalism, obtained from a
3D reduction to a strongly pancakelike confined system in the
$x_{3}-$direction. In the formalism two possible trap interactions are being
considered in the $(x_{1},x_{2})-$plane. The breakup conditions are
investigated in rotating condensate systems, in which we have verified that
the usual TF approximation can successfully describe the density distribution,
as well as it is quite helpful in defining the contours of the condensate,
before and after the breakup. Therefore, this approximation will be used in
order to obtain the necessary rotating frequency in which the fragmentation of
the condensate occurs in two parts.

For a condensed system of particles having mass $m$, confined by a trap
potential $V_{3D}(x_{1},x_{2},x_{3})\equiv
V_{2D}(x_{1},x_{2})+m\omega_{3}^{2}x_{3}^{2}/2$, where $\omega_{3}$ is the
trap frequency along the $x_{3}-$direction, within a rotating frame with
angular frequency $\overline{\Omega}$ about $x_{3}$, where the the angular
momentum operator is $L_{3}={\rm i}\hbar\left(x_{2}\partial/\partial
x_{1}-x_{1}\partial/\partial x_{2}\right)$, the 2D reduction of the original
3D stationary solution of the modified Gross-Pitaevskii (GP) formalism can be
written, for the chemical potential $\mu_{2}\equiv\mu_{3}-\hbar\omega_{3}/2$,
as

$\displaystyle\mu_{2}\psi_{2D}$ $\displaystyle=$
$\displaystyle\left[-\frac{\hbar^{2}}{2m}\nabla^{2}_{2D}+V_{2D}+\overline{\Omega}L_{3}+g_{2D}\left|\psi_{2D}\right|^{2}\right]\psi_{2D},$
(1)

where $\nabla_{2D}^{2}\equiv\frac{\partial^{2}}{\partial
x_{1}^{2}}+\frac{\partial^{2}}{\partial x_{2}^{2}}$ and
$\psi_{2D}\equiv\psi_{2D}(x_{1},x_{2})$ is the 2D wave function normalized to
the number of atoms $N$ (the corresponding density given by
$n_{2D}\equiv|\psi_{2D}|^{2}$), with $V_{2D}\equiv V_{2D}(x_{1},x_{2})$ being
the 2D trap potential. The strength of the nonlinear interaction, $g_{2D}$,
related to the corresponding $g_{3D}$ due to the two-body $s-$wave interaction
$a_{s}$, is given by

$g_{2D}\equiv
g_{3D}\sqrt{\frac{m\omega_{3}}{2\pi\hbar}}\equiv\left(\frac{4\pi\hbar^{2}a_{s}}{m}\right)\sqrt{\frac{m\omega_{3}}{2\pi\hbar}}.$
(2)

The total energy, related to (1) can be written in terms of the 2D vector
${\bf r}$ and its conjugate momentum $\mathbf{p}$, as

$\displaystyle E[\psi_{2D}]$ $\displaystyle=$
$\displaystyle\int\\!\\!\mathrm{d}^{2}\mathbf{r}\left\\{\frac{1}{2m}\psi_{2D}^{*}(\mathbf{p}-\mathbf{\overline{\Omega}}\times
m\mathbf{r})^{2}\psi_{2D}\right.$ $\displaystyle+$
$\displaystyle\left.\left[V_{2D}-\frac{m\overline{\Omega}^{2}r^{2}}{2}\right]|\psi_{2D}|^{2}+\frac{g_{2D}}{2}|\psi_{2D}|^{4}\right\\}.$

With the energy and length units given, respectively, by $\hbar\omega$ and
$l_{\omega}\equiv\sqrt{\frac{\hbar}{m\omega}}$, we have
$({x},{y})=(x_{1}/l_{\omega},x_{2}/l_{\omega})$, the GP formalism can be
written in dimensionless units. For that, we redefine the respective
quantities, appropriately, as $\mu\equiv\mu_{2}/(\hbar\omega)$, $V\equiv
V(x,y)\equiv V_{2D}/(\hbar\omega)$, $\psi\equiv\psi(x,y)\equiv
l_{\omega}\psi_{2D}/\sqrt{N}$, $\Omega\equiv\overline{\Omega}/\omega$ and
$L_{z}=L_{3}/\hbar$. Therefore, from Eqs. (1)-(II), we obtain the
corresponding dimensionless equations, as

$\displaystyle\mu\psi=\left[-\frac{1}{2}\left(\frac{\partial^{2}}{\partial{{x}^{2}}}+\frac{\partial^{2}}{\partial{{y}^{2}}}\right)+V+\Omega
L_{z}+g\left|\psi\right|^{2}\right]\psi,$ (4)

where $\psi$ is normalized to one. In the above, within our units, $\omega$ is
an arbitrary frequency parameter, which is assumed to be perpendicular to the
$z-$direction, with $\lambda\equiv\omega_{3}/\omega$ being the trap aspect
ratio. In this case, the dimensionless nonlinear strength carries the
information about the number of atoms $N$, trap aspect ratio
$\lambda\equiv{\omega_{3}/\omega}$ and two-body $s-$wave scattering length
$a_{s}$:

$g=N\sqrt{8\pi\lambda}\left(\frac{a_{s}}{l_{\omega}}\right).$ (5)

We assume repulsive two-body interaction, with $a_{s}>0$ throughout this
study, due to stability requirements of the condensate cloud, such that $g$ is
always positive. The dimensionless energy, corresponding to (II) is given by

$\displaystyle{\cal E}[\psi]$ $\displaystyle=$
$\displaystyle\int\\!\\!dx\,dy\left\\{\frac{1}{2}\left[\left|{\rm
i}\frac{\partial\psi}{\partial x}-\Omega y\psi\right|^{2}+\left|{\rm
i}\frac{\partial\psi}{\partial y}+\Omega x\psi\right|^{2}\right]\right.$ (6)
$\displaystyle+$
$\displaystyle\left.\left[V(x,y)-\frac{{\Omega}^{2}(x^{2}+y^{2})}{2}\right]|\psi|^{2}+\frac{g}{2}|\psi|^{4}\right\\},$

from where we can observe that the trap potential is going to be effectively
reduced as we increase the rotation frequency. About our full-dimension trap
interaction, $V_{3D}(x_{1},x_{2},x_{3})$, as usual we are assuming that the
trap aspect ratio $\lambda=\omega_{3}/\omega$ is large enough in order to
validate our quasi-2D geometry approach, in consonance with possible
experimental realizations. For that, one can also realize that the number of
particles $N$ together with $\lambda$ and the two-body $s-$wave contact
interaction $a_{s}$ can be conveniently adapted according to the system being
investigated, by adjusting the strength $g$ [Eq. (5)].

Within our aim to study the possibility for the rotation frequency to distort
up to the breakup of the condensate, let us assume non-symmetric 2D trap
potential in the $(x,y)$ plane, given by a quadratic-quartic interaction. For
that, let us fix an harmonic oscillator trap in the $y-$direction, implying
$\omega_{y}=\omega$ in our dimensionless units. For the quartic term, which
impose a more strong confinement, we are going to consider two possibilities:
The first, just $x-$dependent, with the second more symmetric, exchanging the
$x^{4}$ by $r^{4}$. We can write both trap models as given by

$\displaystyle V_{\sigma}({x},{y})$ $\displaystyle=$
$\displaystyle\frac{y^{2}}{2}+\frac{b^{2}}{4}\left[\left(x^{2}-\frac{a^{2}}{b^{2}}\right)^{2}+\sigma\left(r^{4}-x^{4}\right)\right],$
(7)

where $\sigma=0$ ($x-$quartic) and $\sigma=1$ ($r-$quartic) define the two
kind of traps. In both cases, $b^{2}/4$ is the strength of the quartic term.
The minima for $V_{\sigma}(x,y)$ are given by $x\equiv\overline{x}=\pm a/b$
and $y=0$, along $x$ and $y$ directions, respectively.

In the following, we provide a detailed Thomas-Fermi approximation analysis
for the density distributions, which will show up to be quite realistic as
compared with the numerical results, when used to define the threshold
conditions.

###  II.1 Thomas-Fermi approximation.

From Eq. (II), we can derive the TF approximation, when the atomic repulsive
interactions are dominating, by expanding the condensate to a mean radius that
exceeds the mean oscillator length $l_{\omega}$. In this case, the expansion
reduces the radial gradient of the density, with the kinetic energy becoming
negligible relative to the trap and interactions, such that

$\displaystyle{\cal E}_{TF}[\psi]$ $\displaystyle\approx$
$\displaystyle\int\\!\\!dx\,dy\left\\{\left[V_{\sigma}(x,y)-\frac{{\Omega}^{2}r^{2}}{2}\right]|\psi|^{2}+\frac{g}{2}|\psi|^{4}\right\\}.$
(8)

By minimizing this energy with respect to $|\psi|^{2}$ at fixed normalization,
we obtain the TF density, which is a function of $\Omega$ through
$\mu_{TF}\equiv\mu_{TF}(\Omega)$. It is given by

$\displaystyle n(x,y;\Omega)$ $\displaystyle\equiv$
$\displaystyle|\psi|^{2}=\frac{\mu_{TF}}{g}-\frac{1}{g}\left[V_{\sigma}(x,y)-\frac{{\Omega}^{2}r^{2}}{2}\right].$
(9)

We emphasize here that besides the double-well structure provided by
$V_{\sigma}$, the centrifugal part effectively enlarge the distance between
the minima, as can be seen from the term inside square brackets in the above
Eq. (9), where for a fixed $\Omega$, the density $n(x,y;\Omega)$ have two
maxima located at
$(x,y)=\left(\displaystyle\pm\sqrt{\frac{\Omega^{2}+a^{2}}{b^{2}}},0\right)$.
In this case, the density at the center of the cloud, given by

$n_{0}\equiv
n(0,0;\Omega)\equiv\frac{\mu_{TF}-V_{\sigma}(0,0)}{g}=\frac{1}{g}\left[\mu_{TF}-\displaystyle\frac{a^{4}}{4b^{2}}\right],$
(10)

shall decrease by increasing the rotation, due to centrifugal forces, with the
distances between the two density maxima increasing. Within this dynamical
process, $n_{0}$ can be zero, with the creation of a giant vortex at the
center, as well as induce a breakup of the condensate in two clouds, which
will depend on the specific trap potential configuration. So, the conditions
to occur the breakup, or to generate the giant vortex, are established when
the effective chemical potential $\mu_{TF}-a^{4}/(4b^{2})$ becomes zero
($gn_{0}=0$). However, by increasing even more the rotation $\Omega$, the
effective chemical potential becomes negative, with $gn_{0}$ defined by Eq.
(10) representing this effective chemical potential, and not just the density
at the center.

For the $V_{\sigma}$ traps defined in Eq. (7), the centrifugal contribution
can either generate a giant vortex or split the cloud for $\sigma=1$
($r-$quartic) case, in contrast to the $\sigma=0$ ($x-$quartic) case, where
the cloud can split as well, but no giant-vortex phase is generated (as it
will be shown). Therefore, a more relevant criterion is to analyze the contour
of the Thomas-Fermi density in polar coordinates, as a function $r(\theta)$,
since the breakup can easily be identified with the absence of real solutions
at $\theta=\pm\pi/2$ (in view of the maxima distribution along the
$x-$direction). In order to support this TF criterion to define the critical
rotation conditions for the condensate breakup (both kind of potentials), as
well as for the giant-vortex formation (case of $r-$quartic), we rely on the
close agreement between the TF approximation with the numerical GP solution
when defining the condensate contour. The two trap models defined in Eq. (7),
for $\sigma=0$ ($x-$quartic case) and $\sigma=1$ ($r-$quartic case), will be
treated separately in the next subsections.

###  II.2 Case of $x-$quartic confinement

The contour of the cloud condensate, obtained by solving $n(x,y;\Omega)=0$,
before and after the breakup, is similar to what is known in mathematics as
the Cassini Ovals cassini . Just at the breakup limit, the corresponding oval
contour curve can also be identified with the Bernoulli Lemniscate. In polar
coordinates, with
$(x,y)\equiv[x(\theta),y(\theta)]\equiv[r(\theta)\cos{\theta},r(\theta)\sin{\theta}]$,
by applying the contour equation $n(x_{c},y_{c};\Omega)=0$, we obtain the
angular dependence of $r_{c}\equiv r_{c}(\theta)$, which is given by

$4gn_{0}=b^{2}r_{c}^{4}\cos^{4}\theta-2r_{c}^{2}[(1+a^{2})\cos^{2}\theta-(1-\Omega^{2})].$
(11)

The limit $\theta\to\pi/2$, given by (11),
$r_{c}^{2}(\pi/2)=2gn_{0}/(1-\Omega^{2})$, reduces to zero in the breakup,
given by Eq. (10),

$n_{0}(\Omega_{s})=0\rightarrow\mu_{TF}(\Omega_{s})=\frac{a^{4}}{4b^{2}},$
(12)

where $\Omega_{s}$ is the critical rotation frequency to break the condensate
in two separated clouds. The above relation expresses the fact that the
chemical potential, which is a general function of the rotation parameter
$\Omega$ given by Eq. (9), at the breakup (when the density reduces to zero)
is given by $\Omega_{s}$. Therefore, the dependence of the breakup condition
on the condensate density parameter $g$ can only become explicit once a
functional relation becomes established between $\Omega_{s}$ and $g$, as it
will be shown.

The more general solution of Eq. (11), for the contour $r_{c}$, is given by

$\displaystyle b^{2}\,r_{c,\pm}^{2}\cos^{2}\theta$ $\displaystyle=$
$\displaystyle 1+a^{2}-(1-\Omega^{2})\sec^{2}\theta$ $\displaystyle\pm$
$\displaystyle\sqrt{\left[{1+a^{2}-(1-\Omega^{2})\sec^{2}\theta}\right]^{2}+{4b^{2}gn_{0}}},$

where $r_{c,+}(\theta)$ and $r_{c,-}(\theta)$ are the two possible solutions.
For $n_{0}>0$, only one solution is possible, $r^{2}_{c,+}(\theta)$, giving us
the full TF contour for $0\leq\theta<2\pi$. For $n_{0}<0$ (after the breakup),
two possible solutions emerge, with $r_{c,-}(\theta)$ corresponding to inner
and $r_{c,+}(\theta)$ to outer contour points.

Therefore, aside of $\Omega_{s}$ being the critical rotation for breakup, it
is also the threshold minus sign solution for $\Omega>\Omega_{s}$.
Specifically at the critical rotation frequency $\Omega_{s}$, with $\mu_{TF}$
given by Eq. (12), we have the contour delimited by the $+$ sign solution as

$\displaystyle r_{c,+}^{2}(\theta)$ $\displaystyle=$ $\displaystyle
2\frac{\sec^{4}{\theta}}{b^{2}}\left[a^{2}+\Omega_{s}^{2}-(1+a^{2})\sin^{2}\theta\right]$
(14) $\displaystyle\times$
$\displaystyle\Theta\left(\frac{a^{2}+\Omega_{s}^{2}}{1+a^{2}}-\sin^{2}\theta\right),$

where $\Theta(\xi)$ is the Heaviside step function (=1 for $\xi\geq 0$, =0 for
$\xi<0$). The above implies that, at the critical breakup frequency, the
radial contour is limited within a fixed interval. By applying the
normalization in Eq. (9), together with (14), after the radial integration,
for the given trap parameters $a$ and $b$, we obtain the relation between the
critical splitting rotation $\Omega_{s}$ and the density nonlinear parameter
$g$, which is given by

$\displaystyle{gb^{4}}$ $\displaystyle=$
$\displaystyle\int^{\theta_{m}}_{0}d\theta\cos^{4}\theta\,r_{c}^{6}(\theta)=\frac{32}{105}\sqrt{\frac{(a^{2}+\Omega_{s}^{2})^{7}}{1-\Omega_{s}^{2}}},$
(15)

where
$\theta_{m}\equiv\arcsin\left(\sqrt{(a^{2}+\Omega_{s}^{2})/(1+a^{2})}\right)\in[0,\pi/2]$
is fixed by the Heaviside step function in Eq. (14). The parameter $g$, as
defined in (5), is determined by the number of atoms $N$, the two-body
scattering length and the pancake 2D asymmetry $\lambda$. Therefore, given $g$
and the trap parameters $a$ and $b$, we obtain analytically the critical
frequency $\Omega_{s}$ for the condensate to split. Within our purpose to
verify the effect of a rotation in breaking a condensate, we are obviously
considering that the confining asymmetric double-well shaped trap (7) have the
parameter-ratio $a/b$ such that the two minima position are close enough in
order to have a continuous distribution of the condensate when $\Omega=0$.

Figure 1:  Rotation frequency $\Omega$ as function of $gb^{4}$, for the trap
defined by Eq. (7), with $\sigma=0$ and the nonlinear strength $g$ given by
Eq. (5). The curves, for $a=0$ (dashed) and $a=0.5$ (solid), refer to the TF
lower limiting values, $\Omega_{s}$, to split the condensate. The circles
(triangles) are indicating the positions of density-plot illustrations for
$a=0.5$, shown in 4 panels of Fig. 2 (3 panels of Fig. 3), with $gb^{4}=0.05$
($gb^{4}=0.25$), where $\Omega_{s}=$ 0.5616 (0.7640). With units defined in
the text, all parameters are dimensionless.

Figure 2:  The densities $|\psi(x,y)|^{2}$, obtained from the numerical GP
solutions, are shown in the upper panels, for $\Omega=0.3$ (a), $0.4$ (b),
$0.5$ (c) and $0.6$ (d), with the TF approximation boundaries represented by
white-dashed lines. In panels (e) and (f), corresponding to $y=0$ densities
(c) and (d), respectively, we compare TF (dashed curves) with GP (solid
curves) results. The other parameters are $g=500$, $b=0.1$ ($gb^{4}=0.05$) and
$a=0.5$. The effective potential for the panel (d) is illustrated by the
surface plot [with $z\equiv V_{0}(x,y)-0.18r^{2}$]. With the defined units,
all quantities are dimensionless.

In Fig. 1 we are resuming the TF results for $\Omega_{s}$ as a function of
$gb^{4}$, by considering two indicative fixed values of the parameter $a$ (0
and 0.5). Few sample results for the densities, showing how the rotation is
affecting the dynamics of the condensate, are also illustrated in Figs. 2
(with $g=500$) and 3 (with $g=2500$), considering numerical GP results. In
both the cases, we have $a=0.5$ and $b=0.1$, such that the minima of the
quartic interaction are at $x=\pm 5$. The corresponding parameter positions of
these results ($gb^{4}=$0.05 and 0.25, respectively) are also indicated inside
the panel of Fig. 1. The corresponding breakup rotational parameter are
$\Omega_{s}=$ 0.5616 and 0.7640.

In Fig. 2, we show the numerical GP simulations together with the TF
boundaries for the densities, as considering four panels with different
arbitrary values of $\Omega$, from 0.3 till 0.6. This figure illustrates how
the breakup occurs in a condensate trapped by $V_{0}(x,y)$ given by Eq. (7),
induced by the rotation frequency, which is reducing the strength and shifting
the position of the minima of the corresponding effective double-well
potential. For a better quantitative analysis, we have also included the
panels (e) and (f) corresponding to $|\psi(x,0)|^{2}$ of panels (c) and (d).
The effective potential related to panel (d) is also being represented in the
bottom part of the figure by a surface 3D plot. In particular, we consider
remarkable the agreement between the numerical GP computed results with the
corresponding TF approximation estimates, when applied to the contour of the
condensate (inside which we can find most of the density distributions). Such
results obtained by using the TF approximation, as compared with the numerical
computed ones, are indicating that it can be quite useful in order to verify
the exact rotation frequency, $\Omega_{s}$, at which one should expect the
condensate being fragmented in two pieces, once considered fixed the potential
parameters. The results presented in the panels (e) and (f) are particularly
significant to support the TF approximation in order to define the breakup
limiting conditions. In spite of the obvious limitations in describing the
vortex pattern structure, the TF approximation provides an averaged good
agreement with the numerical GP results for the internal averaged density
distributions, remarkably reproducing the corresponding contour borders.

In Table 1, the quantitative deviations between the TF and the GP numerical
results are verified for the chemical potentials corresponding to the four
specific cases shown in Fig. 2. Among the results, we also add the ones
corresponding to the total energy, together with the cases with $\Omega=0$.

Table 1:  The chemical potentials (GP numerical simulations and TF approximation), for $\Omega=0$ and the frequencies considered in Fig. 2 are given in this table, together with the corresponding GP Energy values. The units are $\hbar\omega$ for the energy and chemical potentials, with $\omega$ for the frequency $\Omega$. $\Omega$ | $\mu_{TF}$ | $\mu$ | Energy
---|---|---|---
0 | 7.332 | 7.395 | 4.849
0.3 | 5.994 | 6.469 | 3.962
0.4 | 4.813 | 5.407 | 2.955
0.5 | 3.061 | 3.801 | 1.271
0.6 | 0.343 | 1.215 | -1.490

Figure 3:  Density distributions (left) are shown with the corresponding phase
diagrams (right), for the nonlinear interaction $g=$ 2500 and rotation
frequencies $\Omega=$ 0.7 [panel (a)] and 0.77 [panel (b)], whereas the
breakup critical value is at $\Omega_{s}=$ 0.764. The trap parameters are the
same as in Fig. 2. The TF boundaries are given by dashed lines, with the upper
color bars indicating the respective density and phase levels. With the
defined units, all quantities are dimensionless.

From the GP results shown in Fig. 2, we observe the crescent number of
vortices formed inside the condensate as we increase the rotation. The linear
increasing behavior of the vortices number with the frequency is clearly being
broken near the breakup limit, giving us some indication of the occurrence of
possible hidden vortices at the low-density region. This behavior will be
further discussed in the next section, within our analysis of vortex pattern
distributions.

In order to improve our illustrations of the $x-$quartic confining potential,
with densities in which we can verify more clearly how the vortex pattern
distributions occur before and after the breakup, we include the Fig. 3, where
we consider two specific cases with frequency rotations $\Omega=$ 0.70 (a) and
0.77 (b), with the nonlinear parameter, $g=2500$, which is much larger than in
the case of Fig. 2. However, in both the cases, we are keeping the same trap
parameters, $a=0.5$ and $b=0.1$, with the potential minima at $x=\pm 5$. These
plots, in which a richer pattern structure can be observed, were obtained by
the GP numerical simulations for the densities (left panels), with the
corresponding phase diagrams (right panels). The boundary definitions,
indicated inside the panels, were obtained by using the TF approximation. As
shown in the panel (b), with rotation frequency $\Omega=0.77$, the condensate
is already broken in two fragments. Due to the rotation $\Omega=0.77$, the
positions of the minima of the effective potential move to $x=\pm 9.18$, with
both fragments presenting about the same kind of triangle vortex pattern
distributions. As for $\Omega=0.70$, panel (a), we observe about the same
vortex pattern, with the minima of the effective potential being at $x=\pm 9$.
In both phase diagrams we observe the existence of hidden vortices which are
located in the low-density regions, along the $y-$ axis, close to $x=0$. More
specifically, these hidden vortices are located outside the region defined by
the TF approximation, as they are manifested only by solving the full GP
formalism.

Figure 4: The dependence of the effective chemical potentials [GP and TF
results], given by $\mu-{a^{4}}/{(4b^{2})}$, on the rotation frequency
$\Omega$ for the trap parameters $a=0.5$, $b=0.1$. The TF (GP) results are
indicated, respectively, by the circles (squares), being fitted by solid
(dashed) lines. At $\Omega=0$, $\mu_{0}=$7.332 (TF) and 7.395 (GP). For zero
density, when $\mu=(5/4)^{2}$, we obtain $\Omega_{s}=0.5616$, within the TF
approximation; and 0.590 within the GP result. With the given units, all
quantities are dimensionless.  Figure 5: For the $x-$quartic trap potential
[Eq.(7) with $\sigma=0$], we show by a color-level map (indicating the
necessary split frequency $\Omega_{s}$) the corresponding parametric region
for which a rotation is required for breaking the cloud. The upper white
region corresponds to $\Omega_{s}>0.95$ at which the condensate is likely to
fly apart in the $y-$direction before the split occurrence. Below the white-
dashed line ($\Omega_{s}=0$), the cloud is already spatially separated by the
two-well shaped trap. With given units, all quantities are dimensionless.

By considering the same potential parameters used in Table 1 and in Fig. 2,
with $a=0.5$ and $b=0.1$, we also present in Fig. 4 the corresponding behavior
of the effective chemical potential, given by $\mu-{a^{4}}/{(4b^{2})}$, in
terms of the rotation frequency $\Omega$, which is zero at the condensate
breakup, when $\Omega=\Omega_{s}=0.5616$. For that we consider the TF
approximation results together with numerical solutions obtained from the GP
equation. The results can be well fitted by approximate analytical polynomial
expressions. From these analytical results, we observe the quite small
deviation of the TF critical breakup rotation limits,
$\Delta\Omega_{s}=0.028$. We also observe that the deviation of the TF
chemical potential $\mu_{TF}$, from the GP values $\mu$ becomes larger for
fast rotations, which one could expect due to the fact that, by increasing the
rotation, more vortices are emerging in the condensate, with the effective
confining interaction being diminished in comparison with the kinetic energy
term. This is evidenced by the density plots shown in Fig. 2.

In Fig. 5, a diagrammatic picture shows the relation between the trap
parameters $gb^{4}$ and $a$ for the $x-$quartic potential [Eq. (7) with
$\sigma=0$], indicating the necessary rotation $\Omega_{s}$ to split the
condensate. $\Omega_{s}$ is represented by the color density, considering the
upper-white region $\Omega_{s}>0.95$ the region at which the condensate likely
flies apart before the breakup. The parameter region at which the condensate
is already broken by the potential parameters is shown below the dashed line
(where $\Omega_{s}=0$).

###  II.3 Case of $r-$quartic confinement

By considering the confining potential given by Eq. (7) with $\sigma=1$, we
have a trap with quadratic+quartic radial symmetry in the $(x,y)$ plane, with
the symmetry broken in the $x-$direction, which can be written as

$\displaystyle V_{1}(x,y)$ $\displaystyle=$
$\displaystyle\frac{r^{2}}{2}+\frac{b^{2}}{4}r^{4}-(1+a^{2})\frac{x^{2}}{2}+\frac{a^{4}}{4b^{2}},$
(16)

Different from the previous $x-$quartic trap, in this case we have also a
quartic component in the $y-$direction. The $|x|$ and $|y|$ extreme limits for
non-zero density, obtained by using the Eq. (9) are defined by

$\displaystyle n(x,y)$ $\displaystyle\equiv$
$\displaystyle|\psi|^{2}=\frac{\mu_{TF}}{g}-\frac{1}{g}\left[V(x,y)-\frac{{\Omega}^{2}r^{2}}{2}\right]$
$\displaystyle=$ $\displaystyle
n_{0}-\frac{1}{2g}\left({{y}^{2}-{a^{2}x^{2}}+\frac{b^{2}}{2}r^{4}-\Omega^{2}r^{2}}\right),,$

where $gn_{0}\equiv\left[\mu_{TF}-a^{4}/(4b^{2})\right]$ as also given by Eq.
(10). Analogously as in the $x-$quartic case, we could follow by using polar
coordinates. However, before that, it may be useful to verify the extremes of
non-zero density $n\equiv n(x,y)$. So, by deriving in the $x-$coordinate, we
can identify a minimum at $x=x_{0}=0$ with two maxima at $x=x_{\pm}$ of
(II.3), given by

$\displaystyle x_{\pm}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\frac{\Omega^{2}+a^{2}}{b^{2}}-y^{2}},\;\;x_{0}=0.$
(18)

Similarly, by deriving in the $y-$coordinate (given $x$), we can identify the
existence of two maxima at $y=y_{\pm}$ with a minimum at $y=y_{0}=0$, only
when $\Omega^{2}>1+b^{2}x^{2}$,

$\displaystyle y_{\pm}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\frac{\Omega^{2}-1}{b^{2}}-x^{2}},\;\;y_{0}=0.$ (19)

Combining both results, we verify that for $\Omega^{2}>1$ we can have extrema
(two maxima and one minimum) only in the $x-$direction, located at
$(x_{\pm},0)$ ($\theta=0$), as the maxima obtained by (19) are located outside
the positions where the previously broken condensate parts are located, as we
can see by replacing $x_{\pm}^{2}$ in (19).

Now, let us consider the density contour given by Eq. (II.3) in polar
coordinates, $(x,y)\equiv[x(\theta),y(\theta)]$
$\equiv[r(\theta)\cos\theta,r(\theta)\sin(\theta)]$, such that $r_{c}\equiv
r_{c}(\theta)$ is given by

$4gn_{0}=b^{2}r_{c}^{4}-2r_{c}^{2}\Big{[}(1+a^{2})\cos^{2}{\theta}-(1-\Omega^{2})\Big{]},$
(20) $\displaystyle b^{2}r_{c}^{2}$ $\displaystyle=$
$\displaystyle{(1+a^{2})\cos^{2}\theta-(1-\Omega^{2})}$ $\displaystyle\pm$
$\displaystyle\sqrt{\left[{(1+a^{2})\cos^{2}\theta-(1-\Omega^{2})}\right]^{2}+{4b^{2}gn_{0}}}.$

We should emphasize that the necessary and sufficient condition for the
breakup is that no real solution occurs for $r_{c}^{2}(\pm\pi/2)\equiv
y_{c}^{2}$ (implying density zero along the $y$ axis, not only at $y=y_{c}$).
In the $x-$quartic case ($\sigma=0$) the $n_{0}=0$ condition is enough for the
breakup, with the frequency being limited to $\Omega\leq 1$, as shown.
Nevertheless, differently from $x-$quartic case, as we take $\theta=\pm\pi/2$
in Eq. (20) we still end up with a bi-quadratic equation for the radius, with
the Eq. (II.3) simplified as

$b^{2}y_{c}^{2}={(\Omega^{2}-1)}\pm\sqrt{{(\Omega^{2}-1)^{2}}+{4b^{2}\left[\mu_{TF}-\left(\frac{a^{4}}{4b^{2}}\right)\right]}}.$
(22)

The existence of solutions for this equations implies that, at
$\theta=\pm\pi/2$, the contour is connected. So, next we consider all the
possibilities in (22), where $gn_{0}=\mu_{TF}-a^{4}/(4b^{2})$: If $gn_{0}>0$
only the $+$ sign solution can exist, corresponding to a connected contour in
the TF description. However, if $gn_{0}\leq 0$, the existence of real
solutions (characterizing contour connected) depends on $\Omega$: For
$\Omega<1$, neither the $+$ nor the $-$ solutions are possible, characterizing
the breakup of the cloud. However, for $\Omega>1$, both solutions are
possible, if $0>4b^{2}gn_{0}>-(\Omega^{2}-1)^{2}$, corresponding to the upper
and lower solutions of the connected contour, which is identifying the giant-
vortex formation at the center. The other possibility, when
$4b^{2}gn_{0}<-(\Omega^{2}-1)^{2}$, there is no real solution for Eq. (22),
which will corresponds to the breakup. Summarizing:

$\begin{array}[]{c|l|l}\Omega>1&\displaystyle{-(\Omega^{2}-1)^{2}}<4b^{2}gn_{0}\leq
0&\textrm{Giant vortex}\\\ &\displaystyle{-(\Omega^{2}-1)^{2}}\geq
4b^{2}gn_{0}&\textrm{Breakup}\\\ \hline\cr\Omega\leq 1&gn_{0}(\Omega)\leq
0&\textrm{Breakup}\end{array}$ (23)

The critical values for the rotation parameter $\Omega$ are corresponding to
the transitions among the phases, which are obtained when the equality holds
in the above relations, being defined as $\Omega_{G}$ for the Giant-vortex
formation (which can happen only for $\Omega>1$), and $\Omega_{s}$ the
critical frequency to split the cloud.

Differently from the first trap studied, in the one defined by Eq. (16), the
breakup regime can be achieved for larger values of $b$ since the condensate
will no longer fly apart due to the presence of quartic confinement also in
the $y$-direction.

In the above, we have considered the possible breakup of the condensate (as
well as the giant-vortex formation in the origin) along the $y-$axis
(considering $\theta=\pm\pi/2$, or $x=0$), when by increasing the rotation
(with $\Omega>1$) the quadratic confinement in the $y-$direction is reduced
till the condensate split in two fragments. However, a similar situation could
occur in the $x-$direction, if the asymmetry of the trap potential which is
holding the two fragments can allow minima for the density at $y=0$
$(\theta=0)$ for fixed values of $\pm x$, in Eq. (II.3). This should be given
by the following second-order equation:

$b^{2}x_{c}^{2}={(\Omega^{2}+a^{2})}\pm\sqrt{{(\Omega^{2}+a^{2})^{2}}+{4b^{2}gn_{0}(\Omega)}}.$
(24)

For larger values of $\Omega^{2}$ ($\gg a^{2}$), the Eqs. (22) and (24)
provide the same results for $y_{c}$ and $x_{c}$, respectively.

In the following, let us first consider $n_{0}(\Omega)=0$ in Eq. (II.3), which
give us a critical $\Omega_{c}$. This is the split critical frequency
$\Omega_{s}$ for $\Omega\leq 1$. However, for $\Omega\geq 1$ it is a critical
rotation for a giant vortex to appear, such that $\Omega=\Omega_{G}$. For
$n_{0}(\Omega_{c})=0$, the critical condition gives

$\displaystyle r_{c}^{2}(\theta)$ $\displaystyle=$
$\displaystyle\frac{2(1+a^{2})}{b^{2}}\left[C(\Omega_{c})-\sin^{2}\theta\right]\Theta\left(C-\sin^{2}\theta\right),$
(25) $\displaystyle{\rm with}\;\;\;C\equiv
C(\Omega_{c})\equiv\frac{a^{2}+\Omega_{c}^{2}}{1+a^{2}}.$

In the above, the Heaviside step function is only effective in case
$\Omega_{c}\leq 1$, when $C$ can define an upper limit for the angle $\theta$,
given by $\theta_{m}$. In such a case, we can assume
$C\equiv\sin^{2}\theta_{m}$ for $\Omega_{c}\leq 1$. In the other case, for
$\Omega_{c}\geq 1$, we keep $C$ as defined in Eq. (25).

By applying the normalization condition in Eq. (II.3), together with Eqs. (20)
and (25), we obtain

$\displaystyle gb^{4}$ $\displaystyle=$
$\displaystyle\frac{b^{6}}{24}\int^{\pi}_{0}d\theta r_{c}^{6}$
$\displaystyle=$
$\displaystyle\Theta(1-\Omega)\frac{2\left(1+a^{2}\right)^{3}}{3}\int^{\theta_{m}}_{0}d\theta\left[C-\sin^{2}\theta\right]^{3}+$
$\displaystyle+$
$\displaystyle\Theta(\Omega-1)\frac{2\left(1+a^{2}\right)^{3}}{3}\int^{\pi/2}_{0}d\theta\left[C-\sin^{2}\theta\right]^{3}$
$\displaystyle=$
$\displaystyle\frac{\left(1+a^{2}\right)^{3}}{48}\left\\{-(2\theta_{m})\cos(2\theta_{m})\left[2\cos^{2}(2\theta_{m})+{3}\right]\right.$
$\displaystyle+$
$\displaystyle\left.\left[{\left({5}-\frac{11}{3}\sin^{2}(2\theta_{m})\right){\sin(2\theta_{m})}}\right]\right\\}\Theta(1-\Omega)$
$\displaystyle+$
$\displaystyle\frac{\pi\left(1+a^{2}\right)^{3}}{3}\left\\{\left(C-\frac{1}{2}\right)\left[C^{2}-C+\frac{5}{8}\right]\right\\}\Theta(\Omega-1)$

In the first case ($\Omega\leq 1$), when the integral is limited by
$\theta_{m}$, we can identify $\Omega_{c}$ with the critical splitting
frequency $\Omega_{s}$. And, in the second case, for $\Omega\geq 1$, the upper
limit of the integral is fixed at $\pi/2$, with the critical value
corresponding to the creation of a giant vortex, such that we identify
$\Omega_{c}$ with $\Omega_{G}$, which is within the definition of $C$, given
in Eq. (25). Therefore, in the above we have
$C\equiv{(a^{2}+\Omega_{G}^{2})}/{(1+a^{2})}$.

Some specific cases:
(i) With $\Omega^{2}\leq 1$, we notice that in the non-rotating case
$\Omega=0$, for $a^{2}=1$, we obtain $gb^{4}={2}/{9}$.

(ii) With $\Omega\geq 1$, for the giant-vortex critical frequency, when
$\Omega_{G}^{2}=1$, we obtain $gb^{4}={5\pi(1+a^{2})^{3}}/{48}$ (= 0.639, if
$a=1/2$); and with $\Omega_{G}=2$, $gb^{4}=52.11$ (with $a=1/2$) and 46.27
(with $a=0$).

Next, we consider the breakup for the second case, $\Omega^{2}\geq 1$, which
is obtained by increasing the rotation frequency beyond $\Omega_{G}$, till the
condensate split in two fragments. So, the density at $x=0$ (for any value of
$y$), which is given by

$\displaystyle n(0,y)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{g}\left[\left(\mu_{TF}-\frac{a^{4}}{4b^{2}}\right)-{\frac{b^{2}}{4}y^{4}-\frac{1-\Omega^{2}}{2}y^{2}}\right],$
(27)

becomes zero along the $y-$axis, at a critical frequency
$\Omega=\Omega_{s}>\Omega_{G}$, when $gn(0,y)=0$. This condition is expressed
by $4b^{2}gn_{0}=-(\Omega_{s}^{2}-1)^{2}$ in (23), implying in
$b^{2}y_{c}^{2}=(\Omega_{s}^{2}-1)$ at the critical breakup limit. Next, by
following the same procedure as before, we consider the corresponding contour
$r_{c}(\theta)$. By applying the normalization (II.3), we obtain

$\displaystyle gb^{4}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{3}(1+a^{2})^{3}\int_{0}^{\pi/2}d\theta\cos^{3}\theta\sqrt{\left[\frac{2\Omega^{2}+a^{2}-1}{a^{2}+1}-\sin^{2}\theta\right]^{3}}.$
(28) $\displaystyle=$
$\displaystyle\frac{2(1+a^{2})^{3}}{3}\int_{0}^{1}dz(1-z^{2}){\left[\zeta^{2}+1-z^{2}\right]^{3/2}}\big{|}_{\zeta^{2}\equiv{\frac{2(\Omega^{2}-1)}{a^{2}+1}}}$
$\displaystyle=$
$\displaystyle\frac{(1+a^{2})^{3}}{72}\left[\zeta\left(3\zeta^{4}+22\zeta^{2}+15\right)\right.$
$\displaystyle-$
$\displaystyle\left.3(\zeta^{2}-5)(\zeta^{2}+1)^{2}\cot^{-1}\zeta\right].$

Therefore, in case $\Omega^{2}=1$ we have $\zeta=0$, with
$gb^{4}={5\pi}{(1+a^{2})^{3}}/{48},$ in agreement with the previous result we
have obtained for $\Omega^{2}\leq 1$.

With Fig. 6 we resume our results for the critical frequencies in terms of
$gb^{4}$, considering two values for the parameter $a$ (=0 and 0.5). For a
better display of the lower rotation frequency, $\Omega<1$, the results are
presented in log-scale in the inset. For such cases with $\Omega<1$, we can
only observe breakup of the condensate, with the limiting condition
$\Omega_{c}$ given when the density at the center reduces to zero. To
illustrate this case, we have the density distribution of the condensate shown
in Fig. 7, where we include three panels with sample results, obtained by the
numerical GP results. With $\Omega=0.6$ [panel (a), before the breakup],
$\Omega=0.75$ [panel (b), close to the breakup], and $\Omega=0.85$ [panel (c),
after the breakup]. To further appreciate the utility of the TF approximation
for the analysis of the radial-quartic trap confinement, we are also including
the Fig. 8 with several panels, in which the upper panels [(a), (b) and (c)],
with $\Omega^{2}<1$, are in agreement with the same numerical results already
presented in Fig. 7 (parameters $a=0.5$, $b^{2}=0.02$ and $g=500$). As also
indicated in Fig. 6, no giant vortex can be formed for these low-rotation
frequency cases, but just breakup is possible. In the lower panels of Fig. 8,
we have also TF results with $\Omega^{2}>1$, corresponding to potential
parameters given by $a=0$, $b^{2}=0.18$ and g=2500, for $\Omega=$2.1 [(d),
just below the formation of a giant vortex], 2.8 [(e), after the giant vortex
is formed] and 4.3 [(f), after the breakup]. All these cases are indicated in
the Fig. 6; with $gb^{4}$ given in a log-scale inset panel for lower
frequencies, and with $gb^{4}$ in normal scale shown in the main panel for
higher frequencies.

Figure 6: The condensate critical frequencies $\Omega$ to originate the giant
vortex ($\Omega_{G}$) (GV, dashed lines) and for the breakup ($\Omega_{s}$)
(BU, solid lines) for $a=0$ and $a=0.5$ (line with bullets). The crossed-
circles are indicating the coordinate positions of density plots shown in the
three panels of Fig. 9 and also for the density plots given in the three lower
panels of Fig. 8 (indicated by the dotted lines at $gb^{4}=$16.2 and 81,
respectively). The inset (in log scale) shows the corresponding behavior for
lower $\Omega<1$, with the crossed-squares indicating the coordinate positions
of the density plots shown in the upper panels of Fig. 8 (indicated by the
dotted line at $gb^{4}=0.2$). With the defined units, all quantities are
dimensionless.  Figure 7:  For $r-$quartic potential [Eq (16)], with $g=500$,
$a=0.5$ and $b^{2}=0.02$ ($gb^{4}=0.2$), these three panels, for $\Omega=0.6$
(a), 0.75 (b) and 0.85 (c), are exemplifying the density profiles
corresponding, respectively, to low frequency conditions, before, close to,
and after the breakup. With the defined units, all quantities are
dimensionless.  Figure 8: Densities by using the TF approximation, for the
$r-$quartic potential. The upper panels are for low frequencies, as in Fig. 7,
with $g=500$, $a=0.5$, $b^{2}=0.02$ ($gb^{4}=0.2$). The lower panels are for
high frequencies, when giant vortex and breakup can be verified, with
$g=2500$, $a=0$, $b^{2}=0.18$ ($gb^{4}=81$). The values of $\Omega$ are
indicated inside the panels. With the defined units, all quantities are
dimensionless.

Figure 9: The upper panels show $|\psi|^{2}$, by using the Eq. (16), with
$g=500$, $b^{2}=0.18$ ($gb^{4}=16.2$) and $a=0$, considering $\Omega=1.65$
(a), 2.00 (b) and 2.75 (c). These panels are exemplifying specific points of
Fig. 6, with giant hole in (b) and after the breakup, in (c). The TF contour
are shown by the white-dashed lines, obtained as solution of Eq. (II.3). The
surface plot shows the effective potential corresponding to panel (c)
[$z\equiv V_{1}(x,y)-(2.75)^{2}/2$]. With the defined units, all quantities
are dimensionless.  Figure 10:  The density with vortex distributions (upper
panel), together with the corresponding phase diagram (lower panel), are shown
for $g=2500$, with $a=0$, $gb^{4}=16.2$ and $\Omega=2.75$, as in Fig. 9(c)
[implying $b^{2}\approx 0.08$]. The results are from numerical GP equation,
with the dashed contour lines given by the TF approximation. With the defined
units, all quantities are dimensionless.

When considering higher frequencies, $\Omega>1$, as we increase $\Omega$ (for
given potential parameter), we first observe a giant vortex being created at
the center, for which the limiting $\Omega$ is identified by $\Omega_{G}$,
represented in Fig. 6 by the dashed lines. In order to verify the splitting of
the cloud, we need to increase $\Omega$, even more, such that
$\Omega_{G}<\Omega$, with the limiting given by $\Omega_{s}$, represented in
Fig. 6 by the solid lines. This case with $\Omega>1$ is also illustrated with
three panels in Fig. 9, for $\Omega=$1.65 [panel (a), before the giant vortex
is created], $\Omega=$ 2 [panel (b), where we observe the giant vortex at the
center, before the breakup], and $\Omega=$2.75 [panel (c), after the breakup].
In these three panels, we also show that the TF approximation can give quite
precisely the contour of the cloud, along all the stages of the rotation-
frequency conditions; namely, before the creation of the giant vortex, during
and after that, and even after the cloud is fragmented in two parts. We are
also illustrating this figure with a 3D surface representation of the
effective potential, corresponding to the panel (c).

In Fig. 9, the panels (a), (b) and (c) refer to the three points shown in Fig.
6 for $gb^{4}=16.2$. One may clearly note the centered-hole phase in Fig. 9
(b) and the breakup in Fig. 9 (c). Corresponding to the three specific cases
shown in Fig. 9, together with the $\Omega=0$ case, in Table 2 we have the TF
and GP chemical potentials, with the GP total energy.

Table 2:  The chemical potentials (GP numerical results and TF approximation), for $\Omega=0$ and for the three cases presented in Fig. 9, are shown in this table. The units are $\hbar\omega$ for the energy and chemical potentials, and $\omega$ for $\Omega$. $\Omega$ | $\mu$ | $\mu_{TF}$ | Energy
---|---|---|---
0 | 15.787 | 15.637 | 9.814
1.65 | 3.808 | 1.484 | -0.199
2 | -6.082 | -8.922 | -10.274
2.75 | -52.490 | -62.741 | -57.417

In view that it was clearly shown that the TF approximation works quite well
to define the contour of the condensate, as shown in the three panels of Fig.
9, we found useful also to verify the effect of increasing the density factor
from $g=500$ to $g=2500$, by considering the panel (c) of Fig. 9. For that, in
Fig. 10 we kept the same values $gb^{4}$ used in Fig. 9, with $a=0$ and
$\Omega=2.75$. As shown, a length-scale transformation is observed, which is
also reflected in the corresponding density. The length-scale transformation
is transparent in the Eqs. (II.3), (18) and (19). As one can see clearly, we
have a scale transformation, such that $(\tilde{x},\tilde{y})=(bx,by)$ implies
$gn(x,y)\to[gb^{2}\tilde{n}(\tilde{x},\tilde{y})]$ (with $\int
d\tilde{x}d\tilde{y}\,\tilde{n}(\tilde{x},\tilde{y})=1$).

By considering this scaling transformation, it is possible to be verified that
both results, shown in Fig. 9(c) and Fig. 10, will correspond to the scaling
of all the corresponding physical quantities, as the trap interaction,
chemical potential and energies (which will be multiplied by a factor
$b^{2}$). Within this scaling procedure, the strength of the nonlinear term
will be rescaled, such that $g\to\overline{g}=gb^{4}$. Therefore, the net
effect due to the kinetic energy, which is obtained by solving the full GP
formalism, will be reflected on the relation of this energy with the frequency
parameter $\Omega$, with the increasing vorticity of the system.

Another aspect which may be of interest, is to estimate numerically the free
expansion of the condensate, by switching off the trap, in order to check with
possible experimental realizations. By considering a limited short time
interval, with $t<16$ (units $\omega^{-1})$, we have verified that the vortex
patterns and corresponding numbers are kept about the same during the
expansion, when considering a case with high frequency in which a giant vortex
can be observed at the center. As expected, the density starts being reduced
as the condensed cloud expands, with the size of the giant vortex increasing
at the same proportion. However, we understand that the required time period
for numerical simulations should be much larger, covering an increasing
spatial distribution, in order to draw a more conclusive picture about this
dynamics, such that could be more useful for experimental studies. This is a
perspective theoretical investigation of great interest to be performed within
a full dynamical study in which one could introduce time-dependent parameters
in the trap interaction.

##  III Vortex pattern distribution and the Feynman rule

As observed from the numerical results, the vortex pattern distribution
follows approximately the corresponding space occupied by the condensate
cloud. The number of vortex is related to the angular frequency of the
rotation and the space occupied by a rotating superfluid, which is the well-
known Feynman rule (FR) 1955-Feynman  established when studying the liquid
helium (For a detailed analysis, considering rapidly rotating condensate with
dense array of vortices, see Ref. 2009Fetter ). According to this rule, the
number of vortices $N_{v}$ inside a given area $A$ is linearly proportional to
the rotation frequency $\Omega$, which in our dimensionless units can be
written as $N_{v}=\Omega A/\pi$. As pointed out in 2010-Wen , this rule was
initially formulated for uniform superfluid helium, been intensively studied
both theoretically and experimentally for a BEC trapped in a single harmonic
potential. So, it is of interest to verify the extension of the rule when
considering trap geometries other than symmetric ones. In our case, by
considering the results shown in Fig. 2, we noticed that apparently it fails
if we consider the visible vortices within the full area defined by the TF
boundaries. However, one should realize that may not be the correct way to
apply the rule, in particular when the area is not regular, having low-density
regions, such that the visible vortex distribution becomes not uniform. In
this regards, the discrepancy in the applicability of the FR has been observed
before for the case of double-well trap confinement 2010-Wen , being
attributed to hidden vortices.

In order to clarify that, let us first consider the cases represented in Fig.
3, where the nonlinear term is much larger than in Fig. 2, with $g=2500$.
Notice that, after the breakup, when $\Omega=0.77$ [panel (b)], the densities
of the two fragments are more or less established visibly within two triangle
formats, if one considers the limiting space defined by the TF contours. In
this case, we can visually approximate the total area as the corresponding one
occupied by two triangles, having $A=2\times 90$ (units of $l_{\omega}^{2}$).
With that, by applying the FR just for the visible vortices, we obtain
$N_{v}\approx 0.77\times(180/\pi)\approx 44$, which matches exactly with the
observed visible vortices inside the area defined by the TF boundary. In the
other case, before the breakup, shown by the panel (a) of Fig. 3, with
$\Omega=0.7$, we have 41 visible vortices inside the TF contour limits, also
in accordance with the corresponding FR prediction. So, the assumption
presented in Ref. 2010-Wen , as requiring hidden vortices to enter in the
estimated value of $N_{v}$, apparently fails when considering these specific
results. Let us try to clarify that. By considering the respective phase
diagrams for the panels (a) and (b), we can clearly verify the existence of
hidden vortices in the low-density regions (see, in particular near the $x=0$
positions). However, they are outside the area defined by the TF contour,
where we are applying the Feynman rule. More precisely, they are emerging only
when we consider a more extended area for the densities, which is given by the
GP numerical solution. So, once consistently we define the respective area in
which enough vortices exist, the applicability of the FR is observed. Suppose,
for instance, that we ignore the TF limiting space, considering a circular
area with center at $(x,y)=0$ and arbitrary radius, covering a space such that
we include the hidden vortices near $x=0$. Again we could verify,
approximately, the applicability of the Feynman rule; this time, including the
hidden vortices.

In resume, the possible regularity of the trap can only be used to explain
visual discrepancy in the applicability of the FR if, inside the given
specific region there are hidden vortices, which need to be included. In Figs.
2 and 3 we have about the same confining potential shape. So, for small non-
linearity, as with $g=500$, the rule apparently is requiring hidden vortices
to be verified (present inside the TF limiting contour), but not in the second
case. In conclusion, the visual verification of the FR in a large area
requires an enough dense cloud, which can support larger number of vortices,
as in the second case, when we have the nonlinear parameter $g=2500$
(indicating larger condensed cloud). With enough dense cloud, one can select
an arbitrary specific area inside the condensate to apply the rule.

By considering a minimum area $A=\pi R_{m}^{2}$ to have one vortex inside, the
rule can be expressed by

$\eta_{v}\equiv\frac{N_{v}}{A}=\frac{\Omega}{\pi},\;\;\to\;\;{R_{m}}=\sqrt{\frac{1}{\Omega}},$
(29)

which also implies in about $2R_{m}$ the distance between two vortices. From
this reading, the FR can be approximately verified, not only for the case that
we have a condensate enough dense with $g=2500$ (where we have verified
complete agreement), but even for the case that we have $g=500$, shown in Fig.
2. This is shown by the Table 3, where $R_{m}$ comes from Eq. (29), with $D$
obtained visually from the respective panels. In this way, the applicability
of the FR can also be extended by including low-density regions, where regular
patterns of hidden vortices can be detected, by examining the corresponding
phase diagrams. See, for example, the results shown in the Fig. 10, in which
hidden vortices (not visible in the upper panel) can be identified (outside
the TF contour line) in the corresponding phase diagram.

Table 3: Applicability of the Feynman rule in the given examples of this contribution, Figs. 2-10. $2R_{m}$ is the predicted distance between two vortices. $D$ is a visual approximate distance between the closest vortices (extracted near the minima of the effective potential, with possible $\pm$0.04 error).  Fig. | $Panels$ | $\Omega$ | $2R_{m}$ | $D$ | $g$
---|---|---|---|---|---
2 | (a) | 0.3 | 3.65 | (111Not applicable, as we have only one vortex near each minimum.) | 500
| (b) | 0.4 | 3.16 | 3.00 | 500
| (c) | 0.5 | 2.83 | 2.60 | 500
| (d) | 0.6 | 2.58 | 2.40 | 500
3 | (a) | 0.70 | 2.39 | 2.40 | 2500
| (b) | 0.77 | 2.28 | 2.30 | 2500
7 | (a) | 0.60 | 2.58 | 2.50 | 500
| (b) | 0.75 | 2.31 | 2.25 | 500
| (c) | 0.85 | 2.17 | 2.00 | 500
9 | (a) | 1.65 | 1.56 | 1.46 | 500
| (b) | 2.00 | 1.41 | 1.37 | 500
| (c) | 2.75 | 1.21 | 1.09 | 500
10 | - | 2.75 | 1.21 | 1.20 | 2500

##  IV Conclusions

The parameter conditions for a single ultra-cold rotating Bose-Einstein
condensate to be fragmented, as well as to form giant vortex at the center,
were investigated in this work, by using two kinds of asymmetric pancakelike
external quartic-quadratic potentials. The corresponding GP formalism was
solved by numerical approach, as well as by a detailed analysis using the TF
approximation. The computational results were obtained for low and high
rotational frequencies, describing the usual vortex-pattern distributions. In
this case, the vortex density distributions was found to follow closely the
well-known Feynman rule, once we have appropriately defined the specific
region where enough vortices can be observed.

Of particular relevance, considering that it can be quite useful in
experimental setup analysis, is TF approximation reliability to define the
limiting frequencies for the density distribution of the condensate, even when
considering the quite non-uniform condensate distributions. The space contour
limits of the condensate, as well as the averaged density distributions, are
shown to be well described by the TF approximation, as compared with the GP
numerical results, even after the breakup of the condensate, which happens as
the rotation frequency increases. These results, verified for different
quartic-quadratic confinement, and also by varying the repulsive interaction
parameter, provides some further support to the TF approximation even when
considering single atomic systems with more complex properties.

As expected, the vortex pattern picture is provided by the GP numerical
approach, in which we noticed a quite regular triangular shape structure for
the vortices, when the rotation is enough high, before the condensate breakup.
The vortex distributions, close to their symmetrically distributed maxima
(minima of the effective potential), apparently are not strongly affected by
the breakup conditions when higher rotation are implemented. In these cases,
we have also verified that the visible vortex pattern distributions follow the
well-known Feynman rule, when applying this rule for sample spacial region
with enough number of vortices. With this rule given in terms of the vortex
density, we can define a unit of area in order to apply the FR even to low-
density regions, where hidden vortices are more likely to appear.

The present study for the breakup conditions in rotating condensate was done
for the simplest possible situation in which we consider an asymmetric
quartic-quadratic trap confinement, enough large to control possible higher
rotation conditions, in which the condensate fragments are still kept
confined. The possible relevance of the present investigation relies in the
fact that we are providing quite well defined limiting conditions via TF
approximation (to our acknowledgment, not presented before in the usual
literature), which can be easily followed in actual experimental realizations,
with different condensed atomic samples, trapped in quartic-quadratic
pancakelike potentials.

We understand this kind of study can be quite useful particularly when
designing experimental setups of other more involved investigations, such as
with dipolar single-component BEC systems, binary miscible and immiscible
combinations, etc. As a perspective interesting study, one could further
investigate the confined system by assuming time-dependent rotation frequency
or trap configurations, which could be introduced by using the present
approach. By manipulating the trap configuration as shown, for example, in
this work, one could implement another mechanism to control the condensed
cloud, in addition to the use of Feshbach resonance techniques.

Acknowledgments
We are grateful to Romain Dubessy for his interest and helpful discussion. The
authors thank the Brazilian agencies Fundação de Amparo à Pesquisa do Estado
de São Paulo (FAPESP) [Contracts 2018/02737-4(AA), 2017/05660-0(LT),
2016/17612-7 (AG)], Conselho Nacional de Desenvolvimento Científico e
Tecnológico [Procs. 304469-2019-0(LT) and 306920/2018-2 (AG)] and Coordenação
de Aperfeiçoamento de Pessoal de Nível Superior [Proc. 88887.374855/2019-00
(LB)].

## References

  * (1) A. Maeder, Physics, formation and evolution of rotating stars (Springer, Berlin, 2009).
  * (2) H. Uchida, M. Shibata, T. Yoshida, Y. Sekiguchi, Yuichiro, and H. Umeda, Gravitational collapse of rotating supermassive stars including nuclear burning effects, Phys. Rev. D 96, 083016 (2017).
  * (3) C. Chiappini, U. Frischknecht, G. Meynet, R. Hirschi, B. Barbuy, M. Pignatari, T. Decressin and A. Maeder, Imprints of fast-rotating massive stars in the galactic bulge, Nature bf 472, 454 (2011).
  * (4) A. Bohr and B. R. Mottelson, Physics of rapidly rotating nuclei, Physics Today 32, 6-25 (1979).
  * (5) R. Clark and B. Wadsworth, A new spin on nuclei, PhysicsWorld, (1998).
  * (6) Y. Guo, R. Dubessy, M. G. de Herve, A. Kumar, T. Badr, A. Perrin, L. Longchambon, and H. Perrin, Supersonic rotation of a superfluid: A long-lived dynamical ring, Phys. Rev. Lett. 124 025301 (2020).
  * (7) R. K. Kumar, T. Sriraman, H. Fabrelli, P. Muruganandam, and A. Gammal, Three-dimensional vortex structures in a rotating dipolar Bose-Einstein condensate, J. Physics B: At., Mol. and Opt. Phys. 49, 155301 (2016).
  * (8) R. A. Barankov, Boundary of two mixed Bose-Einstein condensates, Phys. Rev. A 66, 013612 (2002).
  * (9) R. M. Wilson, C. Ticknor, J. L. Bohn, and E. Timmermans, Roton immiscibility in a two-component dipolar Bose gas, Phys. Rev. A 86, 033606 (2012).
  * (10) L. Wen, W. M. Liu, Y. Cai, J. M. Zhang, and J. Hu, Controlling phase separation of a two-component Bose-Einstein condensate, Phys. Rev. A 85, 043602 (2012).
  * (11) R. W. Pattinson, T. P. Billam, S. A. Gardiner, D. J. Mc-Carron, H. W. Cho, S. L. Cornish, N. G. Parker, and N. P. Proukakis, Equilibrium solutions for immiscible two-species Bose-Einstein condensates in perturbed harmonic traps, Phys. Rev. A 87, 013625 (2013).
  * (12) G. Filatrella, B. A. Malomed, and M. Salerno, Domain walls and bubble droplets in immiscible binary Bose gases, Phys. Rev. A 90, 043629 (2014).
  * (13) R. K. Kumar, P. Muruganandam, L. Tomio, and A. Gammal, Miscibility in coupled dipolar and non-dipolar Bose-Einstein condensates, J. Phys. Commun. 1, 035012 (2017).
  * (14) R. K. Kumar, L. Tomio, and A. Gammal, Spatial separation of rotating binary Bose-Einstein condensates by tuning the dipolar interactions, Phys. Rev. A 99, 043606 (2019).
  * (15) A. L. Fetter, Rotating vortex lattice in a Bose-Einstein condensate trapped in combined quadratic and quartic radial potentials, Phys. Rev. A 64, 063608 (2001).
  * (16) E. Lundh, Multiply quantized vortices in trapped Bose-Einstein condensates, Phys. Rev. A 65, 043604 (2002).
  * (17) G. M. Kavoulakis and G. Baym, Rapidly rotating Bose-Einstein condensates in anharmonic potentials, New J. of Phys. 5, 51 (2003).
  * (18) A. L. Fetter, B. Jackson, and S. Stringari, Rapid rotation of a Bose-Einstein condensate in a harmonic plus quartic trap, Phys. Rev. A 71, 013605 (2005).
  * (19) A. L. Fetter, Rotating trapped Bose-Einstein condensates, Rev. Mod. Phys. 81, 647 (2009).
  * (20) A. Recati, F. Zambelli, and S. Stringari, Overcritical rotation of a trapped Bose-Einstein condensate, Phys. Rev. Lett. 86 377 (2001).
  * (21) K. Kasamatsu, M. Tsubota, and M. Ueda, Giant hole and circular superflow in a fast rotating Bose-Einstein condensate, Phys. Rev. A 66, 053606 (2002).
  * (22) P. Engels, I. Coddington, P. C. Haljan, V. Schweikhard, and E. A. Cornell, Observation of long-lived vortex aggregates in rapidly rotating Bose-Einstein condensates, Phys. Rev. Lett. 90, 170405 (2003).
  * (23) U. R. Fischer and G. Baym, Vortex states of rapidly rotating dilute Bose-Einstein condensates, Phys. Rev. Lett. 90, 140402 (2003).
  * (24) V. Bretin, S. Stock, Y. Seurin, J. Dalibard, Fast rotation of a Bose-Einstein condensate, Phys. Rev. Lett. 92, 050403 (2004).
  * (25) A. Aftalion and I. Danaila, Giant vortices in combined harmonic and quartic traps, Phys. Rev. A 69, 033608 (2004).
  * (26) T. P. Simula, A. A. Penckwitt, and R. J. Ballagh, Giant vortex lattice deformations in rapidly rotating Bose-Einstein condensates, Phys. Rev. Lett. 92, 060401 (2004).
  * (27) G. Baym and C. J. Pethick, Vortex core structure and global properties of rapidly rotating Bose-Einstein condensates, Phys. Rev. A 69, 043619 (2004).
  * (28) M. Correggi, F. Pinsker, N. Rougerie, and J. Yngvason, Rotating superfluids in anharmonic traps: From vortex lattices to giant vortices, Phys. Rev. A 84, 053614 (2011).
  * (29) M. Correggi, F. Pinsker, N. Rougerie, and J. Yngvason, Giant vortex phase transition in rapidly rotating trapped Bose-Einstein condensates, Eur. Phys. J. Spec. Top. 217, 183 (2013).
  * (30) S. K. Adhikari, Stable controllable giant vortex in a trapped Bose-Einstein condensate, Laser Phys. Lett., 16 (2019).
  * (31) S. Pandey, H. Mas, G. Drougakis, et al., Hypersonic Bose-Einstein condensates in accelerator rings, Nature 570, 205-209 (2019).
  * (32) E. Lundh, A. Collin and K.-A. Suominen, Rotational states of Bose gases with attractive interactions in anharmonic traps, Phys. Rev. Lett. 92, 070401 (2004).
  * (33) L. D. Carr and C. W. Clark, Vortices in attractive Bose-Einstein condensates in two dimensions, Phys. Rev. Lett. 97, 010403 (2006).
  * (34) H. Sakaguchi and B. A. Malomed, Localized matter-waves patterns with attractive interaction in rotating potentials, Phys. Rev. A 78, 063606 (2008).
  * (35) J. B. Sudharsan, R. Radha, H. Fabrelli, A. Gammal, and B. A. Malomed, Stable multiple vortices in collisionally inhomogeneous attractive Bose-Einstein condensates, Phys. Rev. A 92, 053601 (2015).
  * (36) F. Dalfovo, S. Giorgini, L. P. Pitaevskii and S. Stringari, Theory of Bose-Einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999).
  * (37) A. Aftalion, Vortices in Bose-Einstein condensates, Birkhäuser Boston, c/o Springer Science+Business Media, New York, NY 10013, USA, 2006. Vol. 67 of Progress in Nonlinear Differential Equations and Their Applications, edited by H. Brezis.
  * (38) C.J. Pethick and H. Smith, Bose-Einstein condensation in dilute gases, Cambridge University Press, Cambridge, England, 2008.
  * (39) L. Pitaevskii and S. Stringari, Bose-Einstein condensation and superfluidity, Oxford University Press, Oxford, 2016.
  * (40) S. Inouye, M.R. Andrews, J. Stenger, H.-J. Miesner, D.M. Stamper-Kurn, and W. Ketterle, Observation of Feshbach resonances in a Bose-Einstein condensate, Nature 392, 151 (1998).
  * (41) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * (42) M. Brtka, A. Gammal and L. Tomio, Relaxation algorithm to hyperbolic states in Gross-Pitaevskii equation, Phys. Lett. A 359, 339 (2006).
  * (43) J. Yang, Newton-conjugate-gradient methods for solitary wave computations, J. Comput. Phys. 228, 7007 (2009).
  * (44) M. Cozzini, A. L. Fetter, B. Jackson, and S. Stringari, Oscillations of a Bose-Einstein condensate rotating in a harmonic plus quartic trap, Phys. Rev. Lett. 94, 100402 (2005).
  * (45) R. P. Feynman, Application of quantum mechanics to liquid helium, Progress in Low Temperature Physics, Vol. I, edited by C. J. Gorter (North-Holland, Amsterdam, 1955).
  * (46) E. W. Weisstein, “Cassini Ovals”. From MathWorld–A Wolfram Web Resource.
https://mathworld.wolfram.com/CassiniOvals.html.

  * (47) L. Wen, H. Xiong and B. Wu, Hidden vortices in a Bose-Einstein condensate in a rotating double-well potential, Phys. Rev. A 82 053627 (2010).