Level spacing statistics and spectral correlations in the diffuse van der Waals clusters
Abstract
We present a statistical analysis of eigenenergies and discuss several measures of spectral fluctuations and spectral correlations for the van der Waals clusters of different sizes. We show that the clusters become more and more complex with increase in cluster size. We study nearestneighbour level spacing distribution , the level number variance , and the DysonMehta statistics for various cluster sizes. For large clusters we find that although the BohigasGiannoniSchmit (BGS) conjecture seems to be valid, it does not exhibit true signatures of quantum chaos. However contrasting conjecture of Berry and Tabor is observed with smaller cluster size. For small number of bosons, we observe the existence of large number of quasidegenerate states in lowlying excitation which exhibits the Shnirelman peak in distribution. We also find a narrow region of intermediate spectrum which can be described by semiPoisson statistics whereas the higher levels are regular and exhibit Poisson statistics. These observations are further supported by the analysis of the distribution of the ratio of consecutive level spacings which is independent of unfolding procedure and thereby provides a tool for more transparent comparison with experimental findings than . Thus our detail numerical study clearly shows that the van der Waals clusters become more correlated with the increase in cluster size.
pacs:
03.75.Hh, 31.15.Xj, 03.65.Ge, 03.75.Nt.I Introduction
Weakly bound fewbody systems are being studied since a long time back and have achieved revived interest recently as the physics of such weakly bound systems can be investigated experimentally in ultracold atomic gases (1). Utilizing the Feshbach resonance, the effective interatomic interaction can be changed essentially to any desired values (2); (3). The recent experiments on cold atoms also provide evidence of the existence of large weakly bound clusters. Thus our present study is motivated by the recent experiments on ultracold Bose gas. We treat the threedimensional bosonic cluster with maximum up to Rb atoms interacting through twobody van der Waals potential. Alkali atoms, specially Rb atoms, are good candidates for laser manipulation and to observe BoseEinstein condensate (4). At ultracold temperature the interatomic interaction is fairly well represented by a single parameter , the wave scattering length. For our present system we keep which corresponds to the JILA experiment (4). Thus the system is weakly interacting, and diffuse as the average size of the cluster increases with cluster size. The binding of such body cluster is provided by the twobody van der Waals potential having a short range repulsive core below a cutoff radius and a tail which represents the long range attractive interaction.
The stability of such body clusters, their energetics and various structural properties are recently studied (5). We propose the use of twobody basis function to describe various properties of bosonic clusters. With more than three particles the system becomes more complex as the number of degrees of freedom increases. We have investigated correlations between energies of the and systems and observe the generalized Tjon line (5) for large cluster. Now we consider the spectral statistics and spectral correlation of the atomic clusters of different sizes as these contain rich physics and also plays an pivotal role to establish the universal properties of quantum systems. Berry and Tabor conjectured that the fluctuation property of energy levels of a quantum system whose classical analog is regular, is characterised by Poisson statistics (6). Whereas, the fluctuation property of energy levels of a quantum system whose corresponding classical dynamical system is fully chaotic obeys the BohigasGiannoniSchmit (BGS) conjecture (7). This tells that Gaussian orthogonal ensemble (GOE) or Gaussian unitary ensemble (GUE) or Gaussian sympletic ensemble (GSE) statistics of random matrix theory, depending on time reversal symmetry and rotational symmetry of the system, will describe the fluctuation properties. However this conjecture is often interpreted in another way and the observation of level repulsion in the spectrum is treated as an indication of the nonintegrability of the system. The Poisson distribution implies complete randomness in the relative positions of energy levels as they are completely uncorrelated. On the other hand Wigner distribution implies strong correlation among the energy levels.
Earlier the spectral properties of many different quantum systems like atoms, atomic nuclei, quantum billiards have been studied (8); (9); (10); (11); (12); (13); (14); (15); (16). Also some attempts have been made for noninteracting manybosons and interacting bosonic system (17); (18); (19); (20). Recently we have reported the level spacing distribution of ultracold interacting bosons trapped in a harmonic potential (21); (22); (23). We found intriguing effect of both the interatomic interaction and the trap and observed deviation from the BGS cojecture. In this paper we are interested in similar type of calculation in the van der Waals bosonic clusters. Unlike the BoseEinstein condensate where the external trapping provides the stability of the condensate, the van der Waals clusters are bound due to the van der Waals interaction. In the very dilute condition one may treat it as a uniform Bose gas. Apart from the experimental interest, this kind of systems are also challenging for the following reasons. First, solving the manybody Schrödinger equation itself is a challenging numerical task due to many degrees of freedom and the obvious question is what kind of approximation is to be valid for the description of such clusters. Secondly for large cluster size when the system becomes very much correlated, one may expect Wigner type spectral distribution. However it needs an exhaustive study as level repulsion in the energy spectrum may not always lead to Wigner distribution which signifies chaos. It indicates that one may need to use some deformed GOE type of distribution for the correct description of nonintegrable but nonchaotic system. We propose to study several measures of spectral fluctuations and spectral correlation to determine the degree of influence of the interatomic interaction. This kind of study is also relevant as the statistical fluctuation can be directly observed experimentally in the context of ultracold Bose gases. We calculate nearest neighbour level spacing distribution (NNSD) , the level number variance and the DysonMehta statistics (24) for various cluster sizes. However all these measures require unfolding of the spectrum to remove variation in the density of energy levels in different parts of the spectrum. We can either unfold the spectrum of each member of the ensemble separately and form ensemble averaged NNSD or a single unfolding function can be used for all the members of the ensemble. Depending on the unfolding procedure, the final outcome of NNSD may vary. Moreover suitable unfolding function is not always known a priori and generally is approximated by higher order polynomials. Therefore to verify the outcome of the NNSD, we further analyze the distribution of quotients of successive spacings which does not require any unfolding and is independent of the energy level density.
The paper is organised as follows. In Section II, we introduce the manybody potential harmonic expansion method. Section III discusses the numerical results and Section IV concludes with the summary of our work.
Ii Methodology:Manybody calculation with potential harmonic basis
To study the spectral statistics and different spectral correlations we need to calculate a large number of energy levels of the diffuse Rb cluster. We approximately solve the full manybody Schrödinger equation by our recently developed Potential harmonic expansion method. We have earlier applied it successfully to study different properties of BEC (25); (26); (27); (28); (29); (30); (31) and atomic clusters (5); (32); (33). The methodology has already been described in detail in our earlier works (34); (35); (36). Hence here we describe it briefly for interested readers.
We consider a system of Rb atoms, each of mass and interacting via twobody potential. The timeindependent quantum manybody Schrödinger equation is given by
(1) 
Where is the total energy of the system, is the twobody potential and is the position vector of the th particle. It is usual practice to decompose the motion of a manybody system into the motion of the center of mass and the relative motion of the particles in center of mass frame. In absence of any confinig potential the center of mass behaves as a free particle in laboratory frame and we set its energy as zero. Hence, after elimination of the center of mass motion and using standard Jacobi coordinates, defined as (37); (38); (39)
(2) 
we obtain the equation for the relative motion of the atoms
(3) 
is the sum of all pairwise interactions. Now it is to be noted that Hyperspherical harmonic expansion method is an abinitio tool to solve the manybody Schrödinger equation where the total wave function is expanded in the complete set of hyperspherical basis (37). Although Hyperspherical harmonic expansion method is a complete manybody approach and includes all possible correlations, it is highly restricted to only. But for a diffuse cluster like Rbcluster, only twobody correlation and pairwise interaction are important. Therefore we can decompose the total wave function into twobody Faddeev component for the interacting pair as
(4) 
It is important to note that is a function of twobody separation () only and the global hyperradius , which is defined as . Thus the effect of twobody correlation comes through the twobody interaction in the expansion basis. is symmetric under the exchange operator for bosonic atoms and satisfy the Faddeev equation
(5) 
where is the total kinetic energy operator. In this approach, we assume that when () pair interacts, the rest of the bosons are inert spectators. Thus the total hyperangular momentum quantum number as also the orbital angular momentum of the whole system is contributed by the interacting pair only. Next the th Faddeev component is expanded in the set of potential harmonics (PH) (which is a subset of hyperspherical harmonic (HH) basis and sufficient for the expansion of ) appropriate for the () partition as
(6) 
denotes the full set of hyperangles in the dimensional space corresponding to the interacting pair and is called the PH. It has an analytic expression:
(7) 
is the HH of order zero in the dimensional space spanned by Jacobi vectors; is the hyperangle between the th Jacobi vector and the hyperradius and is given by = . For the remaining noninteracting bosons we define hyperradius as
(8)  
such that . The set of quantum numbers of HH is now reduced to only as for the noninteracting pair
(9)  
(10)  
(11) 
and for the interacting pair , and . Thus the dimensional Schrödinger equation reduces effectively to a four dimensional equation with the relevant set of quantum numbers: Energy , orbital angular momentum quantum number , azimuthal quantum number and grand orbital quantum number for any . Substituting in Eq(4) and projecting on a particular PH, a set of coupled differential equation for the partial wave is obtained
(12) 
where , ,
and .
is a constant and represents the overlap of the PH for
interacting partition with the sum of PHs corresponding to all
partitions (39).
The potential matrix element is given by
(13) 
Here we would like to point out that we did not require the additional shortrange correlation function for Rb clusters as was necessary for dilute BEC. A BEC is designed to be very dilute and hence confined by a harmonic oscillator potential of low frequency ( Hz). The average interatomic separation is thus very large () compared with the range of atomatom interaction (). Moreover the kinetic energy of the atoms is extremely small. Hence the effective interaction for large is controlled by the wave scattering length () (40). This is achieved by the inclusion of the correlation function (35); (36). On the other hand, diffuse van der Waals clusters are weakly bound by the actual interatomic van der Waals potential (of range ), without any confinement. Hence no correlation function is needed. The average interparticle separation is large enough, so that only twobody correlations are expected to be adequate, at least for light clusters.
Iii Results
iii.1 Choice of interaction and calculation of many body effective potential
As pointed earlier we choose the van der Waals potential with a hard core of radius as the interaction potential, = for and = for . For Rb atoms, the value of is 2803 eV (40). The unmanipulated scattering length corresponding to Rbdimer is . We obtain by solving the twobody Schrödinger equation for zeroenergy (36). We adjust the hard core radius in the twobody equation to obtain the dimer scattering length. In the Fig. 1 of Ref. (36) , we see the value of changes from negative to positive passing through an infinite discontinuity as decreases. Each discontinuity corresponds to one extra twobody bound state. We observe that tiny change in across the infinite discontinuity causes to jump from very large positive value to very large negative value. For our present calculation, we tune such that it corresponds to single bound state of the dimer. Thus calculated is for dimer scattering length of Rb atoms. With this set of values of and , we next solve the coupled differential equation [12] by hyperspherical adiabatic approximation (41). In hyperspherical adiabatic approximation, the hyperradial motion is assumed slow compared to hyperangular motion. For the hyperangular motion for a fixed value of , we diagonalize the potential matrix together with the hypercentrifugal term. Thus the effective potential for the hyperradial motion is obtained as a parametric function of . For the ground state of the system we choose the lowest eigenpotential [corresponding eigen column vector being ] as the effective potential. We plot the effective potential as a function of hyperradius , at the dimer scattering length and for various cluster size 3, 5 and 40 in Fig. 1.

(a) =3 

(b) =5 

(c) = 40 
With increase in cluster size the depth of the eigen potential increases sharply which indicates stronger binding of the cluster. The average size of the cluster also increases with increases in . The energy of the cluster is finally obtained by solving the adiabatically separated hyperradial equation in the extreme adiabatic approximation
(14) 
subject to appropriate boundary condition.
In our earlier published work we have reported ground state and few lowlying excitation of Rb cluster with maximum size of . (5). However in the calculation of level statistics and spectral correlation we also need higher multipolar excitations. In our manybody picture the collective motion of the cluster is described by the effective potential. The excited states in this potential are denoted by which corresponds to th radial excitation with th surface mode. Thus corresponds to the ground state and are the different excitations for . To calculate the higher levels with we follow the next procedure. We have noted that for , a large inaccuracy is involved in the calculation of the offdiagonal potential matrix. As the main contribution to the potential matrix comes from the diagonal hypercentrifugal term we disregard the contribution coming from offdiagonal part. Thus we get the effective potential for . Substituting in Eq. (14) we solve for different radial modes and repeat the numerical procedure for various to obtain the higher multipolar excitations.
Before discussing the statistical behavior of the energy spectrum we should discuss how accurate our calculated energy levels are. It is to be noted that the potential harmonic expansion method has been successfully applied in the calculation of collective excitations and thermodynamic properties of trapped bosons (42). For the investigation of thermodynamic properties we need to calculate a large number of energy levels. The calculated critical temperature and the condensate fraction are in good agreement with the experimental results (42). The effect of twobody correlations on thermodynamic properties of trapped bosons is also observed (42). Very recently we have also studied the energetics of diffuse Rb clusters (5) and also compared with the well studied He, Ne, and Ar clusters. Thus the calculated energy levels are accurate for further analysis. We also check for the convergence such that the error is considerably smaller than the mean level spacings.
iii.2 Levelspacing statistics for different cluster sizes
NNSD or distribution is the most common observable which is used to study the short range fluctuation. Now to compare the statistical property of different parts of the spectrum we need to unfold them. By unfolding, the smooth part of the level density is removed, it basically maps the energy levels to another with the mean level density equal to . For our present calculation we use polynomial unfolding of sixth order. We observe that for small cluster size with and , as the effective potential is very shallow, the number of energy levels are very small and not sufficient for the calculation of NNSD. Instead, we also calculate the manybody collective levels including higher order excitations with different . We then unfold each spectrum separately for a specific value of and then form an ensemble having the same symmetry. From the unfolded spectrum we calculate the nearest neighbour spacing as and calculate . is defined as the probablity density of finding a distance between two adjacent levels. Uncorrelated spectra obey the Poisson statistics which gives exponential distribution . Whereas for system with timereversal symmetry, level repulsion leads to the WignerDyson distribution (43).
The distribution of the unfolded spectrum with cluster size is plotted in Fig. 2. We observe that for very small and also for large . In our earlier calculation of Rb diffuse cluster, we have calculated the several lowenergy excitations. We have observed that due to the heavier mass of Rb atom, kinetic energy of Rb clusters is small while the interaction energy is large. It implies that although the system is tightly bound, it is less correlated for smaller . Thus unlike the trapped bosons, the smaller diffuse cluster does not exhibit any degeneracy in the calculation of lowlying excitations. It is reflected in Fig. 2(a) where we observe that for very small .


(a) lowest 22 levels 

(b) 300 levels 400 
The level spacing distribution for higher levels is shown in Fig. 2(b) which indicates that for such small cluster, the energy levels are completely uncorrelated. Though it looks very similar to the Poisson distribution the peak value at is less than . To determine how closely the histogram matches with the Poisson distribution we fit it with Brody distribution (14)
(15) 
where and is the Brody parameter. Depending on the value of the Brody parameter , this distribution interpolates between the Poisson distribution and the Wigner distribution (). Here we found . This implies that there is negligible correlation between the energy levels and the system is very close to regular. Actually for there are only interacting pair and the net attractive interaction is very weak.
Next to study the effect of interatomic interaction we gradually increase the effective interaction. We can vary the effective interaction either by tuning the scattering length or by changing the number of bosons. Here we increase the number of bosons to . It is already known from the earlier study of He cluster that decreases smoothly as a function of which indicates the saturation in the density and predicts liquiddrop behavior in He cluster with larger (5). However diffuse Rb cluster which is the system of our present interest is dilute and less compact which indicates sharp change in with change in cluster size. The average size of the cluster also increases. Thus the cluster with is more tightly bound, stable and more correlated compared with the cluster size . Due to more correlation in the energy spectrum, we can expect the very closely spaced energy levels which leads to the quasidegeneracy. This is reflected in Fig. 3 where we plot the distribution for the lowest 30 levels. The sharp peak in the first bin near clearly exhibits the signarure of quasidegeneracy. This peak is known as Shnirelman peak (44).

(a) 

(b) 
For better understanding of the structure of Shnirelman peak, we plot the same histogram in Fig. 4 as in Fig. 3 in finer details. Reducing the bin size gradually, a huge peak appears in the first bin which demonstrates the existence of global quasidegeneracy. The peak has a finite width which is further associated with Poisson tail. The resolution of the peak is further studied as the integral level spacing distribution (here being the number of levels), normalized to unity. We plot as a function of in Fig. 5. The linear dependence between and is shown in the left most part of Fig. 5 which represents the structure of the Shnirelman peak. Whereas the rightmost steep increase of corresponds to the Poisson tail.

(a) level 

(b) level 
However for the higher levels we observe that the system exhibits pseudointegrability. It is reflected in Fig. 6(a) where we observe the semiPoisson (SP) distribution. For comparison, in the same figure we plot the analytic expression of SP statistics given by (45). We observe the level repulsion at smaller values of (), where and asymptotic decay of is exponential. The SP distribution is observed within a narrow intermediate region between the quasidegenerate regime and the completely integrable regime. distribution for the higher levels are plotted in the Fig. 6(b) which is again very similar to Poisson distribution. We again fit the Brody distribution with the histogram and find the Brody parameter . The observation of SP distribution and increase in the value of Brody parameter clearly manifests the enhanced effect of interatomic correlation with increase in cluster size. However we fail to give any physical reason which causes this SP and Poisson statistics. As pointed earlier, for smaller cluster size only effective potential is not enough to calculate sufficient number of levels for the study of distribution. So the findings of SP statistics may be physically acceptable whose origin is not clear to us or it may be due to overlap of several values.
Thus to get further insight we significantly change the cluster size to where only effective potential is deep enough to support sufficient number of states for calculation of distribution. We plot the distribution in Fig. 7. We observe similarity with Wigner distribution as very small value of near signifies the level repulsion. However the peak at overshoots . The large peak at signifies large accumulation of levels with level spacing . Though we tried to fit the histogram again with Brody distribution we fail to appropriately fit it. Now it is worthy to mention that Guhr and Weidenmüller (46) proposed a modified uniform spectrum in terms of a deformed GOE, which combines uniform, GOE and Poisson. As the distribution of Fig. 7 is quite similar to Fig. 1, Fig. 2, and Fig. 3 of Ref. (46), the use of deformed GOE may be an ideal step for future investigation. As the Fig. 7 does not match with the Wigner distribution we conclude that the Hamiltonian is not chaotic. However the deformed GOE type distribution signifies the system is strictly nonintegrable and exhibits strong interatomic correlation. Thus it is indeed required to calculate the energy level correlation which we discuss in the following section.
iii.3 Energy level correlation
So far we have considered only the NNSD which is commonly used to characterize the shortrange fluctuations in the spectrum. However in order to confirm our findings of the effect of correlation on the spectral properties and to investigate how the correlation gradually builds in with the increase in cluster size which makes the system too complex, we study the long range correlations of the spectrum. The level number variance is the most commonly used observable to characterize correlations between pair of levels. It mainly determines the longrange fluctuations in the spectrum. It is defined as the average variance of the number of levels in the energy interval containing an average number of levels and is calculated as
(16) 
where represents the average over the energy value and determines the number of eigen energy levels below . For the uncorrelated Poisson statistics , whereas for GOE, increases logarithmically with . From the earlier study of level spacing distribution it has been observed that for the system exhibits features which are very close to the noninteracting limit. We have also observed the Poisson distribution in the level statistics of higher levels. However the most interesting observation is the semiPoisson distribution for the intermediate part of the spectrum for . The corresponding is plotted in Fig. 8(a). It approximately increases linearly as which is the value of number variance of SP distribution. Then we plot for higher part of the spectrum in Fig. 8(b). It is approximately proportinal to indicating that the system is correlated but does not exhibit any level repulsion. This further confirms the findings of the Poisson distribution in the distribution. For strongly correlated cluster with we observe that approximately increases logarithmically with [Fig. 9]. This feature is close to GOE results. However there are significant differences between our numerical results and the Wigner surmise. It again indicates that the system does not show full chaos though it exhibits strong nonintegrability.

(a) , 

(b) , 
The other important observable to characterise longrange correlation is statistics (24). Given an energy interval of length , it is defined as the least square deviation of the staircase function from the best straight line fitting it:
(17) 
It is customary to use the average values of . Thus statistics, averaged over energy intervals, measures the deviation of the unfolded spectrum from the equidistant spectrum and hence it gives information on the rigidity of spectrum or spectral stiffness. For uncorrelated Poisson spectra whereas for Wigner spectra . Our
calculated numerical results for is shown in Fig. 10. Though it looks similar to GOE distribution, it is significantly lower than the GOE results which confirms our earlier observation for large cluster. The approach of towards the GOE behavior is valid only upto . The similar kind of observation was made in Fig. 6 of Ref. (46) where deformed GOE behavior is noted in . It indicates that for large cluster size, the levels are strongly correlated. Whereas for smaller cluster (), we observe that distribution gradually approaches to Poisson as we move upward in the spectrum. For a small
intermediate region of the spectrum lies between the GOE and Poisson distribution [Fig. 11(a)] whereas for the upper levels it almost perfectly follows the Poisson distribution [Fig. 11(b)] which indicates that the spectrum has turned soft.

(a) 

(b) 
iii.4 Quotients of successive spacings
Before concluding the paper, we present in this Section, as a test of the observations made in Sections III B, the results of the analysis of the distribution of quotients of successive level spacings [denoted by ], a measure introduced recently, that is independent of the unfolding function and the unfolding procedure (47); (48). Note that in all the analysis presented in Sections III B and III C, we have employed a sixthorder polynomial for the density of levels for unfolding. The P(r) distribution and the related averages allow for a more transparent comparison with experimental results than the traditional level spacing distribution and this measure is particularly important for manybody systems as the theory for the eigenvalue (level) densities for these systems is usually not available. In the recent past, this measure was used in analyzing manybody localization (47); (49); (50); (51) and also in quantifying the distance from integrability on finite size lattices (52); (53). More recently, using it is established conclusively that embedded random matrix ensembles for manybody systems, generated by random interactions in the presence of a meanfield, follow GOE for strong enough twobody interaction (54).
Given a ordered set of the energy levels , the nearest neighbor spacing and the probability distribution of the ratios is subject to normalization . If the system is in integrable domain (described by Poisson NNSD), then the is given by
(18) 
and if the system is chaotic (described by GOE), then the is given by Wignerlike surmise (48),
(19) 
The average value of r, i.e. , is for GOE and is for Poisson. It is also possible to consider . The average value of , i.e. , is for GOE and for Poisson.
Some results for vs for the spectrum of diffuse Rb cluster with the same cluster sizes as above viz. = 3, 5 and 40, are shown in Figs. 12 and 13. Moreover, we have also calculated the averages and and results are given in Table 1. For with levels 122, there is a peak at as seen from Fig. 12a. Similarly for levels 300400, is close to Poisson form as shown in Fig. 13a. These results are consistent with the NNSD results in Figs. 2a and 2b respectively. In addition, the results for and given in Table 1 are also in agreement with these observations. Turning to , with levels 130 the shows peaks at and (see Fig. 12b) and for quantifying this structure, it is necessary to derive that corresponds to Shnirelman peak. Going to levels 4080, it is seen from Fig. 12c that exhibits level repulsion with for but the form of shows clear deviations from the GOE result given by Eq. (19). In order to compare with the conclusion drawn from NNSD in Fig. 6a, it is necessary to derive the formula for for pseudointegrable systems (these systems give semiPoisson form for NNSD). Turning to levels 8501000, it is clearly seen from Fig. 13b that the is close to Poisson and this is in complete agreement with NNSD shown in Fig. 6b. Further, for the curve shows level repulsion and it is closer to GOE than to Poisson (see Fig. 13c). Also, the values of and (shown in Table 1) are close to GOE results. Thus N=40 example exhibits level repulsion as seen in the NNSD result. Combining all these observations, we conclude that the results deduced from NNSD analysis are consistent with those obtained from analysis and thus the unfolding procedure used in Sections III B and III C can be considered to be good.

(a) =3 

(b) =5 

(c) = 40 

(a) =3 

(b) =5 

(c) = 40 
N=3  levels (122)  0.76  1.168 
levels (300400)  0.34  204151  
N=5  levels (130)  0.48  6.18 
levels (4080)  0.64  1.64  
levels (8501000)  0.39  144078  
N=40  levels(160200)  0.76  1.43 
GOE  0.5359  1.75  
Poisson  0.3863 
Iv Conclusions
Study of energy level statistics plays an important role in elucidating the universal properties of quantum systems. Berry and Tabor conjectured that the eigenenergy levels of a quantum system whose classical dynamics shows integrability, must exhibit the fluctuation property as determined by the uncorrelated Poisson statistics. This is in sharp contrast with the BGS conjecture which asserts that the fluctuation property of energy levels of a quantum system whose classical dynamics should exhibit GOE (or GUE or GSE) statistics. However, complicated quantum manybody systems often lie between these two contrasting conjectures.
Thus the purpose of present paper is to consider a relatively complex quantum system whose experimental realization is possible. The van der Waals bosonic cluster is such a quantum system which starts to be more and more complex with increase in cluster size. The above mentioned contrasting conjectures have been examined thoroughly by using various statistical observables like NNSD, level number variance and the spectral rigidity . These observables highlight the short and long range correlation, level repulsion, level clustering and how the features of the above observables crucially depend on the cluster size are also focussed. Our detailed numerical analysis reveals that for smaller cluster when the system is very close to integrability, Berry and Tabor conjecture is followed. For large cluster although we observe similar to BGS conjecture, however deviation occurs. For large clusters the system becomes strongly correlated but does not exhibit true chaos. However our present study reveals that the deformed GOE type of distribution may be suitable for future investigation.
Acknowledgements
This work is supported by the Department of Atomic Energy (DAE), government of India, through Grant No. 2009/37/23/BRNS/1903. SKH acknowledges the Council of Scientific and Industrial Reaserch (CSIR), India for a senior research fellowship through NET (Grant No: 08/561(0001)/2010EMR1). NDC acknowledges financial support from the University Grants Commission (UGC), India [Grant No: F.40425/2011 (SR)]. SKH also acknowledges hospitality of the Maharaja Sayajirao University of Baroda, Vadodara, India during a recent visit for this work.
Footnotes
 Present address: Department of Physics, Presidency University, 86/1 College Street, Kolkata 700073
References
 T. Kramer et. al., Nature (London) 440, 315 (2006).
 E. A. Donley et. al. Nature (London) 417, 529 (2002).
 C. H. Chin, T. Kraemer, M. Mark, J. Herbig, P. Waldburger, H. C. Nagerl, R. Grimm Phys. Rev. Lett. 94, 123201 (2005).
 M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995)
 P. K. Debnath, B. Chakrabarti, T. K. Das and S. Canuto, J. Chem. Phys. 137, 014301 (2012).
 M. V. Berry and M. Tabor, Proc. R. Soc. London, A 356, 375 (1977).
 O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev. Lett. 52, 1 (1984).
 G. Casati, F. ValzGris, and I. Guarneri, Lett. Nuovo Cimento Soc. Ital. Fis. 28, 279 (1980).
 F. Haake, Quantum Signatures of Chaos, (Springer, New York, 2010)
 T. H.Seligman, J. J. M. Verbaarschot, and M. R. Zirnbauer, Phys. Rev. Lett. 53, 215 (1984).
 T. Zimmermann, H.D. Meyer, H. Köppel, and L. S. Cederbaum, Phys. Rev. A 33, 4334 (1986).
 G. Tanner, K. Richter, and J.M. Rost, Rev. Mod. Phys. 72, 497 (2000).
 J. Sakhr and N. D. Whelan, Phys. Rev. A 62, 042109 (2000).
 T. A. Brody et al., Rev. Mod. Phys. 53, 385 (1981).
 V. K. B. Kota, Phys. Rep. 347, 223 (2001).
 J. M. G. Gómez, K. Kar, V. K. B. Kota, R. A. molina, A. Relaño and J. retamosa, Phys. Rep. 499, 103 (2011)
 L. Muñoz, E. Faleiro, R. A. Molina, A. Relaño, and J. Retamosa, Phys. Rev. E 73, 036202 (2006).
 R. J. Leclair, R. U. Haq, V. K. B. Kota and N. D. Chavda, Phys. Lett. A 372, 4373 (2008).
 N. D. Chavda, V. Potbhare, and V. K. B. Kota, Phys. Lett. A311, 331 (2003).
 M. Vyas, V. K. B. Kota, N. D. chavda and V. Potbhare, J. Phys. A, Math. theor. 45, 265203 (2013); arXiv:1010.6054.
 B. Chakrabarti, A. Biswas, V. K. B. Kota, K. Roy and S. K. Haldar, Phys. Rev. A 86, 013637 (2012).
 K. Roy, B. Chakrabarti, A. Biswas, V. K. B. Kota and S. K. Haldar, Phys. Rev. E 85, 061119 (2012).
 K. Roy, B. Chakrabarti and V. K. B. Kota, Phys. Rev. E 87, 062101 (2013).
 M. L. Mehta, Random Matrices (Academic Press, New York, 1991).
 B. Chakrabarti, T. K. Das, and P. K. Debnath, Phys. Rev. A 79, 053629 (2009).
 S. K. Haldar, B. Chakrabarti, and T. K. Das, Phys. Rev. A 82, 043616 (2010).
 P. K. Debnath and B. Chakrabarti, Phys. Rev. A 82, 043614 (2010).
 A. Biswas, T. K. Das, L. Salasnich, and B. Chakrabarti, Phys. Rev. A 82, 043607 (2010).
 A. Biswas, B. Chakrabarti, T. K. Das, and L. Salasnich, Phys. Rev. A 84, 043631 (2011).
 S. K. Haldar, B. Chakrabarti, T. K. Das and A. Biswas, Phys. Rev. A 88, 033602 (2013).
 S. K. Haldar, P. K. Debnath and B. Chakrabarti Eur. Phys. J. D 67, 188 (2013).
 T. K. Das, B. Chakrabarti and S. Canuto J. Chem. Phys. 134, 164106 (2011).
 S. K. Haldar, B. Chakrabarti, and T. K. Das FewBody Systems 53, 283 (2012).
 T. K. Das and B. Chakrabarti, Phys. Rev. A 70, 063601 (2004).
 T. K. Das, S. Canuto, A. Kundu, B. Chakrabarti, Phys. Rev. A 75, 042705 (2007).
 T. K. Das, A. Kundu, S. Canuto, and B. Chakrabarti, Phys. Lett. A 373, 258261 (2009).
 J. L. Ballot and M. Fabre de la Ripelle, Ann. Phys. (N. Y.) 127, 62 (1980).
 M. Fabre de la Ripelle, Ann. Phys. (N. Y.) 147, 281 (1983).
 M. Fabre de la Ripelle, FewBody System 1, 181 (1986).
 C. J. Pethick and H. Smith,BoseEinstein Condensation in Dilute Gases (Cambridge University Press, Cambridge, England, 2001).
 T. K. Das, H. T. Coelho and M. Fabre de la Ripelle, Phys. Rev. C 26, 2281 (1982).
 S. Bhattacharyya, T. K. Das and B. Chakrabarti, Phys. Rev. A 88, 053614 (2013).
 O. Bohigas and M. J. Giannoni, in Mathematical and Computational Methods in Nuclear Physics, Vol. 209 of Lecture Notes in Physics, edited by J. S. Dehesa, J. M. G. Gomez, and A. Polls (Springer, New York, 1984).
 A. I. Shnirelman. Usp. Mat. Nauk 30, 265 (1975); A. I. Shnirelman, addendum in V. F. Lazutkin. KAM Theory and Semiclassical Approximations to eigenfunctions (Springer, Berlins, 1993).
 A. M. GarcíaGarcía, Phys. Rev. E 72, 066210 (2005).
 T. A. Guhr and H. A. Weidenmüller, Ann. Phys. (N.Y.)193, 472 (1989).
 V. Oganesyan, D. A. Huse, Phys. Rev. B 75, 155111 (2007).
 Y. Y. Atas, E. Bogomolny, O. Giraud, G. Roux, Phys. Rev. Lett. 110, 084101 (2013).
 V. Oganesyan, A. Pal, D. A. Huse, Phys. Rev. B 80, 115104 (2009).
 A. Pal, D. A. Huse, Phys. Rev. B 82, 174411 (2010).
 S. Iyer, V. Oganesyan, G. Refael and D. A. Huse, Phys. Rev. B 87, 134202 (2013).
 C. Kollath, G. Roux, G. Biroli, A. M. Läuchli, J. Stat. Mech. P08011 (2010).
 M. Collura, H. Aufderheide, G. Roux, D. Karevski, Phys. Rev. A 86, 013615 (2012).
 N. D. Chavda, V. K. B. Kota, Phy. Lett. A 377, 3009 (2013).