Two-dimensional dipolar gap solitons in free space with spin-orbit coupling

TL;DR

This paper demonstrates the creation and stability of two-dimensional dipolar gap solitons in free space with spin-orbit coupling in dipolar Bose-Einstein condensates, revealing their properties, stability conditions, and dynamics.

Contribution

It introduces stable 2D gap solitons in dipolar spinor BECs with SOC, including analytical and numerical analysis of their stability, shape, and mobility, which was not previously achieved.

Findings

Stable 2D gap solitons with vorticities 0 and 1 are found.

Dipolar interactions enable the formation of isotropic and anisotropic solitons.

Solitons exhibit stable embedded states within spectral bands.

Abstract

We present gap solitons (GSs) that can be created in free nearly two-dimensional (2D) space in dipolar spinor Bose-Einstein condensates with the spin-orbit coupling (SOC), subject to tight confinement, with size $a_{\perp }$, in the third direction. For quasi-2D patterns, with lateral sizes $l\gg a_{\perp }$, the kinetic-energy terms in the respective spinor Gross-Pitaevskii equations may be neglected in comparison with SOC. This gives rise to a bandgap in the system's spectrum, in the presence of the Zeeman splitting between the spinor components. While the present system with contact interactions does not produce 2D solitons, stable gap solitons(GSs), with vorticities $0$ and $1$ in the two components, are found, in quasi-analytical and numerical forms, under the action of dipole-dipole interaction (DDI). Namely, isotropic and anisotropic 2D GSs are obtained when the dipoles are polarized, respectively, perpendicular or parallel to the 2D plane. The GS families extend, as \textit{embedded solitons} (ESs), into spectral bands, a part of the ES branch being stable for isotropic solitons. The GSs remain stable if the competing contact interaction, with the sign opposite to that of the DDI, is included, while the addition of the contact term with the same sign destabilizes the GSs, at first replacing them by breathers, and eventually leading to destruction of the solitons. Mobility and collision of the GSs are studied too, revealing negative and positive effective masses of the isotropic and anisotropic solitons, respectively.

Equations

$$i{\partial }_t{\varphi }_{+}=-{\left(2m\right)}^{-1}{\nabla }^2{\varphi }_{+}+\lambda \left({\partial }_x-i{\partial }_y\right){\varphi }_{-}-\Omega {\varphi }_{+}$$

$$+\left(g|{\varphi }_{+}{|}^2+\tilde{g}|{\varphi }_{-}{|}^2\right){\varphi }_{+}+\kappa {\varphi }_{+}\int R(\mathbf{r}-{\mathbf{r}}^{\prime })(|{\varphi }_{+}({\mathbf{r}}^{\prime }){|}^2+|{\varphi }_{-}({\mathbf{r}}^{\prime }){|}^2)d{\mathbf{r}}^{\prime },$$

$$i{\partial }_t{\varphi }_{-}=-{\left(2m\right)}^{-1}{\nabla }^2{\varphi }_{\pm }-\lambda \left({\partial }_x+i{\partial }_y\right){\varphi }_{+}+\Omega {\varphi }_{-}$$

$$+\left(g|{\varphi }_{-}{|}^2+\tilde{g}|{\varphi }_{+}{|}^2\right){\varphi }_{-}+\kappa {\varphi }_{-}\int R(\mathbf{r}-{\mathbf{r}}^{\prime })(|{\varphi }_{+}({\mathbf{r}}^{\prime }){|}^2+|{\varphi }_{-}({\mathbf{r}}^{\prime }){|}^2)d{\mathbf{r}}^{\prime },$$

$$R_{\mathrm{iso}}(\mathbf{r}-{\mathbf{r}}^{\prime })=1/[{ϵ}^2+(\mathbf{r}-{\mathbf{r}}^{\prime }{)}^2{]}^{3/2},$$

$$R_{\mathrm{aniso}}(\mathbf{r}-{\mathbf{r}}^{\prime })=(1-3{\cos }^2\Theta )/[{ϵ}^2+(\mathbf{r}-{\mathbf{r}}^{\prime }{)}^2{]}^{3/2},$$

$$i{\partial }_t{\varphi }_{+}=\lambda \left({\partial }_x-i{\partial }_y\right){\varphi }_{-}-\Omega {\varphi }_{+}+\left(g|{\varphi }_{+}{|}^2+\tilde{g}|{\varphi }_{-}{|}^2\right){\varphi }_{+}$$

$$+\kappa {\varphi }_{+}\int R(\mathbf{r}-{\mathbf{r}}^{\prime })\left[|{\varphi }_{+}({\mathbf{r}}^{\prime }){|}^2+|{\varphi }_{-}({\mathbf{r}}^{\prime }){|}^2\right]d{\mathbf{r}}^{\prime },$$

$$i{\partial }_t{\varphi }_{-}=-\lambda \left({\partial }_x+i{\partial }_y\right){\varphi }_{+}+\Omega {\varphi }_{-}+\left(g|{\varphi }_{-}{|}^2+\tilde{g}|{\varphi }_{+}{|}^2\right){\varphi }_{-}$$

$$+\kappa {\varphi }_{-}\int R(\mathbf{r}-{\mathbf{r}}^{\prime })\left[|{\varphi }_{+}({\mathbf{r}}^{\prime }){|}^2+|{\varphi }_{-}({\mathbf{r}}^{\prime }){|}^2\right]d{\mathbf{r}}^{\prime },$$

$${\omega }^2={\Omega }^2+{\lambda }^2(p^2+q^2),$$

$$\mu u_{+}=\lambda \left({\partial }_x-i{\partial }_y\right)u_{-}-\Omega u_{+}+\left(g|{\varphi }_{+}{|}^2+\tilde{g}|{\varphi }_{-}{|}^2\right)u_{+}$$

$$+\kappa u_{+}\int R(\mathbf{r}-{\mathbf{r}}^{\prime })\left[|u_{+}({\mathbf{r}}^{\prime }){|}^2+|u_{-}({\mathbf{r}}^{\prime }){|}^2\right]d{\mathbf{r}}^{\prime }.$$

$$\mu u_{-}=-\lambda \left({\partial }_x+i{\partial }_y\right)u_{+}+\Omega u_{-}+\left(g|{\varphi }_{-}{|}^2+\tilde{g}|{\varphi }_{+}{|}^2\right)u_{-}$$

$$+\kappa u_{-}\int R(\mathbf{r}-{\mathbf{r}}^{\prime })\left[|u_{+}({\mathbf{r}}^{\prime }){|}^2+|u_{-}({\mathbf{r}}^{\prime }){|}^2\right]d{\mathbf{r}}^{\prime }.$$

$$N=N_{+}+N_{-}\equiv \int \left(|u_{+}(\mathbf{r}){|}^2+|u_{-}(\mathbf{r}){|}^2\right)d\mathbf{r}.$$

$$u_{+}=U_{+}(r),u_{-}=U_{-}(r)e^{i\theta },$$

$$\mu U_{+}=\lambda \left(\frac{d{U}_{-}}{dr}+\frac{1}{r}{U}_{-}\right)-\Omega U_{+}+\left(g{U}_{+}^2+\tilde{g}{U}_{-}^2\right)U_{+}$$

$$+\kappa U_{+}{\int }_0^{\infty }\varrho (r,r^{\prime })\left[U_{+}^2(r{}^{\prime })+U_{-}^2(r{}^{\prime })\right]r^{\prime }dr{}^{\prime }.$$

$$\mu U_{-}=-\lambda \frac{d{U}_{+}}{dr}+\Omega U_{-}+\left(g{U}_{-}^2+\tilde{g}{U}_{+}^2\right)U_{-}$$

$$+\kappa U_{-}\int \varrho (r,r^{\prime })\left[U_{+}^2(r{}^{\prime })+U_{-}^2(r{}^{\prime })\right]r^{\prime }dr{}^{\prime },$$

$$\varrho (r,r^{\prime })$$

$$k$$

$$\mu =\mp \left(\Omega -\delta \mu \right),0<\delta \mu \ll \Omega ,$$

$$u_{\mp }\approx \left(\lambda /2\Omega \right)\left({\partial }_x\pm i{\partial }_y\right)u_{\pm }.$$

$$\delta \mu \cdot u_{\pm }=\frac{{\lambda }^2}{2\Omega }{\nabla }^2{u}_{\pm }\pm g|u_{\pm }{|}^2{u}_{\pm }\pm \kappa u_{\pm }(\mathbf{r})\int R(\mathbf{r}-{\mathbf{r}}^{\prime })|u_{\pm }({\mathbf{r}}^{\prime }){|}^{2}d{\mathbf{r}}^{\prime },$$

$$\mp i\frac{\partial }{\partial t}{\tilde{u}}_{\pm }=\frac{{\lambda }^2}{2\Omega }{\nabla }_{\pm }^2\tilde{u}\pm g|{\tilde{u}}_{\pm }{|}_{\pm }^2\tilde{u}\pm \kappa {\tilde{u}}_{\pm }(\mathbf{r})\int R(\mathbf{r}-{\mathbf{r}}^{\prime })|{\tilde{u}}_{\pm }({\mathbf{r}}^{\prime }){|}^{2}d{\mathbf{r}}^{\prime },$$

$$N_{\mathrm{at}}\sim 10^{3}N,{\left(x,y\right)}_{\mathrm{phys}}\sim (x,y)\times 20\mu \mathrm{m}.$$

$$\mathbf{P}=i\int \left[\left(\nabla {\varphi }_{+}^{\ast }\right){\varphi }_{+}+\left(\nabla {\varphi }_{-}^{\ast }\right){\varphi }_{-}\right]d\mathbf{r},$$

$${\varphi }_{\pm }(x,y,t=0)=u_{\pm }(x,y)\left\{e^{i\eta x},e^{i\eta y}\right\},$$

$$g<g_{\mathrm{Townes}}^{(\mathrm{iso})}=N_{\mathrm{Townes}}/\left(\Omega N\right),g>g_{\mathrm{Townes}}^{(\mathrm{aniso})}=-N_{\mathrm{Townes}}/\left(\Omega N\right)$$

$${ϵ}_{\mathrm{dd}}\equiv g_{\mathrm{eff}}^{(\mathrm{iso},\mathrm{aniso})}/g,$$

$$\kappa u_{+}(\mathbf{r}=0)\int R_{\mathrm{iso}}(0-{\mathbf{r}}^{\prime })|u_{+}({\mathbf{r}}^{\prime }){|}^{2}d{\mathbf{r}}^{\prime }\equiv g_{\mathrm{eff}}^{(\mathrm{iso})}|u_{+}(\mathbf{r}=0){|}^2{u}_{+}(\mathrm{r}=0),$$

$$\kappa u_{-}(\mathbf{r}=0)\int R_{\mathrm{aniso}}(0-{\mathbf{r}}^{\prime })|u_{-}({\mathbf{r}}^{\prime }){|}^{2}d{\mathbf{r}}^{\prime }\equiv g_{\mathrm{eff}}^{(\mathrm{aniso})}|u_{-}(\mathbf{r}=0){|}^2{u}_{-}(\mathrm{r}=0).$$

Full text

I Introduction

Gap solitons (GSs) are usually defined as self-trapped modes existing in spectral bandgaps of periodic potentials. GSs have been predicted and observed in diverse optical media, such as Bragg gratings Krug , waveguide arrays arrays , and photonic crystals Segev (see also reviews DS-review0 ; DS-review1 ; DS-review2 ), as well as in Bose-Einstein condensates (BECs) trapped in optical lattices BEC ; GS-BEC , and in a plasmonic medium including a lattice potential plasmon . The dynamics of spatial GSs is usually modeled by the nonlinear Schrödinger/Gross-Pitaevskii equations (NLSEs/GPEs) with periodic potentials. GSs in fiber gratings are described by nonlinear coupled-mode equations for counterpropagating waves DS-review0 . In all cases, the underlying spatially periodic structures play a key role, producing spectral bandgaps in which GSs can be created.

A challenging problem in optics and BEC is the creation of stable fundamental and vortical bright solitons in two- and three-dimensional (2D and 3D) geometry old-review ; Dumitru ; recent-review . New possibilities are offered by spin-orbit coupling (SOC) in spinor BEC SOC1 and its counterparts in optics Barcelona ; Sakaguchi2 . In particular, the interplay of the linear SOC with the cubic attractive nonlinearity opens a way for creating 2D ground-state Sakaguchi ; Barcelona and 3D metastable HP solitons in free space, which was previously deemed impossible (1D Kartashov1D and 2D Kartashov2D GSs supported by a combination of lattice potentials and SOC were predicted too).

We aim to demonstrate that SOC offers another unexpected possibility, to create stable 2D GSs in free space, without the use of any periodic potential. The model is formulated in Section II, where we consider the 2D condensate under the action of strong SOC, which makes it possible to neglect the kinetic-energy terms in the corresponding system of coupled GPEs, thus reducing them to a first-order system of coupled-mode equations. The gap in the system’s spectrum is generated by the Zeeman splitting (ZS), which is an essential ingredient of SOC settings Zeeman . However, numerical results demonstrate that the usual contact nonlinearity of any sign fails to build solitons in this 2D system. Our main result is prediction of families of stable isotropic and anisotropic 2D GSs under the action of dipole-dipole interaction (DDI). We note that the realization of SOC in the condensate of dipolar chromium atoms was proposed in Ref. Cr , and elaborated theoretically in subsequent works Xunda .

Using a combination of an analytical approximation, which is available close to edges of the bandgap (similar to the approximation recently developed in Ref. Sherman ), and numerical methods, in Section III we construct families of stationary solutions for isotropic and anisotropic GSs in the condensates composed of dipoles oriented, respectively, perpendicular or parallel to the system’s plane. The isotropic and anisotropic families extend, severally, as embedded solitons (ESs) Jianke1999 across the top and bottom edge of the bandgap into the adjacent Bloch band. Stability of the GSs and ESs is explored by means of systematic direct simulations. Mobility and collisions of the GSs are studied in Section IV, showing that the effective mass is positive for the anisotropic 2D solitons, and negative for the isotropic ones (negative mass is a known property of GSs negative-mass ). In Section V, we consider an extended system, which includes both the dipole-dipole and contact interactions. The conclusion is that the solitons persist as stable modes under the action of competing interactions with opposite signs. If the signs are identical, the addition of the contact interaction tends to destabilize a soliton, replacing it by a breather. Eventually, all self-trapped modes are destroyed. The paper is concluded by Section VI.

II The model

In the scaled form (the corresponding estimates of physical parameters are given below), the coupled GPEs for two components of the spinor wave function, $\phi_{\pm}$, carrying the same magnetic moment and coupled by the SOC of the Rashba type Rashba , are SOCsoliton1 ; SOCsoliton2 ; Sakaguchi ; Sherman ; Xunda ; Guihua2017 :

[TABLE]

[TABLE]

where $m$ is the atomic mass, $g$ and $\tilde{g}$ are strengths of the contact self- and cross-interactions (usually, these strengths are nearly equal in mixtures of different states of the same atomic species Ho ), while $\lambda$, $\Omega>0$, and $\kappa>0$ represent SOC, ZS, and DDI, respectively Cr . Note that the sign of $\kappa$ may be altered by means of a rotating magnetic field Giovanazzi , and ZS may be replaced by the Stark - Lo Surdo splitting in dc electric field.

SOC of the Rashba type is adopted here as it helps to build 2D solitons, while SOC terms of the Dresselhaus type tend to cause delocalization Sherman . Nevertheless, Eqs. (16) and (17), derived below in the limit case when one component is much larger than the other, take a universal form for any combination of the Rashba and Dresselhaus terms. Although, strictly speaking, those asymptotic equations are valid only in small vicinities of the bottom and top edges of the spectral bandgap, see Eq. (14) below, the results reported in the next section [see Fig. 1(b)] clearly demonstrate that the predictions produced by Eqs. (16) and (17) are quite accurate in the entire bandgap, hence the universality implied by those equations remains approximately valid for the full GS families.

In the isotropic setting, with dipoles oriented perpendicular to the $\left(x,y\right)$ plane, the respective repulsive DDI kernel is

[TABLE]

where cutoff $\epsilon$ is determined by the confinement in the transverse dimension Santos ; epsilon ; quadr , while effects of the attractive DDI in the transverse direction are suppressed by the tight confinement, whose strength (trapping frequency) is much larger than chemical potential produced by the effective 2D GPE Tikhonenkov . If the dipoles are polarized parallel to the $\left(x,y\right)$ plane Tikhonenkov , the DDI is anisotropic, with

[TABLE]

where $\Theta$ is the angle between the polarization direction and ${(\mathbf{r}-\mathbf{r^{\prime}})}$. The simplest approximation for the cutoff, adopted here, is sufficient, as the detailed analysis demonstrates that the exact form gives rise to practically the same result Santos ; epsilon ; quadr .

The SOC coefficient relevant to experimental settings with transverse-confinement size $a_{\perp}$ is estimated, in physical units, as $\lambda\gtrsim\hbar^{2}/\left(ma_{\perp}\right)$ Drummond-Pu . It follows from here that the kinetic-energy terms in Eq. (1) may be neglected for all quasi-2D patterns, with lateral sizes $l\gg a_{\perp}$, reducing it to the coupled-mode equations,

[TABLE]

[TABLE]

cf. Ref. Sakaguchi2 in 1D. In a different context, a 2D SOC system with a deep optical-lattice potential and negligible kinetic energy was recently introduced in Ref. flatband .

The linear part of Eqs. (4) and (5) is essentially the same as in the 2D nonlinear Dirac equation (NDE) for bosons in honeycomb lattices Haddad1 , which gives rise to stable states in the presence of an external trap Haddad2 . In the absence of the trap, stable 2D solitons were produced by NDE with the sign-indefinite contact nonlinear terms, $\pm(|\phi_{+}|^{2}-|\phi_{-}|^{2})\phi_{\pm}$, in its two components Kevrekidis . However, this case is not relevant to the BEC spinor wave function, whose components represent atomic states with nearly equal intra- and inter-component scattering lengths Ho . We have checked that Eqs. (4) and (5) with the corresponding contact nonlinearity, $\sim\left(|\phi_{+}|^{2}+|\phi_{-}|^{2}\right)\phi_{\pm}$, fail to create 2D solitons, unlike the system which includes the kinetic energy Sakaguchi . This is explained by the fact that the gradient energy corresponding to the Dirac operator in Eqs. (4) and (5), unlike its Schrödinger counterpart, cannot balance the local cubic nonlinearity. However, it is shown below that GSs are readily produced by the balance of the SOC with the DDI, as well as with a combination of the DDI and contact interaction with opposite signs. The statements concerning the balance are corroborated by analytical results reported below for GSs located close to bandgap edges.

Linearizing Eqs. (4) and (5) for $\phi_{\pm}\sim\exp\left(ipx+iqy-i\omega t\right)$, we derive the dispersion relation,

[TABLE]

with the free-space bandgap of width $2\Omega$, as shown in Fig. 1(a) for $\lambda=1$, $\Omega=10$ (by means of scaling, these values of the ZS and SOC strengths are fixed in the present work). The bandgap will close if the small kinetic-energy terms are kept in Eq. (1), which may cause transformation of the GSs into usual solitons at times much larger than experimentally relevant time scales.

GS solutions to Eqs. (4) and (5) with chemical potential $\mu$ are looked for as $\phi_{\pm}=e^{-i\mu t}u_{\pm}(x,y)$, with $u_{\pm}$ obeying equations

[TABLE]

[TABLE]

The GS is characterized by the norm, which is proportional to the total number of atoms in the binary BEC,

[TABLE]

Furthermore, Eqs. (7) and (8) with the isotropic DDI kernel given by Eq. (2) admit solutions with an exact structure of semi-vortices Sakaguchi ; HP ; Sherman , i.e., isotropic complexes with zero vorticity $S=0$ in component $u_{+}$, and vorticity $S=1$ in $u$:

[TABLE]

where $\left(r,\theta\right)$ are the polar coordinates in the $\left(x,y\right)$ plane, and radial functions $U_{\pm}(r)$ are real. Indeed, the substitution of the semi-vortex ansatz (10) in Eqs. (7) and (8) yields, after performing the angular integration in the DDI terms, to the following equations which include solely the radial coordinate:

[TABLE]

[TABLE]

where the effective radial kernel is

[TABLE]

Here, $E(k)$ is the standard complete elliptic integral of the second kind with modulus $k$.

Numerical results presented in the next section, see Figs. 2(a,b,c) and 3(a1,a2,a3), confirm that Eqs. (7) and (8) with the isotropic kernel indeed give rise to the two-component solitons whose structure precisely conforms to ansatz (10). Furthermore, it is demonstrated below that the anisotropic kernel gives rise, as a matter of fact, to deformed patterns of the same semi-vortex type, with $S=0$ in component $u_{-}$ and $S=-1$ in $u_{+}$, see Figs. 2(d,e,f) and 3(b1,b2,b3) below.

III Families of gap solitons and embedded solitons

III.1 Isotropic and anisotropic solitons near edges of the bandgap

Close to edges of the bandgap of dispersion relation (6), i.e., at

[TABLE]

the two-component problem can be reduced to one of those previously solved for the single component in the semi-infinite gap Pedri ; Tikhonenkov . Close to the bottom (top) edge, Eq. (8) for $u_{-}$ ($u_{+}$) makes it possible to eliminate this component in favor of $u_{+}$ ($u_{-}$):

[TABLE]

In particular, for the isotropic semi-vortex represented by ansatz (10), Eq. (15) amounts to $U_{-}(r)\approx\left(\lambda/2\Omega\right)dU_{+}/dr$.

Substituting relation (15) in Eq. (8) for $u_{\pm}$, one arrives at an equation for the single component:

[TABLE]

where $\nabla^{2}=(\partial_{x}+i\partial_{y})(\partial_{x}-i\partial_{y})\equiv\partial_{x}^{2}+\partial_{y}^{2}$ appears as a square of the SOC operator from Eqs. (4) and (5). It is easy to see that the same $\nabla^{2}$ is produced by squaring the SOC operator which presents a general combination of the Rashba and Dresselhaus terms.

Along with Eq. (16), it is relevant to consider its time-dependent version,

[TABLE]

where $\tilde{u}_{\pm}\left(x,y,t\right)\equiv\exp\left(\pm i\Omega t\right)~{}u\left(x,y,t\right)$. In particular, Eq. (17) is used below to test stability of various solitons generated by Eq. (16).

If the contact nonlinear term $\sim g$ is present in Eq. (16), while the DDI is absent ($\kappa=0$), this equation produces no solitons in the case of the effective self-defocusing, $\pm g<0$ (recall we define $\Omega$ to be positive). In the case of the local self-focusing,

$\pm g>0$, Eq. (16) produces unstable Townes solitons Townes . This argument suggests that the contact interaction of either sign, in the absence of the DDI, cannot support stable solitons in the present system. As mentioned above, this expectation is fully corroborated by numerical results (not shown here in detail).

Proceeding to the opposite case, when the DDI is present, while the contact interaction is absent ($g=0$), we note that, near the bottom edge of the bandgap, Eq. (16) for $u_{+}$, with $\kappa>0$ and the isotropic kernel, $R=R_{\mathrm{iso}}(\mathbf{r}-\mathbf{r^{\prime}})$, gives rise to the ground state in the form an axisymmetric bright soliton with zero vorticity Pedri , while the respective smaller component $u_{-}$, produced by Eq. (15), features vortical structure $\sim re^{i\theta}$ (at small $r$), conforming to the semi-vortex structure of the isotropic GS, as given by Eq. (10). In fact, in this case, the 2D integration in Eq. (16) can be reduced to the radial-only integration, with kernel $R_{\mathrm{iso}}(\mathbf{r}-\mathbf{r^{\prime}})$ substituted by the effective radial one, as done above, in the general form, in Eq. (13).

On the other hand, near the top edge of the bandgap, Eq. (16) for $u_{-}$, with $\kappa>0$ and the anisotropic kernel taken as per Eq. (3), gives rise to the ground state in the form of the 2D anisotropic bright soliton, as previously demonstrated in Ref. Tikhonenkov , while Eq. (15) produces a smaller component, $u_{+}$, in the form of an *anisotropic vortex. *Overall, the anisotropic GS constructed in this form seems as a deformed semi-vortex, as corroborated by numerically exact results displayed below in Fig. 3(b1,b2,b3).

III.2 Numerical findings for generic two-component solitons

The above results predict that the isotropic DDI for dipoles polarized perpendicular to the $\left(x,y\right)$ plane, and the anisotropic DDI for the in-plane polarization, support, severally, stable isotropic GSs [with the semi-vortex structure, as per Eq. (10)] and anisotropic 2D GSs, with $\mu$ taken, respectively, close to the bottom or top edge of the bandgap. These predictions have been corroborated by numerical solutions of Eqs. (7) and (8). The numerical results are collected in Fig. 1(b1,b2) for the system which does not include the contact interaction [$g,\tilde{g}=0$ in Eqs. (7) and (8)]. The results clearly demonstrate that the quasi-analytical approximations remain valid not only close to the edges, but actually across the entire bandgap, and extend, as ESs (see further details below), into the bands.

Numerical solutions of Eq. (8) were produced by means of the squared-operator method Jianke2007 . The scaling invariances of Eqs. (4) and (8) were used to fix $\Omega=10$, $\lambda=1$, and $\kappa=0.1$. Generic results were produced fixing the regularization parameter as $\epsilon=0.5$ (with other reasonable values of $\epsilon$, similar results have been obtained), while the total norm, $N$, was varied as an essential control parameter. The stability of the GS families was identified by means of systematic simulations of the perturbed evolution (the distinction between stable and unstable states could be easily detected, as numerical truncation errors were sufficient to trigger the growth of the instability, if any). As concerns the asymptotic equation (16), its solutions of were produced by applying the imaginary-time method ITM to its time-dependent version given by Eq. (17). Comparison of shapes of stable solitons, as found from Eqs. (7) and (8), and, on the other hand, from the simplified equations (16) and (17), is displayed in Fig. 2.

Branches of isotropic and anisotropic solitons are characterized by $\mu(N)$ dependences displayed in Fig. 1(b1,b2), along with the semi-analytical counterparts of these dependences. It has been thus found that the stable branch of the isotropic GSs extends, across the full bandgap, into the upper Bloch band abutting on the bandgap, as a family of ESs, which may exist, under certain conditions, in spectral bands Jianke1999 . In particular, a model supporting ESs in a 2D system was reported in Ref. 2D-ES . With the increase of $N$, the isotropic-ES branch loses its stability inside of the Bloch band, at $N_{\max}\approx 45$ ($\mu\approx 16.7$). At $N>N_{\max}$, the branch extends indefinitely into the band in an unstable form. On the contrary, the stability of the anisotropic GSs terminates still in the bandgap, at $N_{\max}\approx 28.2$ ($\mu\approx$ $-3.8$), the ES continuation of this branch being fully unstable. The robustness of the solitons in the present system is further attested to by the fact that unstable ones, both isotropic and anisotropic, do not suffer destruction, their vortex component keeping its vorticity: as shown in Fig. 4, unstable solitons commence spontaneous motion, instead of destruction, emitting small amounts of radiation from the vortex component. The evolution of the weakly unstable isotropic ESs (recall all the isotropic GS are stable) does not break their circular symmetry either.

As said above, isotropic GSs are built as semi-vortices [defined as per Eq. (10)], i.e., bound states of zero-vorticity and vortex components, as can be clearly seen in Figs. 3(a1-a3). The extension of the GS solutions in the ES form keeps the semi-vortex shape as well. Although no exact ansatz for a vortex structure is available for anisotropic GSs, Figs. 3(b1-b3) clearly demonstrate that the anisotropic GSs (and their ES extension) feature the shape of deformed semi-vortices. An essential peculiarity of the semi-vortices is that their vortex component carries a relatively small share of the total norm, in comparison with the zero-vorticity counterpart, as shown in Fig. 1(c) for the isotropic and anisotropic semi-vortices alike.

The stable isotropic solitons shown above, with vorticities $(S_{+},S_{-})=(0,1)$ in their large and small components, are fundamental states, as the system cannot produce any state with a simpler structure. However, it is possible to look for more complex modes (excited states). To generate them near the bottom edge of the bandgap, one can use, as a seed, isotropic vortex-solitons solutions of Eq. (16) for the large component, $u_{+}$, with vorticities $S_{+}=\pm 1$, which are known from the study of the single-component dipolar BEC Tikhon2 . Then, Eq. (15) generates vorticity $S_{-}=S_{+}+1$ in the small $u_{-}$ component. We have found that both species of the resulting composite modes, with $(S_{+},S_{-})=(1,2)$ and $(-1,0)$, are unstable against splitting into a pair of fragments (not shown here in detail). In fact, splitting is a common instability mode of vortex solitons old-review ; Dumitru ; recent-review ; Skryabin , although nonlocality may stabilize some of them Kiev ; Tikhon2 .

III.3 Physical parameters for the gap solitons in the Bose-Einstein

condensate

It is relevant to estimate actual parameters of the BEC solitons which can be created according to the results reported above in the scaled form. To this end, we translate the results into physical units corresponding to the experimental realization of SOC Zeeman ; Cr ; Xunda ; Drummond-Pu , taking values of the magnetic moment for 52Cr or 164Dy atoms, and the strength of the transverse trapping potential $\omega_{\perp}\sim 100$ Hz. We thus conclude that the stable quasi-2D solitons may be created with the number of atoms in the range from $N_{\mathrm{at}}\sim 10^{3}$ (near edges of the bandgap)   to $N_{\mathrm{at}}\sim 10^{4}$ (deeper in the bandgap), and physical lateral sizes $l_{\mathrm{phys}}\sim 10$ $\mathrm{\mu}$m. The corresponding relations between the physical quantities and scaled ones displayed in Figs. 1-7 is

[TABLE]

Further, the magnetic field necessary for inducing the appropriate ZS is estimated as $H\sim 0.1-1$ G. In addition to the magnetic realization, the use of a spinor BEC built of small molecules carrying electric dipole moments electric may be feasible too. Another possibility for the realization of the stable 2D solitons predicted above is to use the DDI between moments induced by an external dc magnetic or electric field induced .

IV Mobility and collisions of the isotropic and anisotropic solitons

Soliton mobility in the system under the consideration is a nontrivial issue, as Eqs. (4) and (5) are not Galilean invariant, although they conserve the total momentum,

[TABLE]

and their asymptotic version, which amounts to the single equation (17), remains Galilean invariant. The mobility was tested, in the framework of the system without the contact interactions ($g=\tilde{g}=0$), by applying kick $\eta$ to stable quiescent solitons in the $x$ or $y$ direction, which correspond to simulating Eqs. (4) and (5) with input

[TABLE]

the components of the respective momentum (19) being $P_{x,y}=N\eta$. An example of stable motion of a kicked isotropic GS (semi-vortex) is shown in Fig. 5(a). The velocity of the moving isotropic soliton, $V_{x}$, is displayed in Fig. 5(b) as a function of the kick’s magnitude, $\eta$. It is seen that the mass of kicked isotropic GSs is negative, as they move in the direction opposite to the kick. It is relevant to mention that the negative dynamical mass is a generic feature of GSs known in other settings, including 2D GSs HS ; plasmon .

The effective mass of the isotropic solitons in the $x$ and $y$ directions, defined as $M_{x,y}(N)=P_{x.y}/V_{x,y}$, is displayed in Fig. 5(c1) as a function of $N$ (strictly speaking, the motion makes the soliton slightly anisotropic). Naturally, the mass increases almost linearly with the norm. Close to the bottom edge of the bandgap, where $N$ is very small, the mass is isotropic, $M_{x}\approx M_{y}$, as the corresponding equation (16) is isotropic too. Farther from the bandgap edge, the mass features a weak anisotropy. On the other hand, the anisotropic solitons feature a positive mass, which also grows almost linearly with $N$, as shown in Fig. 5(c2). In the course of its motion, the kicked solitons show very weak transverse displacement, due to the “Magnus force” acting on the small vortical component Magnus . Additional high-accuracy simulations are needed to study the latter effect in detail.

It is natural too to study collisions between stable GSs set in motion by opposite kicks $\pm\eta$. The simulations reveal three generic outcomes of the collisions, for the isotropic and anisotropic solitons alike: merger into a single breather [Figs. 6(a1,a2)], elastic collision [Figs. 6(c1,c2)], and transformation into diffracting quasi-linear beams [Figs. 6(b1,b2)], at, small, large, and intermediate values of $\eta$, respectively.

V Effects of the contact interaction

The local mean-field nonlinearity, induced by interatomic collisions, is always present in the bosonic gas, therefore it is relevant to explore effects of the contact interaction on the soliton families obtained above in the absence of the contact terms in Eqs. (7) and (8). Here, we perform this analysis in the framework of the asymptotic equations (16) and (17), as they are sufficient to capture main effects produced by the local nonlinearity, as shown below.

First, Eq. (16), that includes the contact-interaction term $\sim g$, implies that, due to the possibility of the critical collapse Townes in the same 2D equation, the norm of wave function $u_{\pm}$ cannot exceed a critical value, which is determined by the scaled norm of the Townes’ soliton, $N_{\mathrm{Townes}}\approx 5.85$ Townes , viz., ${N_{\pm}<N_{\mathrm{Townes}}/}\left(\pm g\Omega\right)$. In other words, for given norm $N$, isotropic and anisotropic self-trapped modes exist, respectively, at

[TABLE]

(recall $\Omega=10$ is fixed in this paper).

Furthermore, inside the existence regions (21), systematic simulations of Eq. (17) have revealed an intrinsic stability boundary, $g=$ $g_{\mathrm{cr}}^{(\pm)}\approx\pm 0.25$, such that the isotropic and anisotropic stationary GSs are stable, respectively, at $g<$ $g_{\mathrm{cr}}^{(+)}$ and $g>g_{\mathrm{cr}}^{(-)}$, while in the remaining intervals, g_{\mathrm{cr}}^{(+)}<g<$${g_{\mathrm{Townes}}^{\mathrm{(iso)}}} and {g_{\mathrm{Townes}}^{\mathrm{(aniso)}}<g<}$$g_{\mathrm{cr}}^{(-)}, the GSs are unstable, spontaneously transforming into persistent breathers, as shown in Fig. 7. Thus, we conclude that the stationary GSs exist and remain completely stable when the arbitrarily strong contact interaction is* self-repulsive*, in terms of Eqs. (16) and (17) (which corresponds to $g<0$ and $g>0$ for the isotropic and anisotropic GSs, respectively), being compensated by the effectively attractive DDI. In the limit of very strong local self-attraction, the solitons become very broad, corresponding to $\delta\mu\rightarrow 0$ in terms of Eq. (16), i.e., $\mu\rightarrow-\Omega$ and $\mu\rightarrow+\Omega$, in terms of Figs. 7(c) and (d), respectively. Note an essential difference of these limits from those shown in Figs. 1(b1) and (b2): in the latter case, the soliton’s norm vanishes at the edges of the bandgap, while in the cases displayed in Figs. 7(c) and (d) the soliton branches keep the fixed finite norm, $N_{\pm}=2$.

To explicitly compare strengths of the competing contact interactions and DDI, we define the relative strength,

[TABLE]

where $g_{\mathrm{dd}}^{(\mathrm{iso,aniso})}$ are effective DDI coefficients for the isotropic and the anisotropic GSs, respectively, which are defined, for the definiteness’ sake, at centers of the solitons, as follows:

[TABLE]

[TABLE]

Here, we adopt conditions $|u_{+}(\mathbf{r})|^{2}\gg|u_{-}(\mathbf{r})|^{2}$ and $|u_{-}(\mathbf{r})|^{2}\gg|u_{+}(\mathbf{r})|^{2}$ for the isotropic and anisotropic GSs, respectively, as Figs. 7(c) and (d) clearly demonstrates that these conditions hold at the critical (most essential) points.

The so defined relative DDI strengths (22) are displayed, as functions of $g$, in Fig. 8, for isotropic and anisotropic GSs at $g>0$ and $g<0$, respectively. A natural conclusion is that the GSs suffer the destabilization when the DDI becomes too weak in comparison with the contact interaction.

On the contrary to the setting with the competing contact interaction and DDI, the interplay of the local and nonlocal nonlinear interactions with identical signs leads to the destabilization of the solitons in the present system, as is shown by Figs. 7(c,d). Overall, the present situation qualitatively resembles that reported in Ref. Minhang , in which a binary BEC represented a mixture of two atomic states coupled by a microwave field. In such a setting, nonlocal attraction between the components, mediated by the microwave field, was sufficient for the existence and stability of two-component solitons in the presence of arbitrarily strong local self- and cross- repulsion.

VI Conclusion

The objective of this work is to propose a setting for the creation of stable quasi-2D GSs (gap solitons) and ESs (embedded solitons) in free space. The system is based on the spinor dipolar BEC, whose components are coupled by the SOC (spin-orbit interaction). For quasi-2D states, the kinetic-energy terms in the spinor GPEs are negligible in comparison with SOC, which gives rise to simplified couple-mode equations, with the bandgap provided by the ZS (Zeeman splitting). Stable isotropic and anisotropic 2D solitons were thus found, in the quasi-analytical and numerical forms, for the dipoles polarized perpendicular and parallel to the system’s plane, respectively. Both families continue as ESs into adjacent spectral bands, the isotropic-ES branch being partly stable. Mobility and collisions of the solitons were studied too, concluding that the mass of the isotropic/anisotropic ones is negative/positive. Effects of contact interactions, added to the DDI (dipole-dipole interaction), were studied too, with a conclusion that the stationary GSs persist and remain stable in the presence of the arbitrarily strong local self- and cross-component repulsion, compensated by the effectively attractive DDI.

A challenging possibility is to extend the present analysis from solitons to “quantum droplets” in binary dipolar BEC, in the presence of SOC. “Droplets” stabilized by beyond-the-mean-field effects were recently created in single-component dipolar condensates droplets .

Acknowledgements.

This work was supported, in part, by NNSFC (China) through Grant No. 11575063, and by the joint program in physics between NSF and Binational (US-Israel) Science Foundation through project No. 2015616. B.A.M. appreciates a foreign-expert grant of the Guangdong province (China).