Energy of twodimensional bosons with zerorange interactions
Abstract
We derive an integral equation describing twodimensional bosons with zerorange interactions and solve it for the ground state energy by applying a stochastic diffusion Monte Carlo scheme for up to 26 particles. We confirm and go beyond the scaling predicted by Hammer and Son [Phys. Rev. Lett. 93, 250408 (2004)] in the large limit.
I Introduction
Nonrelativistic twodimensional 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 . The other parameter, the dimer binding energy , can always be set to unity by rescaling the coordinates. This model neighbors the model of onedimensional attractive bosons, which is exactly solvable McGuire , and the model of threedimensional bosons, the zerorange formulation of which is ill defined since its groundstate energy is not bound from below because of the Thomas collapse Thomas .
Bruch and Tjon BruTjo79 have shown that no Thomas collapse occurs in twodimensions. They have found two bound states (ground and first excited) of three twodimensional bosons. The energies of these states have been calculated over the years by various fewbody methods with increasing accuracy BruTjo79 ; AdhDelFre88 ; NieFedJen97 ; NieFedJen99 ; HamSon04 ; KarMal06 ; BroKagKla06 . In particular, the ground state trimer energy is known to be HamSon04 ; KarMal06 .
Calculations for the fourboson system have been done using finiterange potentials, but significant range corrections make the convergence to the zerorange limit problematic Tjo80 ; LimNakAka80 ; VraKil02 . The current best estimate of the tetramer energy is obtained in the zerorange limit by using the effective field theory PlaHamMei04 . Note the fast increase of the binding energy with .
Hammer and Son HamSon04 have studied body droplets in the limit of large and found that their energies increase exponentially with ,
(1) 
and sizes decrease proportionally to . This result relies on the idea of asymptotic freedom and can be obtained in the meanfield 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 corrections starting from the large limit. On the other hand, numerical efforts to approach the large regime from the fewbody side have been impeded by the rapid growth of the configurational space and the necessity of an extrapolating procedure to remove finiterange effects. In particular, finiterange calculations of Blume Blu05 carried out for were not conclusive enough to confirm or disprove Eq. (1). Lee Lee06 calculated energies of up to 10 particles by using lattice effective field theory and observed convergence towards although with a significant errorbar, .
In this paper we derive a manybody integral equation to describe the system directly in the zerorange limit and solve it by using our recently developed stochastic algorithm BazPet17 . We tabulate for all and claim the values and for the first terms in the large expansion
(2) 
Ii Integral equation for N bosons
Consider the system of twodimensional bosons interacting with each other via a zerorange potential characterized by the twobody binding energy (we set ). The body wave function at the energy satisfies the Schrödinger equation
(3) 
supplied with the BethePeierls boundary conditions for each pair in the limit
(4) 
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 being omitted.
We then introduce the decomposition
(5) 
where satisfies Eq. (3), but, in contrast to , is singular only for . So far we do not impose any particular boundary condition on in this limit, but require that be symmetric in its first arguments. For (in this paper we only consider this case) the general form of satisfying the above constraints is (we introduce )
(6) 
where is totally symmetric, is the decaying Bessel function, and . Equation (6) corresponds to a linear combination of states of noninteracting bosons with momenta and a molecule (described by the wave function ) with the binding energy . We have chosen to work in the centerofmass frame where the molecule momentum equals .
The wave function (5) with defined by Eq. (6) is symmetric and satisfies Eq. (3) for arbitrary . This freedom is removed by the boundary condition (4) which, because of the symmetry, can be applied for a single pair, say, . Namely, by considering the limit we observe that the logarithmic divergence in the sum of Eq. (5) comes only from the component with , given explicitly by Eq. (6). The corresponding logarithmic and regular contributions are obtained by substituting in Eq. (6) the asymptotic form . All other components in Eq. (5) contribute only to the regular part of in the limit . Denoting this latter contribution by the boundary condition (4) leads to the equation
(7) 
which, after Fourier transforming, explicitly reads
(8) 
We will call Eq. (II) the STM equation since it is the twodimensional body analog of the threedimensional threebody equation derived by Skorniakov and TerMartirosian STM . The formulation of the body problem based on the STM equation (II) provides a few advantages compared to the one based on Eqs. (3) and (4). Very important for us is the fact that Eq. (II) incorporates the zerorange boundary condition and thus requires no zerorange extrapolation procedure. Another feature is that the STM formulation reduces the configurational space of the problem by one set of singleparticle coordinates giving obvious advantages for small (3 or 4) where Eq. (II) can be solved by deterministic methods. However, this point is not essential for us in this work since we aim at significantly larger . In Sec. III we develop and apply to Eq. (II) a stochastic method based on the diffusion Monte Carlo (DMC) technique.
Iii Stochastic method
In order to solve Eq. (II) 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 BazPet17 . The idea of the method is particularly transparent in the currently discussed bosonic case where the function 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 2dimensional space for which Eq. (II) is the detailedbalance condition and 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. BazPet17 and its supplemental material).
We start with introducing a new function
(9) 
which will actually be the probability density distribution of walkers. The proportionality factor relating and in Eq. (9) is chosen such that be normalizable, as this is not necessarily the case for . 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 shifts the typical momentum of the distribution and influences the convergence rate, but not the final result. Rewriting Eq. (II) in terms of we obtain the integral equation
(10) 
where and are abbreviations for bulky but straightforward expressions which we do not write explicitly; and correspond, respectively, to the first and second integral on the right hand side of Eq. (II).
We then create an initial distribution of walkers in the space and organize an iterative stochastic process of branching and moving them. Namely, at each iteration walkers are subject to two types of operations:

Singleparticle moves: A walker with coordinates is branched on average times and the th coordinate of each of the descendants is moved from to according to the normalized distribution density function . Singleparticle moves of this type are similar to the ones we dealt with in the fermionic case BazPet17 . Here we have to calculate the normalization integral and sample a product of two twodimensional Laplacians, , where 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 is branched times and coordinates and are moved, respectively, from and to and , where the latter are drawn from . In fact, due to the delta function in the second term on the right hand side of Eq. (II) and, therefore, in the function , the value and only the difference has to be sampled. Explicitly, the distribution of is governed by a product of three Laplacians, . The normalization integral is analytical and one can find a very fast rejectionbased 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 , particularly, since the number of pair moves grows proportionally to . 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 STM equation.
Up to the noise related to the probabilistic branching procedure the walker population at a next iteration is given by 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 , which implicitly enters in the weights through the kernels and . Alternatively, and this is what we do in practice, we fix and use as the parameter that controls the walker population and keeps it close to the initial . More precisely, we adjust this parameter during every iteration in such a way that the sum of all branching weights of all walkers equals .
Assuming that a steady equilibrium distribution of walkers is reached and that is fixed at its exact value, it is easy to show that the detailedbalance equation for the above iterative scheme is exactly Eq. (10). With the population correction implemented, we do reach a steady state, but fluctuates around a certain average value . The amplitude of these fluctuations decreases with and converges to the exact value in the limit (see Sec. IV).
Iv Results
The free parameters of our algorithm are and . We find phenomenologically that maximizes the convergence rate. Results presented below are obtained with .
For a given walker population , after thermalization and accumulation of statistics, our algorithm gives the quantity at fixed with a certain statistical uncertainty. We then translate it into the body energy estimate . In order to verify the convergence of towards for a given we run our code with various values of and fit the obtained function to the form , where , , and are fitting parameters. The bestfit values of the exponent for different all belong to the interval . In Fig. 1 we show the quantity for a few representative values of together with the corresponding fits. The energies obtained by using this extrapolating procedure are reported in Tab. 1. The error bars represent the sum of the statistical uncertainty and systematic uncertainty which we define as , 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 . Thus, at approximately the same CPU cost we have a large room for increasing and decreasing the number of iterations, which minimizes the finite systematic error while keeping the statistical uncertainty constant.
3  1.65225(2)  15  8.135(2) 
4  1.9720(1)  16  7.129(4) 
5  2.0745(1)  17  6.232(2) 
6  2.0471(1)  18  5.438(3) 
7  1.9462(1)  19  4.734(2) 
8  1.8070(1)  20  4.119(2) 
9  1.6508(4)  21  3.577(2) 
10  1.4905(2)  22  3.108(4) 
11  1.3345(2)  23  2.694(5) 
12  1.1873(4)  24  2.332(4) 
13  1.0508(3)  25  2.018(4) 
14  9.2596(9)  26  1.748(4) 
We can benchmark our results for and 4 against the most accurate previous calculations. For the trimer we have , in good agreement with HamSon04 and KarMal06 . Our tetramer result is more precise than the previously known value PlaHamMei04 .
Our data for larger are consistent with the value HamSon04 and with the assumption that the function can be expanded in powers of . In Fig. 2 we plot this function together with linear, quadratic, and cubic polynomial fits based, respectively, on the data for , , and . The fitting parameters are, respectively, , , and . Based on these findings we claim and . The fact that the
V Summary and outlook
In this paper we have calculated the energies of twodimensional bosonic droplets with strictly zerorange interactions for . Our results are consistent with the proportionality factor in Eq. (1) and we estimate the leading finite corrections; we claim where and . That turns out to be rather large indicates significant finite effects for as large as 10. Our results should be useful as a reference point for developing a theory beyond Ref. HamSon04 expected, in particular, to describe the droplet dynamics and excitations. The calculation of excited states with our current method is challenging since the solution of the STM equation in this case is a signchanging function.
The exponential scaling of the size and energy raises obvious questions concerning the feasibility of observing such droplets. In the quasitwodimensional geometry for small ratio of the scattering length to the confinement oscillator length the dimer size is exponentially large, PetShl2001 . The body droplet size, in this case proportional to , should be much larger than and smaller than the waist of the laser beams used for the quasitwodimensional confinement. This defines the window of parameters where the droplet is observable; one can see that is not necessarily very small. Given extremely low values of routinely produced in experiments (see, for example, Fattori ; Hulet ; Bourdel ), 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 effectiverange correction to the body energy (as has been done in three dimensions BazPet17 ). One can then estimate the energy correction due to finiterange effects associated with the finite width of the cloud in the quasitwodimensional case. In addition, it should provide a link to fewbody calculations based on realistic finiterange potentials. Finally, we mention that our method can be extended to bosonic mixtures with different masses.
Acknowledgment
We thank G. Astrakharchik, J. Carlson, and D. Lee for useful communications and O. Mustaki for technical assistance. This research was supported by the European Research Council (FR7/20072013 Grant Agreement No. 341197) and from the Pazi Fund. We also acknowledge support by the IFRAF Institute.
References
 (1) J. B. McGuire, Study of Exactly Soluble OneDimensional NBody Problems, J. Math. Phys. 5, 622 (1964).
 (2) L. H. Thomas, The interaction between a neutron and a proton and the structure of H, Phys. Rev. 47, 903 (1935).
 (3) L. W. Bruch and J. A. Tjon, Binding of three identical bosons in two dimensions, Phys. Rev. A 19, 425 (1979).
 (4) S. K. Adhikari, A. Delfino, T. Frederico, I. D. Goldman, and L. Tomio, Efimov and Thomas effects and the model dependence of threeparticle observables in two and three dimensions, Phys. Rev. A 37, 3666 (1988).
 (5) E. Nielsen, D. V. Fedorov, and A. S. Jensen, Threebody halos in two dimensions, Phys. Rev. A 56, 3287 (1997).
 (6) E. Nielsen, D. V. Fedorov, and A. S. Jensen, Structure and Occurrence of ThreeBody Halos in Two Dimensions, FewBody Syst. 27, 15 (1999).
 (7) H.W. Hammer and D. T. Son, Universal Properties of TwoDimensional Boson Droplets, Phys. Rev. Lett. 93, 250408 (2004).
 (8) O. I. Kartavtsev and A. V. Malykh, Universal lowenergy properties of three twodimensional bosons Phys. Rev. A 74, 042506 (2006).
 (9) I. V. Brodsky, M. Yu. Kagan, A. V. Klaptsov, R. Combescot, and X. Leyronas, Exact diagrammatic approach for dimerdimer scattering and bound states of three and four resonantly interacting particles, Phys. Rev. A 73, 032724 (2006).
 (10) J. A. Tjon, Bound states of four identical bosons in two dimensions, Phys. Rev. A 21, 1334 (1980).
 (11) T. K. Lim, S. Nakaichi, Y. Akaishi, and H. Tanaka, Nearthreshold behavior of the groundstate binding energies of the fewatom systems of He and other bosons in two and three dimension, Phys. Rev. A 22, 28 (1980).
 (12) L. Vranješ and S. Kilić, Helium trimers and tetramers in two dimensions and quasitwodimensions, Phys. Rev. A 65, 042506 (2002).
 (13) L. Platter, H.W. Hammer, and U.G. Meißner, Universal Properties of the FourBoson System in Two Dimensions, FewBody Syst. 35, 169 (2004).
 (14) D. Blume, Threshold behavior of bosonic twodimensional fewbody systems, Phys. Rev. B 72, 094510 (2005).
 (15) D. Lee, LargeN droplets in two dimensions, Phys. Rev. A 73, 063204 (2006).
 (16) B. Bazak and D. S. Petrov, FiveBody Efimov Effect and Universal Pentamer in Fermionic Mixtures, Phys. Rev. Lett. 118, 083002 (2017).
 (17) G. V. Skorniakov and K. A. TerMartirosian, Three Body Problem for Short Range Forces. I. Scattering of Low Energy Neutrons by Deuterons, Zh. Eksp. Teor. Fiz. 31, 775 (1956) [Sov. Phys. JETP 4, 648 (1957)].
 (18) D. S. Petrov and G. V. Shlyapnikov, Interatomic collisions in a tightly confined Bose gas, Phys. Rev. A 64, 012706 (2001).
 (19) M. Fattori, C. D’Errico, G. Roati, M. Zaccanti, M. JonaLasinio, M. Modugno, M. Inguscio, and G. Modugno, Atom Interferometry with a Weakly Interacting BoseEinstein Condensate, Phys. Rev. Lett. 100, 080405 (2008).
 (20) J. H. V. Nguyen, P. Dyke, D. Luo, B. A. Malomed, and R. G. Hulet, Collisions of matterwave solitons, Nat. Phys. 10, 918 (2014).
 (21) S. Lepoutre, L. Fouché, A. Boissé, G. Berthet, G. Salomon, A. Aspect, and T. Bourdel, Production of strongly bound K bright solitons, Phys. Rev. A 94, 053626 (2016).