Vortex solitons: Old results and new perspectives

TL;DR

This review discusses the stability, creation, and recent advances in vortex solitons across 2D and 3D systems, highlighting stabilization methods, new schemes involving spin-orbit coupling, and experimental efforts.

Contribution

It provides a comprehensive overview of stabilization techniques for vortex solitons, including novel schemes like spin-orbit coupling and microwave coupling, and introduces the concept of robust quantum droplets.

Findings

Stabilization of vortex solitons via competing nonlinearities and trapping potentials.

Prediction of stable quantum droplets with vortex structures.

Development of new schemes for creating stable 2D and 3D vortex solitons.

Abstract

A review is given of some well-known and some recent results for two- and three-dimensional (2D and 3D) solitons, with emphasis on states carrying embedded vorticity. Unlike typically stable 1D solitons, 2D and 3D ones are vulnerable to instabilities induced by the critical and supercritical collapse, respectively, in the 2D and 3D models with cubic focusing nonlinearity. Vortex solitons are subject to a still stronger splitting instability. For this reason, a central problem is looking for physical settings in which 2D and 3D solitons may be stabilized. The review addresses, first, two well-established topics, viz., the stabilization of vortex solitons by means of competing nonlinearities, or by trapping potentials (harmonic-oscillator and spatially-periodic ones). The former topic includes a new addition, viz., the prediction of robust quantum droplets. Two other topics outline new schemes elaborated for the creation of stable vortical solitons. One scheme relies on the use of spin-orbit coupling (SOC) in binary BEC, making it possible to predict stable 2D and 3D solitons, which couple or mix components with vorticities S=0 and (+/-)1 (semi-vortices (SVs) or mixed modes (MMs), respectively). In this system, the situation is different in the 2D and 3D geometries. In 2D, SOC helps to create a ground state (GS), represented by stable SV or MM solitons, whose norm falls below the threshold value at which the critical collapse sets in. In the 3D geometry, the supercritical collapse does not allow one to create a GS, but metastable solitons of the SV and MM types can be constructed. Another new scheme makes it possible to create stable 2D vortex-ring solitons with arbitrarily high S in a binary BEC with components coupled by microwave radiation. Some other topics are addressed too, such as vortex solitons in dissipative media and attempts to create vortex solitons in experiments.

Tables

Table 1. Table 1: Table I: The middle line shows numerically found minimum values of the norm, N thr subscript 𝑁 thr N_{\mathrm{thr}} , necessary for the stability of 2D vortex QDs (quantum droplets) with winding number S 𝑆 S , as per Ref. GZ-ln . The bottom line presents results produced by the analytical approximation ( 24 ), which is relevant for S ≥ 3 𝑆 3 S\geq 3 .

$S$	$1$	$2$	$3$	$4$	$5$
$N_{\mathrm{thr}}$	$60$	$200$	$510$	$1380$	$3550$
$N_{\mathrm{thr}}^{\left(\mathrm{analyt}\right)}$	n/a	n/a	$486$	$1536$	$3750$

Table 2. Table 2: Table II. β collapse subscript 𝛽 collapse \beta_{\mathrm{collapse}} and β collapse ( analyt ) superscript subscript 𝛽 collapse analyt \beta_{\mathrm{collapse}}^{\mathrm{(analyt)}} are numerically obtained and analytically predicted (see Eq. ( 94 )) values of the contact-interaction strength, β 𝛽 \beta , up to which the fundamental and vortex solitons exist, in the framework of Eq. ( 87 ). β stability subscript 𝛽 stability \beta_{\mathrm{stability}} is the numerically identified stability boundary of the vortex solitons.

$S$	${\beta }_{\mathrm{collapse}}$	${\beta }_{\mathrm{collapse}}^{\left(\mathrm{analyt}\right)}$	${\beta }_{\mathrm{stability}}$	$S$	${\beta }_{\mathrm{collapse}}$	${\beta }_{\mathrm{collapse}}^{\left(\mathrm{analyt}\right)}$	${\beta }_{\mathrm{stability}}$
$0$	$11.8$	$n/a$	$\equiv {\beta }_{\mathrm{collapse}}$	$3$	$132.5$	$130.6$	$41$
$1$	$48.3$	$43.5$	$11$	$4$	$175.5$	$174.1$	$57$
$2$	$89.7$	$87.0$	$28$	$5$	$218.5$	$217.7$	$70$

Equations

$$i\frac{\partial }{\partial z}{u}_j^{(\mathrm{FF})}+{\left(-1\right)}^j\frac{b}{2}{u}_j^{(\mathrm{FF})}$$

$$+\frac{1}{2}\left(\frac{{\partial }^2}{\partial x^2}+\frac{\partial }{\partial y^2}+D_{\mathrm{FF}}\frac{{\partial }^2}{\partial {\tau }^2}\right)u_j^{(\mathrm{FF})}+u^{(\mathrm{SH})}{\left(u_{3-j}^{(\mathrm{FH})}\right)}^{\ast }=0,$$

$$2i\left(\frac{\partial }{\partial z}+ic\frac{\partial }{\partial \tau }\right)u^{(\mathrm{SH})}-q{u}^{(\mathrm{SH})}$$

$$+\frac{1}{2}\left(\frac{{\partial }^2}{\partial x^2}+\frac{\partial }{\partial y^2}+D_{\mathrm{SH}}\frac{{\partial }^2}{\partial {\tau }^2}\right)u^{(\mathrm{SH})}+u_1^{(\mathrm{FF})}{u}_2^{(\mathrm{FF})}=0,$$

$$\tau \equiv t-z/V_{\mathrm{gr}}$$

$$b=0,u_1^{(\mathrm{FF})}=u_2^{(\mathrm{FF})}\equiv u^{(\mathrm{FF})}/\sqrt{2}.$$

$$i\frac{\partial u}{\partial z}+\frac{1}{2}\left(\frac{{\partial }^2}{\partial x^2}+\frac{\partial }{\partial y^2}\right)u+\frac{|u{|}^{2}u}{1+|u{|}^2/U_0^2}=0,$$

$$u\left(x,y,z\right)=\exp \left(ikz+iS\theta \right)U(r),$$

$$\frac{d^{2}U}{d{r}^2}+\frac{1}{r}\frac{dU}{dr}-\frac{S^2}{r^2}U+\frac{2U^3}{1+U^2}=2kU,$$

$$U\sim r^{|S|},$$

$$U\sim r^{-1/2}\exp \left(-\sqrt{2k}r\right)$$

$$i\frac{\partial u}{\partial z}+\frac{1}{2}\frac{{\partial }^{2}u}{\partial {\tau }^2}+\frac{1}{2}\left(\frac{{\partial }^2}{\partial x^2}+\frac{\partial }{\partial y^2}\right)u+|u{|}^{2}u-|u{|}^{4}u=0.$$

$$u_{\mathrm{sol}}(z,\tau )=2e^{ikz}\sqrt{\frac{k}{1+\sqrt{1-16k/3}\cosh \left(2\sqrt{2k}\tau \right)}},$$

$$0<k<3/16.$$

$$u\left(x,y,z,\tau \right)=\exp \left(ikz+iS\theta \right)U(r,\tau ),$$

$$\frac{{\partial }^{2}U}{\partial r^2}+\frac{1}{r}\frac{\partial U}{\partial r}-\frac{S^2}{r^2}U+\frac{{\partial }^{2}U}{\partial {\tau }^2}+2U^3-2U^5=2kU.$$

$$R\approx \left(1/\sqrt{6}\right)\ln \left({\left(3/16-k\right)}^{-1}\right).$$

$$U^4-U^2+\left(k+\frac{S^2}{2r^2}\right)=0.$$

$$r<{\left(r_0\right)}_{\mathrm{TF}}=\frac{S}{\sqrt{1/2-2k}},$$

$$k_{\min }(S)<k<3/16=0.1875$$

$$k_{\min }(1)=0.1487,k_{\min }(2)=0.1619,k_{\min }(3)=0.1700,k_{\min }(4)=0.1769,k_{\min }(5)=\mathrm{0.1806.}$$

$$U\left(R,\delta \right)=C(-1{)}^{S+1}{R}^{-1//2}\exp \left(-\sqrt{-2\mu }R\right)\cos \delta ,$$

$$i\frac{\partial \psi }{\partial t}=-\frac{1}{2}{\nabla }_{3D}^2\psi -|\psi {|}^2\psi +|\psi {|}^3\psi ,$$

$$i\frac{\partial \psi }{\partial t}=-\frac{1}{2}{\nabla }_{2D}^2\psi +\ln \left(|\psi {|}^2\right)|\psi {|}^2\psi .$$

$$H=\frac{1}{2}\int \int \left[{\left|{\nabla }_{2D}\psi \right|}^2+|\psi {|}^4\ln \left(\frac{|\psi {|}^2}{\sqrt{e}}\right)\right]dxdy,$$

$$N_{\mathrm{thr}}^{(\mathrm{analyt})}\approx 6N^4$$

$${\psi }_{\mathrm{pert}}\left(x,y\right)\approx {\left(x+iy\right)}^S-\epsilon \left(x+iy\right)$$

$$\equiv \left(x+iy\right)\left[{\left(x+iy\right)}^{S-1}-\epsilon \right].$$

$$x_{\mathrm{piv}}^{(1+j)}+i{y}_{\mathrm{piv}}^{(1+j)}={\epsilon }^{1/(S-1)},j=1,...,S-1.$$

$${\psi }_{\mathrm{pert}}\left(x,y\right)={\left(x+iy\right)}^S-\epsilon \left(x-iy\right).$$

$$x_{\mathrm{piv}}^{(k)}+i{y}_{\mathrm{piv}}^{(k)}={\epsilon }^{1/(S+1)}{\left|\epsilon \right|}^{2/\left(S^2-1\right)},k=1,...,S+1,$$

$$i\frac{\partial \psi }{\partial t}=-\frac{1}{2}{\nabla }^2\psi +U(x,y,z)\psi -|\psi {|}^2\psi .$$

Full text

List of acronyms: 1D – one-dimensional; 2D – two dimensional; 3D – three-dimensional; BEC – Bose-Einstein condensate; CGLE – complex Ginzburg-Landau equation; CQ – cubic-quintic (nonlinearity), FF – fundamental frequency; GPE – Gross-Pitaevskii equation; GS – ground state; GVD – group-velocity dispersion; HO – harmonic oscillator (potential); HV – hidden vorticity; LHY – Lee-Huang-Yang (correction to the mean-field dynamics of BEC); MF – mean field; MM – mixed mode; NLSE – nonlinear Schrödinger equation; OL – optical lattice; PT – parity-time (symmetry); QD – quantum droplet; SH – second harmonic; SOC – spin-orbit coupling; SV – semi-vortex; TF – Thomas-Fermi (approximation); TS – Townes’ soliton; VA – variational approximation; VAV – vortex-antivortex (composite state in a two-component system); VK – Vakhitov-Kolokolov (stability criterion)

I Introduction: The topic and objectives of the article

I.1 The concept of vortex solitons

Since term “soliton” was coined by Zabuski and Kruskal ZK , it is commonly applied to localized self-trapped modes which spontaneously emerge in diverse media as a result of the balance between the diffraction and/or material dispersion affecting the propagation of linear waves, and self-focusing (self-attractive) material nonlinearity. Theoretical and experimental studies of solitons have grown into a huge research area penetrating many disciplines, in physics, mathematics, and beyond. The largest share of these works have been dealing with one-dimensional (1D) solitons KA ; Peyrard , the dominant areas in current studies of solitons being nonlinear optics KA and Bose-Einstein condensates (BECs) in ultracold atomic gases Brazh ; BEC-sol1 ; Morsch ; PhysicaD ; BEC-sol2 ; BEC-sol3 . The extension of the concept of solitons to two- and three-dimensional (2D and 3D) geometry is a highly nontrivial generalization (multidimensional optical solitons, which are self-trapped (localized) both in the spatial and temporal directions, are often called spatiotemporal solitons old-review , alias “light bullets” Silb ). A well-known problem impeding the studies of multidimensional solitons is that, while a majority of 1D solitons, such as ones described by the Korteweg - de Vries Kruskal , nonlinear Schrödinger equation ZS , sine-Gordon AKNS , and Landau-Lifshitz Borovik ; Mikhailov ; Rodin equations, may be produced as solutions of integrable Zakharov ; Newell ; Ablowitz or nearly integrable 1989 models, only few 2D equations are integrable (most essential among them are Kadomtsev-Petviashvili equations of two types Dryuma , with opposite signs of the dispersion, only one of them, the so-called KP-I equation, giving rise to 2D solitons), and there are no integrable 3D equations which would find any realization in physics. The lack of (nearly) integrable models obviously makes the theoretical study of multidimensional solitons more difficult. Another fundamental issue is that the most common cubic self-focusing nonlinearity tends to create 2D and 3D solitons which are subject to instability driven by the wave collapse (alias blowup) that occurs in the same equations, i.e., spontaneous formation of singularities, driven by the self-focusing effect Askaryan ; Townes , after a finite evolution time Berge ; SulemSulem ; Gadi . The collapse induced by the cubic nonlinearity is critical in 2D, which means that it sets in when the norm of the underlying wave field exceeds a certain finite critical value, and supercritical in 3D, where an initial state with an arbitrarily small norm may blow up due to the collapse. A clear experimental demonstration of the collapse of an optical vortex beam (in water), which is accompanied by spontaneous breaking of the azimuthal symmetry of the collapsing beam, and is obviously relevant to the topic under the consideration, was reported in Ref. Gaeta ).

On the other hand, solitons created in 2D and 3D geometries offer a great variety of new possibilities. Most significant, vorticity may be embedded in both 2D and 3D localized states, thus creating vortex solitons, which is the main subject of the present review. The intrinsic vorticity is characterized by an integer winding number, alias topological charge, $S$ (“spin”), which is defined through a total change of the phase, $\Delta\varphi=2\pi S$, accumulated in the course of a trip along a closed trajectory surrounding the vortex’ pivot (phase singularity), placed at $r=0$, where $r$ is the radial coordinate. The presence of the vorticity implies that, similar to the commonly known asymptotic structure of the Bessel functions, $J_{|S|}(r)$, the local amplitude of the wave field(s) in vortex solitons must vanish $\sim r^{|S|}$ at $r\rightarrow 0$ (here $|S|$ is written as integer $S$ may be negative). Due to this peculiarity, the vortex soliton features a “hole” in the center, i.e., the localized mode is shaped as a ring (or, quite frequently, as a broad annulus) in 2D, and a torus (“donut”) in 3D. Nevertheless, an exception is known: in the 2D model based on the Gross-Pitaevskii equation (GPE) for a BEC with the quintic self-defocusing interaction and an external potential $U=-U_{0}r^{-2}$, $U_{0}>0$, which pulls atoms to the center, stable vortex states are possible in which the amplitude does not vanish at $r\rightarrow 0$, but instead diverges $\sim r^{-1/2}$ pull . This states is a result of the balance between the pull to the center and quintic self-repulsion. Indeed, the absolute values of the wave function in a vortex state with winding number $S\neq 0$ cannot be finite at $r=0$, as its phase, $\varphi=S\theta$ (where $\theta$ is the angular coordinate) is not defined at $r=0$. The usual solution of this problem is resorting to solutions with the amplitude vanishing at $r=0$. However, an alternative solution is possible too, with the amplitude diverging at $r\rightarrow 0$, as in Ref. pull . Indeed, it is commonly known that the Bessel equation, which produces the asymptotic form $\sim r^{|S|}$ of the solution for the vortex’ amplitude, also has the singular solution (given by the Neumann function), $\sim r^{-|S|}$. In the 2D setting, the latter solutions with any $|S|\geq 1$ is unphysical, as it gives rise to a divergent norm. Nonetheless, the above-mentioned specific singularity, $\sim r^{-1/2}$, which is produced by the interplay of the vortex structure with the pulling potential and quintic self-defocusing, is integrable (the corresponding integral norm converges at $r\rightarrow 0$), hence this solution is a physically relevant one.

In 3D, solitons with more sophisticated topological structures are known too, in the form of skyrmions skyrmion1 ; skyrmion2 ; skyrmion3 ; semi1 ; semi2 ; VPG-skyrmion , hopfions, hopfion ; Yasha , knots super ; Sutcliffe ; Radu , etc., which may carry two independent topological charges. In particular, hopfions may be realized in terms of a single complex wave field, which forms a twisted vortical torus, with the phase of the wave function winding both along the torus-forming ring, which determines the overall vorticity of the hopfion, and along the ring in the torus’ cross section, which determines the intrinsic twist of the hopfion Yasha . In addition to physics of nuclear matter and the classical field theory, which are the origin of the Skyrme model and diverse 3D states generated by it skyrmion1 ; skyrmion2 ; skyrmion3 ; Sutcliffe ; Radu ; field-theory , these complex 3D modes find important realizations in ferromagnets ferro ; ferro-Sutcliffe , semiconductors semi1 ; semi2 , superconductors super , and in various configurations of BEC.

In comparison with the fundamental (zero-vorticity) solitons in the same geometries, the creation of 2D and 3D vortex rings and tori is a still more challenging objective, because, in addition to the above-mentioned collapse-driven instability, they are subject to an even stronger azimuthal instability, which tends to break the axially symmetric ring or torus into fragments, each one being, roughly speaking, a fundamental soliton, see illustration below in Fig. 1 Minsk1 ; Minsk2 . Therefore, theoretical prediction and experimental realization of physically relevant settings that may support stable fundamental and, especially, vortex solitons in 2D and 3D geometries is a challenging objective, which has drawn much interest in the course of the two last decades. It was the subject of several reviews, which were published both relatively long ago and more recently, being chiefly focused on selected aspects of the broad topic old-review ; Dum ; me ; NatureRev .

The present article aims to provide a review specifically focused on vortex solitons. Because this theme is a vast one too, the review addresses some selected aspects in detail, and others in a brief form, as a comprehensive review might easily grow to the size of a book. Two first topics selected for the presentation are well-established ones. They rely on the stabilization of vortex states by competing nonlinearities Manolo ; Pego ; nine ; CQ2comp2D ; FWM ; S=1and2chi2chi32D ; S=3and4chi2chi32D ; chi2chi33D , or in trapping potentials 2D ; Sadhan ; Ueda ; Ueda2 ; Dum2D ; Dum3D (3D and 2D harmonic-oscillator (HO) potentials are considered in detail, and spatially periodic lattices, which were a subject of several earlier published reviews, are addressed in a brief form). The former topic includes a new addition, viz., the consideration of the recently created Tarr1 -Ing2 quantum droplets (QDs). Models of QDs make it possible to predict robust vortex solitons with the help of specific competing nonlinearities 3DLHY ; GZ-ln ; necklace ; HS , as well as of necklace-shaped circular chains of QDs, which may carry the angular momentum necklace (unlike fundamental QDs, creation of such modes has not yet been reported).

Two other topics in the review present entirely recent developments. These are schemes for the creation of stable 2D and 3D solitons in models of binary BEC, with its two components interacting by means of the spin-orbit coupling (SOC) we ; Sherman2 ; HP , or through a resonant microwave field Jieli2 . The latter section also briefly outlines results for vortex solitons obtained in other models with effectively nonlocal nonlinearities. The section of the article preceding the concluding one touches upon some other topics which are relevant to the broad theme of vortex solitons, such as settings with spatially modulated local nonlinearities, discrete media, and, in a somewhat more detailed form, nonlinear dissipative and $\mathcal{PT}$-symmetric systems. The same section also briefly mentions the state of the art as concerns experimental creation of vortex solitons (thus far, they were only observed as transient states, in very specific settings).

I.2 A related field: “dark” vortices supported by a

finite background

A broad topic related to bright vortex solitons, i.e., localized states with embedded vorticity, is the prediction and creation of solitons in media with self-repulsive nonlinearities. These are 2D (and 3D) vortex states, embedded in flat modulationally stable background fields. These states may be considered as 2D counterparts of dark solitons, which are well-known solutions of the integrable 1D nonlinear Schrödinger equation ZS . In fact, these delocalized vortex modes are most frequently called “vortices”, in the context of optics and BEC alike. They are well-known objects, studied in detail theoretically and created in experiments, prior to the work on bright vortex solitons. Results obtained for the “dark vortices”, starting from the earliest ones Coullet ; Neu ; Swartz1 , have been collected in many publications dealing with various settings in photonics Swartz2 ; Soskin ; vort-review-early ; Vyslo ; Bessel-beam-DiTrapani ; vort-in-liq-cryst ; vort-in-communications (including exciton-polariton condensates polaritons1 ; polaritons ; polaritons2 , where the vortices exist in dissipative media) and BEC Aftalion ; Cornell-vortex ; Cornell2 ; Cornell3 ; Hall ; TsuKasa ; BEC-exper ; Fetter , as well as in quantum Fermi gases Fermi ; Fermi-review ; see also book black for an overview of the theme. The studies of vortices in this context are also known as “singular optics”, with pivots of the vortices considered as singularities of optical fields; many publications on the latter theme were collected in a special issue of Journal of Optics Editorial (see also book Gbur ).

It is relevant to mention too that vortices are closely related to the concept of the optical angular momentum Allen ; opt-angular-mom ; Bliokh-Nori , and the possibility of storage of the angular momentum in ultracold gases angular-momentum-storage .

As concerns the stability of the “dark vortices”, in the 2D geometry the stability of the vortex with $S=1$ is secured by the conservation of the associated angular momentum, while multiple vortices, with $S\geq 2$, are usually unstable against splitting into $S$ unitary vortices Neu . The instability is qualitatively explained by the fact that two separated vortices with equal topological charges repel each other. In the 3D geometry, even the vortex with $S=1$ may be unstable against spontaneous bending of its pivotal line, as predicted theoretically Svidzinsky ; bending ; bending2 and observed experimentally Dalibard-bending ; Bagnato-bending .

Lastly, it is relevant to mention a 2D two-component system which gives rise to stable bound states with a zero-vorticity component filling the central “hole” of the dark-vortex structure induced in the other component Berloff .

II Earlier results and recent additions to them: competing

nonlinearities and trapping potentials

II.1 Unstable vortex solitons in media with quadratic and saturable

nonlinearities

As said above, the first objective in the studies of multidimensional solitons is securing their stability. One possibility is to consider media with collapse-free nonlinearities. In optics, they may be represented by quadratic (alias $\chi^{(2)}$, or second-harmonic-generating) terms chi1 ; chi2 ; chi3 ; chi4 , which do not lead to collapse in 2D and 3D geometries Rubi ; HaoHe . This circumstance stimulated the experimental creation of fundamental spatiotemporal solitons in a quasi-2D $\chi^{(2)}$ setting (the localization in the third direction was provided by a waveguiding structure, which was necessary, as this direction was “sacrificed” for inducing sufficiently strong group-velocity dispersion (GVD), that was necessary for self-trapping in the temporal direction) Frank1 ; Frank2 . There are two basic forms of the quadratic nonlinearity, which are usually referred to as “Types I and II”. They correspond, respectively, to the degenerate two-wave system and full three-wave one. The latter system is modeled by three $\chi^{(2)}$-coupled equations for the paraxial propagation of two components $u_{1,2}^{\mathrm{(FF)}}$ of the fundamental-frequency (FF) wave, which represent two mutually orthogonal polarizations of light, and amplitude $u^{\mathrm{(SH)}}$ of the second-harmonic (SH) wave chi1 ; chi2 ; chi3 ; chi4 :

[TABLE]

with $j=1,2$, where $\ast$ stands for the complex conjugate, $z$ is the propagation distance, the paraxial-diffraction operator $(1/2)\left(\partial^{2}/\partial x^{2}+\partial^{2}/\partial y^{2}\right)$ acts on functions of transverse coordinates $\left(x,y\right)$,

[TABLE]

is the reduced time, which combines the usual temporal variable, $t$, and the propagation distance, $V_{\mathrm{gr}}$ being the group velocity of the FF’s carrier wave KA , and walkoff coefficient $c$ is proportional to the mismatch between group velocities of the SH and FF waves. Further, real $b$ accounts for the phase birefringence of the SH components, $D_{\mathrm{FF,SH}}$ are group-velocity-dispersion (GVD) coefficients of the FF and SH waves, and real $q$ is the phase mismatch between them, while the diffraction and $\chi^{(2)}$ coefficients are scaled to be $1$. The system creates bright spatiotemporal solitons in the case of the anomalous GVD, i.e., $D_{\mathrm{FF,SH}}>0$. The degenerate Type-I system is obtained from Eqs. (1) and (2) by setting

[TABLE]

Further analysis has demonstrated that the $\chi^{(2)}$ nonlinearity of either type (I or II), while readily maintaining stable fundamental 2D and 3D solitons, cannot stabilize solitons with embedded vorticities $S_{1,2}^{\mathrm{(FF)}}$ and $S^{\mathrm{(SH)}}=S_{1}^{\mathrm{(FF)}}+S_{2}^{\mathrm{(FF)}}$ of the FF and SH components against splitting, even in the 2D settings, for any nonzero values of $S_{1,2}^{\mathrm{(FF)}}$ (in the two-wave system, one has the single value of $S^{\mathrm{(FF)}}$ and, accordingly, $S^{\mathrm{(SH)}}=2S^{\mathrm{(FF)}}$). This negative result was predicted theoretically for the two-wave Dima1 ; Petrov1 ; Petrov2 ; Dima2 and three-wave 3wave ; Herve-3waves systems in 2D (which implies the consideration of $\tau$-independent solutions of Eqs. (1) and (2)), and demonstrated experimentally for the former one Petrov . The instability of bright vortex states makes it possible to realize experimentally their controllable splitting in a set of fundamental solitons, by means of interaction with an additional seed beam, which was taken too as a vortical one stimulated . Nevertheless, in some cases the instability may be very weak Herve-3waves , which may allow one to produce quasi-stable vortex solitons.

Another option is to resort to nonlinearities which impose saturation on the growth of the cubic term at large intensities of the wave field(s). The simplest realization is represented by the following generalized nonlinear Schrödinger equation (NLSE), written in terms of the propagation of an optical wave, with complex amplitude $u\left(z,y;z\right)$ in a bulk waveguide Rypdal :

[TABLE]

where positive constant $U_{0}^{2}$ characterizes the saturation intensity. In fact, one may set $U_{0}=1$, by means of rescaling $u\equiv U_{0}\tilde{u}$, $\left(x,y\right)\equiv\left(\tilde{x},\tilde{y}\right)/U_{0}$, $z\equiv\tilde{z}/U_{0}^{2}$, hence Eq. (5) may be written in the parameter-free form. The nonlinear term in this equation adequately models, in particular, optical properties of warm atomic gases Tikho . The temporal variable does not appear in Eq. (5), as it governs the evolution of the monochromatic wave in the spatial domain KA .

A general form of localized vortex-soliton solutions to Eq. (5) is

[TABLE]

where $\left(r,\theta\right)$ are polar coordinates in the $\left(x,y\right)$ plane, $k>0$ is a real propagation constant, which is a free parameter of the vortex-soliton family, integer $S$ is the vorticity (winding number), and real amplitude $U(r)$ obeys the ordinary differential equation,

[TABLE]

supplemented by the above-mentioned boundary condition at $r\rightarrow 0$,

[TABLE]

and subject to the condition of the exponential localization at $r\rightarrow\infty$:

[TABLE]

(pre-exponential factor $r^{-1/2}$ in Eq. (9) is essentially the same as in the standard asymptotic approximation for cylindrical functions at $r\rightarrow\infty$).

However, both the experiment Tikho and numerical analysis Tikho ; Dima1 have demonstrated that all the vortex solitons generated by Eqs. (6)-(9) are unstable. In particular, the simplest one, with $S=1$, spontaneously splits in two fragments, which are close to fundamental (zero-vorticity) solitons.

II.2 The stabilization of two- and three-dimensional vortex solitons

by cubic-quintic and other competing nonlinearities

II.2.1 Formulation of the model

The NLSE of the cubic-quintic (CQ) type appears if the saturable nonlinearity in Eq. (5) is replaced by its truncated expansion:

[TABLE]

Equation (10) is written in the spatiotemporal  form, including the GVD term (its sign corresponds to the anomalous dispersion, because normal dispersion cannot support bright spatiotemporal solitons KA ), $\tau$ being the same reduced time as in Eqs. (1) and (2). All free coefficients are scaled out in the normalized form of Eq. (10).

The CQ terms, as written in Eq. (10), adequately model the nonlinear response in specific optical media. In particular, the use of this nonlinearity in a bulk waveguide filled by liquid carbon disulfide has made it possible to create robust fundamental (zero-vorticity) 2D spatial solitons Cid-robust . Further, a recently developed technique makes it possible to design desirable coefficients of cubic, quintic, and higher-order nonlinearities in colloidal suspensions of metallic nanoparticles Cid-review . In particular, stable 2D fundamental solitons were created in a medium with an effective quintic-septimal nonlinearity, where the usual cubic term was negligible Cid2 .

The 1D reduction of Eq. (10), for solutions which do not depend on $x$ and $y$, admits well-known exact soliton solutions Pushkarovs ; Enns ,

[TABLE]

which exist and are completely stable in interval

[TABLE]

3D vortex-soliton solutions are produced by the substitution of a generalization of ansatz (6),

[TABLE]

in Eq. (10), where $U(r,\tau)$ is a real even function of $\tau$ obeying an accordingly modified version of Eq. (7):

[TABLE]

In the 3D case, boundary conditions (8) and (9) are supplemented by the condition of the exponential localization of the solution at $|\tau|\rightarrow\infty$: $U\sim\exp\left(-\sqrt{2k}|\tau|\right)$, while being subject at $r\rightarrow 0$ to the same condition (8) as in 2D.

II.2.2 The shape and stability of 2D vortex solitons

Vortex solitons in the 2D version of the NLSE with the CQ nonlinearity were introduced in work Drits . The existence of stable vortex solitons with $S=1$ in this model was revealed by work Manolo .

Numerical solution of Eq. (14) without the temporal derivative produces vortex-soliton radial profiles with the propagation constant belonging to the same interval (12) as in the 1D case, because Eq. (14) takes an asymptotically 1D form at $r\rightarrow\infty$. In the limit of $k\rightarrow 3/16$, the 2D solitons feature a “flat-top” shape in a circular area of large radius,

[TABLE]

It is filled by the solution with a nearly constant amplitude, close to the largest possible value, $U_{\max}=\sqrt{3}/2$ (it is a constant solution of Eq. (14) for $k=3/16$), with a “hole” of radius $r_{0}$ produced by the vorticity around the pivot (placed at $r=0$), see Figs. 1(a) in Manolo and 4-6 in Pego . The hole’s radius may be estimated by means of the Thomas-Fermi (TF) approximation, which neglects derivatives in Eq. (14), thus yielding two solutions: $U=0$, and one produced by equation

[TABLE]

Obviously, relevant roots of Eq. (16), with real $U^{2}>0$, do not exist at

[TABLE]

which determines the radius of the vorticity-drilled “hole” in the TF approximation (at $r<r_{0}$, this approximation admits solely the $U=0$ solution). For $k\rightarrow 3/16$, Eq. (17) yields $r_{0}\approx 2\sqrt{2}S$, which is indeed much smaller than the outer radius (16) of the flat-top-shaped soliton. The TF approximation can be applied to vortex solutions in other contexts too Fetter .

The stability of the 2D vortex solitons was for the first time addressed, in the framework of Eq. (10) (without the time-derivative term) in work Manolo by means of direct simulations. It was found that a part of the vortex-soliton family with $S=1$ is stable, and, furthermore, such solitons, if set in motion, may collide quasi-elastically. In a rigorous form, the stability of the vortex-soliton solutions was investigated in work Pego , through numerical computation of stability eigenvalues for small perturbations added to the stationary states. Thus, stability regions in the existence interval (12) were identified for all values of winding number $S$, in the form of

[TABLE]

(all fundamental solitons with $S=0$ are stable). The numerically found stability boundaries are Pego :

[TABLE]

Thus, with the increase of the winding number from $1$ to $5$ the relative width of the stability regions, with respect to the existence interval (12), drops from $\approx 21\%$ to $\approx 4\%$.

The 2D version of the same NLSE (10) can be realized in optics as the equation for the spatiotemporal propagation in a planar waveguide, which implies dropping coordinate $y$ in Eq. (10), while keeping temporal variable $\tau$. Thus one may define spatiotemporal vortices spatiotemp-vort1 ; spatiotemp-vort2 as modes given by ansatz (6) in which $\left(r,\theta\right)$ are realized as the polar coordinates in the $\left(x,\tau\right)$ plane. The spatiotemporal vortex soliton may be excited in the planar waveguide by an input whose amplitude and phase are patterned according to Eq. (6), taken at $z=0$ spatiotemp-vort1 .

Lastly, it is possible to introduce a two-dimensional NLSE/GPE with specially selected coefficients depending on the radial coordinate, which makes it possible to produce vortex solitons with various more sophisticated shapes NonlinDyn ; NonlinDyn2 .

II.2.3 3D vortex solitons

In the framework of 3D NLSE (10), the stability of vortex solitons, shaped as vortex tori, was investigated in work nine . In this case too, 3D solitons exist in interval (12). The main numerical finding is that they are stable in a part of this interval, defined as in Eq. (18), with $k_{\min}^{\mathrm{(3D)}}(S=1)\approx 0.13$ nine . No stability interval was found for 3D vortex solitons with $S\geq 2$. Unstable vortex tori spontaneously split into fragments, which seem as fundamental ($S=0$) solitons, see Fig. 1.

II.2.4 Stable vortex solitons supported by competing quadratic and

cubic nonlinearities.

As mentioned above, systems of equations (1) and (2), which model the co-propagation of fundamental-frequency and second-harmonic waves, coupled by $\chi^{(2)}$ terms, produce stable solitons solely with zero vorticity, but no stable vortex solitons. Two-component vortex solitons in such a system (reduced to two equations by substitution (4)) may be stabilized if it includes the repulsive cubic self-interaction. In 2D, this was demonstrated for vortex solitons with $S=1$ and $2$, with the relative size of the stability regions $\approx 8\%$ and $\approx 5\%$, respectively S=1and2chi2chi32D . Small stability regions, with the respective size $\approx 3\%$ and $\approx 1.5\%$, were found for vortex solitons with $S=3$ and $4$ as well S=3and4chi2chi32D . Finally, 3D vortex solitons with $S=1$ may also be stable in the system with the competing quadratic and self-defocusing cubic interactions, the corresponding relative size of the stability domain being $\approx 10\%$ , while all the 3D solitons with $S\geq 2$ are unstable, as well as in the CQ model chi2chi33D .

II.2.5 Interactions between multidimensional solitons

Interactions between stable 2D and 3D solitons have also been theoretically investigated. In particular, an effective potential for the interaction between far separated 2D identical solitons with chemical potential $\mu<0$, which carry vorticity $S$, was derived in Ref. interaction :

[TABLE]

where $R$ and $\delta$ are the distance and phase shift between the solitons, and $C$ is a positive constant (the interaction potential was derived in Ref. interaction in a more general form, which applies as well to models based on complex Ginzburg-Landau equations (CGLEs), that include dissipation and gain, see Eq. (101) below). Interactions between vortex solitons were investigated by means of direct simulations too SKA-int1 ; SKA-int2 .

II.3 A new realization of competing nonlinearities: Stable 2D and 3D

vortex quantum droplets with embedded vorticity

II.3.1 Gross-Pitaevskii equations with the Lee-Huang-Yang (LHY)

corrections

Recent theoretical and experimental works with binary BEC have revealed a fascinating possibility of the creation of stable 3D and 2D soliton-shaped objects, which are stabilized by the LHY corrections to the mean-field (MF) theory, which account for the effect of quantum fluctuations around MF states LeeHY . An appropriate system is provided by a binary BEC with self-repulsion in each component (in this case the LHY correction is relevant) and the MF attraction between the components. In this setting, the normalized GPE for identical wave functions of the two components in 3D is Petrov-QD

[TABLE]

where the self-focusing cubic term implies that the cross-attraction is slightly stronger than self-repulsion, while the defocusing quartic term represents the LHY corrections. Similar to other models with competing attractive and repulsive nonlinearities, Eq. (21) readily generates a family of stable fundamental (zero-vorticity) soliton-like modes, which are often called “quantum droplets” (QDs), as they represent self-trapped modes filled by an ultradilute superfluid Petrov-QD ; Astra . The creation of stable 3D fundamental solitons, and collisions between moving ones, by GPE which includes both the LHY term and an additional self-repulsive quintic one was recently addressed in Ref. SKA .

Soon after being predicted, QDs were experimentally created in a binary BEC of 39K atoms. The experimentally observed shape may be nearly two-dimensional, if the condensate is strongly confined in one direction Tarr1 ; Tarr2 , or isotropic, if the trapping potential is weak Ing ; Ing2 .

The reduction of the effective dimension from 3 to 2 by means of a tightly binding potential applied in one direction entails the replacement of the cubic-quartic nonlinearity in Eq. (21) by the single cubic term multiplied by a logarithmic factor Astra :

[TABLE]

The Hamiltonian corresponding to this equation is

[TABLE]

where $e=\allowbreak 2.\,\allowbreak 718...$ is the base of natural logarithms.

II.3.2 Vortex solitons in the 2D setting

Obviously, the nonlinearity in Eq. (22) is attractive at small values of the normalized density, $|\psi|^{2}<1$, for which the logarithm is negative, and it reverses its sign to repulsive at $|\psi|^{2}>1$, i.e., this nonlinearity seems as a competing one. In particular, similar to the above-mentioned model with the CQ nonlinearity, Eq. (22) gives rise to broad 2D solitons with the flat-top shape, featuring a nearly constant value of the density in its intrinsic area determined by minimization of the energy density, as defined by Eq. (23): $n_{\mathrm{flat-top}}=1/\sqrt{e}\approx 0.6065$. In this case, a crude approximation for the radius of the inner “hole” created by the embedded vorticity is provided by the TF approximation, which neglects term $\nabla_{\mathrm{2D}}^{2}\psi$ in Eq. (22): $\left(r_{0}\right)_{\mathrm{TF}}=\sqrt{e/(2-\sqrt{e})}S\approx 2.8S$, cf. Eq. (17). Comparison with numerical findings demonstrates that this approximation produces reasonable results for $1\leq S\leq 4$ GZ-ln .

Numerical and approximate analytical consideration of Eq. (22) has produced completely stable families of fundamental ($S=0$) solitons, and partly stable families of solitary vortices, up to $S=5$ GZ-ln . They chiefly feature the flat-top shape, although become closer to narrow annuli near the stability boundary. Typical examples of the vortex solitons are displayed in Fig. 2. They are stable provided that their norm exceeds a threshold value, $N>N_{\mathrm{thr}}$(actually, it is the double norm, defined for the binary modes with equal components). The vortex solitons are unstable against spontaneous splitting at $N<N_{\mathrm{thr}}$. For $1\leq S\leq 5$, values $N_{\mathrm{thr}}$ are collected in Table I. For $S$ large enough (actually, for $S\geq 3$), the dependence of $N_{\mathrm{thr}}$ on $S$ can be predicted in an approximate analytical form, which is based on comparison of values of the energy (Hamiltonian, see Eq. (23)) of the unsplit vortex and a set of far separated fragments produced by the splitting. The instability does not occur if the splitting would imply the increase of the energy, which yields the following prediction for the stability boundary GZ-ln :

[TABLE]

(in fact, coefficient $6$ is a fitting factor, while the prediction produces scaling $N_{\mathrm{thr}}\sim S^{4}$).

It is relevant to mention that the 2D model, written in an explicitly two-component form, admits vortex states with the hidden vorticity (HV), i.e., identical amplitude profiles and chemical potential of the components but opposite vorticities, $S_{1}=-S_{2}$. It was found in Ref. GZ-ln that the HV modes with $S_{1}=-S_{2}=1$ have their limited stability region too, while they are completely unstable for $S_{1}=-S_{2}\geq 2$.

II.3.3 “Invisible” splitting instability of multiple vortices

As shown above, QD families with multiple vorticity keep their partial stability up to $S=5$. Nevertheless, it was demonstrated in Ref. GZ-ln that all multiple-vortex states, with $S\geq 2$, demonstrate virtually invisible but, strictly speaking, existing structural instability against splitting. This means that a specially selected small perturbation, without initiating any dynamical instability, may split the pivot (phase singularity) of the multiple vortex with topological charge $S$ into sets of $S$ or $S+2$ phase singularities corresponding to unitary vortices, although the splitting remains almost invisible, as it occurs in the broad central “hole” induced by the multiple vorticity, where absolute values of the wave function remain extremely small (in the case of the splitting into $S+2$ phase singularities, one of them belongs to a unitary antivortex, so that the total vorticity is conserved). It is worthy to note that the structural instability is not specific to the model of the LHY type, but is a generic effect, relevant to all models which support dynamically stable higher-order vortex solitons. While this effect seems intuitively obvious, it makes sense to produce it in an explicit analytical form, following Ref. GZ-ln .

To demonstrate the structural instability against splitting, one can place the pivot of the state with $S\geq 2$ at the origin ($x=y=0$), and add a specific perturbation mode carrying its own intrinsic vorticity $s=1$, with small complex amplitude $-\varepsilon$, thus producing a perturbed configuration,

[TABLE]

Phase singularities of unitary vortices, into which the original vortex is split by the perturbation, are identified as zeros of $\left|\psi_{\mathrm{pert}}\left(x,y\right)\right|$.  As seen from Eq. (25), they are located at the origin, $x_{\mathrm{piv}}^{(1)}=y_{\mathrm{piv}}^{(1)}=0$, and at $\left(S-1\right)$ additional points, produced by $\left(S-1\right)$ branches of the root of degree $1/\left(S-1\right)$:

[TABLE]

Another option is to consider a small disturbance with intrinsic vorticity $s=-1$ (instead of $+1$, which was taken above):

[TABLE]

In this case, the pivots are located at points defined by equation $\left(x+iy\right)^{S}=\varepsilon\left(x-iy\right)$, i.e., $\left(x+iy\right)^{S+1}=\varepsilon r^{2}$, which yields a set of $S+2$ roots, viz.,

[TABLE]

plus the central pivot at the origin.

Verification of this consideration by direct simulations of Eq. (22) is displayed in Fig. 3. Strong magnification of the picture observed in the nearly empty “hole” confirms the splitting of the pivot with multiplicity $S$ into $S$ or $S+2$ sets of unitary-vortex pivots, under the action of the small initial perturbation with its intrinsic vorticity $s=+1$ or $s=-1$, respectively (the vortex soliton with $S=1$ does not undergo the splitting). On the other hand, the splitting remains virtually invisible on the normal scale of $\left|\psi\left(x,y\right)\right|$, hence it does not imply actual instability of the solitary states with multiple vorticity.

II.3.4 3D vortex solitons

Numerical solutions of three-dimensional Eq. (21) for solitons with embedded vorticity were presented in Ref. 3DLHY . Being supported by the competing set of the cubic and quartic terms, the solitons again tend to feature a flat-top shape. An analytical estimate for the stability of the 3D vortex soliton against the spontaneous splitting, which is based, as in the 2D case, on the comparison of energies of the unsplit state and the set of far separated fragments, yields scaling $S_{\min}\sim S^{6}$, cf. the 2D result given by Eq. (24). In accordance with this prediction, a conspicuous stability region was found, in a numerical form, for $S=1$, while for $S=2$ only vortex solitons with a very large norm may be stable. Examples of the evolution of stable and unstable 3D QDs with $S=2$ are displayed in Fig. 4.

All the HV states turn out to be unstable in the 3D model against spontaneous splitting, see an example in Fig. 5.

Lastly, it is relevant to mention that the creation of QDs was predicted and experimentally realized in dipolar BECs, making use of the attractive long-range interactions and LHY effect dip1 ; dip2 . However, the theoretical analysis has demonstrated than all QDs with embedded vorticity are unstable in that setting Macri .

II.4 Trapping potentials for vortex solitons

II.4.1 The harmonic-oscillator (HO) trap

The use of trapping potentials, which are natural ingredients of experimentally relevant settings in BEC and optics alike, makes it possible to stabilize both fundamental and vortex solitons. The model is written as the normalized GPE Pit for the mean-field wave function $\psi(x,y,z,t)$:

[TABLE]

Here $U$ is the trapping potential, and the negative sign in front of the cubic term implies self-attraction, as above. Actually, a quasi-two-dimensional trapping potential, i.e., $U=U\left(x,y\right)$, may be sufficient to support stable 3D solitons low-dim ; low-dim2 ; Herve-lowdim .

In optics, time $t$ in Eq. (29) is replaced by the propagation distance, $z$, while the original coordinate $z$ is replaced by the temporal variable $\tau$, cf. Eq. (10). The effective potential in the optical waveguide may be solely two-dimensional, being proportional to the local change of the refractive index, $U\left(x,y\right)\sim-\delta n(x,y)$.

For axially and spherically symmetric potentials, $U=U\left(r,z\right)$, where $\left(r,z,\theta\right)$ is the set of cylindrical coordinates, solutions to Eq. (29) are looked for in the form similar to that defined above in Eq. (13):

[TABLE]

with real chemical potential $\mu$, and real amplitude function $R\left(r,z\right)$ satisfying the stationary equation which is similar to Eq. (14):

[TABLE]

Most relevant for the realization in BEC is the harmonic-oscillator (HO) trapping potential Pit ,

[TABLE]

where $\Omega^{2}$ accounts for anisotropy of the trap, the limits of $\Omega^{2}\gg 1$ and $\Omega^{2}\ll 1$ corresponding, respectively, to nearly 2D (“pancake-shaped” pancake ) and nearly one-dimensional (“cigar-shaped” Randy-NJP ) configurations. Equation (29) can be derived from the Hamiltonian,

[TABLE]

which is a dynamical invariant (conserved quantity) of Eq. (29). Two other dynamical invariants of Eq. (29) are the total norm,

[TABLE]

and $z$-component of the angular momentum,

[TABLE]

whose value in stationary state (30) is

[TABLE]

Note that, unlike the model-specific Hamiltonian, the definitions of the norm and angular momentum, given by Eqs. (34) - (36) or their 2D counterparts are universal.

II.4.2 Lattice (spatially-periodic) trapping potentials

Many experiments in BEC and photonics made use of spatially periodic potentials, which are represented by optical lattices (OLs) for atomic condensates Brazh ; Morsch , by photonic crystals steering the transmission of light in optics PhCr ; JY , or by similar periodic structures acting on exciton-polariton fields in semiconductor microcavities Kriz . In addition to assembling the traditional photonic-crystal structures, similar multi-channel waveguides can be produced by burning multi-core patterns in bulk silica Jena , and, on the other hand, virtual (rewritable) structures can be created in photorefractive materials, illuminating them by counterpropagating laser beams in the ordinary polarization, which interfere linearly and induce an effective lattice, while solitons are built by extraordinarily-polarized beams, which are subject to the action of saturable nonlinearity MSegev .

It is relevant to stress that, while the HO potential, considered above, maintains bound states in the absence of any nonlinearity, just as eigenmodes of this potential in quantum mechanics, no self-trapping is possible in periodic potentials. As the basic model, one can take the 3D GPE with the full 3D OL potential, or its quasi-2D version (as mentioned above, the 2D potential may be sufficient to create stable vortex solitons in the fully 3D form low-dim ; low-dim2 ; Herve-lowdim ):

[TABLE]

Here, $\sigma=1$ or [math], for the full 3D or quasi-2D lattice, respectively, and $g=+1$ and $-1$ corresponds, severally, to the self-focusing and defocusing signs of the cubic nonlinearity. Note that the defocusing nonlinearity may create gap solitons in the presence of the lattice potential Brazh ; Morsch ; Yang-book . In optics, as mentioned above, $t$ is replaced by the propagation coordinate, $z$, while original coordinate $z$ is replaced by the reduced-time variable (3), and solely the 2D potential, with $\sigma=0$ in Eq. 37, is relevant. The sign of the cubic nonlinearity is usually $g=+1$ (self-focusing) in optical media. GPE (37) conserves two dynamical invariants, viz., the total norm/energy, defined as in Eq. (34), and the Hamiltonian,

[TABLE]

cf. Eq. (33).

II.5 Vortex solitons in the HO trap

II.5.1 3D vortex modes and their stability

Before displaying numerical findings, it is relevant to outline approximate analytical results which are available in the present setting. In the linear limit, Eq. (29) is tantamount to the quantum-mechanical Schrödinger equation for the 3D anisotropic HO. In the Cartesian coordinates, the corresponding 3D eigenfunctions are built as

[TABLE]

where $\Phi_{j}$, $\Phi_{k}$ and $\Phi_{l}$ are stationary wave functions of 1D harmonic oscillators with quantum numbers $j,k,l$, which correspond to energy eigenvalues $j+1/2$, $k+1/2$ and $\left(l+1/2\right)\Omega$, respectively, the chemical potential being

[TABLE]

The states which carry over into ones (30) with vorticity $S$ in the nonlinear model are constructed as combinations of factorized wave functions (39) with $l=0$ and $j+k=S$. Since the correction to $\mu$ from the self-attractive nonlinearity is negative, this restriction and Eq. (40) with $l=0$ impose a bound on $\mu$,

[TABLE]

In particular, the eigenfunctions of the linear model, with vorticities $S=1$ and $S=2$, are written, in terms of the cylindrical coordinates, $\left(r,\theta,z\right)$, as

[TABLE]

[TABLE]

Numerical solution of Eq. (31) generates $\mu(N)$ and $E(N)$ dependences for families of vortical trapped modes with $S=1$. They are displayed in Fig. 6 for $\Omega=10$ and $1$, which correspond to the pancake-like and isotropic trapping configurations, respectively (for the cigar-shaped trapped modes, that correspond, e.g., to $\Omega=0.1$, the dependences, which are not shown here, are very close to those for $\Omega=1$ Dum3D ). An obvious feature of the figure is the presence of a largest norm, corresponding to the turning point of the $\mu(N)$ curve, which bounds the existence of the trapped modes. This limitation is caused by the possibility of the supercritical collapse in the 3D model Berge ; SulemSulem ; Gadi : if the norm is too large, the self-attraction cannot be balanced by the gradient part of energy (33). The largest norm in this model was numerically found, as a function $\Omega$, in Ref. Sadhan . Another implication of the possibility of the supercritical collapse is that all the stable modes found in the 3D model are actually metastable ones. They cannot play the role of the GS, which does not exist in the model admitting the supercritical collapse HP .

The stability of stationary solutions produced by Eq. (31) was identified through the computation of eigenvalues $\lambda$ of small perturbations Dum2D ; Dum3D . To this end, a perturbed solution to Eq. (29) is looked for as

[TABLE]

where $\left(u,v\right)$ are eigenmodes of infinitesimal perturbations corresponding to integer values of azimuthal index $L$, and $\ast$ stands for the complex conjugate. The substitution of ansatz (44) in Eq. (29) and linearization lead to the Bogoliubov - de Gennes equations Pit ,

[TABLE]

supplemented by the boundary conditions demanding that $u(r,z)$ and $v(r,z)$ decay exponentially at $r\rightarrow\infty$ and $\left|z\right|\rightarrow\infty$, and decay as $r^{\left|S\pm L\right|}$ at $\rho\rightarrow 0$. The underlying stationary mode is stable if all eigenvalues $\lambda$ produced by numerical solution of Eq. (45) Ueda are pure imaginary. Results of the numerical stability analysis are included in Fig. 6, where stable and unstable portions of the vortex-soliton families are designated.

Note that the stability region of the vortex mode with $S=1$ is essentially smaller than formally predicted by the celebrated Vakhitov-Kolokolov (VK) criterion, $d\mu/dN<0$, which provides a necessary, but generally, not sufficient condition for the stability of self-trapped modes VK ; Berge ; Gadi (it cannot detect instabilities accounted for by complex eigenvalues $\lambda$ – in particular, the splitting instability of the modes with $S\geq 1$). In the region where the VK criterion holds, but the vortices are unstable, they are vulnerable to the splitting instability induced by perturbations with $L=2$ and $3$ in Eq. (44). The vortex modes with $S\geq 2$ were found to be completely unstable Dum2D ; Dum3D .

The predictions for the stability based on the computation of the stability eigenvalues were verified by direct simulations of Eq. (29), starting with stationary modes to which small arbitrary perturbations were added. The robustness of stable vortices is illustrated by Fig. 7, which demonstrates that they absorb the perturbations and clean themselves up.

Lastly, it is relevant to mention the consideration of a binary model with the repulsive cubic nonlinear terms (both intra- and inter-component ones) reported in Ref. Radik . By means of systematic simulations and a semi-analytical approximation based on a finite-mode Galerkin decomposition of the wave functions, nontrivial two-component stationary states and dynamical regimes have been found in that system.

II.5.2 2D trapped vortex modes

The one-component setting

As said above, taking the limit of $\Omega\rightarrow\infty$ in HO potential (32) leads to the reduction of the 3D model to its 2D version, with Eqs. (GPE), (31), (33), (34), (35), and (45) carrying over into their 2D counterparts, and the effective 2D HO potential taking the form of $U\left(x,y\right)=(1/2)\left(x^{2}+y^{2}\right)$. In particular, stationary solutions of the 2D version of Eq. (29) are looked for as $\psi=\exp\left(-i\mu t+iS\theta\right)R(r)$, cf. Eq. (30). In the limit form of the linear equation, the chemical potential is given by the usual eigenvalue for the 2D linear Schrödinger equations with this potential:

[TABLE]

cf. Eq. (41).

The shape and stability of nonlinear 2D modes trapped in the HO potential was studied in a number of works 2D ; 2D2 ; 2D3 ; 2D4 ; Dum2D . The results, produced by the computation of stability eigenvalues in the framework of the 2D Bogoliubov - de Gennes equations (45), demonstrate that the family of the fundamental ($S=0$) trapped modes is stable in its entire existence region,

[TABLE]

where $N_{\mathrm{TS}}$ is the numerically found value of the norm of the *Townes soliton *(TS), which determines the threshold for the onset of the critical collapse in the 2D NLSE Berge ; SulemSulem ; Gadi . In the absence of trapping potentials, the family of TSs is degenerate, as, at all values of $\mu$, they assume the single value of the norm, which is exactly equal to $N_{\mathrm{TS}}$. The degeneracy is a consequence of the specific scaling invariance of the NLSE in 2D Berge ; SulemSulem ; Gadi . The trapping potential introduces a characteristic spatial scale (the usual HO length, or the spatial period, in the case of lattice potentials), which breaks the scaling invariance and thus lifts the degeneracy, making $N$ a function of $\mu$. In fact, Eq. (47) demonstrates that the degeneracy is lifted so that the soliton’s norm falls below the collapse-onset threshold. This circumstance lends the trapped modes with $S=0$ protection against the collapse, i.e., stability, and actually makes them the system’s GS, which did not exist in the absence of the trapping potential. This is a major difference from the stabilization mechanism provided by the trapping potential in 3D, where, as said above, the GS cannot exist (as the supercritical collapse occurs at any value of the norm), only metastability of the trapped modes being possible.

The family of TSs with embedded vorticity $S=1$ is also degenerate in the free space, admitting a single value of the norm,

[TABLE]

Minsk2 (it is relevant to mention that the first generalization of the TS, in the form of higher-order radial states with $S=0$, was introduced in Ref. Yankauskas , soon after the concept of the TS was established; while such states are unstable, they can be stabilized by means of the “management” technique, which makes the coefficient of the cubic nonlinearity a periodically varying function of $t$ VPG-management ). Recently, an analytical approximation was developed, which produces values $N_{\mathrm{TS}}^{(S)}$ for $S\geq 1$ with a good accuracy Jieli2 , see Eq. (94) and Table II below.

The HO trapping potential stabilizes the modes with $S=1$ in interval

[TABLE]

the corresponding stability region in terms of the chemical potential being $1.276\equiv\mu_{\mathrm{cr}}<\mu\equiv\mu_{0}^{\mathrm{(2D)}}(S=1)\equiv 2$ (the right edge of the region is determined by Eq. (46)) Dum2D . Trapped 2D vortices with $S\geq 2$ they remain completely unstable.

The partial stability of the family of trapped vortex modes with $S=1$, predicted through the computation of the corresponding eigenvalues, was corroborated by direct simulations of the perturbed evolution. The simulations have also revealed a noteworthy dynamical regime for the trapped vortices with $S=1$ in interval

[TABLE]

adjacent to one given by Eq. (49). In interval (50), the evolution of the unstable vortex is time-periodic, as shown in Fig. 8: it spits into two fragments which then recombine back into the vortex, keeping the vorticity of the configuration in the course of the cycles, while the splitting orientation slowly rotates in the $\left(x,y\right)$ plane. The vortices with still larger values of the norm, $N>10.30$, also split in two fragments, which, however, fail to recombine. Instead, each one quickly blows up, i.e., collapses.

A similar scenario of the instability development of 3D vortex solitons with embedded vorticity $S=1$, trapped in the three-dimensional HO potential, was revealed by simulations reported in Refs. Ueda ; Ueda2 .

Two-component systems

A noteworthy generalization of the above analysis was performed for a system of two nonlinearly coupled fields, which may be realized as a binary BEC, or as co-propagation of two optical beams in a bulk waveguide Brtka . The respective coupled two-dimensional GPE/NLSE system is

[TABLE]

[TABLE]

In addition to the symmetric states with equal vorticities of both components, this system gives rise to composite modes with the HV (hidden vorticity), formed by the components with equal norms $N/2$ and opposite vorticities, $S_{1,2}=\pm 1$. The HV modes are stable, roughly, in the region of $N<7$, $-1<\eta<0.2$. A related problem is the study of compound states in which one component is fundamental ($S=0$) and the other carries vorticity VPG-mixed .

Another relevant two-component system models the dynamics of BEC in parallel 2D layers, coupled by hopping of atoms between them, with the HO trapping potential acting in each of them Viskol :

[TABLE]

The asymmetry of the binary vortex state is characterized by parameter

[TABLE]

A standard diagram for the symmetry-breaking bifurcation (transition to an asymmetric state) is usually displayed in the form of the $\theta(N)$ dependence. For the model based on Eqs. (53a), it is plotted in Fig. 10, for $\kappa=0.15$. It is seen that, with the increase of the total norm, the symmetric vortex becomes unstable at a critical (bifurcation) point, $N_{\mathrm{cr}}^{(S=1)}$. The bifurcation gives rise to two mutually symmetric branches of stable asymmetric vortex states, characterized by dependence $\pm\theta(N)$, a typical example of which is shown in Fig. 10 (in the figure, only the branch with $\theta>0$ is presented). Numerical findings presented in Ref. Viskol suggest a dependence of $N_{\mathrm{cr}}^{(S=1)}$ on the linear-coupling constant $\kappa$, which may be well fitted by a simple linear relation,

[TABLE]

At still larger values of $N$, the asymmetric states, with one dominating component, are destroyed by instability which is, generally, similar to that outlined above for the single-component model. In Fig. 10, solitons belonging to the dashed portion of the asymmetric branch are subject to the latter instability.

II.6 Vortex solitons supported by a spatially periodic (lattice)

potential

Creation and stabilization of solitons by means of lattice potentials was addressed in many works, see reviews Brazh ; Morsch ; big ; Vyslo1 ; Vyslo2 ; Yukalov ; RMP ; Entropy , therefore these results are presented here in a brief form. The possibility to stabilize 2D solitons with embedded vorticity by means of the lattice potential was first demonstrated in Refs. lattice1 and lattice2 ; lattice3 . As shown in Figs. 11 and 12, these vortex states do not seem as familiar annulus-shaped (alias crater-shaped) objects (cf. Figs. 1(a), 2(a1-a4), and 8(a,c)), but rather as circular chains of local peaks, in terms of the local density. Vorticity $S$ is represented by the phase structure of these modes, which feature the phase circulation of $2\pi S$ corresponding to a round trip along a path encircling the mode’s pivot. It is seen that the basic structure of the lattice vortex with $S=1$ is represented by the set of four density peaks, of two different types: rhombic, with a nearly empty central site, which is displayed in Fig. 11, and a densely packed square-shaped one, which does not include an empty central site, in Fig. 12. These two types of lattice vortices are usually called on-site-centered and off-site centered ones, the latter type usually being essentially more stable.

In addition to the basic on- and off-cite-centered vortex solitons with the natural symmetric shape, the analysis has also revealed a possibility of the existence of asymmetric ones, built as triangular or stretched rhombic and rectangular multi-peak patterns Asymm . A related result is a possibility of replication of a vortex mode, originally created in one potential well of a 2D double-well configuration, by building a twin vortex in the adjacent well replication .

In the limit of a very deep lattice potential, i.e., $\varepsilon\rightarrow\infty$ in Eq. (37), the 2D version of Eq. (37) carries over into the discrete NLSE PGK-book :

[TABLE]

where $\left(m,n\right)$ is the set of discrete coordinates replacing $\left(x,y\right)$, and real $C>0$ is a coefficient of the intersite coupling in the discrete lattice. Stationary states are represented by solutions to Eq. (56) in the form of of

[TABLE]

where the stationary lattice field $u_{m,n}$ is real for fundamental discrete solitons, and complex for vortex solitons discrvort . An example of a numerically generated stable vortex discrete soliton with $S=1$ is displayed in Fig. 13. In discrete models, on-site-centered vortex solitons also tend to be essentially more stable than their off-site-centered counterparts discrvort .

The discrete NLSE supports vortex solitons with higher winding numbers, $S\geq 2$, but they are completely unstable discrvort . Getting back to the NLSE (37), stable higher-order vortex solitons were produced in Ref. HS-EPL . An example of a stable circular-chain soliton with $S=4$, built of $12$ local peaks, is shown in Fig. 14.

Furthermore, the 2D version of Eq. (37) supports more complex supervortex complexes with two independent vorticities, local one $s$ and global $S$ HS-EPL (see also Ref. Radik2 ). As shown in Fig. 15, the supervortex is built as a circular chains of compact local vortices with winding number $s=1$, each squeezed into one cell of the lattice potential, and global vorticity $S$ imprinted onto the chain. Thus, the supervortices with global vorticities $\pm S$ and fixed local one $s$ are* different states*. In particular, for $\varepsilon=10$ in Eq. (37) and $s=1$, the supervortices are stable for $S=\pm 1$ and $\pm 2$, unstable for $|S|>3$, and marginally stable for $S=\pm 3$ discrvort .

Lastly, the full 3D form of the model based on Eq. (37) with the quasi-2D lattice potential, $\sigma=0$, gives rise to stable 3D solitons with embedded vorticities $S=1$ and $2$, which are also built as chains of local soliton-like objects in the $\left(x,y\right)$ plane, each being self-trapped in the transverse direction, see an example in Fig. 16.

III Stabilization of 2D and 3D semi-vortex and mixed-mode solitons by

the spin-orbit coupling (SOC)

This section aims to summarize recent findings which put forward a completely new approach to the stabilization of vorticity-carrying 2D and 3D solitons in free space (without the use of any external potential), realized in models of the binary (pseudo-spinor) atomic BEC with the SOC. The presentation is chiefly based on Refs. we , HP , and Sherman2 .

III.1 The models

A great deal of attention has been lately drawn to the use of ultracold quantum gases, both bosonic and fermionic, as simulators of various fundamental effects that were previously predicted and/or discovered experimentally in much more complex settings of condensed-matter physics simulator . In particular, much interest has been recently attracted to the implementation of the (pseudo-) SOC in atomic BEC, as an efficient emulation of the fundamental SOC in semiconductors, where the direct SOC is induced by the coupling of the electron’s magnetic moment to the magnetic field generated by the intrinsic electrostatic field of the underlying ionic lattice, in the reference frame moving along with the electron Dresselhaus ; Rashba .The experimental implementation of the pseudo-SOC was proposed Campbell and realized in the condensate of 87Rb atoms, using appropriately designed laser illumination and magnetic fields socbec ; socbec2 ; socbec3 ; socbec4 . Parallel to the experiments, many theoretical studies on this topic have been carried out theory-SOC -Fukuoka2 , we ; HP , see also reviews rf:1 -rf:15 . The SOC emulation in the atomic condensate is provided by mapping the spinor wave function of semiconductor electrons into a two-component pseudo-spinor wave function of the binary BEC composed of atoms in two different hyperfine states. Namely, a pair of states of the 87Rb atom, $\left|\psi_{+}\right\rangle=\left|F=1,m_{F}=0\right\rangle$ and $\left|\psi_{-}\right\rangle=\left|F=1,m_{F}=-1\right\rangle$, were used to map the spin-up and spin-down electron’s wave function into them socbec ; socbec2 ; socbec3 ; socbec4 . Thus, the dynamics of the fermionic wave functions of electrons in the semiconductor may be emulated by the mean-field dynamics of the bosonic gas.

The consideration of the interplay of the SOC, which is, essentially, linear mixing between the two components of the spatially inhomogeneous binary BEC, and the intrinsic nonlinearity in the bosonic condensate has made it possible to predict diverse nonlinear patterns strongly affected or created by the SOC, including 1D solitons Konotop -Konotop5 , 2D gap solitons supported by OL potentials gap-sol , and 2D vortices and vortex lattices, in forms specific to the spin-orbit-coupled BEC Fukuoka -Fukuoka10 .

The 2D model of the binary SOC BEC is based on the following system of coupled GPEs for two components of the pseudo-spinor wave function,$\Psi\equiv\left\{\psi_{+},\psi_{-}\right\}$ we :

[TABLE]

where $\sigma_{x,y,z}$ are the Pauli matrices, $\lambda$ is a real coefficient of the SOC of the Rashba type (a combination with the SOC of the Dresselhaus type, which is modeled by combination $\sigma_{x}\partial/\partial x-\sigma_{y}\partial/\partial y$, instead of $-\sigma_{y}\partial/\partial x+\sigma_{x}\partial/\partial y$ in Eq. (62), is not considered here, as it tends to destroy 2D solitons Sherman2 ), $\eta$ is the relative strength of the cross-attraction between the components (cf. Eq. (51a)), while the strength of the self-attraction is normalized to be $1$. Coefficient $1/\lambda$ has the dimension of length, defining a fixed scale which breaks the scale invariance of the NLSE/GPE in the free 2D space and thus makes it possible, as mentioned above, to create stable solitons with norms falling below the collapse threshold.

Stationary solutions of Eq. (62) for 2D solitons with real chemical potential $\mu$ are looked for as $\psi_{\pm}=\exp\left(-i\mu t\right)u_{\pm}\left(x,y\right)$, where complex stationary wave functions are determined by equations

[TABLE]

Dynamical invariants of Eqs. (62) are the same total norm $N$ as defined above, Hamiltonian, and linear momentum:

[TABLE]

[TABLE]

The consistent derivation of the effective 2D SOC model from the full 3D system of GPEs may give rise to the 2D equations with nonpolynomial nonlinearity, as a generalization of the cubic terms in Eq. (62) Wesley . Such a generalized system also creates stable solitons in the 2D free space. Another relevant generalization addresses a model of a dual-core nonlinear coupler in optics, where SOC is emulated by temporal dispersion of the linear inter-core coupling Kart . It was demonstrated that this model gives rise to stable 2D spatiotemporal solitons.

The 3D model is taken here with the SOC of the Weyl type HP :

[TABLE]

where $\lambda$ is again the SOC coefficient, and the 3D matrix vector is ${\boldsymbol{\sigma=}}\left\{\sigma_{x},\sigma_{y},\sigma_{z}\right\}$. The 3D system conserves the norm and linear momentum, along with the Hamiltonian,

[TABLE]

III.2 Stable 2D solitons: quiescent and mobile semi-vortices and mixed

modes

III.2.1 2D semi-vortices

Unlike the models considered above, which give rise to fundamental states with zero vorticity, a vortical component is inherently present in any self-trapped state generated by the nonlinear SOC systems. First, Eq. (62) admits stationary solutions with real chemical potential $\mu$, written in terms of the polar coordinates, $\left(r,\theta\right)$:

[TABLE]

where real functions $f_{1,2}\left(r\right)$ take finite values and have zero derivatives at $r=0$, and feature the following asymptotic form at $r\rightarrow\infty$:

[TABLE]

with constants $F$ and $\delta$. As it follows from Eq. (74), the solutions may be exponentially localized at

[TABLE]

Solutions (73) are built as bound states of a fundamental (zero-vorticity, $S_{+}=0$) soliton in component $\psi_{+}$ and a solitary vortex, with vorticity $S_{-}=1$, in $\psi_{-}$, therefore composite modes of this type are called semi-vortices (SVs) we . The invariance of Eq. (29) with respect to transformation

[TABLE]

gives rise to a conjugate semi-vortex, which is a mirror image of (73), with $\left(S_{+}=0,S_{-}=1\right)$ replaced by $\left(S_{+}=-1,S_{-}=0\right)$:

[TABLE]

Numerically, stable SVs can be readily generated, as solutions to Eq. (29), by means of imaginary-time simulations im-time -im-time4 , starting from the Gaussian input,

[TABLE]

where $A_{1,2}$ and $\alpha_{1,2}>0$ are real constants. A typical example of the SV is displayed, by means of cross sections of its components, in Fig. 17(a).

Further, Fig. 17(b) represents the family of the SVs, showing their chemical potential as a function of the norm. Note that the $\mu(N)$ dependence satisfies the VK criterion, $d\mu/dN<0$, which is the above-mentioned necessary condition for the stability of solitary modes supported by the self-attractive nonlinearity VK ; Berge ; Gadi . The family of the SV solitons exists precisely in the interval of norms (47), which, as said above, should secure their stability against the critical collapse. It is also worthy to note that there is no finite minimum (threshold) value of $N$ necessary for the existence of the SVs in the free space. In the limit of $\mu\rightarrow-\infty$, the vortex component of the SV vanishes, while the fundamental one degenerates into the usual TS, with $N=N_{\mathrm{TS}}$, as shown by means of the dependence of ratio $N_{+}/N$ on $N$ in Fig. 17(c).

Systematic real-time simulations confirm the stability of the whole SV family at $\eta\leq 1$, while they are unstable at $\eta>1$ we , where, however, there is another family of stable solitons in the form of mixed modes (MMs), see below. In fact, the SVs at $\eta\leq 1$ and MMs at $\eta\geq 1$ are the first ever found examples of stable solitons supported by the cubic self-attractive nonlinearity in the free 2D space.

The study of 2D SV solitons was extended to a system with long-range anisotropic cubic interactions, mediated by dipole-dipole forces in a bosonic gas composed of atoms carrying magnetic moments Raymond . An essentially novel feature found in the nonlocal model is spontaneous shift of the pivot of the vortical component ($\psi_{-}$) with respect to its zero-vorticity counterpart ($\psi_{+}$). Those solitons also demonstrate a mobility scenario which is essentially different from the one outlined below for the present model: they respond to an applied kick by drift in the opposite direction (i.e., with an effective negative mass) along a spiral trajectory.

III.2.2 2D mixed modes

Aside from the SVs, the same SOC system (62) gives rise to another type of vorticity-carrying solitons, in the form of MMs, which combine terms with zero and unitary vorticities, $\left(S=0,S=-1\right)$ and $\left(S=0,S=+1\right)$, in the spin-up and spin-down components, $\psi_{+}$ and $\psi_{-}$. Numerically, the MM can be produced by imaginary-time simulations initiated by the following input:

[TABLE]

In fact, the MM may be considered as a superposition of the SV (73) and its mirror image (77). Accordingly, symmetry reflection (76) transforms the MM into itself.

A typical example of the MM and the respective $\mu(N)$ dependence are displayed in Figs. 18(a) and (b). Note that peak positions of the two components, $|\psi_{+}\left(x,y\right)|$ and $|\psi_{-}\left(x,y\right)|$, in this state are separated along $x$, Fig. 18(c) showing the separation ($\mathrm{DX}$) as a function of the norm. For a small amplitude of the vortex component, $A_{2}$, Eq. (79) yields $\mathrm{DX}\approx A_{2}/\left(\alpha_{1}A_{1}\right)$.

The $\mu(N)$ dependence for the MM family shows in Fig. 18(b) that the VK criterion holds in this case too, and, as well as SVs, the MMs do not have any threshold value of $N$ necessary for their existence. The family exists in the interval of $N<\tilde{N}_{\mathrm{TS}}(\eta)=2N_{\mathrm{TS}}/(1+\eta)$, where $N_{\mathrm{TS}}$ is the same critical (TS) norm as in Eq. (47). In the limit of $N\rightarrow\tilde{N}_{\mathrm{TS}}(\eta)$ the vortex components vanish in the MM, and it degenerates into a two-component TS, cf. the above-mentioned degeneration of the SV in the limit of $N\rightarrow N_{\mathrm{TS}}$. Separation $\mathrm{DX}$ between peaks of the two components vanishes in this limit too, see Fig. 18(c).

Direct simulations demonstrate that the MMs are unstable at $\eta<1$, and stable at $\eta\geq 1$, i.e., precisely in the regions where the SVs are, severally, stable and unstable we (exactly at $\eta=1$, both the SV and MM solutions are stable Fukuoka2 ). The stability switch between the SV and MM is explained by the comparison of energy (65) for them at equal values of the norm: the energy is smaller for the SV at $\eta<1$, and for the MM at $\eta>1$ we . Accordingly, the SV and MM realize the system’s stable GS, respectively, at $\eta<1$ and $\eta>1$.

Lastly, it is relevant to mention that 2D solitons of the MM type, additionally stabilized by the LHY correction to the model of the binary condensate with SOC, were considered too SOC-LHY .

III.2.3 Mobility of stable 2D solitons

Although the underlying system (62) conserves the momentum (66), the SOC terms break the Galilean invariance of the original NLSEs. For this reason, generating moving solitons from quiescent ones, which were considered above, is a nontrivial problem. As shown in Ref. we , the system gives rise to the mobility of 2D solitons of the MM type along the $y$ axis, but not along $x$ (note that the essential anisotropy of the MM modes with respect to the $x$ and $y$ directions manifests itself by the splitting of peaks of the two components along the $x$ direction, see Figs. 18(a) and (c)). The respective solutions moving at velocity $v_{y}$ can be looked for as

[TABLE]

(Eq. (80) produces the Galilean transform of the wave functions $\psi_{\pm}$ in Galilean-invariant systems). The substitution of ansatz (80) into Eq. (62) leads to the coupled GPEs in the moving reference frame, which differ from Eqs. (62) by the presence of linear mixing between the two components we :

[TABLE]

(here, $\nabla^{2}\equiv\partial^{2}/\partial x^{2}+\partial^{2}/\partial\left(y^{\prime}\right)^{2}$).

Stationary solutions to equations (81) can be obtained, as well as in the case of Eqs. (62), by means of the imaginary-time-evolution method. In particular, at $\eta=2$, when the GS is represented by the quiescent MM soliton, its moving version, which is displayed in Figs. 19(a,b) for $N=3.1$ and $v_{y}=0.5$, exists and is stable too. As well as its quiescent counterpart, this mode features the mirror symmetry between the profiles of $|\phi_{+}\left(x,y\right)|$ and $|\phi_{-}\left(x,y\right)|$. Figure 19(c) shows the amplitude of the moving soliton, $A=\sqrt{|\phi_{+}(x=0,y^{\prime}=0)|^{2}+|\phi_{-}(x=0,y^{\prime}=0)|^{2}}$, as a function of $v_{y}$. The amplitude monotonously decreases with the growth of the velocity, vanishing at

[TABLE]

i.e., the mobile solitons exist in the limited interval of the velocities we .

The SV solitons may also be made mobile, but in a very narrow interval of velocities – e.g., at $v_{y}<\left(v_{y}\right)_{\max}^{(\mathrm{SV})}\approx 0.03$ for $\eta=0,\lambda=1$, and $N=3.7$, cf. the limit velocity given by Eq. (82) for the MMs. At $v_{y}>0.03$, the imaginary-time solution of Eq. (81) with the SV input converges to stable MM solitons, instead of the SV we .

III.3 3D metastable semi-vortices and mixed modes

III.3.1 Analytical considerations

The creation of metastable 3D solitons in the model based on Eq. (71) can be predicted, starting from evaluation of scaling of different terms in the respective energy functional (72). Assuming that a localized state has characteristic size $L$ and norm $N$, an estimate for the amplitudes of the wave function is $A\sim\sqrt{N}L^{-3/2}$. Accordingly, the three terms in Eq. (72) scale with $L$ as

[TABLE]

with positive coefficients $c_{\mathrm{kin}}$, $c_{\mathrm{soc}}$, and $c_{\mathrm{int}}^{\mathrm{(self/cross)}}$. As shown in Fig. 20, Eq. (83) gives rise to a local minimum of $E_{\mathrm{tot}}(L)$ at finite $L$, provided that

[TABLE]

Although this minimum cannot represent the GS, which formally corresponds to $E_{\mathrm{tot}}\rightarrow-\infty$ at $L\rightarrow 0$ in the collapsed state, as is suggested by Fig. 20 too, the local minimum corresponds to a self-trapped state which should be stable against small perturbations. Condition (84) suggests that metastable 3D solitons may exist when the SOC term is present, while its strength $\lambda$ is not too large, $N$ and $\eta$ being not too large either HP .

Like the 2D system, its 3D counterpart admits the existence of solitons of the SV and MM types. Also similar to the 2D case, for $\eta<1$ the energy of the SV is lower than that for the MM, and vice versa for $\eta>1$. Another similarity to the 2D system is that there is no threshold (minimum norm) necessary for the existence of the metastable 3D solitons, while their existence is bounded by a certain maximum norm, $N<N_{\max}$ HP .

In the numerical form, stationary 3D solitons were produced by means of the imaginary-time method applied to Eq. (71), and their stability against small perturbations was verified in real-time simulations HP . Typical examples of density profiles of stable SV and MM solitons are displayed in Fig. 21. Naturally, the MM states exhibit a more sophisticated profile. Lastly, similar to what is shown above for the 2D system in Fig. 19, it was found that stable 3D solitons can be set in motion in a limited interval of velocities, cf. Eq. (82).

The use of two-component BEC is the simplest possibility to emulated BEC in bosonic gases. More complex realizations were elaborated in models of three-component (spin-1) condensates. In particular, vortex solitons supported by SOC in spin-1 condensates have been numerically constructed too SKA-SOCvort1 ; SKA-SOCvort2 .

IV Giant microwave-coupled vortex rings in binary BEC

IV.1 The model

Another possibility to create stable 2D vortex solitons with arbitrarily large values of the vorticity was recently predicted in the hybrid model based on a system of two GPEs coupled by a microwave field through a magnetic-dipole transition. Two components of the matter waves, which are governed by the coupled GPEs, represent a BEC mixture of two different hyperfine states of the same atomic species, the microwave field being generated by the transition between the levels, which gives rise to the feedback of the matter waves on the microwave radiation. This physical model was introduced, in the 1D form, in Ref. Jieli1 .

The 2D system of the coupled GPEs for two components, $\phi_{\uparrow\downarrow}$, which form a pseudo-spinor wave function of the binary BEC, was derived in Ref. Jieli2 . In the scaled form, the system takes the form

[TABLE]

[TABLE]

Here, the integral terms represent the action of the microwave field on the matter waves, the field itself being replaced by the solution of the respective Poisson equation, produced by means of the 2D Green’s function, $\gamma$ is the strength of the field-matter interaction (which appears as the effective long-range interaction in Eqs. (85) and (86)), and $H_{0}$ is a background magnetic field. Further, $\beta$ is the coefficient of the additional contact (collisional) interaction between the components of the wave function, which may be present too, $\beta>0$ implying the attractive interactions.

For identical components, $\phi_{\downarrow}=\phi_{\uparrow}\equiv\phi\exp\left(-iH_{0}t\right)$, Eqs. (85) and (86) reduce to a single equation,

[TABLE]

where $\mathbf{r}^{\prime}=\left(x^{\prime},y^{\prime}\right)$, and its soliton solutions with vorticity $S$ and chemical potential $\mu$ can be looked for in the usual form,

[TABLE]

cf. Eq. (6). The Hamiltonian of Eq. (87), written in terms of the radial wave function $\Phi_{S}(r)$, defined as per Eq. (88), is

[TABLE]

The following analysis is presented for the fixed normalization of the wave function,

[TABLE]

(i.e., the total norm of the binary state is $1$).

IV.2 Basic results for vortex solitons

In the presence of the self-focusing contact interaction, $\beta>0$, solitons cannot be created by Eq. (87) (in other words, by Hamiltonian (89) beyond a point of the onset of the collapse, $\beta=\beta_{\mathrm{collapse}}$). A simple analytical approximation makes it possible to predict this point quite accurately for all vortex solitons, with $S\geq 1$ Jieli2 . To this end, one should make use of the fact that, according to numerical results, for $S\geq 2$ and $\beta$ large enough, vortex solitons generated by the present model assume the shape of narrow annuli, see an example for $S=5$ in Fig. 22(b). In the radial direction, it may be approximated by the usual 1D soliton ansatz. With regard to normalization (90), it is written as

[TABLE]

where $R$ is the radius of the narrow annulus. The substitution of this approximation in Eq. (89) yields

[TABLE]

Then, the equilibrium value of $R$ is selected as a point of the energy minimum: $dE/dR=0$, i.e.,

[TABLE]

and the analytically predicted collapse point, $\beta_{\mathrm{collapse}}^{\mathrm{(analyt)}}$, is one at which $R_{\mathrm{eq}}^{2}$ vanishes, i.e., the annulus collapses onto the center:

[TABLE]

As seen in Table II, this analytical prediction is very close to its numerically found counterparts at $S\geq 2$, and is rather close even at $S=1$, when the form of the vortex soliton is not actually close to a narrow ring, see Fig. 22(a).

It is worthy to note that prediction (94) does not depend on strength $\gamma$ of the effective long-range interaction, because this coefficient appears only as an overall factor in expression (93), and it remains valid in the limit of $\gamma\rightarrow 0$, i.e., for the usual two-dimensional NLSE with the local cubic self-attractive term,

[TABLE]

To explain this fact, we note that, at the eventual stage of the collapse, when the shrinking 2D vortex ring becomes extremely narrow, Eq. (87) becomes asymptotically tantamount to the local NLSE, therefore the condition for the onset of the collapse is identical in both equations (87) and (95). However, the latter equation has fundamental and vortex-soliton solutions solely at $\beta=\beta_{\max}$, which are completely unstable, while the effective long-range interaction in Eq. (87) creates stable fundamental solitons and vortices for all $S$ (see below). Actually, the approximate analytical result given by Eq. (94) provides an explanation for the numerical findings that were reported in Ref. Minsk2 and many other works. In this connection, it also relevant to mention that, taking into regard normalization condition (90), value $\beta_{\mathrm{collapse}}(S=1)=48.3$ in Table II is tantamount to the collapse threshold for $S=1$ written above in the form of Eq. (48).

The stability of the solitons was identified by real-time simulations of Eqs. (85) and (86) with random perturbations added to the initial conditions (independent perturbations were taken for $\phi_{\uparrow}$ and $\phi_{\uparrow}$, to verify the stability against breaking the equality of the two components) Jieli2 . The fundamental solitons ($S=0$) are stable in their entire existence region, $\beta<\beta_{\mathrm{collapse}}(S=0)\approx 11.8$ (in the present notation, it is tantamount to the TS norm, see Eq. (47)).

For the vortex solitons with $S\geq 1$, the numerical analysis reveals an internal stability boundary, $\beta_{\mathrm{stability}}(S)<\beta_{\mathrm{collapse}}(S)$, the vortices being stable at $\beta<\beta_{\mathrm{stability}}(S)$, see Table 1. In the interval of $\beta_{\mathrm{stability}}(S)<\beta<\beta_{\mathrm{collapse}}(S)$, they are broken by azimuthal perturbations into rotating necklace-shaped clusters of fragments, which resembles the initial stage of the instability development of localized vortices in usual models; however, unlike those models necklace1 ; necklace2 ; necklace3 , the necklace does not expand, remaining confined under the action of the effective nonlocal interaction. Typical examples of the stable and unstable evolution of the vortex solitons in the present model are displayed in Figs. 22.

The stability boundary, $\beta_{\mathrm{stability}}(S)$, can be found in an approximate analytical form too Jieli2 . To this end, the wave function of an azimuthally perturbed vortex ring is approximated by

[TABLE]

cf. Eq. (88). An evolution equation for the modulation amplitude, $A$, is derived by substituting ansatz (96) in Eq. (87) and averaging it in the radial direction:

[TABLE]

(recall $R$ is considered as the constant large radius of the narrow vortex annulus). In this approximation, the stability analysis amounts to studying the modulational stability of the solution with $A=1$ against azimuthal perturbations $\sim\exp\left(ip\theta\right)$, with its own integer winding numbers $p$, in the framework of the effectively one-dimensional NLSE (97). A simple result is that the stability is maintained under condition $p^{2}\geq\left(8/3\right)\left(\beta/8\pi\right)^{2}$, if the term $\sim\beta^{2}$ dominates in Eq. (97). Further, numerical results demonstrate that, as usual, the critical instability corresponds to $p^{2}=S^{2}$ (for instance, the appearance of five fragments in Fig. 22(d) demonstrates that, for $S=5$, the dominant splitting mode has $p=5$). Thus, it is expected that the vortex soliton remains stable at

[TABLE]

On the other hand, the numerically found stability boundary, presented in Table II, follows an empirical formula, $\beta_{\mathrm{stability}}(S)\approx\allowbreak 15S-4$. Thus, the analytical approximation given by Eq. (98) is accurate enough for $S\geq 2$.

It follows from these results that the giant vortex rings, with higher values of $S$, are much more robust than their counterparts with smaller $S$. This feature is opposite to what was previously found in those (few) models which are able to produce stable vortex solitons with $S>1$, cf. Refs. Pego ; Brtka ; 3DLHY ; GZ-ln . It is also relevant to mention that, while the fundamental soliton is the system’s GS at $\beta<\beta_{\mathrm{stability}}(S=0)$, the ground state does not exist at $\beta>\beta_{\mathrm{stability}}(S=0)$, due to the possibility of the collapse. Thus, the vortex solitons with winding number $S$ cannot represent the GS in the case of $\beta_{\mathrm{stability}}(S=0)<\beta<\beta_{\mathrm{stability}}(S)$ (e.g., $11.8<\beta<70$ for $S=5$, as per Table II). Nevertheless, the vortex solitons exist in this region as metastable states, cf. the 3D solitons in the SOC system based on Eq. (71), as illustrated above by Fig. 20. In that system, the vortical solitons of the SV and MM types also exist as metastable states, although the system does not have a GS, due to the possibility of the supercritical collapse HP .

Lastly, in the case of the self-repulsive contact interactions, $\beta<0$, Eqs. (85) and (86) give rise to stable solitons for all values of $S$ (at least, up to $5$) and indefinitely large values of $|\beta|$, as shown in work Jieli2 .

IV.3 A brief outline: vortex solitons in other systems with nonlocal

nonlinearity

The nonlocality is an essential feature of Eqs. (85) and (86). Other forms of long-range interactions naturally occur in many other physically significant models. In optics, the nonlocality naturally emerges in the propagation equation for beams in liquid crystals LC1 ; LC2 , as well as in the case of the thermal nonlinearity thermal . A typical model is based on the following two-dimensional NLSE Briedis :

[TABLE]

cf. Eq. (87), where $\sigma$ is a spatial scale of the long-range interaction.

It was demonstrated that, unlike the self-focusing local cubic (Kerr) nonlinearity, its nonlocal version may support stable vortex solitons, with unitary ($S=1$) and higher-order ($S\geq 2$) vorticities, as well as higher-order states characterized by zeros in the radial direction Briedis . This fact is not surprising, because, in the limit of strong nonlocality, i.e., with $\sigma$ much larger than a characteristic size of solitons, Eq. (99) reduces to an effectively linear equation for $u\left(x,y,z\right)\equiv\tilde{u}\left(x,y,z\right)\exp\left(iN/\left(\pi\sigma^{2}\right)\right)$, with the two-dimensional trapping HO potential:

[TABLE]

where $N$ is again the norm of the field (this linear equation, with the HO potential proportional to $N$, is known as the model of “accessible solitons” Snyder ). Obviously, the linear Schrödinger equation (100) admits a full set of stable states with all values of $S$, as well as other higher-order modes.

Experimentally, the creation of a stable two-component bound state, built of a vortex soliton in one optical beam, and a fundamental soliton in another one, carried by a different wavelength, in a nematic liquid crystal was demonstrated in Ref. Izdeb .

BEC made of dipolar atoms are described by GPEs with nonlocal anisotropic terms which represent dipole-dipole interactions Pfau . Stable 2D and 3D vortex solitons were predicted in various models including such interactions Lashkin ; Tikhon ; SKA-dipolar ; SKA-dipolar2 ; Goran , although no experimental results on this topic have been reported thus far.

V An outline of other topics and experimental studies

V.1 Conservative systems

There are other aspects of the studies of vortex solitons which are not included in the present review. As concerns the theory, a specific possibility is provided by the use of spatially nonuniform self-repulsive nonlinearity, with the local strength growing from the center to periphery at any rate faster than $r^{D}$, where $r$ is the distance from the center, and $D$ the spatial dimension. HS . This model makes it possible to create very robust solitons of diverse types, including fundamental ones and solitary vortices ICFO , further -further7 , as well as sophisticated 3D states – notably, hopfions (vortex tori with inner twist, which feature two independent topological numbers) hopfion . A challenging problem is finding still more general physically relevant conditions for the creation of complex 3D modes, such as skyrmions (which, similar to hopfions, carry two different topological charges) skyrmion1 ; skyrmion2 ; skyrmion3 , monopoles field-theory , linked vortex rings, globally coupled vortex clusters clusters , and others.

Another relevant ramification is the investigation of discrete vortex solitons in lattice media discrvort ; Pelin ; big ; Thaw , which essentially generalize the above-mentioned simplest discrete NLSE (56). Various forms of stable topological modes were created in virtual photonic lattices induced in strongly birefringent photorefractive materials (using the experimental technique proposed in Ref. MSegev ), see original works discrvort1 ; discrvort2 ; discrvort3 ; Jena-vortex and a review in big .

V.2 Dissipative systems

V.2.1 Free-space models

A separate vast research area is the study of multidimensional solitons with embedded vorticity in dissipative media, see, e.g., a review in Ref. Dum-diss and recent theoretical works NNR1 ; NNR2 ; NNR3 ; DimaHS ; DimaThaw2 ; DimaThaw1 and references therein. A basic class of models for the generation of dissipative vortex solitons is based on CGLEs Kramer with the CQ or, more generally, saturable nonlinearity. In particular, a model for the evolution of amplitude $u\left(x,y,z\right)$ of electromagnetic fields in laser cavities Lugiato ; NNR0 ; Mandel may be based on the following form of CGLE, which was first introduced, on phenomenological grounds, by Petviashvili and Sergeev Petvia :

[TABLE]

where positive coefficients $\gamma$, $\mu$, $\varepsilon$, and $\nu$ account, severally, for the linear and quintic losses, cubic gain, and coefficient of the quintic self-defocusing. The presence of the linear loss is necessary to make the zero background stable, thus admitting the existence of stable dissipative solitons. The gain is provided by the cubic term, while the quintic loss is necessary for the overall stability of the model. The presence of the effective viscosity term with $\beta>0$ is necessary for stability of vortex solitons Lucian . This term may represent, e.g., diffusion of free carriers (electrons and holes) generated by the electromagnetic fields in the laser cavity.

Solutions to Eq. (101) for vortex dissipative solitons with winding number $S$ are looked for, in terms of the polar coordinates $\left(r,\theta\right)$, as

[TABLE]

where real $k$ is a propagation constant, $\Phi(r)$ is a real phase, and $A(r)$ is a real amplitude function. Unlike conservative models, such as the one based on the NLSE with the CQ nonlinearity, see Eq. (10), in which solitons exist in the form of continuous families parametrized by $k$, CGLEs give rise only to pairs of dissipative solitons with two discrete values of eigenvalues $k$ for each integer $S$, starting from $S=0$. One of these solutions, which features a larger amplitude, may be stable, as an attractor in the dissipative system, while the additional solution plays the role of a separatrix between attraction basins of the attractors represented by the zero solution and stable dissipative soliton PhysicaDD .

Unlike solitary vortices in conservative media, the phase structure of their dissipative counterparts assumes the spiral shape Lucian . Indeed, the asymptotic form of solution (102) at $r\rightarrow\infty$ is

[TABLE]

with a real radial wavenumber, $q$, and $\lambda>0$ accounting for the exponential localization of the dissipative soliton. In particular, $q=\gamma/\lambda$ in the case of $\beta=0$. The spiral structure of the phase is produced by the combination of terms $S\theta$ and $qr$ in Eq. (103), see an example for $S=1$ in Fig. 23. The model based on CGLE (101) with $\beta>0$ gives rise to stable spiral solitons with vorticities $S=1$ and $S=2$ Lucian .

The 3D version of the CGLE with the CQ nonlinearity also makes it possible to generate stable 3D vortex solitons, in the form of “donuts”, with embedded vorticities $S=1$ and $2$ CGL-3Dvort , provided that the quasi-2D viscosity term, with $\beta>0$, is kept in the 3D version of Eq. (101). An example of stable evolution of an initially perturbed 3D vortex solitons with $S=2$ is displayed in Fig. 24.

V.2.2 Dissipative systems with the trapping potential

In most cases, models of laser cavities do not contain the diffusion term $\sim\beta$ in Eq. (101), because photons are not subject to diffusion. As mentioned above, dissipative vortex solitons (unlike ones with $S=0$) cannot be stable, in the free space, in the absence of this term. However, they can be readily stabilized by trapping potentials Skarka0 . In particular, the respectively modified CGLE,

[TABLE]

where $\Omega^{2}$ is the strength of the HO trapping potential, supports stable vortex modes with $S=1$, while all the states with $S\geq 2$ are unstable against spontaneous splitting, similar to what is demonstrated above for the 2D version of the conservative NLSE/GPE (29) with the same HO trapping potential. Also similar to the conservative model, those vortex modes with $S=1$ which are unstable suffer splitting in the framework of Eq. (104). However, in contrast with the case of Eq. (29), in the latter case fragments produced by the splitting do not perform cycles of periodic recombination-splitting evolution, see Fig. 8, nor are they destroyed by the intrinsic collapse. Instead, they form a stably rotating dipole, as shown in Fig. 25. Similarly, an unstable dissipative vortex soliton with $S=2$ splits in a set of three fragments, which form a steadily rotating “tripole” Skarka0 .

Two-dimensional vortex patterns in dissipative media may be also supported by the linear gain applied in an appropriately defined spatially confined area, such as a ring-shaped one ring (experimentally, a similar structure was created for pumping a spin-orbit-coupled exciton-polariton condensate in a semiconductor microcavity Amo ). In a similar context, many vortex modes have been found in the framework of the two-dimensional CGLE with the CQ nonlinearity, taken in the form similar to Eq. (101), and with a radially-dependent linear-loss coefficient, which has a minimum at $r=0$ Skarka1 :

[TABLE]

Even without the inclusion of the diffusion term, $\beta=0$ (whose physical origin is not straightforward), this model gives rise to a great variety of stable vortex states, such as ordinary stationary vortex solitons, elliptically deformed and crescent-shaped steadily rotating ones, and eccentrically deformed vortices, which combine steady rotation and orbiting around the origin, see an example in Fig. 26.

A great variety of dynamical self-trapped patters, initiated by inputs with embedded vorticity, is produced by CGLE (101) with linear gain, instead of the linear loss, placed around $r=0$, which corresponds to

[TABLE]

with $\gamma_{0}>0$ (cf. Eq. (105)), i.e., the gain area is introduced at $0\leq r^{2}<\gamma_{0}/\gamma_{2}$. As reported in Ref. Skarka2 , this model readily gives rise to spontaneous breaking of the underlying axial symmetry and thus generates rotating and periodically transmuting patterns shaped, roughly, as four-, five-, six-, seven-, eight-, nine-, and ten-pointed stars. While these rotating modes carry the orbital angular momentum, they lose the original vorticity (the phase singularity is expelled from the original pattern).

A recent ramification of studies of vortex solitons under the action of the trapping potential is represented by the model with the uniform linear gain, $\gamma>0$ (and still including the diffusion term), cubic loss ($\varepsilon>0$), and HO trap DimaThaw2 . The model may apply to laser cavities and exciton-polariton condensates:

[TABLE]

(here the diffusion coefficient is denoted $\eta$ instead of $\beta$, cf. Eq. (104)). The decrease of $\eta$ in this model gives rise to a variety of stably rotating multi-vortex states, with $S=1,2,3,4,6,7$, as shown in Fig. 27. For the same values of other parameters and $\eta<$ $0.22$, the model generates a turbulent state, while at $\eta>0.5$ it gives rise to stable axisymmetric vortex with $S=1$.

A similar model was elaborated for the model of the spinor (two-component) exciton-polariton condensate under the action of the SOC DimaHS ,

[TABLE]

where $\beta$ is the SOC strength (note that in exciton-polariton condensates SOC is represented by the second-order differential operators, unlike the first-order ones in the model of the atomic BEC, cf. Eq. (62)), $\alpha$ is the coefficient of the nonlinear interaction between the components, the strength of the HO trapping potential is fixed to be $1$ by rescaling, and $\Omega$ represents the Zeeman splitting between the components (if any).

The analysis of Eqs. (LABEL:HSDima) has demonstrated that, in the absence of the Zeeman terms ($\Omega=0$), the system may give rise to stable vortex-antivortex (VAV) complexes, in the form of

[TABLE]

written in polar coordinates $\left(r,\theta\right)$, with a common amplitude function $\phi(r)$, and stable MM (mixed-mode) states, see Fig. 28(a). The figure demonstrates a small bistability region in which the VAV and MM species coexist as stable ones. In addition to that, Fig. 28(b) shows that the system with $\Omega>0$ may support stable SV (semi-vortex) states, in the form of

[TABLE]

In comparison with the SV defined as per Eqs. (73) and (77), the second-order SOC operator in Eqs. (LABEL:HSDima) imparts winding number $S=2$ to the vortex component of the SV, instead of $S=1$.

V.3 $\mathcal{PT}$-symmetric systems

A specific variety of the dissipative systems is represented by ones subject to the condition of the parity-time ($\mathcal{PT}$) symmetry, i.e., settings with spatially separated and mutually balanced gain and loss elements. They maintain a real spectrum of energy eigenvalues (i.e., avoid instability, which is often driven in dissipative systems by linear gain), provided that the strength of the gain-loss terms does not exceed a certain critical level PT1 ; PT2 . In the combination with nonlinearity, $\mathcal{PT}$-symmetric systems readily create solitons PTsol1 ; PTsol2 . While these solitons are usually considered in the 1D form, they may be found in multidimensional models as well. In particular, in some specific systems, such as one with a $\mathcal{PT}$-symmetric lattice potential Optica , and another one with a specially designed profiles of spatially inhomogeneous nonlinearity and gain-loss terms Luz , stable vortex solitons were predicted too. Stable “dark” vortices in a $\mathcal{PT}$-symmetric medium, supported by a finite-amplitude background field, were studied too Panos-PTvort .

V.4 A brief outline of experimental results

In terms of experiments, the field of multidimensional solitons remains a challenging one. As concerns vortex solitons, which is the central topic of the present review, none of them was experimentally created in a stationary form. Transient modes in the form of semi-discrete “bullets” with embedded vorticity were observed in arrays of optical waveguides with the Kerr nonlinearity, which form a hexagonal lattice in the transverse cross section Jena-vortex , see an illustration in Fig. 29. Spatial 2D solitons with embedded vorticity, which were effectively stabilized, as transient states, by nonlinear (cubic) losses in the medium with the CQ nonlinearity (liquid carbon disulfide), were reported too Cid3 , as shown in Fig. 30. Thus, the current state of the topic of vortex solitons calls for performing new experiments.

VI Conclusion: A summary of the review

This article aims to present a relatively brief review of the broad area of multidimensional solitons, focusing on ones with embedded vorticity. The multidimensional solitons and solitary vortices find most important physical realizations in BEC and nonlinear optics, the crucially important problem being search for physically relevant settings which provide for the stabilization of these states against the collapse and splitting, that tend to destroy the fundamental (zero-vorticity) and vortical solitons, respectively, in media with the most common cubic self-attractive nonlinearity.

The article reviews, first, some well-established topics, which remain relevant to the current studies. One topic is the stabilization of the 3D and 2D states with vorticities $S=0$ and $S\geq 1$ by competing attractive and repulsive nonlinearities. This part includes the newest addition, in the form of the 3D and 2D GPEs with the LHY (Lee-Huang-Yang) corrections, which is related to the recent breakthrough in theoretical and experimental studies of self-trapped states in atomic BEC, viz., QDs (quantum droplets). The other well-established topic deals with the use of trapping HO (harmonic-oscillator) and spatially periodic (lattice) potentials, combined with the cubic self-attraction, for the creation of stable vortex solitons featuring the multi-peak structure.

Two other sections of the article address two recently elaborated schemes for the stabilization of multidimensional vortex solitons in BEC with cubic nonlinearities. One scheme predicts the creation of stable 2D and 3D solitons, which mix terms with $S=0$ and $1$ (SVs, i.e., semi-vortices) or $S=0$ and $\pm 1$ (MMs, mixed modes), in a two-component system which realizes the SOC in the atomic condensates. The second novel scheme predicts giant vortex rings (which may be stable with indefinitely large values of $S$) in 2D two-component BEC coupled by microwave radiation. A remarkable peculiarity of the latter model is that virtually all the essential results, concerning the existence and stability of the vortex solitons in it, can be produced in an approximate analytical form.

The analysis of the two novel schemes reveals a drastic difference between the 2D and 3D settings. In the former case, the stabilization mechanism creates the GS (ground state), that did not exist otherwise. The total norm of the GSs takes values below the threshold necessary for the onset of the critical collapse, driven by the cubic attractive nonlinearity in the 2D space. In the 3D settings, the supercritical collapse does not allow the creation of a GS, but, nevertheless, an appropriate mechanism may create metastable solitons.

The review also summarizes, in a brief form, some other results relevant to the general topic of vortex solitons. In particular, these are relatively old and new findings for self-trapped vortex states in dissipative media. Relevant experimental findings are briefly outlined two.

To conclude, it is relevant to stress that there remain vast room for further theoretical and, especially, experimental studies of vortex solitons and objects related to them in photonics, BEC, and other quickly developing areas of physics.

Acknowledgments

I would like to thank my collaborators in original works that were dealing with various topics related to multidimensional solitons (most of those works are cited in this review): C. B. de Araújo, B. B. Baizakov, V. Besse, O. V. Borovkova, G. Boudebs, C. M. Brtka, W. B. Cardoso, R. Carretero-González, Zhaopin Chen, Zhigang Chen, J. Cuevas, G. Dong, N. Dror, Z. Fan, D. J. Frantzeskakis, A. Gammal, G. Gligorić, L. Hadžievski, X. Jiang, Y. V. Kartashov, P. G. Kevrekidis, V. V. Konotop, H. Leblond, B. Li, Y. Li, V. E. Lobanov, A. Maluckov, D. Mazilu, D. Mihalache, W. Pang, H. Pu, J. Qin, A. S. Reyna, H. Sakaguchi, L. Salasnich, M. Salerno, E. Ya. Sherman, Ya. Shnir, L. Tarruell, F. Toigo, L. Torner, F. Wise, Y.-C. Zhang, and Z.-W. Zhou. I also thank Prof. Victor Pérez-García, the editor of Physica D, for the invitation to write a review article for the journal.

Parts of this work were supported by grants No. 2010239 from the Binational (US-Israel) Science Foundation (BSF), and No. 2015616 from the joint program in physics between the National Science Foundation (US) and BSF.