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

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–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–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, ${\phi }_{\pm }$ [13], see equation (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 equation (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 equation (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 dipole–dipole 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 self- and 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–30].

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–34]. Furthermore, the LHY terms were predicted [35–42] and experimentally demonstrated (in the ultracold gas of dysprosium atoms) [44–49] to stabilize 2D and 3D self-trapped QDs supported by the long-range dipole–dipole interactions. In particular [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 ³⁹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 [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 [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 section 2. Stable MMs and their unstable excited states (which may be weakly unstable) are studied by means of numerically methods in section 3. 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 section 4. 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 section 5, 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 section 6.

In the scaled form, the 2D system of coupled GPEs for the spinor wave function, ${\boldsymbol{\phi }}=({\phi }_{+},{\phi }_{-})$, 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

$$\begin{eqnarray}{\rm{i}}\displaystyle \frac{\partial {\phi }_{+}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{+}+(g| {\phi }_{+}{| }^{2}-\gamma | {\phi }_{-}{| }^{2}){\phi }_{+}+\lambda \hat{D}{\phi }_{-},\\ {\rm{i}}\displaystyle \frac{\partial {\phi }_{-}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{-}+(g| {\phi }_{-}{| }^{2}-\gamma | {\phi }_{+}{| }^{2}){\phi }_{-}-\lambda {\hat{D}}^{* }{\phi }_{+},\end{eqnarray} \tag{ 1 }$$

where the SOC operators may be written in the Cartesian or polar coordinates:

$$\begin{eqnarray}\hat{D} & = & {\partial }_{x}-{\rm{i}}{\partial }_{{y}}\equiv {{\rm{e}}}^{-{\rm{i}}\theta }\left(\displaystyle \frac{\partial }{\partial r}-\displaystyle \frac{{\rm{i}}}{r}\displaystyle \frac{\partial }{\partial \theta }\right),\\ {\hat{D}}^{* } & = & {\partial }_{x}+{\rm{i}}{\partial }_{{y}}\equiv {{\rm{e}}}^{{\rm{i}}\theta }\left(\displaystyle \frac{\partial }{\partial r}+\displaystyle \frac{{\rm{i}}}{r}\displaystyle \frac{\partial }{\partial \theta }\right),\end{eqnarray} \tag{ 2 }$$

while g and γ are effective 2D coupling constant [33]. They are defined so that both are positive or negative, i.e., $g\gamma \gt 0$. For the time being, the LHY terms are not included. Note that the polar form of the SOC operators, if substituted in equation (1) demonstrates that, if any anisotropic solution $\{{\phi }_{+}(r,\theta ),{\phi }_{-}(r,{\theta }_{-})\}$ is found, it generates a family of solutions rotated by an arbitrary angle, ${\theta }_{0}$, in the form of

$$\begin{eqnarray}&&{({\phi }_{+},{\phi }_{-})}_{\mathrm{rotated}}={\phi }_{+}(r,\theta +{\theta }_{0}),\,{{\rm{e}}}^{-{\rm{i}}{\theta }_{0}}{\phi }_{-}(r,\theta +{\theta }_{0}).\end{eqnarray} \tag{ 3 }$$

By means of further rescaling, one may set $\lambda =1$, and replace the above-mentioned values of g and γ by either $\gamma =1$ (which we adopt below in the case of $g,\gamma \gt 0$), while $g\gt 0$ is varied, or $g=-1$, while $\gamma \lt 0$ is varied (which is adopted below in the case of $g,\gamma \lt 0$). Note that g = 0 (no SPM) may also correspond to the model of binary Fermi gases [54].

The energy of the system is

$$\begin{eqnarray}E & = & \displaystyle \int \displaystyle \int \left\{\displaystyle \frac{1}{2}(| {\rm{\nabla }}{\phi }_{+}{| }^{2}+| {\rm{\nabla }}{\phi }_{-}{| }^{2})+\displaystyle \frac{g}{2}(| {\phi }_{+}{| }^{4}+| {\phi }_{-}{| }^{4})-\gamma | {\phi }_{+}{| }^{2}| {\phi }_{-}{| }^{2}\right.\\ & & +\lambda ({\phi }_{+}^{* }\hat{D}{\phi }_{-}-{\phi }_{-}^{* }{\hat{D}}^{* }{\phi }_{+})\}{\rm{d}}{x}{\rm{d}}{y}.\end{eqnarray} \tag{ 4 }$$

In terms of the radial coordinate, $r=\sqrt{{x}^{2}+{y}^{2}}$, the asymptotic form of the confined modes at $r\to \infty$ is looked for as ${\phi }_{\pm }\sim \exp (-{\rm{i}}\mu t-{qr})$, where chemical potential μ is real, while radial wavenumber q may be complex, with $\mathrm{Re}(q)\gt 0$ . The substitution of this in equation (1) and linearization easily leads to the expression for q in terms of μ:

$$\begin{eqnarray}&&q=\pm {\rm{i}}\lambda +\sqrt{2\mu +{\lambda }^{2}}.\end{eqnarray} \tag{ 5 }$$

As it follows from here, the chemical potential satisfying the localization condition, $\mathrm{Re}(q)\gt 0$, takes values

$$\begin{eqnarray}&&\mu \lt -{\lambda }^{2}/2\equiv -1/2\end{eqnarray} \tag{ 6 }$$

(recall we have set $\lambda \equiv 1$). In particular, the limit case of delocalized states, with $\mathrm{Re}(q)=0$, corresponds to $\mu =-1/2$.

In the case of $g,\gamma \gt 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,

$$\begin{eqnarray}&&{N}_{+}={N}_{-}\equiv \int \int | {\phi }_{\pm }(x,y){| }^{2}{\rm{d}}{x}{\rm{d}}{y}.\end{eqnarray} \tag{ 7 }$$

It is relevant to mention that, in terms of the total norm, $N\equiv {N}_{+}+{N}_{-}$, the threshold for the onset of the critical 2D collapse in the framework of equations (1) is

$$\begin{eqnarray}&&{N}_{\mathrm{thr}}(g)=2{(\gamma -g)}^{-1}{N}_{\mathrm{Townes}}\approx 11.70{(\gamma -g)}^{-1},\end{eqnarray} \tag{ 8 }$$

where ${N}_{\mathrm{Townes}}\approx 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,\theta )$, as in [13], i.e.,

$$\begin{eqnarray}{\phi }_{+}(t=0) & = & {A}_{1}{\rm{\exp }}(-{\alpha }_{1}{r}^{2})-{A}_{2}r{\rm{\exp }}(-{\rm{i}}\theta -{\alpha }_{2}{r}^{2}),\\ {\phi }_{-}(t=0) & = & {A}_{1}{\rm{\exp }}(-{\alpha }_{1}{r}^{2})+{A}_{2}r{\rm{\exp }}({\rm{i}}\theta -{\alpha }_{2}{r}^{2}),\end{eqnarray} \tag{ 9 }$$

with ${\alpha }_{\mathrm{1,2}}\gt 0$ and real amplitudes ${A}_{\mathrm{1,2}}$. Further, following [13], it is possible to construct excited states built on top of the MM, starting from ansatz

$$\begin{eqnarray}{\phi }_{+}(t=0) & = & {A}_{1}r{\rm{\exp }}({\rm{i}}\theta -{\alpha }_{1}{r}^{2})-{A}_{2}{r}^{2}{\rm{\exp }}(-2{\rm{i}}\theta -{\alpha }_{2}{r}^{2}),\\ {\phi }_{-}(t=0) & = & {A}_{1}r{\rm{\exp }}(-{\rm{i}}\theta -{\alpha }_{1}{r}^{2})+{A}_{2}{r}^{2}{\rm{\exp }}(2{\rm{i}}\theta -{\alpha }_{2}{r}^{2}).\end{eqnarray} \tag{ 10 }$$

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 equation (1) [13]:

$$\begin{eqnarray}&&{\phi }_{+}={{\rm{e}}}^{-{\rm{i}}\mu t}{f}_{+}(r),\,{\phi }_{-}={{\rm{e}}}^{-{\rm{i}}\mu t+{\rm{i}}\theta }{{rf}}_{-}(r),\end{eqnarray} \tag{ 11 }$$

with chemical potential $\mu \lt -{\lambda }^{2}/2$ and real amplitude functions ${f}_{\pm }(r)$, with nonzero values at r = 0, which decay exponentially, as $\exp (-\sqrt{-2\mu -{\lambda }^{2}}r)$ at $r\to \infty$, see equation (5). Below, we use the input in the form of

$$\begin{eqnarray}&&{\phi }_{+}(t=0)={A}_{1}{\rm{\exp }}(-{\alpha }_{1}{r}^{2}),\quad {\phi }_{-}(t=0)={A}_{2}r{\rm{\exp }}({\rm{i}}\theta -{\alpha }_{2}{r}^{2}),\end{eqnarray} \tag{ 12 }$$

with real amplitudes ${A}_{\mathrm{1,2}}$ and ${\alpha }_{\mathrm{1,2}}\gt 0$, to construct numerically exact SV solutions of equation (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 [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 $\mu$ m. Further, assuming typical values of the transverse confinement length ≃3 $\mu {\rm{m}}$ and scattering length $\sim -0.1\mathrm{nm}$ for the interatomic attraction, we conclude that N = 1 in the present notation is tantamount to $\simeq \,3000$ atoms.

3.1. Stable fundamental MMs

By means of the imaginary-time method [55, 56], MM solutions to equation (1) can be generated from initial guess (9). A typical example of the so produced MM with $(N,g,\gamma )=(20,0.5,1)$ is displayed in figure 1. Note a small separation between maxima of the two components in figure 1(a), which is a characteristic feature of MMs [13].

Zoom In Zoom Out Reset image size

The numerical solution produces a family of stable MMs in the interval of $0\leqslant g\lt \gamma \equiv 1$ (recall $\gamma \equiv 1$ is imposed in equation (1) by rescaling). In precise agreement with equation (8), for fixed $g\lt 1$ stable MMs are found at $N\lt {N}_{\mathrm{thr}}(g)$, and input in the form of equation (9) suffers collapse at $N\gt {N}_{\mathrm{thr}}(g)$. The chemical potential, μ, and energy (4) of the MMs are displayed versus N in figures 2(a) and (b), respectively. In the limit of $N\to 0$, the MM turns into an infinitely broad state with an infinitely small amplitude and the above-mentioned limit value of the chemical potential, $\mu =-1/2$. Note that figure 2(a) demonstrates that the MM family, which is supported by the dominant attractive XPM, satisfies the Vakhitov–Kolokolov (VK) criterion, ${\rm{d}}\nu /{\rm{d}}{N}\lt 0$ , that, in turn, is a well-known necessary condition for the stability of the self-trapped modes [30–57].

Zoom In Zoom Out Reset image size

Mobility of the MMs is a nontrivial issue, as equation (1) is not Galilean invariant. Following the pattern of [13], the mobility was studied by rewriting equation (1) in the reference frame moving with velocity ${\bf{V}}=({V}_{x},{V}_{y})$:

$$\begin{eqnarray}{\rm{i}}\displaystyle \frac{\partial {\phi }_{+}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{+}+{\rm{i}}({\bf{V}}\cdot {\rm{\nabla }}){\phi }_{+}+(g| {\phi }_{+}{| }^{2}-\gamma | {\phi }_{-}{| }^{2}){\phi }_{+}+\lambda \hat{D}{\phi }_{-},\\ {\rm{i}}\displaystyle \frac{\partial {\phi }_{-}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{-}+{\rm{i}}({\bf{V}}\cdot {\rm{\nabla }}){\phi }_{-}+(g| {\phi }_{-}{| }^{2}-\gamma | {\phi }_{+}{| }^{2}){\phi }_{-}-\lambda {\hat{D}}^{* }{\phi }_{+}.\end{eqnarray} \tag{ 13 }$$

Here we only consider the case of ${V}_{x}=0$, because the imaginary-time method could produce stationary solution of equation (13) solely in this case, see [13]. Figures 3(a), (b) present a typical stable solution of equation (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.

Zoom In Zoom Out Reset image size

The numerical solution demonstrates that, for fixing N and g, the amplitude of the moving MM, $| {\phi }_{\pm }^{\max }|$, decays (while it width increases) with the increase of V_y, as shown in figure 3(c). In this panel, the amplitude vanishes at ${V}_{y}^{\max }\approx 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 self-attractive SPM, i.e., $g\lt 0$ [13]. In agreement with a qualitative analysis of equation (13), moving solitons with ${V}_{y}\lt 0$ are mirror images, with respect to the y-axis, of their counterparts with ${V}_{y=0}\gt 0$. Accordingly, the dependence of their amplitude on $| {V}_{y}|$ is the same as shown in figure 3(c).

3.2. Unstable excited MM 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 figure 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.

Zoom In Zoom Out Reset image size

3.3. Stable SVs

In the system with the same signs of the competing SPM (repulsive) and XPM (attractive) terms as above, i.e., $g,\gamma \gt 0$, imaginary-time simulations of equation (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,\gamma \lt 0$ in equation (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\lt {N}_{\mathrm{Townes}}\approx 5.85$. Typical examples of numerical solutions for stable SVs are displayed in figure 5(a), which shows that the size of the SV increases with the increase of $| \gamma |$. The stability of the SVs is corroborated by the fact that their $\mu (N)$ dependence also satisfies the above-mentioned VK criterion, ${\rm{d}}\mu /{\rm{d}}{N}\lt 0$, as seen in figure 5(b).

Zoom In Zoom Out Reset image size

To highlight a relation between the size and structure of the SV and the magnitude of the repulsive XPM coefficient $| \gamma |$, 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

$$\begin{eqnarray}&&{\bar{r}}_{\pm }={\left(\displaystyle \frac{{\int }_{0}^{\infty }| {\phi }_{\pm }(r){| }^{2}{r}^{3}{\rm{d}}{r}}{{\int }_{0}^{\infty }| {\phi }_{\pm }(r){| }^{2}{r}{\rm{d}}{r}}\right)}^{1/2},\quad {F}_{-}=\displaystyle \frac{{N}_{-}}{N}.\end{eqnarray} \tag{ 14 }$$

Figure 5(c) displays the widths as functions of γ for $(N,g)=(5,-1)$ (naturally, the width is larger for the vortical component, ${\phi }_{-}$, which has the 'hole' at its center). The norm share of the vortex component is displayed, as a function of γ, in figure 5(d) for $(N,g)=(4,-1)$ and $(5,-1)$. For the SVs, unlike the MMs, point $| g| =| \gamma | =1$ is not a critical one, as norms of the two SV's components are strongly different. However, figures 5(c), (d) demonstrates that there is a critical value ${\gamma }_{\mathrm{cr}}$ of γ, with the widths diverging and the norm share of the vortex component, ${F}_{-}$, approaching 0.5 at $\gamma \to {\gamma }_{\mathrm{cr}}$. In the same limit, the chemical potential, μ, approaches the above-mentioned limit value, $\mu =-1/2$ (see a typical example in figure 5(e) for N = 5), i.e., at $\gamma ={\gamma }_{\mathrm{cr}}$ the SV family displayed in figure 5 degenerates into an infinitely extended state with an infinitesimal amplitude. Obviously, ${\gamma }_{\mathrm{cr}}$ is a function of N, with ${\gamma }_{\mathrm{cr}}\to 0$ in the linear limit, i.e., at $N\to 0$. Figure 5(f) displays ${\gamma }_{\mathrm{cr}}$ in the interval of $3\lt N\lt {N}_{\mathrm{Townes}}$, the maximum value being ${\gamma }_{\mathrm{cr}}(N\to {N}_{\mathrm{Townes}})\approx 2.15$.

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 self-trapped 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}_{\mathrm{thr}}$, see equation (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 [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, figure 6 in [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 [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–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 ³⁹K atoms.

According to [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 equations (4) and (3) in [33]. However, it is shown in the same work that, in the case of

$$\begin{eqnarray}&&{n}_{2{\rm{D}}}\ll {a}^{-2},\end{eqnarray} \tag{ 15 }$$

where ${n}_{2{\rm{D}}}$ 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 equations (5) and (6) in [33]. A straightforward estimate demonstrates that inequality (15) amounts to the following condition for a characteristic size ${l}_{2{\rm{D}}}$ of a localized 2D pattern:

$$\begin{eqnarray}&&{l}_{2{\rm{D}}}^{2}\gg {{aa}}_{\perp },\end{eqnarray} \tag{ 16 }$$

where ${a}_{\perp }$ 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 [33] (it effectively represents zero-point energy of the Bogoliubov modes on top of the mean-field state):

$$\begin{eqnarray}&&{{ \mathcal E }}_{2{\rm{D}}}=\displaystyle \frac{1}{2}g{({| {\phi }_{+}| }^{2}-{| {\phi }_{-}| }^{2})}^{2}+\displaystyle \frac{{g}^{2}}{8\pi }{({| {\phi }_{+}| }^{2}+{| {\phi }_{-}| }^{2})}^{2}\mathrm{ln}\left(\displaystyle \frac{{| {\phi }_{+}| }^{2}+{| {\phi }_{-}| }^{2}}{\sqrt{e}{n}_{0}}\right),\end{eqnarray} \tag{ 17 }$$

where ${n}_{0}\equiv | {\phi }_{+}^{f}{| }^{2}+| {\phi }_{-}^{f}{| }^{2}$ is the density of the symmetric ($| {\phi }_{+}^{f}{| }^{2}=| {\phi }_{-}^{f}{| }^{2}$) 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

$$\begin{eqnarray}{\rm{i}}\displaystyle \frac{\partial {\phi }_{+}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{+}+\lambda \hat{D}{\phi }_{-}+g(| {\phi }_{+}{| }^{2}-| {\phi }_{-}{| }^{2}){\phi }_{+}+\displaystyle \frac{{g}^{2}}{4\pi }(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\rm{ln}}(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\phi }_{+},\\ {\rm{i}}\displaystyle \frac{\partial {\phi }_{-}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{-}-\lambda {\hat{D}}^{* }{\phi }_{+}+g(| {\phi }_{-}{| }^{2}-| {\phi }_{+}{| }^{2}){\phi }_{-}+\displaystyle \frac{{g}^{2}}{4\pi }(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\rm{ln}}(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\phi }_{-},\end{eqnarray} \tag{ 18 }$$

where the remaining scaling invariance is used to fix ${n}_{0}\equiv 1$.

Numerical solutions of equation (18) of the MM type, which are relevant here because of the opposite signs of the self- and cross-interactions in the terms $\sim g$, were produced by means of the same techniques as used above for the solution of equation (1). Thus, the controlled parameters in this section are a set of $(g,N,\lambda )$.

Numerical solution of equations (18) produces anisotropic (elongated) QD density profiles when $\lambda \ne 0$. The size and area of the QDs are much larger than for the solitons found in above sections. In accordance with equation (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 figure 6. When λ is smaller enough, the QD profile is nearly isotropic, see an example in the last column of figure 6. In the fourth row of figure 6, one can see that the total-density profile, $n({\bf{r}})=| {\phi }_{+}({\bf{r}}){| }^{2}+| {\phi }_{-}({\bf{r}}){| }^{2}$, naturally has the oval shape similar to that of each component.

Zoom In Zoom Out Reset image size

For the comparison of characteristics of the QD with different sets of control parameters, we plot the patterns of $n({\bf{r}})$ corresponding to different $(g,N,\lambda )$ sets in figure 7, choosing the vertical orientation for all the QDs. The figure shows that the area and the anisotropy of the $n({\bf{r}})$ patterns strongly depend on the parameters.

Zoom In Zoom Out Reset image size

In order to further quantify the dependence of the QDs on controlled parameters, we define average horizontal and vertical widths for $n({\bf{r}})$ as

$$\begin{eqnarray*}&&\{{W}_{x},{W}_{y}\}={\left(\displaystyle \frac{1}{N}\int \int \{{x}^{2},{y}^{2}\}n({\bf{r}}){\rm{d}}{x}{\rm{d}}{y}\right)}^{1/2}.\end{eqnarray*}$$

The degree of the anisotropy and the area of the QD can be defined by ${W}_{x,y}$ as

$$\begin{eqnarray}&&\mathrm{Ans}={W}_{y}/{W}_{x},\quad S=\pi {W}_{x}\times {W}_{y}.\end{eqnarray} \tag{ 19 }$$

These two characteristics, as well as the chemical potential of the QD, are displayed, as functions of g, N and λ in figure 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.

Zoom In Zoom Out Reset image size

It is also relevant to stress that, unlike the first and second row of figure 8, in the last row of figure 8 the dependence of the chemical potential, μ, on parameters $(g,N,\lambda )$ is not monotonous. In particular, the non-monotonous dependence $\mu (N)$ implies that the QD family does not satisfy the VK criterion, ${\rm{d}}\mu /{\rm{d}}{N}\lt 0$, although the entire QD family, including the segment with ${\rm{d}}\mu /{\rm{d}}{N}\gt 0$, is found to be completely stable. The violation of the VK criterion as the necessary stability condition [57] 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.

As shown in [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, equation (18) is rewritten as

$$\begin{eqnarray}{\rm{i}}\displaystyle \frac{\partial {\phi }_{+}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{+}+({\lambda }_{{\rm{R}}}\hat{D}-{\rm{i}}{\lambda }_{{\rm{D}}}{\hat{D}}^{* }){\phi }_{-}+g(| {\phi }_{+}{| }^{2}-| {\phi }_{-}{| }^{2}){\phi }_{+}+\displaystyle \frac{{g}^{2}}{4\pi }(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\rm{ln}}(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\phi }_{+},\\ {\rm{i}}\displaystyle \frac{\partial {\phi }_{-}}{\partial t} & = & -\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}{\phi }_{-}-({\lambda }_{{\rm{R}}}{\hat{D}}^{* }+{\rm{i}}{\lambda }_{{\rm{D}}}\hat{D}){\phi }_{+}+g(| {\phi }_{-}{| }^{2}-| {\phi }_{+}{| }^{2}){\phi }_{-}+\displaystyle \frac{{g}^{2}}{4\pi }(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\rm{ln}}(| {\phi }_{+}{| }^{2}+| {\phi }_{-}{| }^{2}){\phi }_{-},\end{eqnarray} \tag{ 20 }$$

where ${\lambda }_{{\rm{R}}}$ and ${\lambda }_{{\rm{D}}}$ are strengths of the Rashba and Dresselhaus SOC. Note that the Dresselhaus-only system, with ${\lambda }_{{\rm{R}}}=0$ and ${\lambda }_{{\rm{D}}}\ne 0$, is made tantamount to its Rashba-only counterpart, with ${\lambda }_{{\rm{R}}}\ne 0$ and ${\lambda }_{{\rm{D}}}=0$, by swapping $x\leftrightarrow y$ [59]. On the other hand, the mixed form of the SOC operators in equation (20) obviously breaks the rotational symmetry defined above by equation (3).

It was found in [16] that, for fixed ${\lambda }_{{\rm{R}}}$, there are critical values, ${({\lambda }_{D}^{{\rm{cr}}})}_{\mathrm{SV}}\lt {({\lambda }_{D}^{{\rm{cr}}})}_{\mathrm{MM}}\lt {\lambda }_{{\rm{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 ${\lambda }_{{\rm{D}}}={\lambda }_{{\rm{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 equation (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 ${\lambda }_{{\rm{R}}}$ and ${\lambda }_{{\rm{R}}}$. To illustrate these findings, we currently fix ${\lambda }_{{\rm{R}}}=1$ and vary ${\lambda }_{{\rm{D}}}$ in the interval of $[0,1]$. Typical examples of the the respective numerical solutions for QD with different values of ${\lambda }_{{\rm{D}}}$ and fixed values of $(g,N)$ are shown in figures 9(a1)–(a4). In particular, the orientation of the QDs tends to be diagonal. In the limit of ${\lambda }_{{\rm{D}}}={\lambda }_{{\rm{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({\bf{r}})$, as

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{1}{N}\int \int | n(x,y)-n(-y,x)| {\rm{d}}{x}{\rm{d}}{y},\end{eqnarray} \tag{ 21 }$$

where $n(-y,x)$ is the $90^\circ$ counterclockwise rotation of $n(x,y).$ Figure 9(b) shows Δ as a function of ${\lambda }_{{\rm{D}}}$ for fixed values of other parameters, $(g,N,{\lambda }_{{\rm{R}}})=(1,2,20)$. In figure 9, the QD profile become practically isotropic at ${\lambda }_{{\rm{D}}}\gt 0.6$.

Zoom In Zoom Out Reset image size

Lastly, direct simulations of the QD evolution in the framework of equation (20) demonstrate that asymmetric QDs are stable (see figure 9(c)), while effectively isotropic ones tend to be unstable against spontaneous drift in the diagonal direction, as shown in figure 9(d).

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 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\lt {N}_{\mathrm{Townes}}$, where ${N}_{\mathrm{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 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 VK 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 [20]. It may be relevant to consider 3D system with the competing signs of the SPM and XPM interactions, as well as MMs in 3D stabilized by the respective LHY terms.

We appreciate valuable discussions with G E Astrakharchik, H Sakaguchi, E Ya Sherman, Zheng Wei, and Pang Wei. This work was supported by NNSFC (China) through Grant No. 11575063, 61471123, 61575041, by the Natural Science Foundation of Guangdong Province, through Grant No. 2015A030313639, by the joint program in physics between NSF and Binational (US-Israel) Science Foundation through project No. 2015616, and by the Israel Science Foundation through Grant No. 1286/17. BAM appreciates a foreign-expert grant of the Guangdong province (China).