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

The energy and structure of a dilute hard-disks Bose gas are studied in the framework of a variational many-body approach based on a Jastrow correlated ground state wave function. The asymptotic behaviors of the radial distribution function and the one-body density matrix are analyzed after solving the Euler equation obtained by a free minimization of the hypernetted chain energy functional. Our results show important deviations from those of the available low density expansions, already at gas parameter values $x\sim 0.001$. The condensate fraction in 2D is also computed and found generally lower than the 3D one at the same $x$.


I. INTRODUCTION
The study of two dimensional (2D) quantum manybody 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 twodimensional 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 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 [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 longrange 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 ofh 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 the case of strong interactions, as a hard-core potential, correlation effects are crucial. Within CBF theory they are included by means of a many-body correlation operator, F (1, 2, . . . , N ), acting on the non-interacting ground state wave function, Φ 0 (1, 2, . . . , N ), [8]: with Φ 0 = 1 for homogeneous bosons. F may be determined according to the variational principle by minimizing the ground state energy.
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, constant). 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 equations) are briefly outlined in Section II; an analysis of the EL asymptotic behaviors is presented in Section III; Section IV presents the results for the HD model; while summary and conclusions are given in Section V.

II. VARIATIONAL THEORY
Given the Jastrow correlated wave function (3), the energy per particle of the 2D homogeneous Bose gas is: ∇ 2 ln f 2 (r 12 ) .
(6) 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), g(r 12 ) = f 2 2 (r 12 ) e N (r12)+E(r12) 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), 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 non-zero, longrange 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), = (2π) 2 n n 0 δ(k) + dr e ik·r (n 1 (r) − n n 0 ) 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, 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 straightfordwardly 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 δ δf (r) +µn dr 1 − f 2 (r) = 0 (14) which requires the use of a Lagrange multiplier µ. The minimization is performed under the healing condition 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µ/h 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 fromh(k) as Assuming thatH(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: and 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 lim k→0 kn(k) = n 0 2 mc h . (26)

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 ofh 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 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 V MC SR /E HN C SR , 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 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 Figure (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 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 Figure (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 many-body 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 Figure (3) at x = 0.01 and x = 0.1. The quantity shown is x r 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 [x r(1 − f EL (r))] and [x rN (r)] to enlighten the asymptotic limits (23) and (22) whose dimensionless values are c/4π and c/2π, respectively. The fact that one limit is twice the other is clearly appreciated in the figure. Notice also that due to the chain process implied by the HNC scheme, N (r) merges to the 1/r law at larger distances then f (r). To illustrate the different asymptotic behaviors, we also show [x r(1 − f SR (r))] which goes quickly to zero. The EL static structure function, S(k), is shown in Figure (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 figure (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 Figure (7) by showing the quantity r[n 1 (r)/(n 0 n)−1], whose dimensionless asymptotic value is c/(8πx). 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 Figure (8). The Figure also contains the low-density prediction [12], n LD 0 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.

VI. APPENDIX
In this appendix we discuss the EL equations for an interacting many-body system in 2D. As in the 3D case, the solution to the optimization equation, can be obtained in momentum space, yielding where t(k) =h 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 diagramas in terms of the induced interaction ω I (r). In momentum space, we have Eqs. (36) to (38) 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. (36). 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, can then be built, leading to whith the integration weights Eqs. (43) turn integrals into algebraic products that can be carried out numerically in a neat and fast way using available standard libraries.