main suppSupplementary References

# Fragmentation and correlations in a rotating Bose-Einstein condensate
undergoing breakup

Sunayana Dutta<EMAIL_ADDRESS>Department of Physics, University of
Haifa, Haifa 3498838, Israel  Haifa Research Center for Theoretical Physics
and Astrophysics, University of Haifa, Haifa 3498838, Israel  Axel U. J. Lode
Institute of Physics, Albert-Ludwig University of Freiburg, Hermann-Herder-
Strasse 3, 79104 Freiburg, Germany  Ofir E. Alon  Department of Physics,
University of Haifa, Haifa 3498838, Israel  Haifa Research Center for
Theoretical Physics and Astrophysics, University of Haifa, Haifa 3498838,
Israel

###### Abstract

The theoretical investigation of rotating Bose-Einstein condensates has mainly
focused on the emergence of quantum vortex states and the condensed properties
of such systems. In the present work, we concentrate on other facets by
examining the impact of rotation on the ground state of weakly interacting
bosons confined in anharmonic potentials computed both at the mean-field level
and particularly at the many-body level of theory. For the many-body
computations, we employ the well-established many-body method known as the
multiconfigurational time-dependent Hartree method for bosons (MCTDHB). We
present how various degrees of fragmentation can be generated following the
breakup of the ground state densities in anharmonic traps without ramping up a
potential barrier for strong rotations. The breakup of the densities is found
to be associated with the acquisition of angular momentum in the condensate
due to the rotation. In addition to fragmentation, the presence of many-body
correlations is examined by computing the variances of the many-particle
position and momentum operators. For strong rotations, the many-body variances
become smaller than their mean-field counterparts, and one even finds a
scenario with opposite anisotropies of the mean-field and many-body variances.
Further, it is observed that for higher discrete symmetric systems of order
$k$, namely three-fold and four-fold symmetry, breakup to $k$ sub-clouds and
emergence of $k$-fold fragmentation take place. All in all, we provide a
thorough many-body investigation of how and which correlations build up when a
trapped Bose-Einstein condensate breaks up under rotation.

##  I Introduction

The successful experimental realization of a rotating Bose-Einstein condensate
(BEC) has paved the way to explore various rich physics of correlated quantum
systems [1, 2, 3, 4, 5]. Butts and Rokhsar [6] first evaluated the wave
function of a rotating BEC using the lowest Landau level approximation with
the help of a Gross-Pitaevskii functional. The rotating ultracold bosonic
gases have led to the investigation of the occurrence of quantized vortices
[6, 3, 1], vortex nucleation [7], emergence of quantum fluctuations [8, 9],
and presence of the fractional quantum Hall effect in weakly interacting
quantum systems [10].

There have been studies where breaking up of fast rotating objects is widely
detected in systems extending from astronomical objects, e.g., galaxies and
supermassive rotating stars [11] to the quantum systems in nuclear physics
[12]. In atomic physics, the emergence of superfluid flow of a rotating
quantum gas has been explored experimentally in Ref. [13], in an anharmonic
potential. Theoretically, the emergence of breakup in a rotating dipolar
condensate was investigated in three-dimension in Ref. [14] and finally, in a
rotating pancake-like asymmetric quartic-quadratic potential in Ref. [15]. A
substantial volume of literature exists corresponding to the rotating BEC,
investigating various rich quantum features by employing Gross-Pitaevskii
mean-field approximation [16, 17, 18]. However, there exists much less studies
of rotating BEC in the many-body domain, see [19, 20, 21, 22, 23, 24], that
explore interesting many-body quantum features like the fragmentation of the
condensate and correlations.

Condensation and fragmentation are widely explored many-body features of BEC
derived from the properties of the one-body reduced density matrix [25, 26,
27, 28, 29, 30]. According to Penrose and Onsager, interacting bosons are said
to be condensed if they have a single macroscopic eigenvalue of the one-body
reduced density matrix [31], and fragmented if there exist two or more
macroscopic eigenvalues [32]. The fragmentation of condensates has been
thoroughly studied for non-rotating systems [33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 29, 30, 46, 47]. However, these many-body features have
relatively been less extensively investigated in the rotating frame [48, 7].

In the regime of “ultrafast rotation”, when the rotation frequency comes
closer to the trapping frequency in case of a harmonic trapping potential, the
system would tend to escape as the centrifugal force would cancel the trapping
force. This problem can be resolved either by introducing an anharmonic term
into the confining potential or by adding anisotropy in the harmonic
potential. Anisotropy makes the confining potential elongated. Ref. [49]
investigated the appearance of vortex rows of BEC in an anisotropic potential
under rotation using Gross-Pitaevskii formalism. Motivated also by these
findings, in this work, we study a system of bosons subject to rotation
confined in various anharmonic trap geometries with discrete rotational
symmetry in two-dimension 2D, namely an elongated trap, a three-fold symmetric
trap and finally four-fold symmetric trap.

There exist several ideas to induce mechanical rotation and coupling of the
internal states in the condensate, such as, phase imprinting (using an
electromagnetic field) as proposed by Williams and Holland in Ref. [50] and
optical spoon stirring method where quantized vortices are observed in a
stirred gaseous condensate of atomic rubidium [1]. Here, we study the physics
of a 2D weakly interacting ultracold bosons confined in different anharmonic
potentials both at the mean-field level and particularly at the many-body
level of theory. Specifically, we investigate how the system transformed from
a fully-condensed state into a fragmented state followed by breaking up of the
condensate density induced by the rotational motion. As variance is a
sensitive probe that characterizes many-body correlations even for fully
condensed systems [51, 52], hence we also analyze the variances of many-
particle operators of the fragmented state in the rotating frame. The
qualitative difference between the mean-field and many-body variances is a
useful tool to explain the nature of many-body correlations. The many-particle
position variance characterizes to what extent a wave-packet spreads or
narrows down, similarly the momentum variance is associated with the size of
the wave-packet in momentum space. Hence, we emphasize on the emergence of
correlations by investigating variances of many-particle operators like the
position, momentum and angular momentum. Interestingly, the fluctuations
present in the system are hardly observed in the angular momentum variance.
Thus, we only present the position and momentum variances in the main text and
give the detail discussion of the angular momentum variance in the
supplemental material.

The time-dependent Gross-Pitaevskii mean-field theory [53, 54] is the most
celebrated theoretical model to investigate the many-particle systems of
ultracold bosonic atoms. However, this method is unable to study fragmentation
and correlations owing to its building via the mean-field ansatz. In this
paper, we employ a well established many-body numerical method named the
multiconfigurational time-dependent Hartree method for bosons (MCTDHB) [55,
56] to accurately solve the Schrödinger equation at the many-body level for
ultracold atoms subject to a rotation. The MCTDHB method is a bosonic version
of the MCTDH family of methods [57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
68, 69, 70, 71, 72] which is able to self-consistently describe the physics
involving the presence of many-body correlations. The main focus of the
applications of MCTDHB has been the emergence of fragmentation of the
condensate, where the one-body reduced density matrix has multiple significant
eigenvalues. For numerical simulation of the results presented in this work,
we use the MCTDH-X software [72, 73, 74, 75]. Finally, the supplemental
material reports the benchmarking of MCTDHB for an exactly solvable many-body
model under rotation and also presents the convergence of the many-body
results of our present work.

##  II Setup and theoretical tools

We consider a system of weakly interacting bosonic atoms in two spatial
dimensions 2D confined in non-spherically symmetric trapping potentials in the
rotating frame. The properties of these trapped bosons can be described by the
(time-dependent) many-body Schrödinger equation. The Schrödinger equation
dealing with a many-boson system is usually solved by employing the mean-field
Gross-Pitaevskii approximation. However, the reduced density matrix involved
in the Gross-Pitaevskii approximation has only a single eigenvalue and it
involves a single basis state and thereby is unable to capture the many-body
features such as fragmentation and correlations.

In MCTDHB, the (time-dependent) optimized one-body basis is used. Here the
basis set and the expansion coefficients in the basis are optimized
variationally [55, 56]. The MCTDHB is a numerically exact method [76] and can
describe both coherent and fragmented condensates. MCTDHB includes the theory
of Gross-Pitaevskii approximation as a special case when only a single one-
body state is considered.

###  II.1 Hamiltonian

The general Hamiltonian of $N$ interacting bosons is given as

$\displaystyle\hat{H}=\sum_{j=1}^{N}\hat{h}(\mathbf{r}_{j})+\sum_{j<k}\hat{W}(\mathbf{r}_{j}-\mathbf{r}_{k}),$
(1)

where the single-particle Hamiltonian

$\displaystyle\hat{h}({\bf r})=\hat{T}(\mathbf{r})+\hat{V}(\mathbf{r})$ (2)

is composed of the kinetic energy and the external potential energy,
respectively. Here, the interaction of ultracold dilute bosonic gases is
considered to be a finite range interaction and modelled by a Gaussian
function [77, 63, 64],
$\hat{W}(\mathbf{r}-\mathbf{r}^{\prime})=\frac{\lambda_{0}}{2\pi\sigma^{2}}e^{-\frac{({\bf
r}-{\bf r^{\prime}})^{2}}{2\sigma^{2}}}$ with $\sigma=0.25$. This avoids the
regularization of the delta contact potential in 2D. The interaction strength
$\lambda_{0}$ is scaled with the number of bosons $N$ as
$\Lambda=\lambda_{0}(N-1)$, where $\Lambda$ is the interaction parameter. One
uses the interaction parameter to define the mean-field regime. In our study,
we work in the units $\hbar=m=1$ and all the quantities are dimensionless. We
also consider three different trapping potentials $\hat{V}({\bf r})$ that we
shall discuss in the next section. The first setup is the elongated trap $-$
that leads to breaking up of the ground state density into two clouds. After
that we move to more complex traps, namely three-fold symmetric and four-fold
symmetric traps to investigate and establish the generality of the results.
Hence, our strategy is to first study an elongated trap, and then a three-fold
symmetric trap, and finally, a four-fold symmetric trap to see what stays
between the two to three and the three to four-fold symmetric traps.

In the rotating frame, the kinetic energy operator is modified and can be
written as,

$\displaystyle\hat{T}({\bf
r})=\frac{1}{2}(\hat{p}_{x}^{2}+\hat{p}_{y}^{2})-\omega_{r}\hat{l}_{z},$ (3)

here $\omega_{r}$ is the rotation frequency and
$\hat{l}_{z}=\hat{x}\hat{p}_{y}-\hat{y}\hat{p}_{x}$ is the angular-momentum
operator.

An alternative way to mimic the rotational effect in the condensate is by
introducing a synthetic gauge field $\mathbf{A}(\mathbf{r})$ as

$\hat{T}(\mathbf{r})=\frac{1}{2}\big{[}-i\nabla_{\mathbf{r}}-q\mathbf{A}(\mathbf{r})\big{]}^{2}.$
(4)

Consider the following general form of the gauge field:

$\mathbf{A}(\mathbf{r})=(ay,bx,0).$ (5)

Then, expansion of Equation (4) leads to,

$\displaystyle\hat{T}(\mathbf{r})=\frac{1}{2}(\hat{p}_{x}-ay)^{2}+\frac{1}{2}(\hat{p}_{y}-bx)^{2}=\frac{1}{2}(\hat{p}_{x}^{2}+\hat{p}_{y}^{2})-(\hat{p}_{x}ay+\hat{p}_{y}bx)+\frac{1}{2}(a^{2}y^{2}+b^{2}x^{2}).$
(6)

For the specific case $b=-a$, Equation (6) becomes

$\hat{T}(\mathbf{r})=\frac{1}{2}(\hat{p}_{x}^{2}+\hat{p}_{y}^{2})-a\hat{l}_{z}+\frac{1}{2}a^{2}(y^{2}+x^{2}).$
(7)

Combining Equations (2), (3) and (7) we have

$\hat{h}({\bf r})=\hat{T}(\mathbf{r})+\hat{V}^{\prime}(\mathbf{r}),$ (8)

where the modified confining potential is
$\hat{V}^{\prime}(\mathbf{r})=\hat{V}(\mathbf{r})+\frac{1}{2}a^{2}\mathbf{r}^{2}$
and $a=\omega_{r}$ corresponds to the rotation frequency of the condensate.

###  II.2 Many-body method

The MCTDHB method employs time-adaptive orbitals to represent the field
operator as a sum of the $M$ time-dependent single-particle states

$\hat{\Psi}(\mathbf{r},t)=\sum_{j=1}^{M}\hat{b}_{j}\phi_{j}(\mathbf{r},t).$
(9)

The ansatz of the MCTDHB wavefunction is

$|\Psi(t)\rangle=\sum_{\vec{n}}C_{\vec{n}}|\vec{n},t\rangle.$ (10)

The summation in Equation (10) runs over all $\binom{N+M-1}{N}$ possible time-
dependent configurations $\vec{n}=(n_{1},...,n_{M})$ with fixed particle
number $N=\sum_{i=1}^{M}n_{i}$. To derive the MCTDHB equations, the time-
dependent variational principle [78, 79, 80] is employed for the ansatz in
Equation (10). Thus, resulting in two-coupled equations of motion – a set of
linear equations for the coefficients $\\{C_{\vec{n}}\\}$ and a set of non-
linear equations for the orbitals $\\{\phi_{j}(\mathbf{r});j=1,...,M\\}$, see
Refs. [56, 55] for details and derivation of the equations of motion. In the
following work, the self-consistent ground state is achieved by relaxing the
system via imaginary-time propagation and is hence determined by the
variational principle. Thus, in the following sections, we omit the time-
dependency from the various quantities and observables that are involved in
the many-body simulations.

###  II.3 Quantities of interest

In this section, we define the quantities of interest, namely the one-body
density, the eigenvalues of the one-body reduced density matrix (RDM), the
expectation value of the angular momentum operator, and finally, the many-
particle variances of the position, momentum, and angular momentum operators.

####  II.3.1 One-body reduced density matrix (RDM), one-body density, and
natural occupations

The one-body RDM of the N-boson state $|\Psi\rangle$ is a hermitian matrix and
is defined as

$\rho^{(1)}(\mathbf{r},\mathbf{r}^{\prime})=\langle\Psi|\hat{\Psi}^{\dagger}(\mathbf{r}^{\prime})\hat{\Psi}(\mathbf{r})|\Psi\rangle=\sum_{k,q}\rho_{kq}\phi^{*}_{k}(\mathbf{r}^{\prime})\phi_{q}(\mathbf{r})$
(11)

in its eigenbasis $\\{\phi_{q}(\mathbf{r})\\}$. The matrix elements
$\rho_{kq}=\langle\Psi|\hat{b}^{\dagger}_{k}\hat{b}_{q}|\Psi\rangle$ represent
the one-body RDM using $M$ orbitals corresponding to the creation
(annihilation) operators $\hat{b}^{\dagger}_{k}$ ($\hat{b}_{q}$). The diagonal
of $\rho^{(1)}(\bf{r},\bf{r}^{\prime})$ is referred to as the one-body density
$\rho(\mathbf{r})$ which is
$\rho(\mathbf{r})=\rho^{(1)}(\mathbf{r},\mathbf{r}^{\prime}=\mathbf{r})$.

The eigenvalues of the one-body RDM are obtained by the diagonalization of
Equation (11) which corresponds to a unitary transformation of the orbitals
$\phi_{q}(\mathbf{r})$ to the natural orbitals $\phi_{j}^{(NO)}(\mathbf{r})$
as

$\frac{\rho^{(1)}(\mathbf{r},\mathbf{r}^{\prime})}{N}=\sum_{j}n_{j}\phi_{j}^{(NO),*}(\mathbf{r}^{\prime})\phi_{j}^{(NO)}(\mathbf{r}).$
(12)

Here, the eigenvalues $n_{j}$ are normalized as $\sum_{j=1}^{M}n_{j}=1$ and,
without loss of generality, they are sorted in magnitude such that $n_{1}\geq
n_{2}\geq...$ throughout this work. The eigenvalues $n_{j}$ are termed natural
occupations and characterize the degree of condensation and fragmentation of
the bosons. Thus, the system with one-body RDM consisting only a single
macroscopically-contributing eigenvalue $n_{1}$ is said to be condensed [31].
When the one-body RDM has $k$ macroscopically-occupied eigenvalues, the system
is referred to as $k$-fold fragmented [34].

####  II.3.2 Angular momentum

For a 2D many-particle systems, there is only a single component of the
angular momentum operator, i.e.,

$\displaystyle\hat{L}_{Z}=\sum_{j=1}^{N}\hat{l}_{z_{j}}=\sum^{N}_{j=1}\frac{1}{i}\left(\hat{x}_{j}\frac{\partial}{\partial
y_{j}}-\hat{y}_{j}\frac{\partial}{\partial x_{j}}\right).$ (13)

Bosonic systems with angular momentum provide rich quantum features beyond
mean-field theory, such as, phantom vortices [81], spatially partitioned
vortices [51], and fragmentation [82, 51, 81]. In the following studies, we
investigate the expectation value of the angular momentum operator per
particle $\frac{1}{N}\langle\Psi|\hat{L}_{Z}|\Psi\rangle$, for three different
confining traps to see the effect of the rotation. In our system, we expect an
intricate dynamics of the angular momentum acquisition under rotation and
build up of correlations.

####  II.3.3 Many-particle variances

The variance of a many-particle observable
$\hat{O}=\sum_{j=1}^{N}\hat{o}(\mathbf{r}_{j})$ can be written as [83]

$\displaystyle\frac{1}{N}\Delta^{2}_{\hat{O}}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\bigg{[}\langle\Psi|\hat{O}^{2}|\Psi\rangle-\langle\Psi|\hat{O}|\Psi\rangle^{2}\bigg{]}=$
(14) $\displaystyle=$
$\displaystyle\frac{1}{N}\bigg{\\{}\sum_{j}n_{j}\int\mathrm{d}\mathbf{r}\>\phi^{*(NO)}_{j}(\mathbf{r})\hat{o}^{2}\phi^{(NO)}_{j}(\mathbf{r})-\Bigg{[}\sum_{j}n_{j}\int\mathrm{d}\mathbf{r}\>\phi^{*(NO)}_{j}(\mathbf{r})\hat{o}\phi^{(NO)}_{j}(\mathbf{r})\Bigg{]}^{2}+$
$\displaystyle+$
$\displaystyle\sum_{jklm}\rho_{jklm}\bigg{[}\int\mathrm{d}\mathbf{r}\phi^{*(NO)}_{j}(\mathbf{r})\hat{o}\phi^{(NO)}_{l}(\mathbf{r})\bigg{]}\bigg{[}\int\mathrm{d}\mathbf{r}\phi^{*(NO)}_{k}(\mathbf{r})\hat{o}\phi^{(NO)}_{m}(\mathbf{r})\bigg{]}\bigg{\\}}.$

Here, the expectation value of $\hat{O}=\sum_{j=1}^{N}\hat{o}(\mathbf{r}_{j})$
is dependent only on the one-body operator, whereas the expectation of
$\hat{O}^{2}$ is a combination of one- and two-body operators
$\hat{O}^{2}=\sum_{i=1}^{N}\hat{o}^{2}(\mathbf{r}_{j})+\sum_{j<k}2\hat{o}(\mathbf{r}_{j})\hat{o}(\mathbf{r}_{k})$.
$\rho_{jklm}$ are the two-particle reduced density matrix elements,
$\rho(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}^{\prime}_{1},\mathbf{r}^{\prime}_{2})=\sum_{jklm}\rho_{jklm}\phi^{*}_{j}(\mathbf{r})\phi^{*}_{k}(\mathbf{r})\phi_{l}(\mathbf{r})\phi_{m}(\mathbf{r})$.
In the following work, to analyze the emergence of many-body correlations in
the rotating condensate, we study the many-particle variances per particle of
the position, momentum, and angular momentum operators [52].

##  III Many-body physics of Bose-Einstein condensates breaking up under
rotation: Results and Discussion

We investigate the impact of rotation on the ground state of weakly
interacting bosonic atoms by dividing the analysis into two main parts
depending on the confining anharmonic potentials. First, the breakup and
fragmentation processes in an elongated trap are investigated and then the
breakup and fragmentation in more complicated traps of discrete spatial
symmetry are explored. In the following sections, we analyze the effect of
rotation on various static properties, namely the ground state energy, one-
body density, natural occupations, expectation value of the angular momentum
operator, and finally, the variances of the many-particle position, momentum,
and angular momentum operators as a function of the rotation frequency
$\omega_{r}$. For the numerical computations, the MCTDH-X implementation of
the MCTDHB theory is employed [72, 73, 74, 75]. In our work, we consider
$N=10$ weakly interacting ultracold bosonic atoms interacting via a Gaussian
interaction [77, 63, 64] with the interaction parameter chosen as
$\Lambda=0.1$ throughout the computations. The grid used to represent the
Hamiltonian Equation (18), see section Sec.II.1, extends from $[-8,8)$ to
$[-8,8)$ and comprises of $128\times 128$ discrete variable representation
(exponential) functions to represent each of the orbitals
${\phi_{j}(\bf{r})}$.

Figure 1: Behaviour of the ground state energy per particle E/N as a function
of the rotation frequency $\omega_{r}$ in an elongated trap computed at the
mean-field (MF) and many-body (MB) levels (with M=1 and M=4 self-consistent
orbitals, respectively). The inset shows the energy difference between the
mean-field and many-body energies. All quantities shown are dimensionless.

###  III.1 Breaking up to two clouds and pathway to two-fold fragmentation

Let us discuss the impact of rotation on the weakly interacting bosons trapped
in an elongated confinement. The elongated trap is an anharmonic potential
that elongates the condensate in the $x$-direction and can be written as

$\displaystyle V(\mathbf{r})=\frac{1}{4}(0.8x^{2}+y^{2})^{2}.$ (15)

For this trap, we consider $M=4$ self-consistent orbitals to investigate the
ground state properties of the condensate. We recomputed these results with
$M=8$ self-consistent orbitals to check the convergence of the system, see the
elaborate discussion in the supplemental material. The rotation frequency
range is $\omega_{r}=[0,2.0]$ in this trap.

Figure 1 shows the behaviour of the ground-state energy per particle E/N, in
the rotating frame computed both at the mean-field and many-body levels.
Initially, E/N remains almost constant for slow rotation. Then, E/N drops
gradually with further increase in $\omega_{r}$ as evident from Figure 1. It
is also observed that the energies computed both at the mean-field and many-
body

Figure 2: The one-body densities per particle are shown for three different
rotation frequencies $\omega_{r}$ in the elongated trap at the many-body
level. M=4 self-consistent orbitals are used. The density computed at the
mean-field (not shown) and many-body levels depict identical features for all
$\omega_{r}$. All quantities shown are dimensionless.

levels practically coincide with each other for all $\omega_{r}$. The inset of
Figure 1 corresponds to the energy difference per particle $\Delta E/N$,
between the mean-field and many-body energies defined as $\Delta
E=E_{MF}-E_{MB}$. The energy difference remains minimum till about
$\omega_{r}=1.2$. Subsequently, $\Delta E/N$ exhibits some structures for an
intermediate range of $\omega_{r}$ around $\omega_{r}\sim 1.3$. Finally, the
energy difference slightly rises from of the order of $10^{-5}$ to $10^{-3}$
for larger rotation frequencies $\omega_{r}\geq 1.6$. The presence of these
structures at the intermediate rotation might suggest that something
interesting is happening at the many-body level. Therefore, we dig deeper than
the energy of the system to see the many-body features.

Now let us discuss the behaviour of the ground-state densities per particle
$\frac{\rho(\mathbf{r})}{N}$ of the rotating condensate confined in the
elongated trap. The densities computed at the many-body level of theory are
shown for three different rotation frequencies $\omega_{r}$ in Figure 2. In
absence of rotation, the density displays a single cloud where all bosons
accumulate in the center of the trap. This behaviour persists with the
inclusion of rotation for slow rotation frequencies, e.g., at
$\omega_{r}=0.5$. Further increase in rotation induces a breaking of the
density into two clouds. The breakup of the ground-state density confined in
the potential given by Equation (34) corresponds to the scenario where the
minimum of the potential is splitted into two parts and shifted with rotation
by creating an effective double-well potential. It is also observed that the
distance between the two densities increases with increase of $\omega_{r}$
[Figure 2(b)-(c)]. We have also computed the density per particle at the mean-
field level and identical features are observed in the density profile. Hence,
it is observed that the densities computed at the mean-field and many-body
levels in the real space show identical pattern. In the density profile, we
show three specific frequencies that correspond to slow, fast, and faster
rotations.

To understand whether the splitting of the density into two clouds [as shown
in Figure 2(c)] for faster rotations leads to fragmentation of the condensate,
we further discuss the dependence of the natural occupations
$\frac{n_{j}}{N}$. In the elongated trap, four natural occupations
corresponding to the four natural orbitals are employed and found to vary with
the rotation frequency $\omega_{r}$ (see Figure 3).

Figure 3: Pathway from condensation to fragmentation in the rotating elongated
trap. two-fold fragmentation is observed with increase in rotation. The two
leading natural occupations $n_{1}/N$ and $n_{2}/N$, are shown as a function
of the rotation frequency $\omega_{r}$. $M=4$ self-consistent orbitals are
used. The third and fourth natural occupations satisfy $n_{3}/N$, $n_{4}/N\leq
10^{-5}$, see the supplemental material. All quantities shown are
dimensionless.

The system remains fully condensed, i.e., $\frac{n_{1}}{N}\sim 1$,
$\frac{n_{2}}{N}\sim\frac{n_{3}}{N}\sim\frac{n_{4}}{N}\leq 10^{-5}$ with the
inclusion of rotation till about $\omega_{r}=1.1$ [see Figure S2(a) of the
supplemental material, which displays the depletion as a function $\omega_{r}$
on a log scale]. As $\omega_{r}$ increases further, the first natural
occupation number $n_{1}/N$, falls of gradually following an increase in
population of the second natural occupation number $n_{2}/N$. The other two
natural occupation numbers remain almost the same as
$\frac{n_{3}}{N}\sim\frac{n_{4}}{N}\leq 10^{-5}$. For faster rotation at
$\omega_{r}=2$, the state becomes essentially fully two-fold fragmented with
natural occupations of $\frac{n_{1}}{N}\approx\frac{n_{2}}{N}\approx 50\%$.
This signifies equally populated two leading natural orbitals, whereas
$\frac{n_{3}}{N}$ and $\frac{n_{4}}{N}$ remain essentially unpopulated. For an
intermediate frequency, $\omega_{r}=1.3$, we observe a scenario where $n_{1}$
shows a deep followed by a peak in $n_{2}$. To understand this feature, we
zoom in at the intermediate points between the rotation frequencies
$\omega_{r}=1.2$ and $\omega_{r}=1.4$ as shown in the inset of Figure 3. A
smooth transition is found, from coherence via loss of coherence to build up
of coherence. This transition might suggest the presence of a resonant-like
behaviour of the interacting bosons in the elongated 2D trap at this specific
rotation frequency.

To intermediately summarize, for weakly interacting bosons confined in an
elongated trap [Equation (34)], inclusion of rotation triggers a transition
from a fully condensed state to a fully two-fold fragmented state with equally
populated two leading natural orbitals. Thus, it can be concluded that the
rotation can be used as a tool to manipulate fragmentation.

Figure 4: Expectation value of angular momentum operator
$\langle\hat{L}_{Z}\rangle/N$, computed at the mean-field (MF) and many-body
(MB) levels [with $M=1$ and $M=4$ self-consistent orbitals, respectively] as a
function of the rotation frequency $\omega_{r}$ for the elongated trap. Actual
data are points. The continuous curves are only to guide the eye. All
quantities shown are dimensionless.

Next, let us discuss the behaviour of the expectation value of the angular
momentum per particle $\frac{\langle\psi|\hat{L}_{Z}|\psi\rangle}{N}$,
computed at the mean-field and many-body levels as a function of the rotation
frequency $\omega_{r}$. It is observed from Figure 4 that the mean-field and
many-body angular momenta exactly coincide each other at all $\omega_{r}$
except for the intermediate rotation frequency around $\omega_{r}=1.3$ [see
the inset of Figure 4]. The angular momentum remains minimum till about
$\omega_{r}=1.1$. For $\omega_{r}=1.2$ (the rotation frequency that
corresponds to the breakup of the density), the rotation generates a state
where significant angular momentum enters the system with
$\langle\hat{L}_{Z}\rangle/N>1$. The angular momentum of the condensate
gradually increases with further increase of the rotation.

As we know, the variance is a sensitive probe of correlations that allows one
to study the quantum fluctuations present in a system [51, 52]. Thus, it would
be interesting to investigate the variance of many-particle operators which
signifies many-body correlations for the fragmented condensate in the rotating
frame. Therefore, we further analyze the impact of rotation on the behaviour
of the many-particle variances of position and momentum operators which are
sensitive to rotation.

Figure 5: Dependence of the the many-particle position and momentum variances
on the rotation in the elongated trap. Shown are (a)
$\frac{1}{N}\Delta^{2}_{\hat{X}}$, (b) $\frac{1}{N}\Delta^{2}_{\hat{Y}}$, (c)
$\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}$, and (d)
$\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}$ as a function of $\omega_{r}$ at the
many-body level (MB) [$M=4$ self-consistent orbitals] and at the mean-field
level (MF) [$M=1$ self-consistent orbital]. All quantities shown are
dimensionless.

Figures 5(a)-(b) display the behaviour of the many-particle position variance
per particle $\frac{1}{N}\Delta^{2}_{\hat{X},\hat{Y}}$ along the $x$ and
$y$-directions respectively, as a function of the rotation frequency
$\omega_{r}$ in the elongated trap. The mean-field and many-body position
variances along the $x$\- direction $\frac{1}{N}\Delta^{2}_{\hat{X}}$,
coincide till about $\omega_{r}=0.9$, indicating the essential absence of
correlations in the system [Figure 5(a)]. Further, the mean-field position
variance monotonically increases with $\omega_{r}$ which signifies spreading
of the density as observed in Figure 2. Also, the many-body position variance
slowly increases till $\omega_{r}=1.2$. But, now, we observe a deep at
$\omega_{r}=1.3$ that corresponds to the loss of coherence in the natural
occupation, see Figure 3. After increasing and reaching a maximum at rotation
frequency of about $\omega_{r}\sim 1.4$, $\frac{1}{N}\Delta^{2}_{\hat{X}}$
starts to decrease which incorporates the emergence of a small amount of
depletion of the condensate [see Figure 5(a)]. For a faster rotation, at
$\omega_{r}=2$, the position variance $\frac{1}{N}\Delta^{2}_{\hat{X}}$
decreases significantly which goes hand in hand with the emergence of
fragmentation of the condensate. Similar behaviour of the position variance
along the $x$-direction is observed in a two-dimensional double well in [84],
albeit without the resonant-like behaviour described above. The position
variance per particle along the $y$-direction
$\frac{1}{N}\Delta^{2}_{\hat{Y}}$ is almost frozen and varies slowly as shown
in Figure 5(b). This might suggest that excitations along the tighter
$y$-directions are practically not involved, at least as far as the position
variance is considered. The $\frac{1}{N}\Delta^{2}_{\hat{Y}}$ computed at the
mean-field level matches that at the many-body level almost for all the
rotation frequencies except of small difference for the intermediate
$\omega_{r}$. As above, convergence of the results is detailed in the
supplemental material.

Further, we can describe the anisotropy of the variance by considering two
facts,
$(1)$ by comparing two quantities, one along the $x$\- direction and the other
along the $y$\- direction, whether they are similar or different.
$(2)$ by comparing these quantities at the mean-field and many-body levels
which demonstrates whether the anisotropy of the many-particle variances are
alike or opposite. Hence, in the elongated trap and for small rotations, it is
observed that

$\frac{1}{N}\Delta^{2}_{\hat{X}}|_{MB}<\frac{1}{N}\Delta^{2}_{\hat{Y}}|_{MB}$,
at MB

$\frac{1}{N}\Delta^{2}_{\hat{X}}|_{MF}<\frac{1}{N}\Delta^{2}_{\hat{Y}}|_{MF}$,
at MF.

Thus, the anisotropy of the many-particle position variance computed at the
mean-field and many-body levels are alike at slow rotation. However, for fast
rotation we find

$\frac{1}{N}\Delta^{2}_{\hat{X}}|_{MB}<\frac{1}{N}\Delta^{2}_{\hat{Y}}|_{MB}$,
at MB

$\frac{1}{N}\Delta^{2}_{\hat{X}}|_{MF}>\frac{1}{N}\Delta^{2}_{\hat{Y}}|_{MF}$,
at MF.

Hence, it indicates that the many-particle position variances display opposite
anisotropy when computed at the mean-field and many-body levels.

Now, we investigate the behaviour of the many-particle momentum variance per
particle $\frac{1}{N}\Delta^{2}_{\hat{P}_{X},\hat{P}_{Y}}$ along the $x$\- and
$y$-directions as a function of rotation frequency $\omega_{r}$ in the
elongated trap [Figure 5(c)-(d)]. It is observed that unlike the position
variance, the mean-field and many-body momentum variances along the
$x$-direction are almost similar, see Figure 5(c). The momentum variances
computed at the mean-field and many-body levels gradually increase from
$\omega_{r}=1.2$ onwards. This corroborates the narrowing of the density lobes
along the $x$-direction in real space. The momentum variance in the
$y$-direction displays a completely different picture, see Figure 5(d). The
mean-field and many-body variances match each other till $\omega_{r}=1.1$.
However, for a faster rotation, the momentum variance computed at the mean-
field deviates from that computed at the many-body level in the $y$-direction.
$\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}$ displays similar behaviour as the
$\frac{1}{N}\Delta^{2}_{\hat{X}}$ as clear from Figure 5(a). Thus, excitations
along the $y$-direction plays a role in the momentum space as a result of
rotation.

Inverse to the position variance, in the case of the momentum variance for
slow rotation we find

$\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}|_{MB}>\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}|_{MB}$,
at MB

$\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}|_{MF}>\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}|_{MF}$,
at MF.

Hence, the anistropy of the momentum variances are alike when computed at the
many-body and mean-field levels of theory for slow rotation. For a faster
rotation however, we find

$\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}|_{MB}>\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}|_{MB}$,
at MB

$\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}|_{MF}<\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}|_{MF}$,
at MF.

This signifies the presence of opposite anisotropy of the momentum variance.
This feature was not found before for static double well [84]. Hence, this
indicates one of the distinct features of the rotation, which makes both the
many-particle position and momentum variances show opposite anisotropies with
respect to the mean-field and many-body levels in the elongated trap. The
rotation localized the position variance in the long direction and the
momentum variance in the narrow direction and this is purely a many-body
effect.

Let us briefly summarize the results so far. The ground state of a Bose-
Einstein condensate in a two-dimensional elongated trap is analyzed in
presence of rotation. Here, the ground state density splits into two clouds
for fast rotation. It is fascinating to observe the effect of rotation on the
condensate as it leads to a transition from condensed state to a fully two-
fold fragmented state in a single well without ramping up a barrier. Thus,
rotation can be used as a probe to manipulate various degrees of fragmentation
which we will further discuss in the next section where higher-order
fragmentations are explored. Further, we observe the presence of an opposite
anisotropy both in the position and momentum variances by comparing the mean-
field and many-body results for fast rotations.

###  III.2 Breaking up to several clouds and pathway to higher-order
fragmentation

So far in the elongated trap, we observe breaking up of the ground state
density into two clouds, followed by the emergence of two-fold fragmentation
and presence of opposite anisotropy both in the position and momentum
variances for faster rotations. Now, it would be fascinating to explore some
more complex potentials with $k$-fold rotational symmetry, explicitly three-
fold and four-fold symmetric potentials, to analyze the impact of rotation on
generic properties. Would the density break into more than two clouds? Is
higher-order fragmentation possible? How would the angular momentum enter into
the condensate? Finally, when correlations set in, who wins, the many-body or
the mean-field variance? Both for the position and momentum variances?

####  III.2.1 Bosons in a three-fold symmetric trap under rotation

In this section, we analyze the impact of rotation on weakly interacting
bosons confined in a three-fold symmetric trap of the form

$\displaystyle
V({\bf{r}})=\frac{1}{5}(x^{2}+y^{2})^{2}+\frac{1}{5}(x^{2}y-\frac{1}{3}y^{3}).$
(16)

In this trap, we consider $M=3$ self-consistent orbitals to compute the ground
state properties of the system. We also computed the following results with
$M=6$ self-consistent orbitals for convergence and checked the consistency of
our results (see the supplemental material). The rotation frequency range is
$\omega_{r}=[0,2.5]$ for the following analysis.

Figure 6: The one-body densities per particle are shown for three different
rotation frequencies $\omega_{r}$ in the three-fold symmetric trap at the
many-body level. M=3 self-consistent orbitals are used. The density computed
at the mean-field (not shown) and many-body levels depict identical features
for all $\omega_{r}$. All quantities shown are dimensionless.

We computed the ground-state energy $E/N$, and found it to display similar
pattern as that for the elongated trap. That is for slow rotation, E/N remains
almost constant and then E/N drops gradually with increase in $\omega_{r}$. In
addition, we found that the mean-field and many-body energies are practically
identical for all $\omega_{r}$. The results are shown in the supplemental
material, see Figure S1(b).

Figure 7: Pathway from condensation to fragmentation in the rotating three-
fold symmetric trap. three-fold fragmentation is observed with increase in
rotation. The three leading natural occupations $n_{1}/N$, $n_{2}/N$ and
$n_{3}/N$, are shown as a function of the rotation frequency $\omega_{r}$.
$M=3$ self-consistent orbitals are used. All quantities shown are
dimensionless.

Figure (6) displays the ground-state densities per particle of a rotating
condensate for three different rotation frequencies $\omega_{r}$, confined in
the three-fold symmetric trap [Equation (35)] computed at the many-body level.
In the absence of rotation, the one-body density depicts a single cloud where
all the bosons accumulate in the center of the trap. However, faster rotations
lead to breakup of the density profile into three clouds. We have also
computed the ground state density at the mean-field level and it shows
essentially identical features as the many-body level density in real space.

To further understand the breaking of the ground-state density, we discuss the
behaviour of the natural occupations $\frac{n_{j}}{N}$ as a function of the
rotation frequency $\omega_{r}$ as shown in Figure 7. It is evident from
Figure (7) that the system remains fully condensed, i.e., $\frac{n_{1}}{N}\sim
1,\frac{n_{2}}{N}\sim\frac{n_{3}}{N}\sim 10^{-6}$ with the inclusion of
rotation till about $\omega_{r}=1.1$. From about $\omega_{r}=1.2$ onwards, the
condensate starts to deplete with gradual decrease in the population of the
first natural orbital followed by corresponding increase in the populations of
the second and third natural orbitals, $\frac{n_{2}}{N}\sim\frac{n_{3}}{N}\sim
10^{-4}$ [see Figure S2(b) in the supplemental material]. For a faster
rotation, say $\omega_{r}=2$, the system transits to a fragmented state with
decrease in the population of the first natural orbital followed by
macroscopic population of the second and third natural orbitals. Further
increase in the rotation to $\omega_{r}=2.5$ leads to a three-fold fragmented
state having the natural occupations $\frac{n_{1}}{N}\sim 60\%$ and
$\frac{n_{2}}{N}\sim\frac{n_{3}}{N}\sim 20\%$ of the first, second, and third
natural orbitals, respectively. It is observed that to achieve a fully three-
fold fragmentation in this three-fold symmetric trap, we would need a much
faster rotation in comparison with the elongated trap. Therefore, in case of a
weakly interacting bosons confined in the three-fold symmetric trap given by
Equation (35), switching on rotation leads to a transition from a fully
condensed state to three-fold fragmentation.

Figure 8: Expectation value of angular momentum operator
$\langle\hat{L}_{Z}\rangle/N$, computed at the mean-field (MF) and many-body
(MB) levels [with $M=1$ and $M=3$ self-consistent orbitals, respectively] as a
function of the rotation frequency $\omega_{r}$ for the three-fold symmetric
trap. All quantities shown are dimensionless.

Let us move to the behaviour of average angular-momentum per particle
$\langle\hat{L}_{Z}\rangle/N$, as a function of the rotation frequency
$\omega_{r}$ for the three-fold symmetric trap computed both at the mean-field
and many-body levels. It is observed from Figure (8) that the mean-field and
many-body angular momenta exactly match each other for all $\omega_{r}$. The
angular momentum remains minimum till the rotation frequency of about
$\omega_{r}=1.1$. For about $\omega_{r}=1.2$, the rotation produces a state
where significant value of angular momentum generates with
$\langle\hat{L}_{Z}\rangle/N>1$ at which the breakup of the density is
observed. The angular momentum gradually increases with further increase in
the rotation of the condensate. Even for strong rotation, the angular momentum
computed at the mean-field and many-body levels coincide each other.
Therefore, we can conclude that, at least for the ground state, the angular
momentum and its variance [see Figure S6(b) in the supplemental material] are
not good indicators for many-body effects. Nonetheless, we stress that the
acquisition of angular momentum in the condensate follows the breakup of the
density, the emergence of depletion, and the eventual fragmentation in the
three-fold symmetric trap.

Figure 9: Dependence of the the many-particle position and momentum variances
on the rotation in the three-fold symmetric trap. Shown are (a)
$\frac{1}{N}\Delta^{2}_{\hat{X}}$ and $\frac{1}{N}\Delta^{2}_{\hat{Y}}$ and
(b) $\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}$ and
$\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}$ as a function of $\omega_{r}$ at the
many-body level (MB) [$M=3$ self-consistent orbitals] and at the mean-field
level (MF) [$M=1$ self-consistent orbital] All quantities shown are
dimensionless.

Now we move to the discussion of the impact of rotation on the behaviour of
the many-particle variances of the position and momentum operators along the
$x$\- and $y$-directions, to further characterize the many-body properties of
the rotating condensate undergoing breakup.

Figure 9(a) displays the behaviour of the many-particle position variance per
particle $\frac{1}{N}\Delta^{2}_{\hat{X},\hat{Y}}$ as a function of the
rotation frequency $\omega_{r}$ in the three-fold symmetric trap computed at
the mean-field and many-body levels. It is observed that the mean-field and
many-body $\frac{1}{N}\Delta^{2}_{\hat{X},\hat{Y}}$ coincide till about
$\omega_{r}=1.1$. For faster rotations, from about $\omega_{r}=1.2$ onwards
the mean-field and many-body position variances deviate. The mean-field
$\frac{1}{N}\Delta^{2}_{\hat{X},\hat{Y}}$ increase monotonously depicting the
spreading and finally the breakup of the density. However, the many-body
$\frac{1}{N}\Delta^{2}_{\hat{X},\hat{Y}}$ first increase and then, after
reaching a maximal value starts decreasing with further increase in
$\omega_{r}$, because of the depletion and eventual fragmentation. One of the
important features of the position variance is that the variance along the
$x$\- and $y$-direction exactly coincide each other for both the mean-field
and many-body regimes, thereby indicating the fact that the three-fold
rotational symmetry of the condensate is preserved even for fast rotations.
Clearly, there is no anisotropy of the variance for all $\omega_{r}$ in this
trap.

Finally, Figure 9(b) shows the behaviour of the momentum variance per particle
$\frac{1}{N}\Delta^{2}_{\hat{P}_{X},\hat{P}_{Y}}$, computed at the mean-field
and many-body levels as a function of $\omega_{r}$. The momentum variances
remain small even with increase in rotation till about $\omega_{r}=1.1$. From
about $\omega_{r}=1.2$ onwards, the momentum variances computed both at the
mean-field and many-body levels slowly start to increase. Here, the mean-field
momentum variance only increases monotonously. The many-body momentum variance
increases, decreases, and again increases but remain much smaller than the
mean-field momentum variance, indicates the presence of depletion and
fragmentation. Similar to the position variance, the momentum variance along
the $x$\- and $y$-directions exactly coincide each other both at the mean-
field and many-body levels, due to the rotational symmetry.

####  III.2.2 Bosons in a four-fold symmetric trap under rotation

Figure 10: The one-body densities per particle are shown for three different
rotation frequencies $\omega_{r}$ in the four-fold symmetric trap at the many-
body level. M=4 self-consistent orbitals are used. The density computed at the
mean-field (not shown) and many-body levels display identical features at all
$\omega_{r}$. All quantities shown are dimensionless.

Now we move to a more complicated system, a four-fold symmetric trap, to show
the stability of the ground-state properties found above for the three-fold
symmetric trap. The potential of the four-fold symmetric trap is given by

$\displaystyle V(\mathbf{r})=\frac{1}{4}(x^{4}+y^{4}).$ (17)

Here, we consider $M=4$ self-consistent orbitals to obtain the ground-state
properties of bosons under rotation in this trap. We also compute the results
with $M=8$ self-consistent orbitals to verify the numerical convergence, see
the supplemental material. The range of rotation frequencies is taken to be
$\omega_{r}=[0,2.0]$ for the following study.

Figure 11: Rotation leads to four-fold fragmentation in the four-fold
symmetric trap. $M=4$ self-consistent orbitals are used. The variation of four
natural occupations $n_{1}/N$, $n_{2}/N$, $n_{3}/N$ and $n_{4}/N$ are shown as
a function of rotation frequency $\omega_{r}$. All quantities computed are
dimensionless.

We computed the ground state energy $E/N$, and found it to display similar
pattern as that for the elongated and three-fold symmetric traps. That is for
slow rotation, E/N remains almost constant and then E/N drops gradually with
increase in $\omega_{r}$. In addition, we found that the mean-field and many-
body energies are practically identical for all $\omega_{r}$. The results are
shown in the supplemental material, see Figure S1(c).

Figure 10 shows the behaviour of the ground-state densities per particle of a
rotating BEC confined in the four-fold symmetric trap for three different
rotation frequencies $\omega_{r}$ at the many-body level. Similar to the
elongated and three-fold symmetric traps, for slow rotation the density
displays a single cloud. With increasing rotation, a deep in the density
emerges [Figure 10(b)] and finally, faster rotations lead to splitting of the
density into four sub-clouds as evident from Figure 10(c).

For a deeper understanding of the many-boson density profile, the behaviour of
the natural occupations $\frac{n_{j}}{N}$ as a function of the rotation
frequency $\omega_{r}$ is analyzed in Figure (11).

Figure 12: Expectation value of angular momentum operator
$\langle\hat{L}_{Z}\rangle/N$, computed at the mean-field (MF) and many-body
(MB) levels [with $M=1$ and $M=4$ self-consistent orbitals, respectively] as a
function of the rotation frequency $\omega_{r}$ for the four-fold symmetric
trap. All quantities shown are dimensionless.

It is found that the system preserves the fully condensed state, with
$\frac{n_{1}}{N}\sim
1,\frac{n_{2}}{N}\sim\frac{n_{3}}{N}\sim\frac{n_{4}}{N}\leq 10^{-6}$, till a
rotation frequency of about $\omega_{r}=1.0$. Further increase in $\omega_{r}$
leads to slow depletion of the condensate. For a faster rotation,
$\omega_{r}=2.0$, four-fold fragmentation of the condensate with finite
population of all the four natural orbitals is observed.

Let us discuss the behaviour of the average angular momentum per particle as a
function of $\omega_{r}$, see Figure 12. The angular momentum remains minimal
for $\omega_{r}=1$. From about $\omega_{r}=1.1$, the rotation produces a state
where significant value of angular momentum generates with
$\langle\hat{L}_{Z}\rangle/N>1$ at which the breakup of the density is
observed. The angular momentum gradually increases with further increase in
the rotation. Even for strong rotation, the angular momenta computed at the
mean-field and many-body levels coincide each other. Finally, it can be
concluded for the four-fold trap as well that, at least for the ground state,
the angular momentum and its variance [see Figure S6(C) in the supplemental
material] do not precisely signifies many-body effects. However, we can
conclude that accumulation of angular momentum in the condensate, the breakup
of the density, emergence of depletion, and eventual fragmentation in the
four-fold symmetric trap are in sync similar to the elongated and three-fold
symmetric traps.

Figure 13: Dependence of the the many-particle position and momentum variances
on the rotation in the four-fold symmetric trap. (a) depicts the position
variances $\frac{1}{N}\Delta^{2}_{\hat{X}}$, $\frac{1}{N}\Delta^{2}_{\hat{Y}}$
and (b) displays momentum variances $\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}$,
$\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}$ at the many-body level (MB) [$M=4$ self-
consistent orbitals] and at the mean-field level (MF) [$M=1$ self-consistent
orbital]. All quantities shown are dimensionless.

Figures 13 displays the behaviour of the position and momentum variances
computed at the mean-field and many-body levels as a function of rotation
frequency $\omega_{r}$ for four-fold symmetric trap along the $x$\- and
$y$-directions. The mean-field and many-body position variances
$\frac{1}{N}\Delta^{2}_{\hat{X},\hat{Y}}$ coincide till about $\omega_{r}=1$.
For faster rotations, from about $\omega_{r}=1.1$ onward, the mean-field
position variance deviates from the many-body position variance and the former
is always larger than the latter. Further, the position variances along the
$x$\- and $y$-directions exactly coincide with each other both at the mean-
field and many-body levels similar to the three-fold symmetric trap, see
Figure 9(a). Side by side, the behaviour of the momentum variance
$\frac{1}{N}\Delta^{2}_{\hat{P}_{X},\hat{P}_{Y}}$ displays similar feature as
that of the position variance, both at the mean-field and many-body levels. In
particular, mean-field is larger than many-body from about $\omega_{r}=1.5$.

Finally, we can conclude that the depletion, angular momentum, and the
position and momentum variances follow hand in hand both in the three-fold and
four-fold symmetric traps. The rotating interacting bosons acquire unique
many-body properties while undergoing breakup in space.

##  IV Concluding remarks

In the present work, we have studied the impact of rotation on the ground
state of weakly interacting bosonic atoms confined in two-dimensional
anharmonic potentials, first in, an elongated trap and, subsequently, in
three-fold and four-fold symmetric traps. Here, the multiconfigurational time-
dependent Hartree method for bosons, which is particularly suitable to
describe many-body properties, is employed to investigate the ground-state
energy, density, the depletion and fragmentation, angular momentum, and
finally, many-particle variances as a function of the rotation frequency to
characterize the correlations present in the system.

In the elongated trap, it is observed that the ground-state density breaks up
into two clouds with rotation. The splitting of the density is followed by the
emergence of the two-fold fragmentation. Interestingly, the ground state
exhibits opposite anistoropy both for the many-particle position and momentum
variances when computed at the many-body and mean-field levels. The rotation
squeezes the position variance in the elongated direction and the momentum
variance in the narrow direction, thereby producing unique correlations.
Finally, a synchronized pattern among the density breakup, eventual
fragmentation, acquisition of angular momentum in the condensate, and many-
particle variances is observed. When the angular momentum sets in the breakup
of density and the many-particle position and momentum variances start to
increase.

For the three-fold and four-fold symmetric traps, the ground-state density
eventually splits into three and four clouds, respectively, with the inclusion
of the rotation. Side-by-side, the rotation leads to transition from condensed
to three-fold and four-fold fragmented condensates, respectively, at the many-
body level of theory. We find that the depletion, the accumulation of angular
momentum, and finally, the increase in the variances of position and momentum
follow hand in hand.

The rotating frame in our work can be viewed as a specific case of synthetic
gauge fields. In future, continuation of this investigation includes the
extension of the study in various synthetic gauge fields. It would be
interesting to explore the ground-state breakup, eventual condensation or
fragmentation, and finally, various correlations of the condensate in presence
of synthetic gauge fields. It would be fascinating to discover many-body
features beyond the capacity of a “simple” rotation.

##  V Acknowledgements

This work is supported by the Israel Science Foundation (ISF) grant no.
1516/19 and by the Austrian Science Foundation (FWF) under grant P-32033-N32.
SD acknowledges Anal Bhowmik for helpful discussions. Computation time at the
High-Performance Computing Center Stuttgart (HLRS) is gratefully acknowledged.

## References

  * Madison  _et al._ [2000] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Vortex formation in a stirred Bose-Einstein condensate, Phys. Rev. Lett. 84, 806 (2000).
  * Raman  _et al._ [2001] C. Raman, J. Abo-Shaeer, J. Vogels, K. Xu, and W. Ketterle, Vortex nucleation in a stirred Bose-Einstein condensate, Phys. Rev. Lett. 87, 210402 (2001).
  * Matthews  _et al._ [1999] M. R. Matthews, B. P. Anderson, P. Haljan, D. Hall, C. Wieman, and E. A. Cornell, Vortices in a Bose-Einstein condensate, Phys. Rev. Lett. 83, 2498 (1999).
  * Haljan  _et al._ [2001] P. Haljan, B. Anderson, I. Coddington, and E. A. Cornell, Use of surface-wave spectroscopy to characterize tilt modes of a vortex in a Bose-Einstein condensate, Phys. Rev. Lett. 86, 2922 (2001).
  * Abo-Shaeer  _et al._ [2002] J. Abo-Shaeer, C. Raman, and W. Ketterle, Formation and decay of vortex lattices in Bose-Einstein condensates at finite temperatures, Phys. Rev. Lett. 88, 070409 (2002).
  * Butts and Rokhsar [1999] D. Butts and D. Rokhsar, Predicted signatures of rotating Bose-Einstein condensates, Nature (London) 397, 327 (1999).
  * Dagnino  _et al._ [2009] D. Dagnino, N. Barberán, M. Lewenstein, and J. Dalibard, Vortex nucleation as a case study of symmetry breaking in quantum systems, Nat. Phys. 5, 431 (2009).
  * Coddington  _et al._ [2003] I. Coddington, P. Engels, V. Schweikhard, and E. A. Cornell, Observation of Tkachenko oscillations in rapidly rotating Bose-Einstein condensates, Phys. Rev. Lett. 91, 100402 (2003).
  * Schweikhard  _et al._ [2004] V. Schweikhard, I. Coddington, P. Engels, V. P. Mogendorff, and E. A. Cornell, Rapidly rotating Bose-Einstein condensates in and near the lowest Landau level, Phys. Rev. Lett. 92, 040404 (2004).
  * Regnault and Jolicoeur [2003] N. Regnault and T. Jolicoeur, Quantum Hall fractions in rotating Bose-Einstein condensates, Phys. Rev. Lett. 91, 030402 (2003).
  * Chiappini  _et al._ [2011] 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 (London) 472, 454 (2011).
  * Bohr and Mottelson [1979] A. Bohr and B. R. Mottelson, Physics of rapidly rotating nuclei, Physics Today 32, 25 (1979).
  * Guo  _et al._ [2020] Y. Guo, R. Dubessy, M. d. 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).
  * Kumar  _et al._ [2016] R. K. Kumar, T. Sriraman, H. Fabrelli, P. Muruganandam, and A. Gammal, Three-dimensional vortex structures in a rotating dipolar Bose-Einstein condensate, J. Phys. B: At., Mol. Opt. Phys. 49, 155301 (2016).
  * Brito  _et al._ [2020] L. Brito, A. Andriati, L. Tomio, and A. Gammal, Breakup of rotating asymmetric quartic-quadratic trapped condensates, Phys. Rev. A 102, 063330 (2020).
  * Aftalion [2007] A. Aftalion,  _Vortices in Bose-Einstein Condensates_ , Vol. 67 (Springer Science & Business Media, 2007).
  * Fetter [2009] A. L. Fetter, Rotating trapped Bose-Einstein condensates, Rev. Mod. Phys. 81, 647 (2009).
  * Tsatsos  _et al._ [2016] M. C. Tsatsos, P. E. Tavares, A. Cidrim, A. R. Fritsch, M. A. Caracanhas, F. E. A. dos Santos, C. F. Barenghi, and V. S. Bagnato, Quantum turbulence in trapped atomic Bose-Einstein condensates, Phys. Rep. 622, 1 (2016).
  * Viefers  _et al._ [2000] S. Viefers, T. Hansson, and S. Reimann, Bose condensates at high angular momenta, Phys. Rev. A 62, 053604 (2000).
  * Reimann  _et al._ [2006] S. Reimann, M. Koskinen, Y. Yu, and M. Manninen, Rotating quantum liquids crystallize, New J. Phys. 8, 59 (2006).
  * Cremon  _et al._ [2013] J. Cremon, G. Kavoulakis, B. Mottelson, and S. Reimann, Vortices in Bose-Einstein condensates: Finite-size effects and the thermodynamic limit, Phys. Rev. A 87, 053615 (2013).
  * Cremon  _et al._ [2015] J. Cremon, A. Jackson, E. Karabulut, G. Kavoulakis, B. Mottelson, and S. Reimann, Rotating Bose-Einstein condensates: Closing the gap between exact and mean-field solutions, Phys. Rev. A 91, 033623 (2015).
  * Beinke  _et al._ [2018] R. Beinke, L. S. Cederbaum, and O. E. Alon, Enhanced many-body effects in the excitation spectrum of a weakly interacting rotating Bose-Einstein condensate, Phys. Rev. A 98, 053634 (2018).
  * Khanore and Dey [2022] M. P. Khanore and B. Dey, The quantum vortex states in extended Bose-Hubbard model: effects of lattice geometries, inter-particle interactions and spatial inhomogeneity, Eur. Phys. J. D 76, 1 (2022).
  * Löwdin [1955] P.-O. Löwdin, Quantum Theory of Many-Particle Systems. I. Physical Interpretations by Means of Density Matrices, Natural Spin-Orbitals, and Convergence Problems in the Method of Configurational Interaction, Phys. Rev. 97, 1474 (1955).
  * Yang [1962] C. N. Yang, Concept of Off-Diagonal Long-Range Order and the Quantum Phases of Liquid He and of Superconductors, Rev. Mod. Phys. 34, 694 (1962).
  * Davidson [1976] E. Davidson,  _Reduced Density Matrices in Quantum Chemistry_ , Vol. 6 (Academic Press: New York, NY, USA, 1976).
  * Coleman and Yukalov [2000] A. J. Coleman and V. I. Yukalov,  _Reduced Density Matrices: Coulson’s Challenge_ , Vol. 72 (Lectures Notes in Chemistry; Springer: Berlin, Germany, 2000).
  * Fischer  _et al._ [2015] U. R. Fischer, A. U. Lode, and B. Chatterjee, Condensate fragmentation as a sensitive measure of the quantum many-body behavior of bosons with long-range interactions, Phys. Rev. A 91, 063621 (2015).
  * Lode [2016] A. U. Lode, Multiconfigurational time-dependent Hartree method for bosons with internal degrees of freedom: Theory and composite fragmentation of multicomponent Bose-Einstein condensates, Phys. Rev. A 93, 063601 (2016).
  * Penrose and Onsager [1956] O. Penrose and L. Onsager, Bose-Einstein Condensation and Liquid Helium, Phys. Rev. 104, 576 (1956).
  * Mueller  _et al._ [2006] E. J. Mueller, T.-L. Ho, M. Ueda, and G. Baym, Fragmentation of Bose-Einstein condensates, Phys. Rev. A 74, 033612 (2006).
  * Girardeau [1962] M. Girardeau, Simple and Generalized Condensation in Many-Boson Systems, Phys. Fluids 5, 1468 (1962).
  * Spekkens and Sipe [1999] R. Spekkens and J. Sipe, Spatial fragmentation of a Bose-Einstein condensate in a double-well potential, Phys. Rev. A 59, 3868 (1999).
  * Streltsov  _et al._ [2004] A. I. Streltsov, L. S. Cederbaum, and N. Moiseyev, Ground-state fragmentation of repulsive Bose-Einstein condensates in double-trap potentials, Phys. Rev. A 70, 053607 (2004).
  * Streltsov and Cederbaum [2005] A. I. Streltsov and L. S. Cederbaum, Properties of fragmented repulsive condensates, Phys. Rev. A 71, 063612 (2005).
  * Alon and Cederbaum [2005] O. E. Alon and L. S. Cederbaum, Pathway from condensation via fragmentation to fermionization of cold bosonic systems, Phys. Rev. Lett. 95, 140402 (2005).
  * Bader and Fischer [2009] P. Bader and U. R. Fischer, Fragmented many-body ground states for scalar bosons in a single trap, Phys. Rev. Lett. 103, 060402 (2009).
  * Fischer and Bader [2010] U. R. Fischer and P. Bader, Interacting trapped bosons yield fragmented condensate states in low dimensions, Phys. Rev. A 82, 013607 (2010).
  * Zhou and Cui [2013] Q. Zhou and X. Cui, Fate of a Bose-Einstein condensate in the presence of spin-orbit coupling, Phys. Rev. Lett. 110, 140407 (2013).
  * Kawaguchi [2014] Y. Kawaguchi, Goldstone-mode instability leading to fragmentation in a spinor Bose-Einstein condensate, Phys. Rev. A 89, 033627 (2014).
  * Song  _et al._ [2014] S.-W. Song, Y.-C. Zhang, H. Zhao, X. Wang, and W.-M. Liu, Fragmentation of spin-orbit-coupled spinor Bose-Einstein condensates, Phys. Rev. A 89, 063613 (2014).
  * Kang and Fischer [2014] M.-K. Kang and U. R. Fischer, Revealing single-trap condensate fragmentation by measuring density-density correlations after time of flight, Phys. Rev. Lett. 113, 140404 (2014).
  * Jen and Yip [2015] H.-H. Jen and S.-K. Yip, Fragmented many-body states of a spin-2 bose gas, Phys. Rev. A 91, 063603 (2015).
  * Fischer and Kang [2015] U. R. Fischer and M.-K. Kang, “photonic” cat states from strongly interacting matter waves, Phys. Rev. Lett. 115, 260404 (2015).
  * Kolovsky [2017] A. R. Kolovsky, Bogoliubov depletion of the fragmented condensate in the bosonic flux ladder, Phys. Rev. A 95, 033622 (2017).
  * Tomchenko [2020] M. Tomchenko, On a fragmented condensate in a uniform bose system, J. Low Temp. Phys. 198, 100 (2020).
  * Sakmann and Kasevich [2016] K. Sakmann and M. Kasevich, Single-shot simulations of dynamic quantum many-body systems, Nat. Phys. 12, 451 (2016).
  * Sánchez-Lotero and Palacios [2005] P. Sánchez-Lotero and J. J. Palacios, Vortices in a rotating Bose-Einstein condensate under extreme elongation, Phys. Rev. A 72, 043613 (2005).
  * Williams and Holland [1999] J. Williams and M. Holland, Preparing topological states of a Bose-Einstein condensate, Nature (London) 401, 568 (1999).
  * Beinke  _et al._ [2015] R. Beinke, S. Klaiman, L. S. Cederbaum, A. I. Streltsov, and O. E. Alon, Many-body tunneling dynamics of Bose-Einstein condensates and vortex states in two spatial dimensions, Phys. Rev. A 92, 043627 (2015).
  * Alon [2019a] O. E. Alon, Analysis of a trapped Bose-Einstein condensate in terms of position, momentum, and angular-momentum variance, Symmetry 11, 1344 (2019a).
  * Gross [1961] E. P. Gross, Structure of a quantized vortex in boson systems, Il Nuovo Cimento 20, 454 (1961).
  * Pitaevskii [1961] L. P. Pitaevskii, Vortex lines in an imperfect Bose gas, Sov. Phys. JETP 13, 451 (1961).
  * Streltsov  _et al._ [2007] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Role of excited states in the splitting of a trapped interacting Bose-Einstein condensate by a time-dependent barrier, Phys. Rev. Lett. 99, 030402 (2007).
  * Alon  _et al._ [2008] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Multiconfigurational time-dependent Hartree method for bosons: Many-body dynamics of bosonic systems, Phys. Rev. A 77, 033613 (2008).
  * Beck  _et al._ [2000] M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, The multiconfiguration time-dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wavepackets, Phys. Rep. 324, 1 (2000).
  * Wang and Thoss [2003] H. Wang and M. Thoss, Multilayer formulation of the multiconfiguration time-dependent Hartree theory, J. Chem. Phys. 119, 1289 (2003).
  * Manthe [2008] U. Manthe, A multilayer multiconfigurational time-dependent Hartree approach for quantum dynamics on general potential energy surfaces, J. Chem. Phys. 128, 164116 (2008).
  * Wang [2015] H. Wang, Multilayer Multiconfiguration Time-Dependent Hartree Theory, J. Phys. Chem. A 119, 7951 (2015).
  * Manthe and Weike [2017] U. Manthe and T. Weike, On the multi-layer multi-configurational time-dependent Hartree approach for bosons and fermions, J. Chem. Phys. 146, 064117 (2017).
  * Manthe [2017] U. Manthe, Wavepacket dynamics and the multi-configurational time-dependent Hartree approach, J. Phys.: Condens. Matter 29, 253001 (2017).
  * Bhowmik  _et al._ [2020] A. Bhowmik, S. K. Haldar, and O. E. Alon, Impact of the transverse direction on the many-body tunneling dynamics in a two-dimensional bosonic Josephson junction, Sci. Rep. 10, 1 (2020).
  * Bhowmik and Alon [2022] A. Bhowmik and O. E. Alon, Longitudinal and transversal resonant tunneling of interacting bosons in a two-dimensional josephson junction, Sci. Rep. 12, 1 (2022).
  * Zanghellini  _et al._ [2003] J. Zanghellini, M. Kitzler, C. Fabian, T. Brabec, and A. Scrinzi, An MCTDHF Approach to Multielectron Dynamics in Laser Fields, Laser Phys. 13, 1064 (2003).
  * Cao  _et al._ [2013] L. Cao, S. Krönke, O. Vendrell, and P. Schmelcher, The multi-layer multi-configuration time-dependent Hartree method for bosons: Theory, implementation, and applications, J. Chem. Phys. 139, 134103 (2013).
  * Cao  _et al._ [2017] L. Cao, V. Bolsinger, S. I. Mistakidis, G. M. Koutentakis, S. Krönke, J. M. Schurer, and P. Schmelcher, A unified ab initio approach to the correlated quantum dynamics of ultracold fermionic and bosonic mixtures, J. Chem. Phys. 147, 044106 (2017).
  * Haxton and McCurdy [2015] D. J. Haxton and C. W. McCurdy, Two methods for restricted configuration spaces within the multiconfiguration time-dependent Hartree-Fock method, Phys. Rev. A 91, 012509 (2015).
  * Miyagi and Madsen [2017] H. Miyagi and L. B. Madsen, Time-dependent restricted-active-space self-consistent-field theory with space partition, Phys. Rev. A 95, 023415 (2017).
  * Lévêque and Madsen [2017] C. Lévêque and L. B. Madsen, Time-dependent restricted-active-space self-consistent-field theory for bosonic many-body systems, New J. Phys. 19, 043007 (2017).
  * Lévêque and Madsen [2018] C. Lévêque and L. B. Madsen, Multispecies time-dependent restricted-active-space self-consistent-field theory for ultracold atomic and molecular gases, J. Phys. B: At., Mol. Opt. Phys. 51, 155302 (2018).
  * Lode  _et al._ [2020] A. U. J. Lode, C. Lévêque, L. B. Madsen, A. I. Streltsov, and O. E. Alon, Colloquium : Multiconfigurational time-dependent Hartree approaches for indistinguishable particles, Rev. Mod. Phys. 92, 011001 (2020).
  * Lode  _et al._ [2019] A. U. J. Lode, M. Tsatsos, E. Fasshauer, R. Lin, L. Papariello, P. Molignini, C. Lévêque, and S. Weiner, MCTDH-X: The time-dependent multiconfigurational Hartree for indistinguishable particles software (2019).
  * Lin  _et al._ [2020] R. Lin, P. Molignini, L. Papariello, M. C. Tsatsos, C. Lévêque, S. E. Weiner, E. Fasshauer, R. Chitra, and A. U. J. Lode, MCTDH-X: The multiconfigurational time-dependent Hartree method for indistinguishable particles software, Quantum Sci. Technol. 5, 024004 (2020).
  * Lode  _et al._ [2021] A. U. J. Lode, O. E. Alon, M. A. Bastarrachea-Magnani, A. Bhowmik, A. Buchleitner, L. Cederbaum, R. Chitra, E. Fasshauer, L. d. Forges de Parny, S. Haldar, _et al._ , MCTDH-X: The Multiconfigurational Time-Dependent Hartree Method for Indistinguishable Particles High-Performance Computation Project, in  _High Performance Computing in Science and Engineering’20_ (Springer, 2021) pp. 21–45.
  * Lode  _et al._ [2012] A. U. Lode, K. Sakmann, O. E. Alon, L. S. Cederbaum, and A. I. Streltsov, Numerically exact quantum dynamics of bosons with time-dependent interactions of harmonic type, Phys. Rev. A 86, 063606 (2012).
  * Doganov  _et al._ [2013] R. A. Doganov, S. Klaiman, O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Two trapped particles interacting by a finite-range two-body potential in two spatial dimensions, Physical Review A 87, 033631 (2013).
  * Dirac [1930] P. A. M. Dirac, Note on Exchange Phenomena in the Thomas Atom, Math. Proc. Cambridge Philos. Soc. 26, 376 (1930).
  * Mott and Frenkel [1934] N. F. Mott and J. Frenkel,  _Math. Gaz._ , Vol. 18 (Clarendon Press, Oxford, 1934) p. 208.
  * Kramer and Saraceno [2007] P. Kramer and M. Saraceno,  _Geometry of the time-dependent variational principle in quantum mechanics_ (Springer, Lecture Notes in Physics, 2007).
  * Weiner  _et al._ [2017] S. E. Weiner, M. C. Tsatsos, L. S. Cederbaum, and A. U. Lode, Phantom vortices: Hidden angular momentum in ultracold dilute Bose-Einstein condensates, Sci. Rep. 7, 40122 (2017).
  * Tsatsos  _et al._ [2010] M. C. Tsatsos, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Fragmented many-body states of definite angular momentum and stability of attractive three-dimensional condensates, Phys. Rev. A 82, 33613 (2010).
  * Klaiman and Alon [2015] S. Klaiman and O. E. Alon, Variance as a sensitive probe of correlations, Phys. Rev. A 91, 063613 (2015).
  * Alon [2019b] O. E. Alon, Variance of a Trapped Bose-Einstein Condensate, J. Phys.: Conf. Ser. 1206, 012009 (2019b).
  * Lode [2014] A. U. J. Lode,  _Tunneling Dynamics in Open Ultracold Bosonic Systems: Numerically Exact Dynamics–Analytical Models–Control Schemes_ (Springer, 2014).
  * Cohen and Lee [1985] L. Cohen and C. Lee, Exact reduced density matrices for a model problem, J. Math. Phys. 26, 3105 (1985).

## Supplemental material to fragmentation and correlations in a rotating Bose-
Einstein condensate undergoing breakup

This Supplemental material provides further many-body analysis that supports
our main results. Sec. S1 benchmarks the multiconfigurational time-dependent
Hartree for bosons (MCTDHB) method [76, 85, 72] in a harmonic interaction
model (HIM) under rotation that gives an exactly solvable many-body model.
This method is employed for the numerical computations of the main results.
Sec. S2 gives the convergence checks of the ground state energy per particle
of our main results with respect to the number of self-consistent orbitals in
different anharmonic traps, namely, elongated, three-fold symmetric, and four-
fold symmetric traps. Sec. S3 illustrates the convergence of natural
occupations in these confining traps with respect to the number of orbitals.
Sec. S4 presents the many-particle variances discussed in the main text, along
with their convergences with respect to the number of self-consistent
orbitals. This section is further divided into two subsections. The first
subsection deals with convergence of the position and momentum variances.
Convergence with respect to the number of self-consistent orbitals of the
expectation value of the angular momenta and angular momentum variance is
elaborated in the second subsection that supports our main results.

##  VI Benchmarking of MCTDHB in the rotating frame: An exactly Solvable many-
body model

In this section, we present a numerical benchmark of the MCTDHB method in the
rotating frame in two dimensions 2D. The many-body Hamiltonian in the rotating
frame can be written as

$\hat{H}(\mathbf{r_{1}},\mathbf{r_{2}},...,\mathbf{r_{N}})=\sum_{j=1}^{N}\hat{h}(\mathbf{r}_{i})+\sum_{j<k}^{N}\hat{W}(\mathbf{r}_{j},\mathbf{r}_{k})-\omega_{r}\hat{L}_{Z},$
(18)

where the single-particle Hamiltonian is
$\hat{h}(\mathbf{r})=-\frac{1}{2}\frac{\partial^{2}}{\partial\mathbf{r}^{2}}+\hat{V}(\mathbf{r})$.
We consider a harmonic interaction model (HIM) in 2D. In this HIM, the
confining trapping potential $\hat{V}(\mathbf{r})$ is a translated 2D harmonic
oscillator of the following form

$\hat{V}(\mathbf{r})=\frac{1}{2}m\omega^{2}[(x-L)^{2}+y^{2}],$ (19)

and the two-body interaction potential
$\hat{W}(\mathbf{r}_{j},\mathbf{r}_{k})$ is also harmonic

$\displaystyle\hat{W}(\mathbf{r}_{j},\mathbf{r}_{k})=\lambda_{0}(\mathbf{r}_{j}-\mathbf{r}_{k})^{2},$
(20)

where the two-body interaction strength $\lambda_{0}$ defines the mean-field
interaction parameter as $\Lambda=\lambda_{0}(N-1)$, $\omega_{r}$ is the
rotation frequency, and
$\hat{L}_{Z}=\sum_{j=1}^{N}\hat{l}_{z_{j}}=\sum^{N}_{j=1}\left(\hat{x}_{j}\hat{p}_{y_{j}}-\hat{y}_{j}\hat{p}_{x_{j}}\right)$
is the many-body angular momentum operator.

Now, the ground state solution of the HIM in the rotating frame is obtained by
solving the many-boson Schrödinger equation given as

$\hat{H}(\mathbf{r})\Psi(\mathbf{r})=E\Psi(\mathbf{r}).$ (21)

The $N$-boson wave function of the ground state in the rotating frame is
obtained as

$\Psi(\mathbf{r}_{1},\mathbf{r}_{2},...,\mathbf{r}_{N})=\bigg{(}\frac{m\Omega}{\pi}\bigg{)}^{\frac{N-1}{4}}\bigg{(}\frac{m\omega}{\pi}\bigg{)}^{\frac{1}{4}}e^{+i\beta\sum_{j=1}^{N}y_{j}}e^{-\frac{1}{2}m\omega\sum_{j=1}^{N}(x_{j}-x_{0})^{2}}e^{-\frac{1}{2}m\omega\sum_{j=1}^{N}y_{j}^{2}},$
(22)

where the dressed frequency is

$\Omega=\sqrt{\omega^{2}+\frac{2\lambda_{0}N}{m}},$ (23)

and the translation and momentum parameters are

$x_{0}=\frac{1}{1-(\frac{\omega_{r}}{\omega})^{2}}L$ (24)

and

$\beta=m\omega_{r}x_{0}=m\omega_{r}\frac{1}{1-(\frac{\omega_{r}}{\omega})^{2}}L.$
(25)

Thus, we can evaluate the ground state properties, such as, the energy,
densities, and reduced density matrices from the many-boson wavefunction
$\Psi$ in Equation (22).

As we know, the ground state energy per particle in the HIM without rotation
is given as [86]

$\frac{E(\omega_{r}=0)}{N}=\bigg{[}(N-1)\Omega+\omega\bigg{]}=\bigg{[}(N-1)\sqrt{\omega^{2}+\frac{2\lambda_{0}N}{m}}+\omega\bigg{]}.$
(26)

Finally, in the rotating frame the energy per particle is obtained as

$\frac{E(\omega_{r})}{N}=\frac{E(\omega_{r}=0)}{N}-\frac{1}{2}m\omega^{2}\frac{(\frac{\omega_{r}}{\omega})^{2}}{1-(\frac{\omega_{r}}{\omega})^{2}}L^{2}.$
(27)

Now we are going to evaluate the expectation values and variances of many-
particle observables, explicitly the position, momentum, and angular momentum
operators. The expectation values of the many-particle position operators are
given as

$\frac{\langle\hat{X}\rangle}{N}=x_{0}=\frac{1}{1-(\frac{\omega_{r}}{\omega})^{2}}L,~{}~{}~{}~{}~{}~{}~{}\frac{\langle\hat{Y}\rangle}{N}=0,$
(28)

and the expectation values of the many-particle momentum operators are given
as

$\frac{\langle\hat{P}_{X}\rangle}{N}=0,~{}~{}~{}~{}~{}~{}~{}\frac{\langle\hat{P}_{Y}\rangle}{N}=\beta=m\omega_{r}x_{0}=m\omega_{r}\frac{1}{1-(\frac{\omega_{r}}{\omega})^{2}}L.$
(29)

Similarly, the expectation value of the many-particle angular momentum
operator reads

$\frac{\langle\hat{L}_{Z}\rangle}{N}=x_{0}\beta=m\omega_{r}\Bigg{[}\frac{L}{1-(\frac{\omega_{r}}{\omega})^{2}}\Bigg{]}^{2}.$
(30)

The many-particle position variances per particle are obtained as,

$\frac{1}{N}\Delta^{2}_{\hat{X}}=\frac{1}{N}\Delta^{2}_{\hat{Y}}=\frac{1}{2m\omega},$
(31)

and the many-particle momentum variances per particle are the following

$\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}=\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}=\frac{m\omega}{2}.$
(32)

It is evident that the many-particle variances of the position and momentum
operators in the rotating frame boil down to those for the non-rotating case.

However, the many-particle angular momentum variance per particle is affected
by the rotation and takes on the following form

$\displaystyle\frac{1}{N}\Delta^{2}_{\hat{L}_{Z}}=x_{0}^{2}\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}+\beta^{2}\frac{1}{N}\Delta^{2}_{\hat{X}}=\frac{1}{2}m\omega\frac{\bigg{[}1+(\frac{\omega_{r}}{\omega})^{2}\bigg{]}}{\bigg{[}1-(\frac{\omega_{r}}{\omega})^{2}\bigg{]}^{2}}L^{2}.$
(33)

With a solvable many-boson model in the rotating frame, we can now proceed and
benchmark MCTDHB in the rotating frame. In the following simulations, we work
in the units $\hbar=m=\omega=1$. To compute the ground state of MCTDHB in the
rotating frame, we consider $N=10$ weakly interacting bosons with repulsive
interaction parameter $\Lambda=-0.1$. The box extensions of the confining
potential are chosen to be $\bf{=[-12,12)\times[-12,12)}$ with $\bf{128\times
128}$ DVR functions to represent each of the orbitals. We translated the
confining potential by $L=2$ in Equation (19).

In Table 1, we show the numerical convergence of the energy per particle
$\frac{E}{N}$. Exact analytical versus numerical MCTDHB ground state energies
are obtained for three different rotation frequencies $\omega_{r}$. The
numerical convergence of the many-body energy per particle is achieved by
using $M=6$ self-consistent orbitals. Thus, we conclude that for all
$\omega_{r}$, the ground state energy of the system is converged.

$\omega_{r}$ | ${\frac{E}{N}}_{numerical}$ | ${\frac{E}{N}}_{analytical}$
---|---|---
$0$ | 0.89372 | 0.89372
$0.5$ | 0.22705 | 0.22705
$0.8$ | -2.66183 | -2.66183
Table 1:  Benchmarking of the ground state many-body energy E/N with respect to the rotation frequency $\omega_{r}$. $M=6$ self-consistent orbitals are used. The ground state energies shown are in dimensionless units. $\omega_{r}$ | ${\frac{1}{N}\langle\hat{X}\rangle}_{numerical}$ | ${\frac{1}{N}\langle\hat{X}\rangle}_{analytical}$ | ${\frac{1}{N}\langle\hat{Y}\rangle}_{numerical}$ | ${\frac{1}{N}\langle\hat{Y}\rangle}_{analytical}$
---|---|---|---|---
0 | 2.00000 | 2.00000 | 0.00000 | 0.00000
0.5 | 2.66666 | 2.66666 | 0.00000 | 0.00000
0.8 | 5.55555 | 5.55555 | 0.00000 | 0.00000
$\omega_{r}$ | ${\frac{1}{N}\Delta^{2}_{\hat{X}}}_{numerical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{X}}}_{analytical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{Y}}}_{numerical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{Y}}}_{analytical}$
0 | 0.50000 | 0.50000 | 0.50000 | 0.50000
0.5 | 0.50000 | 0.50000 | 0.50000 | 0.50000
0.8 | 0.50000 | 0.50000 | 0.50000 | 0.50000
Table 2:  Benchmarking of the expectation values and variances of the many-
particle position operators along the $x$\- and $y$\- directions with respect
to the rotation frequency $\omega_{r}$ computed with $M=6$ self-consistent
orbitals. All the quantities are shown in dimensionless units.

Table 2 shows the comparison of the exact analytical and numerical MCTDHB
results of the many-particle position variances for various rotation
frequencies $\omega_{r}$. The numerical convergence of MCTDHB results computed
with $M=6$ self-consistent orbitals are clearly evident from the table.

$\omega_{r}$ | ${\frac{1}{N}\langle\hat{P}_{X}\rangle}_{numerical}$ | ${\frac{1}{N}\langle\hat{P}_{X}\rangle}_{analytical}$ | ${\frac{1}{N}\langle\hat{P}_{Y}\rangle}_{numerical}$ | ${\frac{1}{N}\langle\hat{P}_{Y}\rangle}_{analytical}$
---|---|---|---|---
0 | 0.00000 | 0.00000 | 0.00000 | 0.00000
0.5 | 0.00000 | 0.00000 | 1.33333 | 1.33333
0.8 | 0.00000 | 0.00000 | 4.44444 | 4.44444
$\omega_{r}$ | ${\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}}_{numerical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}}_{analytical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}}_{numerical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}}_{analytical}$
0 | 0.50000 | 0.50000 | 0.50000 | 0.50000
0.5 | 0.50000 | 0.50000 | 0.50000 | 0.50000
0.8 | 0.50000 | 0.50000 | 0.50000 | 0.50000
Table 3:  Benchmarking of the expectation values and variances of the many-particle momentum operators along the $x$\- and $y$\- directions with respect to the rotation frequency $\omega_{r}$ obtained with $M=6$ self-consistent orbitals. All the quantities are shown in dimensionless units. $\omega_{r}$ | ${\frac{1}{N}\langle\hat{L}_{Z}\rangle}_{numerical}$ | ${\frac{1}{N}\langle\hat{L}_{Z}\rangle}_{analytical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{L}_{Z}}}_{numerical}$ | ${\frac{1}{N}\Delta^{2}_{\hat{L}_{Z}}}_{analytical}$
---|---|---|---|---
0 | 0.00000 | 0.00000 | 2.00000 | 2.00000
0.5 | 3.55555 | 3.55555 | 4.44444 | 4.44444
0.8 | 24.6913 | 24.6913 | 25.3086 | 25.3086
Table 4:  Benchmarking of the expectation values and variances of the many-
particle angular momentum operator with respect to the rotation frequency
$\omega_{r}$ computed with $M=6$ self-consistent orbitals. All the quantities
are shown in dimensionless units.

Table 3 and 4 correspond to benchmarking of the many-particle momentum and
angular momentum variances for various rotation frequencies computed with
$M=6$ self-consistent orbitals. As the many-particle variances of the momentum
and angular momentum are very sensitive quantities, it is generally more
difficult to achieve numerical convergence of the analytical and numerical
results. However, from these tables, it is evident that the MCTDHB results are
converged with exact analytical results.

##  VII Convergence of the ground state energy

In this section, we investigate and report the convergence of the many-
particle ground state energy $\frac{E}{N}$ for the three confining potentials,

Figure 14: Many-body ground state energy per particle $\frac{E}{N}$ as a
function of the rotation frequency $\omega_{r}$. (a) for the elongated trap
$M=8$ self-consistent orbitals are used to show the convergence of the $M=4$
results, (b) for the three-fold symmetric trap $M=6$ self-consistent orbitals
are used to check the convergence of $M=3$ results, and finally in (c) $M=8$
self-consistent orbitals are used to check the convergence of $M=4$ results of
the four-fold symmetric trap. The inset shows the energy difference per
particle as defined in Equation (37). The quantities shown are dimensionless.

namely the elongated trap

$\displaystyle V(\mathbf{r})=\frac{1}{4}(0.8x^{2}+y^{2})^{2},$ (34)

the three-fold symmetric trap

$\displaystyle
V({\bf{r}})=\frac{1}{5}(x^{2}+y^{2})^{2}+\frac{1}{5}(x^{2}y-\frac{1}{3}y^{3}),$
(35)

and, finally, the four-fold symmetric trap defined as

$\displaystyle V(\mathbf{r})=\frac{1}{4}(x^{4}+y^{4}).$ (36)

The convergence of the many-body energies reported in the main text with $M$
self-consistent orbitals is here verified by taking $2M$ self-consistent
orbitals and comparing the respective results.

Figure 14 shows the behavior of the ground state energy $\frac{E}{N}$ as a
function of the rotation frequency $\omega_{r}$ for different numbers of self-
consistent orbitals for three confining potentials given by Equations (34),
(35), and (36). The inset in all the three panels of Figure 14 display the
energy difference between the energies computed for two different orbital
numbers $M$,

$\displaystyle\frac{\Delta E}{N}=\frac{E_{M}}{N}-\frac{E_{2M}}{N}.$ (37)

Figure 14(a) corresponds to the elongated trap [Equation (34)]. The
convergence of energy is checked for two different orbitals $M=4,8$. We find
that the two energy curves $\frac{E}{N}$ computed with $M=4,8$ orbitals fall
on top of each other for all the frequencies. It is evident that our results
converges with respect to the number of orbitals for elongated trap. Similar
convergence is observed in case of the three-fold symmetric trap with $M=3,6$
orbitals in Figure 14(b) [see Equation (35)] and the four-fold symmetric trap
[Equation (36)] with $M=4,8$ orbitals as evident from Figure 14(c). We
conclude that the ground state energies are fully converged with respect to
the number of orbitals in all the three confining potentials at all rotation
frequencies.

##  VIII Convergence of the natural occupations

The behavior of the natural occupations $n_{j}/N$ for the confining potentials
defined in Equations (34), (35), and (36), as a function of the rotation
frequency $\omega_{r}$ are shown in Figures 15 for different orbital numbers
$M$.

Figure 15:  Convergence of the natural occupations with respect to the number
of natural orbitals. (a) for the elongated trap with M=4,8 self-consistent
orbitals, (b) for the three-fold symmetric trap with M=3,6 self-consistent
orbitals, and (c) for the four-fold symmetric trap with M=4,8 self-consistent
orbitals. The inset in panels (b) and (c) depict the variation of the smallest
natural occupations with rotation.

It is observed that all the natural occupations are fully converged for the
three confining potentials. The inset in Figure 15(b) shows the variation of
the smallest natural occupations $\frac{n_{4}}{N}$, $\frac{n_{5}}{N}$, and
$\frac{n_{6}}{N}$ with the rotation frequency $\omega_{r}$ for the three-fold
symmetric trap. All the three natural occupations stays minimal with
$\frac{n_{4}}{N}\sim\frac{n_{5}}{N}\sim\frac{n_{6}}{N}\sim 10^{-6}$ for all
$\omega_{r}$. Similarly, the inset of Figure 15(c) depicts the natural
occupations $\frac{n_{5}}{N}\sim 10^{-6}$,
$\frac{n_{6}}{N}\sim\frac{n_{7}}{N}\sim\frac{n_{8}}{N}\sim 10^{-7}-10^{-6}$
for the four-fold symmetric trap. The convergence of the natural occupations
is evident in all three confining traps.

##  IX Convergence of the position, momentum, and angular momentum variances

The variance of a many-particle operator $\hat{A}=\sum_{j=1}^{N}\hat{a}_{j}$
per particle is defined as

$\displaystyle\frac{1}{N}\Delta^{2}_{\hat{A}}=\frac{1}{N}\bigg{[}\langle\Psi|\hat{A}^{2}|\Psi\rangle-\langle\Psi|\hat{A}|\Psi\rangle^{2}\bigg{]}.$
(38)

The detail discussion is included in the main text. Here, we report the
convergence of many-particle variances of the position, momentum and angular
momentum operators [76, 86] for all the three confining potentials in the
rotating frame with respect to the number of orbitals.

Figure 16:  Convergence of the many-particle position variances (a)-(b) and
momentum variances (c)-(d) computed for the elongated trap with M=4,8 self-
consistent orbitals along the $x$\- and $y$\- directions.

###  IX.1 The position and momentum variances

This section presents the convergence of many-particle position variances
$\frac{1}{N}\Delta^{2}_{\hat{X}}$, $\frac{1}{N}\Delta^{2}_{\hat{Y}}$ and many-
particle momentum variances $\frac{1}{N}\Delta^{2}_{\hat{P}_{X}}$,
$\frac{1}{N}\Delta^{2}_{\hat{P}_{Y}}$ with respect to the number of orbitals
$M$ for the three confining potentials.

Figures 16(a)-(b) show the convergence of the many-particle position variances
and Figures 16(c)-(d) correspond to the many-particle momentum variances for
the elongated trap [Equation (34)] with $M=4,8$ orbitals.

Figure 17: Convergence of the many-particle position (a) and momentum
variances (c) with M=3,6 self-consistent orbitals for the three-fold symmetric
trap. Panels (b) and (d) correspond to many-particle position and momentum
variances with M=4,8 self-consistent orbitals for the four-fold symmetric
trap.

Figures 17(a),(c) correspond to the many-particle position and momentum
variances for three-fold symmetric trap [Equation (35)]. Similarly, Figures
17(b),(d) depict the many-particle position and momentum variances for the
four-fold symmetric trap [Equation (36)].

It is evident that the convergence with respect to the number orbitals of the
many-particle variances of the position and momentum operators for all the
three confining potentials is excellent.

###  IX.2 Angular momentum properties

Before we discuss the angular momentum variance, we wish to show that the
expectation value of the angular momentum operator discussed in the main text
converges at the many-body level. Then, we proceed to the angular momentum
variance, and augment the text by its analysis.

Figure 18:  Convergence of the expectation values of the angular momentum
operator per particle. (a) for the elongated trap with M=4,8 self-consistent
orbitals, (b) for the three-fold symmetric trap with M=3,6 self-consistent
orbitals and (c) for the four-fold symmetric trap with M=4,8 self-consistent
orbitals.

Figure (18) proves the convergence of the expectation values of angular
momentum operator $\hat{L}_{Z}=\sum_{j=1}^{N}\hat{l}_{z_{j}}$ per particle
$\frac{1}{N}\langle\Psi|\hat{L}_{Z}|\Psi\rangle$, with respect to the number
of orbitals $M$ for the three confining

Figure 19: Convergence of the angular momentum variance. (a) for the elongated
trap with M=4,8 self-consistent orbitals, (b) for the three-fold symmetric
trap with M=3,6 self-consistent orbitals and (c) for the four-fold symmetric
trap with M=4,8 self-consistent orbitals. The insets of (a), (b), and (c)
demonstrate the angular momentum variances at the mean-field and many-body
levels that are essentially the same in our rotating traps.

potentials; the elongated trap [Equation (34)], the three-fold symmetric trap
[Equation (35)] and the four-fold symmetric trap [Equation (36)].

Finally, we discuss the convergence of the many-particle angular momentum
variance $\frac{1}{N}\Delta^{2}_{\hat{L}_{Z}}$ in all the three confining
traps. It is evident from Figure (19) that the angular momentum variances are
fully converged with respect to the number of orbitals for all the three
confining potentials at all $\omega_{r}$. The insets of Figure (19)
demonstrate that in the present rotating traps and interaction strengths the
mean-field and many-body angular momentum variances essentially coincide.