Formation and dynamics of self-bound droplets in dipolar molecular condensate

Abstract

We study self-bound quantum droplets in the regime dominated by microwave-induced non-axisymmetric dipole−dipole interactions, using the extended Gross−Pitaevskii equation with the Lee−Huang−Yang corrections. We identify the existence region through numerical simulations and employ an anisotropic Gaussian−super-Gaussian variational ansätz to capture the characteristic density profile of the droplets, with a Gaussian profile along the narrow x direction and super-Gaussian profiles in the extended (y, z) plane. Within this variational framework, we characterize the self-binding, spatial localization, and density-compression properties of the droplets and find good agreement between the variational predictions and the numerical results. Collisions between droplets moving along different directions reveal a strong directional dependence, with outcomes ranging from quasi-elastic rebound and merger to fragmentation. In addition, we explore the rotational dynamics of a single self-bound droplet about all three Cartesian axes, revealing rich and controllable three-dimensional rotational dynamics. Together, these results demonstrate how non-axisymmetric dipolar interactions provide versatile means for controlling the translational, collisional, and rotational dynamics of self-bound quantum droplets.

Graphical abstract

Keywords

quantum droplets / dipole−dipole interactions / BEC

1 Introduction

Ultracold dipolar molecules [1–4] are emerging as a promising testbed for many-body physics [5, 6], quantum simulations [7–9], and quantum data processing [10, 11]. A major challenge in these settings arises from severe collisional losses due to short-range molecular encounters [12–17]. Recent advances in microwave shielding have helped to significantly suppress these losses, enabling the creation of long-lived molecular Bose−Einstein condensates (BECs) [18–27].

Building on this progress, it has been demonstrated that dual-frequency microwave-field shielding [28–35] enables highly tunable strong long-range dipole−dipole interactions (DDIs), thus providing a control level beyond that achievable in BEC of magnetic atoms [36–42]. In the dual-field dressing scheme, the ellipticity of the \sigma field introduces an non-axisymmetric component of the DDI, proportional to \sim {\sin }^2\theta \cos 2\varphi, which breaks the axial symmetry and makes the interaction direction-dependent in the transverse plane [28, 43]. The linearly polarized \pi field provides an additional control channel for the usual axially symmetric component. By tuning the microwave amplitudes and detunings along a physically allowed dressing trajectory, the axially symmetric component can be compensated while a non-axisymmetric interaction remains [44].

The recent experimental realization [45, 46] of self-bound molecular quantum droplets (QDs) in the BEC of NaCs molecules with strong DDI [46] has put forward a new setup for the work with self-trapped states in the polar-molecule BEC [47], similar to QDs in atomic BECs [48–53]. The recent work by Baillie [54] established a comprehensive framework for microwave-dressed dipolar quantum droplets, systematically mapping their ground-state (GS) structures and properties across the two interaction parameters, {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)} and {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}, using Gaussian variational and extended Gross−Pitaevskii (GPE) approaches. Here, {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)} and {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} characterize the relative strengths of the axisymmetric and non-axisymmetric DDI, respectively. Building on this framework, we concentrate on the purely non-axisymmetric limit, {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}=0 and {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\ne 0, where the droplets acquire strongly anisotropic, flattened density profiles. In this regime, we refine the variational description by adopting an anisotropic Gaussian−super-Gaussian ansätz and extend the analysis beyond equilibrium properties to the dynamics of self-bound droplets, including their directional collision behavior and rotational motion about the three Cartesian axes.

The subsequent presentation is organized as follows. Section 2 introduces the theoretical model. Section 3 analyzes the ground-state properties using the Gaussian−super-Gaussian variational ansätz and numerical simulations. Section 4 addresses head-on collisions between identical QDs and the rotational response of a single QD. Section 5 summarizes the main results and discusses directions for future work.

2 The model

We consider BECs of electric dipolar bialkali molecules [45], which are dressed by two microwave fields [28, 55, 56]. These are a linearly polarized field oscillating along the Z-axis (the \pi-field), and an elliptically polarized one rotating in the (X,Y) plane (the \sigma-field), respectively, see a schematic illustration of the control fields in Fig. 1. The interaction between the dipole molecules is controlled by a strongly repulsive short-range hardcore contact interaction and the long-range DDIs with the potential [22, 28]

(1)U_{\mathrm{d}\mathrm{d}}(\mathit{R})={\displaystyle \frac{3G_s}{4\pi R^3}}[{ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}(1-3{\cos }^2\theta )+\sqrt{3}{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}{\sin }^2\theta \cos 2\varphi ],

where \theta is the angle between \mathit{R} and Z-axis, and \varphi is the angle between the projection of \mathit{R} onto the (X,Y) plane and the X-axis, m and \hslash being the mass of the molecule and Planck constant, respectively. Coefficient G_s=4\pi {\hslash }^2{a}_s/m is the contact-interaction strength, with a_s being the s-wave scattering length.

The first term of Eq. (1) represents the usual axially symmetric DDI between dipolar molecules, which can be tuned by the combined effects of the two control microwaves [45] with a relative strength {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}=a_{\mathrm{d}0}/a_s, where a_{\mathrm{d}0} is the effective DDI scattering length. The second term of Eq. (1) represents the non-axially symmetric DDI between the dipolar molecules, which is controlled by the ellipticity of the \sigma-field [28, 45] with a relative strength {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=a_{\mathrm{d}2}/a_s, where a_{\mathrm{d}2} is the respective DDI scattering length. Here, both {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)} and {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} can be tuned over a wide range of positive, zero and negative values.

The angular structure of Eq. (1) can be naturally understood in terms of the spherical harmonics Y_{lm}(\theta ,\varphi ). In particular, the first term is proportional to the spherical harmonic Y_{20}(\theta ,\varphi ), corresponding to the m=0 component. This term preserves rotational symmetry about the Z-axis and represents the conventional dipole−dipole interaction between polarized molecules. In contrast, the second term is proportional to the real linear combination Y_{2,2}(\theta ,\varphi )+Y_{2,-2}(\theta ,\varphi ), corresponding to the m=\pm 2 components. This term explicitly breaks the axial symmetry and introduces a four-fold anisotropy in the transverse (X,Y) plane.

The evolution of the system’s wave function \mathrm{\Psi } is governed by the extended GPE (alias eGPE) [54],

(2)\begin{array}{rl}\mathrm{i}\hslash \frac{\mathrm{\partial }}{\mathrm{\partial }T}\mathrm{\Psi }=& [-\frac{{\hslash }^2}{2m}{\mathrm{\nabla }}_{XYZ}^2+G_s{\left|\mathrm{\Psi }\right|}^2\\ & +\int U_{\mathrm{d}\mathrm{d}}(\mathit{R}-{\mathit{R}}^{\mathrm{\prime }}){\left|\mathrm{\Psi }({\mathit{R}}^{\mathrm{\prime }})\right|}^2\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}+{\mathrm{\Gamma }}_{\mathrm{Q}\mathrm{F}}{\left|\mathrm{\Psi }\right|}^3]\mathrm{\Psi }.\end{array}

Here, the effect of quantum fluctuations (QFs) is represented by the strength

(3){\mathrm{\Gamma }}_{\mathrm{Q}\mathrm{F}}={\displaystyle \frac{32G_s{a}_s^{3/2}}{3\sqrt{\pi }}}\text{Re}\{{\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)})\},

of the Lee−Huang−Yang (LHY) term [54], with

(4){\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)})={\displaystyle \int {\displaystyle \frac{\mathrm{d}{\mathrm{\Omega }}_k}{4\pi }}{\left[1+{\displaystyle \frac{{\tilde{U}}_{\mathrm{d}\mathrm{d}}(\mathit{k})}{G_s}}\right]}^{5/2}.}

Here, {\mathrm{\Omega }}_k is the solid angle in k-space, and {\tilde{U}}_{\mathrm{d}\mathrm{d}}(\mathit{k}) is the Fourier transform of U_{\mathrm{d}\mathrm{d}}(\mathit{R}), is given by

(5){\tilde{U}}_{\mathrm{d}\mathrm{d}}(\mathit{k})=G_s[{ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}\left({\displaystyle \frac{3k_Z^2}{k^2}}-1\right)-\sqrt{3}{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}{\displaystyle \frac{k_X^2-k_Y^2}{k^2}}].

In the numerical evaluation of this Fourier transform, a spherical cutoff of radius R_c is introduced to regularize the DDI integral, which gives rise to a cutoff function s(k{R}_c), the full momentum-space potential {\tilde{U}}_{\mathrm{d}\mathrm{d}}(\mathit{k}) being recovered in the limit of R_c\to \mathrm{\infty }. Further details of the calculation of {\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}) are provided in Appendix A.

The total particle number N is defined by wave function \mathrm{\Psi } as

(6)N={\displaystyle \int |\mathrm{\Psi }(\mathit{R}){|}^2\mathrm{d}\mathit{R}.}

By means of rescaling,

(7)R=l_{0}r,T=t_{0}t,\mathrm{\Psi }=l_0^{-3/2}\psi ,{\displaystyle \frac{m{l}_0^2}{\hslash t_0}}=1,

with \psi denoting the corresponding dimensionless wave function. To obtain self-bound stationary states, we consider solutions of the form

(8)\psi (\mathit{r},t)=\phi (\mathit{r}){\mathrm{e}}^{-\mathrm{i}\mu t},

where \mu is the chemical potential and \phi (\mathit{r}) is a time-independent wave function. Substituting Eqs. (7) and (8) into Eq. (2) yields the stationary eGPE,

(9)\begin{array}{rl}\mu \phi =& -\frac{1}{2}{\mathrm{\nabla }}_{xyz}^2\phi +g_s|\phi {|}^2\phi \\ & +\phi \int u_{\mathrm{d}\mathrm{d}}(\mathit{r}-{\mathit{r}}^{\prime })|\phi ({\mathit{r}}^{\prime }){|}^2\mathrm{d}{\mathit{r}}^{\prime }+{\gamma }_{\mathrm{Q}\mathrm{F}}|\phi {|}^3\phi .\end{array}

Here, the dimensionless strength of the contact interaction is g_s=4\pi a_s/l_0, while the LHY coefficient is {\gamma }_{\mathrm{Q}\mathrm{F}}=4{g_s}^{5/2}\text{Re}\{{\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)})\}/(3{\pi }^2). The dimensionless dipolar interaction potential is given by u_{\mathrm{d}\mathrm{d}}(\mathit{r})=3g_{\mathrm{s}}/(4\pi r^3)[{ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}(1-3{\cos }^2\theta )+\sqrt{3}{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}{\sin }^2\theta \cos 2\varphi ]. Obviously, the integral \int |\psi (\mathit{r}){|}^2\mathrm{d}\mathit{r}=N.

The energy functional associated with the stationary eGPE is [54]

(10)E=E_{\mathrm{p}}+E_{\mathrm{M}\mathrm{F}}+E_{\mathrm{L}\mathrm{H}\mathrm{Y}}+E_{\mathrm{d}\mathrm{d}},

where

(11)E_{\mathrm{p}}={\displaystyle \int \left[{\displaystyle \frac{1}{2}}|{\mathrm{\nabla }}_{xyz}\phi (\mathit{r}){|}^2\right]\mathrm{d}\mathit{r},}

(12)E_{\mathrm{M}\mathrm{F}}={\displaystyle \int \left[{\displaystyle \frac{g_s}{2}}|\phi (\mathit{r}){|}^4\right]\mathrm{d}\mathit{r},}

(13)E_{\mathrm{L}\mathrm{H}\mathrm{Y}}={\displaystyle \int \left[{\displaystyle \frac{2}{5}}{\gamma }_{\mathrm{Q}\mathrm{F}}|\phi (\mathit{r}){|}^5\right]\mathrm{d}\mathit{r},}

(14)E_{\mathrm{d}\mathrm{d}}={\displaystyle \frac{1}{2}}{\displaystyle \int \mathrm{d}\mathit{r}{\displaystyle \int \mathrm{d}{\mathit{r}}^{\mathrm{\prime }}{u}_{\mathrm{d}\mathrm{d}}(\mathit{r}-{\mathit{r}}^{\mathrm{\prime }})|\phi (\mathit{r}){|}^2|\phi ({\mathit{r}}^{\mathrm{\prime }}){|}^2.}}

To quantify the spatial localization of the QD, we define the effective volume [57]

(15)V={\displaystyle \frac{{\left({\displaystyle \int |\phi (\mathit{r}){|}^2\mathrm{d}\mathit{r}}\right)}^2}{{\displaystyle \int |\phi (\mathit{r}){|}^4\mathrm{d}\mathit{r}}}}.

In the following, all quantities are presented in dimensionless units. We choose l_0=2 μm, yielding the characteristic time scale t_0=m{l}_0^2/\hslash \simeq 9.83\mathrm{m}\mathrm{s}. Consequently, r=1 and t=1 correspond to physical scales of 2 μm and 9.83\mathrm{m}\mathrm{s}, respectively. The physical density is related to the dimensionless density through |\mathrm{\Psi }{|}^2=l_0^{-3}|\psi {|}^2, such that |\psi {|}^2=1 corresponds to the density of 1.25\times {10}^{11}{\mathrm{c}\mathrm{m}}^{-3}.

In what follows below we set {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}=0, to focus on effects of the non-axisymmetric DDIs. This condition can be experimentally achieved by tuning the amplitude (Rabi frequency) and the detuning of the \pi-field, while keeping the \sigma-field parameters fixed [45].

In Section 3, we first fix the effective scattering length a_s and vary the particle number N and the relative strength {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} of the non-axisymmetric DDI. In Section 4, we consider a second idealized parameter cut in which a_{\mathrm{d}2} and N are fixed while a_s is varied.

3 Ground-state (GS) properties

3.1 Formation of QDs

In this section, we adopt the representative value a_s=2100a_0, which is taken from the experimental data in Ref. [45], where a_0 is the Bohr radius, the corresponding value of the coefficient in Eq. (9) being g_s=4\pi a_s/l_0=0.7. The GS-QD state was produced by Eq. (9) by means of the imaginary-time method [58]. Stability of the QDs was tested by direct simulations of their perturbed evolution (not shown here in detail), confirming that all QDs are stable, as might be expected. Typical examples of the QDs for different values of {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} are shown in Figs. 2 (a1)–(a3), with N=10000 and {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=1.5, 1.9, and 4, respectively. As {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} increases, the QD becomes denser at intermediate interaction strengths and gradually develops a more anisotropic shape at larger {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}.

Figure 2(b) shows the existence region of self-bound QDs for {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}=0 and a_s=2100a_0. Stable self-bound states exist only in the orange region above the boundary (black line), whereas the grey region corresponds to unbound states where no bound-state solution is found. Self-bound states are not supported for sufficiently small particle numbers N or small values of {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}. In particular, the limit of {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=0 does not support self-bound states because the effective mean-field interaction becomes purely repulsive (g_s>0). At N=1000, the threshold found in Fig. 2(b) is consistent with the result reported by Baillie [54]. Notably, in Fig. 2(b) confirms the stability of the QDs against decay.

As shown in Fig. 2, under the condition {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}\equiv 0 and {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\ne 0, the droplets exhibit a highly anisotropic, elliptical pancake-like density profile, characterized by a narrow width along the x direction and an extended, nearly flat-top density distribution in (y,z) plane. In Ref. [54], an anisotropic Gaussian variational ansätz was employed to characterize the GS properties of microwave-dressed dipolar quantum droplets in the general case, where both {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)} nd {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} are nonzero. While this ansätz provides a useful description of the overall anisotropic structure, it does not fully capture the pronounced flat-top feature of the density profile in the purely non-axisymmetric regime considered here. To better account for this characteristic droplet shape, we introduce an anisotropic Gaussian–super-Gaussian variational ansätz:

(16){\phi }_{\mathrm{V}\mathrm{A}}(\mathit{r})=A{\mathrm{e}}^{-{\left(\frac{x}{w_x}\right)}^2}{\mathrm{e}}^{-\left[{\left(\frac{y}{w_y}\right)}^{{\alpha }_y}+{\left(\frac{z}{w_z}\right)}^{{\alpha }_z}\right]},

where w_x, w_y, and w_z are the variational widths (w_{x,y,z}>0), while {\alpha }_y and {\alpha }_z control the flatness of the density profile along the y and z directions, respectively, with {\alpha }_{y,z}\ge 2. The amplitude A is determined by the particle number N and the variational parameters.

Substituting Eq. (16) into the energy functional in Eq. (10), the total energy is expressed as a function of variational parameters (w_x,w_y,w_z,{\alpha }_y,{\alpha }_z). The kinetic, contact-interaction, and LHY contributions are evaluated analytically from the corresponding integrals of the trial state. The dipolar contribution is evaluated numerically using FFT convolution with the same spherical cutoff as in the imaginary-time calculations. The variational parameters are then determined by numerically minimizing the total energy using the \mathrm{f}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{h} routine [59].

In Fig. 3, the dashed lines represent the variational profiles. Overall, the variational results agree well with the imaginary-time profiles over the main density region, with a slight overestimation of the low-density tails. The variational results yield {\alpha }_y,{\alpha }_z well above 2, confirming the pronounced flat-top character of the density profiles in the (y,z) plane and supporting the validity of the Gaussian−super-Gaussian ansätz for describing the droplet structure.

For all considered interaction strengths, the QDs exhibit pronounced directional dependence, with the largest spatial extent along the y direction. As {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} increases, the density profiles first broaden and then become more localized, indicating a nonmonotonic evolution of the QD size. This observation motivates a quantitative analysis of the QD properties. We use chemical potential \mu and total energy E to characterize the cohesion of the self-bound state, and the effective volume V and peak density n_{\mathrm{p}} to quantify its spatial localization and density compression. We further compare these numerical results with the variational approximation (VA) introduced above to assess the capability of the Gaussian−super-Gaussian ansätz in describing the QD properties. Firstly, we examine how these quantities depend on {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} and N.

As shown in Figs. 4(a) and (b), the solid and dotted curves represent the numerical and VA results, respectively. For fixed particle numbers N=5000 and N=10000, both the chemical potential \mu and the total energy E exhibit a nonmonotonic trend with increasing {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}.

The two quantities first decrease and reach minima near {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\approx 3.2, indicating the strongest self-binding within the sampled range. For larger {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}, both quantities gradually increase, signaling weakened effective binding. The minima occur at nearly the same value for the two particle numbers shown here, suggesting only a weak N dependence in this parameter window; a denser scan would be needed to establish the precise dependence of the minimum position on N.

Besides the energetic properties, it is also instructive to examine how the density and size of the QD evolve. As shown in Fig. 4(c), the effective volume V first increases with {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}, reaches a maximum near {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\approx 1.9, and then decreases at larger anisotropy. A correlated behavior is observed in the peak density n_{\mathrm{p}} shown in Fig. 4(d), which first decreases and subsequently increases at stronger anisotropy, reflecting initial expansion of the QD followed by compression at large {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}.

To further characterize the basic properties of self-bound QDs, we next select a representative {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=2 and investigate their dependence on particle number N. As shown in Figs. 5(a) and (b), the chemical potential \mu decreases while the total energy E becomes increasingly negative with increasing N, reflecting enhanced self-binding. The variation is rapid at small particle number and gradually slows at around N=10000, approaching an approximately linear dependence within the explored range.

Importantly, the condition is satisfied throughout the explored parameter range, consistent with the Vakhitov−Kolokolov (VK) stability criterion for self-bound states [60]. This behavior suggests a crossover from a dilute regime at small particle numbers to a denser and more weakly compressible QD state at larger N.

At the same representative {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=2, we next examine how the density and size properties evolve with N. As shown in Fig. 5 (c), the peak density n_{\mathrm{p}} with N. Meanwhile, the effective volume V in Fig. 5(d) increases rapidly at small particle number and then grows more slowly as N keeps growing, suggesting that additional particles contribute to expansion of the QD, rather than increase of its density. This behavior is consistent with the flat-top regime of self-bound quantum droplets [48].

Having discussed the basic properties of the QD, we now focus on its geometric anisotropy induced by the non-axisymmetric dipole-dipole interaction. To characterize the anisotropic deformation, we define the anisotropy ratios

(17){\eta }_{zy}=W_z/W_y,

(18){\eta }_{xy}=W_x/W_y,

where the effective lengths along the three spatial directions are defined as

(19)W_x={\displaystyle \frac{{\left({\displaystyle \int |\phi (x,0,0){|}^2\mathrm{d}x}\right)}^2}{{\displaystyle \int |\phi (x,0,0){|}^4\mathrm{d}x}}},

(20)W_y={\displaystyle \frac{{\left({\displaystyle \int |\phi (0,y,0){|}^2\mathrm{d}y}\right)}^2}{{\displaystyle \int |\phi (0,y,0){|}^4\mathrm{d}y}}},

(21)W_z={\displaystyle \frac{{\left({\displaystyle \int |\phi (0,0,z){|}^2\mathrm{d}z}\right)}^2}{{\displaystyle \int |\phi (0,0,z){|}^4\mathrm{d}z}}}.

This definition is convenient for the present flat-top droplets: in the Thomas−Fermi limit, where the central density is approximately constant within the half-width R_i^{\mathrm{T}\mathrm{F}} and vanishes outside it, the above definition gives W_i=2R_i^{\mathrm{T}\mathrm{F}} [57]. Thus, W_i/W_j directly measures the corresponding Thomas−Fermi aspect ratio.

Here, W_y is chosen as the reference length since the QD generally has the largest spatial extent along the y direction. Therefore, smaller values of {\eta }_{zy} and {\eta }_{xy} indicate stronger anisotropy, corresponding to a more elongated density distribution along y.

Figures 6(a) and (b) show the anisotropy ratios as functions of {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} and N, respectively. Since increasing {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} strengthens the anisotropic dipolar interaction, while increasing N enhances the overall nonlinear interaction energy, the QD is expected to become increasingly elongated along the y direction. As a result, both {\eta }_{zy} and {\eta }_{xy} are generally expected to decrease with increasing {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} or N, reflecting enhanced geometric anisotropy. Overall, the numerical results follow this tendency.

However, a slight deviation from this trend appears near {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=1.9, where {\eta }_{zy} exhibits a small local maximum, corresponding to a relatively thicker density profile along the z direction. This behavior may originate from the competition between the anisotropic DDI and the stabilizing effects of the contact and LHY interactions, which can temporarily weaken the effective compression along the z direction.

Overall, the VA results agree well with the numerical results over the explored parameter ranges, demonstrating that the present Gaussian−super-Gaussian ansätz captures well the main properties of the QDs.

3.2 Contact interaction tuning at fixed a_d2

To further clarify the role of the contact interaction in the formation of self-bound quantum droplets, we consider the dependence of the droplet properties on the s-wave scattering length a_s while keeping the dipolar length a_{\mathrm{d}2} fixed. Although the independent experimental control of these two interaction parameters may require appropriate tuning of the microwave dressing fields, a reasonable schemes for their independent tunability have been discussed in Ref. [61]. Here, we therefore treat a_s and a_{\mathrm{d}2} as independently controllable parameters and focus on the theoretical role of the contact interaction by varying a_s at a fixed a_{\mathrm{d}2}. This provides a direct way to isolate how the short-range contact interaction affects the formation and properties of the self-bound droplets.

Specifically, we fix the non-axisymmetric dipolar length at a_{\mathrm{d}2}=4200a_0 and set N=2000, while varying a_s. This defines an idealized fixed-a_{\mathrm{d}2} path within the present eGPE model. For {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}=0, substituting G_s=4\pi {\hslash }^2{a}_s/m and {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=a_{\mathrm{d}2}/a_s into Eq. (1) gives

(22)\begin{array}{rl}{U}_{\mathrm{d}\mathrm{d}}(\mathit{R})=& \frac{3\sqrt{3}{G}_s{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}}{4\pi R^3}{\sin }^2\theta \cos 2\varphi \\ =& \frac{3\sqrt{3}{\hslash }^2{a}_{\mathrm{d}2}}{m{R}^3}{\sin }^2\theta \cos 2\varphi .\end{array}

Thus, the non-axisymmetric mean-field DDI kernel automatically unchanged along this path, while the contact-interaction coefficient G_s and the LHY coefficient {\gamma }_{\mathrm{Q}\mathrm{F}} vary with a_s.

Figure 7 summarizes the dependence of the QD properties on a_s. Panels (a1)–(a4) show the density distributions in the central x=0 plane, corresponding to the (y,z) plane, for a_s=2100a_0, 1600a_0, 800a_0, and 50a_0, respectively. Panel (b) shows the corresponding effective volume V as a function of a_s.

Over the parameter range shown, decreasing a_s produces a continuous compression of the QD, as indicated by the decreasing effective volume and increasing peak density. These trends indicate increasing compression and suggest a possible approach to collapse for a_s\to 0. The stationary states obtained by ITP for the parameters shown in Fig. 7 remain stable in direct simulations for a sufficiently long time. Note that, in this figure, as a_s becomes very small, the increasingly compressed QDs require progressively finer spatial resolution, and the numerical results become sensitive to the grid size; we therefore terminate the numerical scan the a_s at a small finite value.

4 Dynamics of QDs

4.1 Collision dynamics

To investigate the dynamical properties of self-bound QDs, we consider head-on collisions between two identical QDs with N=10000, {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=2, and a_s=2100a_0. The initial states are prepared by imprinting opposite momenta onto two spatially separated stationary QDs.

For a head-on collision along a coordinate axis j=x,y,z, the initial state is written once in the form [62]

(23)\begin{array}{rl}\psi (\mathit{r},t=0)=& {\phi }_1(\mathit{r}-{\mathit{r}}_0){\mathrm{e}}^{\mathrm{i}{\mathit{\zeta }}_j\cdot \mathit{r}}\\ & +{\phi }_2(\mathit{r}+{\mathit{r}}_0){\mathrm{e}}^{-\mathrm{i}{\mathit{\zeta }}_j\cdot \mathit{r}},\end{array}

where \mathit{r}=(x,y,z), {\mathit{r}}_0=r_{0j}{\hat{\mathit{e}}}_j is the initial displacement of one droplet, and {\mathit{\zeta }}_j={\zeta }_j{\hat{\mathit{e}}}_j is the corresponding kick vector. The parameters are r_{0x}=12 for the x collision, r_{0y}=128 for the y collision, and r_{0z}=64 for the z collision. The initial displacements are chosen to avoid significant density overlap at t=0.

For collisions along the x direction, the density cross sections in the (x,t)-plane, shown in Figs. 8(a)–(d), clearly indicate that the collision dynamics strongly depend on the kick strength, {\zeta }_x. For weak kicks with {\zeta }_x=0.39, shown in Fig. 8(a), the QDs approach each other and subsequently reverse their motion before substantial density overlap occurs, producing a quasi-elastic rebound with only weak deformation. For stronger kicks with {\zeta }_x=0.79, shown in Fig. 8(b), the QDs merge upon collision and subsequently form a single self-bound QD.

For a still larger kick with {\zeta }_x=1.87, shown in Fig. 8(c), the collision still produces a merged central droplet, but two weak emitted branches also appear and propagate away from the collision region. At the largest kick considered here, {\zeta }_x=3.34, shown in Fig. 8(d), the central merged component is almost suppressed after the collision, and the final density distribution is dominated by two outgoing emission branches.

The (y,t)- and (z,t)-plane density cross sections in Figs. 9(a1)–(a2) and Figs. 9(b1)–(b2) show that the collisions are likewise strongly inelastic in both transverse directions. For the weaker kicks {\zeta }_y=0.39 in the y direction and {\zeta }_z=0.20 in the z direction, shown in Figs. 9(a1) and (b1), respectively, the QDs deform and fragment before substantial overlap. For the stronger kicks {\zeta }_y=0.79 and {\zeta }_z=0.39, shown in Figs. 9(a2) and (b2), respectively, the QDs first overlap and then break into fragments. No long-lived self-bound remnant is observed in either direction.

4.2 Rotational dynamics

We investigate the rotational response of a single self-bound QD to a time-dependent anisotropic DDI. We obtain a stationary QD for N=5000, {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}=0, {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}=2, and a_s=2100a_0. The stationary state serves as the initial state for real-time evolution with a continuously rotating anisotropic DDI about one of the three Cartesian axes.

The rotation angle evolves as \beta (t)=\mathrm{\Omega }t, where \mathrm{\Omega } is the angular velocity. The time-dependent orientation of the anisotropic DDI is implemented by transforming the Fourier-space momentum vector \mathit{k} according to

(24){\mathit{k}}_{\beta }^{\{x,y,z\}}={\mathsf{T}}_{\{x,y,z\}}\left[\beta (t)\right]\mathit{k},

where {\mathsf{T}}_{\{x,y,z\}} denotes the rotation matrix about the corresponding Cartesian axis, and k_{\beta ,j}^{\{x,y,z\}} is the jth component of {\mathit{k}}_{\beta }^{\{x,y,z\}} with j=x,y,z.

The corresponding transformed components are

(25)\begin{array}{rl}{\mathit{k}}_{\beta }^x=& \left(k_{\beta ,x}^x,k_{\beta ,y}^x,k_{\beta ,z}^x\right)\\ =& (k_x,k_y\cos \beta -k_z\sin \beta ,k_y\sin \beta +k_z\cos \beta ),\end{array}

(26)\begin{array}{rl}{\mathit{k}}_{\beta }^y=& \left(k_{\beta ,x}^y,k_{\beta ,y}^y,k_{\beta ,z}^y\right)\\ =& (k_x\cos \beta -k_z\sin \beta ,k_y,k_x\sin \beta +k_z\cos \beta ),\end{array}

(27)\begin{array}{rl}{\mathit{k}}_{\beta }^z=& \left(k_{\beta ,x}^z,k_{\beta ,y}^z,k_{\beta ,z}^z\right)\\ =& (k_x\cos \beta -k_y\sin \beta ,k_x\sin \beta +k_y\cos \beta ,k_z).\end{array}

Substituting the transformed momentum components into Eq. (5) gives

(28)\begin{array}{rl}{\tilde{U}}_{\mathrm{d}\mathrm{d}}(\mathit{k},t)=& G_s[20{ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}\left(\frac{3{\left(k_{\beta ,z}^{\{x,y,z\}}\right)}^2}{k^2}-1\right)\\ & -\sqrt{3}{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\frac{{\left(k_{\beta ,x}^{\{x,y,z\}}\right)}^2-{\left(k_{\beta ,y}^{\{x,y,z\}}\right)}^2}{k^2}].\end{array}

The wave function is propagated using the same split-step Fourier method as in the collision simulations. The angular velocities are chosen as {\mathrm{\Omega }}_x={\mathrm{\Omega }}_z=2.5\pi \times {10}^{-3} and {\mathrm{\Omega }}_y=1.25\pi \times {10}^{-2}. To reduce boundary effects, the numerical box is enlarged in the yz, xz, and xy planes for rotations about the x, y, and z axes, respectively.

The resulting density evolution is shown in Fig. 10. The three rows correspond to rotations about the x, y, and z axes, respectively, while the columns show \beta =0^{\circ }, {45}^{\circ }, and {90}^{\circ }. The density isosurface is plotted at n=0.1n_{\mathrm{m}\mathrm{a}\mathrm{x}}. As \beta increases, the changing orientation of the anisotropic interaction induces a corresponding deformation of the self-bound QD, which remains localized throughout the evolution.

5 Conclusion

We have investigated self-bound quantum droplets (QDs) supported by non-axisymmetric microwave-dressed dipole−dipole interactions (DDIs) in a polar-molecule Bose−Einstein condensate by solving the extended Gross−Pitaevskii equation (eGPE) including the Lee−Huang−Yang (LHY) correction. We find that stable self-bound QDs can exist even when the conventional axially symmetric DDI is completely suppressed, and identify their existence region in the (N,{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}) parameter plane. To characterize their strongly anisotropic and flattened density profiles, we introduce an anisotropic Gaussian−super-Gaussian variational ansätz. The optimized variational parameters accurately capture the characteristic widths and flat-top features of the QDs, with only small deviations in the low-density tails. The variational density profiles are in good overall agreement with the numerical eGPE results.

We further characterize the ground-state properties of the QDs through the chemical potential \mu, total energy E, effective volume V, peak density n_{\mathrm{p}}, and anisotropy ratios as functions of the particle number N and the non-axisymmetric dipolar interaction strength {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}. The results reveal pronounced anisotropic deformation, nonmonotonic dependences of the droplet properties on {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}, and increasingly strong self-binding with increasing N. The corresponding variational predictions show good agreement with the numerical simulations, demonstrating that the Gaussian−super-Gaussian ansätz provides a useful description of the equilibrium properties of these highly anisotropic QDs.

We have also explored the dynamical properties of the QDs through head-on collisions between identical droplets along different directions and through their rotational dynamics under reorientation of the microwave-controlled DDI. The collision outcomes exhibit a pronounced directional dependence. Collisions along the x axis lead to quasi-elastic rebound and merger, whereas collisions along the y and z axes are destructive and result in fragmentation without a long-lived self-bound remnant. In addition, direct simulations of the rotational dynamics show that a self-bound QD can dynamically follow the reorientation of the non-axisymmetric DDI, undergoing orientation-dependent deformation while maintaining its spatial localization. These results demonstrate that non-axisymmetric microwave-dressed DDIs provide versatile control over the equilibrium structure, collisional behavior, and three-dimensional rotational dynamics of self-bound quantum droplets.

As an extension of the present work, it would be interesting to investigate collective excitations [63–65] of non-axisymmetric QDs, because the pronounced anisotropic deformation, found here, breaks rotational symmetry and may split and mix low-lying shape modes. It would also be important to clarify specific mechanisms which account for outcomes of the direction-dependent collisions, reported in this work. Actually, the different outcomes suggest action of direction-dependent coupling between translational motion, internal deformation, and fragmentation. A systematic mode-resolved analysis of the corresponding energy-transfer pathways would require a dedicated dynamical study (Refs. [66–68]), which is a subject for further work.

6 Appendix A: The derivation of {\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)})

The calculation of {\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}) is derived by Baillie in Ref. [54]. In this appendix, we briefly reproduce the derivation, to facilitate understanding of properties of this function.

By introducing a spherical cutoff of radius R_c, the Fourier transform of the dipole−dipole interaction given in Eq. (5) can be rewritten as

(A1)\begin{array}{c}\begin{array}{rl}{\tilde{U}}_{\mathrm{d}\mathrm{d}}^{(R_c)}(\mathit{k})=& {\tilde{U}}_{\mathrm{d}\mathrm{d}}(\mathit{k})s(k{R}_c)\\ =& G_s[{ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}\left(3{\cos }^2{\theta }_k-1\right)\\ & -\sqrt{3}{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}{\sin }^2{\theta }_k\cos 2{\phi }_k]s(k{R}_c),\end{array}\end{array}

where {\theta }_k and {\phi }_k denote the polar and azimuthal angles of the wave vector \mathit{k}, respectively, and the cutoff function s(k{R}_c) is

(A2)\begin{array}{c}s(k{R}_c)=1+3(k{R}_c{)}^{-2}\cos (k{R}_c)-3(k{R}_c{)}^{-3}\sin (k{R}_c).\end{array}

In the limit of R_c\to \mathrm{\infty }, one has s(k{R}_c)\to 1, so that {\tilde{U}}_{\mathrm{d}\mathrm{d}}^{R_c}(\mathit{k})\to {\tilde{U}}_{\mathrm{d}\mathrm{d}}(\mathit{k}).

To evaluate the angular average, the integrations over {\theta }_k can be carried out analytically, yielding

(A3)\begin{array}{c}\begin{array}{rl}{\mathcal{Q}}_5(x,0)=& \frac{5(x-1{)}^3}{16\sqrt{3x}}[\ln (\sqrt{1+2x}-\sqrt{3x})-\frac{\ln (1-x)}{2}]\\ & +\frac{1}{16}(11+4x+9x^2)\sqrt{1+2x}.\end{array}\end{array}

The remaining dependence on {\varphi }_k is represented by an auxiliary function,

(A4)\begin{array}{c}\begin{array}{rl}I({\varphi }_k)=& {\left(1-\frac{2}{\sqrt{3}}{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\cos 2{\varphi }_k\right)}^{5/2}\\ & \times {\mathcal{Q}}_5\left(\frac{\sqrt{3}{ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}+{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\cos 2{\varphi }_k}{\sqrt{3}-2{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}\cos 2{\varphi }_k},0\right).\end{array}\end{array}

The full coefficient {\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)}) is then obtained as

(A5)\begin{array}{c}\text{Re}\{{\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)})\}={\displaystyle \frac{2}{\pi }}{\displaystyle {\int }_0^{\pi /2}\mathrm{d}{\varphi }_k|\text{Re}\{I({\varphi }_k)\}|,}\end{array}

(A6)\begin{array}{c}\text{Im}\{{\mathcal{Q}}_5({ϵ}_{\mathrm{d}\mathrm{d}}^{(0)},{ϵ}_{\mathrm{d}\mathrm{d}}^{(2)})\}={\displaystyle \frac{2}{\pi }}{\displaystyle {\int }_0^{\pi /2}\mathrm{d}{\varphi }_k|\text{Im}\{I({\varphi }_k)\}|.}\end{array}

In parameter regions where {\mathcal{Q}}_5 is complex, only its real part is retained in the extended GPE. The so obtained values of Re[{\mathcal{Q}}_5] for different values of {ϵ}_{\mathrm{d}\mathrm{d}}^{(2)} and {ϵ}_{\mathrm{d}\mathrm{d}}^{(0)}=0 are collected in Table A1.

References

[1] R. V. Krems , Cold controlled chemistry, Phys. Chem. Chem. Phys. 10(28), 4079 (2008)

[2] M. Baranov , Theoretical progress in many-body physics with ultracold dipolar gases, Phys. Rep. 464(3), 71 (2008)

[3] K. K. Ni , S. Ospelkaus , M. H. G. De Miranda , A. Peer , B. Neyenhuis , J. J. Zirbel , S. Kotochigova , P. S. Julienne , D. S. Jin , and J. Ye , A high phase-space-density gas of polar molecules, Science 322(5899), 231 (2008)

[4] S. A. Moses , J. P. Covey , M. T. Miecnikowski , D. S. Jin , and J. Ye , New frontiers with quantum gases of polar molecules, Nat. Phys. 13(1), 13 (2017)

[5] T. Lahaye , C. Menotti , L. Santos , M. Lewenstein , and T. Pfau , The physics of dipolar bosonic quantum gases, Rep. Prog. Phys. 72(12), 126401 (2009)

[6] M. A. Baranov , M. Dalmonte , G. Pupillo , and P. Zoller , Condensed matter theory of dipolar quantum gases, Chem. Rev. 112(9), 5012 (2012)

[7] A. Micheli , G. K. Brennen , and P. Zoller , A toolbox for lattice spin models with polar molecules, Nat. Phys. 2(5), 341 (2006)

[8] L. D. Carr , D. DeMille , R. V. Krems , and J. Ye , Cold and ultracold molecules: Science, technology, and applications, New J. Phys. 11(5), 055049 (2009)

[9] E. Altman , K. R. Brown , G. Carleo , L. D. Carr , E. Demler , C. Chin , B. DeMarco , S. E. Economou , M. A. Eriksson , K. M. C. Fu , M. Greiner , K. R. A. Hazzard , R. G. Hulet , A. J. Kollár , B. L. Lev , M. D. Lukin , R. Ma , X. Mi , S. Misra , C. Monroe , K. Murch , Z. Nazario , K.K. Ni , A. C. Potter , P. Roushan , M. Saffman , M. Schleier-Smith , I. Siddiqi , R. Simmonds , M. Singh , I. B. Spielman , K. Temme , D. S. Weiss , J. Vučković , V. Vuletić , J. Ye , and M. Zwierlein , Quantum simulators: Architectures and opportunities, PRX Quantum 2(1), 017003 (2021)

[10] D. DeMille , Quantum computation with trapped polar molecules, Phys. Rev. Lett. 88(6), 067901 (2002)

[11] P. Rabl , D. DeMille , J. M. Doyle , M. D. Lukin , R. J. Schoelkopf , and P. Zoller , Hybrid quantum processors: molecular ensembles as quantum memory for solid state circuits, Phys. Rev. Lett. 97(3), 033003 (2006)

[12] K. K. Ni , S. Ospelkaus , D. Wang , G. Quemener , B. Neyenhuis , M. H. G. De Miranda , J. L. Bohn , J. Ye , and D. S. Jin , Dipolar collisions of polar molecules in the quantum regime, Nature 464(7293), 1324 (2010)

[13] G. Quéméner and P. S. Julienne , Ultracold molecules under control!, Chem. Rev. 112(9), 4949 (2012)

[14] J. L. Bohn , A. M. Rey , and J. Ye , Cold molecules: Progress in quantum engineering of chemistry and quantum matter, Science 357(6355), 1002 (2017)

[15] X. Ye , M. Guo , M. L. González-Martınez , G. Quéméner , and D. Wang , Collisions of ultracold ²³Na⁸⁷Rb molecules with controlled chemical reactivities, Sci. Adv. 4(1), eaaq0083 (2018)

[16] R. Bause , A. Christianen , A. Schindewolf , I. Bloch , and X. Y. Luo , Ultracold sticky collisions: Theoretical and experimental status, J. Phys. Chem. A 127(3), 729 (2023)

[17] A. V. Gorshkov , P. Rabl , G. Pupillo , A. Micheli , P. Zoller , M. D. Lukin , and H. P. Buchler , Suppression of inelastic collisions between polar molecules with a repulsive shield, Phys. Rev. Lett. 101(7), 073201 (2008)

[18] T. Karman and J. M. Hutson , Microwave shielding of ultracold polar molecules, Phys. Rev. Lett. 121(16), 163401 (2018)

[19] L. Lassablière and G. Quéméner , Controlling the scattering length of ultracold dipolar molecules, Phys. Rev. Lett. 121(16), 163402 (2018)

[20] L. Anderegg , S. Burchesky , Y. Bao , S. S. Yu , T. Karman , E. Chae , K. K. Ni , W. Ketterle , and J. M. Doyle , Observation of microwave shielding of ultracold molecules, Science 373(6556), 779 (2021)

[21] A. Schindewolf , R. Bause , X. Y. Chen , M. Duda , T. Karman , I. Bloch , and X. Y. Luo , Evaporation of microwaveshielded polar molecules to quantum degeneracy, Nature 607(7920), 677 (2022)

[22] X. Y. Chen , A. Schindewolf , S. Eppelt , R. Bause , M. Duda , S. Biswas , T. Karman , T. Hilker , I. Bloch , and X. Y. Luo , Field-linked resonances of polar molecules, Nature 614(7946), 59 (2023)

[23] J. Lin , G. Chen , M. Jin , Z. Shi , F. Deng , W. Zhang , G. Quéméner , T. Shi , S. Yi , and D. Wang , Microwave shielding of bosonic NaRb molecules, Phys. Rev. X 13(3), 031032 (2023)

[24] N. Bigagli , C. Warner , W. Yuan , S. Zhang , I. Stevenson , T. Karman , and S. Will , Collisionally stable gas of bosonic dipolar ground state molecules, Nat. Phys. 19(11), 1579 (2023)

[25] F. Deng , X. Y. Chen , X. Y. Luo , W. Zhang , S. Yi , and T. Shi , Effective potential and superfluidity of microwave-shielded polar molecules, Phys. Rev. Lett. 130(18), 183001 (2023)

[26] X. Y. Chen , S. Biswas , S. Eppelt , A. Schindewolf , F. Deng , T. Shi , S. Yi , T. A. Hilker , I. Bloch , and X. Y. Luo , Ultracold field-linked tetratomic molecules, Nature 626(7998), 283 (2024)

[27] J. Dutta , B. Mukherjee , and J. M. Hutson , Universality in the microwave shielding of ultracold polar molecules, Phys. Rev. Res. 7(2), 023164 (2025)

[28] T. Karman , N. Bigagli , W. Yuan , S. Zhang , I. Stevenson , and S. Will , Double microwave shielding, PRX Quantum 6(2), 020358 (2025)

[29] N. Bigagli , W. Yuan , S. Zhang , B. Bulatovic , T. Karman , I. Stevenson , and S. Will , Observation of Bose-Einstein condensation of dipolar molecules, Nature 631(8020), 289 (2024)

[30] W. Yuan , S. Zhang , N. Bigagli , H. Kwak , C. Warner , T. Karman , I. Stevenson , and S. Will , Extreme loss suppression in an ultracold molecular gas with widely tunable dipolar interactions, Science 393(6817), 1226 (2026)

[31] Z. Shi , Z. Huang , F. Deng , W. J. Jin , S. Yi , T. Shi , and D. Wang , Bose−Einstein condensate of ultracold sodiumrubidium molecules with tunable dipolar interactions, Nat. Phys. 22(9), 1467 (2026)

[32] F. Deng , X. Hu , W. J. Jin , S. Yi , and T. Shi , Two- and many-body physics of ultracold molecules dressed by dual microwave fields, Nat. Commun. 16(1), 11219 (2025)

[33] T. A. Cardinale , T. Bland , and S. M. Reimann , Exploring supersolids of single-microwave shielded molecules via exact and mean-field theories, Commun. Phys. 9(1), 191 (2026)

[34] W. J. Jin , F. Deng , S. Yi , and T. Shi , Bose−Einstein condensates of microwave-shielded polar molecules, Phys. Rev. Lett. 134(23), 233003 (2025)

[35] G. Li , X. Jiang , B. Liu , Z. Chen , B. A. Malomed , and Y. Li , Two-dimensional anisotropic vortex quantum droplets in dipolar Bose−Einstein condensates, Front. Phys. (Beijing) 19(2), 22202 (2024)

[36] I. Ferrier-Barbut , H. Kadau , M. Schmitt , M. Wenzel , and T. Pfau , Observation of quantum droplets in a strongly dipolar Bose gas, Phys. Rev. Lett. 116(21), 215301 (2016)

[37] M. Schmitt , M. Wenzel , F. Bottcher , I. Ferrier-Barbut , and T. Pfau , Self-bound droplets of a dilute magnetic quantum liquid, Nature 539(7628), 259 (2016)

[38] H. Kadau , M. Schmitt , M. Wenzel , C. Wink , T. Maier , I. Ferrier-Barbut , and T. Pfau , Observing the Rosensweig instability of a quantum ferrofluid, Nature 530(7589), 194 (2016)

[39] A. Boudjemâa , Fluctuations and quantum self-bound droplets in a dipolar Bose−Bose mixture, Phys. Rev. A 98(3), 033612 (2018)

[40] G. Semeghini , G. Ferioli , L. Masi , C. Mazzinghi , L. Wolswijk , F. Minardi , M. Modugno , G. Modugno , M. Inguscio , and M. Fattori , Self-bound quantum droplets of atomic mixtures in free space, Phys. Rev. Lett. 120(23), 235301 (2018)

[41] C. R. Cabrera , L. Tanzi , J. Sanz , B. Naylor , P. Thomas , P. Cheiney , and L. Tarruell , Quantum liquid droplets in a mixture of Bose−Einstein condensates, Science 359(6373), 301 (2018)

[42] L. Chomaz , I. Ferrier-Barbut , F. Ferlaino , B. Laburthe-Tolra , B. L. Lev , and T. Pfau , Dipolar physics: A review of experiments with magnetic quantum gases, Rep. Prog. Phys. 86(2), 026401 (2023)

[43] P. Xu and G. Chen , Effective potential and scattering length of shielded polar molecules, Phys. Rev. A 112(4), 043306 (2025)

[44] C. J. Ho , J. Dutta , B. Mukherjee , J. M. Hutson , and M. R. Tarbutt , Tuning interactions between static-field-shielded polar molecules with microwaves, Phys. Rev. Res. 8(2), 023087 (2026)

[45] S. Zhang , W. Yuan , N. Bigagli , H. Kwak , T. Karman , I. Stevenson , and S. Will , Observation of self-bound droplets of ultracold dipolar molecules, Nature 651(8106), 601 (2026)

[46] T. Langen , J. Boronat , J. Sánchez-Baena , R. Bombın , T. Karman , and F. Mazzanti , Dipolar droplets of strongly interacting molecules, Phys. Rev. Lett. 134(5), 053001 (2025)

[47] M. Schmidt , L. Lassabliere , G. Quéméner , and T. Langen , Self-bound dipolar droplets and supersolids in molecular Bose−Einstein condensates, Phys. Rev. Res. 4(1), 013235 (2022)

[48] D. S. Petrov , Quantum mechanical stabilization of a collapsing Bose-Bose mixture, Phys. Rev. Lett. 115(15), 155302 (2015)

[49] D. S. Petrov and G. E. Astrakharchik , Ultradilute low-dimensional liquids, Phys. Rev. Lett. 117(10), 100401 (2016)

[50] R. N. Bisset , R. M. Wilson , D. Baillie , and P. B. Blakie , Ground-state phase diagram of a dipolar condensate with quantum fluctuations, Phys. Rev. A 94(3), 033619 (2016)

[51] A. R. P. Lima and A. Pelster , Quantum fluctuations in dipolar Bose gases, Phys. Rev. A 84(4), 041604 (2011)

[52] I. Ferrier-Barbut , M. Schmitt , M. Wenzel , H. Kadau , and T. Pfau , Liquid quantum droplets of ultracold magnetic atoms, J. Phys. At. Mol. Opt. Phys. 49(21), 214004 (2016)

[53] M. Guo and T. Pfau , A new state of matter of quantum droplets, Front. Phys. (Beijing) 16(3), 32202 (2021)

[54] D. Baillie , Symmetry and self-bound droplets in dipolar molecular gases, Phys. Rev. Res. 8(2), 023219 (2026)

[55] W. Zhang , H. Liu , F. Deng , K. Chen , S. Yi , and T. Shi , Supersolid phases in ultracold gases of microwave shielded polar molecules, (2025)

[56] A. Schindewolf , J. Hertkorn , I. Stevenson , M. Ciardi , P. Groß , D. Wang , T. Karman , G. Quéméner , S. Will , T. Pohl , and T. Langen , Colloquium: Strongly dipolar molecular Bose−Einstein condensates: From few- to many-body physics, Rev. Mod. Phys. 98(3), 031002 (2026)

[57] Z. Zhao , G. Li , Z. Chen , H. B. Luo , B. Liu , B. A. Malomed , and Y. Li , Bulk modulus of three-dimensional quantum droplets, Phys. Rev. A 113(4), 043315 (2026)

[58] J. Yang , Nonlinear Waves in Integrable and Non-Integrable Systems, SIAM, Philadelphia, (2010)

[59]

[60] N. G. Vakhitov and A. A. Kolokolov , Stationary solutions of the wave equation in a medium with nonlinearity saturation, Radiophys. Quantum Electron. 16(7), 783 (1973)

[61] S. Singh and C. H. Greene , Efimov effect in ultracold microwave-shielded polar molecules, Phys. Rev. Lett. 136(23), 233402 (2026)

[62] Y. Hu , Y. Fei , X. L. Chen , and Y. Zhang , Collisional dynamics of symmetric two-dimensional quantum droplets, Front. Phys. (Beijing) 17(6), 61505 (2022)

[63] F. Wächtler and L. Santos , Ground-state properties and elementary excitations of quantum droplets in dipolar Bose−Einstein condensates, Phys. Rev. A 94(4), 043618 (2016)

[64] L. Chomaz , R. M. W. Van Bijnen , D. Petter , G. Faraoni , S. Baier , J. H. Becher , M. J. Mark , F. Wächtler , L. Santos , and F. Ferlaino , Observation of roton mode population in a dipolar quantum gas, Nat. Phys. 14(5), 442 (2018)

[65] I. Ferrier-Barbut , M. Wenzel , F. Bottcher , T. Langen , M. Isoard , S. Stringari , and T. Pfau , Scissors mode of dipolar quantum droplets of dysprosium atoms, Phys. Rev. Lett. 120(16), 160402 (2018)

[66] D. Baillie , R. M. Wilson , R. N. Bisset , and P. B. Blakie , Self-bound dipolar droplet: A localized matter-wave in free space, Phys. Rev. A 94, 021602(R) (2016)

[67] G. Ferioli , G. Semeghini , L. Masi , G. Giusti , G. Modugno , M. Inguscio , A. Gallemı , A. Recati , and M. Fattori , Collisions of self-bound quantum droplets, Phys. Rev. Lett. 122(9), 090401 (2019)

[68] S. K. Adhikari , Statics and dynamics of a self-bound matter-wave quantum ball, Phys. Rev. A 95(2), 023606 (2017)

Rights & permissions

Higher Education Press