Quantum correlations and degeneracy of identical bosons in a 2D harmonic trap

We consider a few number of identical bosons trapped in a 2D isotropic harmonic potential and also the $N$-boson system when it is feasible. The atom-atom interaction is modelled by means of a finite-range Gaussian interaction. The spectral properties of the system are scrutinized, in particular, we derive analytic expressions for the degeneracies and their breaking for the lower-energy states at small but finite interactions. We demonstrate that the degeneracy of the low-energy states is independent of the number of particles in the noninteracting limit and also for sufficiently weak interactions. In the strongly interacting regime, we show how the many-body wave function develops holes whenever two particles are at the same position in space to avoid the interaction, a mechanism reminiscent of the Tonks-Girardeau gas in 1D. The evolution of the system as the interaction is increased is studied by means of the density profiles, pair correlations and fragmentation of the ground state for $N=2$, $3$, and $4$ bosons.


I. INTRODUCTION
The problem of a particle trapped in a harmonic trap is one of the best known quantum systems. Going from a single particle to a system composed of N interacting particles is, however, far more involved. Interestingly, recent advances in ultracold atomic gases have opened the possibility of studying systems of a few atoms, either fermions or bosons, trapped in potentials of different kind [1][2][3][4].
For the bosonic case, there are important results in 1D where the fermionization of the bosonic gas was preconized by Tonks and Girardeau [5] for the case of infinitely repulsive bosons and later confirmed experimentally in ultracold atomic gases [6,7]. There are many works studying fermionization in 1D, for instance, in optical lattices [8], in few-atom mixtures [9][10][11], for attractive interactions [12] and for few dipolar bosons [13]. In other cases, the focus are quantum correlations [14,15], its effects in mixtures of distinguishable and identical particles [16] and analytic ansatz to capture the physics in all interaction regimes [17].
The case of two particles with contact interactions was considered in one, two and three dimensions in Ref. [18]. There, they obtained semi-analytic results finding the energies and wave functions as solution of transcendental equations. More general cases of few-body systems have been studied mostly in 3D, see Ref. [19] and references therein.
In 2D, semi-analytical approximate solutions to the case of two bosons with finite range interactions have been presented in Ref. [20]. Other 2D works include two and three-body exact solutions for fermions and bosons with contact interaction [21], fast-converging numerical methods for computing the energy spectrum for a few bosons [22], the study of finite-range effects [23,24] and universality [25,26], condensation in trapped fewboson systems [27], and interacting few-fermions systems [28,29].
In this paper, we study the properties for N = 2, 3, and 4 identical bosons interacting through a finite-range interaction confined in a 2D isotropic harmonic trap by means of direct diagonalization of the Hamiltonian. We analyze the properties of the system as we increase the strength of the interaction, going from the noninteracting regime to the strongly interacting one. In Sect. II, we present the many-body Hamiltonian, including the two-body Gaussian-shaped interaction potential considered. We split the center-of-mass and relative parts of the Hamiltonian making use of Jacobi coordinates. In Sect. III, we consider first the noninteracting Hamiltonian and discuss in some details the degeneracies present in the many-body spectrum. In Sect. IV, we focus on the effect of interactions on the many-body spectrum of the system. In Sect. V, we discuss the correlations which build in the ground state as the interaction is increased. Finally, the conclusions and summary are presented in Sect. VI.

II. THE N -BOSON HAMILTONIAN
We consider a system of N identical bosons of mass m trapped by an isotropic harmonic potential. The manybody Hamiltonian in first quantized form reads In usual ultracold atomic gases experiments, the atomatom interactions are well approximated by a contact potential. In our case, we use a finite-size Gaussian potential, where g and s characterize the strength and range of the interaction, respectively. Both parameters are considered to be tuneable. For instance, g can be varied by means of a suitable Feshbach resonance. In the limit of s going to zero, we recover a contact interaction with strength g. Regardless of N , we can split the Hamiltonian in two parts, H = H cm + H r , using Jacobi coordinates, x i , k = 1, ... , N − 1 .
The center-of-mass part and relative part of the total Hamiltonian read with the definitions M ≡ N m and µ ≡ m/2. The interaction only appears in the relative part and takes the formṼ ( r 1 , ... , r N −1 ) ≡ N i<j V x i ( R, r k , ... , r N −1 ) − x j ( R, r k , ... , r N −1 ) .
As a consequence, the change in the energy spectrum with increasing the interaction through g or changing the range s will come from a change in the energy associated to H r .

A. Second-quantized N -boson Hamiltonian
Our numerical method to study the excitation spectrum will consist in truncating the Hilbert space of the N -boson system such that the particles can populate only the first M single-particle eigenstates. We label the single particle states, ψ i (x, y), and their corresponding eigenenergies, ǫ i = n x + n y + 1, with an index i = 1, ... , M running through the pair of quantum numbers n x and n y . With this truncation, the second quantized Hamiltonian readsĤ WhereĤ 0 andĤ int correspond to the single particle and interaction terms, where The explicit analytical form of these integrals, V i,j,k,l , is provided in Appendix C. The operatorâ † i (â i ) creates(destroys) a particle in the single particle mode i, a † i |n 1 , ... , n M = √ n i + 1 |n 1 , ... , n i + 1, ... , n M , a i |n 1 , ... , n M = √ n i |n 1 , ... , n i − 1, ... , n M . (10) They satisfy bosonic commutation relations, [â i ,â † j ] = δ i,j . We introduce the Fock basis, where |vac ≡ |0, ... , 0 is the vacuum state and, as we consider a fixed number of particles N , the quantum numbers n i verify The dimension of the Fock space is which, for N = 2, 3, and 4, gives, respectively,

III. DEGENERACIES IN THE NONINTERACTING LIMIT
In this section, we will discuss the degeneracies present in the system in absence of interactions. First, we consider the two-boson case, in which the analysis is simpler, and then we shall explain the main degeneracies for N bosons.

A. The two-boson system
In the noninteracting case, g = 0, for the two-boson system, we can write down the Hamiltonian in second quantization, splitting the center of mass and the relative motion. Using from now on harmonic oscillator units, ω for energy and /(mω) for length, we have, in polar coordinates,Ĥ =Ĥ cm +Ĥ r =n cm +n r + 2, whereĤ cm =n cm + 1,Ĥ r =n r + 1. Therefore, we have a 2D harmonic oscillator for each part of the Hamiltonian. The corresponding eigenstates can be labelled as |n cm , m cm , n r , m r , namely, n cm |n cm , m cm , n r , m r = n cm |n cm , m cm , n r , m r , n r |n cm , m cm , n r , m r = n r |n cm , m cm , n r , m r , L z,cm |n cm , m cm , n r , m r = m cm |n cm , m cm , n r , m r , L z,r |n cm , m cm , n r , m r = m r |n cm , m cm , n r , m r , whereL z,cm andL z,r are the third component of the center-of-mass orbital angular momentum and the relative orbital angular momentum, respectively, expressed in units of . However, those four quantum numbers have a restriction imposed by the symmetry of the wave function under the exchange of particles. The full wave function in polar coordinates for R and r reads χ ncm,mcm,nr,mr (R, r, ϕ R , ϕ r ) = The L k n (x) are the associated Laguerre polynomials defined as and N n,m (α) is a normalization constant, The wave function corresponding to the center of mass is symmetric under the exchange of particles, because R and ϕ R remain unchanged upon exchanging particles 1 and 2, since R = 1 2 ( x 1 + x 2 ). However, the relative wave function is symmetric or antisymmetric depending on the quantum number m r . We have defined the relative coordinate as r = x 1 − x 2 , therefore the angle ϕ r changes to ϕ r + π and, due to the form of the wave function, see Eq. (17), a factor (−1) mr appears. For this reason, only the states with m r = even can describe the twoboson system. This implies that n r must also be an even With the previous possible quantum numbers, we can determine the degeneracy for each energy level. We define the excitation energy number as the excitation energy per energy unit, N E ≡ E − E 0 . Then, the degeneracy for a given value of N E (see Appendix A) is where ⌊N E /2⌋ indicates the floor function of N E /2. The previous equation is valid for spinless bosons, which is the case considered in this work. However, for fermions and bosons with spin, the spatial antisymmetric states should be considered. The degeneracy for those states (see Appendix A) is Notice that the total degeneracy is given by [30],

Unperturbed energy states
We are also interested in knowing how many states have m r = 0 for each energy level, because these states are the ones that do not feel a zero-range interaction. For a finite but small range, these states are also expected to remain almost unperturbed for the considered range of interaction strengths. The number of states in each energy level such that their energy should not change significantly with a small Gaussian width (see Appendix A) is B. N -boson system The procedure described above in order to compute the degeneracy is not valid for systems with more than two bosons. The reason is that we cannot label the symmetric (neither the antisymmetric) states under the exchange of a pair of particles using the previous quantum numbers. The symmetry of the relative Jacobi coordinates, defined in Eq. (3), under the exchange of two particles is not well defined. An alternative way for counting the degeneracy is by making use of the Fock basis introduced in the previous section, Eq. (11). Those states are eigenstates of H 0 , i.e.,Ĥ The ground state of a system of N identical spinless bosons in a 2D isotropic harmonic potential is always non-degenerate. In particular, for the noninteracting case, it corresponds to a state with all the bosons populating the non-degenerate single-particle ground state, i.e., the state |N, 0, ... , 0 . For any higher energy level of this system, labelled with From left to right, if we have reached d max NE , one of the degenerate states is the one with N E bosons in the single-particle states with excitation energy, E sp exc = E sp − E sp 0 = 1. Therefore, we have N ≥ N E bosons. From right to left, if we have N ≥ N E bosons, we have reached the maximum degeneracy because having less bosons would not allow us to have the previous discussed state, which is degenerate. Adding more bosons would not increase the number of degenerate states, since it is impossible to introduce new states with the same energy as the previous ones. This is due to the finite ways of decomposing N E as a sum of positive integers, without considering the order, that is, the number of partitions p(N E ) [31,32].
Therefore, the degeneracy of the first N E + 1 energy levels is independent of the number of particles N for any N ≥ N E . In Table II, we give the low-energy states with their corresponding energies, excitation energy numbers and degeneracies for a system of N bosons. In Table III, we give d max NE for the first values of N E . Computing the maximum degeneracy is analogous to computing the number of partitions of the integer N E where there are n + 1 different kinds of part n for n = 1, 2, 3, ..., [33] and we can obtain it from its generating function, and also, 1  1  2  N+2  2  2  6  N+3  3  3  14  N+4  4  5  33  N+5  5  7  70  N+6  6  11  149   TABLE III: Energy, excitation energy number, number of partitions of the excitation energy number and maximum degeneracy for the low-energy levels of a system of N noninteracting identical bosons in a 2D isotropic harmonic potential. The maximum degeneracy, d max N E , is equal to the degeneracy of the level NE if and only if N ≥ NE (see the text for explanation).
where P L(k) are the planar partitions of k [34]. Notice that the number of partitions is a lower bound of the maximum degeneracy, and the equality would hold for non-degenerate singleparticle states, e.g. for the 1D case.

IV. ENERGY SPECTRA
Our numerical method consists in the direct diagonalization of the truncated second-quantized Hamiltonian, as described in Sec. II A. We will consider systems with N = 2, 3, and 4 bosons. Direct diagonalization provides the energy spectrum of the Hamiltonian in the truncated space. In particular, we have used the ARPACK implementation of the Lanczos algorithm to obtain the lower part of the many-body spectrum.

A. Two-boson energy spectrum
In Fig. 1, we show the low-energy spectrum for the system of two interacting identical bosons in the harmonic trap. In the figure, we compare results obtained with three different values of s = 0.1, 0.5 and 1. In all cases, the energy spectrum has a number of common features.
First, in the spectrum, there are the states discussed in Sect. III A 1, which are essentially insensitive to the interaction. In the zero range limit, these are basically states with non-zero relative angular momentum, which do not feel the contact interaction [20]. With finite interactions but for a small range, s = 0.1 and 0.5, they remain mostly flat for g up to 20. For s = 1, their energy increases slightly with g, deviating from the zero range prediction.
Second, the ground state energy increases linearly with g for small values of g. Up to first order perturbation theory, the energy is given by . However, the ground state energy tends to saturate as g is increased. For smaller values of s, this saturation takes place at smaller values of g.
Third, there are the energies coming from the relative part of the Hamiltonian with the center of mass at the ground state, i.e. n cm = 0. The ground state is one of these states and there is one state of this type in each energy level with an even N E in the noninteracting limit.
Finally, the spectrum also contains center-of-mass excitations [18], which are easily recognized as constant energy shifts independent of g with respect to states with n cm = 0.
For comparison, we depict also the approximate values of [20] in panel (a) of Fig. 1. As reported in Ref. [20], their approximate solution -which is not variational -starts to deviate from the exact numerical results at values of g ≃ 4. The approximation gives, however, a fairly good overall picture of the low-lying two-particle spectrum.

B. Degeneracy for the interacting two-boson system
We can label the states with three quantum numbers. Two are the ones corresponding to the center of mass, n cm and m cm , and the other is a new quantum number, ν r , that labels the nondegenerate eigenstates of the relative part of the Hamiltonian. We can write those states as where χ ncm,mcm ( √ 2, R, ϕ R ) is given in Sect. III and f νr (r) is the relative wave function, that depends on g and s.
The other states that are in the spectrum are the unperturbed ones (almost unaffected by the interaction). Their degeneracy is given in Sec. III. The states of Eq. (30), for a given ν r , are degenerate with degeneracy given by the 2D harmonic oscillator of the center-of-mass part, i.e., their degeneracy is n cm + 1. From each noninteracting energy level with even N E , a state with a new ν r arises, and its center-of-mass excitations appear in higher energy levels with degeneracy n cm + 1, too. To sum up, the ground state is nondegenerate. The first excited state is two-degenerate and the two states are the two possible center-of-mass excitations of the ground state. The third noninteracting energy manifold (6 states with E(g = 0) = 4) splits in three: 1) three center-of-mass excitations of the ground state, 2) two unperturbed states and, 3) the new relative state with quantum numbers n cm = 0, m cm = 0 and ν r = 1 with E(g = 2) = 4.21.
We give the degeneracy and the quantum numbers of the low-energy states in Table IV. C. Three and four-boson energy spectra Our exact diagonalization scheme allows us to obtain the lowest part of the many-body spectrum for systems of up to 4 bosons with good accuracy, up to values of g ≃ 20. In Fig. 2 we report the ground state energy for ncm nr mcm mr νr E(g = 0) E(g = 2) dint(g = 2) 0 -0 0 1 2 2  TABLE IV: Quantum numbers, energy in the noninteracting limit, energy at g = 2 and degeneracy, for the low-energy levels of a system of two interacting identical bosons trapped in a 2D isotropic harmonic potential. The energies are in units of ω and the ones with g = 2 correspond to a vertical cut in Fig. 1 panel (b), s = 0.5. N = 3 and N = 4 bosons compared with a simple meanfield variational ansatz using the following wave function, and finding the optimum α * that minimizes the energy As expected, this mean-field ansatz captures well the behaviour of the ground state of the system for small values of g. For g ≃ 2, however, we already observe substantial deviations, with the meanfield prediction overstimating the ground state energy considerably. In particular, as we will see below, the system develops strong beyond-mean-field correlations as g is increased.
In addition, we introduce a two-body-correlated variational many-body ansatz of Jastrow type [35], close, and some times even below, the exact diagonalization ones. To improve the latter, one needs to enlarge the Hilbert space (larger M ) to get a slightly lower upper bound. In principle, the exact diagonalization procedure for a given M provides an upper bound for the ground state and each excited state. In the next section, we explain the physical interpretation of the variational parameters and discuss how well the ansatz captures the physics of the problem. The low-energy spectrum for N = 3 and N = 4 at smaller values of g is fairly similar. This is not unexpected as the degenerate manifolds are the same irrespective of the number of particles, see Sect. III B. The first excited state is a center-of-mass excitation, the Kohn mode, as seen clearly in the excitation spectra shown in Fig. 3.
Even for g up to 16, the low-energy spectra for N = 3 and N = 4 are quite similar. The overall picture is qualitatively the same for both cases, although for N = 4 there are extra levels crossing. In Fig. 3 panel (b), there is a level that starts crossing the highest energy level depicted at g ≃ 3. This line in the spectrum comes from the fourth excited level in the noninteracting limit and is also expected to appear for systems with more particles, e.g. N = 5. It arises from the existence of a degenerate kind of states that are found only for N ≥ 4, as it is explained in Sect. III.

D. Degeneracy for the interacting three and four-boson systems
One major difference for more than two particles, is that we do not find states not affected by the interaction. Moreover, the degeneracy of the eigenfunctions of the relative part of the Hamiltonian is not 1. Therefore, the states cannot be uniquely characterized by ν r . However, we can identify the states that are center-of-mass excitations of lower energy states. In Fig. 3, in both panels, for example, for g = 1, we know the degeneracy of all the energy levels and we can identify them. The ground state is nondegenerate. As we have said before, the first excited state is a center-of-mass excitation, with degeneracy 2. The second excited state decomposes in three states corresponding to the next center-of-mass excitations of the ground state, there are two degenerate states corresponding to a relative excitation, and finally a different relative excitation. The third excited energy level in the noninteracting limit splits when g is increased in the next center-of-mass excitations of the states of the previous level, i.e., four center-of-mass excitations of the ground state, four center-of-mass excitations of the previous two-degenerate relative excited states, and two more degenerate states corresponding to two center-ofmass excitations of the single-degenerate relative energy level that appeared in the second excited state when g was increased. Moreover, there are two pairs of different relative excited states that split from the noninteracting third energy level. This behaviour is the same independently of N for g sufficiently small, for instance, for N = 4 up to g = 3, where there is the previous discussed crossing of levels.

E. N -boson energies up to first order in perturbation theory
Using the analytic expressions of the integrals of the interaction that are given in Appendix C, we compute the energies of the first three energy levels in first order perturbation theory. For the ground state of the system, the energy is given by The next level has energy The third energy level splits in three, in the way that is discussed in the previous section that is also valid for N bosons, with energies The similarity in the energy difference, E − E 0 , for the case of N = 3 and N = 4 plotted in Fig. 3 for a small g can be understood using the previous expressions. The corresponding excitation energies are, in this approximation, In the first two cases, Eq. (39) and Eq. (40), we recover the first and the second center-of-mass excitations that are red solid lines in Fig. 3. The presence of the factor N in the quantity E 22 − E 0 , see Eq. (41), explains why the slope of the green dashed lines is slightly bigger in absolute value for N = 4, panel (b), than for N = 3, panel (a), in Fig. 3 for g ≃ 0. This effect would be notorious when comparing the spectrum for two very different numbers of particles. Finally, we also see that the second term in E 23 − E 0 is proportional to N , but in that case, for s small, the second term becomes negligible. Therefore, the blue dotted lines are very close to the red solid lines in the spectra for g ≃ 0, as we have used s = 0.5. In the zero-range limit, this approximation gives E 23 (s → 0) = E 21 (s → 0). As the perturbative correction affects only the relative motion, the corrections to E 0 , E 1 and E 21 are equal.

V. INTERACTIONS AND QUANTUM CORRELATIONS
As seen in the previous section, the ground state energy of the system for N = 2, 3 and 4 tends to saturate as we increase the strength of the atom-atom interactions. This saturation starts to occur for values g for which the mean-field variational ansatz starts to deviate from the exact results. This reminds of a similar effect found in 1D systems, where the ground state evolves from meanfield to Tonks-Girardeau gas as the interaction strength is increased. In the Tonks-Girardeau limit, the atoms do avoid completely the atom-atom contact interaction by building strong correlations which in 1D are easily understood from the Bose-Fermi mapping theorem [36]. In 2D, no such mapping exist. However, we expect that the system should build suitable correlations to avoid the interaction, resulting in a saturation of the energy for increasing g.
For the ground state, besides the exact diagonalization method, we have also made use of a correlated variational ansatz, Eq. (33), to enlighten the discussion. The energies and properties associated to this variational ansatz are evaluated by means of Monte-Carlo methods (standard Metropolis algorithm). The physical meaning of the variational parameters is quite transparent. α directly affects the overall size of the cloud. The two-body Jastrow correlations are parameterized by a and b. Two limiting cases are illustrative. If the system is fully condensed we will have a = 0, while a = 1 would correspond to building a zero of the wave function whenever two atoms are at the same position. b affects the two-body correlation length. Thus, we expect the following behavior: for val- ues of g ≃ 0 we should have a = 0 (b is thus irrelevant) and α close to 1. For increasing g, α decreases to avoid the interaction by simply putting the atoms apart. As we increase g, two-body correlations build in, a = 0 and α should stop decreasing as the correlation is more efficient to separate the atoms. Let us first discuss the density profile of the clouds, see Appendix B for definitions. In Fig. 4 we show the density profile, normalized to unity, depending on the radial coordinate X = x 2 + y 2 , computed with our exact diagonalization procedure. Due to the symmetry of the trap, the density profile of the ground state does not have angular dependence, see Appendix B, Eq. (B6). In panels (a), (b) and (c) we show results for N = 2, 3, and 4. In all cases, with the same value of s = 0.5. We compare densities obtained for different values of g.
Irrespective of N we observe a number of common features. For g = 0, the system has a Gaussian density profile which, as g is increased, evolves into a profile with a flat region for X ≤ 1 at g ≃ 16. As N is increased, the size of the inner plateau increases, thus tending towards an homogeneous density.
The quality of our variational approach is seen in Fig. 5. We compare density profiles obtained with the exact diagonalization procedure with those obtained variationally by means of Eq. (33). As seen in the figure, the variational wave function provides a fairly accurate representation of the density profile for N = 2, 3, and 4. In particular, it captures well the appearance of the plateau.
The effect of increasing the interaction among the atoms is manifold. As we have seen above, the density profile is modified and the gas becomes close to homogeneous in the inner part of the trap. This change in the density is however accompanied by a change in the correlations present in the system. Actually, the gas goes from a fully condensed state to a largely fragmented one as we increase the interaction. In Fig. (6), we depict how the condensed fraction for N = 2, 3, and 4 decreases when increasing the interaction strength. For the same value of g, the fragmentation in the system is larger for larger number of particles. The most populated eigenstate of the one-body density matrix (natural orbit), is found to have the approximate form, using the |n x , n y basis, and its wave function reads This natural orbit is a superposition of the two first single-particle states of the 2D harmonic oscillator with zero angular momentum, m = 0, the state |n = 0, m = 0 and the state |n = 2, m = 0 , thus the wave function has no angular dependence. For the noninteracting case, C 0 = 1 and C 1 = 0, since the particles condense in the ground state of the harmonic oscillator. When the interaction is increased, C 0 becomes smaller and C 1 starts to increase. In Fig. 7, we plot the wave function of Eq. (44) using the corresponding values of C 0 and C 1 computed for N = 2, 3, 4 and different values of the interaction strength g. The advent of correlations beyond mean-field ones should also become apparent when computing twoparticle correlations. In particular, we can evaluate the probability of finding two particles at given positions. For simplicity we consider one of them at the origin and the second one at a distance X. The probability density of finding a particle in the space once we have fixed a particle at the center is given by η(X)/ρ(0) and is normalized to unity (see Appendix B). Without interactions, the pair correlation function is proportional to the density, since the probability density for finding a particle in a particular place is not correlated with the positions of the others, see Eq. (B20). In Fig. 8, we show how η(X)/ρ(0) evolves with increasing the interaction for the systems with N = 2, 3, and 4 bosons. In all cases, the central peak gets smaller when increasing the interaction, being fairly close to zero for g ≃ 16. This is in line with the fact that the atoms build correlations to avoid the interaction, e.g. as g is increased the probability of finding two atoms at the same location decreases. In between, next to the center of the trap, the function is uniform. When the in-teraction is strong there is a minimum at the position of the first atom, the probability density η(X)/ρ(0) develops a maximum corresponding to the preferred distance between particles. Increasing the number of bosons, this maximum shifts towards larger distances.

VI. SUMMARY AND CONCLUSIONS
In this work, we have studied systems of a few number of bosons trapped in an isotropic 2D harmonic trap interacting by a finite-range Gaussian potential.
First, we have explored in detail the noninteracting case, paying particular attention to the degeneracies of the excitation spectrum of the system. In particular, for the N -boson case, we have explained how to compute the degeneracy of the low-energy states which is independent of the number of particles.
By means of a direct diagonalization of the Hamiltonian in a truncated space, we have studied the interacting system and we have computed the low-energy spectra for N = 2, 3, and 4 bosons. We have also proposed a variational ansatz with two-body correlations which provides an accurate description of both the energy and the structure of the ground state in the full range of interaction considered. Center-of-mass and relative excitations are clearly identified in the spectrum. As the interaction is increased, we have shown how the ground state and all low lying states tend to saturate as a function of the interaction strength.
The effect of increasing the interaction on the ground state is twofold. On one side, the density at the center of the trap decreases becoming almost flat in the bulk of the gas, with the cloud thus becoming larger. On the other side, the atoms develop strong two-body correlations to avoid the interaction. This is achieved by building holes in the many-body wave function whenever two atoms are at the same position, as is clearly seen in the computed pair correlations and also on the explicit zeros introduced in our variational wave function. This mechanism is similar to the one present in the Tonks-Girardeau gas in 1D and is also responsible for the observed saturation of the energies of the system as we increase the interaction strength. Finally, the onset of correlations in turn produces fragmentation on the one-body density matrix, which has been shown to increase with the number of particles.
Thus, we have shown that our exact diagonalization method allows one to study interacting bosonic systems in 2D. We are presently implementing this method for spin-orbit coupled bosonic systems.
for sending his notes on this topic. We also would like to show our gratitude to N. L. Harshman for discussions about degeneracy and a careful reading of the manuscript. The authors acknowledge financial support by grants 2014SGR-401 from Generalitat de Catalunya and FIS2014-54672-P from the MINECO (Spain). P.M. is supported by a FI grant from Generalitat de Catalunya and B.J.-D. is supported by the Ramón y Cajal program.
Appendix A: Computation of degeneracies in the noninteracting limit

The two-boson system
We compute the degeneracy of each energy level depending on the excitation energy number, N E = E − E 0 , for the two-boson system with the possible states labelled using the quantum numbers of Eq. (20). First, we fix the excitation energy number, N E , and consider it to be even. Then, the values that n r can take are n r = 0, 2, ... , N E , so n r = 2k with k = 0, 1, ... , N E /2. Since we have n cm + n r = N E , for each value of n r there is the corresponding n cm . Now, we count the number of states with a given n r with excitation energy number N E taking into account the degeneracy due to the quantum numbers m cm and m r , that is, Therefore, we have to sum over k to find the degeneracy. The sum goes from k = 0 to k = N E /2 if N E is even and to k = (N E − 1)/2 if N E is odd, which can be generalized using the floor function, summing from k = 0 to k = ⌊N E /2⌋. The degeneracy is The previous equation, Eq. (A2), for N E even is and for N E odd is For the spatial fermionic states, which are the ones with m r = odd and antisymmetric upon exchanging particles 1 and 2, we compute the degeneracy analogously, using that n r = odd, a. Unperturbed energy states We are also interested in knowing the number of states in each energy level with m r = 0. We compute this number of states subtracting from the total number of degenerate states, d b NE , the ones with m r = 0, that is, where we have used Eq. (A2). As before, we can separate the case with N E even, and the case with N E odd, Appendix B: Computation of the density profile, the pair correlation function and the condensed fraction 1. The density profile a. First-quantized density operator For a system of N particles, the density operator in first quantization, normalized to unity, is defined aŝ Therefore, the density profile for a given state of a system of N identical bosons, Ψ( x 1 , ... , x N ), would be In particular, for a two-boson system in 2D, the previous equation reduces to We compute the density profile for the general interacting case, in the harmonic trap, for the ground state of the system as (B4) where we have made use of the explicit form of the manybody wave function of the ground state, This way of writing the wave function of the ground state is equivalent to separate the center-of-mass part from the relative part. Using the change of variables r = x − x 2 and polar coordinates in Eq. (B4), we express the density as We have used that 2π 0 dϕ e A cos ϕ+B sin ϕ = 2π I 0 where I 0 is a modified Bessel function. Notice that, as we would expect, in Eq. (B6) we have demonstrated that the density only depends on the radial coordinate X ≡ x 2 + y 2 , and we can rewrite that equation as This result is valid not only for our Gaussian-shaped potential but also for any potential dependent only on the modulus of the relative coordinate. In these other cases, the explicit form of the interaction defines the relative wave function f (r). In the noninteracting case, we can compute the integral analytically, by substituting the explicit form of f 0 (r), and we recover the known result, where ϕ 0 (X) is the wave function of the single-particle ground state of the 2D harmonic oscillator. The previous result, ρ 0 (X) = |ϕ 0 (X)| 2 , is also valid for the case of N noninteracting bosons in the 2D harmonic potential, since the many-body wave function factorizes, Ψ 0 ( x 1 , x 2 ..., x N ) = ϕ 0 ( x 1 ) ... ϕ 0 ( x N ). We recover the previous result replacing the factorized wave function into Eq. (B2).
b. Second-quantized density operator For our numerical computations, we make use of the second-quantized form of the density operator, For a state written in our Fock basis, Eq. (11), as where the index k labels each state of the basis, |k = |n 1 , ... , n M , the density profile is computed as where ψ i ( x) are the single-particle eigenstates of the 2D harmonic oscillator.

The pair correlation function
The pair correlation operator, normalized to unity, for a system of N particles readŝ from which we obtain the pair correlation function for a state of the N -boson system, Ψ( x 1 , ... , x N ), as For the particular case of the ground state of two bosons in 2D, we have where Ψ ( x, x ′ ) is the corresponding wave function, Eq. (B5). For the noninteracting case, in the harmonic trap, we know the function of the relative part, Eq. (B9).
In that case, the pair correlation function is The last result is also valid for the system of N bosons, because then we can factorize, Ψ 0 ( x 1 , x 2 ..., x N ) = ϕ 0 ( x 1 ) ... ϕ 0 ( x N ), and replace the wave function into Eq. (B15) in order to find the same result. Now, we fix one particle at the origin, and compute the function (B18) Notice that the previous function depends only on the radial coordinate X ≡ x 2 + y 2 , so we can write Again, for the noninteracting case we have an analytical expression for the previous function, that reads and is proportional to the density, Eq. (B10). The probability density of finding a particle in the space once we have found a particle at the origin is given by the quantity η(X)/ρ(0). We verify its normalization to unity in the general case, where we have used that all the particles are identical, Eq. (B2) and Eq. (B15).

a. Second-quantized pair correlation operator
The second-quantized form of the pair correlation operator iŝ (B22) For a state written in our Fock basis, Eq. (11), as where the index k labels each state of the basis, |k = |n 1 , ... , n M , the pair correlation function is computed as η( x, x ′ ) = 1 N (N − 1) M i,j,p,q=1 where ψ i ( x) are the single-particle eigenstates of the 2D harmonic oscillator.

The condensed fraction
The degree of condensation is characterized using the one-body density matrix, where, i, j = 1, ... , M . Diagonalizing this matrix, its eigenvalues n i are computed, which are the occupations of the corresponding singe-particle eigenstates |φ i . The state |Ψ is fully condensed when |Ψ = |φ 1 ⊗N and then, the one-body density matrix has only a single nonzero eigenvalue, n 1 = 1. If there is fragmentation in the system, the highest eigenvalue n 1 < 1, due to the normalization, M i=1 n i = 1.

4A
is a confluent hypergeometric function of the second kind that we have expressed in series and Q ∈ N is defined as The next step is computing the integral in Eq. (C5) by replacing the explicit form of I x (x ′ ), Eq. (C8). First, we notice that depending on the parity of the integrand, the integral will be zero since we integrate in a symmetric interval. The possible situations are I xx ′ = 0 n x1 + n x2 + n x3 + n x4 odd I xx ′ = 0 n x1 + n x2 + n x3 + n x4 even.
In the second case, we compute the integral replacing again the Hermite polynomials by their series representation and substituting (C8) into (C5), The expression is analogous for I yy ′ and all the sums that appear are finite and have few terms when n xi are small. Now, knowing the form of I xx ′ and I yy ′ we have V i,j,k,l . Moreover, many of the integrals are zero V i,j,k,l = 0 4 i=1 n xi odd or 4 i=1 n yi odd V i,j,k,l = 0 4 i=1 n xi even and 4 i=1 n yi even, (C14) and we also take profit from the symmetries of I xx (n x1 , n x2 , n x3 , n x4 ), which verifies I xx (n x1 , n x2 , n x3 , n x4 ) = I xx (n x4 , n x2 , n x3 , n x1 ) = I xx (n x1 , n x3 , n x2 , n x4 ) = I xx (n x4 , n x3 , n x2 , n x1 ).
(C15) Therefore, we are computing four integrals at the same time.