{'fraction_non_alphanumeric': 0.040965340666635056, 'fraction_numerical': 0.01676070361765682, 'mean_word_length': 4.238172109772756, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Novel phases in rotating Bose condensed gas: vortices and quantum correlation

We present the exact diagonalization study of rotating Bose-condensed gas interacting via finite-range Gaussian potential confined in a quasi-2D harmonic trap. The system Hamiltonian matrix is diagonalized in given subspaces of quantized total angular momentum to obtain the lowest-energy eigenstate employing the beyond lowest-Landau-level approximation. In the co-rotating frame, the quantum mechanical stability of angular momentum states is discussed for the existence of phase transition between the stable states. Thereby analyzing the von Neumann entanglement entropy and degree of condensation provide the information about quantum phase correlation in the many-body states. Calculating the conditional probability distribution, we further probe the internal structure of quantum mechanically stable and unstable states. Much emphasis is put on finding the spatial correlation of bosonic atoms in the rotating system for the formation and entry of singly quantized vortices, and then shaping into canonical polygons with and without a central vortex at the trap center. The results are summarized in the form of a movie depicting the vortex patterns having discrete p-fold rotational symmetry with p=2,3,4,5,6.

Here we examine the ground state properties of harmonically trapped quasi-2D Bose-Einstein condensate hosting upto seven vortices with repulsive finite-range Gaussian interaction. Our aim in the present work, is to examine as to how a quasi-two-dimensional Bosecondensed gas with finite-range particle-particle Gaussian interaction responds to an externally impressed rotation in the angular momentum regime 0 ≤ L z ≤ 5N , where N is the number of bosons. The sign of the Gaussian interaction potential is taken to be positive corresponding to repulsive interaction [8]. In order to analyze the many-body ground state, exact diagonalization of the n × n many-body Hamiltonian matrix has been performed in given subspaces of (quantized) total angular momentum with Hilbert dimensionality n ∼ 10 5 (e.g., for L z = N = 16, n = 384559). In constructing the manybody basis states we go beyond the lowest Landau level approximation so as to include the single-particle basis states with radial quantum number n r ≥ 0 and angular momentum quantum number m of either sign [36,37]. The degree of condensation and von Neumann entanglement entropy of angular momentum states are calculated as a measure of many-body quantum correlation. To analyze the internal structure of stable as well as unstable states, the conditional probability distribution is also obtained. Beyond the mean-field approximation, we explore the novel phases of vortex states in subspaces of quantized total angular momentum and their transition to stable vortex patterns (configurations) is examined. In particular, we examine the critical angular velocity corresponding to stable vortex patterns. In response to the externally impressed rotation, singly quantized vortices are observed and at high angular momentum L z , the system shapes up into canonical polygons with or without a vortex located at the center of the harmonic trap. This paper is organized as follows. In Sec. II, we present the theoretical model of a rotating Bose system and describe the many-body Hamiltonian with repulsive Gaussian interaction. In Sec. III, we present the numerical results for N = 16 bosons and discuss the various physical quantities of interest (which are experimentally accessible) to examine in detail the many-body ground state properties of the system in quasi-2D harmonic trap. Finally in Sec. IV, we briefly summarize our results of the present study. Few important supplemental derivations and computational scheme for exact diagonalization calculation are deferred to the Appendixes.

II. MODEL SYSTEM
We consider a system of N interacting spinless bosonic atoms each of mass M , trapped in an external harmonic oscillator potential V (r) = 1 2 M ω 2 r r 2 + ω 2 z z 2 with radial ω r and axial ω z frequencies, given r = x 2 + y 2 the radial distance from trap center. The harmonic trap is assumed to be highly anisotropic with λ z ≡ ω z /ω r ≫ 1 yielding an effectively quasi-two-dimension (quasi-2D) system with x-y rotational symmetry. The system is rotated about the z-axis with trap angular velocity Ω ≡ Ωê z . We chosehω r as the unit of energy and a r = h/M ω r as the corresponding unit length. Introducing Ω ≡ Ω/ω r (≤ 1) as the dimensionless angular velocity and L z (scaled byh) being the z projection of the total angular momentum operator, the manybody Hamiltonian in the co-rotating frame is given as The first two terms in the Hamiltonian (1) correspond to the kinetic and potential energies respectively. The third term U (r i , r j ) arises from the atom-atom interaction, is modelled by the Gaussian potential with parameter σ (scaled by a r ) being the effective range of the potential [52,[58][59][60][61][62] where g 2 = 4πa s /a r is a measure of the strength of effective two-body interaction with a s being s-wave scattering length taken to be positive (a s > 0) so that the effective interaction is repulsive. In the limit σ → 0, the normalized Gaussian potential in Eq. (2) reduces to the zero-range δ-function potential g 2 δ (r − r ′ ) used in earlier studies [4]. For a given value of σ, the s-wave scattering length in the parameter g 2 = 4πa s /a r is adjusted in such a way that N a s /a r becomes relevant to the experimental value [4,36]. In addition to being physically more realistic, the finite-range Gaussian potential (2) is expandable within a finite subspace of single-particle basis functions [63] and hence computationally more feasible compared to the zero-range δ-function potential [52,[58][59][60].

III. NUMERICAL RESULTS AND DISCUSSION
The results are presented for an interacting system of N = 16 Bose atoms of 87 Rb confined in quasi-2D harmonic trap, with radial frequency ω r = 440π Hz corresponding to the trap length a r = h/M ω r = 0.727µm and the aspect ratio λ z (≡ ω z /ω r ) = √ 8 [64]. Using Davidson iterative algorithm [65], exact diagonalization of the n × n Hamiltonian matrix is carried out beyond lowest Landau level approximation [36,37] separately for each of the subspaces of total angular momentum in the regime 0 ≤ L z ≤ 5N . Details of the diagonalization scheme and the beyond lowest Landau level approximation used here, have been presented in Appendix A. It is to be noted that for a many-body system (under consideration here), the characteristic energy scale for the interaction is determined by the dimensionless parameter, N a s /a r , where a s is the s-wave scattering length. Owing to increasing dimensionality of the Hilbert space with increase in number of atoms N , the computation (exact diagonalization) becomes impractical beyond a few hundred atoms, for instance, n = 384559 for L z = N = 16. We, therefore, vary s-wave scattering length a s in our calculations to achieve a suitable value of (N a s /a r ) ≤ 1, relevant to experimental situation [4]. Accordingly, the parameters of two-body Gaussian potential in Eq. (2) have been chosen as: interaction range σ = 0.1 (in units of a r ) and s-wave scattering length a s = 1000a 0 , where a 0 = 0.0529 nm is the Bohr radius. The corresponding value of the interaction parameter g 2 (= 4πa s /a r ) turns out to be 0.9151 leading to (N a s /a r ) ∼ 1 in the moderately interacting regime [36,37]. We begin by examining the response of Bose-condensed system to the external rotation.

A. Stability of quantum states
The system of N bosons confined in a harmonic trap is subjected to an externally impressed rotation along z-axis with angular velocity Ω (≡ Ω/ω r ). In the corotating frame, the simultaneous eigenstate of Hamiltonian and total angular momentum minimizes the energy Ψ H lab − ΩL z Ψ at zero temperature, to become the ground state of the system. This is equivalent to minimizing H lab with respect to Ψ, subject to the constraint that the system has angular momentum expectation value L z with angular velocity Ω identified as the corresponding Lagrange multiplier. Here Ψ is the variationally obtained N -body ground state associated with L z (please see Appendix A). Thereby for the stable states in the rotating frame, one may obtain a number of successive critical angular velocities Ω c L i z ; i = 1, 2, 3, . . . for the angular momentum states at which a particular excited state becomes the ground state of the system. However for unstable angular momentum states such a critical angular velocity doesn't exist, which therefore remain as the excited states. Thus, the critical angular velocity Ω c L i z is the one, beyond which the higher total angular momentum state L i z becomes lower in energy in the rotating frame compared to the lower angular momentum state L (i−1) z < L i z and is given by [36,52]: where E lab (L z ) is the variationally obtained N -body state energy for the total angular momentum state L z in the non-rotating frame.   1 shows the total angular momentum L z versus the angular velocity Ω, for the many-body ground state, and is called the L z − Ω stability curve. There exist a series of critical angular velocities Ω c L i z < 1; i = 1, 2, 3, . . . , at which the ground state of the rotating Bose-condensed system changes its quantized total angular momentum abruptly. For N = 16 bosons with Gaussian repulsion (2), we observe a step-like structure in the L z − Ω stability curve ( Fig. 1) with plateaus at L i z = 0, 16,26,28,32,35,36,48,52,60,71,76. These L i z ground states corresponding to critical angular velocities Ω c L i z , are referred to as quantum mechanically stable states [28][29][30]. Thus, for rotating bosons interacting via finite-range Gaussian potential, one may recovers the step-like structure in the stability curve, well known for the Bose-condensed gas interacting via zero-range (δfunction) potential [49].
Initially, when the Bose-condensed system is nonrotating, the angular momentum state L z = 0 corresponds to the ground state of the system. As the system is subjected to an external rotation, higher angular momentum states (L z > 0), which minimize the free energy, become the ground state of the system. As observed in Fig.1, the state with total angular momentum L z = 0 continues to be the ground state of the system till angular velocity Ω c L 1 z = N , where the state with total angular momentum L z = N , becomes the new ground state in the co-rotating frame. The system thereby undergoes a quantum phase transition from a non-rotating state to the single vortex state, seen as vertical line (a "jump" in ground state total angular momentum) in the stability curve. At angular velocities beyond Ω c L 1 z = N , the ground state angular momentum has a plateau L z = N till angular velocity reaches Ω c L 2 z , where the second jump in ground state angular momentum takes place. Thereafter, one observes a sequence of plateaus and jumps. The rotational angular velocity Ω (≡ Ω/ω r ) grows from lower limit zero to the upper limit Ω = 1 beyond which the condensate will become centrifugally unstable. We notice that with increasing Ω < 1, the width of the plateaus and the hight of the jumps in the stability curve, in general, decreases which may plausibly be attributed to condensate depletion, as discussed later in the section.
The critical angular velocities Ω c L i z are calculated for the onset of different vortex states and undergoes quantum phase transition. At an angular velocity Ω < Ω c L i z , the nucleation of the first vortex state with angular momentum L i z begins which stabilizes at Ω = Ω c L i z and becomes the new ground state of the system. For angular velocities beyond Ω c L i z , the vortex state L i z continues to be the ground state of the system till the nucleation of the next vortex state begins. The angular momentum states L i z = 26, 28, 32, 36, 48, 52, correspond to off-centred vortex states (without a central vortex), as will be seen when we examine their internal structure in later sections. The single-particle angular momentum ℓ z = L z /N is not an integer multiple ofh. In contrast, for axisymmetric multi-vortex states, the total angular momentum L z is quantized in integer units and hence the topological discreteness of the vortices. We next correlate these observations with our results on the many-body ground states for their quantum correlation in terms of von Nuemann entropy as well as degree of condensation.

B. Many-body quantum correlation
An interesting and useful measure of the quantum phase correlation in the many-body ground state of a confined system is provided by the von Neumann entanglement entropy [53][54][55][56][57][58], an extension to the classical Gibbs entropy, defined as whereρ is the single-particle reduced density operator obtained from the many-body ground state wavefunction.
In terms of eigenvalues {λ µ } and eigenfunctions {χ µ (r)}, the single-particle reduced density operatorρ takes the form as in Eq. (B2), for details please see Appendix B.
The von Neumann entropy is, then, evaluated explicitly as in different subspaces of quantized total angular momentum L z , as shown in Fig. 2 for N = 16 bosons. It is clear from Fig. 2 that for a (harmonically confined) rotating Bose gas, the von Neumann entropy S 1 has a tendency to grow with angular momentum L z , but the whole curve shows a slightly oscillatory behavior. Further, the curve of entropy S 1 versus L z exhibits a series of local minima L S1 z . For instance, the local minima appear at angular momenta L S1 z = 0, 16,24,26,28,30,36,39,44,48,52,55,60,66,71,76 as shown in Fig. 2. It is interesting to note that the entropy S 1 plotted in Fig. 2, present features similar to that of stability curve in Fig. 1. In fact, most of these local minima (obtained from S 1 ) correspond to the stable angular momentum L i z -states of ground state Ψ, characterized by critical angular velocities Ω c L i z ; i = 1, 2, 3 . . . (see Fig. 1). Therefore, the von Neumann entropy S 1 plotted in Fig. 2, can also be used to study the L z -Ω stability curve Fig. 1 [55], for zero-range (δ-function) interaction potential in LLL approximation.
The von Neumann entanglement entropy S 1 defined in Eq. (5) also provides information about the degree of condensation C d defined in Eq. (C1), which is sensitive to the loss of macro-occupation, see Appendix C. To ascertain the relation between C d and S 1 , we compare corresponding plots of both these quantities with respect to angular momentum L z , as shown in Fig. 3 and 2, respectively. It is observed from Figures that the value of C d approaches to unity and the value of entropy S 1 reduces to zero, for a perfect Bose-Einstein condensate [66], since all bosons occupy one and the same mode i.e. λ 1 = 1 and λ µ = 0 for µ = 2, 3, . . . . The angular momentum state L z = 0 corresponds to a fairly condensed state with S 1 ∼ 0. As the condensate depletes, with more than one eigenvalue {λ µ } becoming non-zero, the entropy S 1 increases. We observe that lesser the value of S 1 , greater will be the quantum many-body correlation between the particles thus enhancing the condensation. In the single vortex state corresponding to L z = N = 16, the degree of condensation C d has a larger value and there is a dip in the entropy S 1 . For other vortex states too, we have found the dip in S 1 and corresponding peak in C d at the stable angular momentum states associated with the critical angular velocities of the system.
For repulsive Bose-condensed gas rotating with angular velocity Ω c L i z smaller than the trapping frequency ω r , there exists stable states containing the various vortex configurations. The system corresponding to these stable vortex states is well described by a stationary state with some finite non-zero vorticity m 1 , identified as the angu- I. The many-body ground state energies E lab (L i z ) in the laboratory frame, the values of critical angular velocity Ωc L i z ; i = 1, 2, 3, . . . and the p-fold rotational symmetry of each stable vortex state for N = 16 bosons in given subspaces of 0 ≤ Lz ≤ 5N . The single-particle quantum number m1 corresponding to largest eigenvalue λ1 of the SPRDM, the degree of condensation (C d ) and von Neumann entropy (S1) is also reported. Here, the interaction strength g 2 = 0.9151 and range σ = 0.1 of the Gaussian potential (2). lar momentum quantum number of the most dominant single-particle state χ 1 (r) corresponding to the largest condensate fraction λ 1 (macroscopic eigenvalue) of the SPRDM Eq.(B3). For each stable vortex state in the regime 0 ≤ L z ≤ 5N with N = 16 bosons, Table I lists values of the following: (i) the lowest energy E lab (L i z ) of the N -body states in the laboratory frame (ii) critical angular velocity Ω c L i z ; i = 1, 2, 3, . . . (iii) respective discrete rotational symmetry of the vortex states denoted by integer p, (iv) along with the vorticity defined by single-particle angular momentum quantum number m 1 of the SPRDM, (v) the degree of condensation C d and (vi) the von Neumann entanglement entropy S 1 of harmonically confined rotating Bose-condensate, corresponding to the stable states as observed in L z -Ω stability curve (see Fig. 1). In order to gain an insight we examine the spatial correlation by looking into the internal structure of the condensate for various stable as well as unstable vortex states.
C. Internal structure of vortex states Circular symmetries of the rotating Hamiltonian often hide the internal structure in the exact many-body state. The pair correlation function, a standard tool in many-body physics, is often used to analyze the internal structure of exact many-body state. Here, the pair correlation function is conditional probability distribution (CPD) which provides a test for the presence of spatial correlation among bosons. The CPD P (r, r 0 ) is defined as the probability of finding a particle at r, given the presence of another at r 0 and is computed by the expression [37-40, 67, 68]: Here |Ψ = ν C ν |Φ ν is the many-body ground state obtained through exact diagonalization and r 0 = (x 0 , y 0 ) is the reference point. The details of CPD in finite system are sensitive to the selection of this reference point r 0 [69]. Usually, the reference point r 0 is given a value of the order of spatial extent of the confined system, so that the effects of pair repulsion between atoms do not affect the vortex structure significantly. For large bosonic systems, CPD show the vortices clearly by giving r 0 a significantly large value. However in case of systems with small number of particles N (or depleted condensate); in order to obtain optimum information, the reference point r 0 is taken as a position where the density of the system is maximum. The reference point r 0 in our calculation is taken to be located on the x-axis i.e. r 0 = (r 0 , 0) for a system of rotating Bose-condensed gas of N = 16 particles. For slow rotating regime (0 ≤ L z ≤ 2N ) where the condensate is large, we choose r 0 = (3, 0). However for moderately to rapidly rotating regime (2N < L z ≤ 5N ) where the condensate depletes due to centrifugal effect, we choose r 0 = (1.5, 0). We have carefully selected the position of reference point r 0 to ensure that it never coincides with the vortex position (low density region). It is observed that a small displacement of the reference point r 0 alters the pair correlation considerably for unstable states. The reference point thus acts as an infinitesimal symmetry breaking perturbation on the two-body correlation function (spontaneous symmetry breaking). We have reported here the calculation for quantum mechanically stable ground states in the rotating frame (corresponding to critical angular velocities Ω c L i z , i = 1, 2, 3, . . . ) for harmonically confined rotating Bose gas in the regime of 0 ≤ L z ≤ 5N . It appears that N = 16, is sufficiently large system to give insight into vortex nucleation. Calculating the CPD using Eq.(6) for quantized angular momentum L z -values, one may study the internal structure (spatial correlation) of N -body quantum states. CPD is indeed an observable quantity in experiments which indicates the tendency of a system's evolution. Unlike the usual density distribution, which is circularly/cylindrically symmetric under rotational invariant confinement, the CPD plots exhibit asymmetric behavior and reflect the intrinsic density distribution of bosons for each successive critical angular velocity Ω c L i z . The respective CPD plots for N = 16 are presented in Fig. 4 which depicts the condensate density isosurface, viewed along the axis of rotation (z-axis). The value of quantized total angular momentum L z is marked in the upper left corners of each CPD plot. Knowing the stability of vortex states from L z -Ω stability curve (see Fig. 1), we can extract valuable information about the vortex patterns (configurations) and their formation. The angular momentum values associated with the Sta- A striking feature of the rotating Bose-condensate is their lack of full rotational symmetry. Instead they have only p-fold symmetry with p = 2, 3, 4, 5 and 6 as shown in CPD plots of Fig. 4. Thus CPD can also be regarded as a measure of the possible spontaneous symmetry breaking. A succession of discontinuous, symmetry-changing transitions follow at higher values of angular momenta L z , with each new stable state corresponding to a differ-ent configuration of vortices. To elucidate, we have also produced a movie featuring the CPD contours of manybody states which provide information about evolution of the system with different values of L z such as the formation of various pattern of vortices along with its discrete p-fold rotational symmetry (see the Supplemental Material [70]). The gradual entry of vortex in the rotating system can clearly be seen, first changing to a central single-vortex state at L z = N = 16, and then to higher vortex states. When the angular momentum is increased beyond L z > N , additional vortices appear in the form of regular structures as shown in Fig. 4. As their number grows they may organise themselves in a triangular Abrikosov lattice [16,17,71].
At L z = 26, 28 and 32, the second vortex is appeared in the confined system as single-particle quantum number m 1 = 1 (vorticity associated with first centred vortex state) changes to m 1 = 2 with the formation of 2-fold vortex pattern. Out of these two-fold symmetry vortex states observed in L z − Ω stability curve shown in Fig. (1), the L z = 28 state with larger plateau size found to be more stable phase coherent state corroborated by the lower value of S 1 and higher value of C d . We further observe that in the angular momentum regime 22 ≤ L z ≤ 63 for N = 16, no vortex is formed along the axis of the trap as seen in Fig. 4 (also refer to the [70]). It may be noted that the stable angular momentum state L z = 60 corresponds to 5-fold rotational symmetry of pentagonal vortex pattern without a central vortex at the trap center. On the other hand the L z = 71 and 76 represent the 5-fold and 6-fold rotational symmetry respectively, of pentagonal and hexagonal vortex patterns with a central vortex. Unlike the six-fold rotational symmetry of hexagonal vortex pattern identified as a precursor to thermodynamically stable triangular vortex lattice, the vortex pattern with five-fold rotational symmetry found to exist in confined system will not survive in the thermodynamic limit. Looking beyond the mean-field approximation it is important to mention here that, though we do not observe the linear configuration of three vortices as suggested in [72], our exact diagonalization result bears the signature of triangular vortex lattice [71]. We observe that with increasing angular momentum L z , the system expands radially due to centrifugal force and consequently it becomes less dense. The centrifugal effect is prominent for the multi-vortex states but becomes fairly small near the trap center. There is an upper limit on the rotational angular velocity in a harmonic trap. This limit is the trap frequency ω r above which the center-ofmass of the condensate destabilizes [73]. In the vicinity of high angular velocity Ω < ∼ ω r , the centrifugal force influences the shape of the condensate. Notice that the non-rotating L z = 0 state is more compact relative to the expanded higher angular momentum states.
In Fig. 5, we have plotted three largest eigenvalues λ 1 > λ 2 > λ 3 of the SPRDM obtained from the ground state of the system versus the quantized total angular momentum in the regime 0 ≤ L z ≤ 5N . It can be observed from the figure that the nucleation of the first centered vortex (at L z = N ) in the rotating Bosecondensate, does not occur through a smooth entrance of the vortex. Indeed, the system passes through a correlated meta-stable angular momentum state at L z = 12, where the system achieves critically re-arrangement in configuration space, as evident by the von Neumann entropy (S 1 ) being maximum (see Fig. 2) and the two largest condensate fractions λ 1 (= 0.4543), λ 2 (= 0.4059) becoming comparable [74]. Thus the system is preparing to undergo a critical phase transition from L z = 11 to L z = 12 state and the corresponding transition of singleparticle quantum number m 1 = 0 to m 1 = 1 takes place, associated with largest condensate fraction λ 1 . A similar pattern is observed for the entry of the second vortex at angular momentum state L z = 21 with λ 1 ∼ λ 2 . The system is undergo a critical change in spatial correlation from L z = 21 to L z = 22 state and the corresponding transition of single-particle quantum number m 1 related to λ 1 , that is m 1 = 1 to m 1 = 2 takes place. It is also observed that for a particular vortex state, the vorticity m 1 corresponding to largest condensate fraction λ 1 of SPRDM (B2) always attains the same value, irrespective of different values of quantum number n. For example, one may refer to Table II where the single quantized vortex state has the same value of m = 1 for quantum numbers n = 1 and 3. We found that the total angular momentum states with L z = 12, 22, 33, 43 etc., where the quantum jump of single-particle angular momentum m 1 takes place, correspond to the unstable states that provide a physical mechanism for the entry of vorticity into the condensate (refer to Table III). For better understanding, we focus on the angular momentum regime 0 ≤ L z ≤ N , in which the first vortex nucleates. The respective CPD plots calculated using Eq.(6), are shown in Fig. 6, which depicts the isosurface density profile viewed along the axis of rotation (z-axis). The non-rotating state with angular momentum L z = 0 has Gaussian density profile along the axis of the trap as seen in Fig. 6(a). The density profile of this state has complete x-y rotational symmetry of the harmonic trap. The maximum number of bosons are confined around the trap center and occupy the non-rotating lowest singleparticle state with angular momentum quantum number m 1 = 0 corresponding to largest condensate fraction λ 1 (the macroscopic eigenvalue of the SPRDM). It may also be noted from Table I that in the absence of external rotation the vortex-less state for L z = 0 is the most stable phase coherent state with the maximum value of degree of condensation C d = 0.9857 and minimum value of von Neumann entropy S 1 = 0.0890. For an axially symmetric harmonic potential confining N (≫ 1) bosons, the CPD plots may be directly related to the position of the vortex with respect to the trap center. The regime 0 < L z < N , in above plots of Fig. 6, describes an off-center unstable vortex state, for instance see Figs. 6(b) to 6(h). The total angular momentum state L z = N corresponds to a singly quantized vortex for which the axis of rotation coincides with the trap center as seen in the Fig. 6(i). This state is identified by a high condensate fraction λ 1 (see Fig. 5) with the single-particle quantum number m 1 = 1 as shown in Table I, for which there is a peak in C d curve (Fig. 3) and a corresponding dip in S 1 (Fig. 2). The orientation of vortex can be represented as a measure of the possible spontaneous breaking of rotational symmetry. In our system of harmonically confined N = 16 rotating bosons, we observe that the circular symmetry is broken spontaneously, beyond L z = 11 as shown in Fig. 6(d).
The entry of first vortex is clearly seen to be from the periphery of the trap (L z = 0). It is gradually moving inwards to its center (L z = N ) with increasing total angular momentum L z .

IV. SUMMARY AND FUTURE WORK
In summary, we have studied the novel phases of a rotating system of N = 16 bosonic atoms confined in a quasi-2D harmonic trap and interacting via repulsive finite-range Gaussian potential. The many-body ground states were obtained using exact diagonalization method within subspaces of total angular momenta in the regime with 0 ≤ L z ≤ 5N , which goes beyond the mean-field approximation. By minimizing the energy of system in the co-rotating frame at zero-temperature, first, we have discussed the rotational critical angular velocity associated with the quantized total angular momentum introduced as L z − Ω stability curve. From this we were able to provide information about the mechanical stability of angular momentum states, most notably the quantum phase transition between the stable vortex states. In order to analyze the criterion for the existence of Bose-condensate in the rotating system considered, we have calculated the first order correlation as embodied in single-particle reduced density matrix. Many-body quantum correlation measured in terms of degree condensation and von Neumann entanglement entropy as a function of L z , provided the direct evidence of quantum mechanically phase coherent states in consistent with the stability curve. Due to robust nature of quantum phase coherence in the rotating system, the stable vortex states corresponding to the successive critical angular velocities has been found with lower value of entanglement entropy and higher value of degree of condensation. Thereby internal structure of stable and unstable sates has been investigated by calculating the conditional probability distribution (CPD) of angular momentum states in the co-rotating frame. It has been observed that the mechanism of formation and entrance of the singly quantized vortices in the considered system is not smooth, indeed it passes through a critically re-arrangement in the configuration space, advocated by the analysis of internal structures in the unstable stable vortex states where the two largest condensate fractions become comparable. The stable vortex states with upto seven singly quantized vortices has been studied possessing the p-fold discrete rotational symmetry with p = 2, 3, 4, 5, 6. The central vortex initially nucleated at L z = N = 16 state is peculiarly reappeared at L z = 71 and 76 states forming a pentagon vortex pattern with vorticity m 1 = 6 and a hexagon vortex patterns with with vorticity m 1 = 7, respectively. In addition to the existence of five-fold rotational symmetry of vortex pattern in the confined system, our exact diagonalization results bear the signature of six-fold rotational symmetry precursor to thermodynamically stable triangular vortex lattice. The complete exact diagonalization results for the model system considered here are summarised in a Table III. To make the physical picture of condensate density isosurface more richer, we have also presented the CPD contour plots in the form of a movie [70] depicting the formation of various vortex configurations from one stable state to the other with rotation. Finally, we emphasize that the ground state properties, namely, critical angular velocity, degree of condensation, von Neumann entropy and CPD, analyzed in the present work are all measurable quantities and hence experimentally relevant. In future work we would like to investigate the effect of finite-range Gaussian interaction on the quantum mechanical stability of the angular momentum states and discuss how the range of interaction can be used as a probe to further explore the internal structure of vortex states.
Appendix A: Computational scheme Since the system being studied here is rotationally invariant in the x-y plane, the z-component of the total angular momentum L z is a good quantum number leading to block diagonalization of the N -body Hamiltonian matrix in subspaces of quantized total angular momentum [36]. To obtain the many-body eigenstates, we carry out exact diagonalization of the Hamiltonian matrix in different subspaces of L z with inclusion of lowest as well as higher Landau levels of the single-particle basis states in constructing the N -body basis states [36,37,40]. The associated Hilbert space may be restricted to the space spanned by the normalized single-particle basis functions u n,m (r, θ) u nz (z) ≡ u n,m,nz (r) spanning the quasi-2D plane, where (n, m, n z ) is set of single-particle quantum numbers.

The single-particle basis states
For the construction of N -body basis states, the singleparticle basis functions are chosen to be the eigenfunctions of the non-interacting single-particle Hamiltonian identified as the harmonic oscillator Hamiltonian with an externally impressed rotation about the z-axis given the single-particle angular momentum ℓ z . The eigensolutions for H Ω sp u n,m,nz (r) = ǫ n,m,nz u n,m,nz (r), in dimensionless form, are known to be: ǫ n,m,nz = (n + 1 − mΩ) + λ z (n z + 1/2) u n,m,nz (r) = 1 2 {n − |m|} ! λ z /π 3 1 2 {n + |m|} ! 2 nz n z ! e −(r 2 +λz z 2 )/2 × e imθ r |m| L |m| 1 2 (n−|m|) r 2 H nz λ z z 2 (A2) where n = 2n r + |m| with the radial quantum number n r and the single-particle angular momentum quantum number m. Here L |m| 1 2 (n−|m|) r 2 is the associated Laguerre polynomial and H nz λ z z 2 the Hermite polynomial. With λ z ≫ 1, the system is taken to be quasi-2D and hence there is practically no excitation along the relatively stiffer z-axis. We, therefore, set n z = 0 in Eq. (A2) implying that all the particles occupy the lowest-energy state u 0 (z) = (λ z /π) 1/4 e −λz z 2 /2 of z co-ordinate degree of freedom, thereby reducing Eq. (A2) to ǫ n,m = (n + 1 − mΩ) + λ z /2 , with n = 2n r + |m| Restricting to n r = 0 and taking m ≥ 0 corresponds to the lowest Landau level (LLL) approximation. For harmonically confined interacting Bose gas, most of the mean field as well as the many-body calculations have used only the LLLs. There has been some studies [36,37] wherein pointed out that for slow rotating [75] and moderately to strongly interacting [76] bosons, it becomes necessary to consider the single-particle states beyond LLL approximation with n r ≥ 0 and m taking positive as well as negative values in constructing the many-body basis states.
The single-particle basis functions thus chosen are used to construct the N -particle basis states {Φ ν } in the manybody variational trial wavefunction Ψ (r 1 , r 2 , . . . , r N ) of the system for given value of total angular momentum L z . Thus, the many-body variational wavefunction for bosons, is a linear combination of the N -body basis states {Φ ν } given by the expansion: where {C ν } are the variational parameters to be determined.

Construction of N-body basis states
The N -body basis states {Φ ν (r 1 , r 2 , . . . , r N )}, for the given subspace of L z , are constructed as the symmetrized products a finite number of normalized single-particle basis functions u j (r) ≡ u n,m (r) defined in Eq. (A3): where P permutes the co-ordinates of N -particles and ν j is the occupancy of jth single-particle basis function u j (r). The many-body index ν ≡ (ν 0 , ν 1 , . . . , ν j , . . . , ν k ), labeling the N -body basis function Φ ν (r 1 , r 2 , . . . , r N ), stands for a set of single-particle quantum numbers {j ≡ (n, m)} and their respective occupancies {ν j }. In the second-quantized notation, the Bose field operator is expanded in terms of single-particle basis functions asψ (r) = jb j u j (r) whereb j (b † j ) is the annihilation(creation) operator for the Bose quanta in state j, obeying commutation rules. In occupation-number representation, the N -body basis function in Eq. (A5) for a given subspace of L z is written in second-quantized form as: with k j=0 ν j = N and k j=0 m j ν j = L z where j = (n j , m j ). With these constraints on {ν j , m j }, only the most important Fock states (spanning the active Fock space) constructed from the full single-particle basis set, with a given L z are included. This procedure reduces the dimension of the Hamiltonian matrix to the order of 10 2 − 10 5 , amenable to diagonalization on a small workstation of 3.20 GHz processor.

Diagonalization of the Hamiltonian
The many-body Hamiltonian (1) in terms of bosonic field operators, is written asĤ lab = ψ † (r) H spψ (r) dr + 1 2 ψ † (r)ψ † (r ′ ) Uψ (r ′ )ψ (r) dr dr ′ . Expanding the Bose field operators in terms of single-particle basis statesψ (r) = jb j u j (r), the Hamiltonian is obtained asĤ are the usual bosonic creation (annihilation) operators, i|H sp |j and i, j|U |k, l are the one-body and the two-body matrix-elements respectively over the single-particle basis functions chosen. Once the active Fock states are constructed as in (A6), we calculate the matrix elements and subsequently diagonalize the Hamiltonian matrix using Davidson algorithm [65]. The variational parameters {C ν } for the ground state wavefunction Ψ = ν C ν Φ ν are determined by minimizing E {Ψ} ≡ Ψ|Ĥ lab |Ψ − E lab Ψ|Ψ with respect to Ψ. The Lagrange multiplier E lab is identified as the variational energy for the ground state.
The above expression is written in terms of single-particle basis functions u n (r) with quantum number n ≡ (n, m).
Being hermitian, this can be diagonalized to give where χ µ (r) ≡ n c µ n u n (r) and µ λ µ = 1 with 1 ≥ λ 1 ≥ λ 2 ≥ · · · λ µ ≥ · · · ≥ 0. The {λ µ } are the eigenvalues, ordered as above, and {χ µ (r)} are the corresponding eigenvectors of the SPRDM (B2); each µ defines a fraction of the BEC. According to the Penrose-Onsager criterion, "B.E. condensation is said to be present whenever the largest eigenvalue of the single-particle reduced density matrix is an extensive rather than an intensive quantity" [77]. Conveniently, a fragmented state is defined as a state for which there exits more than one macroscopically large eigenvalue of the SPRDM [78]. The fragmented ground state can as well be viewed as a quantum mechanical average over the symmetry broken single condensed states [79].   For the quasi-2D system, where z degree of freedom in the SPRDM has been traced out, Eq. (B2) reduces to where index µ = 1, 2, 3, . . . are the new labels for the single-particle basis quantum numbers corresponding to each fragment of the condensate. We further have found from our calculation that the eigenvector corresponding to a given fragment has the form For a particular total angular momentum L z -state, every fraction λ µ of the SPRDM is characterized by a value of single-particle angular momentum quantum number m µ . It is observed that m µ corresponding to a condensate fraction always takes the same value, irrespective of quantum number n. For example, one may refer to Table II where the single quantized vortex state have the same value of m µ = 1 for quantum numbers n = 1 and 3 corresponding to µ = 1. For a many-body system, the realization of BEC corresponds to a single eigenvalue λ 1 becoming significantly larger than the remaining eigenvalues, i.e. λ 1 ≫ λ (µ≥2) ∀ µ, and identified as the largest condensate fraction. The corresponding eigenvector χ 1 (r) = n c 1 n u n (r), takes the role of the macroscopic order parameter of the system, in the mean-field Gross-Pitaevskii description. The quantum number m 1 , corresponding to the largest fractional amplitude c 1 n 2 of single-particle function u n (r), is recognized as vorticity of the condensate.
Appendix C: Degree of condensation It is to be noted that the usual definition of condensation for a macroscopic system, given by the largest eigenvalue λ 1 of the SPRDM, is not appropriate for systems with small number of particles being studied here [45,58,61]. For example, in the absence of condensation, there is no macroscopic occupation of a single quantum state and all levels are equally occupied, such a definition would imply a condensate though with small magnitude. To circumvent this situation, one introduces a quantity which is sensitive to the loss of macroscopic occupation called the degree of condensation defined as µ=2 λ µ is the mean of the rest of eigenvalues. It can be seen that the degree of condensation, defined as in Eq. (C1), approaches zero for equal eigenvalues, as one would expect. In Fig. 3, we present the variation of degree of condensation C d with total angular momentum L z for N = 16 bosons. TABLE III. Complete exact diagonalization results: For N = 16 bosons in given subspaces of total angular momentum 0 ≤ Lz ≤ 5N , the lowest eigenenergy E lab (Lz) in units ofhωr of the states in the laboratory frame, the values of critical angular velocity Ωc L i z with p-fold rotational symmetry of stable vortex states, the largest three eigenvalues λ1 > λ2 > λ3 and the corresponding single-particle quantum numbers (n1, m1) , (n2, m2) and (n3, m3) in the SPRDM (B2) for the ground state of the rotating BEC. Also given are the values of degree of condensation C d (C1) and von Neumann entropy S1 (5) of the many-body quantum states. Here, interaction parameter g 2 = 0.9151 and range σ = 0.1 of the repulsive Gaussian interaction potential (2). See supplementary material [70] for the internal structure of each tabulated angular momentum Lz-state (includes the stable and unstable states) which clearly shows the nucleation and entry of the vortices alongwith its arrangement to the various stable vortex configurations in a harmonically confined system.