Structured heterosymmetric quantum droplets

We predict that the Lee-Huang-Yang effect makes it possible to create stable quantum droplets (QDs) in binary Bose-Einstein condensates with a heterosymmetric or heteromultipole structure, i.e., different vorticities or multipolarities in their components. The QDs feature ﬂat-top shapes when either chemical potential, μ 1 , 2 , of the droplet approaches an edge of a triangular existence domain in the ( μ 1 ,μ 2 ) plane. QDs with different vorticities of their components are stable against azimuthal perturbations, provided that the norm of one component is large. We also present multipole states in which the interaction with a strong fundamental component balances the repulsion between poles with opposite signs in the other component, leading to the formation of stable bound states. Extended stability domains are obtained for dipole QDs; tripole ones exist but are unstable, while quadrupoles are stable in a narrow region. The results uncover the existence of much richer families of stable binary QDs in comparison to states with identical components. DOI


I. Introduction
Binary Bose-Einstein Condensates (BECs) offer a unique possibility, facilitated by the use of the Feshbach-resonance technique [1][2][3], to equilibrate intra-component repulsion and inter-component attraction.It was predicted [4,5] that the balance of the weak residual attractive mean-field nonlinearity and additional self-repulsion, induced by the Lee-Huang-Yang (LHY) effect [6], i.e., a correction to the mean-field dynamics produced by Bogoliubov fluctuations, makes it possible to build self-trapped three-and two-dimensional (3D and 2D) states, called quantum droplets (QDs), as they are filled by a nearly incompressible quantum liquid.A unique asset of QDs is their stability against the collapse, that destabilizes 2D and 3D solitons of the mean-field type in diverse physical media [7].The prediction was followed by experimental creation of QDs, with both nearly-2D (oblate) [8,9] and isotropic 3D [10,11] shapes, in mixtures of two different atomic states in 39 K and in an attractive mixture of 41 K and 87 Rb atoms [12].Still earlier, it was shown that a similar mechanism stabilizes QDs in a single-component gas of dipolar atoms with long-range attraction between them [13][14][15][16][17][18][19].QDs have been also predicted in Bose-Fermi mixtures [20,21], as well as in the form of supersolids [22][23][24].
The LHY effect strongly depends on the density of states of the Bogoliubov modes, which makes it essentially different in different spatial dimensions [5].As a result, the corresponding nonlinearity in the modified Gross-Pitaevskii equations (GPEs) is quartic self-repulsion in 3D [4], while in 2D it amounts to a cubic term multiplied by a logarithmic factor [5] (in the 1D limit, it is represented by a quadratic self-attraction term [5]).The logarithmic multiplier implies switching of the 2D nonlinearity from self-attraction to repulsion with the increase of the condensate's density, which makes it possible to create stable states.In addition to the ground-state QDs, interest was drawn to QD vortices.While they are unstable in the dipolar BEC [18], stability areas were found for 3D QDs with vorticities 1 m  and 2 m  [25].The LHY-modified GPE in 2D produces stable QDs, at least, up to 5 m  [26][27][28].
To date, studies of QDs with an intrinsic topological structure were focused on binary states with equal vorticities, are possible too, but they feature a strongly reduced stability area [27]).Identical components were also assumed in the recently introduced ring-shaped droplet clusters, which may be quasi-stable states too [29,30].Admitting hetero-symmetric states, with or different symmetries in the components, may essentially expand the variety of modes existing in the system, cf.1D states with unequal components [31].In the experiment, different vorticities can be imprinted onto different components of the binary condensate by laser beams which carry the respective vorticities, coupled to different atomic states which represent the two components [32][33][34][35][36][37][38].In nonlinear optics, a similar method made it possible to create two-component vortex solitons with 1 2 m m  in photorefractive media [39].
Another possibility to create stable states with unequal vorticities in the two components is to use a trapping potential in cubic nonlinear media [40].It is relevant too to mention states composed of more than two components with different vorticities [41].
Also predicted [42] and observed [43,44] were binary heterosymmetric solitons with one fundamental component and another one carrying a dipolar structure, as well as "propellers" with a rotating dipolar structure [45].Similar to what is mentioned above for binary states with different vorticities, such complexes can be created by means of a pair of laser beams coupling to different atomic states, one beam being structureless, and the other one carrying a dipole or quadrupole transverse profile.
In this work we aim to introduce and explore hetero-symmetric QDs with different vorticities or different multipolarities (fundamental/dipole or fundamental/quadrupole) in the two components.We identify stability domains of such complexes in the plane 1 2 ( , )   of chemical potentials of the components.Quite naturally, the domain shrink with the increase of the vorticity or multipole's order.
The model, based on coupled two-dimensional GPEs with the LHY corrections, is introduced in Section II.Results demonstrating the existence and stability of the hetero-symmetric QDs, obtained by means of numerical and analytical methods, are summarized in Section III.The paper is concluded by Section IV.

II. The model
As mentioned above, in two dimensions the evolution of the two-component wave function 1,2 ( , , ) x y t  of a binary BEC is governed by coupled GPEs [5,26] with the cubic terms multiplied by logarithmic factors that represent the LHY corrections: Here the wave functions and coordinates   , x y are measured in units of 1/2 0 n and , respectively, where 0 n is the equilibrium density of the 2D QD [5], s a and a  being the scattering length of atomic collisions and the transverse confinement scale.Here the strength of the difference nonlinearity, , is scaled to be 1 , while the relative strength  of the LHY terms is kept as a free parameter, which takes generic values 1  , in the scaled notation.In is found that increase of  leads to practically linear growth of the size of the existence and stability domains, therefore we hereafter fix 1   .The system conserves total U and partial norms 1,2 U (which are proportional to numbers of atoms in the two components): as well as linear and angular momenta and the Hamiltonian, in which the last term represents the LHY-modified interaction.
To produce stationary droplet solutions of Eq. ( 1), we used the Newton iterative method.The evolution of QDs in time was simulated by means of split-step fast Fourier method applied to Eq. (1).

III. Results
First, we address axisymmetric QDs with different vorticities 1,2 m and chemical potentials 1,2 0   of their components.In polar coordinates   , r  , the corresponding solutions are sought for as where 1,2 ( ) w r are real radial profiles of the wave function.The sta- bility of the stationary solutions is analysed below by considering perturbed solutions, with azimuthal index n and growth rate  .The substitution of this in Eq. ( 1) and derivation of the corresponding Bogoliubov -de Gennes equations, by linearization of the GPEs for small perturbations , leads to an eigenvalue problem for the growth rate:    Representative profiles of QDs with 1 2 ( , ) (2, 0) m m  and (2,1) are shown in Figs.1(a,b) (the complex of the former type, with vorticity carried by a single component, may be naturally called a half-vortex [46] or semi-vortex [47]).The increase of either 1  or 2  typi- cally leads to growth of amplitudes of both components, which feature broad shapes (solid lines in the plots), due to the effect of the LHY-induced logarithmic factor in the nonlinearity in Eq. ( 1).Even if one component has zero vorticity, it develops a ring shape with a pronounced minimum at the centre.Decrease of 1,2  at fixed 2,1  leads to gradual vanishing of 2,1  .However, the limit case in which one component vanishes, is not adequately modelled by Eq. (1) [31].
In contrast to the case of QDs with identical components [5,27], in our case total norm U is a non-monotonous function of the chemical potentials, as shown in Figs.1(c,d) for Existence domains for QDs and the results of the corresponding linear-stability analysis, based on numerical solution of Eq. ( 5), are summarized in Fig. 2(a-c).Vortex QDs exist in triangular regions in the 2 1 ( , )   plane.Close to the top and right edges of the triangles,

2
 and 1  vanish, respectively, for any combination of 1 m and 2 m .At the left edge of the triangular regions, both components acquire flat-top shapes, with different amplitudes for 1 2    , hence one ex- pects to find a stability region close to this border of the existence domain.
In addition to the numerical results, the left edge of the triangular existence area, identical for all combinations of 1,2 m , can be found in an exact analytical form.To this end, in the flat-top regime, with derivatives of 1,2  neglected in the lowest approximation, Eq. ( 1) is combined with the condition of the conservation of the formal Hamiltonian, h , in the stationary version of Eq. ( 1) in the same flat-top limit [ h is given by the last two terms in the Hamiltonian density corresponding to Eq. ( 3)], cf.Ref. [25].A straightforward calculation, based on the latter principle, yields an exact result, viz., a parametric form of the left edge of the triangle, which completely coincides with its numerically found counterpart displayed in Fig. 2: where total density (the latter value is a numerical solution of an algebraic equation).While in Fig. 2(a,c) it is possible to produce numerical data close to the triangle's vertex at 1,2 0,   in Fig. 2(b) this is challenging because the existence domain strongly shrinks.
By and large, vortex QDs are prone to azimuthal instabilities with low azimuthal perturbation indices n , which, however, may be suppressed by the LHY effect [26][27][28].Lines of different colours in Fig. 2(a-c) indicate borders at which the instabilities with different values of the perturbation azimuthal index n in Eq. ( 4) disappear ( 0 n  does not produce instability).Dependencies of the real part of the perturbation growth rates , re  , on 1  for fixed 2  are presented in Fig. 3, using the same colour coding for n as in Fig. 2.These results predict a possibility to observe different dynamical scenarios by applying specific perturbation modes, corresponding to the respective values of n , to unstable QDs.In addition to the prediction of the exact border of the stability domain, the dependencies re 1 ( )   allow us to clearly identify the azimuthal index of the fastest growing perturbation mode.A final conclusion is that QDs are stable in crescentshaped domains s in Fig. 2, where re 0   for all n .We thus find that QDs with larger vorticities are vulnerable to a broader set of azimuthal perturbations.While the mode with The most destructive perturbation for the QD modes shown here has 2 n  .Nevertheless, for QDs of all the types considered here a stability domain s is identified in Figs.2(a-c), getting narrower for higher vorticities, see Figs. 2(a) and (c).Note that the stability domains cover a variety of QDs with all values of shares of the components' norms, 1,2 / U U , as they extend from the left vertex (with 2 0 w  ) to the bottom one (with 1 0 w  ).This result substantially expands the class of stable 2D QDs found in previous works.We have found similar but narrower stability domains for QDs with vorticities 1 2 , m m up to 5 and 1 2 m m  .Generally, the possibility of having narrow but nonvanishing stability areas for vortex solitons with 1 m  is a known feature of 2D models with completing nonlinearities [48].We have also found that QDs with component vorticities of opposite signs, such as Stable and unstable evolution of QDs, produced by direct simulations of Eq. ( 1), is displayed in Fig. 4. Unstable modes break into sets of fundamental (zero-vorticity) QDs, whose components typically carry different norms and fly away in tangential directions, to conserve the angular momentum.The breakup, induced by different perturbation eigenmodes, is shown in Figs.4(a-d).Further, examples of the stable evolution of QDs with ) m m  are displayed in Figs.4(e,f).Even in the presence of considerable perturbations, such QDs keep their original shape over indefinitely long times.
Another basic finding is that stable hetero-symmetric QDs can be built with different multipolarities, rather than vorticities, in the two components.In such complexes one component is a set of several poles with opposite signs (local QDs), which do not escape under the action of mutual repulsion, as they are nonlinearly coupled to the other (structureless) component, whose wave function does not feature a sign-changing pattern.The simplest representative of this class of composite QDs is the (semi-) dipole shown in Fig. 5(a).It is produced by Eq. ( 1) in the form of w .The corresponding initial guess for the dipole component was taken as , 0   .Dependences of the norm on chemical potentials for such QDs are qualitatively similar to those for the vortices, see Figs. 1(c,d).The existence domain for the semidipole QDs is shown in Fig. 2(d).The domain is substantially narrower than its counterpart for the axisymmetric vortex droplets.The dipole component of this composite mode vanishes at the right edge of the triangular existence domain, while the structureless component vanishes at the upper edge.The left edge corresponds to the structure in which both components develop broad shapes.In the latter case, the distance between poles of the dipole component, 1  , substantially increases, while component 2  transforms into two weakly overlapping in-phase fragments.Therefore, such QDs may be considered as a bound state of two fundamental QDs, whose 1  and 2  components are juxtaposed so as to be out-of-phase and in- phase, respectively.Dipole QDs are stable in a limited part of their existence domain with a complex shape, labelled s in Fig. which is bounded by the red line.It close to the bottom vertex of the existence triangle, in contrast to the stability domains for the vortices displayed above.The cause of the instability of the broad modes is elongation of the nodal line between the poles and emergence of a pattern resembling a quasi-1D dark soliton.When the nodal line becomes long enough, a snake instability sets in (which is typical for such states [49,50]), accompanied by oscillations of amplitudes of the two poles of the dipole.An example of such instability, which is actually rather weak, is presented in Fig. 5(b), while the evolution of a stable semi-dipole is shown in Fig. 5(a).Note also that an additional, very narrow, stability domain for semi-dipoles is stretched along the entire right edge of the existence triangle in Fig. 2(d).In the latter case, a weak dipole component is stabilized by the stronger structureless one.
We have also found QDs with a more complex structure.These are unstable tripoles, with one component built of three poles with alternating signs, set along a line, and (semi-) quadrupoles, with the quadrupole structure carried by a single component, see an example in Fig. 5(c).The existence domain for the semi-quadrupoles practically coincides with that for dipoles [see Fig. 2(d) above], while the stability takes place in a narrow area adjacent to the right edge of the existence triangle.The evolution of a stable semi-quadrupole is shown in Fig. 5(c), while Fig. 5(d) illustrates the decay of an unstable solution, in which four separate poles fuse into a broad pattern.

IV. Conclusion
In this work, we have uncovered a ramified variety of structured hetero-symmetric QDs in binary BECs, with different vorticities or multipolarities of the two components.In particular, these states include semi-vortices and semi-multipoles, with the respective structure carried by a single component, as well as hetero-vortical complexes, with unequal nonzero vorticities in both components.The stability of these states against azimuthal and symmetry-breaking perturbations is an example of novel phenomena offered by the LHY (Lee-Hung-Yang) corrections to the mean-field dynamics.In two-dimensional settings, the LHY effect adds the logarithmic factor to the cubic self-and cross-interactions.Undoing rescaling that led to Eq. ( 1), an estimate for the number of atoms in the QD with relevant physical parameters is , cf.Refs.[26,27].Challenging extensions of this work may be the exploration of three-dimensional composite modes, including more complex topologically organized ones, such as skyrmions [51][52][53][54][55]. Another challenging extension may be investigation of the hetero-symmetric and hetero-multipole modes, both two-and three-dimensional ones, in the framework of the full multi-body quantum theory, rather than reducing the quantum effects to the LHY corrections, cf.Ref. [56].We also anticipate the relevance of the application of the hetero-symmetry concept to complexes hold in trapping potentials, cf.Refs.[46,57]. with
the left and bottom vertices of the triangle, at which the density of either component vanishes, can be obtained from Eq. (6).In particular, the coordinates of

Fig. 2 .
Fig. 2. Stability and existence domains for QDs with different vorticities in the two components (a)-(c) and (semi-) dipole QDs (d) in the 2 1 ( , )   plane.QDs exist in triangular regions bounded by lines with black dots [the analytical form of the left edge is given by Eq. (6)].Color lines with dots in (a)-(c) are borders of domains where the instability with azimuthal indices n for vortex QDs disappear upon the increase of 1  ( 1 n  -red, 2 n  -green, 3 n  -blue, 4 n  -cyan, 5 n  -magenta).QDs are fully stable in domains marked by s , located close to the left edge of the triangle.For dipoles, the stability domain s in (d) is bounded by the red line with dots.