Two-dimensional solitons and quantum droplets supported by competing self- and cross-interactions in spin-orbit-coupled condensates

We study two-dimensional (2D) matter-wave solitons in spinor Bose-Einstein condensates (BECs) under the action of the spin-orbit coupling (SOC) and opposite signs of the self- and cross-interactions. Stable 2D two-component solitons of the mixed-mode (MM) type are found if the cross-interaction between the components is attractive, while the self-interaction is repulsive in each component. Stable solitons of the semi-vortex type are formed in the opposite case, under the action of competing self-attraction and cross-repulsion. The solitons exist with the total norm taking values below a collapse threshold. Further, in the case of the repulsive self-interaction and inter-component attraction, stable 2D self-trapped modes, which may be considered as quantum droplets (QDs), are created if the beyond-mean-field Lee-Huang-Yang (LHY) terms are added to the self-repulsion in the underlying system of coupled Gross-Pitaevskii equations. Stable QDs of the MM type, of a large size with an anisotropic density profile, exist with arbitrarily large values of the norm, as the LHY terms eliminate the collapse. The effect of the SOC term on characteristics of the QDs is systematically studied. We also address the existence and stability of QDs in the case of SOC with mixed Rashba and Dresselhaus terms, which makes the density profile of the QD more isotropic. Thus, QDs in the spin-orbit-coupled binary BEC are for the first time studied in the present work.


I. INTRODUCTION
Producing stable bright solitons and solitary vortices in two-and three-dimensional (2D and 3D) free space remains a challenging problem in several areas of physics, such as nonlinear optics and Bose-Einstein condensates (BECs). It is well known that the ubiquitous cubic self-attractive nonlinearity readily predicts multidimensional solitons which, however, are made completely unstable by the collapse. Diverse methods have been proposed to stabilize 2D and 3D solitons, solitary vortices, and more complex self-trapped states. In particular, stable 2D optical solitons have been predicted and created in media with saturable [1], quadratic [2], cubic-quintic [3] and nonlocal nonlinearities [4,5], which do not give rise to the collapse. Stable quasi-2D matter-wave solitons were predicted in dipolar BECs with long-range interactions [6][7][8][9][10][11][12].
It was predicted too that 2D and 3D solitons can be stabilized in spinor (two-component) BECs with the help of the spin-orbit coupling (SOC) of the Rashba type [13][14][15][16][17][18][19][20][21]. There are two types of multidimensional solitons in this system, viz., semi-vortices (SVs) and mixed modes (MMs). The SVs are characterized by an exact ansatz with vorticities (S + = 0, S − = +1) [or its flipped counterpart, with (S + = −1, S − = +0)] in its two components, φ ± [13], see Eq. (12) below. MMs do not admit a representation in the form of an exact ansatz, but they can be constructed from an input which combines terms with vorticities (S + = 0; −1) and (S − = 0; +1) in the two components [13], as per Eq. (9), see below. Still more complex self-trapped modes can be found as excited states built on top of the MMs, starting from ansatz given below by Eq. (10). Excited states can be also built on top of the SV, given by an exact ansatz with vorticities (S + = +1, S − = +2) in the two components [13], or its flipped counterpart, with (S + = −2, S − = −1). In the previously studied systems, the excited states were found to be completely unstable, unlike the fundamental MMs and SVs [13].
Recently, we have found that SOC can also support stable 2D gap solitons in the free-space spinor BEC with dipoledipole interactions (possibly combined with contact interactions, which cannot create such solitons by themselves) in the presence of the Zeeman splitting between the components [22]. In that case, the bandgap in the system's spectrum, which is populated by the gap solitons, is a result of the interplay of the SOC and Zeeman splitting. In nonlinear optics, SOC can be emulated by dispersive coupling between parallel cores in twin planar waveguides, which makes it possible to predict stable spatiotemporal solitons ("light bullets") in the system [23].
Thus, the two-component SOC systems offer a considerable potential for the creation of stable self-trapped modes in 2D and even 3D [20] settings . Thus far, the studies of these systems dealt with the case when both the self-and cross-component nonlinear interactions were attractive. On the other hand, the relative sign of the interactions in the binary BECs can be switched by means of the Feshbach resonance [24,25]. In particular, the opposite signs of the self-and cross-interactions make it possible to predict one-dimensional solitons of the symbiotic type, which may exist only in the two-component form [26,27].
One objective of this work is to consider 2D composite solitons in the spinor BEC with SOC and competing selfand cross-interaction terms with opposite signs. We find that, in the case of self-repulsion and cross-attraction, stable MMs exist below a critical value of the total norm, provided that the cross-attraction is stronger. Above the critical norm, the MMs are destroyed by the collapse. In the opposite case of the self-attraction and cross-repulsion, stable SVs exist with the total norm smaller than the critical value corresponding to the Townes solitons [28].
Another possibility for the prediction of stable multidimensional self-trapped states in the form "quantum droplets" (QDs) is offered by taking into account the Lee-Huang-Yang (LHY) corrections [31] to the coupled mean-field Gross-Pitaevskii equations (GPEs) in binary BEC. The self-repulsive LHY terms may compensate the contact attraction between the components, which, otherwise, would destabilize the solitons due to the possibility of the collapse [32][33][34]. Furthermore, the LHY terms were predicted [35][36][37][38][39][40][41][42] and experimentally demonstrated (in the ultracold gas of dysprosium atoms) [44][45][46][47][48][49] to stabilize 2D and 3D self-trapped QDs supported by the long-range dipole-dipole interactions. In particular, Ref. [35] developed a general formalism for adding small corrections to the mean-field equations. A very recent profoundly important experimental finding is that 3D QDs, with very low densities, can be also created via contact isotropic interactions in the gas of 39 K bosonic gas [50].
The second objective of the present work is to consider the role of the LHY terms in the 2D spinor condensate combining SOC with the contact self-repulsion in each component and cross-attraction between them, or vice versa (i.e., competing self-and cross-cubic terms). We find that MMs with arbitrarily large norms exist as stable modes when the LHY correction terms are present. The MMs in this case feature a large size and an anisotropic density profile, which may be oriented in any direction in the system's 2D plane, if the SOC is chosen the pure Rashba or pure Dresselhaus form (we demonstrate that the interplay of mixed Rashba-Dresselhaus SOC with the LHY terms also supports stable MM states, but the effective isotropy is broken in that case). In the course of the analysis, we take into account the effective renormalization of the strength of the LHY correction under the action of SOC, that was demonstrated in Ref. [51], in which, however, the formation of QDs in spin-orbit-coupled BEC was not addressed. In this connection, it relevant to mention that effects of another linear interaction between the two components in the spinor BEC, viz., the Rabi coupling, was considered in Ref. [52], including the effect of this interaction on the formation of the QDs under the action of the LHY terms.
The paper is structured as follows. The model and some analysis of the system are introduced in Sec. II. Stable MMs and their unstable excited states (which may be weakly unstable) are studied by means of numerically methods in Sec. III. Stable SVs, which are formed in case of the competing self-attraction and cross-repulsion, are also considered in that section. QDs of the MM type in the spin-orbit-coupled BEC are addressed in Sec. IV. This is followed by the consideration of the interplay of the mixed Rashba-Dresselhaus SOC with the LHY terms QDs of the MM type in the more general SOC system, combining the Rashba and Dresselhaus couplings, are considered in Sec. V, where it is found that QDs of the MM type persist and may be stable even when the Rashba and Dresselhaus terms are equally mixed, while the mean-field solitons (in the absence of the LHY corrections) do not exist in the latter case. The paper is concluded by Sec. VI.

II. THE MODEL
In the scaled form, the 2D system of coupled GPEs for the spinor wave function, φ = (φ + , φ − ), with competing self-and cross-interaction local cubic terms (alias terms representing self-and cross phase-modulation, i.e., SPM and XPM, respectively, in terms of nonlinear optics [53]), and SOC of the Rashba type, with strength λ, is written as where the SOC operators may be written in the Cartesian or polar coordinates: while g and γ are effective 2D coupling constant [33]. They are defined so that both are positive or negative, i.e., gγ > 0. For the time being, the LHY terms are not included. Note that the polar form of the SOC operators, if substituted in Eq. (1) demonstrates that, if any anisotropic solution {φ + (r, θ), φ − (r, θ − )} is found, it generates a family of solutions rotated by an arbitrary angle, θ 0 , in the form of By means of further rescaling, one may set λ = 1, and replace the above-mentioned values of g and γ by either γ = 1 (which we adopt below in the case of g, γ > 0), while g > 0 is varied, or g = −1, while γ < 0 is varied (which is adopted below in the case of g, γ < 0). Note that g = 0 (no SPM) may also correspond to the model of binary Fermi gases [54].
The energy of the system is In terms of the radial coordinate, r = x 2 + y 2 , the asymptotic form of the confined modes at r → ∞ is looked for as φ ± ∼ exp (−iµt − qr), where chemical potential µ is real, while radial wavenumber q may be complex, with Re(q) > 0. The substitution of this in Eq. (1) and linearization easily leads to the expression for q in terms of µ: As it follows from here, the chemical potential satisfying the localization condition, Re(q) > 0, takes values (recall we have set λ ≡ 1). In particular, the limit case of delocalized states, with Re(q) = 0, corresponds to µ = −1/2. In the case of g, γ > 0, which corresponds to the competition between the self-repulsion and cross-attraction (this sign combination may give rise to symbiotic solitons in the 1D geometry without SOC [26,27]), the system does not support SVs, but it readily creates MMs, with equal norms in the two components, It is relevant to mention that, in terms of the total norm, N ≡ N + + N − , the threshold for the onset of the critical 2D collapse in the framework of Eqs. (1) is where N Townes ≈ 5.85 is the norm of single-component Townes soliton [28]- [30]. This value, as the boundary between stable and collapsing MM states, is completely corroborated by numerical results. The initial ansatz for the fundamental MM is adopted, in polar coordinates (r, θ), as in Ref. [13], i.e., with α 1,2 > 0 and real amplitudes A 1,2 . Further, following Ref. [13], it is possible to construct excited states built on top of the MM, starting from ansatz These ansätze are used below to obtain numerically exact solutions for fundamental and excited MMs. Unlike MMs, self-trapped states in the form of SVs correspond to an exact ansatz, which is compatible with Eq. (1) [13]: with chemical potential µ < −λ 2 /2 and real amplitude functions f ± (r), with nonzero values at r = 0, which decay exponentially, as exp − −2µ − λ 2 r at r → ∞, see Eq. (5). Below, we use the input in the form of with real amplitudes A 1,2 and α 1,2 > 0, to construct numerically exact SV solutions of Eq. (1). It is relevant to stress that the norm of the wave function, N , defined above, is proportional to the number of atoms in the condensate, but is not identical to it. To estimate the actual number of atoms, we notice that, according to Ref. [16], the comparison of scaled 2D equations (1) with the underlying system of 3D GPEs written in physical units shows that the unit length in the scaled equations equations corresponds to the physical length about 1 µm. Further, assuming typical values of the transverse confinement length ≃ 3 µm and scattering length ∼ −0.1 nm for the interatomic attraction, we conclude that N = 1 in the present notation is tantamount to ≃ 3000 atoms.

III. NUMERICAL RESULTS FOR THE MEAN-FIELD SYSTEM (WITHOUT THE LHY TERMS)
A. Stable fundamental mixed modes By means of the imaginary-time method [55,56], MM solutions to Eq. (1) can be generated from initial guess (9). A typical example of the so produced MM with (N, g, γ) = (20, 0.5, 1) is displayed in Fig. 1. Note a small separation between maxima of the two components in Fig. 1(a), which is a characteristic feature of MMs [13].
The numerical solution produces a family of stable MMs in the interval of 0 ≤ g < γ ≡ 1 [recall γ ≡ 1 is imposed in Eq. (1) by rescaling]. In precise agreement with Eq. (8), for fixed g < 1 stable MMs are found at N < N thr (g), and input in the form of Eq. (9) suffers collapse at N > N thr (g). The chemical potential, µ, and energy (4) of the MMs are displayed versus N in Figs. 2(a) and (b), respectively. In the limit of N → 0, the MM turns into an infinitely broad state with an infinitely small amplitude and the above-mentioned limit value of the chemical potential, µ = −1/2. Note that Fig. 2(a) demonstrates that the MM family, which is supported by the dominant attractive XPM, satisfies the Vakhitov-Kolokolov (VK) criterion, dν/dN < 0, that, in turn, is a well-known necessary condition for the stability of the self-trapped modes [57]- [30].
Mobility of the MMs is a nontrivial issue, as Eq. (1) is not Galilean invariant. Following the pattern of Ref. [13], the mobility was studied by rewriting Eq. (1) in the reference frame moving with velocity V = (V x , V y ): Here we only consider the case of V x = 0, because the imaginary-time method could produce stationary solution of Eq. (13) solely in this case, cf. Ref. [13]. Figures 3(a,b) present a typical stable solution of Eq. (13) for N = 20 and g = 0.5. Note weakly anisotropic deformation of the 2D profile caused by the steady motion in the y direction.
The numerical solution demonstrates that, for fixing N and g, the amplitude of the moving MM, |φ max ± |, decays (while it width increases) with the increase of V y , as shown in Fig. 3(c). In this panel, the amplitude vanishes at V max y ≈ 0.72, 1.52, and 2.83, for (P, g) = (15, 0.5), (15, 0.4), and (20, 0.5), respectively, and the MMs do not exists at V y exceeding these limit values. Thus, the MM's mobility range expands with the increase of the norm and decrease of g, up to g = 0. The asymmetry of the moving solitons remains weak even when V y is close to the limit value. It is relevant to mention that the mobility remains restricted to final values of V y also in the system with the selfattractive SPM, i.e., g < 0 [13]. In agreement with a qualitative analysis of Eq. (13), moving solitons with V y < 0 are mirror images, with respect to the y-axis, of their counterparts with V y=0 > 0. Accordingly, the dependence of their amplitude on |V y | is the same as shown in Fig. 3(c).

B. Unstable excited mixed-mode states
By means of the squared-operator method [58], excited MMs, initiated by ansatz (10), have been found too. Numerical simulations demonstrate that all such excited states are unstable. Typical examples of their shape and evolution, shown in Fig. 4, demonstrate that, with the increase of g, the mode's shape becomes sharper and essentially less unstable, transforming from a "trefoil" in (a) to a 12-petal "camomile" in (b) to a 6-vane "windmill" in (c). In particular, the latter pattern seems nearly stable.

C. Stable semi-vortices
In the system with the same signs of the competing SPM (repulsive) and XPM (attractive) terms as above, i.e., g, γ > 0, imaginary-time simulations of Eq. (1), starting from the SV ansatz (12), converge to the MM solution, clearly demonstrating that the system supports no SVs in this case. SV solutions are readily generated in the opposite case, with attractive SPM and repulsive XPM terms, i.e., g, γ < 0 in Eq. (1). Fixing in this case, as said above, g = −1 by means of rescaling, we find stable SVs with the total norm taking values in the natural interval of N < N Townes ≈ 5.85. Typical examples of numerical solutions for stable SVs are displayed in Fig. 5(a), which shows that the size of the SV increases with the increase of |γ|. The stability of the SVs is corroborated by the fact that their µ(N ) dependence also satisfies the above-mentioned VK criterion, dµ/dN < 0, as seen in Fig. 5(b).
To highlight a relation between the size and structure of the SV and the magnitude of the repulsive XPM coefficient |γ|, which competes with the attractive SPM, we define effective widths of the two components and the share of the norm in the vortex component, respectively, as that there is a critical value γ cr of γ, with the widths diverging and the norm share of the vortex component, F − , approaching 0.5 at γ → γ cr . In the same limit, the chemical potential, µ, approaches the above-mentioned limit value, µ = −1/2 [see a typical example in Fig. 5(e) for N = 5], i.e., at γ = γ cr the SV family displayed in Fig. 5 degenerates into an infinitely extended state with an infinitesimal amplitude. Obviously, γ cr is a function of N , with γ cr → 0 in the linear limit, i.e., at N → 0. Figure 5(f) displays γ cr in the interval of 3 < N < N Townes , the maximum value being γ cr (N → N Townes ) ≈ 2.15.

IV. MIXED MODES UNDER THE ACTION OF THE LEE-HUANG-YANG TERMS.
In the case of the repulsive SPM, the LHY theory, which introduces effects of quantum fluctuations around the mean-field state, gives rise to post-mean-field corrections which are crucially important for the stabilization of selftrapped QDs by means of the effective higher-order self-repulsion [32]- [48]. In particular, it removes the collapse induced by the attractive XPM at N = N thr , see Eq. (8). It is relevant to mention that the LHY correction to the mean-field energy is itself affected by the presence of SOC in the spinor BEC, due to the degeneracy of atomic zero-point-energy spectrum imposed by the SOC, as shown in Ref. [51]. As a result, the interplay with SOC gives rise to an effective renormalization of the strength of the original LHY correction, in the case of moderately strong SOC. If the SOC is too strong, Fig. 6 in Ref. [51] suggests that the LHY correction may qualitatively change its form. In fact, the mean-field part of the model may also undergo a dramatic change in the case of strong SOC, replacing the regular solitons by gap solitons, as shown in Ref. [22]. Thus, both the mean-field and LHY parts of the model adopted in the present work are relevant (the LHY part is displayed below), provided that the SOC does not grow too strong. Further, the model assumes effective locality of the LHY correction to the mean-field energy. This assumption is definitely corroborated, for the nondipolar binary condensates, by the previous theoretical results [32][33][34], as well as by the very recent experimental work [50] reporting the creation of stable QDs supported by purely contact (local) interaction in the gas of 39 K atoms.
According to Ref. [33], the LHY terms, which are added to the coupled GPEs, are derived from the expression for the LHY energy. In the general case, this expression is quite complex, as shown by Eqs. (4) and (3) in Ref. [33]. However, it is shown in the same work that, in the case of where n 2D is the peak value of the two-dimensional atomic density and a is the scattering length responsible for the self-repulsion of each component, the energy can be cast in an essentially simpler form given by Eqs. (5) and (6) in Ref. [33]. A straightforward estimate demonstrates that inequality (15) amounts to the following condition for a characteristic size l 2D of a localized 2D pattern: where a ⊥ is the length of the transverse confinement of the quasi-2D condensate. Condition (16) definitely holds for any non-collapsing soliton. In this case, the nonlinear terms in the coupled GPEs can be derived from the respective part of the system's energy density, which was obtained in Ref. [33] (it effectively represents zero-point energy of the Bogoliubov modes on top of the mean-field state): where n 0 ≡ |φ f flat state corresponding to this expression, and e is the base of the natural logarithms.
Finally, coupled GPEs in the LHY form are derived by the variation of the total energy corresponding to density (17), as where the remaining scaling invariance is used to fix n 0 ≡ 1. Numerical solutions of Eq. (18) of the MM type, which are relevant here because of the opposite signs of the selfand cross-interactions in the terms ∼ g, were produced by means of the same techniques as used above for the solution of Eq. (1). Thus, the controlled parameters in this section are a set of (g, N, λ).
Numerical solution of Eqs. (18) produces anisotropic (elongated) QD density profiles when λ = 0. The size and area of the QDs are much larger than for the solitons found in above sections. In accordance with Eq. (3), in the presence of the Rashba-type SOC, the anisotropic soliton may be oriented in any direction in the (x, y) plane. To illustrate this property, typical examples of the QDs with vertical, diagonal, and horizontal directions are displayed in the first three columns of Fig. 6. When λ is smaller enough, the QD profile is nearly isotropic, see an example in the last column of Fig. 6. In the fourth row of Fig. 6, one can see that the total-density profile, n(r) = |φ + (r)| 2 + |φ − (r)| 2 , naturally has the oval shape similar to that of each component.
For the comparison of characteristics of the QD with different sets of control parameters, we plot the patterns of n(r) corresponding to different (g, N, λ) sets in Fig. 7, choosing the vertical orientation for all the QDs. The figure shows that the area and the anisotropy of the n(r) patterns strongly depend on the parameters.
In order to further quantify the dependence of the QDs on controlled parameters, we define average horizontal and vertical widths for n(r) as The degree of the anisotropy and the area of the QD can be defined by W x,y as These two characteristics, as well as the chemical potential of the QD, are displayed, as functions of g, N and λ in Fig.  8. It clearly demonstrates that: (i) the increase of the nonlinearity strength, g, reduces the anisotropy and effective area of the QD; (ii) the increase of the total norm, N , reduces the anisotropy and increases the effective, respectively; (iii) opposite to the effect of N , the increase of the SOC strength, λ, can increase the anisotropy and decreases the effective area, respectively. It is also relevant to stress that, unlike the first and second row of Fig. 8, in the last row of Fig. 8 the dependence of the chemical potential, µ, on parameters (g, N, λ) is not monotonous. In particular, the non-monotonous dependence µ(N ) implies that the QD family does not satisfy the VK criterion, dµ/dN < 0, although the entire QD family, including the segment with dµ/dN > 0, is found to be completely stable. The violation of the VK criterion as the necessary stability condition [57]- [30] is explained by the fact that the QD modes are supported not by the purely attractive interaction between the two components, but actually by its combination with the self-repulsive LHY terms.

V. QUANTUM DROPLETS UNDER THE ACTION OF MIXED RASHBA-DRESSELHAUS SOC
As shown in Ref. [16], the addition of the Dresselhaus terms to the SOC dominated by the Rashba coupling tends to delocalize the solitons in the framework of the mean-filed GPE system (which does not include the LHY terms). Therefore, the existence of stable 2D solitons under the action of the mixed Rashba-Dresselhaus SOC is a challenging issue. In this case, Eq. (18) is rewritten as where λ R and λ D are strengths of the Rashba and Dresselhaus SOC. Note that the Dresselhaus-only system, with λ R = 0 and λ D = 0, is made tantamount to its Rashba-only counterpart, with λ R = 0 and λ D = 0, by swapping x ↔ y [59]. On the other hand, the mixed form of the SOC operators in Eq. (20) obviously breaks the rotational symmetry defined above by Eq. (3). It was found in Ref. [16] that, for fixed λ R , there are critical values, (λ cr D ) SV < (λ cr D ) MM < λ R , of the Dresselhaus-SOC strength, above which the SVs and MMs, respectively, disappear through delocalization. In particular, this implies that the mean-field solitons never exist at λ D = λ R , i.e., for equal strengths of both SOC interactions.
In this section we aim to briefly explore the possibility of the existence of stable QDs in the framework of Eq. (20), which includes the LHY terms, in the case of different strengths of the Rashba and Dresselhaus interactions. Numerical solutions demonstrate that, in this LHY-affected system, QDs can be found with arbitrary values of λ R and λ R . To illustrate these findings, we currently fix λ R = 1 and vary λ D in the interval of [0,1]. Typical examples of the the respective numerical solutions for QD with different values of λ D and fixed values of (g, N ) are shown in Fig. 9(a1-a4). In particular, the orientation of the QDs tends to be diagonal. In the limit of λ D = λ R , the QDs do not suffer delocalization, unlike their mean-field counterparts (see above), but rather become isotropic. To quantify the latter property, we define the diagonal anisotropy, in terms of the total density of the two components, n(r), as where n(−y, x) is the 90 • counterclockwise rotation of n(x, y). Figure 9(b) shows ∆ as a function of λ D for fixed values of other parameters, (g, N, λ R ) = (2, 20, 1). In Fig. 9, the QD profile become practically isotropic at λ D > 0.6. Lastly, direct simulations of the QD evolution in the framework of Eq. (20) demonstrate that asymmetric QDs are stable [see Fig. 9(c)], while effectively isotropic ones tend to be unstable against spontaneous drift in the diagonal direction, as shown in Fig. 9(d).

VI. CONCLUSION
The first objective of this work is to study 2D matter-wave solitons in SOC (spin-orbit-coupled) two-component BEC with opposite signs of the cubic SPM and XPM interactions, unlike the recently studied system in which stable MM (mixed-mode) and SV (semi-vortex) solitons are created by entirely attractive interactions. In the present case, stable MMs were found provided that the XPM attraction is stronger than the SPM repulsion. These modes exist with the total norm falling below a critical value corresponding to the onset of the collapse. Excited states of the MMs were found too, being completely unstable (in some cases, their instability may be weak). In the opposite case of the attractive SPM and repulsive XPM, the system supports stable SVs with the total norm N < N Townes , where N Townes is the well-known critical norm corresponding to the Townes solitons.
Another setting studied in this work is one with the repulsive SPM and attractive XPM interactions, including the post-mean-field LHY (Lee-Huang-Yang) terms, which eliminate the collapse and help to create stable QDs (quantum drops) of the MM type with arbitrarily large values of the norm. In particular, the QDs feature a large size and an anisotropic (elliptical) density profile, that may be oriented in any direction, when the LHY terms play an essential role, and the SOC is taken in the Rashba form. The entire family of the QDs is stable in this case, even if it does not satisfy the Vakhitov-Kolokolov criterion. The growth of the SOC strength leads to the increase the anisotropy degree of the QD and decrease of its size. The QD family remains partly stable under the action of the mixed Rashba-Dresselhaus SOC.
A challenging possibility is to extend the present analysis to 3D systems. It was recently found that SOC can support metastable solitons of both SV and MM types in the 3D spinor BEC with fully attractive cubic interactions