# Use of two-body correlated basis functions with van der Waals interaction to study the shape-independent approximation for a large number of trapped interacting bosons.

###### Abstract

We study the ground state and the low-lying excitations of a trapped Bose gas in an isotropic harmonic potential for very small () to very large () particle numbers. We use the correlated two-body basis functions and the shape-dependent van der Waals interaction in our many-body calculations. We present an exhaustive study of the effect of inter-atomic correlations and the accuracy of the mean-field equations considering a wide range of particle numbers. We calculate the ground state energy and the one-body density for different values of the van der Waals parameter . We compare our results with those of the modified Gross-Pitaevskii results, the correlated Hartree hypernetted-chain equations (which also utilize the two-body correlated basis functions), as well as of the Diffusion Monte Carlo for hard sphere interactions. We observe the effect of the attractive tail of the van der Waals potential in the calculations of the one-body density over the truly repulsive zero-range potential as used in the Gross-Pitaevskii equation and discuss the finite-size effects. We also present the low-lying collective excitations which are well described by a hydrodynamic model in the large particle limit.

###### Keywords:

Bose-Einstein condensate Many-body physics Collective excitations Shape-independent approximation###### pacs:

03.75.Hh 31.15.Ja 03.65.Ge 03.75.Nt^{†}

^{†}journal: J Low Temp Phys

∎

## 1 Introduction

Laboratory realization of gaseous Bose-Einstein condensates (BEC) Anderson (); Bradley (); Davis () and subsequent experiments characterizing low-lying collective excitations Jin (); Jin1 () have prompted various theoretical investigations of these systems Stringari (); Dalfovo1 (); Hu (); Andras (); Fil (); JP (). These theoretical calculations include the trapped Fermi gas and BEC in a shallow trap. Unlike phenomena such as superfluidity of liquid helium, the atomic vapor is very dilute where the fundamental interactions are characterized by the -wave scattering length . With this type of interactions the system is quite easy for theoretical understanding. The standard theory for BEC in a dilute atomic vapor uses the Gross-Pitaevskii (GP) equation Dalfovo () which is a nonlinear Schrödinger equation with the inter-atomic interaction characterized by the -wave scattering length . The exact shape of the inter-particle interaction and inter-atomic correlation are ignored in this picture. The elementary excitations of a BEC in a harmonically confined dilute Bose gas have been studied using the GP equation Dalfovo2 ().

The GP theory is the most popular tool for the description of weakly interacting bosons. However, in recent experiments the number of trapped atoms varies from just a few to . Therefore, it has become imperative to study the effect of realistic interaction, inter-atomic correlation and the accuracy of the GP equation. The astonishing success of the mean-field theory lies in the fact that the gas parameter ( is the atom number density) is small. However, nowadays, utilizing the Feshbach resonances, can be tuned by changing the magnetic field. Particle-particle correlation becomes important at large and the accuracy of the GP theory needs a deeper study. The mean-field description uses an effective mean-field potential obtained by the shape-independent pseudopotential approximation (SIA). The SIA implies that the calculated ground state energy remains unchanged irrespective of the shape of the interaction potential. As the hard sphere interaction is completely characterized by a single parameter, it is beyond the capacity of characterizing the universal behavior of the ground state properties. Instead, we utilize a realistic interaction having a controllable parameter .

The SIA has been addressed in different context. For homogeneous systems, Cowell et al Cowell () have shown that the SIA fails when , i.e., different potentials having the same scattering length leads to different ground state energies. In the case of an inhomogeneous system with few atoms in the trap, the SIA is less valid for tight confinement. For a homogeneous Bose gas, Giorgini et al. Giorgini () find a small dependence of the ground state energy on the exact shape of the two-body potential when . For larger or for stronger confinement the SIA becomes less applicable. The validity of the SIA has also been addressed for both weakly and strongly interacting BECs. The ground state energy of a trapped BEC is calculated by the diffusion Monte Carlo (DMC) approach using different two-body potentials that generate identical Blume (). It is seen that different potentials produce indistinguishable total ground state energies for a small gas parameter. Whereas for larger , inter-atomic correlations play an important role and SIA becomes invalid and quantum corrections to the ground state energy are required Blume (). The local density versus correlated basis approaches for trapped bosons have been examined and the beyond the GP approximation has been prescribed Poll (); Poll1 ()

In the present manuscript, as stated before, we report the ground state properties of weakly interacting systems with a wide range of using the two-body correlated basis functions. Although the ground state properties have been addressed in different theories as described above, we do not find any exhaustive study which keeps the effect of inter-atomic correlations, uses the realistic van der Waals interaction and treats the real experimental situation where the number of particles is finite, varying from very few to a quite large number. Thus, the first key question of our study is to test the accuracy of the mean-field theory for large and to study the SIA for a wide range of . As the modified GP (MGP) equation accounts for quantum fluctuations, we also check how closely our correlated basis functions reproduce the ground state energy obtained with the MGP. We also compare our calculated results with the correlated Hartree HNC results Poll () which utilizes hard sphere bosons. We observe that the MGP substantially improves the GP results and our correlated many-body energies are in good agreement with the MGP. All the earlier calculations in this direction consider a truly finite number of particles and use the standard short-ranged two-body potential. Here we use the realistic van der Waals potential with a short-ranged hard core and a long-ranged attractive tail characterized by the parameter . Comparison with the DMC calculation is also made for few bosons with a repulsive hard sphere interaction to justify the accuracy of our two-body correlated basis function in the dilute regime since the DMC is essentially exact. Thus, we are further assured to tackle real experimental situations which consider quite a large number of bosons and, hence, is beyond the scope of the DMC.

We investigate the effect of the long attractive tail (by changing , adjusting the cutoff radius , such that they produce the same scattering length) on the ground state energy for a wide range of . For , we observe good agreement between our many-body results and the mean-field results. This establishes the applicability of the SIA for a wide range of . In addition to the ground state energy, we also calculate the effect of the long-range attractive tail of the van der Waals interaction on the one-body density for a small and a large number of atoms in the trap. Although we observe that the one-body density profile is almost independent of the choice of the long-range attractive tail of the van der Waals potential, the peak value of the density substantially differ from the mean-field results. We next analyze the finite size effect on the one-body density profile.

Studies of collective excitations at both low and high energies at the large- limit is quite interesting for the following reason. The excitations at high energies are expected to be of single-particle nature. However, a laboratory BEC is a strongly inhomogeneous system due to the external trap and may significantly differ from the uniform Bose gas at low energies where only phonons are present. Thus the study of the transition from low-lying collective excitations to single-particle excitations by the correlated many-body method using the realistic interaction is itself interesting. In the present many-body calculation, we can in principle calculate all the low-lying and high-lying collective excitations. However, as mentioned in the formalism in the next section and in the result section, we solve the coupled differential equations by the hyperspherical adiabatic approximation (HAA) and consider the lowest eigen potential as the many-body effective potential. HAA with the lowest eigen potential is very accurate for the ground and few low-lying collective excitations. However, for the high-lying excitations the effect of higher order eigen potential will come in the picture and one should use the coupled adiabatic approximation. Thus for our present calculation we report only on the few low-lying excitations where the effect of higher order eigen potentials can be safely ignored.

The paper is organized as follows. In Sec. 2 we discuss the correlated potential harmonic basis for a large particle number. Sec. 3 mainly considers ground-state properties of the systems and comparison with the mean-field results is provided. The validity of the shape-independent approximation is presented in Sec 4. Sec. 5 deals with the calculation of collective excitations at low energies. Sec. 6 concludes with a summary.

## 2 Formalism

### 2.1 The correlated potential harmonic expansion method

In the present work we calculate the ground state energy and the low-lying collective excitation frequencies of a dilute BEC for a large number of bosons () using the potential harmonics expansion method (PHEM) Fabre83 (). The most fundamental feature of the PHEM is that the two-body correlations are dominant in the many-body system and the two-body Faddeev component of the -body wave function is a function of the two-particle relative separation and a global length called hyperradius . We have already successfully utilized the two-body correlated basis functions for the description of dilute BECs TKD (). So in the present description we only point out the essential components of the PHEM for the completeness and clarity of the manuscript. For a detailed formulation we refer the reader to Refs. Das1 (); Das2 (); PLA2009 ().

In the PHEM we expand the two-body Faddeev component corresponding to the interacting pair of bosons, in a condensate with bosons, in terms of the potential harmonic (PH) basis as

(1) |

where corresponds to the full set of hyperangles in the -th partition while and are the orbital angular momentum of the system and its projection. In the hyperspherical coordinate, the variables are characterized by the hyperradius (, being the Jacobi vectors describing the relative motion) and hyperangles Ballot (). However, for the potential harmonic expansion method of a weakly interacting BEC, as the only the two-body correlations are dominating, we assume that when the pair of atoms interact, the rest of the atoms are inert spectators Das1 (). As all the degrees of freedom coming from the inert spectators are frozen, the number of quantum numbers becomes effectively four irrespective of the number of bosons. These are the orbital angular momentum , the azimuthal , the grand orbital angular momentum , and the energy quantum number. The closed analytic expression for the PH, can be found in Ref. Fabre83 (). It is indicated in our earlier works that the above expansion is in general very slow as the lowest order PH is a constant and does not represent the strong short-range repulsion of the inter-atomic interaction. Therefore, we introduced a short-range correlation function which is obtained as the zero-energy solution of the two-body Schrödinger equation

(2) |

with the chosen two-body potential PLA2009 (), and corresponds to the appropriate -wave scattering length as described in next section. After including the correlation function in the PHEM basis, we call it the correlated potential harmonic expansion method (CPHEM). The expansion Eq. (1) now takes the form

(3) |

Substituting this expansion (Eq. (3)) in the many-body Schrödinger equation, one gets a set of coupled differential equations (CDE) in . The coupling potential matrix element is given by

(4) | |||||

where , and are the Jacobi polynomial, its norm and weight function respectively, with and . The CDEs are solved using the hyperspherical adiabatic approximation (HAA) HAA (). The HAA basically reduces the whole dimensional problem to an effective one-dimensional problem. In the HAA, the coupled potential matrix along with the diagonal hypercentrifugal repulsion is diagonalized to get an effective potential as the lowest eigen value of the matrix for a particular value of .

The basic length scale for a harmonic oscillator trap of frequency is . For the typical experimental BEC, is of the order of . However, the effective potential in hyperspherical space due to the hypercentrifugal repulsion together with the harmonic oscillator trap has a minimum at about . As an example, with , the minimum of the effective potential will be near , which is almost times larger than the typical range of the inter-atomic interaction. This shows that for such a typical case, the entire contribution to in the integral in Eq. (4) comes from an extremely narrow interval of -integration (). The integral in Eq. (4) also varies rapidly within this narrow interval because of the following reason. The integrand contains the Jacobi polynomial and its weight function Abram (). For large , both and change very rapidly with respect to . varies from zero at to a maximum of at and then rapidly reaching a value about of the peak value at . Although the peak value is extremely large for large , partial cancellation results from the factor Abram (). Thus, any standard quadrature to evaluate the integrand in Eq. (4) gives essentially zero for . Usually we solve this problem by splitting the interval into gradually increasing subintervals and evaluating the integral in each subinterval using a 32-point Gauss-Legendre quadrature. This permits us to evaluate for up to 15000 with an accuracy of one part in TKD-PRA-2007 ().

### 2.2 Extension to

As pointed out earlier, the experimental BEC treats up to atoms in the trap. But the numerical code mentioned above can treat only up to 15000 atoms which is far from the experimental situation. To circumvent the problem and extend the correlated many-body technique to quite a large number of atoms, we recently made a direct mathematical transformation Sofianos (), which transforms the PHEM into a two-variable integro-differential equation. With this transformation, the Jacobi polynomial is replaced with the associated Laguerre polynomial. In our initial attempt we applied the CPHEM using the Laguerre polynomial (CPHEL) for the order of atoms for the ground state Sofianos (). We utilize a mathematical relation to transform the Jacobi polynomials into the associated Laguerre polynomials Abram (). An outline of the derivation, including derivation of the relations between the Jacobi and the associated Laguerre polynomials in the limit are as follows. Starting from the mathematical relation Abram ()

(5) |

interchanging and , and using the relation Abram ()

(6) |

we obtain

(7) |

Substituting and , we get

(8) |

This relation was used in evaluating appearing in the CDEs for large Sofianos (). In this limit, the weight function of the Jacobi polynomial transforms as

(9) |

In the limit , the last factor becomes . Hence for large ,

(10) |

This has the correct functional form for the weight function of the associated Laguerre polynomial . Substituting equations (8) and (10) in equation (4), and using the explicit expression of the norm of the Jacobi polynomial Abram (), we obtain

(11) |

where , is the hard-core radius of our chosen realistic van der Waals potential and

(12) | |||||

where .

## 3 The ground-state energy for atoms in an external trap

Throughout our calculation we keep the system parameters fixed at values which correspond to the JILA trap Anderson (). The mass amu, the trap frequency Hz and the scattering length Bohr. However, for comparison with the DMC results with repulsive a hard sphere interaction, we also compute the ground-state energy for few bosons with a larger scattering lengths Bohr and Bohr. As a unit of length we choose the oscillator unit (o.u.) of length and the energy unit as the harmonic oscillator energy . For the mean-field GP equation, the two-body potential is chosen as the zero-range potential Dalfovo (). This potential is shape-independent and completely ignores the dependence of energy on the scattering amplitude. We choose the realistic van der Waals interaction

(13) |

with o.u. for Rb atoms Pethick (). For a given value of , is calculated by looking at the zero-energy solution of the two-body Schrödinger equation Eq. (2), where is the van der Waals potential Pethick (). Its asymptotic form quickly attains from which is determined Pethick (). We choose the value of o.u. which corresponds to . It is to be noted that although is almost four times smaller than , they are of same order. Smaller values of (corresponding to a larger number of nodes in ) are not chosen to avoid the presence of the many-body bound states and clustering. It should also be noted that for hard-sphere scattering, the effective range increases linearly with , as , whereas for the van der Waals potential is determined from by Gao ()

(14) |

where and . The calculated value of for our present work is o.u. which is comparable with but greater than the value of , as expected.

As pointed earlier, the lowest eigen potential is treated as the many-body effective potential which describe the collective phenomena of a dilute BEC. This is also in good agreement with experimental situation. As at the ultra-cold temperature all the individual atoms in the condensate lie within a single de-Broglie wavelength, the condensate is treated as a single lump of quantum stuff. However, before calculating the ground-state energies, it is indeed required to check the convergence for atoms. With the increase in the number of particles , the effective interaction increases, the condensate becomes more repulsive. The condensate density is pushed out as its average radius increases sharply. So higher may be needed for convergence of the ground state for larger . In Fig. 1 we plot as a function of for and with and observe a very fast convergence. Although after , all the graphs appear to overlap completely as seen in the Fig. 1, still we noticed that the minimum of the effective potential decreases very slowly (which is not visible in the figure) as increases. This is consistent with the Rayleigh-Ritz principle. So, throughout our calculation we fix and calculate the ground state energy per particle. In Table 1, we compare our many-body results with the TF, GP and MGP results for different diluteness of the condensate given by . We evaluate the gas parameter at the center of the trap which can be directly expressed in terms of relevant parameters as Dalfovo ()

(15) |

Note that the value of is deliberately kept to justify the usage of the two-body correlated basis functions. Comparison with the MGP is needed for better justification as the MGP includes the correction due to quantum fluctuations.

TF | GP | MGP | CPHEL | HNC | ||
---|---|---|---|---|---|---|

1.90 | 2.42 | 2.43 | 2.43 | 2.43 | ||

4.76 | 5.04 | 5.08 | 5.19 | 5.04 | ||

11.96 | 12.10 | 12.25 | 12.67 | 12.20 | ||

30.05 | 30.12 | 30.66 | 31.67 | 30.48 | ||

75.49 | 75.52 | 77.48 | 79.48 | 76.85 |

We calculate the GP energy by solving the standard GP energy functional

(16) |

and the MGP energy is calculated by solving

(17) | |||||

The additional term in Eq. (17) basically adds quantum corrections to the mean-field effective potential. The TF energy is calculated by using a simple analytic expression Dalfovo (). All the results are presented in Table 1.

We observe that our many-body ground state energy is almost indistinguishable from those of both the GP and MGP for small , as expected. However, for larger , we observe that our many-body results are consistently closer to the MGP results. Note that the TF results are always lower than both the GP and many-body results as the kinetic energy term is completely ignored in the TF limit. We calculate the correlation energy PRA2014 ()

as a measure of the small deviation of the mean-field GP and the TF results from our many-body results.

We present our results in Fig. 2 and observe that both and converge to the same very small value at the very large particle-number limit. This reaffirms that the ground state properties of a dilute BEC in the truly thermodynamic limit should be correctly described by the TF equation. In the same table (Table 1), we also compare our results with the correlated Hartree HNC results (solution of eq. (17) of Ref Poll ()) which utilizes hard sphere bosons. We find that our many-body energies are larger than the HNC results. The presence of a hard-core part in the van der Waals interaction produces an excluded volume. Therefore, for the same , our many-body method should produce larger clusters than the GP and the HNC, and hence, has a larger contribution from the harmonic trapping potential. Of course there is also the effect of the attractive tail of the van der Waals potential which tries to reduce the energy. However, in the dilute condition it does not have a significant contribution which can be detected experimentally. Therefore, the combined effect of a hard core and an attractive tail along with the presence of a centrifugal repulsion term in the coupled differential equations make the CPHEL ground state energy more repulsive. As expected, this term will have a larger contribution for larger as shown in Table 1.

Before addressing the validity of the SIA for large , it is instructive to have an estimation of the accuracy of our two-body correlated basis function method. Thus we compare the CPHEM results with the available DMC results Blume (), which are essentially exact. For completeness we also include the GP and the MGP results. We have presented the total ground state energies for various numbers of bosons interacting via a repulsive hard-sphere potential both for weak and strong interaction in Table 2. For Bohr, the CPHEM results are in good agreement with the DMC results which guarantees the applicability of two-body correlated basis function in the dilute regime. It further assures us to carry forward the calculations for large where the DMC approach fails and our correlated basis functions approach works. It is also noted from Table 2 that the relative difference between the DMC and the CPHEM results is almost of the order of for Bohr. Although we are not interested in the strong interaction in the present work, we also report results for larger , for completeness. The relative difference between the DMC and the CPHEM results is again of the order of for Bohr while for = 10000 Bohr the relative difference is of the order of or more. We may conclude that for stronger interactions, the inclusion of higher body correlations may be required.

DMC | GP | MGP | CPHEM | ||
---|---|---|---|---|---|

100 | 3 | 4.510 | 4.510 | 4.510 | 4.506 |

5 | 7.534 | 7.534 | 7.534 | 7.531 | |

10 | 15.153 | 15.153 | 15.153 | 15.143 | |

20 | 30.640 | 30.638 | 30.639 | 30.622 | |

1000 | 3 | 4.603 | 4.600 | 4.602 | 4.575 |

5 | 7.835 | 7.826 | 7.834 | 7.775 | |

10 | 16.426 | 16.383 | 16.426 | 16.280 | |

20 | 35.475 | 35.293 | 35.497 | 35.180 | |

10000 | 3 | 5.553 | 5.329 | 5.611 | 5.055 |

5 | 10.577 | 9.901 | 10.722 | 9.471 | |

10 | 26.22 | 23.61 | 26.84 | 23.072 | |

20 | 66.9 | 57.9 | 68.5 | 50.678 |

## 4 The shape-independent approximation

To conclusively address the issue of validity of the SIA, a detailed study of the ground-state energies over a wide range of the parameter is required. Although some papers TKD-PRA-2008 (); Esry1 (); Hau (); Esry2 (); Blume () present interesting discussions of this issue, none of them considers the wide range of and none use the realistic inter-atomic interaction. Considering few tens of atoms in the trap, it is shown that the many-body results approach closer to the mean-field results when the number of particles in the trap increases. It can be intuitively understood that tthe system becomes more classical with increasing . In fact, the SIA is valid for a truly finite number of atoms (few tens) and in the extremely dilute condition. However, our present calculation starts with only a few atoms and goes up to the order of millions of atoms in the trap. This requires a deeper and thorough study of the validity of the SIA for such a wide range of . Using the van der Waals interaction instead of the hard-sphere interaction, it is very easy to study the effect of the long-range attractive tail in the ground-state properties of the BEC and to observe its universal behavior.

We tune the parameter and by changing the cutoff radius , we fix the scattering length to 100 Bohr, which mimics the JILA experiment. The calculated many-body results together with the GP and MGP results are presented in Table 3. Again we observe that the many-body ground state energy is almost indistinguishable for up to a few thousands of atoms. For larger , the ground-state energies deviate slightly and the relative change is almost negligible. To have an estimation, we consider the relative change for two extreme choices of . When the value of for these two choices was varied by a factor of , the ground-state energy for the largest choice of the number of atoms () varied by a factor of = . This indicates that the ground-state energy is almost insensitive to the exact value of for a wide range of within the given density parameter . Thus, the ground-state energy considerably satisfies the SIA.

GP | MGP | ||||||
---|---|---|---|---|---|---|---|

3 | 1.512 | 1.512 | 1.512 | 1.512 | 1.512 | 1.511 | 1.511 |

5 | 1.519 | 1.519 | 1.519 | 1.519 | 1.519 | 1.514 | 1.515 |

10 | 1.522 | 1.522 | 1.522 | 1.522 | 1.522 | 1.517 | 1.519 |

20 | 1.548 | 1.547 | 1.547 | 1.547 | 1.547 | 1.533 | 1.535 |

100 | 1.678 | 1.677 | 1.677 | 1.676 | 1.676 | 1.652 | 1.653 |

2.435 | 2.434 | 2.432 | 2.432 | 2.430 | 2.424 | 2.43 | |

5.199 | 5.198 | 5.194 | 5.188 | 5.182 | 5.08 | 5.08 | |

12.70 | 12.68 | 12.66 | 12.64 | 12.62 | 12.10 | 12.25 | |

31.75 | 31.68 | 31.64 | 31.60 | 31.56 | 30.12 | 30.66 | |

79.80 | 79.56 | 79.42 | 79.31 | 79.22 | 75.52 | 77.48 |

Next we compare the one-body density for various parameters and for a wide range of . Although we have reported some results on the one-body density for smaller and have observed an appreciable effect of finite size, in the present work we are interested in the thermodynamic limit. The one-body density is a key quantity as it contains information regarding one-particle aspect of the condensate and can be indirectly measured in the interferometry experiments. We define it as the probability density of finding a particle at a distance from the centre of mass of the condensate JCP2010 ()

(18) |

where is the full many-body wave function and the integral over the hypervolume excludes the variable . In Fig. 3 we present the calculated one-body density for atoms in the trap and compare with the GP results. The calculated one-body density for various values of perfectly match with the GP results which confirms our previous observations. In Fig. 4 we plot the one-body density for atoms (panel (a)) and for atoms (panel (b)) in the trap, and compared with the GP results. The density profiles calculated from the GP equation have the same qualitative features, however, disagreement remains in the peak value of the density distribution as well as in the extension of the density profile, which will be discussed later. We plot the enlarged profile of the one-body density near the peak [panel(c) for and panel(d) for ] and near the tail part [panel(e) for and panel(f) for ] for various parameters. It is seen that all the many-body results calculated for various and having the same scattering length are almost indistinguishable and may not be detected in an experiment. This means that in the dilute regime , the condensate is well described by the single parameter .

The disagreement between the many-body results and the GP results at the peak value needs additional discussion for which we plot the one-body density in Fig. 5, having the same value of the effective repulsive interaction , but different choices of and . The actual two-body attraction is determined by the integration . In a many-body calculation which uses the van der Waals potential having a long attractive tail , the net effective interaction is more attractive than in the GP case. Being more repulsive, the GP treatment lowers the central density and expands the density distribution. It is clearly seen that for the choice of and = 10000, the many-body results perfectly agree with the GP and MGP results at the peak value. Thus keeping constant, one can make the effect of repulsion stronger by appropriately choosing and . By gradually increasing , the effect of the repulsive interaction increases in our many-body calculation, the central peak gradually shifts downward and extended outward, while the mean-field results are independent of the separate choices of the number of atoms and scattering length. This supports the observed effect of using finite number of atoms in the many-body calculations, which has been discussed earlier in many contexts.

(a) | (b) |

(c) | (d) |

(e) | (f) |

## 5 Collective excitations at low energy

It was already pointed out that the low-energy collective excitations provide valuable information about the interaction, while the high-lying excitations are of single particle nature and are useful for the study of statistical properties. In our present picture, the collective motion of the condensate in the hyperradial space takes place in the effective potential . The ground state in this well gives the ground-state energy of the condensate corresponding to and . We use the notation for the energy (in o.u.) of the -th radial excitation of the -th surface mode. The monopole frequency is defined as the lowest hyperradial excitation corresponding to the breathing mode () and is calculated as . For we get the surface modes which can be calculated as hyperradial excitations in the eigen potential corresponding to different values of . However, large inaccuracy relating to the calculation of the off-diagonal potential matrix elements for and difficulties with the slow converge present challenges. For large , we observed that the diagonal hypercentrifugal term is very large which contributes most to the potential matrix. Hence we disregard contributions to the off-diagonal matrix elements and construct the effective potential in the hyperradial space for non-zero orbital angular momentum. In Fig. 6 we plot several breathing mode excitation frequencies like (monopole), (second breathing mode), (third breathing mode) and (fourth breathing mode) as functions of . In the description of the breathing mode frequencies at the large limit, it is important to calculate and compare frequencies of the different modes using the hydrodynamic (HD) model Stringari (); Dalfovo ()

(19) |

where and are the angular momentum quantum number and number of nodes in the radial solution, respectively.

(a) monopole | (b) second breathing mode |

(c) third breathing mode | (d) fourth breathing mode |

For quite a large , when the dimensionless parameter is large, the kinetic energy term in the ground state GP equation becomes negligibly small compared to mean-field term and one gets the TF approximation Dalfovo (). In the same limit the eigen frequencies are calculated using the hydrodynamic (HD) equation of superfluidity (Eq. (19)). So the comparison of the many-body results to the HD prediction is well justified. Note that the HD equation (Eq. (19)) depends only on and and not explicitly on . Thus the effect of finite size correction does not appear here and the many-body results should coincide with HD results in the true thermodynamic limit although small deviation may exist for finite size system. In Fig. 6, all the breathing mode frequencies saturate at the large limit. The asymptotic values for the several breathing modes are presented in Table 4 and compared with the HD results.

CPHEL | HD | |
---|---|---|

2.236 | 2.236 | |

4.468 | 3.742 | |

6.702 | 5.196 | |

8.936 | 6.633 |

It is seen that the HD prediction is very accurate for the description of the lowest excitation of large systems. However, we observe gradually increasing deviations as we go to higher modes. The slow but smooth increase in with increase in is visible in all four panels of Fig. 6, which basically manifests the finite size effect. The deviation from the HD results for larger can be attributed to the following reason. In the calculation of for we have taken the contributions from the diagonal matrix element and neglected the off-diagonal matrix elements as we face numerical difficulty as discussed earlier. Though their contribution is practically insignificant compared to diagonal part, however, this approximation may not be true for higher modes.

Another aspect of the many-body calculation is to calculate the high-lying excitations which are of single particle nature. However, as mentioned earlier, we calculate the excitation spectrum by using the eigen potential and for the present calculation we consider as the lowest eigen potential. This approximation is quite justified for the ground state and the low-lying excitations as shown in Fig. 7. For high-lying excitations, the effect of higher eigen potential may come in the picture and for the accurate calculation of the high-lying excitations we must use the coupled adiabatic channels. Therefore, for the present manuscript we leave the calculation of high-lying excitations to future work.

## 6 Summary and conclusion

We have calculated the ground-state energy, the one-body density and the low-lying excitation frequencies for a very large number of trapped bosons, which is close to the real experimental situation. We utilized the two-body correlated basis functions and used the van der Waals interaction as the inter-atomic interaction. The many-body method described in this manuscript can reveal realistic features of the trapped bosons. The most convenient and widely used tool in this direction is the mean-field GP equation which basically ignores inter-atomic correlation and uses the simple contact interaction. In that respect our many-body approach is a few steps ahead of the mean-field approach as by keeping all possible two-body correlations one can expect to address the beyond the mean-field effects. On the other hand the diffusion Monte-Carlo method is the exact many-body technique. However, due to computational difficulties it can handle only up to few hundreds of bosons in the trap, which is far from the real experimental situation. Our many-body method keeps only two-body correlations and can handle as many as atoms in the trap. The effect of only two-body correlations is relevant as the higher-body correlations are almost negligible in a dilute BEC.

In the first part of our calculation we applied the many-body approach for the calculation of ground-state energies. Comparison with the mean-field results through the correlation energies with respect to the GP and TF results demonstrate that the BEC looses its many-body effects and becomes more classical at a truly large-particle limit. This can be understood from the fact that for a large effective repulsion (large with ), particles are far apart from each other and the effect of interaction becomes small. The TF approximation well describes such a situation. Our present calculation also deals with the wide range of particle numbers and present an exhaustive study of the validity of the shape-independent approximation. However, another fundamental motivation of the present work is to study the collective excitations. We calculate several excited modes of breathing mode frequencies and compare with the HD model. We observe that the asymptotic value () of the lowest breathing mode (monopole frequency) exactly matches with the HD prediction. Whereas the higher breathing modes () have the same qualitative nature as , the asymptotic values are a bit higher than that of HD prediction. For smaller , the finite size effects exist and the breathing modes show slow and smooth increase until it approaches the asymptotic value (). We conclude that the low-lying collective excitations are well described by the HD model at . Our present calculation is an exhaustive study of both the static and dynamic behavior of trapped bosons. However, in the present work we strictly confine our attention to diluteness of the order of as we keep the effect of two-body correlations. Thus we leave the study of higher density regime to future works.

###### Acknowledgements.

BC acknowledges the financial support of Department of Science and Technology (DST), Govt. of India, under a Major Research Project [Sanc. No. SR/S2/CMP-126/2012].## References

- (1) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995)
- (2) C. C. Bradley et. al. Phys. Rev. Lett. 75, 1687 (1995)
- (3) K. K. Davis et. al. Phys. Rev. Lett. 75, 3969 (1995)
- (4) D. S. Jin, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 77, 420 (1996)
- (5) D. S. Jin et al, Phys. Rev. Lett. 78, 764 (1997)
- (6) S. Stringari, Phys. Rev. Lett. 77, 2360 (1996)
- (7) F. Dalfovo, S. Stringari, Phys. Rev. A 53, 2477 (1996)
- (8) H. Hu, G. Xianlong, and X. J. Liu, Phys Rev A 90, 013622 (2014)
- (9) A. Csordas and Z. Adam, Phys. Rev. A 74, 035602 (2006)
- (10) D. V. Fil and S. I. Shevchenko, Phys. Rev. A 64, 013607 (2001)
- (11) J. P. Martikainen, Phys. Rev. A 63, 043602 (2001)
- (12) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999)
- (13) F. Dalfovo, S. Giorgini, M. Guilleumas, L. Pitaevskii, and S. Stringari, Phys. Rev. A 56, 3840 (1997)
- (14) S. Cowell, H. Heiselberg, I. E. Mazets, J. Morales, V. R. Pandharipande, and C. J. Pethick, Phys. Rev. Lett. 88, 210403 (2002)
- (15) S. Giorgini, J. Boronat and J. Casulleras, Phys. Rev. A 60, 5129 (1999)
- (16) D. Blume and C. H. Greene, Phys. Rev. A 63, 063601 (2001)
- (17) A. Fabrocini and A. Polls, Phys. Rev. A 60, 2319 (1999)
- (18) A. Fabrocini and A. Polls, Phys. Rev. A 64, 063610 (2001)
- (19) M. Fabre de la Ripelle, Ann. Phys. (N.Y.) 147, 281 (1983)
- (20) B. Chakrabarti, and T. K. Das, Phys. Rev. A 81, 015601 (2010); A. Biswas, T. K. Das, L. Salasnich, and B. Chakrabarti, Phys. Rev. A 82, 043607 (2010); S. Goswami, T. K. Das, and A. Biswas, Phys. Rev. A 84, 053617 (2011); S. Bhattacharyya, T. K. Das, and B. Chakrabarti, Phys. Rev. A 88, 053614 (2013); S. K. Halder, B. Chakrabarti, T. K. Das, and A. Biswas, Phys. Rev. A 88, 033602 (2013)
- (21) T. K. Das, and B. Chakrabarti, Phys. Rev. A 70, 063601 (2004)
- (22) T. K. Das, S. Canuto, A. Kundu, and B. Chakrabarti, Phys. Rev. A 75, 042705 (2007)
- (23) T. K. Das, A. Kundu, S. Canuto, and B. Chakrabarti, Phys. Lett. A 373, 258 (2009)
- (24) J. L. Ballot, and M. Fabre de la Ripelle, Ann. Phys. (N.Y.) 127, 62 (1980)
- (25) T. K. Das, H. T. Coelho, and M. Fabre de la Ripelle, Phys. Rev. C 26, 2281 (1982)
- (26) M. Abramowitz, and I. A. Stegun, Handbook of mathematical functions, National Institute of Standards and Technology, USA (1964)
- (27) T. K. Das, S. Canuto, A. Kundu and B. Chakrabarti, Phys. Rev. A 75, 042705 (2007)
- (28) S. A. Sofianos, T. K. Das, B. Chakrabarti, M. L. Lekala, R. M. Adam, and G. J. Rampho, Phys. Rev. A 87, 013608 (2013)
- (29) C. J. Pethick and H. Smith, Bose-Einstein condensation in dilute gases, Cambridge University Press, England (2001)
- (30) B. Gao, Phys. Rev. A 58, 4222 (1998)
- (31) M. L. Lekala, B. Chakrabarti, G. J. Rampho, T. K. Das, S. A. Sofianos, and R. A. Adam, Phys. Rev. A 89, 023624 (2014)
- (32) B. Chakrabarti and T. K. Das, Phys. Rev. A 78, 063608 (2008)
- (33) B. D. Esry, Phys. Rev. A 55, 1147 (1997)
- (34) T. Haugset and H. Haugerud, Phys. Rev. A 57, 3809 (1998)
- (35) B. D. Esry and C. H. Greene, Phys. Rev. A 60, 1451 (1999)
- (36) A. Biswas, B. Chakrabarti and T. K. Das, J. Chem. Phys. 133, 104502 (2010)