Three-dimensional hybrid vortex solitons

We show, by means of numerical and analytical methods, that media with a repulsive nonlinearity which grows from the center to the periphery support a remarkable variety of previously unknown complex stationary and dynamical three-dimensional solitary-wave states. Peanut-shaped modulation profiles give rise to vertically symmetric and antisymmetric vortex states, and novel stationary hybrid states, built of top and bottom vortices with opposite topological charges, as well as robust dynamical hybrids, which feature stable precession of a vortex on top of a zero-vorticity base. The analysis reveals stability regions for symmetric, antisymmetric, and hybrid states. In addition, bead-shaped modulation profiles give rise to the first example of exact analytical solutions for stable three-dimensional vortex solitons. The predicted states may be realized in media with a controllable cubic nonlinearity, such as Bose-Einstein condensates.


I. INTRODUCTION
Self-trapping of three-dimensional (3D) confined modes (solitons or, more properly, solitary waves) in optics [1][2][3], Bose-Einstein condensates (BECs) [4][5][6], ferromagnetic media [7], superconductors [8], semiconductors [9], baryonic matter [10], and general field theory [11,12] is a fundamental problem of nonlinear physics. An apparent condition is that an attractive, or self-focusing, nonlinearity is necessary for the creation of localized states; however, the attractive cubic nonlinearity simultaneously gives rise to collapse [13] of localized modes in higher-dimensional settings and, additionally, to strong azimuthal modulational instability of states with intrinsic vorticity [14], thus making the search for stable 3D fundamental and topological solitons in materials with the cubic (Kerr) nonlinearity a challenging issue.
Various methods have been elaborated over the years, chiefly in the theoretical form, to remedy this situation and stabilize 3D solitary waves, fundamental and vortical ones alike. As outlined in detail in the reviews [1,2] (see also the more recent work [15]), stabilization may be achieved by a higher-order quintic self-defocusing nonlinearity, provided that the underlying physical setting gives rise to such terms. Another possibility is offered by periodic (lattice) potentials [1][2][3]. In particular, a 2D potential may be sufficient for the stabilization of 3D solitons, as well as for the stabilization against supercritical collapse [16]. In addition, it is also possible to stabilize 3D fundamental solitons by means of "nonlinearity management" (time-periodic sign-changing modulation of the nonlinearity coefficient), which should be combined, at least, with a 1D lattice potential [17]. The use of nonlocal nonlinearities may also help to stabilize 3D localized modes [19]. Lastly, it is relevant to mention a very recent result concerning 2D localized modes created by the self-focusing cubic nonlinearity in the free space: while a common belief was that they might never be stable, it has been demonstrated in Ref. [20] that mixed vortex-fundamental modes in a system of two coupled GP equations modeling the spin-orbit-coupled BEC can be stable in the 2D free space. This unexpected result is explained by the fact that the norm of the mixed modes takes values below the well-known 2D-collapse threshold [13].
Unlike the above-mentioned methods, the use of spatially inhomogeneous cubic nonlinearity does not yield stabilization of 3D solitons [3]. In the 2D setting, a nonlinearity subject to a smooth spatial modulation cannot stabilize solitons either [18]. Stabilization of 2D fundamental solitons (but not vortex states) is possible by means of various spatial modulation profiles with sharp edges [21]. For this reason, most of previous studies of solitons in inhomogeneous nonlinearity landscapes have been performed in 1D settings, chiefly for periodic modulation patterns [22].
A radically different approach was recently put forward and elaborated in Refs. [23] and [24]: a repulsive, or defocusing, nonlinearity, whose local strength grows from the center to the periphery, as a function of radius r at any rate faster than r 3 , can readily induce self-trapping of robust localized modes, which are stable not only to weak, but also to strong perturbations (although these solutions are far from those in integrable models, we call them "solitons", as commonly adopted in the current literature when dealing with stable self-trapped modes). In BECs, the necessary spatial modulation of the nonlinearity may be induced by means of the tunable Feshbach resonance, controlled by magnetic [25] and/or optical [26] fields, created with appropriate inhomogeneous profiles [27]. The required magnetic field patterns can be provided by magnetic lattices of various types [28], while the optical-intensity profiles can be painted by laser beams in 3D geometries [29]. In addition to fundamental solitons, landscapes with a growing repulsive nonlinearity were shown to support topological states in the form of vortex-soliton tori, which can exhibit gyroscopic precession under the action of an external torque [24] (precession of a tilted vortex was earlier considered in a different setting in Ref. [30]).
So far, only the simplest 3D vortex solitons were addressed in the framework of the setting based on the spatially modulated strength of the self-repulsion. The possibility of the existence of more complex vorticity-carrying 3D structures remains unexplored. In this context, it should be stressed that the creation of stable structures carrying several topological dislocations is a complex challenge. Previously, such entities were found mostly in the form of vortex-antivortex pairs and vortex arrays in settings with a reduced dimensionality, such as superconductors [31,34], pancake-shaped atomic Bose-Einstein condensates [32,35], and exciton-polariton condensates [33]. To the best of our knowledge, no examples of 3D solitons with coaxial vortex lines threading several objects forming a complex state, or with the topological charge changing along the axis of the soliton, have been reported.
In this work, our analysis reveals that 3D media with a repulsive nonlinearity growing from two symmetric minima to the periphery make it possible to create complex but, nevertheless, stable static and dynamical self-trapped topological modes, in the form of fundamental and vortical dipoles, stationary vortex-antivortex hybrids, and precessing hybrids built as a vortex sitting on top of a zero-vorticity mode. These are remarkable, novel species of 3D localized modes, which have not been reported before in any other systems. The very existence of the stationary vortex-antivortex solitons and precessing vortex-fundamental hybrids is an unexpected finding, because the topology of such states is different in their top and bottom sections. All these previously unknown static and dynamical states are supported by the nonlinearity-modulation profile, which is obtained from the spherical configuration by a deformation in the axial (vertical) direction.
The basic model is introduced in Section II, where we also give a number of analytical results, which can be obtained in spite of the apparent complexity of the system. These include the Thomas-Fermi approximation (TFA) for families of vortex modes, an approximate description of the dipole (antisymmetric) modes in terms of quasi-1D dark solitons embedded into the ordinary symmetric states, and an approximation which explains the existence of stationary vortex-antivortex hybrids. Results of systematic numerical analysis are reported in Section III, including families of stationary antisymmetric and vortex-antivortex hybrid modes, as well as dynamical (steadily precessing) vortexfundamental hybrids. A comprehensive stability analysis is presented too, along with simulations of the spontaneous evolution of unstable states. The work is concluded by Section IV. In the Appendix, we additionally present stable analytical solutions for 3D vortex solitons in a model with a bead-shaped spatial modulation profile, which is the first example of any system admitting exact solutions of this type, thus providing a direct proof of their existence.

II. THE MODELS AND ANALYTICAL RESULTS
A. The general formulation Our system is described by the single-component nonlinear Schrödinger/Gross-Pitaevskii (NLS/GP) equation in the 3D space for the wave function ψ(r, t): where Laplacian ∇ 2 acts on coordinates r = {x, y, z}, and σ(r) > 0 represents the local strength of the repulsive nonlinearity, which must grow at r → ∞ faster than r 3 . Dynamical invariants of Eq. (1) are the norm and Hamiltonian, where the (generally, complex) spatial wave function satisfies the equation While the simplest 3D vortex solitons have been obtained in spherically-symmetric nonlinearity landscapes, such as the one with σ(r) = exp r 2 /2 [24], here our objective is to a show that a deformation of this nonlinearity profile, lending it two local minima, allows us to produce novel species of robust stationary and precessing 3D topological modes. To this end, the spherically symmetric modulation pattern is shifted by the distance ±d/2 along the z axis, and the so produced profiles are stitched together in the midplane, z = 0: with ρ 2 ≡ x 2 + y 2 . This profile keeps the cylindrical symmetry and, accordingly, the z-component of the field's angular momentum, which is the third dynamical invariant of the model, in addition to N and H, where * stands for the complex conjugation. The steep anti-Gaussian profile, adopted in Eq. (3), is not a necessary feature of the model. As mentioned above, the necessary condition for the existence of 3D solitons, which follows from the normalizability of the wave function, is that σ(r) must grow faster than r 3 [23]. The modulation profile (3) is adopted here as it makes it possible to obtain families of stationary vortex modes in an almost exact analytical form, by means of the TFA, thus supporting numerical findings.

B. Symmetric self-trapped vortices and the Thomas-Fermi approximation (TFA)
Among the complex stable modes reported below, the simplest species are confined vortex states, carrying an integer topological charge S, which are looked for, in the cylindrical coordinates, as where Φ is a real function. As follows from Eqs. (4) and (5), the angular momentum of the vortex is M = SN . Below, such modes, with identical vorticities S in the top and bottom parts of the peanut-shaped nonlinearity landscape, are denoted as S/S (the definitions of "top" and "bottom" are arbitrary here, as Eqs. (2) and (3) are obviously invariant with respect to z → −z). The shape of the simplest symmetric vortices and fundamental solitons (S = 0) can be approximated by means of the TFA, which neglects z-and ρ-derivatives in Eq. (2), and is usually relevant in the case of a strong repulsive nonlinearity [5,36,37]: Here the first line represents the hole at the center of the vortex state (see panels marked 1/1 in the top rows of Figs. 1 and 2). Families of self-trapped modes are characterized by dependence N (µ), which can be obtained from Eq. (6) in an approximate analytical form: For S = 0 (the fundamental mode), Eq. (7) reduces to a simple linear dependence, N The constant slope dN/dµ given by the latter expression is, actually, an asymptotically exact result at large µ for any S. Figures 3(a,b) show that, while the TFA predictions for N (µ) may be shifted from their numerically found counterparts, the asymptotic slope is indeed predicted exactly.
Our stability analysis for vortex modes (5), as well as for other stationary modes featuring the cylindrical symmetry, which are considered below, was carried out by numerically solving the linearized equations for small perturbations. Perturbed solutions are sought for as where ǫ is an infinitesimal amplitude of the perturbation, k is its integer azimuthal index, and δ (S, µ, k) is a (generally, complex) instability growth rate. Substitution of expression (8) into Eq. (1) and the linearization gives rise to the eigenvalue problem for δ represented by the following equations: The stability condition is Re {δ (S, µ, k)} = 0, which must hold for all eigenvalues at given values of S and µ.

C. Dipole (antisymmetric) modes
The vortex and fundamental modes can be twisted in the vertical direction, which lends them an antisymmetric (dipole) structure along the z axis, as depicted in the left and middle panels in the bottom rows of Figs. 1 and 2. Dipole modes have been previously studied in diverse 1D and 2D settings [38], including vortex dipoles created in a common plane [31][32][33]. In 3D, such dipole structures can be approximately described by assuming that a quasi-1D dark soliton is embedded into an originally symmetric 3D mode around its midplane (z = 0), as suggested in a different context in Ref. [6]. In particular, for the fundamental states (S = 0) approximated by the TFA expression (6), the respective antisymmetric solution can be easily found from Eq. (2), assuming that the width of the dark soliton in the z direction is much smaller than the intrinsic scale of the TFA mode, i.e., µ is large enough: For the vortex states, a similar approximation is available too, but its applicability condition does not hold around the inner hole of the vortex. Solution (10) corresponds to a gap which cleaves the antisymmetric mode around z = 0, as shown in the bottom row of Figs. 1 and 2. The width of the gap does not depend on ρ, implying that the gap is nearly flat, which is well corroborated by numerical results, see the left panel in the bottom row of Fig. 2. Solution (10) makes it possible to calculate the difference between the norm of the symmetric state and its antisymmetric counterpart. Indeed, Eqs. (6) and (10) yield As shown in Fig. 3(d), this prediction is quite accurate.

D. Hybrid modes
Completely novel species of stationary 3D modes are hybrids of the S/ − S type, which combine vortex states with opposite signs and equal norms in the top and bottom sections of the peanut-shaped structure, as shown in Figs. 1 and 2. Unlike the symmetric and antisymmetric vortices introduced above, the hybrids cannot feature axisymmetric density distributions. A central question is whether the vortex-antivortex hybrids exist as stationary modes and, if they do, whether they can be stable. To address this issue, a stationary solution may be looked for in an approximate form as assuming that φ + (ρ, z) and φ − (ρ, z) rapidly vanish, respectively, at z < 0 and z > 0, so that the two vortical components form a sharp domain wall close to z = 0. Substituting ansatz (12) in Eq. (2), and using the rotatingwave approximation, one arrives at a system of nonlinearly coupled equations, Note that in the right-hand sides of this equation the cross-phase-modulation coefficient is twice as large as its selfphase-modulation counterpart. This is typical for systems which give rise to solutions in the form of sharp domain walls between states with different wave numbers, linear or azimuthal ones [37,39]. Although Eq. (13) is axisymmetric, as the angular coordinate θ does not appear in it, the superposition of the two vortices in Eq. (12) breaks the isotropy of the pattern in the midplane: |φ (ρ, θ, z = 0)| 2 = 4φ 2 0 (ρ) cos 2 (Sθ), where φ + (ρ, z = 0) = φ − (ρ, z = 0) ≡ φ 0 (ρ). The latter pattern is close to the numerically found midplane structures, as can be seen in the right column of Fig. 1.
It is relevant to stress that, unlike the vortical modes of the S/S type considered above, the vortex-antivortex hybrids cannot be classified as symmetric or antisymmetric species, with respect to the top and bottom sections of the "peanut" profile. Indeed, a rotation of an hybrid state by angle π/2 about the vertical axis is effectively tantamount to adding a phase shift of π between the top vortex and the bottom antivortex.
Another novel type of hybrid modes, which is studied by means of direct simulations below, is one of the S = 1/0 type. In this case, the ansatz in the form of the superposition of the vortical (S = 1) and fundamental (S = 0) modes in the top and bottom sections of the system [cf. Eq. (12)] does not lead to a self-consistent approximation. In this situation simulations reveal robust dynamical regimes, with the vortex precessing on top of the fundamental soliton, as illustrated in Fig. 4. Our simulations show that, in suitable parameter regions, such spontaneously established dynamical states survive over indefinitely long evolution times (far exceeding t = 100).

A. Stationary modes and their stability
Stationary solutions for the basic types of 3D confined modes that are defined above were obtained as solutions of Eq. (2), with the modulation function (3), by means of the Newton's method. The stability of the so generated families of different modes was studied by means of a numerical solution of eigenvalue problem (9), and verified by direct simulations of perturbed evolution of the modes that were performed with the help of the split-step algorithm.
As indicated above, the solution families are naturally represented by dependences N (µ), which are collected in Fig. 3 for two values of d in Eq. (3), viz., d = 3 in (a,b), and d = 5 in (c). The plots distinguish stable and unstable families, and include the analytical results presented above, viz., the prediction of the TFA for the symmetric modes of the S = 0/0 and S = 1/1 types (see Eq. (7)). In addition, the norm difference between the symmetric and antisymmetric S = 0/0 states, as predicted analytically by Eq. (11), is presented, together with its numerically computed counterpart, in panel (d) for d = 0.
Typical examples of all stationary modes are displayed in Fig. 1, their shapes being additionally illustrated by means of vertical cross sections in Fig. 2. Antisymmetric 3D modes of the 0/0 and 1/1 types seem as built of two oblate fundamental solitons or vortices "levitating" on top of each other. Symmetric 0/0 and 1/1 states, which feature "peanut"-like shapes, transform into solutions reported in Ref. [24] with the decrease of separation d between the nonlinearity minima.
A salient finding is the existence of the stationary hybrid modes, stable and unstable examples of which are shown, respectively, for S = 1/ − 1 and S = 2/ − 2. Cross sections of the hybrids are displayed in the right column of Fig.  2, along the nodal directions in the midplane (z = 0). Such a choice of the presentation is required because, as noted above, the hybrid modes are axially asymmetric, in contrast to the isotropic ones of types 0/0 and 1/1.
As concerns the stability of the modes, all branches in Fig. 3 satisfy the anti-Vakhitov-Kolokolov criterion, dN/dµ > 0, which is a necessary (but, generally, not sufficient) condition for the stability of self-trapped states supported by repulsive nonlinearities [40] (the Vakhitov-Kolokolov criterion per se, dN/dµ < 0, is a necessary condition for the stability of solitons in media with attractive nonlinearities [13,41]). While the families of the symmetric modes of the S = 0/0 and S = 1/1 types were found to be completely stable, only small segments [the bold black ones in Fig.  3(a,b)] of their antisymmetric counterparts are stable too.
The stability-instability transition for the antisymmetric 0/0 and 1/1 states at small values of µ is additionally illustrated by Fig. 5, which displays the instability growth rates, δ r ≡ Re(δ), as functions of µ and azimuthal index k (limited to k ≤ 5), see Eqs. (8) and (9). In particular, an unusual peculiarity is that, for the antisymmetric state of the 0/0 type, the dominant instability mode for small µ corresponds to k = 1 [the red curve in Fig. 5(a)], while zero-vorticity states are normally destabilized solely by perturbations with k = 0 [1]. These instability eigenvalues are complex, hence the respective dynamics is oscillatory (see below).
Another important finding is a large stability region of the hybrid modes with S = 1/ − 1, as shown, in Fig. 3(c), for d = 5 in Eq. (3). It is worthy to note that this stability region strongly depends on d: a detailed analysis reveals that the vortex-antivortex hybrids are completely unstable at d ≤ 4, when the vortex and antivortex constituents of the hybrid are relatively strongly pressed onto each other, and a stability region appears at d > 4, being µ ≤ 15.8, i.e., N < 394.9, at d = 4.5, and µ ≤ 13.5, i.e., N < 329.3, at d = 5. Thus, it is worthy to note that the size of the stability region does not grow monotonously with the increase of d. Typical examples of the evolution of perturbed modes, of those types which may be unstable [they are marked by circles on red branches in Fig. 3(a,b,c)], are displayed in Fig. 6. In all the cases, the evolution keeps initial values of the norm and angular momentum (4). In particular, weakly unstable antisymmetric (dipole) modes with S = 0/0 and S = 1/1, which are taken close to the boundary of the stability region [see Figs. 3 (a,b) and 5], feature only small oscillations of their amplitude, while keeping their dipole structure and vorticity (in the case of S = 1/1). That is, the regions of effective stability for the dipole modes are actually larger than the rigorously defined bold black segments on the respective N (µ) curves in Figs. 3(a,b). On the other hand, at greater values of N , stronger instability destroys the dipole structure, tending to transform the antisymmetric modes into their symmetric counterparts, as shown in the top and middle rows of Fig. 6.
A remarkable feature of the instability-induced evolution (well corroborated by the simulations) is that the vortical structure survives in the course of the spontaneous transformation of the unstable dipole mode of the 1/1 type into its stable symmetric counterpart (see the middle row in Fig. 6). As concerns unstable hybrids, they, quite naturally, exhibit spontaneous annihilation of the vortex with antivortex, thus gradually transforming themselves into symmetric zero-vorticity (fundamental) states, as seen in the bottom row of Fig. 6. On the other hand, stable hybrid solitons do not show any conspicuous shape transformations even at t > 10 3 , and even in the presence of strong initial perturbations.
As indicated above, hybrids with S = 1/0, built of a vortex placed on top of a fundamental mode, cannot form a stationary state. Nevertheless, direct simulations show, as shown in Fig. 4, that the hybrids of this type readily self-trap in a dynamical form, with the vortex performing periodic precession above the zero-vorticity base. The respective initial configuration was constructed by juxtaposing the top and bottom components taken as respective parts of the symmetric vortex and fundamental states, with S = 0/0 and S = 1/1, which were preliminarily generated, for equal values of the chemical potential, in the same trapping configuration. A systematic numerical analysis shows that such robust dynamical regime is observed in a broad parametric area, provided that d is not too small, namely, d ≥ d min ≈ 4.8.

IV. CONCLUSIONS
Using a systematic numerical analysis and a range of analytical approximations, we have discovered several previously unknown species of self-trapped complex 3D field states, that are supported by the local strength of a repulsive cubic nonlinearity growing from two local minima to the periphery, along the axial and radial directions alike. We have shown that the corresponding axisymmetric "peanut"-shaped 3D nonlinearity-modulation profiles support families of vortex states, which are both symmetric and antisymmetric with respect to the top-bottom reflection. The same system gives rise to a novel species of stable stationary top-bottom vortex-antivortex hybrids, which was not reported previously in any 3D setting, to the best of our knowledge. Another newly found species of self-trapped robust dynamical hybrid states exhibits stable precession of a top vortex above a bottom fundamental mode. In addition, we showed (in the Appendix) that systems with "bead"-shaped 3D modulation profiles produce the first example of exact analytical solutions for stable 3D vortex solitons. Settings of such type may be realized in media that allow a local control of the cubic self-repulsive nonlinearity by means of external fields. In particular, this is possible in Bose-Einstein condensates, using the Feshbach resonance controlled by appropriately designed nonuniform magnetic or optical fields. The latter settings suggest a physical realization of the predicted self-trapped modes.