Multidimensional hybrid Bose-Einstein condensates stabi- lized by lower-dimensional spin-orbit coupling

Y. V. Kartashov,1,2 L. Torner,1,3 M. Modugno,4,6 E. Ya. Sherman,5,6 B. A. Malomed,7 and V. V. Konotop8 1ICFO-Institut de Ciències Fotòniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain 2Institute of Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow, 108840, Russia 3Universitat Politècnica de Catalunya, 08034, Barcelona, Spain 4Department of Theoretical Physics and History of Science, University of the Basque Country UPV/EHU, 48080 Bilbao, Spain 5Department of Physical Chemistry, The University of the Basque Country UPV/EHU, 48080 Bilbao, Spain 6IKERBASQUE Basque Foundation for Science, 48013 Bilbao, Spain 7Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering, and Centre for Light-Matter Interaction, Tel Aviv University, 69978 Tel Aviv, Israel 8Departamento de Física and Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, Edifício C8, Lisboa 1749-016, Portugal


I. Introduction and model
Unlike one-dimensional (1D ) settings, where stable soliton states exist in diverse systems, stability is a major issue for multidimensional self-trapped localized states [1,2]. In 2D and 3D settings with the ubiquitous cubic self-attraction, instability of fundamental solitons is driven, respectively, by critical and supercritical collapse [3][4][5], while vortex solutions are subject to azimuthal self-splitting instability [2]. Several stabilization mechanisms for multidimensional solitons have been elaborated theoretically. In a recent landmark advance, stabilization of Bose-Einstein condensates (BECs) against collapse [6,7] by Lee-Huang-Yang (LHY) quantum corrections to the mean-field interactions [8] has been demonstrated experimentally, leading to the creation of fundamental (zero-vorticity) solitary states in the form of "quantum droplets" in BECs with dipole-dipole [9][10][11] and with contact [12][13][14][15] interactions. Regarding "droplets" with embedded vorticity, they are unstable in the former case [16], while LHY-stabilized 3D [17] and 2D [18,19] vortical droplets have been predicted in BECs with contact nonlinearity.
Another possibility for the creation of stable solitons in atomic BECs was predicted in the model of spinor condensates with cubic attraction and spin-orbit coupling (SOC) between the two components [20][21][22] (see [23] for a similar result for optical solitons and [24] for the stabilization of solitons in BECs by Zeeman lattices). Previously, it was assumed that the cubic self-attraction in the 2D geometry always leads to critical collapse when the norm of the wave function, U , exceeds a critical value, cr U , while any input decays at cr U U  , hence the soliton solution existing at cr U U  is unstable [3,4]. In Refs. [20][21][22] it was shown that SOC may change the situation at cr U U  , creating stable 2D solitons, which represent the otherwise missing ground state in the system. In 3D, the supercritical collapse (which in the presence of SOC was studied in [25]) does not let SOC create the ground state, although 3D soliton solutions stable against small perturbations have been predicted [26].
In this paper, we explore the possibility to create stable 2D and 3D self-trapped states supported by lower-dimensional SOC (1D and 2D, respectively). Such settings suggest the existence of new species of multidimensional solitons, which strongly differ from the solitons supported by the full SOC. Their stability is a particularly challenging issue, as reduction of the dimensionality of the support structure makes it harder to secure the stability, cf. the study of fundamental and vortex 2D and 3D solitons supported, respectively, by the 1D and 2D lattice potentials [2,[27][28][29]. The use of lower-dimensional SOC settings offers a crucial advantage for the experimental creation of solitons, as the majority of experimental realizations of SOC were performed in the 1D geometry [30][31][32], with 2D schemes implemented in a few works [33,34], but not yet in 3D.
The evolution of a spinor BEC in space of dimension D is modelled by coupled Gross-Pitaevskii equations (GPEs) for a spinor wavefunction [30][31][32]: Here , , admitting a difference of the crossc ( ) g and selfs ( ) g interactions ( s c g g  corresponds to Manakov's symmetric system [37]). Localized solutions of Eq. (1) are characterized by the norm, U    .
The analysis for 2D and 3D solitons, which are supported, respectively, by the reduced one-and two-dimensional SOC, is presented below in Sections II and III. Results of a systematic numerical investi-gation are combined with analytical findings, which are based on the consideration of the spectrum of the linear version of Eq. (1), as well as on prediction of the stability of the multidimensional solitons produced by the well-known Vakhitov-Kolokolov (VK) criterion. In addition, essential results for 2D solitons are obtained by means of the variational approximation (VA). Results are summarized in Section IV, which also outlines possibility of their experimental implementation.

II. 2D system with 1D spin-orbit coupling
First, we consider the 2D spinor BEC, with  36]. Then, one may expect that such a Zeeman lattice stabilizes 2D solitons [28,29]. However, it is more relevant to start the analysis with identifying the linear spectrum of Eq. (1) for the chemical potential  that can be obtained by substitution , where C is constant spinor and x,y k are wavenumbers of the respective excitation, into the linearized version of Eq. (1). A straightforward calculation yields the spectrum which consists of two branches, The critical nature of the 2D collapse in the system with 0   (without SOC) is underlain by the fact that both the kinetic and nonlinear-interaction energies scale as 2 x,y k , while the total norm takes a single value cr U U  for all solitons. The breaking of the scaling (conformal) invariance by SOC (with 0   ) and lifting the norm's degeneracy may create a stability range for 2D solitons at Under the action of the self-attraction, solitons that have the form of , with complex spinor  T 1 2 ( ( , ), ( , )) w x y w x y W , bifurcate upon variation of chemical potential  from the minimum of the lower branch of the linear dispersion relation given by Eq. (2),  ,2)]. These 2D solitons, featuring the dipole structure in the first component, are drastically different from the self-trapped 2D states ("semi-vortices" and "mixed modes" [20,21]) supported by the full 2D SOC. In addition to the different symmetry, the latter states never show striped patterns. In a sense, the 2D solitons displayed in Fig. 2 are hybrids, which combine the 2D stability with structural features resembling those found in 1D solitons, cf. Refs. [23,24].   Fig. 1(b) is non-monotonous. Also, in contrast to the solitons supported by the full SOC [21,22,25], but similar to 1D states maintained by SOC with the Zeeman lattice [24], the present 2D solitons exist above a threshold value of the norm, th min ( ) U U  , which vanishes solely at 2 2   , being a nonmonotonous function of  and  . This is illustrated by Fig. 1(c), which compares ( ) U  dependencies for the solitons with the singlepeak 2 ( 2 )    and striped structures. The transition between these species occurs via shapes strongly elongated along y [ Fig. 2(c)], clearly indicating hybrid nature of such states. An example of a single-peaked soliton is presented in Fig. 2(d). In all cases, irrespective of the ratio of  and 2  , at    the norm approaches that of the Townes soliton [38], Townes 5.85 U  , which is the single possible value for (unstable) 2D solitons in the absence of SOC [3,4].
To explain the vanishing of th U at the point of transition between the 2D solitons with single-peak and striped structures, , we note that, at this point, the expansion of dispersion relation (2) near the origin ( ) For small deviation of the chemical potential from the minimum, min we explore the bifurcation of solitons from the lower branch of the linear spectrum, i.e., from the state T (0,1) . Accordingly, we look for a stationary solution in the form of x y  is the slowly varying amplitude. Applying the multiple-scale expansion (see Appendix A for the details of derivation), we find that the amplitude solves a stationary equation, which is written in variables Ref. [39]). This equation gives rise to an exact scaling relation, , which satisfies the necessary stability condition given by the VK criterion, / 0 U     [3,4,40]. Soliton solutions of Eq. (3) can be predicted by means of VA [41] based on the Gaussian ansatz with norm U : The VA yields strongly anisotropic relations for parameters of the ansatz: . In all cases, the stability exactly follows the VK criterion, the branches with / 0 U     and / 0 U     [the black and red ones in Figs. 1(b,c)] being stable and unstable, respectively. Thus, 1D SOC stabilizes almost the entire soliton family, except for small segments with . It is very plausible that the stable soliton realizes, for given norm, the ground state of the 2D system, although rigorous proof of this conjecture requires additional analysis. Stable propagation of a perturbed 2D soliton is displayed in Fig. 2(e). Unstable broad solitons transform into much narrower stable ones with the same norm, as shown in Fig. 2(f).
Our results remain valid for non-Manakov nonlinearity, with

III. 3D system with 2D Rashba spin-orbit coupling
The expected single-peak and striped 3D solitons again bifurcate from minima of the   branch in Eq. (5). Solitons have the form of   x y plane. The difference between vorticities of the two components of  , m and 1 m  , is imposed by the 2D SOC. Here we address the modes with 1 m   , i.e., zero vorticity in the dominant second component, which minimizes the total energy, thus having the best chance to be stable.
The families of 3D solitons at s,c 1 g  are displayed in Fig. 4(a). Their amplitude vanishes, and size , like in the 2D case. However, their norm vanishes at    , i.e., 3D solitons have no existence threshold. At 2 2   the striped structure manifests itself in concentric amplitude-phase modulation in the ( , ) x y plane [ Fig. 5(a)], while the soliton is elongated along the z -axis. The interplay of the vorticity in the 1  component and striped radial modulation builds a complex phase distribution, which makes the stable 3D solitons radically different from those supported by the full 3D SOC, cf. Ref. [26]. Similar to what is said above about the 2D solitons, the present modes may be considered as hybrids combining 3D and 1D features. With the increase of  , they become more localized and the radial modulation gradually disappears [ Fig. 5(b)].
As seen in Fig. 4(a), the 3D system with 2D SOC produces striped solitons with a non-monotonous dependence ( ) U  , which includes a VK-stable segment with  Fig. 4(b)], at which the striped structure is replaced by the single-peak one, and ( ) U  dependence abruptly changes, again becoming monotonous [ Fig. 4(a)]. Thus, in contrast to the 2D case, the norm of the 3D solitons does not vanish at 2 2   , where striped shape is replaced by a single-peak one, while, as in the 2D case, the stability exactly follows the VK criterion, / 0 U     , even if only the striped solitons are stable in 3D. Similar results are obtained for s c g g  . Unstable 3D solitons are destroyed by the fast collapse, as shown in the top row of Fig. 6. These high-amplitude solitons belong to segments of the ( ) U  curves with / 0 U     , to the left of the stability domain. The evolution of a stable 3D soliton, belonging to the / 0 U     branch, is displayed in the bottom row of Fig. 6. In that case, initial perturbations excite only weak oscillations of the soliton's amplitude. Strong compression of 3D solitons, which are stable against small perturbations, may initiate their supercritical collapse, therefore, they are metastable states, separated by an - and - dependent barrier from the collapse [42]. , showing the shape of 3D solitons, and the respective amplitude and phase profiles in the 0 z  plane (the left, central, and right columns, respectively) for (a) 6.85 In a recent experiment [43] with one-dimensional SOC, the ratio 2 /   has been modified from 0 to 2.6 , which, according to the above results, allows the existence of stable 2D single-peak and stripe soliton states with the - dependent modulation length 3 m   . The strength of the nonlinear interaction in 2D, written in physical units, is , where M is the atomic mass, s a the scattering length (the self-attraction corresponds to 0 s a  ), and z a the zaxis confinement length. For typical values of s a and z a this yields the Townes norm corresponding to 3 10  atoms. Note that for twodimensional SOC [33], the ratio 2 /   is limited by 0.2, and for maintaining a necessary suppression effect of SOC on the 3D collapse one has to limit the number of atoms to

IV. Conclusions
In summary, we have shown that even a lower-dimensional spinorbit-coupling than that of the embedding space can stabilize soliton states in spinor BECs with intrinsic attraction, when it is combined with Zeeman splitting. On physical grounds, the stabilization arises from the nonparabolic dispersion of the system and, importantly, holds in a broad parameter region, which complies with the Vakhitov-Kolokolov criterion. We also found that the multidimensional states feature hybridization of 1D and 2D/3D structural properties. The results are important for the creation of stable 2D and 3D solitons in BECs, as spin-orbit coupling is currently experimentally available solely in low-dimensional forms.  . Looking for a stationary solution of (A1) in the form , where  depends on the coordinates only, we obtain