Energy of N two-dimensional bosons with zero-range interactions

We derive an integral equation describing $N$ two-dimensional bosons with zero-range interactions and solve it for the ground state energy $B_N$ by applying a stochastic diffusion Monte Carlo scheme for up to 26 particles. We confirm and go beyond the scaling $B_N\propto 8.567^N$ predicted by Hammer and Son [Phys. Rev. Lett. {\bf 93}, 250408 (2004)] in the large-$N$ limit.


I. INTRODUCTION
Nonrelativistic two-dimensional bosons with zerorange interactions is one of the simplest and therefore fundamentally important models in modern science. It has essentially only one parameter -the number of particles N . The other parameter, the dimer binding energy B 2 , can always be set to unity by rescaling the coordinates. This model neighbors the model of onedimensional attractive bosons, which is exactly solvable [1], and the model of three-dimensional bosons, the zerorange formulation of which is ill defined since its groundstate energy is not bound from below because of the Thomas collapse [2].
Bruch and Tjon [3] have shown that no Thomas collapse occurs in two-dimensions. They have found two bound states (ground and first excited) of three two-dimensional bosons. The energies of these states have been calculated over the years by various fewbody methods with increasing accuracy [3][4][5][6][7][8][9]. In particular, the ground state trimer energy is known to be B 3 = 16.5226874B 2 [7,8].
Calculations for the four-boson system have been done using finite-range potentials, but significant range corrections make the convergence to the zero-range limit problematic [10][11][12]. The current best estimate of the tetramer energy B 4 = 197.3(1)B 2 is obtained in the zerorange limit by using the effective field theory [13]. Note the fast increase of the binding energy with N .
Hammer and Son [7] have studied N -body droplets in the limit of large N and found that their energies increase exponentially with N , and sizes decrease proportionally to r −N/2 = 0.3417 N . This result relies on the idea of asymptotic freedom and can be obtained in the mean-field approximation by introducing a logarithmic running of the coupling constant as a function of the droplet size. The proportionality coefficient in Eq. (1) is currently unknown and there is no theory which systematically accounts for finite-N corrections starting from the large-N limit. On the other hand, numerical efforts to approach the large-N regime from the few-body side have been impeded by the rapid growth of the configurational space and the necessity of an extrapolating procedure to remove finite-range effects.
In particular, finite-range calculations of Blume [14] carried out for N ≤ 7 were not conclusive enough to confirm or disprove Eq. (1). Lee [15] calculated energies of up to 10 particles by using lattice effective field theory and observed convergence towards B N /B N −1 = r although with a significant errorbar, r = 8.3 (6). In this paper we derive a many-body integral equation to describe the system directly in the zero-range limit and solve it by using our recently developed stochastic algorithm [16]. We tabulate B N for all N ≤ 26 and claim the values c 1 = −2.06(4) and c 2 = −8(2) for the first terms in the large-N expansion

II. INTEGRAL EQUATION FOR N BOSONS
Consider the system of N two-dimensional bosons interacting with each other via a zero-range potential characterized by the two-body binding energy B 2 = κ 2 0 (we set = m = 1). The N -body wave function at the energy E satisfies the Schrödinger equation supplied with the Bethe-Peierls boundary conditions for each pair i, j in the limit r i → r j Ψ ∝ ln(κ 0 |r i − r j |e γ /2).
The condition (4) substitutes the interaction potential by relating the logarithmically diverging part of the wave function and its regular (constant) part, terms proportional to higher powers of |r i − r j | being omitted. We then introduce the decomposition arXiv:1711.02345v2 [cond-mat.quant-gas] 20 Feb 2018 where φ(r 1 , ..., r N −2 ; r N −1 , r N ) satisfies Eq. (3), but, in contrast to Ψ, is singular only for r N −1 → r N . So far we do not impose any particular boundary condition on φ in this limit, but require that φ(r 1 , ..., r N −2 ; r N −1 , r N ) be symmetric in its first N − 2 arguments. For E < 0 (in this paper we only consider this case) the general form of φ satisfying the above constraints is where F is totally symmetric, K 0 is the decaying Bessel function, and κ(q 1 , ..., 6) corresponds to a linear combination of states of M non-interacting bosons with momenta q 1 , ..., q M and a molecule (described by the wave function K 0 ) with the binding energy κ 2 . We have chosen to work in the center-of-mass frame where the molecule momentum equals − M j=1 q j . The wave function (5) with φ defined by Eq. (6) is symmetric and satisfies Eq. (3) for arbitrary F . This freedom is removed by the boundary condition (4) which, because of the symmetry, can be applied for a single pair, say, N − 1, N . Namely, by considering the limit r N → r N −1 we observe that the logarithmic divergence in the sum of Eq. (5) comes only from the component with i = N − 1, j = N , given explicitly by Eq. (6). The corresponding logarithmic and regular contributions are obtained by substituting in Eq. (6) the asymptotic form K 0 (x) = − ln(xe γ /2). All other components in Eq. (5) contribute only to the regular part of Ψ in the limit r N −1 → r N . Denoting this latter contribution by F(r 1 , ..., r N −1 ) the boundary condition (4) leads to the equation which, after Fourier transforming, explicitly reads Equation (8) is the two-dimensional N -body analog of the three-dimensional three-body equation derived by Skorniakov and Ter-Martirosian [17]. The formulation of the N -body problem based on Eq. (8) provides a few advantages compared to the one based on Eqs. (3) and (4). Very important for us is the fact that Eq. (8) incorporates the zero-range boundary condition and thus requires no zero-range extrapolation procedure. Another feature is that this formulation reduces the configurational space of the problem by one set of single-particle coordinates giving obvious advantages for small N (3 or 4) where Eq. (8) can be solved by deterministic methods. However, this point is not essential for us in this work since we aim at significantly larger N . In Sec. III we develop and apply to Eq. (8) a stochastic method based on the diffusion Monte Carlo (DMC) technique.

III. STOCHASTIC METHOD
In order to solve Eq. (8) we adopt the stochastic method developed by us for calculating binding energies and other characteristics of clusters consisting of up to four identical fermions interacting resonantly with another atom of different mass [16]. The idea of the method is particularly transparent in the currently discussed bosonic case where the function F (q 1 , ..., q M ) is symmetric and can be assumed positive for the ground state. The solution procedure is then based on a stochastic process (diffusion of walkers) in the 2M -dimensional space {q 1 , ..., q M } for which Eq. (8) is the detailedbalance condition and F is (proportional to) the probability density distribution. We now briefly describe the procedure in the bosonic case (for more technical details applicable in general see Ref. [16] and its supplemental material).
We start with introducing a new function which will actually be the probability density distribution of walkers. The proportionality factor relating F and f in Eq. (9) is chosen such that f be normalizable, as this is not necessarily the case for F . Another requirement for the form of this proportionality factor is to make the sampling and branching (see below) tasks analytical and, therefore, fast. The parameter D shifts the typical momentum of the distribution and influences the convergence rate, but not the final result. Rewriting Eq. (8) in terms of f we obtain the integral equation where K and L are abbreviations for bulky but straightforward expressions which we do not write explicitly; K and L correspond, respectively, to the first and second integral on the right hand side of Eq. (8).
We then create an initial distribution of N w walkers in the space {q 1 , ..., q M } and organize an iterative stochastic process of branching and moving them. Namely, at each iteration walkers are subject to two types of operations: • Single-particle moves: A walker with coordinates q 1 , ..., p i , ..., q M is branched on average W i = d 2 q i K(q 1 , ..., q M ; p i ) times and the i-th coordinate of each of the descendants is moved from p i to q i according to the normalized distribution density function K(q 1 , ..., q M ; p i )/W i . Single-particle moves of this type are similar to the ones we dealt with in the fermionic case [16]. Here we have to calculate the normalization integral and sample a product of two two-dimensional Laplacians, where v is half of the walker's total initial momentum. Both tasks are analytical and very fast. This move is repeated over all walkers and over all coordinates.

• Pair moves:
A walker with coordinates q 1 , ..., p i , ..., p j , ..., q M is branched W i,j = d 2 q i d 2 q j L(q 1 , ..., q M ; p i , p j ) times and coordinates i and j are moved, respectively, from p i and p j to q i and q j , where the latter are drawn from L(q 1 , ..., q M ; p i , p j )/W i,j . In fact, due to the delta function in the second term on the right hand side of Eq. (8) and, therefore, in the function L, the value q i + q j = −q 1 − ... − p i − ... − p j − ... − q M and only the difference q = q i − q j has to be sampled. Explicitly, the distribution of q is governed by a product of three Laplacians, The normalization integral is analytical and one can find a very fast rejection-based sampling procedure. This move is repeated over all walkers and over all pairs of coordinates.
We emphasize that fast branching and sampling is essential for the efficient implementation of this method for large N ∼ M , particularly, since the number of pair moves grows proportionally to M 2 . In this context we note that the problem of several ideal bosons resonantly interacting with an impurity is technically easier as there are no pair terms in the corresponding integral equation.
Up to the noise related to the probabilistic branching procedure the walker population at a next iteration is given by M i=1 W i + i<j W i,j summed over all walkers. If this number is systematically larger or smaller than the population at the previous iteration, we can correct the value of E, which implicitly enters in the weights W through the kernels K and L. Alternatively, and this is what we do in practice, we fix E = −1 and use κ 0 as the parameter that controls the walker population and keeps it close to the initial N w . More precisely, we adjust this parameter during every iteration in such a way that the sum of all branching weights of all walkers equals N w .
Assuming that a steady equilibrium distribution of walkers is reached and that κ 0 is fixed at its exact value, it is easy to show that the detailed-balance equation for the above iterative scheme is exactly Eq. (10). With the population correction implemented, we do reach a steady state, but κ 0 fluctuates around a certain average value κ 0 (N w ). The amplitude of these fluctuations decreases with N w and κ 0 (N w ) converges to the exact value in the limit N w → ∞ (see Sec. IV).

IV. RESULTS
The free parameters of our algorithm are D and N w . We find phenomenologically that 2 D 6 maximizes the convergence rate. Results presented below are obtained with D = 4.
For a given walker population N w , after thermalization and accumulation of statistics, our algorithm gives the quantity κ 0 (N w ) at fixed E = −1 with a certain statistical uncertainty. We then translate it into the Nbody energy estimate B N (N w )/B 2 = 1/κ 2 0 (N w ). In order to verify the convergence of B N (N w ) towards B N for a given N we run our code with various values of N w and fit the obtained function B N (N w ) to the form B N (N w ) = B N + cN −γ w , where B N , c, and γ are fitting parameters. The best-fit values of the exponent γ for different N all belong to the interval (0.5, 1). In Fig. 1 we show the quantity 1 − B N (N w )/B N for a few representative values of N together with the corresponding fits. The energies obtained by using this extrapolating procedure are reported in Tab. I. The error bars represent the sum of the statistical uncertainty and systematic uncertainty which we define as |B N (N w,max ) − B N |, i.e., the deviation of the energy calculated with the largest walker pool from the extrapolated value. We should note here that the intrinsic time scale of our stochastic process is related only to the dimension of the configurational space and we benefit from relatively short thermalization and correlation times, at most a few hundreds of iterations for largest N . Thus, at approximately the same CPU cost we have a large room for increasing N w and decreasing the number of iterations, which minimizes the finite-N w systematic error while keeping the statistical uncertainty constant.
We can benchmark our results for N = 3 and 4 against the most accurate previous calculations. For the trimer we have B 3 /B 2 = 16.5225 (2), in good agreement with Our data for larger N are consistent with the value r = 8.567 [7,18] and with the assumption that the function ln B N r −N /B 2 can be expanded in powers of 1/N . In Fig. 2

V. SUMMARY AND OUTLOOK
In this paper we have calculated the energies of twodimensional bosonic droplets with strictly zero-range interactions for N ≤ 26. Our results are consistent with the proportionality factor r = 8.567 in Eq. (1) and we estimate the leading finite-N corrections; we claim B N /B 2 ≈ r N exp(c 1 + c 2 /N ) where c 1 = −2.06(4) and c 2 = −8 (2). That c 2 turns out to be rather large indicates significant finite-N effects for N as large as 10. Our results should be useful as a reference point for developing a theory beyond Ref. [7] expected, in particular, to describe the droplet dynamics and excitations. The calculation of excited states with our current method is challenging since the solution in this case is a sign-changing function.
The exponential scaling of the size and energy raises obvious questions concerning the feasibility of observing such droplets. In the quasi-two-dimensional geometry for small ratio of the scattering length a < 0 to the confinement oscillator length l the dimer size is exponentially large, ∝ l exp( π/2l/|a|) [20]. The N -body droplet size, in this case proportional to l exp( π/2l/|a| − N ln √ r), should be much larger than l and smaller than the waist of the laser beams used for the quasi-two-dimensional confinement. This defines the window of parameters where the droplet is observable; one can see that N is not necessarily very small. Given extremely low values of a routinely produced in experiments (see, for example, [21][22][23]), it is realistic to observe droplets of a few tens of particles. The preparation of a droplet can be realized, for example, by simultaneously switching off of the radial confinement and sweeping the scattering length from a positive to negative value near a zero crossing.
Our method can be generalized to include the leadingorder effective-range correction to the N -body energy (as has been done in three dimensions [16]). One can then estimate the energy correction due to finite-range effects associated with the finite width of the cloud in the quasi-two-dimensional case. In addition, it should provide a link to few-body calculations based on realistic finite-range potentials. Finally, we mention that our method can be extended to bosonic mixtures with different masses.