# Multidimensional solitons: Well-established results and novel findings

###### Abstract

A brief review is given of some well-known and some very recent results obtained in studies of two- and three-dimensional (2D and 3D) solitons. Both zero-vorticity (fundamental) solitons and ones carrying vorticity are considered. Physical realizations of multidimensional solitons in atomic Bose-Einstein condensates (BECs) and nonlinear optics are briefly discussed too. Unlike 1D solitons, which are typically stable, 2D and 3D ones are vulnerable to the instability induced by the occurrence of the critical and supercritical collapse, respectively, in the same 2D and 3D models that give rise to the solitons. 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 one well-established topic, viz., the stabilization of the 3D and 2D states, with and , trapped in harmonic-oscillator (HO) potentials, and another topic which was developed very recently: the stabilization of 2D and 3D free-space solitons, which mix components with and (semi-vortices and mixed modes), in a binary system with the spin-orbit coupling (SOC) between its components. In both cases, the generic situations are drastically different in the 2D and 3D geometries. In the former case, the stabilization mechanism creates a stable ground state (GS, which was absent without it), whose norm falls below the threshold value at which the critical collapse sets in. In the 3D geometry, the supercritical collapse does not allow to create a GS, but metastable solitons can be constructed.

List of acronyms: GPE – Gross-Pitaevskii equation; GS – ground state; HO – harmonic-oscillator (potential); MM – mixed mode; NLSE – nonlinear Schrödinger equation; SOC – spin-orbit coupling; SV – semi-vortex; VA – variational approximation; VK – Vakhitov-Kolokolov (stability criterion)

## I Introduction

Self-trapping of multidimensional (two- and three-dimensional, 2D and 3D) localized modes in nonlinear dispersive/diffractive media, which are usually categorized as solitons, has been a topic of great interest in many areas of physics, as can be seen from reviews general-reviews ()-RMP (). In particular, this area of the theoretical and experimental studies is highly relevant to nonlinear optics, as well as to nonlinear dynamics of matter waves in Bose-Einstein condensates (BECs). Other essential realizations of multidimensional solitons are known in models of ferromagnets ferro (), superconductors super (), and semiconductors semi (), as well as in nuclear physics nuclear (), the classical-field theory hopfion (); field-theory () (including famous skyrmions skyrmion (), which may be considered as 3D solitons), and other fields.

In addition to the their significance to fundamental studies, 2D and 3D solitons offer potential applications. In particular, spatiotemporal optical solitons, alias “light bullets” Silb (), may be used as data carriers in all-optical data-processing schemes Wagner (), and 3D matter-wave solitons are expected to provide a basis for highly precise interferometry interferometry ().

Unlike 1D solitons, which are usually stable objects KA (); Michel (), stability is a major issue for their 2D and 3D counterparts. Indeed, the most common cubic self-focusing nonlinearity, which can easily support 2D and 3D formal solitons solutions, gives rise to the wave collapse (critical collapse in 2D and supercritical collapse in 3D) collapse (), which completely destabilizes the respective soliton families. In particular, the first ever example of solitons which was introduced theoretically in nonlinear optics, viz., the Townes solitons Townes (), i.e., 2D self-trapped modes supported by the cubic self-focusing nonlinearity (in fact, these solutions were discovered before term “soliton” had been coined) are subject to the instability caused by the critical collapse, therefore they have never been observed in the experiment. Multidimensional solitons with embedded vorticity, alias 2D vortex rings and 3D vortex tori, which also exist in cubic media (in particular, there are higher-order versions of Townes solitons with embedded vorticity Minsk ()), are vulnerable to a still stronger splitting instability initiated by modulational perturbations in the azimuthal direction; usually, the modulational instability splits the vortex ring into a few fragments, which resemble fundamental (zero-vorticity) solitons and suffer the destruction via the intrinsic collapse or decay into small-amplitude waves (“radiation”) general-reviews ().

Thus, the stabilization of multidimensional fundamental and vortex solitons in physically relevant models is an issue of great significance to theoretical studies, and creation of the solitons in the respective physical settings is a challenge to the experiment. This paper offers a short review of the broad area of theoretical investigations of 2D and 3D solitons, selecting two topics for a relatively detailed presentation, one well-established, and one very recent. The former one, considered in Section II, relies upon the stabilization of both fundamental (zero-vorticity) solitons and ones with embedded vorticity by an external harmonic-oscillator (HO) trapping potential, axially or spherically symmetric one. This topic has been elaborated in detail in the course of the past 20 years 2D (); Sadhan (); Ueda (); Dum2D (); Dum3D (). A recently developed theme, presented in Section III, is the creation of stable solitons in the 2D and 3D free space (without the help from any external potential), which mix zero-vorticity and vortical components, by means of linear terms that account for the spin-orbit coupling (SOC) in binary BEC. This possibility for 2D and 3D solitons was discovered in Refs. we () and HP (), respectively.

Any stabilization mechanism must secure the solitons against the critical collapse in 2D settings, and against the supercritical collapse in the 3D geometry. This objective is achieved differently in 2D and 3D models, as explained below. In the former case, the stabilizing factor breaks the specific scale invariance, which underlies the critical collapse, and transforms any originally unstable fundamental soliton in the single-component model, or a mixed (fundamental-vortical) soliton in the binary SOC system, into a ground state (GS), which does not exists in the unstabilized 2D model. This is possible because the critical collapse sets in when the norm of the mode exceeds a certain threshold value, and the solitons get stabilized by pushing their norm below the collapse threshold. Vortex solitons with in the 2D single-component model may be stabilized too, as excited states, rather than the GS, provided that their norms are not too large either. On the other hand, in the 3D setting the threshold for the onset of the supercritical collapse is absent, hence a true GS cannot be created. Nevertheless, the stabilization of 3D fundamental solitons and vortical ones with in the single-component model, and of 3D mixed modes in the SOC system, is possible in the form of metastable states, which are robust against small perturbations.

## Ii Multidimensional solitons in trapping potentials

### ii.1 The basic model

A generic method for the stabilization of both fundamental and vortical solitons is the use of trapping potentials, the corresponding model being based on the following nonlinear Schrödinger equation (NLSE), alias the Gross-Pitaevskii equation (GPE) Pit (), for the atomic mean-field wave function , where (or in the 2D geometry) and are appropriately scaled spatial coordinates and time:

(1) |

Here, is the 3D or 2D Laplacian, is the trapping potential, and the negative sign in front of the cubic term implies that the nonlinearity is self-attractive. The trapping potential in the 3D setting may be quasi-two-dimensional, i.e., low-dim ().

In the application to optics, the evolution variable in Eq. (1) is replaced by the propagation distance (), while the longitudinal coordinate is replaced by the local time, , where is the group velocity of the electromagnetic carrier waves KA (). In this case, the trapping potential, which defines a waveguiding channel in the bulk medium by means of the local modulation of the refractive index, , may be only quasi-two-dimensional, as it is virtually important to impose modulation of the refractive index providing the trapping in the temporal direction. The 2D version of Eq. (1), with the propagation distance replacing , models the evolution of spatial light beams in the bulk medium with transverse coordinates .

Axially symmetric potentials (including spherically isotropic ones), , where is the set of cylindrical coordinates, admit looking for solutions to Eq. (1) in the form of vortex modes,

(2) |

with real chemical potential and integer vorticity . Real amplitude function satisfies the stationary equation,

(3) |

supplemented by the boundary conditions (b.c.), at and , and at (assuming ).

Relevant to the experiments with BEC is the HO trapping potential Pit (),

(4) |

where accounts for the anisotropy of the trap, determining the aspect ratio, , between trapping lengths in plane and along axis . Accordingly, the cases of and corresponding, respectively, to nearly two-dimensional (“pancake-shaped”, see Ref. pancake () and references therein) and nearly one-dimensional (“cigar-shaped” Randy-NJP ()) condensates. Equation (1) conserves the energy,

(5) |

norm (scaled number of atoms in the BEC), or the total energy in optics),

(6) |

and -component of the angular momentum,

(7) |

with standing for the complex conjugate (in the case of the isotropic trap, , all the three components of the angular momentum are conserved). For the stationary state (2) with a definite value of , relation holds.

In this section, the basic results for the structure and stability of trapped solitons and solitary vortices are summarized for the HO potential (4), following, chiefly, Refs. Dum2D () and Dum3D (). On the other hand, spatially periodic (lattice) potentials also offer an efficient tool for the stabilization of 2D and 3D fundamental and vortex solitons, as was predicted in various settings lattice (); low-dim (), and demonstrated experimentally for 2D optical solitons with embedded vorticity in photonic lattices photorefr (), and for 2D plasmon-polariton solitons in microcavities with a lattice structure 30 ().

### ii.2 Stationary 3D modes: analytical and numerical results

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

(8) |

where , and are stationary wave functions of 1D harmonic oscillators with quantum numbers , which correspond to energy eigenvalues , and , respectively, the full chemical potential being

(9) |

The states which go over into ones (2) with vorticity in the nonlinear model are constructed as combinations of factorized wave functions (8) with and . Since the correction to from the self-attractive nonlinearity is negative, this restriction and Eq. (9) with impose a bound on ,

(10) |

In particular, the eigenfunctions of the linear model, with vorticities and , are

(11) |

(12) |

To apply the variational approximation (VA) Progress () to the nonlinear model, note that Eq. (3) can be derived from Lagrangian

(13) |

A natural ansatz for the trapped vortices is

(14) |

where amplitude , as well as the radial and vertical widths, and , are free parameters, norm (6) of this ansatz being (below, replaces as one of variational parameters). The trapped modes with are shaped as vortex tori (“doughnuts”), see a typical example below in Fig. 4. This shape is (qualitatively) correctly captured by ansatz (14).

Further results produced by the VA are given here for and , as all vortices with are unstable, see below. The substitution of ansatz (14) with in Lagrangian (13) yields . The first variational equation, , yields , which predicts that the vortex may only exist if it is narrow enough in the radial direction, . Two other equations, , yield , .

Numerical solution of the variational equations generates dependences for families of the fundamental () and vortical () trapped modes which are displayed, severally, in Figs. 1(a) and 2(a) for , , and , which correspond to the pancake-like, isotropic, and cigar-shaped configurations, respectively, while the corresponding curves are displayed in Figs 1(b) and 2(b). The same dependences, obtained from numerical solution of stationary equation (3), are displayed too, demonstrating the accuracy of the VA. It is seen that the properties of the modes trapped in the isotropic () and cigar-shaped () HO potentials are quite close, being strongly different from those is the pancake trap (). Another obvious feature of the figure is the presence of a largest norm, which corresponds to the turning point of the curve and bounds the existence of the trapped modes. This limitation is caused by the presence of the supercritical collapse in the 3D model collapse (): if the norm is too large, the trend to the collapse, driven by the cubic self-attraction, cannot be balanced by the quantum pressure (the gradient part of energy (5)). The largest norm in this model was numerically found, as a function , 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.

### ii.3 Stability of trapped 3D fundamental and vortex modes

The stability of stationary solutions produced by Eq. (3) is identified through the computation of (generally, complex) eigenvalues of small perturbations. To this end, a perturbed solution to Eq. (1) is looked for as

(15) |

where are eigenmodes of infinitesimal perturbations corresponding to integer values of azimuthal index . The substitution of ansatz (15) in Eq. (1) and linearization lead to the Bogoliubov - de Gennes equations Pit (),

(16) |

supplemented by the boundary conditions demanding that and decay exponentially at and , and decay as at . The stability of the underlying stationary mode is secured by the condition that all eigenvalues produced by numerical solution of Eq. (16) Ueda () must be pure imaginary. The findings are included in Figs. 1 and 2, where stable and unstable portions of the solution families are distinguished. The results of the stability analysis are collected in the diagram which displays stability regions in the plane of in Fig. 3.

The stability of the fundamental trapped states () precisely complies with the Vakhitov-Kolokolov (VK) criterion VK (); collapse (), , which is a necessary but, in the general case, not sufficient stability condition. On the other hand, the stability region of the vortex mode with is essentially smaller than formally predicted by this criterion. In the region with the VK criterion still holds for the vortices, but they are actually unstable, they are vulnerable to the above-mentioned splitting instability induced by perturbations with and in Eq. (15). The vortex modes with were found to be completely unstable.

A noteworthy feature revealed by Figs. 1 and 2 is that the smallest and largest maximum values of , up to which the trapped modes with or remain stable, are achieved, respectively, at the largest and smallest values of . In particular, the quasi-1D “tubular” vortical modes with Luca (), which are generated by the model at , are tightly confined in the corresponding cigar-shaped trap, that suppresses their destabilization by the splitting.

As said above, at the 3D setting shrinks into its 2D counterpart, where trapped modes with have a finite stability region, , see below. Therefore, taking into regard that Eq. (10) yields for , the lower stability border () for in Fig. 3 becomes, at large , asymptotically linear and parallel to the upper existence border, . On the other hand, the stability region for the fundamental trapped modes expands at , because the entire fundamental family is stable in the 2D limit, see below too.

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

Lastly, it is relevant to stress that stable matter-wave solitons, created experimentally in condensates of Li Randy-NJP () and Rb Rb85 () atoms loaded in cigar-shaped traps, may be well described by the above solutions for . On the other hand, stable solitons were also observed in the post-collapse state of the Rb condensate trapped in the HO potential with , which is much closer to the isotropic configuration.

### ii.4 2D fundamental and vortex trapped modes

As said above, taking the limit of in HO potential (4) leads to the shrinkage of the 3D model into its 2D simplification, with Eqs. (GPE), (3), (5), (6), (7), and (16) carrying over into their 2D counterparts, and the 2D HO potential taking the form of . In particular, stationary solutions of the 2D version of Eq. (1) are looked for as , cf. Eq. (2), where is realized as the chemical potential of the 2D model, i.e., its limit form corresponding to the linear equation is

(17) |

cf. Eq. (10).

The shape and stability of the 2D trapped modes was studied in a number of works 2D (); Dum2D (). The results, produced by the computation of the stability eigenvalues in the framework of the 2D version of Eq. (16), are summarized in Fig. 5, which demonstrates that the family of the fundamental () trapped modes is entirely stable in its existence region,

(18) |

where is the numerically found value of the norm of the Townes soliton, which determines the threshold for the onset of the critical collapse in the 2D NLSE collapse () (the variational method makes it possible to produce an analytical approximation, , with a relative error Anderson ()). In the free space (in the absence of any trapping potential), the family of the Townes solitons is degenerate, in the sense that they all have the single value of the norm, which is exactly equal to , independently of the soliton’s chemical potential, . The degeneracy is a consequence of the specific scaling invariance of the NLSE in 2D collapse (). The trapping potential introduces a characteristic spatial scale (the standard HO length, in dimensional units, or merely the spatial period, in the case of lattice potentials), which breaks the scaling invariance and thus lifts the degeneracy, making a function of . In fact, Eq. (18) 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 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 is possible at any value of the norm), and only metastability of the trapped modes is possible.

The family of Townes solitons with embedded vorticity is degenerate too in the free space, being pinned to the single value of the norm, Minsk (). As seen in Fig. 5, the HO trapping potential stabilizes them in interval

(19) |

the corresponding stability interval in terms of the chemical potential being (the right edge of the interval is determined by Eq. (17)). As for the vortices with , they remain completely unstable, as in the 3D model.

The stability of the trapped vortex modes with , predicted through the computation of the corresponding eigenvalues, was corroborated by direct simulations of the perturbed evolution, see a typical example in Fig. 6.

The direct simulations also reveal an interesting dynamical regime of the trapped vortices with in interval

(20) |

adjacent to one given by Eq. (19). In this interval, the evolution of the unstable vortex is regular (time-periodic), as shown in Fig. 7: it spits into two fragments which then recombine back into the vortex. The splitting-recombination cycles recur periodically, keeping the vorticity of the configuration. In this dynamical regime, the trapped vortices may be considered as semi-stable modes.

The vortices with still larger values of the norm, , also split into two fragments, which, however, fail to recombine. Instead, each one quickly blows up, i.e., collapses. Note that the corresponding effective collapse threshold, , is close to the double collapse threshold of the fundamental solitons, , which is consistent with the observation that the two fragments tend to suffer the intrinsic collapse if they do not recombine.

A noteworthy generalization of the above analysis was performed for a system of two nonlinearily coupled fields, which may be realized as a binary BEC, or as the copropagation of two optical beams in a bulk waveguide Brtka (). The respective coupled 2D GPE/NLSE system is

where is the relative strength of the attraction () or repulsion () between the components. This system gives rise to composite modes with hidden vorticity, formed by the components with equal norms and opposite vorticities, . Such modes are stable, roughly, in the parameter region with . |

(22) |

where , are the Pauli matrices, is a real coefficient of the SOC of the Rashba type (the Dresselhaus coupling, with combination , instead of in Eq. (22), is not considered here, as it tends to destroy 2D solitons, rather than to create them Sherman ()), is the relative strength of the cross-attraction between the components (cf. Eq. (LABEL:system)), under the condition that the strength of the self-attraction is normalized to be , and is the strength of the Zeeman splitting. Coefficient has the dimension of length, defining a fixed scale which breaks the above-mentioned scale invariance of the NLSE/GPE in the free 2D space. This effect makes it possible to create stable solitons with norms falling below the collapse threshold, see below.

Stationary solutions of Eq. (22) for 2D solitons with real chemical potential are looked for as , where complex stationary wave functions are determined by equations

(23) | |||

(24) |

Dynamical invariants of Eqs. (22) are the total norm, energy, and linear momentum:

(25) |

(26) |

(27) |

It is relevant to note that a solution to the linearized version of Eq. (22), in the form of , gives rise to two branches of the spectrum,

(28) |

which is gapless in the case of , and features a gap, , if the Zeeman splitting is included.

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. (22) Wesley (). This 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 the 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 (“light bullets”).

The 3D model of the SOC system will be taken as per Ref. HP ():

(29) |

where is again the SOC coefficient, and the 3D matrix vector is . The 3D system conserves the norm and momentum, given by 3D counterparts of Eqs. (25) and (27), and energy

(30) | |||

The Zeeman splitting is not included into the 3D system, but it may be added to it.

### iii.2 Stable 2D solitons: standing and mobile semi-vortices and mixed modes

#### iii.2.1 2D semi-vortices (at )

First, we note that Eq. (22) with (in the absence of the Zeeman splitting) admits stationary solutions with chemical potential , written in terms of the polar coordinates:

(31) |

where real functions take finite values and have zero derivatives at , and feature the following asymptotic form at :

(32) |

with constants and . As it follows from Eq. (32), the solutions are exponentially localized, as solitons, at

(33) |

Solutions (31) are built as bound states of a fundamental (zero-vorticity, ) soliton in component and a solitary vortex, with vorticity , in , therefore composite modes of this type are called semi-vortices (SVs) we (). The invariance of Eq. (1) with respect to transformation

(34) |

gives rise to a semi-vortex which is a mirror image of (31), with replaced by :

(35) |

Numerically, stable SVs were generated, as solutions to Eq. (1) by means of imaginary-time simulations im-time (), starting from the Gaussian input,

(36) |

where and are real constants. A typical example of the SV is displayed in Fig. 8(a).

Further, Fig. 8(b) represents the family of the SVs, showing their chemical potential as a function of the norm, cf. Fig. 5(a). Note that the dependence satisfies the VK criterion, , which is the above-mentioned necessary condition for the stability of solitary modes supported by the self-attractive nonlinearity. The family of the SV solitons exists precisely in the interval of norms (18), 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 necessary for the existence of the SVs in the free space. In the limit of , the vortex component of the SV vanishes, while the fundamental one degenerates into the usual Townes soliton, with , as shown by means of the dependence of ratio on in Fig. 8(d).

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

The wave form (36) was used not only as the input for the imaginary-time simulations, but also as an ansatz for the application of the VA. The substitution of the ansatz into expression (26) for the energy yields

(37) |

while norm (25) of the ansatz is Then, values of amplitudes and inverse squared widths , of the ansatz are predicted by the minimization of with respect to the variational parameters, , under the constraint that is kept constant. Figure 8(c) displays the comparison of the so predicted amplitude and the maximum value of obtained from the numerical solution at . It is seen that the accuracy of the VA is quite reasonable.

Very recently, 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 a new parameter of the family of SV solitons, viz., a spontaneous shift of the pivot of the vortical component () with respect to its zero-vorticity counterpart (). Those solitons also demonstrate an essentially different mobility scenario from the one outlined below for the present model: they respond to an applied kick by drift in the opposite direction (with an effective negative mass) along a spiral trajectory.

#### iii.2.2 2D mixed modes (at )

Aside from the SVs, the same SOC system (22) gives rise to another type of vorticity-carrying solitons, in the form of MMs, which combine terms with zero and unitary vorticities, and , in the spin-up and spin-down components, and . Numerically, the MM can be produced by imaginary-time simulations initiated by the following input, which may also serve as the variational ansatz:

(38) |

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

A typical example of the MM state is displayed in Fig. 9(a). Note that peak positions of the two components, and , in this state are separated along , Fig. 9(d) showing the separation () as a function of the norm. For a small amplitude of the vortex component, , Eq. (38) yields .

The dependence for the MM family is displayed in Fig. 9(b), which shows that the VK criterion holds in this case too, and, as well as SVs, the MMs do not have any threshold value of necessary for their existence. The family exists in the interval of , where is the critical norm (18) corresponding to the Townes solitons. In the limit of the vortex components vanish in the MM, and it degenerates into a two-component Townes soliton, cf. the degeneration of the SV in the limit of , which was shown in Figs. 8(b,d). The separation between peaks of the two components vanishes in this limit too, see Fig. 9(d).

The VA for the MM solutions is generated by the insertion of expression (38), as the variational ansatz, into energy functional (26), which yields

(39) |

the total norm of the ansatz being . Numerical solution of the respective variational equations, which minimize energy (39) under the constraint of the fixed norm, produces results which are compared in Fig. 9(c) with their numerically found counterparts. The VA for the MM solution is less accurate than for its VA counterpart, cf. Fig. 8(c), but it is usable too.

Direct simulations demonstrate that the MMs are unstable at , and stable at , i.e., precisely in the regions where the SVs are, severally, stable and unstable (exactly at , both the SV and MM solutions are stable Fukuoka2 ()). This stability switch between the SV and MM is explained by the comparison of energy (26) for them at equal values of the norm: the energy is smaller for the SV at , and for the MM at we () (similar to what is shown below in Fig. 13(a) for the 3D system). Accordingly, the SV and MM realize the system’s stable GS, respectively, at and , while in the opposite case each soliton species represents an excited state, instead of the GS, being therefore unstable.

#### iii.2.3 Effects of the Zeeman splitting

The above considerations did not include the Zeeman splitting, . In the presence of moderately strong splitting, with , the consideration of Eq. (22), similar to that which produced the asymptotic form (32) of the SV soliton at , demonstrates that the SVs exist at

(40) |

cf. existence region (33) at . The asymptotic form of the SV solution, affected by , is more complex than one given by Eq. (32); namely, the radial decay rate and wavenumber in Eq. (32) are replaced as follows:

(41) |

Strong Zeeman splitting, with , replaces existence condition (40) by . In this case, the SV solitons keep the asymptotic form (41) in the semi-infinite interval (40) of the chemical potentials, while, in the adjacent finite interval appearing in this case,

(42) |

the SV soliton exhibits a more dramatic change of its asymptotic shape: the radial wavenumber vanishes, the decay of the solution at being provided by the following exponential factor,

(43) |

As concerns the switch of the stability between the SV and MM at , which was demonstrated above for , it is pushed by to larger values of the relative strength of the cross-attraction, Sherman (). An explicit result can be obtained in the limit of large , when Eq. (23) demonstrates that the chemical potential is close to :

(44) |

The spin-down component, , is vanishingly small in this limit, hence stationary equation (24) simplifies to

(45) |

where is approximately replaced by , pursuant to Eq. (44). Then, the substitution of approximation (45) into Eq. (23) leads to the following equation for :

(46) |

By itself, Eq. (46) is tantamount to the stationary version of the 2D NLSE, which gives rise to the Townes solitons in , while Eq. (45) generates a small vortex component of the SV complex in . Due to the presence of factor in Eq. (46) and the smallness of , the norm (25) of the corresponding SV complex is

(47) |

where is the standard norm of the Townes soliton (18), the last term being a contribution from the small vortex component given by Eq. (45). Equation (47) demonstrates that the total norm of the SV soliton is (slightly) smaller than the collapse threshold, , hence the SV soliton is protected against the collapse, still realizing the stable GS of the system.

The lowest-order approximation presented here does not give rise to terms including the cross-attraction coefficient, from Eq. (22), which implies that, in the limit of large , all solitons belong to the SV type, irrespective of the value of , in accordance with the above-mentioned fact that the increase of leads to conversion of the MMs into SVs Sherman ().

#### iii.2.4 Mobility of stable 2D solitons

Although the underlying system (22) conserves the momentum (27), 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 along the axis, but not along . The respective solutions moving at velocity can be looked for as

(48) |

(Eq. (48) produces the Galilean transform of the wave functions in Galilean-invariant systems). The substitution of ansatz (48) into Eq. (22) leads to the coupled GPEs in the moving reference frame, which differ from Eqs. (22) by the presence of linear mixing between the two components (here , and is set):

(49) |

The same linear mixing can be imposed onto the 2D settings by a GHz wave coupling the two underlying atomic states radio (), hence the linear mixing by itself represents a relevant addition to the model.

Stationary solutions to equations (49), written in the moving reference frame, can be obtained, as well as in the case of Eqs. (22), by means of the imaginary-time evolution method. In particular, at , when the GS is represented by the quiescent MM soliton, its moving version, which is displayed in Figs. 8(a,b) for and , exists and is stable too. As well as its quiescent counterpart, this mode features the mirror symmetry between the profiles of and . Figures 10(a,b) display a typical example of a stable moving MM, and Fig. 10(c) shows the amplitude of the moving soliton, , as a function of . The amplitude monotonously decreases with the growth of the velocity, vanishing at

(50) |

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

The SV solitons may also be made mobile, but in a very narrow interval of velocities – e.g., at for , and . At , the imaginary-time solution of Eq. (49) with the SV input converges to stable MM solitons, instead of the SV.

### iii.3 3D metastable solitons

#### iii.3.1 Analytical considerations

The creation of metastable 3D solitons in the model based on Eq. (29) can be predicted starting from scaling evaluation of different terms in the respective energy functional (30). Assuming that a localized state has characteristic size and norm , an estimate for the amplitudes of the wave function is . Accordingly, the three terms in Eq. (30) scale with as

(51) |

with positive coefficients , , and . As shown in Fig. 11, Eq. (51) gives rise to a local minimum of at finite , provided that

(52) |

Although this minimum cannot represent the GS, which formally corresponds to at in the collapsed state, as is suggested by Fig. 11 too – in fact, this means that the system has no true GS), it corresponds to a self-trapped state which is stable against small perturbations. Condition (52) suggests that the metastable 3D solitons may exist when the SOC term is present, while its strength is not too large, and being not too large either.

More accurate semi-analytical results can be obtained by means of the VA, which is based on the following ansatz for the 3D version of the SV mode, written in cylindrical coordinates , cf. its 2D counterpart (36):

(53) |

(54) |

with real parameters , and , . In particular, amplitudes and account for possibility of having even and odd terms in the dependence of the wave functions of the vertical coordinate, . Further, the ansatz for the 3D version of the MM soliton may be chosen as a superposition, with mixing angle , of the SV (53) and its flipped counterpart, cf. Eq. (35) in the 2D setting:

(55) |

The substitution of ansätze (53)-(55) into the full energy (30) and minimizing it with respect to the free parameters under the constraint of the fixed total norm produces results summarized in Fig. 12, in which (meta)stable 3D solitons are predicted to exist in the shaded areas. It is seen that the solitons indeed exist, provided that , and are not too large, in agreement with the rough prediction of Eq. (52). In other words, for fixed and the stable 3D solitons exist in a finite interval of the norm,

(56) |

which is qualitatively similar to the situation for the single-component model stabilized by the HO trapping potential, cf. Figs. 1 and 2.

As shown in Fig. 13(a), for the VA-predicted energy of the SV is lower than that for the MM, and vice versa for , similar to what was found in the 2D model (22) with . Furthermore, red squares in Fig. 13(b) depict the VA prediction for the soliton’s chemical potential, , plotted as a function of norm for and . These curves clearly demonstrate that, also similar to what was found in the 2D model (cf. Figs. 8(b) and 9(b)), there is no threshold (minimum norm) necessary for the existence of the stable 3D solitons, which exist up to a , see Eq. (56).

While the stability of the upper branch is consistent with the VK criterion, , the lower branch in Fig. 13(b), which does not satisfy the VK criterion, represents solitons corresponding to the energy maximum (instead of the minimum) of the blue dashed curve in Fig. 11, therefore the lower-branch solitons are definitely unstable. In the limit of , they carry over into the known unstable 3D solitons of the single GPE Silberberg ().