Two-dimensional quantum droplets in binary quadrupolar condensates

In 2015, Petrov had proposed a scheme for suppression of the collapse in the binary (two-species) Bose–Einstein condensate (BEC) in the three-dimensional (3D) space with $g_{12}^{2} \gt g_{11}g_{22}$, where g₁₂ is the strength of the inter-species attraction, while g₁₁ and g₂₂ are strengths of the intra-species self-repulsion [1]. The stabilization of the system against the collapse is provided by the Lee-Huang-Yang (LHY) effect, i.e. corrections to the mean-field BEC dynamics induced by quantum fluctuations [2]. The analysis had predicted the arrest of the collapse and formation of stable quantum droplets (QDs) filled by the ultra-diluted quantum fluid. Following the prediction, QDs have been created in a mixture of two hyperfine atomic states of ³⁹K, with the quasi-2D [3, 4] and fully 3D [5] shapes, as well as in a heteronuclear mixture of ²³Na and ⁸⁷Rb [6] atoms. The use of the LHY effect has also made it possible to create QDs in single-component dipolar bosonic gases of dysprosium [7] and erbium [8] atoms. The further analysis has demonstrates that the LHY corrections take different forms in 1D, 2D, and 3D configurations [9].

Studies of QDs have drawn much interest in very different contexts [10–57] . Various phenomena, such as the spin–orbit coupling [58], supersolidity [59, 60], and metastability [61] can be emulated in BECs under the action of the LHY corrections. Further, stable semi-discrete vortex QDs were predicted in arrays of coupled 1D cigar-shaped traps [62], and the symmetry breaking of QDs has been studied in a dual-core trap [63].

In particular, an interesting possibility is to create stable solitary vortices, i.e. QDs with embedded vorticity (traditional quantized vortices in BEC are eddy states supported by a flat background, i.e. they are 2D dark solitons [64], while solitary vortices resemble bright matter-wave solitons). Such modes are often subject to the azimuthal modulational instability that develops faster than the collapse, splitting the vortices into fragments. Therefore, stability conditions for vortex QDs become increasingly stricter with the increase of their winding number (alias the topological charge). The possibility for the stabilization of the vortex modes has been revealed by studies of BEC systems under the action of spatially periodic and quasiperiodic lattice potentials [65–77]. Recently, stable QDs in 2D square lattices [78] and ring-shaped QDs with explicit and hidden vorticity in a radially periodic potential [79, 80] have been predicted ('hidden vorticity' means as a vortex-antivortex bound state in a two-component system with zero total angular momentum [74]). In the free space, stable 2D matter-wave solitons were predicted in microwave-coupled binary condensates [81], and solutions for stable semi-vortices, i.e. 2D [82] and 3D [83] spin–orbit-coupled two-component matter-wave solitons carrying vorticity in one component, are known too.

In dipolar BECs, QDs have been observed in the form of 3D self-bound states. In the 2D geometry, the coherence of QDs was explored in [84]. In [85], the ground-state properties and intrinsic excitations of QDs in dipolar BECs were studied in detail. The previous studies demonstrate that QDs in dipolar BECs feature strong anisotropy of their density profile in the free space, but they do not carry vorticity. In [86], isotropic vortex QDs with the dipoles polarized parallel to the vortical pivot were constructed and found to be completely unstable. Thus, creating stable vortex QDs in BECs with dipole–dipole interactions (DDI) remains a fundamental problem. Recently, solutions for 2D fundamental QDs in dipolar BEC with the dipoles oriented perpendicular to the system's plane have been produced [87]. Furthermore, Li et al have proved that dipolar BEC offers a unique possibility to predict stable 2D anisotropic vortex QDs with the vortex axis oriented perpendicular to the polarization of dipoles [88]. Hence, modifying the attractive nonlinearity by introducing DDI is a promising direction for the stabilization of vortex QDs in BEC. The objective of the present work is to demonstrate that the stabilization is also possible with the help of long-range quadrupole-quadrupole interactions (QQI) between particles in BEC. Previously, the formation of solitons in QQI-coupled BEC was addressed in models that did not include the LHY terms [89–93], which made it more difficult to stabilize the self-trapped states in the BEC.

Electric quadrupoles may be built as tightly bound pairs of dipole and antidipoles directed perpendicular to the system's (x, y) plane, i.e. along the z-axis. By applying an external gradient electric field along the z-axis the electric QQI (alias EQQI) can be induced. The potential of the interaction between two electric quadrupoles in this configuration is written in the polar coordinates $\left( r,\theta \right)$ as

$$\begin{equation} U_{QQ}^{\left(\mathrm{elect}\right)}\left(r,\theta \right) = \frac{4}{3}Q^{2}r^{-5}\left( 3-5\cos ^{2}\theta \right) ,\end{equation} \tag{ 1 }$$

where Q is the quadrupole moment, r is the distance between the quadrupoles, and θ is the angle between the vector connecting the quadrupoles and the line connecting the dipole and antidipole inside the quadrupole. Averaging potentials (1) over the angular range $0\unicode{x2A7D} \theta \lt 2\pi$ yields the respective mean values:

$$\begin{equation} \frac{1}{2\pi }\int_{0}^{2\pi }U_{QQ}^{\left(\mathrm{elect}\right)}\left( r,\theta \right) \mathrm{d}\theta = -2Q^{2}r^{-5}<0.\end{equation} \tag{ 2 }$$

According to equation (2), the EQQI represent the attractive interaction.

Relatively large electric quadrupole moments are featured by some small molecules [94] and by bound states in the form of alkali diatoms [95]. These particles can be used to create BEC in which EQQI plays an essential role. In [89], the creation of quadrupolar matter-wave solitons in the 2D free space was predicted, with a conclusion that, in the presence of EQQI, the solitons feature a higher mass and stronger anisotropy than their DDI-maintained counterparts, for the same environmental parameters. Possibilities of building 2D (quasi-) discrete matter-wave solitons composed of quadrupole particles trapped in deep isotropic and anisotropic optical-lattice potentials was demonstrated in [90, 91], respectively. In [96], soliton solutions of the mixed-mode and semi-vortex types were constructed in binary quadrupolar BECs including spin–orbit coupling.

Magnetic quadrupoles can be built as tightly bound dipole-antidipole pairs directed along a particular axis (x). They can be polarized parallel to each other by a gradient magnetic field also directed along the x-axis. The potential of the magnetic quadrupole–quadrupole interactions (MQQI) is

$$\begin{equation} U_{QQ}^{\left(\mathrm{magn}\right)}\left(r,\theta \right) = \frac{4}{3}Q^{2}r^{-5}\left( 3-30\cos ^{2}\theta +35\cos ^{4}\theta \right) ,\end{equation} \tag{ 3 }$$

cf equation (1). Unlike EQQI, which is attractive on the average, for MQQI the mean value of potential (3) corresponds to repulsion:

$$\begin{equation} \frac{1}{2\pi }\int_{0}^{2\pi }U_{QQ}^{\left(\mathrm{magn}\right)}\left( r,\theta \right) \mathrm{d}\theta = \frac{3}{2}Q^{2}r^{-5}>0,\end{equation} \tag{ 4 }$$

cf equation (2).

In the usual mean-field approximation, which is represented by the Gross–Pitaevskii (GP) equation for the single-particle wave function ψ [97], BEC cannot maintain localized states with the help of MQQI, unless an external trapping potential is used [98]. However, the above-mentioned proposal to create QDs in BEC by means of the LHY correction to the mean-field approximation [1, 9] suggests a possibility to predict self-trapped objects supported by the QQI. Note that the LHY terms take different forms for different spatial dimensions. In 3D, this term in the corrected GP equation is proportional to ${|}\psi |^{3}\psi$, with the sign representing self-repulsion [1]. In 2D, the combined mean-field-LHY term is

$$\begin{equation} \mathrm{2D~local~term}~ = \mathrm{const}\cdot |\psi |^{2}\ln \left(|\psi |^{2}/n_{0}\right)\cdot \psi ,\end{equation} \tag{ 5 }$$

where $\mathrm{const}$ is a positive coefficient, and n₀ is the reference density [9] (its typical value is ${\sim}10^{14}$ cm⁻³). Expression (5) implies the attraction at $|\psi |^{2} \lt n_{0}$ , and repulsion at $|\psi |^{2} \gt n_{0}$. Finally, in 1D, the LHY term is proportional to $-{|\psi |\psi }$, with the sing corresponding to the self-attraction, which helps to build solitons and QDs in the effectively one-dimensional condensate [9].

In this work, we produce 2D solutions for self-bound states, including ones with embedded vorticity, in magnetic quadrupolar BEC, taking the corresponding LHY correction into account. Following the original proposal [1, 9] which derived the correction for the two-component BEC, we also address the two-component system. In this connection, it is relevant to stress that the LHY correction to the GP equation for the dipolar BEC has the same form (in particular, being ${\sim}{| }\psi |^{3}\psi$ in 3D) as in the case of the BEC with contact-only atomic interactions [99, 100]. Accordingly, the same is true for the BEC featuring the long-range MQQI in the combination with the local nonlinearity.

We here consider the 2D model with coordinates $\left( x,y\right)$, including a gradient magnetic field oriented along the x direction. This field is induced by a tapered solenoid, as shown schematically in figure 1. In this configuration, the MQQI, defined as per equation (3), represents self-repulsion, while the local interaction, combining the mean-field and LHY terms according to equation (5) may account for attraction. The self-trapping of QDs is achieved through competition between these terms.

Zoom In Zoom Out Reset image size

Fundamental QDs and quasi-isotropic vortical ones are produced by the analysis presented below. Effects of the norm and strength of the MQQI are systematically studied. The rest of the paper is structured as follows. The model is introduced in section 2. Numerical findings and some analytical estimates are summarized in section 3. The dynamics of the QDs are the focus of section 4. The paper is concluded in section 5.

We aim to consider two-component QDs in the magnetic quadrupolar BEC, which is made effectively two-dimensional by the application of strong confinement in the transverse direction, with confinement size $a_{\perp }\sim 1\ \mathrm{\mu}$m. The orientation of individual quadrupoles is fixed as in figure 1. The system's dynamics are governed by the 2D coupled GP equations, produced by the reduction of the full 3D GP system. The 2D equations include the LHY terms, as derived in [9]. In the scaled form, they are written as [88, 99, 100]

$$\begin{align} i\frac{\partial \psi _{+}}{\partial t} & = -\frac{1}{2}\nabla ^{2}\psi _{+}+\left( {g}_{11}{\left\vert \psi _{+}\right\vert }^{2}+{g}_{12}{ \left\vert \psi _{-}\right\vert }^{2}\right) \psi _{+}+{\gamma }\psi _{+} \left[ \left\vert \psi _{+}\left(\mathbf{r}\right)\right\vert ^{2}+\left\vert \psi _{-}\left( \mathbf{r}\right)\right\vert ^{2}\right] \nonumber\\ &\quad \times \ln \left[ \left\vert \psi _{+}\left(\mathbf{r}\right)\right\vert ^{2}+\left\vert \psi _{-}\left(\mathbf{r}\right)\right\vert ^{2}\right] +\kappa \psi _{+}\left(\mathbf{r}\right)\iint \mathrm{d}\mathbf{r^{^{\prime} }}R\left( \mathbf{r}-\mathbf{ r^{^{\prime} }}\right) \left[ \left\vert \psi _{+}\left( \mathbf{r^{^{\prime} }} \right) \right\vert ^{2}+\left\vert \psi _{-}\left( \mathbf{r^{^{\prime} }} \right) \right\vert ^{2}\right] ,\end{align} \tag{ 6 }$$

$$\begin{align} i\frac{\partial \psi _{-}}{\partial t} & = -\frac{1}{2}\nabla ^{2}\psi _{-}+\left( {g}_{22}{\left\vert \psi _{-}\right\vert }^{2}+{g}_{12}{ \left\vert \psi _{+}\right\vert }^{2}\right) \psi _{-}+{\gamma }\psi _{-} \left[ \left\vert \psi _{-}\left(\mathbf{r}\right)\right\vert ^{2}+\left\vert \psi _{+}\left( \mathbf{r}\right)\right\vert ^{2}\right] \nonumber\\ &\quad \times \ln \left[ \left\vert \psi _{-}\left(\mathbf{r}\right)\right\vert ^{2}+\left\vert \psi _{+}\left(\mathbf{r}\right)\right\vert ^{2}\right] +\kappa \psi _{-}\left(\mathbf{r}\right)\iint \mathrm{d}\mathbf{r^{^{\prime} }}R\left( \mathbf{r}-\mathbf{ r^{^{\prime} }}\right) \left[ \left\vert \psi _{-}\left( \mathbf{r^{^{\prime} }} \right) \right\vert ^{2}+\left\vert \psi _{+}\left( \mathbf{r^{^{\prime} }} \right) \right\vert ^{2}\right] ,\end{align} \tag{ 7 }$$

where $\psi _{\pm }$ are wave functions of the two components, with density $\left\vert \psi _{+}(\mathbf{r})\right\vert ^{2}+\left\vert \psi _{-}( \mathbf{r})\right\vert ^{2}$ normalized to the above-mentioned reference density n₀ (cf equation (5)), r = {x,y} is the set of coordinates measured in units of $a_{\perp }$ (i.e. essentially, measured in microns), $\nabla ^{2} = \partial _{x}^{2}+\partial _{y}^{2}$, the unit of time t is ${\sim}1$ ms, assuming that the mass of the particles is ${\sim}100$ proton masses, while $g_{11,22} \gt 0$ and $g_{12} \lt 0$, which are proportional to the respective scattering lengths of the inter-particle collisions, are strengths of the contact intra-component repulsion and inter-component attraction, respectively. Further, coefficients γ and $\kappa \sim Q^{2}/a_{\perp }^{3}$ (see [90]) represent the LHY correction [9] and MQQI, respectively. The MQQI kernel is

$$\begin{equation} R\left( \mathbf{r}-\mathbf{r^{^{\prime} }}\right) = \frac{3-30\cos ^{2}\theta +35\cos ^{4}\theta }{\left[ b^{2}+\left( \mathbf{r}-\mathbf{r^{^{\prime} }} \right) ^{2}\right] ^{5/2}},\end{equation} \tag{ 8 }$$

cf equation (3), where $\cos ^{2}\theta = \left( x-x^{^{\prime} }\right) ^{2}/(\mathbf{r}-\mathbf{r^{^{\prime} }})^{2}$ and $b\sim a_{\perp }$ is the cutoff, which is determined by the above-mentioned transverse confinement size, that is set to be 1, as said above.

Following [9], we focus the consideration on the symmetric system, with $\psi _{+} = \psi _{-}\equiv \psi /\sqrt{2}$ and ${g} _{11} = {g}_{22} = -{g}_{12}\equiv g$, thus reducing equations (6) and (7) to a single equation in the finally scaled form:

$$\begin{equation} i\frac{\partial }{\partial t}\psi = -\frac{1}{2}\triangledown ^{2}\psi +\ln \left( \left\vert \psi \right\vert ^{2}\right) \cdot \left\vert \psi \right\vert ^{2}\psi +\kappa \psi \int R\left( \mathbf{r}-\mathbf{r^{^{\prime} } }\right) \left\vert \psi \left( \mathbf{r^{^{\prime} }}\right) \right\vert ^{2}\mathrm{d} \mathbf{r^{^{\prime} }},\end{equation} \tag{ 9 }$$

where κ remains the single control parameter, while b = 1 is fixed in equation (8), in accordance with what is said above. Dynamical invariants of equation (9) is the total norm,

$$\begin{equation} N_{\text{2D}} = \int \left\vert \psi \left(\mathbf{r}\right)\right\vert ^{2}\mathrm{d}\mathbf{r,}\end{equation} \tag{ 10 }$$

which is proportional to the number of atoms in the BEC, and the total energy, that includes the mean-field and LHY terms:

$$\begin{equation} E = \frac{1}{2}\int \left[\left\vert \triangledown \psi \right\vert ^{2}\mathbf{+} \left\vert \psi \right\vert ^{4}\ln \left( \frac{\left\vert \psi \right\vert ^{2}}{\sqrt{e}}\right) +\kappa \int R\left( \mathbf{r}-\mathbf{r^{^{\prime} }} \right) ^{2}\left\vert \psi \left( \mathbf{r^{^{\prime} }}\right) \right\vert ^{2}\left\vert \psi \left( \mathbf{r}\right) \right\vert ^{2}\mathrm{d}\mathbf{r^{^{\prime} }}\right]\mathrm{d}\mathbf{r}.\end{equation} \tag{ 11 }$$

Moving states also conserve the integral momentum, $\mathbf{P} = {i}\int \psi \nabla \psi ^{\ast }\mathrm{d}\mathbf{r}$, where $\ast$ stands for the complex conjugate. Strictly speaking, the anisotropy present in equations (8) and (9) breaks the conservation of the angular momentum in the present setting. In fact, the QDs in our system are found to be quasi-isotropic (as shown below), therefore the angular momentum is approximately conserved.

The units of coordinates and time in equation (9) are tantamount, as said above, to $1\ \mathrm{\mu }$m and 1 ms, respectively, while undoing the rescaling demonstrates that the value of the scaled norm N = 100 corresponds to ${\sim}5000$ particles in the underlying BEC.

We look for stationary solutions to equation (9), which represents stationary QDs with the chemical potential µ in the usual form,

$$\begin{equation} \psi \left( \mathbf{r},t\right) = \phi \left( \mathbf{r}\right) \mathrm{e}^{-\mathrm{i}\mu t},\end{equation} \tag{ 12 }$$

where $\phi \left( \mathbf{r}\right)$ is a wave function which is real for fundamental states, and complex for vortical ones. Stationary QD solutions were produced by means of the well-known imaginary-time-propagation (ITP) method [101–103] applied to equation (9).

The strength of the MQQI, κ, and the total norm, N, are used as control parameters in the present analysis. To generate stationary states, the ITP method was initiated with the initial guess

$$\begin{equation} \phi ^{\left( S\right) }\left( x,y\right) = Ar\exp \left( -\alpha r^{2}-iS\theta \right) ,\end{equation} \tag{ 13 }$$

where α and A are real constants, r and θ are 2D polar coordinates, and integer S is vorticity (topological charge, alias winding number, S = 0 corresponding to the ground state). Further, stability of the stationary QD states was verified by direct simulations of the perturbed real-time evolution. The simulations were run by dint of the split-step Fourier-transform algorithm, adding random noise at the 1% amplitude level to the input, taking the input as $\psi = [1+\varepsilon ^{^{^{\prime} }}f(x,y)]\cdot \phi ^{\left( S\right) }\left( x,y\right)$, where $\varepsilon ^{^{^{\prime} }} = 1\%$ and $f(x,y)$ is a random function with the value range $[0,1]$, which can be generated by the standard 'rand' routine.

3.1. Fundamental QDs

Typical examples of numerically constructed stable quasi-isotropic fundamental QDs with S = 0 are displayed in figure 2. The parameters are $(N,\kappa ) = (1000,0.1)$ in (a) and (1000, 0.5) in (b). The stability of these QDs was confirmed by direct simulations of their perturbed evolution, as shown in figures 2(c) and (d). Similar to QDs in binary bosonic gases with contact-only interactions, the quasi-isotropic QDs in the quadrupolar BEC generally feature flat-top density profiles. In this case, one can apply the Thomas-Fermi (TF) approximation, neglecting the gradient term in energy (11). The corresponding approximation yields

$$\begin{equation} E_{\mathrm{TF}} = \frac{1}{2}\left( \kappa \varepsilon n^{2}+n^{2}\ln \left( \frac{n}{\sqrt{e}}\right) \right) A_{S},\end{equation} \tag{ 14 }$$

where $n = \left\vert \psi \right\vert ^{2}$, $A_{S} = N/n$ is the effective area of the flat-top QD, and

$$\begin{equation} \varepsilon \equiv \int \mathrm{d}\mathbf{r}R\left( \mathbf{r}\right) \approx 3.2,\end{equation} \tag{ 15 }$$

with $R\left( \mathbf{r}\right)$ taken as per equation (8), is an effective measure of the repulsive nonlocality. Then, one can obtain the equilibrium density n_e and the effective area A_S of the quasi-isotropic QD from the energy-minimum condition $\mathrm{d}E/\mathrm{d}n = 0$, which yields

$$\begin{align} A_{S} & = N\exp \left( \kappa \varepsilon +1/2\right),\end{align} \tag{ 16 }$$

$$\begin{align} n_{e} & = \exp \left( -\kappa \varepsilon -1/2\right) .\end{align} \tag{ 17 }$$

The chemical potential corresponding to the equilibrium density (17) is

Zoom In Zoom Out Reset image size

$$\begin{equation} \mu _{e} = \kappa \varepsilon n_{e}+n_{e}\ln n_{e}.\end{equation} \tag{ 18 }$$

Further, as a characteristics of the QDs, we define their effective area:

$$\begin{equation} A_{\text{eff}} = \frac{\left( \int \left\vert \psi \right\vert ^{2}\mathrm{d}\mathbf{r} \right) ^{2}}{\int \left\vert \psi \right\vert ^{4}\mathrm{d}\mathbf{r}}.\end{equation} \tag{ 19 }$$

Comparing the results for κ = 0.1 and κ = 0.5 in figure 2 , we conclude that the effective area is larger in the latter case, as is also seen in figure 3(a₁).

Zoom In Zoom Out Reset image size

For the fixed norm, the dependence of the effective area A_eff on the MQQI strength κ is shown in figure 3(a₁), which indicates that the effective area exponentially grows the increase of κ. Here, the solid blue and dotted red lines show, respectively, results of the numerical simulations and the analytical results, produced by equation (16), the two curves being almost identical.

Next, we show the peak density, I_p ($I_{p} = \left\vert \psi \right\vert _{_{\max }}^{2}$), and the chemical potential µ of the stable quasi-isotropic QDs as functions of the MQQI strength κ in figures 3(a₂) and (a₃), respectively. Figure 3(a₂) shows that I_p decreases with κ, which is also in agreement with equation (17), in which the peak density decays exponentially with the increase of the MQQI strength.

In figures 3(b₁)–(b₃), we fix the MQQI strength κ, to address the dependence of the effective area A_eff, peak density I_p , and chemical potential µ on the total norm N. Similar to the 2D anisotropic QDs in dipolar BECs, the effective area of the quasi-isotropic QDs grows linearly with the increase of the total norm, see equation (16), which is a natural feature of flat-top modes. In figure 3(b₂), the peak density saturates at $I_{p}\approx 0.448$, if N is sufficiently large, which corroborates that the superfluid filling the QDs is incompressible, as demonstrated first in [1]. According to equation (17), the theoretical prediction for the equilibrium density is about 0.440, which is very close to the numerical results (see the red dashed line in figure 3(b₂)). The $\mu (N)$ curves in figure 3(b₃) feature a negative slope, i.e. $\mathrm{d}\mu /\mathrm{d}N \lt 0$, satisfying the Vakhitov-Kolokolov (VK) criterion [104], which is a necessary stability condition for self-trapped modes maintained by any self-attractive nonlinearity. Figure 3(b₃) shows that the chemical potential saturates at $\mu \approx -0.217$ when the norm is large. The TF-predicted value given by equation (18) is $\mu _{e}\approx -0.220$ (see the red dashed line in figure 3(b₃)). As shown in figure 3, all the numerical results agree well with the analytical predictions provided by equations (16)–(18).

3.2. QDs with embedded vorticity

The present system maintains stable vortex QDs (VQDs) as well. A typical example of a numerically constructed one with winding number S = 1, produced by input (13) with S = 1, is displayed in figure 4, for N = 1000, and κ = 0.1. The stability of these VQDs was confirmed by direct simulations of its perturbed evolution, see figure 4(c), where 1% random noise is added to the input.

Zoom In Zoom Out Reset image size

To systematically study characteristics of the VQD families, we define their aspect ratio and average angular momentum as follows:

$$\begin{align} A_{l} & = \frac{W_{y}}{W_{x}},\end{align} \tag{ 20 }$$

$$\begin{align} \bar{L}_{z} & = \int \frac{\phi ^{\ast }\hat{L}_{z}\phi }{N}\mathrm{d}\mathbf{r},\end{align} \tag{ 21 }$$

where $W_{x,y}$ and $\bar{L}_{z}$ represent the effective lengths of the VQDs in the $x/y$-direction, and the angular-momentum operator, respectively:

$$\begin{align} W_{x} &\equiv \frac{\left( \int \left\vert \phi \left( x,y = 0\right) \right\vert ^{2}\mathrm{d}x\right) ^{2}}{\int \left\vert \left( x,y = 0\right) \right\vert ^{4}\mathrm{d}x},\end{align} \tag{ 22 }$$

$$\begin{align} W_{y} &\equiv \frac{\left( \int \left\vert \phi \left( x = 0,y\right) \right\vert ^{2}\mathrm{d}y\right) ^{2}}{\int \left\vert \left( x = 0,y\right) \right\vert ^{4}\mathrm{d}y},\end{align} \tag{ 23 }$$

$$\begin{align} \hat{L}_{z} & = -i\left( x\partial y-y\partial x\right) .\end{align} \tag{ 24 }$$

In particular, $A_{l} \lt 1$ indicates that the VQDs manifest anisotropy with the elongation along the x-direction, while $A_{l} = 1$ implies their isotropy [88].

The aspect ratio A_l and average angular momentum $\bar{L}_{z}$ of stable VQDs with N = 500, defined as per equations (20)–(24), are displayed in figures 5(a₁) and (a₂), respectively, as functions of the MQQI strength κ, varying in interval $0 \lt \kappa \lt 1$ Note that, according to figure 3(a₁), the effective area of QDs exponentially grows with the increase of κ, making the calculations more difficult. For this reason, we here present the results for κ < 1. Figures 5(b₁) and (b₂) display the VQD characteristics A_l and $\bar{L}_{z}$ as functions of the norm N, with κ = 0.1. In panel 5(a₁), the aspect ratio attains a maximum value 1.1742 and approaches 1.0677 at $\kappa \rightarrow 1$. Figure 5(a₂) depicts the respective orbital momentum, whose average value is 0.9955. In panels 5(b₁) and (b₂), the aspect ratio and orbital momentum have values of $A_{l}\approx 1.0921$ and $\bar{L}_{z}\approx 0.9957$, respectively. These findings indicates that the system with the repulsive MQQI maintains stable quasi-isotropic VQDs, whose aspect ratio is slightly larger than 1. In this connection, it is relevant to mention that for the fundamental QD modes considered above, values of the aspect ratio are also approximately equal to 1, as the density distribution in those modes seem to be nearly axisymmetric (isotropic).

Zoom In Zoom Out Reset image size

Simulations of in elastic collisions between moving stable VQDs with S = 1 (see figure 10 below) demonstrate that they may merge into a single breather with angular momentum (21) close to $\bar{L}_{z} = 2$. This fact suggests that stable VQDs with S = 2 may exist too. A typical example of such a vortex state, produced by input equation (13) with S = 2, is displayed in figure 6, with the same parameters as in figure 4. The respective values given by equations (20) and (21) for this double VQD are $A_{l} = 0.9718$ and $\bar{L}_{z}\approx 1.9824$ at $t = 10~000$.

Zoom In Zoom Out Reset image size

The family of VQDs with S = 2 are characterized by dependences of their chemical potential µ on κ and N, as shown in figures 7(a) and (b). Next, we fix the characteristic value of the MQQI strength, κ = 0.3, to study a relation between the chemical potential µ of VQDs with S = 2 and their norm N. Figure 7(a) shows that the $\mu (N)$ curve has a negative slope, satisfying the VK criterion. It is worthy to note that, for the current fixed values S = 2 and κ = 0.3, the VQDs are stable in a finite interval, $950 \lt N \lt 1250$, see the red line in figure 7(a), while the black dotted segments of the curves represent unstable states, for the same parameters. In figure 7(b), with N = 1000, the stable VQDs with S = 2 populate the area of κ < 0.45.

Zoom In Zoom Out Reset image size

We have also constructed VQD states with higher vorticities, $S\unicode{x2A7E} 3$, as shown in figure 8. Typical examples of QDs with S = 3 and S = 4, with parameters $(N,\kappa ) = (5000,0.1)$ and $(N,\kappa ) = (15~000,0.1)$, are displayed in figures 8(a₁)–(a₅) and (b₁)–(b₅), respectively. Figures 8(a₁),(b₁) and (a₅),(b₅), respectively, display the initial (t = 0) and final ($t = 10~000$) phase patterns of for S = 3 and 4. Figures 8(a₂)–(a₄) and (b₂)–(b₄) show the density patterns at t = 0, 5000, and $10~000$. The simulations demonstrate that those QDs are stable and feature quasi-isotropic shapes. According to equation (20), we have found $A_{l}\approx 1.0308$ in figure 8(a₄) and $A_{l}\approx 1.1338$ in figure 8(b₄). Further, according to equation (21), we have calculated the average angular momentum of these QDs at $t = 10~000$, $\bar{L}_{z}\approx 2.9633$ in figure 8(a₄) and $\bar{L} _{z}\approx 3.8869$ in figure 8(b₄), which exhibit a small deviation from the corresponding initial values of the angular momentum.

Zoom In Zoom Out Reset image size

It is well known that stable QDs in binary BECs can be set in motion by opposite kicks with magnitude ±η applied along the x- or y -direction, which suggests to consider collisions between two QDs moving in the opposite directions [105, 106]. The kicks can be readily applied by optical pulses shone through the tapered solenoid.

Here, we address the collisions between the QDs moving in the x direction. The corresponding initial states are constructed as

$$\begin{equation} \psi \left( x,y,t = 0\right) = \phi \left( x-x_{0},y\right) \mathrm{e}^{-\mathrm{i}\eta x}+\phi \left( x+x_{0},y\right) \mathrm{e}^{\mathrm{i}\eta x},\end{equation} \tag{ 25 }$$

where $\phi \left( x\pm x_{0},y\right)$ represents the stationary quasi-isotropic QDs initially centered at $x = \mp x_{0}$. It is necessary to choose x₀ large enough so that the two QDs are initially placed far from each other, avoiding any overlap.

4.1. Collisions between fundamental QDs

First, we consider the collisions between moving fundamental QDs (S = 0). In [63], where collisions between two-component QDs have been considered in the model with contact interactions, the simulations demonstrated a trend to inelastic outcomes produced by the collisions between QDs in the in-phase configuration. In [42], the collisional dynamics has been studied for two symmetric QDs with equal intra-species scattering lengths and equal densities of both components. In [88], the collisions between moving 2D anisotropic vortex QDs have been studied, demonstrating the formation of bound states with a vortex-antivortex-vortex structure.

We have produced typical results for the collisions generated by input (25) with ($N,x_{0},\kappa$) = ($500,128,0.05$). The results are shown in figure 9. Figures 9(a₁)–(a₆) show the result of the strongly inelastic collision between moving QDs set in motion by kicks ±0.1. The result is fusion of the two fundamental QDs into a quadrupole breather. It stretches periodically in the x and y directions, resembling the dynamics featured by a liquid drop [42]. When the initial kicks increase to ±0.12, quasi-elastic outcomes of the collision are observed. In this case, figures 9(b₁)–(b₆) show that the colliding droplets separate into another pair of droplets, slowly moving in the perpendicular direction, i.e. the collision leads to the deflection of the direction of motion by 90^o. Similar results were mentioned in [42, 87] . When the norm is fixed, the outcome of the collision is mainly determined by the velocity, i.e. by kick η [42, 63, 107]. If the collision speed exceeds a critical value, an inelastic collision takes place. A typical example of an inelastic collision of two quasi-isotropic fundamental QDs is shown in figures 9(c₁)–(c₆). After the collision, they split into several fragments. Note that the norm and the MQQI strength (both of which affect the QD's size) also affects the outcome of the collision. We do not consider these effects in detail here.

Zoom In Zoom Out Reset image size

4.2. Collision between moving vortex quantum droplets (VQDs)

Next, we consider collisions between moving VQDs. For the input given by equation (25) with ($N,x_{0},\kappa ,S$) = ($500,64,0.05,1$), several typical outcomes of the collision are presented in figure 10. In figures 10(a₁)–(a₆), the collision with a small relative speed, initiated by small kicks (here, they are ±0.01), demonstrates that the VQDs bounce back, keeping their original vorticities, as a result of the contactless collision. According to equation (21), angular momenta of the two droplets in figure 10(a₆) are $\bar{L}_{z,\mathrm{ right}} = \bar{L}_{z,\mathrm{left}}\approx 1.005$. It is worth noting that the mutual rebound does not occur in the in the case of the collision of fundamental QDs, as seen in figure 9. By gradually increasing the kicks to ±0.04, the colliding VQDs merge and eventually split into a couple of non-vortical fragments, see figures 10(b₁)–(b₆). The simulations demonstrate a trend to the merger of the colliding VQDs into a vortex breather. In particular, this outcome, corresponding to the initial kicks ±0.06, is observed in figures 10(c₁)–(c₆). The emerging breather absorbs 86.34% of the total norm in figure 10(c₆). The angular momentum of the breather is $\bar{L} _{z}\approx 1.6860$, which is relatively close to 2. The VQDs being quasi-isotropic, the average angular momentum is approximately conserved, as mentioned above. This result suggests the existence of stable VQDs with vorticity S = 2. As demonstrated above, stable stationary VQDs withe S = 2 are indeed available.

Zoom In Zoom Out Reset image size

The purpose of this work is to establish the stability and characteristics of 2D quasi-isotropic QDs in magnetic quadrupolar BECs. The system is modeled by the coupled GP including the LHY terms, which represent the correction to the mean field theory produced by quantum fluctuations. Stable quasi-isotropic QDs of the fundamental and vortex types (with topological charges $S\unicode{x2A7D} 4$) have been produced by means of imaginary-time simulations. The effect of the norm and strength of the MQQIs, κ, on the QDs were studied in detail. It was thus found that the QD's size increases with the increase of κ. For stable QDs, the dependence between the chemical potential and total norm obeys the VK criterion. The dependence of the effective size and peak density of the QDs on κ and norm N have been studied too. The results of the numerical simulation are consistent with the analytical predictions produced by the TF approximation. Typically, the predicted stable QDs have the size in the range of $20-50\,\mathrm{\mu}$m, with the number of particles in the range of $10^{4}-10^{5}$, with density ${\sim}10^{14}$ cm⁻³.

Collisions between moving QDs have been also studied systematically, for the fundamental and vortex droplets alike. For moving fundamental QDs, the strongly inelastic, quasi-elastic, and inelastic outcomes of the collision have been observed. For moving vortex QDs, a noteworthy results is the contactless rebound of the colliding droplets.

The present analysis can be extended further. In particular, it is relevant to look for essentially anisotropic QDs, especially higher-order vortical ones, in the framework of the present model. A challenging option is to seek for stable vortex QDs in the full three-dimensional setting including MQQI.

This work was supported by the Natural Science Foundation of GuangDong Province through Grant Nos. 2024A1515030131 and 2021A1515010214, GuangDong Basic and Applied Basic Research Foundation through Grant No. 2023A1515110198, National Natural Science Foundation of China through Grant Nos. 12274077 and 11905032, the Research Fund of the Guangdong-Hong Kong-Macao Joint Laboratory for Intelligent Micro-Nano Optoelectronic Technology through Grant No. 2020B1212030010, and Israel Science Foundation through Grant No. 1695/22.

No new data were created or analysed in this study.