Dipolar droplets of strongly interacting molecules

Dipolar droplets of strongly interacting molecules Tim Langen,1, 2 Jordi Boronat,3Juan Sánchez-Baena,3Raúl Bombín,4Tijs Karman,5and Ferran Mazzanti3 15. Physikalisches Institut and Center for Integrated Quantum Science and Technology, Universität Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany 2Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria∗ 3Departament de Física, Universitat Politecnica de Catalunya, Campus Nord B4-B5, 08034 Barcelona, Spain 4Institut des Sciences Moléculaires (ISM), Université de Bordeaux, 351 Cours de la Libération, 33405 Talence, France 5Institute for Molecules and Materials, Radboud University, Nijmegen, The Netherlands (Dated: July 16, 2024) We simulate a molecular Bose-Einstein condensate in the strongly dipolar regime, observing the existence of self-bound droplets, as well as their splitting into multiple droplets by confinementinduced frustration. Our quantum Monte Carlo approach goes beyond the limits of the established effective mean-field theories for dipolar quantum gases, revealing small droplets produced by strong dipolar interactions outside known stable regimes. The simulations include realistic molecular interactions and therefore have direct relevance for current and future experiments. Introduction. Ultracold dipolar Bose-Einstein condensates (BECs) are known to exhibit many exotic quantum phases [1,2], ranging from self-bound droplets [3,4] to droplet supersolids [5–7] and exotic supersolid patterns [8–10]. A key role in the formation, stability, and dynamics of these phases is played by repulsive quantum fluctuations, which can stabilize the dipolar systems against a mean-field collapse [11,12]. In the theoretical modelling, these fluctuations are usually taken into account using an effective mean-field description based on the extended Gross-Pitaevskii equation (eGPE) [13–15]. So far, the experimental study of dipolar gases has mostly been based on magnetic BECs, formed by weakly dipolar lanthanide atoms [2]. In these cases, the effective description yields — except for exceptional cases [16–18] — a good description of the physics observed. However, with the recent dramatic advances [19] in preparing strongly dipolar molecular gases close to or in the quantum regime [20–27] it can be anticipated that this description will soon reach its limits [10,28]. Moreover, on a more fundamental level, it is expected that strongly interacting molecular gases will share important similarities with solid helium [10], where vacancy-induced forms of supersolidity have been investigated for decades [29]. It is thus essential to explore the limits of the current effective mean-field approaches and extend the theoretical description of dipolar quantum matter into the strongly interacting limit. Here, we use quantum Monte Carlo simulations to investigate the existence of dipolar droplets in a molecular BEC. Using realistic interaction potentials for microwave-shielded molecules, we observe that ∗[email protected];[email protected] r (103 a0) r (r0) V(r,θ)/E0 r (r0) r (103 a0) r (r0) r (103 a0) δr = 0.1 δr = 0.5 δr = 1 θ = 0 θ = π/4 θ = 54.7º θ = π/3 θ = 2π/5 θ = π/2 Figure 1. Molecular interaction potentials for single frequency microwave dressing of NaK molecules with a Rabi frequency of Ω/ℏ= 2π×10 MHz and varying relative detuning δr[30]. The dressing establishes a shielding core located at around 1000 a0, which prevents molecular two-body losses. Here, a0denotes the Bohr radius and r0the dipolar length, as defined in the text. By changing the dressing parameters, the potentials can be tuned and various interaction parameters can be realized, potentially resulting in different many-body states of the ultracold molecular gas. self-bound droplets persist in the strongly interacting regime — well outside the regime of validity of meanfield descriptions — for realistic interaction parameters achievable in experiments. Moreover, we observe the formation of metastable droplet assemblies upon confinement-induced frustration, outlining a path towards molecular supersolids. Taken together, this highlights the prospects of molecular BECs to explore novel regimes of dipolar quantum matter. Molecular interactions. Ultracold molecular gases are known to be strongly affected by twoand three-body losses at ultracold temperatures. These losses can, for example, occur through chemical reactions that are alarXiv:2407.09391v2 [cond-mat.quant-gas] 15 Jul 2024 2 δrC6(10−3)deff (D) as(a0)add (a0)ϵdd Ncrit (I) 0.10 0.9080 0.7813 -9558.3 3579.9 -0.375 7 (II) 0.50 1.5471 0.7023 -3644.6 2892.5 -0.794 12 (III) 0.75 2.7026 0.6282 -1666.1 2314.0 -1.389 31 (IV) 1.00 5.0097 0.5552 -596.6 1807.8 -3.030 135 (V) 1.25 9.3089 0.4905 -110.1 1411.0 -12.820 400 Table I. Summary of the molecular interaction parameters considered in this work. We use NaK molecules and a fixed Rabi frequency of Ω/ℏ= 2π×10 MHz. By changing the relative detuning δr=|δ|/Ω, the C6coefficient of the interaction, effective dipole moment deff ,s-wave scattering length asand mean-field dipolar length add, and dipolar parameter ϵdd =add/ascan be tuned in experiments. The critical particle number Ncrit is determined through our Monte Carlo simulations and marks the transition between self-bound droplets and regular gases. Notice the critical value of C6required for a two-body bound state to appear is Ccrit 6∼0.531 ×10−3. For C6> Ccrit 6no two-body bound state is formed, while droplets still exist. lowed for some species [31], or through the formation of short-lived complexes that can be lost from the trapping potential confining the molecules [32]. In contrast to ultracold atoms it is thus imperative to shield molecules from these losses at short range. A very versatile method to achieve this is microwave dressing [33–37], which has been used to realize both degenerate Fermi gases [21] and, very recently, also a BEC of dipolar molecules [27]. An added benefit of such shielding is the tunability of molecular interactions through a simple change of microwave power and frequency. Moreover, as the resulting interaction potentials are formed from microwave-dressed rotational states of the molecules, they are known with high precision and thus ideally suited for numerical simulations, such as the ones presented here. This is in stark contrast to magnetic atoms [2], where the inter-atomic scattering is chaotic and the shape of the interaction potential is thus highly uncertain [16]. In this Letter, we use a recently proposed analytical approximation to model the microwave-dressed molecular interactions [30]. This approximation has been shown to agree very well with the full molecular interaction potentials over the parameter range relevant in typical experiments. We note that our simulations can easily be adapted to describe more complex shielding potentials such as, for example, the recently introduced double microwave dressing [27], DC shielding [38], or combinations of MW and DC fields [10,39]. Expressed in spherical coordinates and for circularly polarized microwaves, the interaction potential of two molecules separated by a distance vector ris given by V(r) = C6 r6sin2θ(1 + cos2θ) + C3 r3(3 cos2θ−1) ,(1) with θthe declination angle. In this expression, C3= d2/[48 πϵ0(1 + δ2 r)] and C6=d4/[128 π2ϵ2 0Ω(1 + δ2 r)3/2] are functions of the relative detuning δr=|δ|/Ω. Here, δand Ωare the detuning and Rabi frequency of the microwave field, respectively. Furthermore, ddenotes the permanent dipole moment in the molecular frame and ϵ0is the vacuum permeability. The resulting effective dipole moment in the lab frame is given by deff =d/p12(1 + δ2 r), and oriented in the z-direction, perpendicular to the microwave polarization. Notice that the dipolar contribution ∼C3/r3has a reversed —anti-dipolar — form [2,40,41]. The corresponding change in sign has a direct influence on the final shape of the system as we discuss below. As a concrete example, in this work we use parameters for bosonic NaK molecules [42], with a dipole moment of 2.72 D, mass m= 62 amu and rotational constant B/ℏ= 2π×2.089 GHz. For the Rabi frequency, we use Ω/ℏ= 2π×10 MHz, which is well-achievable in typical experiments [21]. These parameters describe an interaction potential where bound states are located far from threshold, which is expected to minimize losses in experiments [24]. Example plots of the interaction potential are depicted in Fig. 1for different declination angles θand different interaction strengths. The plots highlight, in particular, the large size of the repulsive shielding core ∼1000 a0, which is in stark contrast to the usual contact interaction potential used to describe ultracold atoms. Besides the effective dipole moment deff , changes in interaction strength are reflected in the s-wave scattering length as, which is expected to be the only relevant scattering parameter in a dilute bosonic gas at ultracold temperatures. Relating this last quantity to the potential parameters requires solving the low-momentum limit of the scattering T-matrix as in Ref. [43], or equivalently the long-range behavior of the s-wave component of the wave function of the E= 0 two-body problem. For the given potential, these calculations are more involved than usual, due to the coupling between modes l, l±2and l±4coming from the purely repulsive part of the interaction. The resulting values of as, and related to this, dipolar parameters ϵdd =add/as, are given in Table I for typical parameters used in this work. Here, add =m d2 eff /[12πϵ0ℏ2]denotes the dipolar length. 3 In view of the large shielding core and strong dipolar interactions of the molecules, it is an interesting question whether the scattering length remains the relevant low-energy parameter that fully characterizes the interactions, whether renormalization of its strength is required as in dysprosium atoms [16,44], or whether a new approach is necessary altogether [45–47]. In this work, we therefore provide both purely experimental parameters such as detuning and Rabi frequency, in parallel with the conventional scattering length. Simulations. Based on the given realistic interaction potential, in the following we analyze a system of N indistinguishable dipolar bosonic molecules in three dimensions described by the Hamiltonian H=−ℏ2 2m N X j=1 ∇2 j+X i<j V(rij),(2) where V(r)is the two-body interaction potential given in Eq. (1). We use the Path Integral Ground State (PIGS) algorithm [48,49] to simulate the ground state of the system. It starts from a variational wave function that is propagated in imaginary time to eliminate all the components that are orthogonal to the true ground state, provided the propagation is carried out over a large enough imaginary time, and a good short-time approximation of the time-evolution operator is used. Notice that this is an ab initio method that leads to the exact ground-state solution for a given Hamiltonian. Remarkably and in contrast to previous works [16,50–52], in the current case an accurate description of the Hamiltonian is known and given by the expressions above, thus allowing for a robust simulation of the molecular BEC’s behavior. The fact that the Hamiltonian is known, makes quantum Monte Carlo methods, such as PIGS, an optimal choice to study this system. In our simulations, we use dimensionless quantities expressed in terms of the dipolar length r0=mC3/ℏ2 and energy E0=ℏ2/mr2 0. Note the difference of 1/3in the numerical pre-factor of r0when comparing to the alternative dipolar length add, which is used in the meanfield description of dipolar gases. In our units convention, C3is always 1, while the changes in the interaction parameters modify the C6parameter. Two-body bound states could appear when the repulsive component of the interaction is reduced below a critical value C6≈0.531 ×10−3, which is, however, outside the set of parameters considered here. In practice, we start the simulation by setting up the system in a tightly confining trap in order to guarantee all particles are close enough, and subsequently let the system equilibrate. Once this has been achieved, the trap is removed and the system is left to evolve freely. When present in the simulations, droplets are 2 1 0 1 2 1 2 3 Position (r0) 0 Column Density (r0 ) -1 1500 1000 500 2 1 0 -1 -2 2 1 0 1 2 x (r0)y (r0) z (r0) ---- (a) Self-bound droplet Gas Critical particle number Ncrit Interaction strength (b) 0.12 relative detuning δr 0.21 0.37 0.61 0.91 1.29 (I) (II) (III) (IV) (V) Figure 2. (a) Single self-bound dipolar droplet in the strongly interacting regime for δr= 0.75 and N= 300 molecules. Due to the anti-dipolar interactions between the molecules, the density distribution is strongly prolate. The column densities have been obtained from the average of a large number of configurations, normalized to the total number of molecules. The inset shows a single snapshot of a droplet. (b) Phase diagram for the transition from self-bound droplets to a normal gas phase, as a function of relative detuning δr. For reference, we also show the interaction strength parameterized by the dipolar parameter ϵ−1 dd =as/add, as it is commonly used to characterize the mean-field limit [16]. Parameters of the white data points (I) to (V) are summarized in Table I. recognized by both their energetic and structural properties: being self-bound states, the droplets total energy remains negative, while at the same time particles stay close to each other and do not spread over the whole simulation box. These two are clear signatures of droplet formation and allow for a good estimation of the critical atom number required for one or more droplets to form. Results. In Fig. 2a we show example results for a simulation with δr= 0.75 and N= 300 particles, leading 4 to a stable dipolar droplet. Due to the specific form of the interaction, which has a fully repulsive short-ranged part combined with the sign-reversed dipolar term, the system tends to arrange in the plane perpendicular to the polarization of the dipoles. The observed droplet thus shows a prolate shape, with a radial size that is significantly larger than its minor axis. This shape is reminiscent of other, previously predicted anti-dipolar droplets in the mean-field regime [41,53,54]. Note that as the interaction potentials do not support a two-body bound state for the parameters considered, the droplets formed are the true manifestation of a many-body selfbound state [28]. A radial integration of the total density shows that the local gas parameter at the center of the system is n|as|3≈4.36, well beyond the reach of mean-field theory [55–57]. Remarkably, this value yields a mean interparticle spacing of n−1/3∼50 nm, which is comparable to the size of the shielding core (cf. Fig. 1). This suggests that the short range details of the interaction, which are neglected in a mean-field approach, play an important role in the description of the droplet. The fact that both length scales — the size of the shielding core and the mean interparticle spacing — are comparable could also be expected to lead to threedimensional structure formation, where the repulsively interacting molecules would spontaneously arrange in a crystal with the same lattice spacing [39]. However, we do not find any evidence for this in our simulations for the parameters investigated. Finally, the observed density ∼6×1021 m−3in the droplet is higher than the density in the molecular BEC and comparable to the densities observed in droplets of magnetic atoms [1]. This emphasizes the need to reduce three-body losses in experimental implementations to a value that is also comparable to that in atoms, where loss rate coefficients L3= 1.3×10−41 m6s−1have been determined [3]. Next, we systematically map out the behavior of the system for different numbers of particles Nand interaction parameters. Our results are shown in Fig. 2b. The parameters for particular values of Ncrit are also summarized in Table I. Similar to the behavior in magnetic atoms, where a subtle balance between dipolar interactions, contact interactions, and quantum fluctuations is required for stable droplets to form [1], we observe a sharp transition between a parameter regime where self-bound droplets are prevalent, and a gaseous BEC. Our results confirm suggestions [10] indicating that for increasing dipolar interaction strength the system supports droplets of decreasing particle number. Specifically, we find critical particle numbers as low as Ncrit = 7for δr= 0.10. The corresponding droplets are characterized by strongly attractive scattering lengths, which is in stark contrast to the usual mean-field description, 2 1 0 -1 -2 -2 -1 0 1 2 3 1.5 1.0 0.5 0 -0.5 -1.0 -2.5 x (add) -3 -2 -1 0 1 2 3 (a) (b) y (add) x (add)y (add) z (add) z (add) -3 -2 -1 0 1 2 -2 -1 0 1 Figure 3. Snapshots of droplets formed from simulations containing (a) N= 90 and (b) N= 300 particles for δr= 0.75, equilibrated in tight confining traps that frustrate the single droplet ground state. The droplets remain in a metastable state after the trap is removed. where such parameters lead to a collapse of the system. Notably, the low particle numbers required are well within the range of the moderate-sized condensates that have been achieved in experiments to date [27]. For larger detunings of the microwave field, corresponding to weaker dipoles, the critical particle number increases, eventually approaching values of several hundred particles, which is similar to values previously observed in weakly-dipolar dysprosium gases [10,16]. Next we compare the results of the quantum Monte Carlo simulations to the predictions of the standard eGPE description, which requires the inclusion of a LeeHuang-Yang (LHY) term to describe the quantum fluctuations that stabilize the droplets. Due to the reversed sign of the dipolar term in the interaction, however, the LHY correction term has to be rederived (see Appendix). The resulting critical atom numbers, obtained as the minimum number of particles required to produce a negative energy, are far below the ones obtained from the PIGS simulations. For δr= 1.25,1.5,1.75 and 2we find N(eGPE) crit = 9,31,69 and 193, respectively. The large discrepancies here are not surprising, considering the large central density of the droplets, together with the fact that with the interaction of Eq. (1) and the conditions considered, the LHY correction acquires an imaginary part that can be as large as a 50%. This again indicates that the system is well beyond a meanfield regime, even when the leading quantum fluctuation terms are taken into account. Despite the fact that in its ground state the system forms a single self-bound droplet, different metastable configurations can be realized by placing it in a more strongly confining trap. When the oscillator lengths associated to the radial frequencies are smaller than the radial size of the ground state droplet, frustration occurs and the system splits, producing two or more nucleation centers depending on the total number of particles and the trap frequencies employed in the initial equilibration [58]. We observe a critical maximum size that each of these droplets can afford, favouring the splitting 5 into more droplets when the total number of particles increases. An example of that is shown in Fig. 3, where two snapshots of equilibrated configurations containing N= 90 and N= 300 particles are presented. The droplets formed in this way are well isolated, with little to no exchange of particles between them. Furthermore, despite not being the true ground state, they still have a very large and negative energy, thus indicating that the relaxation times required for them to evolve into the true ground state droplet greatly exceed the time scales accessible in the simulation. These metastable configurations are therefore expected to be easily realized and observed in experiments, when initial tight confining traps are used to thermalize the system. Conclusion. We have observed self-bound dipolar droplets in a realistic quantum Monte Carlo simulation of a molecular BEC. In particular, the results suggest that single droplets and droplet assemblies, which are well known from magnetic quantum gases [1], exist in molecular BECs over a wide range of interactions parameters and strengths. In certain aspects, strongly dipolar systems can thus behave similarly to weakly dipolar systems in the mean-field limit, as well as to other systems with competing interactions [50,59–61]. However, the parameters and stability regimes for the formation of many-body states can be fundamentally different, as highlighted by the observation of droplets even for strongly attractive short range interactions. We note that, within the scope of the present investigation, the system’s behavior still shows a strong dependence on the value of the scattering length, despite being well outside the universality regime, so that other scattering parameters must also be relevant. More work along the lines of Refs. [44–47] is needed to explore this question in further detail. While the droplets realized in this work do not overlap, it will be interesting to investigate whether supersolid states also persist in the strongly interacting regime [10]. Certainly, the reduced number of particles per droplet makes PIGS method appropriate to study the possible BEC-supersolid transition, as it has been done recently for dysprosium system containing hundreds of atoms [62]. Furthermore, we expect that by suitably engineering the confinement, a cigar shaped molecular BEC could be turned into a self-organized stack of strongly interacting layers with strong intraand interlayer dipolar couplings [53], which is known to give rise to additional many-body phenomena [63–65]. Our work establishes quantum Monte Carlo simulations as a tool to realistically model current and future experiments with ultracold molecular gases. We anticipate that systematic comparisons of such simulations, mean-field models, and experiments will greatly advance our understanding of dipolar quantum matter in the near future. Note added.— During the preparation of this manuscript, we have become aware of the first experimental observation of droplets in a molecular BEC, including ones that are stable at negative scattering lengths [66]. We thank Jens Hertkorn, Thomas Pohl, Matteo Ciardi, and Sebastian Will for discussions, as well as Phillip Groß and Lukas Leczek for a careful reading of the manuscript. T.L. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 949431), from Carl Zeiss Foundation, and the Austrian Science Fund (FWF) 10.55776/COE1. R.B., J.S-B., F.M. and J.B. acknowledge financial support from Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 (Spain) under Grant No. PID2020-113565GB-C21 and from AGAUR-Generalitat de Catalunya Grant No. 2021SGR-01411. R.B. acknowledges funding from ADAGIO (Advanced ManufacturIng Research Fellowship Programme in the Basque – New Aquitaine Region) MSCA COFUND Post-Doctoral fellowship programme (Grant agreement No. 101034379). [1] F. Böttcher, J.-N. Schmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, Reports on Progress in Physics 84, 012403 (2020). [2] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. LaburtheTolra, B. L. Lev, and T. Pfau, Reports on Progress in Physics 86, 026401 (2022). [3] M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, Nature 539, 259 (2016). [4] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Phys. Rev. X6, 041039 (2016). [5] F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051 (2019). [6] L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019). [7] L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Phys. Rev. X 9, 021012 (2019). [8] Y.-C. Zhang, T. Pohl, and F. Maucher, Phys. Rev. A 104, 013310 (2021). [9] J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. Res. 3, 033125 6 (2021). [10] M. Schmidt, L. Lassablière, G. Quéméner, and T. Langen, Phys. Rev. Res. 4, 013235 (2022). [11] A. R. P. Lima and A. Pelster, Phys. Rev. A 84, 041604 (2011). [12] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Phys. Rev. Lett. 116, 215301 (2016). [13] F. Wächtler and L. Santos, Phys. Rev. A 93, 061603 (2016). [14] D. Baillie, R. M. Wilson, R. N. Bisset, and P. B. Blakie, Phys. Rev. A 94, 021602 (2016). [15] H. Saito, Journal of the Physical Society of Japan 85, 053001 (2016). [16] F. Böttcher, M. Wenzel, J.-N. Schmidt, M. Guo, T. Langen, I. Ferrier-Barbut, T. Pfau, R. Bombín, J. Sánchez-Baena, J. Boronat, and F. Mazzanti, Phys. Rev. Res. 1, 033088 (2019). [17] D. Petter, G. Natale, R. M. W. van Bijnen, A. Patscheider, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 122, 183401 (2019). [18] R. Bombín, V. Cikojević, F. Mazzanti, and J. Boronat, Phys. Rev. A 109, 033312 (2024). [19] T. Langen, G. Valtolina, D. Wang, and J. Ye, Nature Physics 20, 702 (2024). [20] L. De Marco, G. Valtolina, K. Matsuda, W. G. Tobias, J. P. Covey, and J. Ye, Science 363, 853 (2019). [21] A. Schindewolf, R. Bause, X.-Y. Chen, M. Duda, T. Karman, I. Bloch, and X.-Y. Luo, Nature 607, 677 (2022). [22] M. Duda, X.-Y. Chen, A. Schindewolf, R. Bause, J. von Milczewski, R. Schmidt, I. Bloch, and X.-Y. Luo, Nature Physics 19, 720 (2023). [23] J. Cao, H. Yang, Z. Su, X.-Y. Wang, J. Rui, B. Zhao, and J.-W. Pan, Phys. Rev. A 107, 013307 (2023). [24] N. Bigagli, C. Warner, W. Yuan, S. Zhang, I. Stevenson, T. Karman, and S. Will, Nature Physics 19, 1579 (2023). [25] J. J. Park, Y.-K. Lu, A. O. Jamison, and W. Ketterle, Nature Physics 19, 1567 (2023). [26] J. Lin, G. Chen, M. Jin, Z. Shi, F. Deng, W. Zhang, G. Quéméner, T. Shi, S. Yi, and D. Wang, Phys. Rev. X13, 031032 (2023). [27] N. Bigagli, W. Yuan, S. Zhang, B. Bulatovic, T. Karman, I. Stevenson, and S. Will, Nature (2024), 10.1038/s41586-024-07492-z. [28] W.-J. Jin, F. Deng, S. Yi, and T. Shi, “Bose-einstein condensates of microwave-shielded polar molecules,” (2024), arXiv:2406.06412 [cond-mat.quant-gas]. [29] S. Balibar, Nature 464, 176 (2010). [30] F. Deng, X.-Y. Chen, X.-Y. Luo, W. Zhang, S. Yi, and T. Shi, Phys. Rev. Lett. 130, 183001 (2023). [31] P. S. Żuchowski and J. M. Hutson, Phys. Rev. A 81, 060703 (2010). [32] R. Bause, A. Christianen, A. Schindewolf, I. Bloch, and X.-Y. Luo, The Journal of Physical Chemistry A 127, 729 (2023). [33] L. Lassablière and G. Quéméner, Phys. Rev. Lett. 121, 163402 (2018). [34] T. Karman and J. M. Hutson, Phys. Rev. Lett. 121, 163401 (2018). [35] T. Karman and J. M. Hutson, Phys. Rev. A 100, 052704 (2019). [36] T. Karman, Phys. Rev. A 101, 042702 (2020). [37] L. Anderegg, S. Burchesky, Y. Bao, S. S. Yu, T. Karman, E. Chae, K.-K. Ni, W. Ketterle, and J. M. Doyle, Science 373, 779 (2021). [38] K. Matsuda, L. De Marco, J.-R. Li, W. G. Tobias, G. Valtolina, G. Quéméner, and J. Ye, Science 370, 1324 (2020). [39] A. V. Gorshkov, P. Rabl, G. Pupillo, A. Micheli, P. Zoller, M. D. Lukin, and H. P. Büchler, Phys. Rev. Lett. 101, 073201 (2008). [40] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Reports on Progress in Physics 72, 126401 (2009). [41] S. B. Prasad, T. Bland, B. C. Mulkerin, N. G. Parker, and A. M. Martin, Phys. Rev. Lett. 122, 050401 (2019). [42] K. K. Voges, P. Gersema, M. Meyer zum Alten Borgloh, T. A. Schulze, T. Hartmann, A. Zenesini, and S. Ospelkaus, Phys. Rev. Lett. 125, 083401 (2020). [43] A. Macia, J. Sánchez-Baena, J. Boronat, and F. Mazzanti, Phys. Rev. Lett. 117, 205301 (2016). [44] R. Ołdziejewski and K. Jachymski, Phys. Rev. A 94, 063638 (2016). [45] S. Yi and L. You, Phys. Rev. A 63, 053607 (2001). [46] C. Ticknor, Phys. Rev. Lett. 100, 133202 (2008). [47] D.-W. Wang, New Journal of Physics 10, 053005 (2008). [48] A. Sarsa, K. Schmidt, and W. Magro, The Journal of Chemical Physics 113, 1366 (2000). [49] R. Rota, J. Casulleras, F. Mazzanti, and J. Boronat, Phys. Rev. E 81, 016707 (2010). [50] F. Cinti and M. Boninsegni, Phys. Rev. A 96, 013627 (2017). [51] Y. Kora and M. Boninsegni, Journal of Low Temperature Physics 197, 337 (2019). [52] M. Boninsegni, Results in Physics 31, 104935 (2021). [53] K. Mukherjee, M. N. Tengstrand, T. A. Cardinale, and S. M. Reimann, Phys. Rev. A 108, 023302 (2023). [54] W. Kirkby, T. Bland, F. Ferlaino, and R. N. Bisset, SciPost Phys. Core 6, 084 (2023). [55] S. Giorgini, J. Boronat, and J. Casulleras, Phys. Rev. A60, 5129 (1999). [56] F. Mazzanti, A. Polls, and A. Fabrocini, Phys. Rev. A 67, 063615 (2003). [57] A. Macia, F. Mazzanti, J. Boronat, and R. E. Zillich, Phys. Rev. A 84, 033625 (2011). [58] M. Wenzel, F. Böttcher, T. Langen, I. Ferrier-Barbut, and T. Pfau, Phys. Rev. A 96, 053630 (2017). [59] Y. H. Liu, L. Y. Chew, and M. Y. Yu, Phys. Rev. E 78, 066405 (2008). [60] G. Pupillo, A. Micheli, M. Boninsegni, I. Lesanovsky, and P. Zoller, Phys. Rev. Lett. 104, 223002 (2010). [61] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359, 301 (2018). [62] R. Bombín, F. Mazzanti, and J. Boronat, New Journal of Physics 26, 013052 (2024). [63] C. Trefzger, C. Menotti, and M. Lewenstein, Phys. Rev. Lett. 103, 035304 (2009). [64] C.-C. Huang and W.-C. Wu, Phys. Rev. A 82, 053612 (2010). [65] L. Du, P. Barral, M. Cantara, J. de Hond, Y.-K. Lu, 7 and W. Ketterle, Science 384, 546 (2024). [66] S. Zhang, N. Bigagli, W. Yuan, B. Bulatovic, H. Kwak, T. Karman, I. Stevenson, and S. Will, “Exploration of dipolar quantum phases with a degenerate gas of NaCs molecules.” Talk at DAMOP 2024. [67] S. Yi and L. You, Phys. Rev. A 61, 041604 (2000). [68] R. Schützhold, M. Uhlmann, Y. Xu, and U. R. Fischer, Int. J. Mod. Phys. B 20, 3555 (2006). [69] A. R. P. Lima and A. Pelster, Phys. Rev. A 86, 063609 (2012). [70] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, Phys. Rev. A 74, 013623 (2006). [71] F. Mazzanti, A. Polls, and A. Fabrocini, Phys. Rev. A 71, 033615 (2005). 8 Appendix: The extended Gross-Pitaevskii equation for the molecular potential The minimization of the energy functional with respect to the single-particle wave function within a mean field approximation yields the non-linear wave equation µψ(r) = −ℏ2∇2 2m+Zdr′VGPE(r−r′)|ψ(r′)|2+Hqu(r)ψ(r),(A.1) which determines the equilibrium state of the condensate for a given chemical potential µ. The pseudopotential VGPE(r)is given by VGPE(r) = gδ(r) + C3 r3(3 cos2θ−1) ,(A.2) where g=4πℏ2a mis the usual contact copuling constant, with athe s-wave scattering length. This pseudopotential is built such that its scattering properties computed in the first-order Born approximation reproduce those of the full potential of Eq. (1) [45,67]. On the other hand, the term Hqu(r) = 32 3gra3 πQ5−r0 3a|ψ(r)|3(A.3) accounts for the effect of quantum fluctuations, which provide a crucial stabilization mechanism for the otherwise collapsing system. The factor 3dividing the argument of the Q5function stems from our definition of the dipole length, which is three times larger than the usual definition add =mC3/(3ℏ2). Notice also that, compared to the beyond-mean-field correction for a standard dipole-dipole interaction, which we refer to as HDDI qu (r), the expression in Eq. (A.3) features an extra minus sign in the argument of the Q5 function to account for the reversed sign of the dipolar term of Eq. (1) [11,68,69]. Remarkably, this reversed sign significantly changes the real and imaginary contributions of the Q5term. We illustrate this in Fig. 4, where we show the real and imaginary parts of both Hqu and HDDI qu computed for the homogeneous system with a/r0= 0.11 (i.e. δr= 1.75). As it can be seen from the figure, the LHY correction for the reversed DDI shows a smaller real part, together with a significantly larger and unphysical imaginary part, which renders the application of Eq. (A.1) unreliable to describe the present system in the regime of parameters considered. In order to compute the critical numbers provided in the main text, we solve Eq. (A.1), and subsequently compute numerically the Hankel transform of the wave function [70,71] to take advantage of the cylindrical symmetry of the problem. The critical number is determined by obtaining the particle number for which the energy of the ground state configuration becomes negative, indicating the presence of a self-bound stable solution. 0 5 10 15 20 25 nr3 0 0 100 200 300 400 500 Hqu/E0 Re{Hqu/E0} Im{Hqu/E0} Re{HDDI qu /E0} Im{HDDI qu /E0} Figure 4. Real (dashed) and imaginary (dot-dashed) parts of the quantum fluctuations correction to the chemical potential Hqu for the potential of Eq. A.2 (blue) and the DDI pseudopotential (red) as a function of the homogeneous density n. The parameters are a/r0= 0.11 (δr= 1.75)