Two-dimensional quantum droplets in binary quadrupolar condensates

We study the stability and characteristics of two-dimensional (2D) quasi-isotropic quantum droplets (QDs) of fundamental and vortex types, formed by binary Bose–Einstein condensate with magnetic quadrupole–quadrupole interactions (MQQIs). The magnetic quadrupoles are built as pairs of dipoles and antidipoles polarized along the x-axis. The MQQIs are induced by applying an external magnetic field that varies along the x-axis. The system is modeled by the Gross–Pitaevskii equations including the MQQIs and Lee-Huang-Yang correction to the mean-field approximation. Stable 2D fundamental QDs and quasi-isotropic vortex QDs with topological charges S⩽4 are produced by means of the imaginary-time-integration method for configurations with the quadrupoles polarized parallel to the system’s two-dimensional plane. Effects of the norm and MQQI strength on the QDs are studied in detail. Some results, including an accurate prediction of the effective area, chemical potential, and peak density of QDs, are obtained in an analytical form by means of the Thomas-Fermi approximation. Collisions between moving QDs are studied by means of systematic simulations.


I. INTRODUCTION
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 2 12 > g 11 g 22 , where g 12 is the strength of the interspecies attraction, while g 11 and g 22 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 39 K, with the quasi-2D [3,4] and fully 3D [5] shapes, as well as in a heteronuclear mixture of 23 Na and 87 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 .Various phenomena, such as the spinorbit 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][66][67][68][69][70][71][72][73][74][75][76][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 Ref. [84].In Ref. [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 Ref. [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][90][91][92][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 (r, θ) as 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 ≤ θ < 2π yields the respective mean values: According to Eq. ( 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 Ref. [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 Refs.[90] and [91], respectively.In Ref. [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 cf. Eq. (1).Unlike EQQI, which is attractive on the average, for MQQI the mean value of potential (3) corresponds to repulsion: cf. Eq. (2).
In the usual mean-field approximation, which is represented by the Gross-Pitaevskii (GP) equation for the singleparticle wave function ψ [97], BEC cannot maintain localized states with the help of MQQI, unless an external trapping potential is used Ref. [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 |ψ| 3 ψ, with the sign representing self-repulsion [1].In 2D, the combined mean-field-LHY term is where const is a positive coefficient, and n 0 is the reference density [9] (its typical value is ∼ 10 14 cm −3 ).Expression (5) implies the attraction at |ψ| 2 < n 0 , and repulsion at |ψ| 2 > n 0 .Finally, in 1D, the LHY term is proportional to −|ψ|ψ, 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 ∼ |ψ| 3 ψ 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 (x, y), including a gradient magnetic field oriented along the x direction.This field is induced by a tapered solenoid, as shown schematically in Fig. 1.In this configuration, the MQQI, defined as per Eq. ( 3), represents self-repulsion, while the local interaction, combining the mean-field and LHY terms according to Eq. ( 5) may account for attraction.The self-trapping of QDs is achieved through competition between these terms.
FIG. 1: The system's scheme.Quadrupoles are built as tightly-bound dipole-antidipole pairs, polarized along the x axis by means of the gradient magnetic field (with the local strength being a linear function of x), cf.Ref. [90].This field can be imposed by the tapered solenoid, as shown in the figure.The expression of the magnetic field at any point xp belonging to the x-axis of the tapered solenoid is produced by the integration of dBx where r1 and r2 are radii of the upper and lower bottom surfaces of the solenoid, h is the length of the solenoid, I is the current, and n is the number of turns of the coil.Red and green rectanges represent two components of the bosonic mixture.
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 Sec.II.Numerical findings and some analytical estimates are summarized in Sec.III.The dynamics of the QDs are the focus of Sec.IV.The paper is concluded in Sec.V.

II. THE MODEL
We aim to consider two-component QDs in the magnetic quadrupolar BEC, which is made effectively twodimensional by the application of strong confinement in the transverse direction, with confinement size a ⊥ ∼ 1 µm.The orientation of individual quadrupoles is fixed as in Fig. 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 Ref. [9].In the scaled form, they are written as [88,99,100] where ψ ± are wave functions of the two components, with density |ψ + (r)| 2 + |ψ − (r)| 2 normalized to the abovementioned reference density n 0 (cf.Eq. ( 5)), r = {x,y} is the set of coordinates measured in units of a ⊥ (i.e., essentially, measured in microns), ∇ 2 = ∂ 2 x + ∂ 2 y , the unit of time t is ∼ 1 ms, assuming that the mass of the particles is ∼ 100 proton masses, while g 11,22 > 0 and g 12 < 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 κ ∼ Q 2 /a 3 ⊥ (see Ref. [90]) represent the LHY correction [9] and MQQI, respectively.The MQQI kernel is cf. Eq. ( 3), where cos 2 θ = (x − x ′ ) 2 /(r − r ′ ) 2 and b ∼ a ⊥ is the cutoff, which is determined by the above-mentioned transverse confinement size, that is set to be 1, as said above.
Following Ref. [9], we focus the consideration on the symmetric system, with ψ + = ψ − ≡ ψ/ √ 2 and g 11 = g 22 = −g 12 ≡ g, thus reducing Eqs. ( 6) and ( 7) to a single equation in the finally scaled form: where κ remains the single control parameter, while b = 1 is fixed in Eq. ( 8), in accordance with what is said above.Dynamical invariants of Eq. ( 9) is the total norm, which is proportional to the number of atoms in the BEC, and the total energy, that includes the mean-field and LHY terms: Moving states also conserve the integral momentum, P = i ψ∇ψ * dr, where * stands for the complex conjugate.
Strictly speaking, the anisotropy present in Eqs. ( 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 Eq. ( 9) are tantamount, as said above, to 1 µm and 1 ms, respectively, while undoing the rescaling demonstrates that the value of the scaled norm N = 100 corresponds to ∼ 5000 particles in the underlying BEC.
We look for stationary solutions to Eq. ( 9), which represents stationary QDs with the chemical potential µ in the usual form, where φ (r) 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][102][103] applied to Eq. ( 9).

III. STATIONARY SOLUTIONS FOR THE QUANTUM DROPLETS (QDS)
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 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 splitstep Fourier-transform algorithm, adding random noise at the 1% amplitude level to the input, taking the input as ψ = [1 + ε ′ f (x, y)] • φ (S) (x, y), where ε ′ = 1% and f (x, y) is a random function with the value range [0, 1], which can be generated by the standard "rand" routine.Typical examples of numerically constructed stable quasi-isotropic fundamental QDs with S = 0 are displayed in Fig. 2. The parameters are (N, κ) = (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 Figs.2(c,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 where n = |ψ| 2 , A S = N/n is the effective area of the flat-top QD, and with R (r) taken as per Eq. ( 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 dE/dn = 0, which yields ) n e = exp (−κε − 1/2) . ( The chemical potential corresponding to the equilibrium density (17) is Further, as a characteristics of the QDs, we define their effective area: Comparing the results for κ = 0.1 and κ = 0.5 in Fig. 2, we conclude that the effective area is larger in the latter case, as is also seen in Fig. 3(a 1 ).Here, the solid blue line shows results of the numerical simulations, and the dotted red line shows the analytical approximation, see Eq. ( 16), (17), and ( 18), for A eff , Ip and µ, respectively.
For the fixed norm, the dependence of the effective area A eff on the MQQI strength κ is shown in Fig. 3(a 1 ), 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 Eq. ( 16), the two curves being almost identical.
Next, we show the peak density, I p (I p = |ψ| 2 max ), and the chemical potential µ of the stable quasi-isotropic QDs as functions of the MQQI strength κ in Figs.3(a 2 ) and (a 3 ), respectively.Figure 3(a 2 ) shows that I p decreases with κ, which is also in agreement with Eq. ( 17), in which the peak density decays exponentially with the increase of the MQQI strength.
In Figs.3(b 1 ) -3(b 3 ), 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 Eq. ( 16), which is a natural feature of flat-top modes.In Fig. 3(b 2 ), the peak density saturates at I p ≈ 0.448, if N is sufficiently large, which corroborates that the superfluid filling the QDs is incompressible, as demonstrated first in Ref. [1].According to Eq. ( 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 Fig. 3(b 2 )].The µ(N ) curves in Fig. 3(b 3 ) feature a negative slope, i.e., dµ/dN < 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 3 ) shows that the chemical potential saturates at µ ≈ −0.217 when the norm is large.The TF-predicted value given by Eq. ( 18) is µ e ≈ −0.220 [see the red dashed line in Fig. 3(b 3 )].As shown in Fig. 3, all the numerical results agree well with the analytical predictions provided by Eqs. ( 16), (17) and (18).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 Fig. 4, for N = 1000, and κ = 0.1.The stability of these VQDs was confirmed by direct simulations of its perturbed evolution, see Fig. 4(c), where 1% random noise is added to the input.
To systematically study characteristics of the VQD families, we define their aspect ratio and average angular momentum as follows: where W x,y and Lz represent the effective lengths of the VQDs in the x/y-direction, and the angular-momentum operator, respectively: In particular, A l < 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 Lz of stable VQDs with N = 500, defined as per Eqs.( 20)-( 24), are displayed in Figs.5(a 1 ) and 5(a 2 ), respectively, as functions of the MQQI strength κ, varying in interval 0 < κ < 1 Note that, according to Fig. 3(a 1 ), 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 1 ) and (b 2 ) display the VQD characteristics A l and Lz as functions of the norm N , with κ = 0.1.In panels 5(a 1 ), the aspect ratio attains a maximum value 1.1742 and approaches 1.0677 at κ → 1. Figure 5(a 2 ) depicts the respective orbital momentum, whose average value is 0.9955.In panels 5(b 1 ) and (b 2 ), the aspect ratio and orbital momentum have values of A l ≈ 1.0921 and Lz ≈ 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).
Simulations of in elastic collisions between moving stable VQDs with S = 1 (see Fig. 10 below) demonstrate that they may merge into a single breather with angular momentum (21) close to Lz = 2.This fact suggests that stable VQDs with S = 2 may exist too.A typical example of such a vortex state, produced by input Eq. ( 13) with S = 2, is displayed in Fig. 6, with the same parameters as in Fig. 4. The respective values given by Eqs.(20) and (21) for this double VQD are A l = 0.9718 and Lz ≈ 1.9824 at t = 10000.According to Eq. ( 20), we have found A l ≈ 1.0308 in Fig. 8(a 4 ) and A l ≈ 1.1338 in Fig. 8(b 4 ).Further, according to Eq. ( 21), we have calculated the average angular momentum of these QDs at t = 10000, Lz ≈ 2.9633 in Fig. 8(a 4 ) and Lz ≈ 3.8869 in Fig. 8(b 4 ), which exhibit a small deviation from the corresponding initial values of the angular momentum.

IV. DYNAMICS OF THE QUANTUM DROPLETS (QDS)
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 ψ (x, y, t = 0) = φ (x − x 0 , y) e −iηx + φ (x + x 0 , y) e iηx , where φ (x ± x 0 , y) represents the stationary quasi-isotropic QDs initially centered at x = ∓x 0 .It is necessary to choose x 0 large enough so that the two QDs are initially placed far from each other, avoiding any overlap.

A. Collisions between fundamental QDs
First, we consider the collisions between moving fundamental QDs (S = 0).In Ref. [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 Ref. [42], the collisional dynamics has been studied for two symmetric QDs with equal intra-species scattering lengths and equal densities of both components.In Ref. [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 , κ) = (500, 128, 0.05).The results are shown in Fig. 9. Figures 9(a 1 ) -(a 6 ) 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, Figs.9(b 1 ) -(b 6 ) 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 Refs.[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 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.

FIG. 3 :
FIG.3:The first column (a1-a3): the effective area (A eff ), peak density (Ip), and chemical potential (µ) of the fundamental quasi-isotropic QDs as functions of the strength of κ with the norm N = 500.The second column (b1-b3): the dependence of A eff , Ip and µ on the norm, N , with κ = 0.1.Here, the solid blue line shows results of the numerical simulations, and the dotted red line shows the analytical approximation, see Eq. (16),(17), and (18), for A eff , Ip and µ, respectively.

B
FIG. 4: (a) and (b): The density and phase patterns of a stable quasi-isotropic VQD (vortex QD) produced by Eq. (13) with winding number S = 1 and (N, κ) = (1000, 0.1).(c) Simulations of the perturbed evolution of the VQDs shown in panel (a) with 1% random noise added to the input.

FIG. 5 :
FIG. 5: Panels (a1) and (a2) display the aspect ratio A l and average angular momentum Lz of stable quasi-isotropic VQDs with N = 500, respectively, as functions of the MQQI strength κ.Panels (b1) and (b2) display the values, (A l , Lz), for VQDs as functions of the norm N , with κ = 0.1.The red dashed line designate the unit levels.

Figure 7 (
a) shows that the µ(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 < N < 1250, see the red line in Fig.7(a), while the black dotted segments of the curves represent unstable states, for the same parameters.In Fig.7(b), with N = 1000, the stable VQDs with S = 2 populate the area of κ < 0.45.

FIG. 7 :
FIG. 7: (a) and (b): The chemical potential of the stable VQDs with S = 2 versus N (at κ = 0.3) and κ (at N = 1000), respectively.The red solid and black dotted segments of the curves represent, respectively, stable and unstable states.