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

###### Abstract

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.

Key-words: Rashba spin-orbit coupling, Dresselhaus spin-orbit
coupling, semi-vortices, mixed modes, Lee-Huang-Yang corrections, quantum
droplets.

## 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 Segev1994 (), quadratic Mihalache2006 (), cubic-quintic Mihalache22006 () and nonlocal nonlinearities Mihalache32006 (); Skupin (), which do not give rise to the collapse. Stable quasi-2D matter-wave solitons were predicted in dipolar BECs with long-range interactions Pedri2005 (); Nath2008 (); Tikhonenkov2008 (); Yongyao2013 (); Tikhonenkov22008 (); Raghunandan2015 (); Jiasheng ().

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 SVS1 (); SVS2 (); Bingjin2017 (); SVS3 (); SVSNJP (); SVS4 (); Guihua2017 (); Yongchang (); Yongping (). 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 [or its flipped counterpart, with ] in its two components, SVS1 (), 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 and in the two components SVS1 (), 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 in the two components SVS1 (), or its flipped counterpart, with . In the previously studied systems, the excited states were found to be completely unstable, unlike the fundamental MMs and SVs SVS1 ().

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 we2017 (). 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 SOCoptics ().

Thus, the two-component SOC systems offer a considerable potential for the creation of stable self-trapped modes in 2D and even 3D Yongchang () 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 FR1 (); FR2 (). 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 symbio1 (); symbio2 ().

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 Townes ().

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 LHY () 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 Petrov2015 (); QDterm (); Sadhan-nondip (). Furthermore, the LHY terms were predicted Lima2011 (); Schutzhold2006 (); Saito2016 (); Oldziejewski2016 (); Bisset2016 (); Wachtler2016 (); Pastukhov2017 (); Cinti2017 () and experimentally demonstrated (in the ultracold gas of dysprosium atoms) Chomaz2016 (); Kadau2016 (); Ferrier2016 (); Schmitt2016 (); Baillie2016 (); Edler2017 () to stabilize 2D and 3D self-trapped QDs supported by the long-range dipole-dipole interactions. In particular, Ref. Schutzhold2006 () 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 Leticia ().

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. Zhengwei2013 (), 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. Luca (), 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 Agrawal ()), and SOC of the Rashba type, with strength , is written as

(1) |

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

(2) |

while and are effective 2D coupling constant QDterm (). They are defined so that both are positive or negative, i.e., . 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 is found, it generates a family of solutions rotated by an arbitrary angle, , in the form of

(3) |

By means of further rescaling, one may set , and replace the above-mentioned values of and by either (which we adopt below in the case of ), while is varied, or , while is varied (which is adopted below in the case of ). Note that (no SPM) may also correspond to the model of binary Fermi gases Fermi ().

The energy of the system is

(4) |

In terms of the radial coordinate, , the asymptotic form of the confined modes at is looked for as , where chemical potential is real, while radial wavenumber may be complex, with . The substitution of this in Eq. (1) and linearization easily leads to the expression for in terms of :

(5) |

As it follows from here, the chemical potential satisfying the localization condition, , takes values

(6) |

(recall we have set ). In particular, the limit case of delocalized states, with , corresponds to .

In the case of , 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 symbio1 (); symbio2 ()), the system does not support SVs, but it readily creates MMs, with equal norms in the two components,

(7) |

It is relevant to mention that, in terms of the total norm, , the threshold for the onset of the critical 2D collapse in the framework of Eqs. (1) is

(8) |

where is the norm of single-component Townes soliton Townes ()-Fibich (). 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 , as in Ref. SVS1 (), i.e.,

(9) |

with and real amplitudes . Further, following Ref. SVS1 (), it is possible to construct excited states built on top of the MM, starting from ansatz

(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 Eq. (1) SVS1 ():

(11) |

with chemical potential and real amplitude functions , with nonzero values at , which decay exponentially, as at , see Eq. (5). Below, we use the input in the form of

(12) |

with real amplitudes and , to construct numerically exact SV solutions of Eq. (1).

It is relevant to stress that the norm of the wave function, , 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. SVS3 (), 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 nm for the interatomic attraction, we conclude that in the present notation is tantamount to atoms.

## Iii Numerical results for the mean-field system (without the LHY terms)

### iii.1 Stable fundamental mixed modes

By means of the imaginary-time method ITP1 (); ITP2 (), MM solutions to Eq. (1) can be generated from initial guess (9). A typical example of the so produced MM with 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 SVS1 ().

The numerical solution produces a family of stable MMs in the interval of [recall is imposed in Eq. (1) by rescaling]. In precise agreement with Eq. (8), for fixed stable MMs are found at , and input in the form of Eq. (9) suffers collapse at . The chemical potential, , and energy (4) of the MMs are displayed versus in Figs. 2(a) and (b), respectively. In the limit of , the MM turns into an infinitely broad state with an infinitely small amplitude and the above-mentioned limit value of the chemical potential, . Note that Fig. 2(a) demonstrates that the MM family, which is supported by the dominant attractive XPM, satisfies the Vakhitov-Kolokolov (VK) criterion, , that, in turn, is a well-known necessary condition for the stability of the self-trapped modes VaKo ()-Fibich ().

Mobility of the MMs is a nontrivial issue, as Eq. (1) is not Galilean invariant. Following the pattern of Ref. SVS1 (), the mobility was studied by rewriting Eq. (1) in the reference frame moving with velocity :

(13) |

Here we only consider the case of , because the imaginary-time method could produce stationary solution of Eq. (13) solely in this case, cf. Ref. SVS1 (). Figures 3(a,b) present a typical stable solution of Eq. (13) for and . Note weakly anisotropic deformation of the 2D profile caused by the steady motion in the direction.

The numerical solution demonstrates that, for fixing and , the amplitude of the moving MM, , decays (while it width increases) with the increase of , as shown in Fig. 3(c). In this panel, the amplitude vanishes at , , and , for , , and , respectively, and the MMs do not exists at exceeding these limit values. Thus, the MM’s mobility range expands with the increase of the norm and decrease of , up to . The asymmetry of the moving solitons remains weak even when is close to the limit value. It is relevant to mention that the mobility remains restricted to final values of also in the system with the self-attractive SPM, i.e., SVS1 (). In agreement with a qualitative analysis of Eq. (13), moving solitons with are mirror images, with respect to the -axis, of their counterparts with . Accordingly, the dependence of their amplitude on is the same as shown in Fig. 3(c).

### iii.2 Unstable excited mixed-mode states

By means of the squared-operator method SOM (), 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 , 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.

### iii.3 Stable semi-vortices

In the system with the same signs of the competing SPM (repulsive) and XPM (attractive) terms as above, i.e., , 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., in Eq. (1). Fixing in this case, as said above, by means of rescaling, we find stable SVs with the total norm taking values in the natural interval of . 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 dependence also satisfies the above-mentioned VK criterion, , 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

(14) |

Figure 5(c) displays the widths as functions of for (naturally, the width is larger for the vortical component, , which has the “hole” at its center). The norm share of the vortex component is displayed, as a function of , in Fig. 5(d) for and . For the SVs, unlike the MMs, point is not a critical one, as norms of the two SV’s components are strongly different. However, Figs. 5(c,d) demonstrates that there is a critical value of , with the widths diverging and the norm share of the vortex component, , approaching at . In the same limit, the chemical potential, , approaches the above-mentioned limit value, [see a typical example in Fig. 5(e) for ], i.e., at the SV family displayed in Fig. 5 degenerates into an infinitely extended state with an infinitesimal amplitude. Obviously, is a function of , with in the linear limit, i.e., at . Figure 5(f) displays in the interval of , the maximum value being .

## 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 self-trapped QDs by means of the effective higher-order self-repulsion Petrov2015 ()-Baillie2016 (). In particular, it removes the collapse induced by the attractive XPM at , 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. Zhengwei2013 (). 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. Zhengwei2013 () 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. we2017 (). 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 Petrov2015 (); QDterm (); Sadhan-nondip (), as well as by the very recent experimental work Leticia () reporting the creation of stable QDs supported by purely contact (local) interaction in the gas of K atoms.

According to Ref. QDterm (), 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. QDterm (). However, it is shown in the same work that, in the case of

(15) |

where is the peak value of the two-dimensional atomic density and 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. QDterm (). A straightforward estimate demonstrates that inequality (15) amounts to the following condition for a characteristic size of a localized 2D pattern:

(16) |

where 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. QDterm () (it effectively represents zero-point energy of the Bogoliubov modes on top of the mean-field state):

(17) |

where is the density of the symmetric () flat state corresponding to this expression, and 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

(18) |

where the remaining scaling invariance is used to fix .

Numerical solutions of Eq. (18) of the MM type, which are relevant here because of the opposite signs of the self- and cross-interactions in the terms , 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 .

Numerical solution of Eqs. (18) produces anisotropic (elongated) QD density profiles when . 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 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, , 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 corresponding to different sets in Fig. 7, choosing the vertical orientation for all the QDs. The figure shows that the area and the anisotropy of the 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 as

The degree of the anisotropy and the area of the QD can be defined by as

(19) |

These two characteristics, as well as the chemical potential of the QD, are displayed, as functions of , and in Fig. 8. It clearly demonstrates that: (i) the increase of the nonlinearity strength, , reduces the anisotropy and effective area of the QD; (ii) the increase of the total norm, , reduces the anisotropy and increases the effective, respectively; (iii) opposite to the effect of , 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 is not monotonous. In particular, the non-monotonous dependence implies that the QD family does not satisfy the VK criterion, , although the entire QD family, including the segment with , is found to be completely stable. The violation of the VK criterion as the necessary stability condition VaKo ()-Fibich () 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. SVS3 (), 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

(20) |

where and are strengths of the Rashba and Dresselhaus SOC. Note that the Dresselhaus-only system, with and , is made tantamount to its Rashba-only counterpart, with and , by swapping Sun2017 (). 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. SVS3 () that, for fixed , there are critical values, , 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 , 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 and . To illustrate these findings, we currently fix and vary in the interval of . Typical examples of the the respective numerical solutions for QD with different values of and fixed values of are shown in Fig. 9(a1-a4). In particular, the orientation of the QDs tends to be diagonal. In the limit of , 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, , as

(21) |

where is the counterclockwise rotation of . Figure 9(b) shows as a function of for fixed values of other parameters, . In Fig. 9, the QD profile become practically isotropic at .

## 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 , where 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 Yongchang (). 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.

###### Acknowledgements.

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. B.A.M. appreciates a foreign-expert grant of the Guangdong province (China).## References

- (1) Segev M, Valley G C, Crosignani B, DiPorto P, and Yariv A, 1994 Phys. Rev. Lett. 73, 3211.
- (2) Mihalache D, Mazilu D, Malomed B A, Lederer F, Crasovan L -C, Kartashov Y V, and Torner L 2006 Phys. Rev. E 74, 047601.
- (3) Mihalache D, Mazilu D, Lederer F, Kartashov Y V, Crasovan L -C, Torner L, and Malomed B A 2006 Phys. Rev. Lett. 97, 073904.
- (4) Mihalache D, Mazilu D, Lederer F, Malomed B A, Kartashov Y V, Crasovan L -C, and Torner L 2006 Phys. Rev. E 73, 025601(R).
- (5) Skupin S, Bang O, Edmundson D, and Królikowski W 2006 Phys. Rev. E 73, 066603.
- (6) Pedri P, and Santos L 2005 Phys. Rev. Lett. 95 200404.
- (7) Nath R, Pedri P, and Santos L 2008 Phys. Rev. Lett. 101, 210402.
- (8) Tikhonenkov I, Malomed B A, and Vardi A 2008 Phys. Rev. A 78, 043614.
- (9) Li Y, Liu J, Pang W, and Malomed B A 2013 Phys. Rev. A 88, 053630.
- (10) Tikhonenkov I, Malomed B A, and Vardi A 2008 Phys. Rev. Lett. 100, 090406.
- (11) Raghunandan M, Mishra C, Lakomy K, Pedri P, Santos L, and Nath R 2015 Phys. Rev. A 92, 013637.
- (12) Huang J, Jiang X, Chen H, Fan Z, Pang W, Li Y 2015 Front. Phys. 10, 100507.
- (13) Sakaguchi H, Li B, Malomed B A 2014 Phys. Rev. E 89, 032920.
- (14) Jiang X, Fan Z, Chen Z, Pang W, Li Y, and Malomed B A 2016 Phys. Rev. A, 93, 023633.
- (15) Liao B, Li S, Huang C, Luo Z, Pang W, Tan H, Malomed B A, Li Y 2017 Phys. Rev. A 96, 043613.
- (16) Sakaguchi H, Sherman E Y, and Malomed B A 2016 Phys. Rev. E 90, 032202.
- (17) Sakaguchi H, Malomed B A 2016 New J. Phys. 18 105005.
- (18) Gautam S, and Adhikari S K 2017 Phys. Rev. A 95, 013608.
- (19) Chen G, Liu Y, Wang H 2017 Commun. Nonlinear Sci. Numer. Simulat. 48, 318.
- (20) Zhang Y, Zhou Z, Malomed B A, and Pu H 2015 Phys. Rev. Lett. 115, 253902.
- (21) Zhang Y, Mossman M E, Busch T, Engels P, and Zhang C 2016 Front. Phys. 11, 118103.
- (22) Li Y, Liu Y, Fan Z, Pang W, Fu S, and Malomed B A 2017 Phys. Rev. A 95, 063613.
- (23) Kartashov Y V, Malomed B A, Konotop V V, Lobanov V E, and Torner L 2015 Opt. Lett. 40, 1045.
- (24) Papp S B, Pino J M, and Wieman C E 2008 Phys. Rev. Lett. 101, 040402.
- (25) Zhang P, Naidon P, and Ueda M, 2009 Phys. Rev. Lett. 103, 133202.
- (26) Pérez-García V M, and Belmonte Beitia J 2005 Phys. Rev. A 72, 033620.
- (27) Adhikari S K, 2005 Phys. Lett. A 346, 179.
- (28) Chiao R Y, Garmire E, and Townes C H 1964 Phys. Rev. Lett. 15, 1056.
- (29) Bergé L 1998 Phys. Rep. 303, 259.
- (30) Fibich G The Nonlinear Schrödinger Equation: Singular Solutions and Optical Collapse (Springer: Cham, 2015).
- (31) Lee T. D, Huang K, and Yang C N 1957 Phys. Rev. 106, 1135.
- (32) Petrov D S 2015 Phys. Rev. Lett. 115, 155302.
- (33) Petrov D S, and Astrakharchik G E 2016 Phys. Rev. Lett. 117, 100401.
- (34) Adhikari S K 2017 Phys. Rev. A 95, 023606
- (35) Schützhold R, Uhlmann M, Xu Y, and Fischer U R 2006 Int. J. Mod. Phys. B 20, 3555.
- (36) Lima A R P and Pelster A 2011 Phys. Rev. A 84, 041604(R).
- (37) Saito H 2016 J. Phys. Soc. Jpn. 85, 053001.
- (38) Bisset R N, Wilson R M, Baillie D, and Blakie B P 2016 Phys. Rev. A 94, 033619.
- (39) Wächtler F, and Santos L 2016 Phys. Rev. A 93, 061603(R).
- (40) Ołdziejewski and R, Jachymski K 2016 Phys. Rev. A 94, 063638.
- (41) Pastukhov V 2017 Phys. Rev. A 95, 023614.
- (42) Cinti F, and Boninsegni M 2017 Phys. Rev. A 96, 013627.
- (43) Adhikari S K 2017 Laser Phys. Lett. 14, 025501 (2017).
- (44) Chomaz L, Baier S, Petter D, Mark M J, Wächtler F, Santos L, and Ferlaino F 2016 Phys. Rev. X 6, 041039.
- (45) Kadau H, Schmitt M, Wenzel M, Wink C, Maier T, Ferrier-Barbut I, and Pfau T 2016 Nature (London) 530, 194.
- (46) Ferrier-Barbut I, Kadau H, Schmitt M, Wenzel M, and Pfau T 2016 Phys. Rev. Lett. 116, 215301.
- (47) Schmitt M, Wenzel M, Böttcher F, Ferrier-Barbut I, and Pfau T 2016 Nature (London), 539, 259.
- (48) Baillie D, Wilson R M, Bisset R N, and Blakie P B 2016 Phys. Rev. A 94, 021602(R).
- (49) Edler D, Mishra C, Wächtler Nath F R, Sinha S, and Santos L 2017 Phys. Rev. Lett 119, 050403.
- (50) Cabrera C R , Tanzi L, Sanz J, Naylor B, Thomas P, Cheiney P, and Tarruell L 2017 ePrint, arXiv:1708.07806.
- (51) Zheng W, Yu Z, Cui X, and Zhai H 2013 J. Phys. B: At. Mol. Opt. Phys. 46 134007.
- (52) Cappellaro A, Macrí T, Bertacco G F, and Salasnich L, 2017 Sci. Rep. 7, 13358.
- (53) Agrawal G P Nonlinear Fiber Optics (Academic Press: San Diego, 1995).
- (54) Astrakharchik G E, Boronat J, Casulleras J, and Giorgini S 2004 Phys. Rev. Lett. 93, 200404.
- (55) Chiofalo L M, Succi S, and Tosi P M 2000 Phys. Rev. E 62, 7438.
- (56) Yang J, and Lakoba T I 2008 Stud. Appl. Math. 120, 265.
- (57) Vakhitov M, and Kolokolov A 1973 Radiophys. Quantum Electron. 16, 783.
- (58) Yang J, and Lakoba T I 2007 Stud. Appl. Math. 118, 153.
- (59) Sun F, Ye J, and Liu W M 2017 New J. Phys. 19, 063025.