A new form of liquid matter: Quantum droplets

Abstract

This brief review summarizes recent theoretical and experimental results which predict and establish the existence of quantum droplets (QDs), i.e., robust two- and three-dimensional (2D and 3D) self-trapped states in Bose–Einstein condensates (BECs), which are stabilized by effective self-repulsion induced by quantum fluctuations around the mean-field (MF) states [alias the Lee–Huang–Yang (LHY) effect]. The basic models are presented, taking special care of the dimension crossover, 2D → 3D. Recently reported experimental results, which exhibit stable 3D and quasi-2D QDs in binary BECs, with the inter-component attraction slightly exceeding the MF self-repulsion in each component, and in single-component condensates of atoms carrying permanent magnetic moments, are presented in some detail. The summary of theoretical results is focused, chiefly, on 3D and quasi-2D QDs with embedded vorticity, as the possibility to stabilize such states is a remarkable prediction. Stable vortex states are presented both for QDs in free space, and for singular but physically relevant 2D modes pulled to the center by the inverse-square potential, with the quantum collapse suppressed by the LHY effect.

Graphical abstract

Keywords

quantum droplet / Bose–Einstein condensate / Lee–Huang–Yang correction / votex state

1 Introduction

The theoretical and experimental work with multidimensional (two- and three-dimensional, 2D and 3D) solitons, i.e., self-trapped modes originating from the balance between nonlinear self-attraction of wave fields and their self-expansion driven, at the linear level, by diffraction and dispersion, is more difficult in comparison to the commonly known concept of 1D solitons [1–4]. On the one hand, 2D and 3D solitons offer many options inaccessible in 1D – in particular, the creation of 2D and 3D self-trapped states with the topological charge, that represents intrinsic vorticity [5], and more sophisticated 3D states, such as monopoles [6], skyrmions or hopfions with two independent topological charges [7–9], and knots [10]. On the other hand, identification of physically relevant models in which multidimensional solitons, both structureless fundamental ones and higher-order states featuring topological structures, are stable, and, eventually, creation of such objects in the experiment, are challenging, because the most common cubic self-attractive nonlinearity is only able to build completely unstable solitons in 2D and 3D, on the contrary to the ubiquitous one-dimensional nonlinear-Schrödinger (NLS) solitons, which are extremely robust objects, realizing the ground state (GS) of the respective model [11]. The fundamental cause of the instability is the fact that 2D and 3D NLS equations with the cubic self-attraction give rise, respectively, to the critical and supercritical collapse, i.e., a trend to develop singular solutions through catastrophic self-compression of the input [12]. Therefore, an issue of profound interest is elaboration of physically relevant settings which admit stabilization of 2D and 3D self-trapped localized modes, fundamental and topologically structured ones alike [1–4].

Among several approaches to resolving this issue, an arguably most successful one was theoretically proposed [13] and experimentally realized [14–26] quite recently. It relies upon the effect of quantum fluctuations as a correction to the mean-field (MF) dynamics of Bose–Einstein condensates (BECs), which was first predicted long ago by Lee, Huang, and Yang (LHY) [35], and was proposed to be used as a mechanism for the stabilization of 3D self-trapped states in Ref. [13]. In this context, soliton-like states, which are called quantum droplets (QDs), form due to the MF attraction, which may be provided either by contact (local) inter-component attraction in binary BEC [13–18], or by long-range dipole-dipole interactions (DDIs) in a single-component condensate of atoms carrying permanent magnetic moments [19–26]. The collapse of the droplets, which would take place in the MF approximation, is arrested by the LHY effect, which is effectively represented by local quartic self-repulsive terms in the respective NLS equations [that are usually called Gross–Pitaevskii equations (GPEs), in the application to BEC]. The competition between the MF attraction and LHY repulsion maintains a superfluid state whose density (taking very low values) cannot exceed a certain maximum, thus making it incompressible. This is a reason why this quantum macroscopic state is identified as a fluid, and localized states filled by it are called “droplets". In addition to the 3D QDs, droplets in the effectively 2D setting, which may exist under the action of tight confinement in one direction, imposed by an external potential, were theoretically predicted [42] and experimentally created [14, 15] too.

The objective of this article to provide a brief review of basic theoretical and experimental results on the theme of QDs. First, basic models, in the form of GPEs with the LHY corrections, are introduced in Section 2, in the full 3D form. Reductions of the equations to the 2D and 1D cases are also presented.

Next, Sections 3 and 4 provide a relatively detailed account of recent experimental findings. These are stable QDs in a binary BEC composed of two different hyperfine states of the same atomic species, as well as in a heteronuclear mixture (Section 3), and in single-component condensates of magnetic-dipolar atoms (Section 4). In Section 3, we also briefly address collisions between 3D droplets in the binary condensate, which were experimentally studied very recently [17].

Sections 5 and 6 report theoretical results. Because the current theoretical literature on QDs is vast, while the size of this article is limited, in these two sections we chiefly focus on most recent theoretical predictions of stable three- and two-dimensional QDs with embedded vorticity. Although they have not yet been reported in the experiment, vortical QDs promise to realize a variety of novel features. Basic results for stable “swirling" (vortical) 3D and 2D droplets, with unitary and multiple topological charges, are summarized in Section 5. Also included are results for vortex QDs in a semi-discrete 2D system [36]. Another type of confined vortex modes is considered in Section 6: effectively two-dimensional ones in the binary BEC pulled to the center by the inverse-square potential, see Eq. (30) below. In the framework of the MF theory, all stationary states in this setting are destroyed by the quantum collapse. However, the LHY terms suppresses the collapse and help one to create an otherwise missing GS, as well as stable vortex states [37]. These solutions are singular at r\to 0, but, nevertheless, physically relevant ones, as their norm (the number of atoms in the condensate) converges.

Section 7 concludes the article. In that section, we briefly mention related topics which are not considered in detail in this review, and discuss directions for the further development of studies in this area.

2 Theoretical models of quantum droplets

2.1 Models of QDs in three, two, and one dimensions

The energy density of a condensed Bose–Bose mixture in the MF approximation is

\begin{array}{}\text{(1)}& {\epsilon }_{\text{MF}}=\frac{1}{2}{g}_{11}{n}_1^2+\frac{1}{2}{g}_{22}{n}_2^2+g_{12}{n}_1{n}_2\mathrm{,}\end{array}

where g_{11}, g_{22}, and g_{12} are, respectively, the intra- and inter-species coupling constants characterizing the interaction between atoms, and n_j is the density of the j-th component of the mixture. In the case when intra-species interactions are repulsive, g_{\mathrm{11,22}}>0, the mixture is miscible, in the framework of the MF theory, if the intra-species repulsion dominates over the inter-species interaction, \sqrt{g_{11}{g}_{22}}>\left|g_{12}\right| [38]. Otherwise, the system is immiscible (at g_{12}>\sqrt{g_{11}{g}_{22}}), or it collapses (at ), if the attraction between the two species, accounted for by , is stronger than the effective single-species repulsion.

The celebrated LHY correction to the MF density (1), which is the leading term in the beyond-MF energy density, originating from the zero-point energy of the Bogoliubov excitations around the MF state [35], takes the following form, as derived by Petrov [13]:

\begin{array}{}\text{(2)}& {\epsilon }_{\text{LHY}}=\frac{128}{30\sqrt{\pi }}g{n}^2\sqrt{n{a}_s^3}\mathrm{,}\end{array}

where a_s is the s-wave scattering length, n is the density of both components, assuming that they are equal, g=4\pi {\hslash }^{2}m/a_s is the corresponding coupling constant, and m is the atomic mass. The LHY term makes sense only for a_s>0, i.e., repulsive intra-species interactions.

In this article, we concentrate on the consideration of symmetric modes in the binary BEC, with equal components of the pseudo-spinor wave function. Asymmetric states were considered too, in 1D [39], 2D [40], and 3D [41] cases alike. The asymmetry essentially affects stability of various modes. In particular, it tends to strongly destabilize vortex modes with unequal topological charges in the two components [40, 41]. Further, QDs in heteronuclear binary BEC [18] are always strongly asymmetric, due to their nature.

In a dilute condensate, the LHY term, \propto n^{5/2} in Eq. (2) is, generally, much smaller than the MF ones, \propto n^2 in Eq. (1), hence the LHY correction is negligible. However, when the binary condensate is close to the equilibrium point, at which the MF self-repulsion in each component is nearly balanced by the attraction between the components, the LHY term becomes essential or even dominant. The result is the spontaneous formation of robust QDs, due to equilibrium between the effective residual MF attraction (assuming, as said above, equal wave functions of the two components) and the LHY repulsion, which stabilizes the droplets against collapsing [13].

In this vein, by defining

\begin{array}{}\text{(3)}& \delta g\equiv g_{12}+\sqrt{g_{11}{g}_{22}}\mathrm{,}\end{array}

for g_{\mathrm{11,22}}>0 and , the latter condition corresponding to the inter-species attraction, Ref. [13] addressed the regime with

The resulting LHY-amended GPE for the wave function of both components in the symmetric 3D system can be written as

\begin{array}{}\text{(5)}& \text{i}{\partial }_{\tilde{t}}\varphi =\left(-\frac{1}{2}{\nabla }_{\tilde{\mathit{r}}}^2-3|\varphi {|}^2+\frac{5}{2}|\varphi {|}^3-\tilde{\mu }\right)\varphi \mathrm{,}\end{array}

where \tilde{\mathit{r}}, \tilde{t}, \tilde{\mu } are, respectively, rescaled coordinates, time, and chemical potential, as defined in Ref. [13], and the quartic self-repulsive term, (5/2)|\varphi {|}^3, corresponds to the LHY energy density given by Eq. (2). The suppression of the collapse in Eq. (5) is guaranteed by the fact that, for large values of the local density, n, the quartic self-repulsion dominates over the cubic self-attraction, that drives the onset of the collapse in the MF theory.

Equation (5) generates a family of stationary 3D-isotropic QDs, which are looked for as

\begin{array}{}\text{(6)}& \varphi \left(\tilde{\mathit{r}}\mathrm{,}\tilde{t}\right)={\varphi }_0\left(\tilde{r}\right)\mathrm{,}\end{array}

where \tilde{r} is the radial coordinate. Radial profiles of the 3D droplets, obtained as solutions of Eq. (5), as well as the respective energy per particle,

\begin{array}{}\text{(7)}& \tilde{E}/\tilde{N}=\frac{1}{2}{\left[{\displaystyle {\int }_0^{\infty }}{\varphi }_0^2(\tilde{r}){\tilde{r}}^2\text{d}\tilde{r}\right]}^{-1}{\displaystyle {\int }_0^{\infty }}\epsilon (\tilde{r}){\tilde{r}}^2\text{d}\tilde{r}\end{array}

(factor 1/2 takes into account the contribution of both mutually symmetric components into the total number of atoms), the particle emission threshold -\tilde{\mu }, and the spectrum of frequencies {\tilde{\omega }}_l of the droplet’s surface modes, with angular-momentum quantum number l, are displayed in Fig. 1. One can find that, in a wide range of {(\tilde{N}-{\tilde{N}}_c)}^{1/4}, where {\tilde{N}}_c\approx 18.65 is the stability boundary (the QDs are unstable at ), all excitation modes cross the threshold for sufficient small \tilde{N}, which means that the modes, excited by initial perturbations on the surface of the 3D droplet, are depleted by emission of small-amplitude waves, in terms of Eq. (5) [13].

Beyond-MF effects in lower dimensions, 2D and 1D, have also drawn much interest. The respective models were derived in Ref. [42], starting from the full 3D setting and including tight transverse confinement, imposed by an external potential acting in one direction, to induce the 3D \to 2D reduction, or in two directions, to impose reduction 3D \to 1D. In the 2D case, the energy density of the symmetric pseudo-spinor condensate, with equal densities of the “top" and “bottom" spinor components, n_{\uparrow \uparrow }=n_{\downarrow \downarrow }\equiv n, and equal scattering lengths, a_{\uparrow \uparrow }=a_{\downarrow \downarrow }\equiv a, was derived in the form of

\begin{array}{}\text{(8)}& {\epsilon }_{\text{2D}}=\frac{8\pi n^2}{{{\displaystyle \ln }}^2(a_{\uparrow \downarrow }/a)}\left[\ln (n/{\left(n_0\right)}_{\text{2D}})-1\right]\mathrm{,}\end{array}

where {\left(n_0\right)}_{\text{2D}}={\left(2\pi \right)}^1\exp \left(-2\gamma -3/2\right){\left(a{a}_{\uparrow \downarrow }\right)}^{-1}\ln (a_{\uparrow \downarrow }/a) is the equilibrium density of each component (\gamma \approx 0.5772 is the Euler’s constant). The corresponding LHY-amended GPE for the common wave function of both components reads

\begin{array}{}\text{(9)}& \text{i}{\partial }_t\psi =-\frac{1}{2}{\nabla }^2\psi \text{}+\text{}\frac{8\pi }{{{\displaystyle \ln }}^2(a_{\uparrow \downarrow }/a)}\ln \left(\frac{|\psi {|}^2}{\sqrt{e}{\left(n_0\right)}_{\text{2D}}}\right)|\psi {|}^2\psi \mathrm{.}\end{array}

The increase of the local density from small to large values leads to the change of the sign of the logarithmic factor in Eq. (8). As a result, the cubic term is self-focusing at small densities, initiating the spontaneous formation of QDs, and defocusing at large densities, thus arresting the transition to the collapse, and securing the stability of 2D QDs.

As shown in Fig. 2(a), the density dependence of the energy per particle, calculated in the framework of the corresponding many-body theory by means of the diffusion Monte Carlo (DMC) simulations [43], converges toward the analytical result given by Eq. (8), with the decrease of 1/\ln (a_{\uparrow \downarrow }/a).

For the 1D setting, the analysis performed in Ref. [42] had yielded the following effective energy density:

\begin{array}{}\text{(10)}& {\epsilon }_{\text{1D}}=\delta g\cdot n^2-4\sqrt{2}{(gn)}^{3/2}/(3\pi )\mathrm{,}\end{array}

the corresponding equilibrium density being {\left(n_0\right)}_{\text{1D}}=8g^3/(9{\pi }^2\delta g^2). The respective LHY-amended GPE features a combination of the usual MF cubic nonlinearity and a quadratic term, representing the LHY corrections in the 1D geometry:

\begin{array}{}\text{(11)}& \text{i}{\partial }_t\psi =-\frac{1}{2}{\partial }_{xx}\psi \text{}+\text{}\delta g\cdot |\psi {|}^2\psi \text{}-\text{}\frac{\sqrt{2}}{\pi }{g}^{3/2}|\psi |\psi \mathrm{.}\end{array}

Note that, in Eq. (11), LHY-induced quadratic term is self-focusing, on the contrary to the defocusing sign of the quartic term in the 3D equation (5). Because the most interesting results for QDs are obtained in the case of the competition between the residual MF term and its LHY-induced correction [44], in the 1D case the relevant situation is one with \delta g>0, when the residual MF self-interaction is repulsive, in contrast with the residual self-attraction adopted in the 3D setting, as mentioned above.

The energy per particle (10) is plotted, versus the rescaled particle number, and for different values of \delta g/g, in Fig. 2(b). Similar to the 2D case, for sufficiently small \delta g/g, the scaled data produced by the DMC method for the many-body theory, is in a good agreement with the Bogoliubov approximation given by Eq. (10).

2.2 Dimensional crossover for quantum droplets

The quasi-2D and quasi-1D description of QDs outlined above is valid for extremely strong transverse confinement. An estimate for experimentally relevant parameters [14–16] yields a respective estimate for the confinement size L_z\ll l_{\text{healing}}\sim 30 nm, where l_{\text{healing}} is the healing length in the condensate [37], while the actual value of L_z used in the experiment is \sim 0.6 \text{μ}m. For this reason, the dimension crossover \text{3D}\to \text{2D} requires a more careful consideration. In particular, for a loosely confined (“thick") quasi-2D layer of the condensate it may be relevant to consider the 2D version of Eq. (5), keeping the quartic LHY term, which is, strictly speaking, relevant in the full 3D space [37].

A detailed consideration of the dimensional crossover was presented in Ref. [45], which, similar to Ref. [42], addressed the binary BEC with repulsive intra-species interactions, g_{11}\approx g_{22}>0, and inter-species attraction, . In addition, periodic boundary conditions were imposed in the vertical direction. As above, the system is tuned to be close to the balance condition, defined as per Eqs. (3) and (4). In this case, the effective quasi-2D energy density representing the LHY effect is defined by the integration of the original 3D expression,

\begin{array}{}\text{(12)}& e_{\text{LHY}}^{\text{(2D)}}(\xi )=\underset{r\to 0}{{\displaystyle \lim }}\frac{\partial }{\partial r}\left[\frac{1}{2}r{\displaystyle \sum _{q_z}}\text{}\text{}{\displaystyle \int }{\text{d}}^2{\mathit{q}}_{\perp }{\text{e}}^{\text{i}\mathit{q}\mathit{r}}({\epsilon }_{bmq}\text{}-\text{}{A}_{\mathit{q}})\right]\mathrm{,}\end{array}

where {\epsilon }_{\mathit{q}}=\sqrt{q^4+2\xi q^2}, A_{\mathit{q}}=q^2+\xi, the summation is performed with respect to discrete wavenumbers in the vertical direction, while {\mathit{q}}_{\perp } is the 2D wave vector in the horizontal plane. A crucially important parameter which appears here is

\begin{array}{}\text{(13)}& \xi =(g_{11}{n}_{11}+g_{22}{n}_{22})/{\epsilon }_0\mathrm{,}\end{array}

with {\epsilon }_0=[{\hslash }^2/(2m)](2\pi /L_z{)}^2. It determines the ratio of the MF energy to the transverse-confinement energy.

For small \xi \ll 1, which is implied by the tight transverse confinement, an approximate calculation of the energy density (12) yields

\begin{array}{}\text{(14)}& e_{\text{LHY}}^{\text{(2D)}}(\xi )=\frac{\pi }{4}{\xi }^2\left[\text{ln}\xi +\text{ln}(2{\pi }^2)+\frac{1}{2}+\frac{{\pi }^2\xi }{3}\right]\mathrm{.}\end{array}

For the \text{3D}\to \text{1D} crossover, one can define the effective LHY energy density similarly, cf. Eq. (12):

\begin{array}{}\text{(15)}& e_{\text{LHY}}^{\text{(1D)}}(\xi )=\underset{r\to 0}{{\displaystyle \lim }}\frac{\partial }{\partial r}\left[\frac{1}{2}r{\displaystyle \sum _{q_x{q}_y}}{\displaystyle \int }\text{d}{q}_z{\text{e}}^{\text{i}\mathit{q}\mathit{r}}({\epsilon }_{\mathit{q}}-A_{\mathit{q}})\right]\mathrm{.}\end{array}

For small ξ, an approximate result is

\begin{array}{}\text{(16)}& e_{\text{LHY}}^{\text{(1D)}}(\xi )=-\frac{2\sqrt{2}}{3}{\xi }^{3/2}+c_2{\xi }^2+c_3{\xi }^3\mathrm{,}\end{array}

where c_2\simeq 3.06 and c_3\simeq 3.55.

The dependence of ratios

\begin{array}{}\text{(17)}& s^{\text{(2D,1D)}}=e_{\text{LHY}}^{\text{(2D,1D)}}/e_{\text{LHY}}^{\text{(3D)}}\end{array}

on ξ, where the 3D LHY energy density is e_{\text{LHY}}^{\text{(3D)}}=16\sqrt{2}\pi {\xi }^{5/2}/15, are demonstrated in Fig. 3. For small values of ξ, the approximate expression matches the numerically exact one well, while for large ξ, the ratios naturally approach 1.

Another approach to calculating the beyond-MF corrections at the dimensional crossover, which is based on the pioneering work [46], was elaborated in Ref. [47]. In that work, a one-component weakly interacting Bose gas satisfying the diluteness condition, \sqrt{n{a}_s^3}\ll 1, is assumed to be confined in one or two directions by a box potential with length l_{\perp } and periodic boundary conditions. The beyond-MF corrections to the energy density with the box potential, denoted by \mathrm{\Delta }{E}_{\text{2D,1D}}, and, in addition, \mathrm{\Delta }{E}_{{}_{\text{1D}}}^h for the \text{3D}\to \text{1D} confinement imposed by the harmonic-oscillator (HO) confinement (rather than by the box), which represent the dimensionality reduction, are displayed, as functions of

\begin{array}{}\text{(18)}& \kappa \equiv n{a}_s{l}_{\perp }\end{array}

in Fig. 4.

3 Experimental observations of two-component quantum droplets (QDs)

3.1 Oblate (quasi-two-dimensional droplets)

The creation of stable QDs in a BEC mixture of two Zeeman states of ³⁹K atoms, namely, |\uparrow \rangle =|F\mathrm{,}{m}_F\rangle =|\mathrm{1,}-1\rangle and |\downarrow \rangle =|\mathrm{1,0}\rangle, where F is the total angular momentum and m_F is its projection was reported in Ref. [14]. The potassium mixture is characterized by intra- and inter-species scattering length a_{\uparrow \uparrow }, a_{\downarrow \downarrow }, and a_{\uparrow \downarrow }. The residual MF interaction is proportional to the effective scattering length

\begin{array}{}\text{(19)}& \delta a=a_{\uparrow \downarrow }+\sqrt{a_{\uparrow \uparrow }{a}_{\downarrow \downarrow }}\end{array}

[cf. Eq. (3)], which identifies the boundary between repulsive (\delta a>0) and attractive () regimes. The interaction strengths can be tuned, via the Feshbach resonance (FR) [48], by an external magnetic field, B, as shown in Fig. 5(a). The MF energy of the mixture is proportional to \delta a, while the LHY correction scales with the intra-species scattering lengths a_{11}, a_{22}. Atoms creating the QD are loaded in a plane of a vertical blue-detuned lattice potential to compensate for gravity and a vertical red-detuned optical dipole trap, which provides a horizontal radial confinement. The experiment started with a sufficiently large magnetic field, B\approx 57.3 G, which, via the Feshbach resonance, corresponds to \delta a\approx 7a_0, where a_0 is the Bohr radius. In this case, the state of the binary BEC superfluid is miscible. Numbers of atoms in each component can be measured by means of the Stern–Gerlach separation in the course of free expansion of the gas, after the trapping potential was switched off. Subsequently, the magnetic field is ramped down, to drive the mixture into the attractive regime with , in which the radial confinement is simultaneously switched off, letting the atoms move freely in the horizontal plane. Figure 5(b) shows the evolution of typical images at different moments of time t, following the removal of the horizontal radial confinement, but keeping the vertical lattice potential. In the case of \delta a>0, the mixture features a gas-like expansion under an overall repulsive MF interaction, while for the attractive regime (with ), in which quantum fluctuations, i.e., the LHY effect, start to dominate. As a result, a stable self-trapped two-component droplet was observed in Ref. [14]. On the other hand, the same Fig. 5(b) demonstrates collapse occurring in a single-component condensate in the attractive regime, dramatically different from the behavior of the binary condensate.

Changes of radial size {\sigma }_r and peak density n_0 of the experimentally created QDs with the variation of the number of atoms, N, are shown in panel (c) of Fig. 5. In the attractive regime, both {\sigma }_r and n_0 remain approximately constant at large N, as expected for a liquid state. The existence of the mixture droplets require a minimum number of atom, N_c, below which a liquid-to-gas transition takes place, and the atomic cloud expands. Further, Fig. 5(d) shows strong dependence of the onset of the liquid-to-gas transition on magnetic field B. The critical number N_c increases with magnetic field B, corresponding to the attenuation of the effective attraction in the mixture.

A more accurate investigation of the liquid-to-gas transition of the binary BEC was reported in Ref. [15], by means of a similar experiment in the mixture of two different atomic states in the potassium condensate. The results also confirm the existence of QDs in the oblate (quasi-2D) configuration. Moreover, it was found that traditional matter-wave bright-soliton states, filled by the gaseous phase, and QDs, filled by the ultradilute superfluid, coexist in a bistable regime, providing an insight into the relation between these two kinds of self-trapped states.

3.2 Three-dimensional (isotropic) droplets

Following the original proposal by Petrov [13], QDs in the full 3D space were created in a weakly interacting binary condensate of ³⁹K [16]. The mixture is composed of two hyperfine states of ³⁹K. A cross dipole potential created by three red-detuned laser beams, and an optical levitating potential were employed in the experiment. The set of perpendicular beams was used to prepare the condensate, while the later element helped to make the residual confinement in all directions negligible. The residual scattering length of the mixture, \delta a=a_{12}+\sqrt{a_{11}{a}_{22}} [cf. Eq. (19)], decreases with the external magnetic field, vanishing at B_c=56.85 G. As Fig. 6(a) shows, at , i.e., in the case of , the binary condensate may be either a QD or an LHY gas, the two phases being separated by a critical number of atoms, N_c. When released from the external dipole trap in the attractive regime, the condensate with N>N_c keeps a constant size, i.e., it demonstrates a well-defined QD, while, in the case of it expands like a gas, as might be expected. In Fig. 6(a), the average size of the atomic cloud is \sigma =\sqrt{{\sigma }_x{\sigma }_z}, where {\sigma }_{x\mathrm{,}z} are half-widths of the density profile in the x and z directions, at the level of 1/\sqrt{e}.

The evolution of the condensate’s size σ in the QD phase at B=56.54G, as well as the total atom number, N, and the population ratio of two atomic states in the mixture, N_1/N_2, are presented in Fig. 6(b). As seen in middle panel, N rapidly drops during a few milliseconds, because of losses induced by three-body (3B) inelastic collisions, and eventually attains the critical values, N_c, at t_c=7 ms, where the liquid-to-gas phase transition takes place. The size, \sigma (t), exhibits a nearly constant value within the time interval 2 ms t_c, as seen in the top panel, which confirms the establishment of a QD. Afterwards, it expands as gas. In the course of the liquid-to-gas transition, the total number of atoms, N, and the population ratio, N_1/N_2, remain constants.

The measurements for σ, N, and N_1/N_2 at critical point N_c were extended to different values of magnetic field B, as shown in Fig. 6(c). The droplet’s size σ and critical atom number N_c increase with the increase of the magnetic field. The colored area in top panel corresponds to the theoretical prediction for σ in the range of norms N_c\le N\le 2N_c. The anisotropy measure, {\sigma }_x/{\sigma }_z-1, remains zero for different magnitudes of B, revealing that the droplet is a spherical isotropic one. The dependence of the critical atom number, N_c, on B, shows good agreement with values theoretically in Ref. [13] for the metastable and stable (dashed and solid lines, respectively) self-trapped states solution.

3.3 Collisions between quantum droplets

Collision of moving classical droplets may lead to their merger into a single one, provided that the surface tension is sufficient to absorb the kinetic energy of the colliding pair. Otherwise, the colliding droplets separate into two or more ones after the collision [49]. Similar phenomena in collisions of two component QDs were recently experimentally demonstrated in Ref. [17]. That work exhibits two different outcomes of the collision, i.e., merger and separation (passage).

For given magnetic field, there exists a critical velocity v_c, such that the colliding QDs merge at , and separate at v>v_c. Typical examples of the evolution of the colliding droplets are displayed in Fig. 7. In the case of , as shown in panel (a), the distance between the droplets decreases and finally stays being equal to zero, which implies the merger, as depicted in panel (b). By contrast, if the case of v>v_c, the kinetic energy of the moving QDs overcomes the surface tension, driving the separation after the collision, as shown in panel (d). As a consequence, the distance between the separating droplets increases, see panel (e). Further, it is shown in panels (c) and (f) that, in both cases, the total atom number decreases due to the strong 3B loss in the system [16].

Figure 7(h) presents a summary of results of experimentally observed collision in the plane of the rescaled atom number \tilde{N} and velocity \tilde{v}, as produced in Ref. [17]. The critical velocity {\tilde{v}}_c, which is the boundary between the merger (red diamonds) and passage (blue squares), exhibits different dependences on the number of atoms at small and large \tilde{N}, due to different energy scales dominating in these cases. In the regime of incompressibility at large \tilde{N}, the surface energy dominates, while the bulk and gradient energies may be negligible. Therefore, v_c is proportional to {\tilde{N}}^{-1/6}. In the opposite case of small \tilde{N}, the bulk energy has to be taken into account, because it cannot be separated from the surface energy. In this case, the consideration of the energy balance yields {\tilde{v}}_c\propto \sqrt{2\left|{\tilde{E}}_{\text{drop}}\right|/\tilde{N}}. The corresponding numerical results for the droplet collision, both with- and without the 3B loss, is displayed in Fig. 7(i). It is seen that the numerical simulation with 3B loss is in good agreement with the experimental results.

To probe the timescale of the collisions for various \tilde{N}, the velocity of moving droplets, v, is set to be slightly larger than {\tilde{v}}_c, to ensure that the separation takes place after the collision. The dependence of the timescale of the collision, \tilde{\tau }, on \tilde{N} is displayed in panel (j), revealing that, in the liquid-like regime at large \tilde{N}, the separation corresponds to longer timescales, because in this limit the colliding pair forms a single cloud in an excited state for a certain time interval, and they separate afterwards. Results of the corresponding numerical simulation, without the 3B loss, are shown in insets of Fig. 7(j).

3.4 Droplets in a heteronuclear bosonic mixture

As outlined above, QDs were first created in mixtures of two different spin states of ³⁹K atoms. The mechanism stabilizing two-component QDs applies as well to mixtures of different atomic species. Experimentally, this possibility was realized in Refs. [18] and [50], using a binary condensate of ⁴¹K and ⁸⁷Rb atoms. The results confirm the existence of stable droplets in the regime of relatively strong inter-species attraction, and expansion of the mixture in the case when the attraction is too weak. In particular, stable QDs were observed with the ratio of atom numbers N_{\text{K}}/N_{\text{Rb}}\approx 0.8, which is consistent with the results predicted by means of the analysis based on Ref. [13]. The heteronuclear droplets were demonstrated to have a lifetime \sim 10 ms, much longer than ones created in the binary condensate of ³⁹K, that, as mentioned above, was determined by 3B losses. The substantially longer lifetime offers one an opportunity to gain insight in intrinsic properties of the QDs, such as the observation of self-evaporation.

4 Single-component QDs in dipolar condensates

The theoretical and experimental findings summarized above demonstrate the possibility of the creation of stable droplets in binary BECs, based on the competition of the cubic nearly-balanced attraction between the two components and self-repulsion in each of them, and the additional quartic LHY-induced self-repulsion, see Eq. (5). Still earlier experiments had produced robust QDs in single-component dipolar BECs made of dysprosium [19–25] and erbium [26] atoms. Generating droplets in this setting is possible with the attraction provided by the long-range dipole-dipole interaction (DDI), and the stabilizing repulsion induced by the contact interaction, including the LHY term. The dipolar BEC are characterized by the scattering length a_s of the contact interaction, and the effective DDI length, a_{dd}. Accordingly, the interplay between the DDI and the contact interactions is controlled by parameter

\begin{array}{}\text{(20)}& {\epsilon }_{dd}=a_{dd}/a_s\mathrm{.}\end{array}

4.1 Quantum droplets in the condensate of dysprosium

In the experiments reported in Refs. [19, 20], isotope ¹⁶⁴Dy with dipolar length a_{dd}\simeq 131a_0, where a₀ is the Bohr radius, is employed to create BEC. The background scattering length of the contact interactions is a_{\text{bg}}=92a_0, which was modulated by many FRs. When {\epsilon }_{dd} is close to 1 [see Eq. ( 20)], the MF contact interactions and DDI nearly balance each other, making the contribution from the beyond-MF LHY effect crucially important for the creation of QDs.

In Ref. [20], a stable BEC containing ～10⁴ dysprosium atoms was created, by tuning the magnetic field to B_{\text{BEC}}\sim 6.962 G. The condensate was loaded into a radially symmetric, pancake-shaped trap with HO frequencies ({\nu }_x\mathrm{,}{\nu }_y\mathrm{,}{\nu }_z)=(\mathrm{46,44,133}) Hz, in the presence of the external magnetic field in the z direction, along which atomic magnetic dipoles are polarized, as shown in Fig. 8(a). Subsequently, the magnetic field was ramped down to a value at which a_s\approx a_{\text{bg}}, resulting in an angular roton instability. Thus, the condensate evolved into a set of N_d droplets, ranging between 2 and 10, which arranged themselves into a triangular structure. As shown in Fig. 8(b), N_d shows a linear dependence on the total number of atoms number, with N/N_d\approx 1750. The droplets strongly repel each other, maintaining distance d\approx 3.0 or 3.3 \text{μ}m for N_d=2 and N_d>2, respectively.

The spatial density distribution of the dysprosium condensate is characterized by its Fourier transform, S(k), which features a local maximum at k=2\pi /d\approx 2.5 \text{μ}m, where k=\sqrt{k_x^2+k_y^2}. The spectral weight,

\begin{array}{}\text{(21)}& \text{SW}={\displaystyle \sum _{k=1.5\text{μ}{\text{m}}^{-1}}^{5\text{μ}{\text{m}}^{-1}}}S(k)\mathrm{,}\end{array}

accounts for the strength of the structured states. It is subject to normalization S{W}_{\text{BEC}}=1 for the entire condensate. To explore properties of the spectral weight, a set of experimental data was collected, as shown in Fig. 8(d). The dysprosium condensate was generated close to the FR at a_s\approx a_{dd}, and then ramped down to a target value of the magnetic field, 6.860 G, with a constant speed, see the red arrow in Fig. 8(c). This was followed by a waiting stage, lasting for 20 ms. Then, the magnetic field was increased back to the high value at which the BEC was originally created, see the green arrow in Fig. 8(c). In the course of the experiment, atomic samples were imaged in situ, and the corresponding spectral weights were calculated as per Eq. (21), see Fig. 8(d). A well-defined hysteresis was thus observed, comparing the stages of the reduction of the magnetic field and return back to the original value, indicating that the system features bistability in the transition region.

To reveal the nature of the droplets observed in Ref. [20], they were trapped in a waveguide [19], which imposes a prolate cigar-like shape with aspect ratio \lambda \simeq 8. The magnetic field was ramped from B_{\text{BEC}}=6.962 G to B_1=6.656 G during 1 ms. When quenching up {\epsilon }_{dd}, the system, instead of collapsing, forms a metastable state composed of droplets whose number is 1\le N_d\le 6. For the case of N_d\ge 2, an average separation between the droplets was measured to be d=2.5\text{μ}\text{m}. The lifetime of these droplets is on order of hundreds of milliseconds, which is much larger than that of the two-component QDs supported by contact interactions. In addition, following the quench of the magnetic field, the expanding droplets overlap, as their size become comparable to or larger than the distance between them, which leads to the appearance of interference fringes. The observation of the fringes indicates that the individual droplets are phase-coherent objects, which was also observed in one-dimensional droplet arrays [27]. In this connection, it is relevant to mention that global phase coherence has been demonstrated in the supersolid state of matter, which was recently realized in several experiments [28–34].

4.2 Quantum droplets in the condensate of erbium

Erbium is another atomic species with permanent magnetic moment, which is appropriate for the realization of dipolar BEC. Isotope ¹⁶⁶Er was employed in Ref. [26] to investigate the BEC-QD crossover. The value of the respective background scattering length a_{\text{bg}} is comparable to the effective dipolar length, a_{dd}=65.5a_0, which makes it easy to realize condition {\epsilon }_{dd}\simeq 1, see Eq. (20). These atoms also feature a convenient set of FRs at ultra-low magnetic field values. The dependence of the scattering length a_s on B is demonstrated in Fig. 9(A). As shown in the upper inset, one can easily enter the range of {\epsilon }_{dd}>1, in which QDs may be observed.

The condensate was prepared at B=1.9 G, corresponding to a_s=81a_0. Following the evaporative cooling procedure, the magnetic field was reduced to 0.8 G [corresponding to a_s=67a_0, via the FR], with the atomic magnetic moments polarized along the weak-trapping axis. Finally, B was ramped down to a target value in the course of time t_r, which is followed by wait time t_h. Then, an absorption image of the gas was taken, after time-of-flight (TOF) t_{\text{TOF}}. Figures 9(b–d) display typical absorption images of the density profiles for t_r=10 ms (quenching), t_h=6 ms, t_{\text{TOF}}=27 ms, and different values of a_s. In particular, in the case of {\epsilon }_{dd}>1, as shown in Fig. 9(e), the density distribution is close to that predicted by the Thomas–Fermi (TF) approximation, which neglects the kinetic-energy term in GPE. The distribution of thermal atoms, see dotted lines in Fig. 9(e), is different from that in the central core, and remains mainly unaffected by the change of a_s. Collective oscillations of the coherent gas cloud is intimately related to the origin of the stabilization mechanism. In this work, the axial mode, which is the lowest-lying excitation in the system above the dipolar mode, was experimentally studied for both adiabatic and nonadiabatic ramps of the magnetic field. For both ramps, the results highlight a qualitative agreement with the theoretical predictions including the LHY term, revealing the fact that the LHY correction plays an essential role in stabilizing the system.

As said above, quantum fluctuations are expected to stabilize the system and help forming the droplets. However, 3B losses favor lower densities. The interplay between quantum fluctuations and 3B losses in the BEC-to-QD crossover is of great interest. Numbers of atoms in both the central-core (N_{\text{core}}) and thermal (N_{\text{th}}) components are shown, as a function of a_s, in Figs. 10(a) and (b), respectively, following the action of the nonadiabatic and adiabatic ramps. Both cases show a similar evolution. When the magnetic field is ramped down, the number of atoms in the central core remains constant for {\epsilon }_{dd}>1, then drops dramatically around {\epsilon }_{dd}\simeq 1, and finally curves up at lower a_s. In contrast to that, N_{\text{th}} shows weak dependence on a_s, confirming a picture in which dynamics of the thermal and condensed components are uncoupled. The observed evolution of N_{\text{core}} matches well with the theoretical calculation including the LHY correction (solid lines), but deviates from the one performed in the absence of the LHY term (dashed lines). The time evolution of N_{\text{core}} for various a_s in the droplet regime is displayed in Fig. 10(c), where N_{\text{core}} shows fast decay in the interval of ms, indicating that atoms are ejected from the high-density core through 3B losses. The steepness of this fast decay critically depends on a_s. The mean in situ density \overline{n} of the high-density component in the BEC-QD crossover is extracted with the help of the general 3B-loss relation, see further details in Ref. [26]. As shown in Fig. 10(d), one can see that the mean density \overline{n} attains a maximum at the threshold, with a_s\simeq a_{dd}, showing a quantitative agreement with numerical simulations including the LHY correction.

As a self-trapped state, the QD is expected to demonstrate its characteristic in the expansion regime. Typical examples of the TOF evolution of width {\sigma }_x of the high-density core are shown in Fig. 10(e). It is observed that, in the case of {\epsilon }_{dd}>1, the atomic cloud exhibits clear slowing-down of the expansion dynamics. The expansion velocity v_x is extracted by fitting the data to {\sigma }_x(t_{\text{TOF}})=\sqrt{{\sigma }_{x\mathrm{,0}}^2+v_x^2{t}_{\text{TOF}}^2}. The dependence of velocity v_x of the expansion of the high-density core on a_s is demonstrated in Fig. 10(f). In the droplet regime v_x gets a minimum at about 56a_0 ({\epsilon }_{dd}\sim 1.17), and grows when a_s gets far away from this minimum point. This behavior cannot be explained by the MF theory. On the other hand, simulations with the LHY correction reproduce the results produced by experimental measurements, see the solid line in Fig. 10(f).

5 Theoretical results: Stable quantum droplets with embedded vorticity

5.1 Three-dimensional vortex rings

Quantum droplets observed in experiments outlined above are fundamental modes, which do not carry any vorticity. It is natural to expect that vortex (alias spinning) modes may offer an opportunity to study more sophisticated properties of the QD state of matter [5]. Thus far, all studies of QDs with embedded vorticity were performed solely in the theoretical form. In particular, in Ref. [5] it was demonstrated that QD solutions with embedded vorticity exist in the model of the single-component dipolar condensate, but they all are unstable, hence physically irrelevant.

On the other hand, it was found that models of binary condensates with contact interactions, based on systems of LHY-amended GPEs readily give rise to stable vortex states. In particular, Ref. [41] addressed the 3D system for the two-component wave function {\psi }_{\mathrm{1,2}}. In the scaled form, the system takes the form of

\begin{array}{}\text{(22)}& \begin{array}{l}\text{i}\frac{\partial {\psi }_1}{\partial t}=-\frac{1}{2}{\nabla }^2{\psi }_1\text{}+\text{}({\left|{\psi }_1\right|}^2\text{}+\text{}{g}_{\text{LHY}}{\left|{\psi }_1\right|}^3){\psi }_1\text{}-\text{}g{\left|{\psi }_2\right|}^2{\psi }_1\mathrm{,}\\ \text{i}\frac{\partial {\psi }_2}{\partial t}=-\frac{1}{2}{\nabla }^2{\psi }_2\text{}+\text{}({\left|{\psi }_2\right|}^2\text{}+\text{}{g}_{\text{LHY}}{\left|{\psi }_2\right|}^3){\psi }_2\text{}-\text{}g{\left|{\psi }_1\right|}^2{\psi }_2\mathrm{,}\end{array}\end{array}

where the strength of the cubic self-repulsion in each component is scaled to be 1, while g>0 is the relative strength of the inter-component attraction. The LHY repulsion is characterized by coefficient g_{\text{LHY}}\simeq (128/3)\sqrt{2/\pi a_s^{3/2}}, where a_s is the intra-component scattering length. The corresponding stationary solutions for vortex droplets with chemical potentials {\mu }_{\mathrm{1,2}} and integer topological charges m_{\mathrm{1,2}} of the components are looked for, in cylindrical coordinates (\rho \mathrm{,}\theta \mathrm{,}z), as

\begin{array}{}\text{(23)}& {\psi }_{\mathrm{1,2}}=u_{\mathrm{1,2}}(\rho \mathrm{,}z)\exp (\text{i}{m}_{\mathrm{1,2}}\theta -\text{i}{\mu }_{\mathrm{1,2}}t)\end{array}

with real stationary wave functions u_{\mathrm{1,2}} obeying equations

\begin{array}{}\text{(24)}& \begin{array}{l}\mu u_{\mathrm{1,2}}+\frac{1}{2}\left(\frac{{\partial }^2}{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial }{\partial \rho }+\frac{{\partial }^2}{\partial z^2}-\frac{m_{\mathrm{1,2}}^2}{{\rho }^2}\right)u_{\mathrm{1,2}}\\ -\left(u_{\mathrm{1,2}}^2+g_{\text{LHY}}{u}_{\mathrm{1,2}}^3\right)u_{\mathrm{1,2}}+g{u}_{\mathrm{2,1}}^2{u}_{\mathrm{1,2}}=0.\end{array}\end{array}

Stability regions for the 3D QDs with embedded vorticity, i.e., vortex rings, which are symmetric with respect to the two components, with {\mu }_1={\mu }_2 and m_1=m_2=1, are shown in Fig. 11, both in the (\mu \mathrm{,}g) and (N\mathrm{,}g) planes. As demonstrated in Fig. 11(a), the vortex rings are stable in region , where {\mu }_{\text{co}} is the cutoff value of the chemical potential, below which no droplet can be found. Actually, {\mu }_{\text{co}} corresponds to the indefinitely broad QDs with a flat-top shape and diverging integral norm (number of atoms), the existence of {\mu }_{\text{co}} being a general property of self-trapped states in models with competing nonlinearities, such as the well-studied cubic-quintic combination [52, 53]. The stability interval of μ expands as g increases, due to the fact that the value of {\mu }_{\text{co}} decreases faster than the stability boundary {\mu }_{\text{st}}. The stability domain in the (N\mathrm{,}g) plane is shown in Fig. 11(b). Value N_{\min } indicates the minimum number of atoms necessary for the formation of a spinning droplet. The 3D vortex droplets are stable at N>N_{\text{st}}, while In the region of they exist but are unstable. The stability-boundary value N_{\text{st}} rapidly increases with the decrease of g, and diverges at g=g_{\min }\approx 1.3, below which system (22) does not maintain stable spinning states.

For the symmetric state with double vorticity, of m_1=m_2=2, a narrow stability region was obtained (not shown here in detail). It may happen that stable higher-order vortices with m_{\mathrm{1,2}}\equiv m\ge 3 exist too in this 3D system, but they were not found. This fact may be explained by scaling

\begin{array}{}\text{(25)}& N_{\text{th}}^{\text{(3D)}}\sim m^6\end{array}

for the minimum number of atoms necessary for the existence of stable vortex modes in this 3D model [40], as this very steep scaling makes it very difficult to create such stable modes with m\ge 3. On the other hand, QDs with hidden vorticity, i.e., as defined in Ref. [54], antisymmetric states with opposite vorticities in the two components,

\begin{array}{}\text{(26)}& m_1=-m_2=-\mathrm{1,}\end{array}

are completely unstable in the framework of Eqs. (22).

5.2 Two-dimensional vortex rings and necklaces

5.2.1 Basic results

In the framework of the reduction of the LHY-amended GPE system to the 2D form, see Eq. (9), families of stable QDs with embedded vorticity, both explicit (identical in both components) and hidden, defined as per Eq. (26), were explored in Ref. [40]. The extension of Eq. (9) for two components, {\psi }_{\pm }, is

\begin{array}{}\text{(27)}& \begin{array}{ll}\text{i}{\partial }_t{\psi }_{\pm }=\hfill & -\frac{1}{2}{\nabla }^2{\psi }_{\pm }+\frac{4\pi }{g}\left(|{\psi }_{\pm }{|}^2-|{\psi }_{\mp }{|}^2\right){\psi }_{\pm }\hfill \\ \hfill & +(|{\psi }_{\pm }{|}^2+|{\psi }_{\mp }{|}^2){\psi }_{\pm }\ln (|{\psi }_{\pm }{|}^2+|{\psi }_{\mp }{|}^2)\mathrm{,}\hfill \end{array}\end{array}

where g>0 is the coupling constant. It was found that Eq. (27) gives rise to stable 2D vortex-ring solutions with explicit vorticity up to m_1=m_2\equiv m=5. An essential finding is that such QDs with embedded vorticity are stable above a certain threshold value, N_{\text{th}}^{\text{(2D)}}, of the number of atoms, which scales as

\begin{array}{}\text{(28)}& N_{\text{min}}^{\text{(2D)}}\sim m^4\end{array}

with the increase of m, cf. the steeper scaling in the 3D model (25). The possibility to find stable 2D vortices with m\ge 3 is a consequence of the relatively mild scaling in Eq. (28), in comparison with Eq. (25).

Vortex QDs in a similar 2D model, with GPEs including both the LHY correction and spin-orbit coupling between the two components, were considered too [55]. In the latter case, QDs are states of the mixed-mode type, in terms of Ref. [56], with each component including terms with vorticities m=0 and +1 or -1.

5.2.2 Semidiscrete vortices

Further, a semidiscrete 2D system, which is constructed as an array of quasi-1D cigar-shaped waveguides for QDs, was introduced in Ref. [36]. In the 1D limit, the LHY correction to the 1D GPE is a quadratic term with the attraction sign [42], on the contrary to the repulsive quartic LHY term in Eq. (5) and the alternating attraction-repulsion sign of the logarithmic factor in the 2D equation (9). The scaled form of the GPE system for this system includes the MF cubic self-repulsion terms, competing with the LHY-induced quadratic self-attraction:

\begin{array}{}\text{(29)}& \begin{array}{ll}i\partial {\psi }_j=\hfill & -\frac{1}{2}{\partial }_{zz}{\psi }_j-\frac{C}{2}\left({\psi }_{j-1}-2{\psi }_j+{\psi }_{j-1}\right)\hfill \\ \hfill & +g|{\psi }_j{|}^2{\psi }_j-|{\psi }_j|{\psi }_j\mathrm{,}\hfill \end{array}\end{array}

where C>0 is the hopping rate which provides linear coupling between adjacent between cores (waveguides), the strength of the quadratic LHY term is normalized to be one, and g>0 is the strength of the cubic self-repulsion. This system gives rise to stable vortical semidiscrete QDs, with the winding number (embedded vorticity) up to m=5. Among them, there are two different species of stable semidiscrete vortex QD with m=1, of the on-site-centered and inter-site-centered types, with the vorticity pivot located, respectively, at a lattice site or between two sites, see examples in Figs. 12(c, d)]. Stable inter-site-centered vortices are found only when the hopping rate between adjacent cores, C, is very small, and they do not exist with m\ge 2. On the other hand, stable on-site-centered vortices were found for arbitrary values of C, with embedded vorticities 1\le m\le 5, i.e., in the same range as indicated above for the 2D continuum model based on Eq. (27).

5.2.3 Necklace clusters

Necklace patterns, built as ring-shaped clusters of solitons, usually appear to be unstable patterns, because interactions between adjacent soliton split the necklaces into sets of separating solitons. Nevertheless, Ref. [57] has demonstrated a possibility to construct robust necklace clusters, composed of fundamental (zero-vorticity) 2D droplets, in the framework of the model based on Eq. (27). The cluster is built of n identical QDs placed on a ring of radius R, with phase difference 2\pi m/n between adjacent QDs, where m is the overall vorticity imprinted onto the cluster. The evolution of the cluster is driven by the initial radius, R, and vorticity m, ranging from contraction to rotation or expansion. As a result, the necklaces which realize an energy minimum feature remarkable robustness, while they are strongly unstable in the framework of the 2D GPE with usual cubic nonlinearity.

Unlike the imprinted vorticity considered above, an alternative way to construct vortex clusters carrying angular momentum is offered by the ground state of a rotating trapped binary BEC, as recently demonstrated in Ref. [58]. Such a system, with attractive inter-species and repulsive intra-species interactions, is confined in a shallow HO harmonic trap with an additional repulsive Gaussian potential placed at the center. If the LHY correction is taken into account, it helps to stabilize the patterns. Numerical simulations have produced rotating necklace-like patterns composed of a few local vortices with topological charge m=1. Traces of these patterns persist in the expanding condensate if it is released into free space when the weakly confining HO trapping potential is switched off.

6 Two-dimensional vortex modes trapped in a singular potential

6.1 Formulation of the problem

The quantum collapse is a well-known peculiarity in quantum mechanics: nonexistence of the GS in 3D and 2D linear Schrödinger equations with attractive potential

\begin{array}{}\text{(30)}& U(r)=-\left(U_0/2\right)r^{-2}\mathrm{,}\end{array}

where U_0>0 is the strength of the pull to the center [59]. In 3D, the collapse occurs when U₀ exceeds a finite critical value, while in 2D the collapse happens at any U_0>0. In both 3D and 2D cases, the potential represents attraction of a particle, carrying a permanent electric dipole moment, to a central charge [60]. In 2D, the same potential ( 30) may be realized as attraction of a magnetically polarizable atom to an electric current (e.g., an electron beam) directed transversely to the system’s plane, or the attraction of an electrically polarizable atom to a uniformly charged transverse thread.

A fundamental issue is regularization of the setting, aiming to create a missing GS. A solution was proposed in Ref. [60], replacing the 3D linear Schrödinger equation by the GPE for a gas of dipole particles pulled to the center by potential (30) and stabilized by repulsive contact interactions, while the long-range DDI between the particles amount to a renormalization of the contact interaction [60]. It was thus found that, in the framework of the MF approximation, the 3D GPE creates the missing GS for arbitrarily large U₀. Further, it was demonstrated that, in terms of the many-body quantum theory, the GS, strictly speaking, does not exist in the same setting, but the interplay of the pull to the center and contact repulsion gives rise to a metastable state, separated from the collapsing state by a tall potential barrier [61].

The situation is more problematic in 2D, as the usual cubic nonlinearity, which represents the contact repulsion in the MF approximation, is not strong enough to create the GS. The problem is that the MF wave function, \psi (r), produced by GPE, gives rise to the density, {\left|\psi (r)\right|}^2, diverging \sim r^{-2} at r\to 0, in 3D and 2D alike. In terms of the integral norm,

\begin{array}{}\text{(31)}& N=\underset{r_{\text{cutoff}}\to 0}{{\displaystyle \lim }}\left[{\left(2\pi \right)}^{D-1}{\displaystyle {\int }_{r_{\text{cutoff}}}^{\infty }}{\left|\psi (r)\right|}^2{r}^{D-1}\text{d}r\right]\mathrm{,}\end{array}

where D=3 or 2 is the dimension, the density singularity \sim r^{-2} is integrable in 3D, while it gives rise to a logarithmic divergence in 2D,

\begin{array}{}\text{(32)}& N\sim \ln \left(r_{\text{cutoff}}^{-1}\right)\mathrm{.}\end{array}

The analysis of GPE demonstrates that a self-repulsive nonlinear term stronger than cubic, i.e., |\text{Ψ}{|}^{\alpha -1}\text{Ψ} with \alpha >3, gives rise to the density with singularity |\text{ψ}(r{)|}^2\sim r^{-4/(\alpha -1)}. Therefore, any value \alpha >3 provides convergence of the 2D integral norm.

Thus, a solution for the regularization of the 2D setting may be offered by the quintic defocusing nonlinearity [60], with \alpha =5. It accounts for three-body repulsive interactions in the bosonic gas [62], although the realization of this feature is the fact that three-particle collisions give rise to losses, kicking out particles from the condensate [63]. On the other hand, the LHY-induced quartic self-repulsive term in Eq. (5) may also be used for the stabilization of the 2D setting under the action of potential (30) [37], provided that the confinement in the transverse direction is realistic (not extremely tight). As mentioned above, in this case the LHY effect is accounted for by the quartic term, added to the effectively two-dimensional GPE. In this connection, it is relevant to note that, if the “fully 2D" equation (9), corresponding to the ultra-tight confinement, is insufficient to create a GS with a convergent norm in 2D. Indeed, in this case the analysis yields a density singularity {\left|\psi \right|}^2\sim r^{-2}/\ln \left(r^{-1}\right) at r\to 0, hence the 2D integral (31) is still diverging, although extremely slowly, N\sim \ln \left(\ln \left(r_{\text{cutoff}}^{-1}\right)\right), cf. Eq. (32).

Thus, the relevant two-dimensional LHY-amended GPE equation, including potential (30), takes the following form in the scaled notation [37]:

\begin{array}{}\text{(33)}& \begin{array}{ll}\text{i}\frac{\partial \psi }{\partial t}=\hfill & -\frac{1}{2}\left(\frac{{\partial }^2\psi }{\partial r^2}+\frac{1}{r}\frac{\partial \psi }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2\psi }{\partial {\theta }^2}\right)\hfill \\ \hfill & -\frac{U_0}{2r^2}\psi +\sigma {|\psi |}^2\psi +|\psi {|}^3\psi \mathrm{,}\hfill \end{array}\end{array}

which is written in polar coordinates \left(r\mathrm{,}\theta \right), coefficient \sigma accounting for the residual MF nonlinearity. In fact, the rescaling makes it possible to set \sigma =\pm 1 or 0. The case of \sigma =0 corresponds to the (nearly) exact cancellation between the intra-component repulsion and inter-component attraction, while all the nonlinearity is represented by the LHY-induced quartic term, cf. Ref. [64]. While Ref. [37] addressed Eq. (33) in its general form, we here focus on the most fundamental case of \sigma =0.

6.2 Analytical considerations

Stationary solutions to Eq. (33) with chemical potential and integer vorticity m are looked for as

\begin{array}{}\text{(34)}& \psi \left(r\mathrm{,}t\right)=\text{exp}\left(-\text{i}\mu t+\text{i}m\theta \right)u(r)\mathrm{,}\end{array}

with real radial function satisfying the equation

\begin{array}{}\text{(35)}& \mu u=-\frac{1}{2}\left(\frac{{\text{d}}^{2}u}{\text{d}{r}^2}+\frac{1}{r}\frac{\text{d}u}{\text{d}r}+\frac{U_m}{r^2}u\right)+u^4\mathrm{,}\end{array}

where, as said above, we set \sigma =0, and define a renormalized potential strength,

\begin{array}{}\text{(36)}& U_m\equiv U_0-m^2\mathrm{.}\end{array}

Simple corollaries of Eq. (35) are scaling relations which show the dependence of the solution on μ:

\begin{array}{}\text{(37)}& u(r;\mu )={(-\mu )}^{1/3}u\left({(-\mu )}^{-1/2}r;\mu =-1\right)\mathrm{,}\end{array}

\begin{array}{}\text{(38)}& N(\mu )={\left(-\mu \right)}^{-1/3}N(\mu =-1)\mathrm{.}\end{array}

A convenient substitution,

\begin{array}{}\text{(39)}& \psi \left(r\mathrm{,}\theta \mathrm{,}t\right)\equiv r^{-2/3}\phi \left(r\mathrm{,}\theta \mathrm{,}t\right)\mathrm{,}u(r)\equiv r^{-2/3}\chi (r)\mathrm{,}\end{array}

transforms Eqs. (33) and (35) into

\begin{array}{}\text{(40)}& \begin{array}{ll}\text{i}\frac{\partial \phi }{\partial t}=\hfill & -\frac{1}{2}\left(\frac{{\partial }^2}{\partial r^2}-\frac{1}{3r}\frac{\partial }{\partial r}+\frac{U_0+4/9}{r^2}+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\theta }^2}\right)\phi \hfill \\ \hfill & +\sigma \frac{|\phi {|}^2\phi }{r^{4/3}}+\frac{|\phi {|}^3\phi }{r^2}\mathrm{,}\hfill \end{array}\end{array}

\begin{array}{}\text{(41)}& \mu \text{χ}=-\frac{1}{2}\left(\frac{{\text{d}}^2\text{χ}}{\text{d}{r}^2}\text{}-\text{}\frac{1}{3r}\frac{\text{dχ}}{\text{d}r}\text{}+\text{}\frac{U_m+4/9}{r^2}\text{χ}\right)\text{}+\text{}\sigma \frac{{\text{χ}}^3}{r^{4/3}}\text{}+\text{}\frac{{\text{χ}}^4}{r^2}\mathrm{.}\end{array}

The expansion of the solution to Eq. (41) at r\to 0 yields

\begin{array}{}\text{(42)}& \text{χ}(r)={\left[\frac{1}{2}\left(U_m+\frac{4}{9}\right)\right]}^{1/3}\left(1+\frac{2\mu }{3U_m}{r}^2\right)+O(r^4)\mathrm{,}\end{array}

which is valid for U_m>0. In the interval of

(the meaning of this interval is explained below), the quadratic term in Eq. (42) is replaced, as the leading correction, by

where const remains indefinite, in terms of the expansion at r\to 0. Exactly at U_m=0, Eq. (42) is replaced by

\begin{array}{}\text{(45)}& \text{χ}(r;U_m=0)={\left(\frac{2}{9}\right)}^{1/3}\left(1+\frac{3\mu }{4}{r}^2\text{ln}\frac{r_0}{r}\right)+\mathrm{...}\mathrm{,}\end{array}

where constant r₀ is also indefinite. In all the cases, it follows from Eq. (39) that the singular form of the density at r\to 0 is

\begin{array}{}\text{(46)}& u^2(r)\approx {\left[\frac{1}{2}\left(U_m+\frac{4}{9}\right)\right]}^{2/3}{r}^{-4/3}\mathrm{,}\end{array}

with which the 2D norm (31) converges.

The asymptotic form of the solution, given by Eq. (42), is meaningful if it yields \text{χ}(r)>0 [otherwise, the derivation of Eq. (35) from Eq. (33) is irrelevant], i.e., for U_m>0, as well as for weakly negative values of the effective strength of the central potential belonging to interval (43). In the limit of U_0+4/9\to +0, Eq. (41) gives rise to an asymptotically exact solution:

\begin{array}{}\text{(47)}& \begin{array}{ll}\text{χ}(r;U_m+4/9\to 0)=\hfill & \frac{\sqrt{3}\mathrm{\Gamma }(1/3)}{\pi }{\left[-\frac{\mu }{4}\left(U_m+\frac{4}{9}\right)\right]}^{1/3}\hfill \\ \hfill & \cdot r^{2/3}{K}_{2/3}\left(\sqrt{-2\mu }r\right)\mathrm{,}\hfill \end{array}\end{array}

where \mathrm{\Gamma }(1/3)\approx 2.68 is the value of the Gamma-function, and K_{2/3} is the modified Bessel function of the second kind. The substitution of this expression in Eqs. (39) and (31) yields the respective value of the norm,

\begin{array}{}\text{(48)}& N(U_m+4/9\to 0)=\frac{{\mathrm{\Gamma }}^2(1/3)}{\sqrt{3}}\frac{{\left(U_m+4/9\right)}^{2/3}}{{\left(-2\mu \right)}^{1/3}}\mathrm{,}\end{array}

which agrees with scaling relation (38).

The counter-intuitive finding that the bound state may exist under the combined action of the defocusing quartic nonlinearity and effectively repulsive potential in interval (43), which was first reported in Ref. [60], is explained in detail in Refs. [37] and [65]). This property is specific for singular bound states (which are physically relevant ones, as they produce the convergent norm).

In the limit of r\to \infty, the asymptotic form of the solution to Eq. (41) is

\begin{array}{}\text{(49)}& \text{χ}(r)\approx {\chi }_0{r}^{1/6}\text{exp}\left(-\sqrt{-2\mu }r\right)\mathrm{,}\end{array}

where {\text{χ}}_0 is an arbitrary constant. A global picture of the nonlinear modes is produced by the TF approximation, which neglects derivatives in Eq. (41):

In the limit of r\to 0, Eq. (50) yields the same exact value of \text{χ}(r=0)={\left[\left(U_m+4/9\right)/2\right]}^{1/3} as given by Eq. (42). On the other hand, the TF approximation predicts a finite radius r₀ of the mode, neglecting the exponentially decaying tail at r\to \infty, cf. Eq. (49).

The TF approximation makes it possible to calculate the corresponding N(\mu ) dependence:

\begin{array}{}\text{(51)}& N_{\text{TF}}(\mu )=2\pi {\displaystyle {\int }_0^{r_0}}{\left[r^{-2/3}{\text{χ}}_{\text{TF}}(r)\right]}^{2}r\text{d}r=C\frac{U_m+4/9}{{(-\mu )}^{1/3}}\end{array}

with C\equiv \pi {\displaystyle {\int }_0^1}{\left(x^{-2}-1\right)}^{2/3}x\mathrm{d}x\approx 3.80, which complies with the exact scaling relation (38). The TF approximation is quite accurate for sufficiently large values of U_m. For instance, at U_m=10 and \mu =-1, a numerically found value of the norm is N_{\text{num}}\approx 41.05, while its TF-predicted counterpart is N_{\text{TF}}\approx 39.68.

6.3 Vortices

Usually, the presence of integer vorticity m\ge 1 implies that the amplitude vanishes at r\to 0 as r^m, which is necessary because the phase of the vortex field is not defined at r=0. However, the indefiniteness of the phase is also compatible with the amplitude diverging at r\to 0. In the linear equation, this divergence has the asymptotic form of the standard singular Bessel’s (alias Neumann’s) cylindrical function, Y_l(r)\sim r^{-m}, which makes the respective 2D state unnormalizable for all m\ge 1. However, in the present system, similar to Ref. [60], Eqs. (39), (42), (44), and (45) demonstrate that the interplay of the central potential and quartic nonlinearity reduces the divergence of the amplitude function to the level of r^{-2/3}, for any m, thus maintaining the normalizability of the states under the consideration.

Stationary solutions of Eq. (41) are not essentially different for m=0 (the GS) and m\ge 1 (vortices). A real difference is revealed by the analysis of their stability. Computation of stability eigenvalues for modes of small perturbations and direct simulations of the perturbed evolution demonstrate that all GSs, including those in the “counter-intuitive" interval (43), are completely stable [37].

The situation is different for the vortices. The analysis of the linearized equations for small perturbations leads to an exact result: they are stable at

\begin{array}{}\text{(52)}& U_0\ge {\left(U_0\right)}_{\text{crit}}^{(m)}=\frac{7}{9}\left(3m^2-1\right)\mathrm{,}\end{array}

while at the vortices are unstable against a perturbation eigenmode which drives the vortex’ pivot out of the central position. This prediction, including the particular values of {\left(U_0\right)}_{\text{crit}}^{(m)}, was accurately corroborated by numerical computation of stability eigenvalues, as well as by direct simulations of perturbed evolution for m=1 and 2 [37].

An example of radial profile \text{χ}(r) of the singular vortex mode, with U_0=1.53, m=1 and a finite norm, is displayed in Fig. 13(a). This value of U₀ is chosen in the instability region, close to its boundary predicted by Eq. (52), {\left(U_0\right)}_{\text{crit}}^{(1)}\approx 1.56. In accordance with the prediction provided by the calculation of eigenmodes of small perturbations, the pivot of the unstable vortex escapes from the central position, slowly moving away along a spiral trajectory, as shown in Fig. 13(b).

Eventually, the pivot is ousted to periphery, thus effectively converting the original unstable vortex into a stable GS with zero vorticity and its center of mass located at the origin, x=y=0. In the course of the simulations, a large part of the initial norm is consumed by an absorber installed at the edge of the simulation domain, to emulate losses due to outward emission of small-amplitude matter waves, in the indefinitely extended system. In particular, the evolution of the unstable vortex displayed in Fig. 13 leads to its transformation into a residual GS with norm N=1.88\approx 41\% of the initial value.

The spontaneous transformation of the vortex mode into the GS implies decay of the mode’s angular momentum. In the extended system, the momentum would be lost with emitted matter waves, while in the present setting it is gradually eliminated by the edge absorber, as shown by the simulations in Ref. [37].

On the other hand, the simulations demonstrate that perturbed vortices with m=1 and m=2 remain completely stable at, respectively, U_0>14/9 and U_0>77/9, in agreement with Eq. (52) [37].

7 Conclusion

The scenarios for the creation of stable 3D QDs with the help of corrections to the MF dynamics induced by quantum fluctuations around the MF states (the LHY effect), proposed in Ref. [13] and experimentally realized in binary BEC in Refs. [14–18], have made a crucially important contribution to the long-standing problem of the making of stable 2D and 3D soliton-like modes. A similar mechanism was also experimentally implemented in the single-component condensate with long-range interactions between atomic magnetic moments [19–26]. In the experiment, full 3D QDs [41], as well as oblate nearly-2D ones [40], have been created as the GS (ground state), i.e., without embedded vorticity. On the other hand, recent theoretical analysis has predicted stable vortical QDs, in the fully 3D and reduced 2D forms, with the unitary and multiple vorticities alike. A very recent addition to the analysis has demonstrated the existence of stable 2D vortex modes pulled to the center by the inverse-square potential (30), in which the quantum collapse is suppressed by the LHY effect [37]. These experimental and theoretical results are summarized in the present review.

There are other directions of the current work on this topic (chiefly, theoretical ones) which are not included in this brief review, but should be mentioned: two-component QDs with the linear Rabi mixing between the components [66], quasi-1D QDs and the spectrum of excitations in them [44, 67], “quantum balls" stabilized by the repulsive three-body quintic interaction acting in combination with the LHY correction and collisions between such “balls" [19, 68], QDs composed of Bose–Fermi mixtures [70–73], QDs in periodic systems [74–77], supersolid crystals built of QDs (see a recent review [78]), miscibility-immiscibility transitions in dipolar QDs [82, 83], and others. Moreover, recent non-perturbative analysis of the strong beyond-mean-field interactions provide extension of the work to the area where the perturbative LHY model is not valid [23, 79, 80].

As concerns relevant directions for the continuation of the work, the creation of the theoretically predicted stable 3D and nearly-2D droplets with embedded vorticity is a challenging aim [84]. There are also interesting possibilities for the development of the theoretical analysis, such as the further consideration of interactions of QDs, as suggested by Ref. [17], and a possibility of forming their mutually-orbiting bound states, precession of 3D droplets with embedded vorticity, considered as gyroscopes (cf. Ref. [81]), the motion of QDs in external potentials.

References

[1] B. A. Malomed, D. Mihalache, F. Wise, and L. Torner, Spatiotemporal optical solitons, J. Optics B 7(5), R53 (2005); B. A. Malomed, D. Mihalache, F. Wise, and L. Torner, On multidimensional solitons and their legacy in contemporary atomic, molecular and optical physics, J. Phys. At. Mol. Opt. Phys. 49(17), 170502 (2016)

[2] B. A. Malomed, Multidimensional solitons: Well-established results and novel findings, Eur. Phys. J. Spec. Top. 225(13–14), 2507 (2016)

[3] D. Mihalache, Multidimensional localized structures in optical and matter-wave media: A topical survey of recent literature, Rom. Rep. Phys. 69, 403 (2017)

[4] Y. Kartashov, G. Astrakharchik, B. Malomed, and L. Torner, Frontiers in multidimensional self-trapping of nonlinear fields and matter, Nat. Rev. Phys. 1(3), 185 (2019)

[5] B. A. Malomed, Vortex solitons: Old results and new perspectives, Physica D 399, 108 (2019)

[6] M. W. Ray, E. Ruokokoski, S. Kandel, M. Möttönen, and D. S. Hall, Observation of Dirac monopoles in a synthetic magnetic field, Nature 505(7485), 657 (2014)

[7] E. Radu and M. S. Volkov, Stationary ring solitons in field theory — knots and vortons, Phys. Rep. 468(4), 101 (2008)

[8] K. Tiurev, T. Ollikainen, P. Kuopanportti, M. Nakahara, D. S. Hall, and M. Möttönen, Three-dimensional skyrmions in spin-2 Bose–Einstein condensates, New J. Phys. 20(5), 055011 (2018)

[9] Y. V. Kartashov, B. A. Malomed, Y. Shnir, and L. Torner, Twisted toroidal vortex-solitons in inhomogeneous media with repulsive nonlinearity, Phys. Rev. Lett. 113(26), 264101 (2014)

[10] I. I. Smalyukh, Review: Knots and other new topological effects in liquid crystals and colloids, Rep. Prog. Phys. 83(10), 106601 (2020)

[11] V. E. Zakharov, S. V. Manakov, S. P. Novikov, and L. P. Pitaevskii, Solitons: The Inverse Scattering Method, Moscow: Nauka Publishers, 1980; New York: Consultants Bureau, 1984 (English translation)

[12] G. Fibich, The Nonlinear Schrӧdinger Equation: Singular Solutions and Optical Collapse, Heidelberg: Springer, 2015

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

[14] C. 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)

[15] P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, Bright soliton to quantum droplet transition in a mixture of Bose–Einstein condensates, Phys. Rev. Lett. 120(13), 135301 (2018)

[16] 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)

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

[18] C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, Observation of quantum droplets in a heteronuclear bosonic mixture, Phys. Rev. Research 1(3), 033155 (2019)

[19] 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)

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

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

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

[23] 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, Dilute dipolar quantum droplets beyond the extended Gross–Pitaevskii equation, Phys. Rev. Research 1(3), 033088 (2019)

[24] I. Ferrier-Barbut, M. Wenzel, F. Bӧttcher, 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)

[25] I. Ferrier-Barbut, M. Wenzel, M. Schmitt, F. Böttcher, and T. Pfau, Onset of a modulational instability in trapped dipolar Bose–Einstein condensates, Phys. Rev. A 97(1), 011604 (2018)

[26] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wachtler, L. Santos, and F. Ferlaino, Quantum-fluctuation-driven crossover from a dilute Bose–Einstein condensate to a macrodroplet in a dipolar quantum fluid, Phys. Rev. X 6(4), 041039 (2016)

[27] M. Wenzel, F. Böttcher, T. Langen, I. Ferrier-Barbut, and T. Pfau, Striped states in a many-body system of tilted dipoles, Phys. Rev. A 96(5), 053630 (2017)

[28] L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Observation of a dipolar quantum gas with metastable supersolid properties, Phys. Rev. Lett. 122(13), 130405 (2019)

[29] F. Böttcher, J. N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Transient supersolid properties in an array of dipolar quantum droplets, Phys. Rev. X 9(1), 011051 (2019)

[30] 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, Long-lived and transient supersolid behaviors in dipolar quantum gases, Phys. Rev. X 9(2), 021012 (2019)

[31] M. Guo, F. Böttcher, J. Hertkorn, J. N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, The low-energy goldstone mode in a trapped dipolar supersolid, Nature 574(7778), 386 (2019)

[32] L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Supersolid symmetry breaking from compressional oscillations in a dipolar quantum gas, Nature 574(7778), 382 (2019)

[33] G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, Excitation spectrum of a trapped dipolar supersolid and its experimental evidence, Phys. Rev. Lett. 123(5), 050402 (2019)

[34] J. Hertkorn, F. Böttcher, M. Guo, J. N. Schmidt, T. Langen, H. P. Büchler, and T. Pfau, Fate of the amplitude mode in a trapped dipolar supersolid, Phys. Rev. Lett. 123(19), 193002 (2019)

[35] T. D. Lee, K. Huang, and C. N. Yang, Eigenvalues and eigenfunctions of a Bose system of hard spheres and its lowtemperature properties, Phys. Rev. 106(6), 1135 (1957)

[36] X. Zhang, X. Xu, Y. Zheng, Z. Chen, B. Liu, C. Huang, B. A. Malomed, and Y. Li, Semidiscrete quantum droplets and vortices, Phys. Rev. Lett. 123(13), 133901 (2019)

[37] E. Shamriz, Z. Chen, and B. A. Malomed, Suppression of the quasi-two-dimensional quantum collapse in the attraction field by the Lee–Huang–Yang effect, Phys. Rev. A 101(6), 063628 (2020)

[38] V. P. Mineev, The theory of the solution of two near-ideal Bose gases, Zh. Eksp. Teor. Fiz. 67, 263 (1974) [Sov. Phys. - JETP 40, 132 (1974)

[39] T. Mithun, A. Maluckov, K. Kasamatsu, B. Malomed, and A. Khare, Modulational instability, inter-component asymmetry, and formation of quantum droplets in one-dimensional binary Bose gases, Symmetry (Basel) 12(1), 174 (2020)

[40] Y. Li, Z. Chen, Z. Luo, C. Huang, H. Tan, W. Pang, and B. A. Malomed, Two-dimensional vortex quantum droplets, Phys. Rev. A 98(6), 063602 (2018)

[41] Y. V. Kartashov, B. A. Malomed, L. Tarruell, and L. Torner, Three-dimensional droplets of swirling superfluids, Phys. Rev. A 98(1), 013612 (2018)

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

[43] S. Pilati, J. Boronat, J. Casulleras, and S. Giorgini, Quantum Monte Carlo simulation of a two-dimensional Bose gas, Phys. Rev. A 71(2), 023605 (2005)

[44] G. E. Astrakharchik and B. A. Malomed, Dynamics of one-dimensional quantum droplets, Phys. Rev. A 98(1), 013631 (2018)

[45] P. Zin, M. Pylak, T. Wasak, M. Gajda, and Z. Idziaszek, Quantum Bose–Bose droplets at a dimensional crossover, Phys. Rev. A 98, 051603(R) (2018)

[46] N. M. Hugenholtz and D. Pines, Ground-state energy and excitation spectrum of a Systemv of interacting bosons, Phys. Rev. 116(3), 489 (1959)

[47] T. Ilg, J. Kumlin, L. Santos, D. S. Petrov, and H. P. Büchler, Dimensional crossover for the beyond-mean-field correction in Bose gases, Phys. Rev. A 98(5), 051604 (2018)

[48] C. D’Errico, M. Zaccanti, M. Fattori, G. Roati, M. Inguscio, G. Modugno, and A. Simoni, Feshbach resonances in ultracold 39 K, New J. Phys. 9(7), 223 (2007)

[49] K. L. Pan, C. K. Law, and B. Zhou, Experimental and mechanistic description of merging and bouncing in head-on binary droplet collision, J. Appl. Phys. 103(6), 064901 (2008)

[50] A. Burchianti, C. D’Errico, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, A dual-species Bose–Einstein condensate with attractive interspecies interactions, Cond. Matter 5(1), 21 (2020)

[51] A. Cidrim, F. E. A. dos Santos, E. A. L. Henn, and T. Macri, Vortices in self-bound dipolar droplets, Phys. Rev. A 98(2), 023618 (2018)

[52] Kh. I. Pushkarov, D. I. Pushkarov, and I. V. Tomov, Self-action of light beams in nonlinear media: Soliton solutions, Opt. Quant. Lectr. 11(6), 471 (1979)

[53] Kh. I. Pushkarov and D. I. Pushkarov, Soliton solutiuons in some non-linear Schrӧdinger-like equations, Rep. Math. Phys. 17(1), 37 (1980)

[54] M. Brtka, A. Gammal, and B. A. Malomed, Hidden vorticity in binary Bose–Einstein condensates, Phys. Rev. A 82(5), 053610 (2010)

[55] Y. Li, Z. Luo, Y. Liu, Z. Chen, C. Huang, S. Fu, H. Tan, and B. A. Malomed, Two-dimensional solitons and quantum droplets supported by competing self- and cross-interactions in spin–orbit-coupled condensates, New J. Phys. 19(11), 113043 (2017)

[56] H. Sakaguchi, B. Li, and B. A. Malomed, Creation of two-dimensional composite solitons in spin–orbit-coupled selfattractive Bose–Einstein condensates in free space, Phys. Rev. E 89(3), 032920 (2014)

[57] Y. V. Kartashov, B. A. Malomed, and L. Torner, Metastability of quantum droplet clusters, Phys. Rev. Lett. 122(19), 193902 (2019)

[58] M. N. Tengstrand, P. Stürmer, E. Ö. Karabulut, and S. M. Reimann, Rotating binary Bose–Einstein condensates and vortex clusters in quantum droplets, Phys. Rev. Lett. 123(16), 160405 (2019)

[59] L. D. Landau and E. M. Lifshitz, Quantum Mechanics: Nonrelativistic Theory, Moscow: Nauka Publishers, 1974

[60] H. Sakaguchi and B. A. Malomed, Suppression of the quantum-mechanical collapse by repulsive interactions in a quantum gas, Phys. Rev. A 83(1), 013607 (2011)

[61] G. E. Astrakharchik and B. A. Malomed, Quantum versus mean-field collapse in a many-body system, Phys. Rev. A 92(4), 043632 (2015)

[62] F. K. Abdullaev, A. Gammal, L. Tomio, and T. Frederico, Stability of trapped Bose–Einstein condensates, Phys. Rev. A 63, 043604 (2001)

[63] E. A. Burt, R. W. Ghrist, C. J. Myatt, M. J. Holland, E. A. Cornell, and C. E. Wieman, Coherence, correlations, and collisions: What one learns about Bose–Einstein condensates from their decay, Phys. Rev. Lett. 79(3), 337 (1997)

[64] N. B. Jørgensen, G. M. Bruun, and J. J. Arlt, Dilute fluid governed by quantum fluctuations, Phys. Rev. Lett. 121(17), 173403 (2018)

[65] H. Sakaguchi and B. A. Malomed, Singular solitons, Phys. Rev. E 101(1), 012211 (2020)

[66] A. Cappellaro, T. Macrì, G. F. Bertacco, and L. Salasnich, Equation of state and self-bound droplet in Rabi-coupled Bose mixtures, Sci. Rep. 7(1), 13358 (2017)

[67] M. Tylutki, G. E. Astrakharchik, B. A. Malomed, and D. S. Petrov, Collective excitations of a one-dimensional quantum droplet, Phys. Rev. A 101, 051601(R) (2020)

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

[69] S. Gautam and S. K. Adhikari, Self-trapped quantum balls in binary Bose–Einstein condensates, J. Phys. At. Mol. Opt. Phys. 52(5), 055302 (2019)

[70] X. Cui, Spin–orbit-coupling-induced quantum droplet in ultracold Bose–Fermi mixtures, Phys. Rev. A 98(2), 023630 (2018)

[71] B. J. DeSalvo, K. Patel, J. Johansen, and C. Chin, Observation of a degenerate Fermi gas trapped by a Bose–Einstein condensate, Phys. Rev. Lett. 119(23), 233401 (2017)

[72] S. K. Adhikari, Fermionic bright soliton in a boson–fermion mixture, Phys. Rev. A 72(5), 053608 (2005)

[73] S. K. Adhikari, A self-bound matter-wave boson–fermion quantum ball, Laser Phys. Lett. 15(9), 095501 (2018)

[74] Z. Zhou, X. Yu, Y. Zou, and H. Zhong, Dynamics of quantum droplets in a one-dimensional optical lattice, Commun. Nonlinear Sci. Numer. Simul. 78, 104881 (2019)

[75] L. Dong, W. Qi, P. Peng, L. Wang, H. Zhou, and C. Huang, Multi-stable quantum droplets in optical lattices, Nonlinear Dyn. 102, 303 (2020)

[76] Y. Zheng, S. Chen, Z. Huang, S. Dai, B. Liu. Liu, Y. LI, and S. Wang, Quantum droplets in two-dimensional optical lattices, Front. Phys. 16(2), 22501 (2021)

[77] B. A. Malomed, The family of quantum droplets keeps expanding, Front. Phys. 16, 22504 (2021)

[78] F. Bottcher, J. N. Sčhmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, New states of matter with fine-tuned interactions: Quantum droplets and dipolar supersolids, arXiv: 2007.06391 (2020)

[79] C. Staudinger, F. Mazzanti, and R. E. Zillich, Self-bound Bose mixtures, Phys. Rev. A 98(2), 023633 (2018)

[80] S. Gautam and S. K. Adhikari, Limitation of the Lee–Huang–Yang interaction in forming a self-bound state in Bose–Einstein condensates, Ann. Phys. 409, 167917 (2019)

[81] R. Driben, Y. V. Kartashov, B. A. Malomed, T. Meier, and L. Torner, Soliton gyroscopes in media with spatially growing repulsive nonlinearity, Phys. Rev. Lett. 112(2), 020404 (2014)

[82] R. N. Bisset, L. A. Peña Ardila, and L. Santos, Quantum droplets of dipolar mixtures, arXiv: 2007.00404 (2020)

[83] J. C. Smith, D. Baillie, and P. B. Blakie, Quantum droplet states of a binary magnetic gas, arXiv: 2007.00366 (2020)

[84] Z. Lin, X. Xu, Z. Chen, Z. Yan, Z. Mai, and B. Liu, Two-dimensional vortex quantum droplets get thick, Commun. Nonlinear Sci. Numer. Simulat. 93, 105536 (2021)

Rights & permissions

Higher Education Press 2021