Finite-range effects in ultradilute quantum drops

In the first experimental realization of dilute Bose-Bose liquid drops using two hyperfine states of $^{39}$K some discrepancies between theory and experiment were observed. The standard analysis of the data using the Lee-Huang-Yang beyond mean-field theory predicted critical numbers which were significantly off the experimental measurements. Also, the radial size of the drops in the experiment proved to be larger than expected from this theory. Using a new functional, which is based on quantum Monte Carlo results of the bulk phase incorporating finite-range effects, we can explain the origin of the discrepancies in the critical number. This result proves the necessity of including finite-range corrections to deal with the observed properties in this setup. The controversy on the radial size is reasoned in terms of the departure from the optimal concentration ratio between the two species of the mixture.


INTRODUCTION
In the last years, there has been an increasingly high interest in understanding dilute, ultracold quantum Bose-Bose mixtures. The focus on this study increased dramatically after the theoretical proposal by Petrov [1] on the formation of self-bound liquid drops. These liquid drops are stabilized by beyond-mean-field effects and can appear in mixtures of two Bose-Einstein condensates with repulsive intraspecies and attractive interspecies interactions. The drops originate from a delicate balance between the collapsed state, predicted by mean-field (MF) theory, and the repulsive character of the first beyond mean-field term (Lee-Huang-Yang -LHY-). The same perturbative theoretical scheme predicts self-binding in low-dimensional mixtures [2,3] and dipolar systems [4,5]. Recently, these predicted quantum drops have been observed in several experiments [6][7][8][9] and they resemble the well-known liquid Helium drops [10,11]. However, the inner density in the Bose-Bose drops is about five orders of magnitude smaller than in 4 He [10,11]. Therefore, these new quantum drops extend the realm of the liquid state to much lower densities than any previous existing classical or quantum liquid.
In the two labs [6,8] where the drops have been observed, the Bose-Bose mixture is composed of two hyperfine states of 39 K. In the first experiment by Cabrera et al. [6], the drops are harmonically confined in one of the directions of space whereas in the second one by Semeghini et al. [8] the drops are observed in free space. This difference in the setup makes that in the first case the drops are not spherical like in the second experiment. This also affects the critical number, that is, the minimum number of atoms required to get a self-bound state. The measured critical numbers differ significantly between the two labs due to the different shape of the drops, the ones in the confined case being smaller than in the free case. In both works, the experimental results for the critical number are compared with the MF+LHY theory. The agreement between this theory and the drops produced in free space is quite satisfactory in spite of the large errorbars of the experimental data that hinder a precise comparison. However, in the confined drops of Ref. [6], where the critical numbers are significantly smaller than in the free case, the theoretical predictions do not match well the experimental data.
Ultradilute liquid drops, which require beyond-mean field corrections to be theoretically understood, offer the perfect benchmark to explore possible effects beyond MF+LHY theory [12] which usually play a minute role in the case of single-component gases [13,14]. Indeed, several theoretical studies [15][16][17][18][19] indicate a strong dependence of the equation of state of the liquid on the details of the interatomic interaction, even at very low densities. This essentially means it is already possible to achieve observations outside the universal regime, in which all the interactions can be expressed in terms of the gas parameter na 3 , with a the s-wave scattering length. The first correction beyond this universality limit must incorporate the next term in the scattering series, that is the effective range r eff [20,21], which in fact can be quite large in these drops and in alkali atoms in general [22,23].
Motivated by experiments with quantum drops, we have investigated the self-bound quantum mixture composed of two hyperfine states of 39 K using nonperturbative quantum Monte Carlo (QMC) methods. Direct QMC simulations [24] of finite particle-number drops, as produced in experiments, would serve as a great test of mean-field theory but, unfortunately, this is not yet achievable because of the large number of particles in realistic drops (N > 10 4 ). Yet, the problem can be addressed in the Density Functional Theory (DFT) spirit, relying on the Hohenberg-Kohm-Sham 2nd theorem [25], which guarantees that a density functional exists that matches exactly the ground-state solution. To build a functional for the quantum Bose-Bose mixture, we have carried out calculations in bulk conditions using the diffusion Monte Carlo (DMC) method, an exact QMC technique applicable to systems at zero temperature. Using that functional we can access to energetics and structure of the liquid drops in the same conditions as in the experiment. We focus on the data obtained by Cabrera et al. [6] in the confined setup since it is in that case where discrepancies between MF+LHY theory were observed.
The rest of the paper is organized as follows. In Sec. II, we introduce the theoretical methods used for the study and discuss the way in which the density functional is built. Sec. III comprises the results obtained for the bulk liquid using the available scattering data of the 39 K mixture. The inclusion of the effective range parameters in the interaction model allows for a better agreement with the measured critical numbers. Finally, we summarize the most relevant results here obtained and derive the main conclusions of our work.

METHODS
We study a mixture of two hyperfine states of 39 K bosons at zero temperature. The Hamiltonian of the system is where V (α,β) (r iαj β ) is the interatomic potential between species α and β. The mixture is composed of N = N 1 + N 2 atoms, with N 1 (N 2 ) bosons of type 1 (2). The potentials are chosen to reproduce the experimental scattering parameters, and we have used different model potentials to investigate the influence of the inclusion of the effective range. The microscopic study has been carried out using a second-order DMC method [26], which allows for an exact estimation of the ground-state of the mixture within some statistical errors. DMC solves stochastically the imaginary-time Schrödinger equation using a trial wave function as importance sampling which guides the diffusion process to regions of expected large probability. In the present case, we used a trial wavefunction built as a product of Jastrow factors [27], where the two-particle correlation functions f (r) are The function f 2b is the solution of the two-body problem for a specific interaction model, R 0 is a variational parameter, and L = (N/ρ) 1/3 is the size of the simulation box. There is a weak dependence of the variational energy on R 0 , and it has been kept as R 0 = 0.9L/2 for all the cases. A careful analysis of imaginary time-step dependence and population size bias has been carried out, keeping both well under the statistical error. The timestep dependence is well eliminated for ∆τ = 0.2×ma 2 11 / and the population bias by using n w = 100. Our simulations are performed in a cubic box with periodic boundary condition, using a number of particles N . The thermodynamic limit is achieved by repeating calculations with different particle numbers; we observe that, within our numerical precision, the energy per particle converges at N ≈ 600 for the range of magnetic fields here considered.
Within density functional theory (DFT), we seek for a many-body wave function built as a product of singleparticle orbitals, (4) These single-particle wave functions, which in general are time-dependent, are obtained by solving the Schrödingerlike equation [28], where V ext is an external potential acting on the system and E int is an energy per volume term that accounts for the interparticle correlations. The differential equation (5) is solved by propagating the wave function ψ with the time-evolution operator To this end, we have implemented a three-dimensional numerical solver based on the Trotter decomposition of the time evolution operator [29,30] with second-order accuracy in the timestep ∆t, as follows with K and V the kinetic and potential terms in Eq. 5.

RESULTS
In order to go beyond the MF+LHY density functional we have carried out DMC calculations of the bulk liquid. In the mixture of 39 K under study, we call the state |F, m F = |↓ = |1, 0 as component 1, and the state |F, m F = |↑ = |1, −1 as component 2. In Fig. 1, we show the energy per particle of the 39 K mixture as a function of the density, using three different sets of potentials in the Hamiltonian (1): i) Hard-core interactions (HCSW) with diameter a ii , i = 1, 2, for the repulsive intraspecies interaction, and a square-well potential with range R = a 11 and depth V 0 for the interspecies potential. The three potentials reproduce the s-wave scattering lengths for the three channels, ii) POT1 stands for a set of potentials which reproduces both the s-wave scattering lengths and effective ranges of the three interacting pairs of the 39 K mixture. To model the interactions, we have chosen a square-well square barrier potential [31] for the 11 channel, a 10-6 Lennard-Jones potential [32] for the 22 channel, and a square-well potential of range R and depth V 0 [13] in the 12 channel, iii) POT2 also reproduce both the s-wave scattering lengths and effective ranges, by using a sum of Gaussians in the 11 channel, a 10-6 Lennard-Jones potential in the 12 channel, and finally a soft-sphere square well in the 22 channel.
In all cases, the attractive interatomic potential does not support a two-body bound state. We have obtained the s-wave scattering length and effective range of the potentials using standard scattering theory [20,21].
We compare our DMC results to the MF+LHY theory, which can be compactly written as [1] assuming the optimal concentration of particles from mean-field theory, N 1 /N 2 = a 22 /a 11 . The energy per particle E 0 /N at the equilibrium density of the MF+LHY approximation ρ 0 and ρ 0 itself is In Fig. 1, we report DMC results for the equation of state corresponding to a magnetic field B = 56.337 G, one of the magnetic fields used in experiments. We show the convergence of the results on the number of particles in the simulation for the particular case of POT1 set of potentials. As we can see, the convergence is achieved with N = 600. We have repeated this analysis for all the potentials and, in all the magnetic field range explored, we arrive to convergence with similar N values. We have investigated the dependence on the effective range by repeating the calculation using the HCSW and POT2 potentials. As it is clear from Fig. 1, only when both scattering parameters, the s-wave scattering length and the effective range, are imposed on the model potentials we get an approximate universal equation of state, mainly around the equilibrium density. The equation of state so obtained shows a significant and overall decrease of the energy compared to the MF+LHY prediction, with a correction that increases with the density. Instead, using the HCSW potentials, which only fulfill the s-wave scattering lengths, the energies obtained are even above the MF+LHY prediction. A similar behavior has been previously shown to hold in symmetric (N 1 = N 2 ) Bose-Bose mixtures [16].
Equations of state of the bulk mixture, for the seven values of the magnetic field used in the experiments (B = 56.230 G to B = 56.639 G), are shown in Fig. 2. The DMC results are calculated using the model POT1, but the differences with the other set POT2 are not significant. In all cases, we take the mean-field prediction for the optimal ratio of partial densities ρ 1 /ρ 2 = a 22 /a 11 . We have verified in several cases that this is also the concentration corresponding to the ground state of the system in our DMC calculations, i.e., the one that gives the minimum energy at equilibrium. The DMC results are compared with the MF+LHY equation of state (8).
Overall, a reduction of the magnetic field, or equivalently an increase in |δa| = a 12 + √ a 11 a 12 , leads to an increase of the binding energy compared to the MF+LHY approximation. This happens clearly due to the influence of the large experimental effective range, since in the limit of zero range one would observe overall repulsive beyond-LHY terms (see also Fig. 1 and Ref. [16]). DMC energies for the 39 K mixture are well fitted using the functional form as it can be seen in Fig. 2 3. Dependence of the critical atom number on the magnetic field. Full circles are predictions using the QMC functional within DFT with the interaction potentials which reproduce both a and r eff . Diamond points are data from the experiment [6]. Empty points show the prediction using the QMC functional with the HCSW model potentials.  Fig. 3 in comparison with the experimental results of Ref. [6]. To make the comparison reliable, we have included the same transversal confinement as in the experiment. In particular, theoretical predictions are obtained within DFT, using a Gaussian ansatz φ = exp −r 2 /(2σ the points at B = 56.23 G, B = 56.337 G and B = 56.453 G in Fig. 3). The observed decrease of N c leads our theoretical prediction closer to the experimental data in a significant amount and in all the δa range, clearly showing the significant influence of the effective range on the N c values. Experiments on quantum droplets were performed either in the harmonic trap [6] or in a free-drop setup [8]. Predictions of N c for these two geometries are given in Tables II and III, using MF+LHY and QMC functionals. The absolute difference of predicted N c values between the two functionals are about 1000 atoms. On the other hand, the relative difference is much higher in the harmonically-trapped system because the presence of an external trap significantly reduces N c .
A second observable measured in experiments is the size of the drops. The radial size of a N = 15000 drop for different values of the magnetic field was reported in Ref. [6]. In Fig. 4, we compare the experimental values with different theoretical predictions. We observe a slight reduction in size using QMC functionals, compared to MF+LHY theory, which is a consequence of the stronger binding produced by inclusion of finite range interactions. Since the experimental data go to the opposite direction, it means that drops size can not be explained solely in terms of the non-zero effective range. One possible explanation for this clear disagreement could be a deviation from the optimal relative number of particles, which can occur in non-equilibrated drops or when one of the components has a large three-body recombination coefficient. Let us define x = N 2 /N 1 a 22 /a 11 . Then, x = 1 stands for the optimal relative particle number, i.e., the concentration corresponding to the ground-state of the system. We have investigated the behavior of both the MF+LHY and QMC functionals under variations in x, and both predict a decrease in the drop size proportional to the deviation from x = 1. Using the MF+LHY functional, we have obtained the x values that fit the experimental size for every B (Fig. 4). We report the result of this analysis in Fig. 5; notice that there is a symmetry on x and so only its absolute deviation from one is important. This result clearly shows the sensitive dependence of drop structural properties on the relative atom number.
As we can see in Fig. 5, the value for x becomes 1 (optimal value) when the drop composed by 15000 particles is studied at the highest magnetic field. This can be understood if we observe that the critical number for this magnetic field matches approximately this number of atoms (see Fig.3 7. Dependence of the radial size σr on the number of particles. The size is obtained from the variational ansatz, since close to the critical atom number the density profile in the radial direction is well approximated by a Gaussian. In both functionals, it is assumed that the relative concentration is optimal N2/N1 = a11/a22. QMC functional includes the correct finite-range r eff through POT. 1 set of potentials, Fig.  1. is larger than the critical number (lower B in Fig. 4) x departs from one. This can be better understood if one calculates the drop phase diagram as a function of x. The result is plotted in Fig. 6. As the number of particles is approaching the critical one, the range of possible values of x, which support a drop state, is reducing. This is a supporting fact that drops close to the critical atom number observed in the experiment fulfill the condition x = 1. On the other hand, there is an increasing range of relative particle concentrations for which a drop can emerge as the number of particles increases.
Close to the critical atom number, the density profile of a drop can change drastically depending on the functional. We illustrate this effect in Fig. 7 for a magnetic field B = 56.337 G. In the figure, we show the dependence of the radial size on the number of particles, with the same harmonic confinement strength as in one of the experiments [6]. We observe a substantial difference between the MF+LHY and QMC functional results, mainly when N approaches the critical number N c . DISCUSSION Experiment in Ref. [6] showed a significant disagreement between the measured data and the MF+LHY perturbative approach. In order to determine the possible origin of these discrepancies we have pursued a beyond MF+LHY theory which incorporates explicitly the finite range of the interaction. To this end, we have carried out DMC calculations of the bulk liquid to estimate accurately its equation of state. We have observed that the inclusion in the model potentials of both the s-wave scattering length and the effective range produces a rather good universal equation of state in terms of these pair of parameters. Excluding the effective range, significant differences are obtained from these universal results. This relevant result points to the loss of universality in terms of the gas parameter in the study of these dilute liquid drops.
Introducing the DMC equation of state into the new functional, following the steps which are standard in other fields, such as DFT in liquid Helium [10], we derive a new functional that allows for an accurate study of the most relevant properties of the drops. In particular, we observe that the inclusion of finite range effects reduces the critical atom number in all the magnetic field range approaching significantly the experimental values. On the other hand, our QMC functional is not able to explain the clear discrepancy between theory and experiment about the size of the drops. We attribute this difference to the dramatic effect on the size that small shifts on the value of x produce. Our analysis provides a reasonable explanation of this feature: above the critical atom number the window of stability of the drops increases from the single point x = 1 to a range of values that, in absolute terms, grow with the number of particles. With the appropriate choice of x, one can obtain agreement with the experiment.
The drops produced in the different setup of Ref. [8] are spherical since all magnetic confinement is removed. The corresponding critical numbers in this case are larger than in the confined setup [6] and MF+LHY theory accounts reasonably well for the observed features. We have applied our formalism also to this case and the corrections are not zero but relatively less important than in the case analyzed here.
We acknowledge very fruitful discussions with Leticia Tarruell, César Cabrera, and Julio Sanz. This work has been supported by the Ministerio de Economia, Industria y Competitividad (MINECO, Spain) under grant No. FIS2017-84114-C2-1-P. V. C. acknowledges financial support from STSM Grant COST Action CA16221.