Ground-state properties of a dilute homogeneous Bose gas of hard disks in two dimensions

Departament d’Electrónica, Enginyeria i Arquitectura La Salle, Pg. Bonanova 8, Universitat Ramon Llull, E-08022 Barcelona, Spain Department d’Estructura i Constituents de la Matèria, Diagonal 645, Universitat de Barcelona, E-08028 Barcelona, Spain Dipartimento di Fisica “E. Fermi,” Università di Pisa, and INFN, Sezione di Pisa, Via Buonarroti, 2 I-56100 Pisa, Italy sReceived 20 October 2004; published 10 March 2005 d


I. INTRODUCTION
The study of two-dimensional ͑2D͒ quantum many-body systems is a subject of considerable interest both from the theoretical and the experimental points of view.Atomic helium monolayers have attracted much attention in these respects ͓1͔, and recently Bose-Einstein condensation ͑BEC͒ in a strongly asymmetrical, quasi-two-dimensional trap has been achieved ͓2,3͔.In a homogeneous 2D gas the one-body density matrix acquires a power law decay below some low critical temperature, T Ͻ T c ͓4͔ and, at T = 0, its r → ϱ limit differs from zero.This behavior defines a true condensate and BEC can occur in such a limit.The modified density of states in confined geometry, as it is the case for atoms in harmonic traps, makes BEC appear even at finite temperature in 2D ͓5,6͔.
In a recent paper ͓7͔, referred to as I hereafter, we have examined the structure of the ground state of a 3D homogeneous gas of bosons, interacting through both hard and soft core potentials.The study was carried out to enlighten the role of the interaction induced correlations along the density of the system and the possible occurrence of universality in the dependence of different properties on the gas parameter, x = na 3 ͑n being the particle density and a the s-wave scattering length͒.The calculations were performed within the correlated basis functions ͑CBF͒ theory, using optimized Jastrow correlated wave functions ͓8͔.The results were in excellent agreement with the existing diffusion Monte Carlo ͑DMC͒ ͓9͔ ones, obtained by the exact stochastic solution of the many-body Schrödinger equation.Moreover, the reliability ͑and limitations͒ of the analytical expansions in x of the energy per particle and of the condensate fraction were clearly assessed ͓10͔.
In this work we use the same variational approach to study the T = 0 ground state of a homogeneous gas of N → ϱ hard disks ͑HD͒ of radius a.We pay particular attention to the long-range structure of the wave function, which is in turn intimately related to the long-range order of the density matrix and to the condensate fraction.
At present, there are not published DMC results for these systems in the literature.However, based on the 3D case experience we believe that the variational approach may provide a very accurate description of the 2D ground state properties and serve as a test of validity for the available low density expansions ͑LDE͒.This fact may be of relevance also to the study of BEC in 2D harmonic traps, since LDE are related to the Gross-Pitaevskii equation and its modification ͓11͔.LDE for the energy per particle in terms of the 2D gas parameter, x = na 2 , have been derived by several authors, starting from the leading order ͓12͔: where the energy is in units of ប 2 /2ma 2 .An interesting question, that we do not address in this paper, is the insurgence and the degree of universality in 2D.In the 3D Bose gas a dependence in the energy on the shape of the potential appears already at x ϳ 0.001 ͓7,9,13͔, and breaks down earlier for quantities other than the energy, as the short range structure of the distribution function and the condensate fraction.A similar study for the 2D gas would be of great interest by itself and also in view of the analysis of the quasi-2D BEC experiments in terms of the Gross-Pitaevskii equation.
In I a Jastrow correlation operator was considered, where the two-body correlation function, f͑r ij ͒, depends on the interparticle distance, r ij = ͦr i − r j ͦ, vanishes for r ij ഛ a and goes to unity at large distances.Expectation values of operators on the Jastrow correlated wave function correspond to multidimensional integrals ͑we consider only operators dependent on the spatial coordinates͒ which can be either directly evaluated by Monte Carlo type integrations ͑variational Monte Carlo, VMC͒ or expanded in cluster diagrams.The hypernetted chain ͑HNC͒ set of integral equations allows for summing an infinite number of diagrams, but it can be only approximately solved, since the class of the elementary diagrams is not summable in a closed way ͓14͔.The approximation amounting to disregard the elementary diagrams ͑HNC/0͒ is expected to be reliable in the low density regime, as shown in I.
In this paper we consider a system of N spinless bosons of mass m on a surface S, described by the Hamiltonian where V͑r͒ is a two-body, symmetric hard-disk potential, and the radius a corresponds, for this potential, to the 2D scattering length.We consider the system in the thermodynamic limit ͑N and S → ϱ, keeping the density, n = N / S, con-stant͒.
Minimization of the energy with respect to f͑r͒ provides the optimal Jastrow factor, which can be obtained through the solution of the Euler-Lagrange ͑EL͒ equation ͓15,16͔, ␦E͓f͔ / ␦f = 0. We adopt here the HNC/0 energy functional to solve the EL equation.The resulting correlation, f EL ͑r͒, shows a very long range structure ͓f EL ͑r → ϱ͒ → 1−␣ / r, with ␣ constant͔ not easily accessible to VMC, due to the limited size of the simulation box.As a consequence we have also used a parametrized shorter range correlation factor, suitable to be used in VMC and whose parameters are variationally fixed.The comparison between the HNC/0 and VMC calculations with this short range correlation provides a check of the accuracy of the truncation in the cluster expansion.Diffusion Monte Carlo is also likely to suffer from the long range behavior of the wave function, as these fine structures, related to collective effects, are hardly distinguishable by numerical, finite size simulations.
The plan of the paper is as follows: the correlated variational theory and its implementation ͑HNC and Euler equa-tions͒ are briefly outlined in Sec.II; an analysis of the EL asymptotic behaviors is presented in Sec.III; Sec.IV presents the results for the HD model; while summary and conclusions are given in Sec.V.

II. VARIATIONAL THEORY
Given the Jastrow correlated wave function ͑3͒, the energy per particle of the 2D homogeneous Bose gas is This expression simplifies for the HD potential ͑5͒, since the radial distribution function ͑RDF͒, vanishes inside the core of the potential and only the kinetic part contributes to the energy.
The HNC equations provide a procedure to evaluate g͑r͒, where N͑r͒ and E͑r͒ are functions representing the sum of the nodal and the elementary diagrams, respectively ͓14͔.
The function E͑r͒ is an input to the HNC equations, which can be solved once a choice for it has been done.The simplest possible approximation corresponds to set E͑r͒ =0 ͑HNC/0͒.This apparently drastic truncation is, however, justified at low densities since the elementary diagrams, due to their high connectivity, do not appreciably contribute at the low densities relevant to BEC experiments.Otherwise, the energy and the RDF can be stochastically evaluated by Monte Carlo sampling of the corresponding many-body integrals ͑7͒.The explicit calculation can be performed by using the standard Metropolis algorithm ͓17͔ ͑see Ref. ͓18͔ for a detailed description of Monte Carlo methods͒.
Connected to the RDF is the static structure function ͑SSF͒, S͑k͒,

͑9͒
often used in the analysis of ͑and extractable, within some approximations, from͒ scattering experiments in condensed matter physics ͑e.g., neutron scattering off liquid helium ͓19͔͒.
The Bose condensate is linked to the nonzero, long-range order of the one-body density matrix ͑OBDM͒, In fact, n 1 ͑r 11 Ј → ϱ͒ / n = n 0 , where n 0 is the condensate fraction.The depletion of n 0 with respect to unity is an unmistakable indication of interparticle interactions ͑and, as a consequence, of correlations͒.The Fourier transform of the OBDM provides the momentum distribution ͑MD͒, As for the distribution function, the correlated OBDM can be computed by using HNC theory ͓20͔.In fact, n 1 ͑r͒ can be expressed in terms of new nodal and elementary functions as where N ww ͑r͒ is the solution of a generalized HNC equation.Again, setting E ww ͑r͒ = 0 these equations can be solved in the HNC/0 approximation.An appropriate choice of the correlation factor is essential for the effectiveness of the variational approach.As stated in the Introduction, the best choice is the one satisfying the Euler-Lagrange equation, ␦E͓g͔ ␦g͑r͒ = 0, ͑13͒ which has been written in terms of the RDF rather than f͑r͒, since in the HNC/0 scheme there is a one-to-one correspondence between these two quantities.The correlation function can then be obtained by inversion of the HNC equations.The EL equations are solved both in configuration and momentum space, as discussed at length in I for the 3D hard-sphere case.The theory can be straightforwardly applied to the 2D hard-disks gas and it is briefly outlined in the Appendix.The correlation function produced by the solution of the EL equation shows a long-range structure that is discussed in the next section.Finite size Monte Carlo techniques have difficulties to correctly deal with this long range behavior and correlations healing to unity inside the simulation box are used.In this respect, we have also adopted a parametrized shorter-range correlation, f SR ͑r͒, obtained by minimizing the energy computed at the two-body order of the cluster expansion, g 2B ͑r͒ = f 2 ͑r͒, constrained by a normalization condition

͑14͒
which requires the use of a Lagrange multiplier .The minimization is performed under the healing condition f SR ͑r ജ d͒ = 1, while is fixed by imposing fЈ͑d͒ = fЉ͑d͒ = 0.The healing distance d is taken as a variational parameter to minimize the many-body energy.
For the particular case of the hard-core potential, the short range correlation, f SR ͑r͒, vanishes at r ഛ a, while where =−4m / ប 2 and J 0 and Y 0 are Bessel functions of the first and second kind, respectively.In the d → ϱ limit, f SR coincides with the exact solution of the 2D zero energy Schrödinger equation, f͑r͒ ϰ ln͑r͒ ͓21͔.

III. ASYMPTOTIC BEHAVIORS
Applying to the 2D gas the sum-rule analysis of Ref. ͓8͔, it can be shown that the SSF at low momenta has a linear dependence, similar to the 3D gas: where c is the sound velocity in the medium.In terms of the RDF and of the correlation function, it corresponds to a long range behavior, different, however, from the 3D case where ͑g 3D −1͒ ϰ 1/r −4 and ͑f 3D −1͒ ϰ 1/r −2 at large r values.
In fact, the 2D Fourier transform for a function, h͑r͒, with circular symmetry can be written as while h͑r͒ is obtained from h ˜͑k͒ as Assuming that H ˜͑k͒ϵkh ˜͑k͒ is a well behaved function at the origin, one finds where only the k = 0 values of the function and of its even derivatives enter in the expansion.
The RDF is obtained from S͑k͒ through so one readily obtains, from Eqs. ͑16͒ and ͑19͒, the asymptotic limit By inverting the HNC/0 equations ͑8͒ one finds for the nodal and correlation functions the following limits: showing that in 2D, correlations have longer range than in 3D.
The long range structure of the OBDM is derived from the previous expressions and the HNC equations.In fact, given the structure ͑23͒ we obtain for N ͑r͒ and for the OBDM The momentum distribution has the same longwavelength limit shown in 3D ͓22͔, namely

IV. RESULTS
In this section we present and analyze results for the energy, radial distribution function, static structure function and one-body density matrix of the hard-disks gas.We have used the optimized Jastrow wave function obtained from the solution of the Euler-Lagrange equations and the short-range correlation of Eq. ͑15͒, mainly to establish a comparison between the HNC/0 and the VMC approaches.In the following, dimensionless quantities will be used: energies and distances will be given in units of ប 2 /2ma 2 and a, respectively.
Several corrections to E LO have been proposed in the literature.Kolomeisky and Straley ͓23͔ used renormalization group techniques to study the ground state of dilute Bose systems as a function of the space dimensionality.They found a general expression, valid for strong interactions when x → 0, that simplifies for hard disks to Cherny and Shanenko ͓24͔ derived an alternative expansion, in the parameter u satisfying the equations: where ␥ = 0.577. . . is the Euler constant.These authors also gave an expansion of u in terms of ␦, allowing us to write the series ͑28͒ in the form

͑30͒
Table I reports the energy per particle as a function of x in the EL approximation ͑EL/HNC͒, the variational Monte Carlo calculation starting from f SR ͑r͒ ͑SR/VMC͒ and the HNC approach with the same correlation function ͑SR/ HNC͒.We remind that the HNC/0 approximation is used everywhere but in the VMC.The results of the x expansions previously discussed are also reported.
The comparison between SR/VMC and SR/HNC shows that the influence of the missing elementary diagrams on the energy is less than 1%, except at the highest value x = 0.1 ͑ϳ2.3% ͒.This gives us confidence that the variational principle is mostly satisfied within our HNC/0 calculations, providing the hierarchy: E exact ഛ E EL ഛ E SR .Only at x = 0.1 these inequalities do not numerically hold.The cases x =10 −5 and x =5ϫ 10 −2 can be considered to fulfill the inequalities if we take into account the numerical accuracy associated to the calculation at these quantities.If we estimate the contribution of the elementary diagrams on E EL by scaling it by the ratio E SR VMC / E SR HNC , we obtain E EL / N͑x = 0.1͒ = 0.9075, restoring all the inequalities.
Lieb ͓25͔ pointed out that a lower bound to the exact energy is given by …͔, ͑31͒ and that E exact / E LO → 1 when x → 0. Both the variational energies ͑EL and SR͒ comply with condition ͑31͒ and seem to tend to E exact when x goes to zero.Table I gives also the healing distance, d, of the SR correlation in units of a. d increases when x → 0, the Lagrange multiplier decreases and the energy goes to zero.Therefore, f SR ͑r͒ x→0 can be approximated by its = 0 limit, which coincides with the zero energy limit of the two-body Scrödinger equation, These results are also shown in Fig. ͑1͒, where the variational scaled energies per particle ͑EL, SR/VMC, and SR/ HNC͒ are compared with the E KS and E CS estimates.All energies have been divided by E LO in order to stress the deviations from the low-density limit.The limit is approached by E EL and E SR from above when x decreases, although it has not been yet fully reached at x =10 −5 .The differences between the variational and the low density energies are still visible, even in the density range relevant to BEC experiments ͓3͔.
E KS does not satisfy the x = 0 lower bound ͑31͒ and, starting from x Ϸ 0.005, becomes higher than the variational upper bound provided by E EL and E SR .E CS satisfies the low density limit, but lies above the variational upper bounds at any value of x.E CS ͑u͒ is always larger than E CS ͑␦͒ and the difference increases drastically along x.Notice that Eq. ͑29͒ does not have solution at x ജ 0.0369.
The EL optimization procedure does not significantly affect the energies obtained with f SR ͑r͒, since the energy is dominated by the short range structure of the HD potential, requiring g͑r͒ to vanish inside the core.The effects of the long range structure of the EL correlation are, instead, clearly evident in the behavior of the radial distribution function.
The EL RDF is shown in Fig. ͑2͒ for different values of x.At low x, g EL ͑r͒ is a monotonically increasing function of the distance.However, it develops a local maximum close to the core radius at densities x Ͼ 0.01.This is a genuine manybody effect induced by the strong correlations at high density.The same behavior was found in I for the 3D Bose gas.As expected, the correlation hole is more pronounced at larger densities.
The long range limit of g EL ͑r͒ is shown in Fig. 3 at x = 0.01 and x = 0.1.The quantity shown is xr 3 (g͑r͒ −1) whose dimensionless asymptotic limit is This ratio is smaller at x = 0.1, implying that the sound velocity increases with x, as expected.Consistent with the previous figure, the asymptotic limit is reached faster at larger densities.Figure 4 gives the EL correlation and the nodal function, N͑r͒, at x = 0.001.We show ͓xr(1− f EL ͑r͒)͔ and ͓xrN͑r͔͒ to enlighten the asymptotic limits ͑23͒ and ͑22͒ whose dimensionless values are c /4 and c /2, respec-TABLE I. Energy per particle for the hard disk model, as a function of x.E EL is obtained by solving the HNC/0 EL equation; E VMC and E SR are the VMC and HNC/0 energies with the SR correlation having a healing distance d, in units of the scattering length.⑀ VMC is the statistical error in the variational Monte Carlo calculation.The last four columns give the energies in different low density expansions ͑see text͒.
The EL static structure function, S͑k͒, is shown in Fig. 5.At low densities the SSF reaches the asymptotic value, S͑k͒ → 1, already at k ϳ 1.As in the RDF case, the overshooting of the SSF at the highest density, x = 0.1, is a consequence of the correlations.The linear regime of S͑k͒ around the origin is appreciable, although the calculation of the ratio S͑k͒ / k shows deviations from a constant value already at low k.
In Fig. 6 we plot the one-body density matrix, n 1 ͑r͒, in the EL approach for x = 0.01, 0.005, and 0.001.The asymptotic limit of n 1 ͑r͒ defines the value of the condensate fraction, which decreases when the gas parameter increases.The asymptotic value is reached faster when x increases.The detailed long-range behavior ͑25͒ is presented in Fig. 7 by showing the quantity r͓n 1 ͑r͒ / ͑n 0 n͒ −1͔, whose dimensionless asymptotic value is c / ͑8x͒.Even if the speed of sound increases with x, the value of this limit is dominated by the presence of the gas parameter in the denominator and the overall quantity increases when x decreases.Finally, the EL and SR/VMC condensate fractions, n 0 ͑x͒, are reported in Fig. 8.The figure also contains the low-density prediction ͓12͔ n 0 LD appears to sensibly overestimate the condensate fraction.The EL and SR/VMC condensate fractions are very similar except for the largest value of x reported in the figure, where the contribution of the elementary diagrams could be important.This fact indicates that the value of n 0 is not very much affected by the inclusion of a long-range structure into the correlation function.However, the use of f EL ͑r͒ is crucial to approach this value in a proper way, that is, to satisfy Eq. ͑25͒.Also reported in the figure is the condensate fraction of the 3D system of hard spheres, taken from I. At fixed x, the 2D condensate fraction is smaller than the 3D one, indicating that correlations in the 2D system are stronger, thus promoting more particles outside the zero-momentum state.

V. SUMMARY AND CONCLUSIONS
In this work we have analyzed the energy and structure of a homogeneous gas of bosons in two dimensions interacting via a hard disk potential whose core radius equals its corresponding 2D scattering length.We have adopted a variational many-body approach, based on a Jastrow correlated ground state wave function.The expectation values have been computed both in the framework of the hypernetted chain theory ͑within the HNC/0 approximation͒ and with the variational Monte Carlo method.Two types of correlation functions have been used: ͑i͒ a long range one, obtained by the free minimization of the HNC/0 ground state energy and preserving the correct asymptotic behaviors of the wave function; ͑ii͒ a short range one, to be used in the Monte Carlo sampling and providing a check of the accuracy of the cluster expansion.
By comparing with the VMC results, the accuracy of the HNC/0 energies are better than 1%, except at the highest density, x = 0.1, where the error is still less than 3%.The EL minimization lowers the energy with respect to the SR correlation by ϳ1.5%.We do not expect further large reductions from a complete DMC calculation.The low density expansions start to severely deviate from the variational results already at x ϳ 0.001, and the most accurate of them appears to be the Cherny and Shanenko expansion in terms of the parameter ␦.However, their use for estimating corrections to the 2D Gross-Pitaevskii equation, especially in the large gas parameter regime, seems questionable.
Finally, the condensate fractions lies well below the values predicted by the low density theories and also below the results for the three dimensional gas of hard spheres at the same gas parameter.
We conclude that the variational theory is a powerful and reliable tool to study dilute systems, also in 2D.Moreover, the homogeneous gas results may be used in a local density type approximation ͓11͔ to analyze bosons in two dimensional harmonic traps.␦E͓g͔ ␦g͑r͒ = 0, ͑A1͒ can be obtained in momentum space, yielding where t͑k͒ = ប 2 k 2 /2m is the free-particle energy spectrum and V ph ͑k͒ is the particle-hole interaction.The latter can be expressed in configuration space and reads, disregarding the contribution of elementary diagrams,

͑A3͒
in terms of the induced interaction I ͑r͒.In momentum space, we have Equations ͑A2͒-͑A4͒ are to be solved simultaneously.This can be done starting from a suitable choice for g͑r͒, performing its FT to get S͑k͒, evaluating I ͑k͒ and V ph ͑r͒, and then deriving a new S͑k͒ with the help of Eq. ͑A2͒.This procedure is iterated until the difference between two consecutive iterations is as small as desired.
Up to this point there are no formal differences between the 2D and 3D cases.The main deviation is the way in which the Fourier transforms are carried out.In 2D and for a general function, f͑r͒, the FT and its inverse read where n is the ͑constant͒ density and J 0 ͑x͒ the zero order Bessel function of the first kind.One way to implement these transformations is to use a finite box in configuration and momentum spaces of length L and K, respectively.The additional conditions f͑r = L͒ = 0 and f͑k = K͒ = 0 lead to a discretized set of coordinates and momenta, related to the zeros j of J 0 ͑x͒ through the relations with j , ␣ =1,2, ... ,N, N being the total number of points in the grids.A Gauss integration rule based on series expansion in Bessel functions and the orthogonality relation, with the integration weights Equations ͑A9͒ turn integrals into algebraic products that can be carried out numerically in a neat and fast way using available standard libraries.

FIG. 1 .
FIG. 1. Scaled energy per particle of the HD gas as a function of x.Solid circles: EL/HNC; open squares: SR/HNC; open diamonds: SR/VMC; stars, pluses, and crosses: low density expansions, KS, CS ͑u͒ , and CS ͑␦͒ , respectively.Notice that diamonds are hardly distinguished from the solid circles.

FIG. 7 .
FIG. 7. Long range behavior of the EL one-body density matrix at several values of x.Solid line: x = 0.01; dashed line: x = 0.005; dotted line: x = 0.001.