Multipoles and vortex multiplets in multidimensional media with inhomogeneous defocusing nonlinearity

We predict a variety of composite quiescent and spinning two- and three-dimensional (2D and 3D) self-trapped modes in media with a repulsive nonlinearity whose local strength grows from center to periphery. These are 2D dipoles and quadrupoles, and 3D octupoles, as well as vortex-antivortex pairs and quadruplets. Unlike other multidimensional models, where such complex bound states either do not exist or are subject to strong instabilities, these modes are remarkably robust in the present setting. The results are obtained by means of numerical methods and analytically, using the Thomas-Fermi approximation. The predicted states may be realized in optical and matter-wave media with controllable cubic nonlinearities.


I. INTRODUCTION
ticular, we demonstrate that there is a bifurcation which destabilizes 2D dipoles, replacing them by stable VAPs, and another bifurcation, which replaces isotropic 2D VAQs by stable anisotropic modes of the same type (AVAQs). Thus, the present model offers the simplest setting in which stable multipoles and vortex-antivortex multiplets can be constructed in the multidimensional space. In this respect, it may be compared to field-theory models, where soliton complexes are created in sophisticated multi-component settings [7,35].

II. THE GOVERNING EQUATIONS
The basic model is represented by the scaled D-dimensional Gross-Pitaevskii equation (GPE) which governs the evolution of the mean-field wave function, ψ(r, t), in BEC: Here t is time, ∇ 2 is the 3D or 2D Laplacian, and σ(r) represents the spatially growing local strength of the repulsive nonlinearity. Following Ref. [30], we here adopt the steep modulation profile, σ(r) = exp(r 2 /2r 2 0 ), where we set r 0 ≡ 1 by means of straightforward rescaling. As well as in Ref. [30], milder profiles, characterized by σ(r) growing as r α (with α > D, as said above) are possible too, but the anti-Gaussian one makes it possible to present results in a compact form. Dynamical invariants of Eq. (1) are the norm, the angular momentum, and the Hamiltonian: In the 2D reduction of the model, the vectorial angular momentum (3) is replaced by its single component, The 2D version of Eq. (1), with time substituted by the propagation distance (z), applies as well to the spatialdomain evolution of the amplitude of electromagnetic waves in a bulk optical waveguide [1] with the self-defocusing cubic nonlinearity. In that case, the 2D reduction of Eq. (2) determines the total power of the optical beam, while Eq. (5) is the beam's orbital angular momentum [43].
Stationary states with real chemical potential µ > 0 (in optics, −µ is the propagation constant) are looked for in the usual form, where the (generally, complex) spatial wave function obeys the equation in the 2D setting, with angular coordinate ϕ. The 3D stationary equation is wherel 2 is the usual angular part of the 3D Laplacian.

III. MULTIPOLES: NUMERICAL RESULTS AND THE THOMAS-FERMI APPROXIMATION (TFA)
Stationary solutions of Eq. (7) and Eq. (8) for multipole modes were obtained by means of the imaginary-time method [44], which is capable to generate not only the ground state but, under special conditions, higher-order modes too [45], and also by means of the modified squared-operator method [46]. The respective inputs, emulating the multipoles sought for, were built as combinations of identical Gaussians with separated centers and alternating signs. The stability of the so generated modes was examined via direct simulations of the perturbed evolution, using the standard pseudospectral split-step fast-Fourier-transform method. 2D simulations were performed in the domain of size (12π) 2 , covered by a mesh of 512 2 points. For the 3D simulations of the octupoles, the split-step algorithm was used too, in the domain of size (12π) 3 , with the mesh of 192 3 points. In particular, the course of the evolution, under the influence of the small numerical noise due to the discrete nature of the mesh, the octupole with µ = 2 kept its 3D shape virtually invariant for t = 200, which exceeds 10 respective diffraction lengths.
Numerically found shapes of the 2D dipole (m = 1) and quadrupole (m = 2) modes, which are displayed in Figs. 1(a) and (b), suggest that these modes may be approximated, in the simplest form, by a real ansatz, The substitution of this into Eq. (9), and the use of the usual mean-field approximation, cos 3 (mϕ) → (3/4) cos (mϕ), yields the following radial equation: The application of the TFA to Eq. (10) implies neglecting the derivative terms, which yields the simple result: This approximation predicts that maxima of the amplitude are located at distance from the origin Amplitude maxima of the numerical solutions, found for µ = 2, are located at distances r  (12) and (17) (see below) predict, for the same case (µ = 2), r (max) dipole = 1.13, r (max) quadrupole = 1.8, and r (max) octupole = 2.75. Naturally, the agreement improves for larger µ, i.e., larger N , as the TFA is more accurate for stronger nonlinearity (the increase of µ from 2 to 10 leads to a decrease of the relative difference between the numerical solutions and their TFA counterparts by a factor 2).
Families of the modes are characterized by the norm (total power) as a function of µ. The TFA based on Eq. (11) produces the following dependence: where ρ ≡ r 2 /2, and the integral can be easily computed numerically. In the limit of µ → ∞, Eq. (13) yields an asymptotically constant slope, Dependences (13) for m = 1 and m = 2 are displayed in Fig. 2(a), along with their counterparts produced by full numerical solutions of Eq. (7).
In the 3D geometry, the typical shape of the octupoles, displayed in Fig. 1(c), suggests that the corresponding angular structure may be approximated by spherical harmonics with quantum numbers l = 3, m = ±2 [47]: where the constant C is not essential here, ϑ is the second angular coordinate, andl 2 Y 32 = 12Y 32 . Thus, the ansatz for the octupole is adopted as For the application of the mean-field approximation to the angular dependence in the cubic term of Eq. (8), The eventual results produced by the 3D version of the TFA are cf. Eq. (11), cf. Eq. (12), and cf. Eq. (13), with cf. Eq. (14). The plot produced by Eq. (18) is displayed in Fig. 2(b), along with its counterpart generated by the full numerical solution. An essential finding revealed by the systematic numerical analysis of the 2D setting is that the increase of µ (and N ) leads to destabilization of the 2D dipole mode at a critical point, by a bifurcation, which gives rise to a new mode in the form of a VAP, which will be presented in the next section. Further, careful analysis demonstrates that both 2D quadrupoles and 3D octupoles are, strictly speaking, entirely unstable solutions. Nevertheless, they are robust objects, provided that their norms are relatively small, as the instability is very weak in that case. For the quadrupoles, the simulations reveal that the robustness interval is µ < µ max 3, corresponding to N < N max 1.55, see Fig. 2(a). In this interval, initially perturbed quadrupoles maintain their shape as long as the simulations were run (with the evolution times measured in many dozens of characteristic dispersion times), exhibiting small-amplitude oscillations between the unperturbed shape and a coexisting stable isotropic VAQ mode (the latter one is presented in detail below). For the octupoles a similarly defined robustness border is located at µ max 2, corresponding to a low value of the 3D norm, N 0.03.  2. (Color online) Lines labeled "Dipole", "Quadrupole", and "Octupole" show, respectively, the 2D and 3D norms, N , vs. µ for these modes, as obtained from numerical computations. The lines marked by "TF" display the predictions for the same dependences produced by the Thomas-Fermi approximation, as per Eqs. (13) and (18). The short solid segment of the N (µ) curve for the dipole-mode family designates a completely stable subfamily (cf. Fig. 5 below). The application of the torque to the stable 2D dipole mode, i.e., the multiplication of the respective stationary solution by the phase factor with real constants p and x 0 [32], readily sets it in persistent rotation, thus creating a robust "propeller" mode, cf.
Ref. [14], which preserves the original two-peak shape. The spinning regime is illustrated by a set of snapshots in Fig.  3, where the rotation period is T = 3.9. Rotation of the dipole was simulated up to t = 100, amplitude oscillations of the dipole in the course of rotation being limited to 2% of its initial value. On the other hand, the application of the torque to the quadrupoles and octupoles does not generate persistently rotating states, in accordance with the above conclusion that these species do not represent fully stable modes.

IV. VORTEX-ANTIVORTEX PAIRS (VAPS)
In addition to multipoles, bound pairs of vortices with opposite signs of angular momenta were investigated. A typical example of the VAP and the respective bifurcation diagram are displayed in Figs. 4 and 5, respectively. As seen in panel (a) of Fig. 4, the amplitude profile of the VAP is quite similar to that of the dipole mode, but panel (b) of the same figure demonstrates the phase structure of the stationary VAP, which its real dipole-mode counterpart does not have. It is exactly the phase structure which makes it possible to identify the new mode as the vortex-antivortex bound state.
In accordance with what is said above, the stability of the modes under the consideration was identified by means of systematic simulations of their perturbed evolution. As shown in Fig. 5, the VAP is stable after the bifurcation point, viz., at 1.64 < µ < 4.06, 1.38 < N < 10.35, then it is destabilized by oscillatory (in t) perturbations in a finite region adjacent to the stability area, 4.06 < µ < 8.06, 10.35 < N < 30.76, which is followed by the eventual restabilization of the VAP at µ > 8.06, N > 30.76. Stable VAPs can be easily set in spinning motion, similar to what is shown above for dipoles.
The emergence of the VAPs determines the instability of the dipole modes. In the interval adjacent to the critical point Eq. (20), 1.64 < µ < 3, the dipole features a weak instability, periodically oscillating between its unperturbed shape and the coexisting stable VAP, whose amplitude profile is close to the dipole's one (not shown in detail here). A similar dynamical regime, featuring remittent shape revivals of a weakly unstable dipole state, was observed in a model with a nonlocal nonlinearity [15]. At µ > 3, the dipole's instability gets stronger, leading to its spontaneous transformation into a fundamental isotropic state (not shown in detail either).  Fig. 2(a)], and for the (partly stable) isotropic and (fully unstable) anisotropic vortex-antivortex quadruplets (VAQs). As usual, stable and unstable families of solutions are indicated by solid and dotted lines, respectively. The dashed-dotted and dashed lines refer, severally, to the Thomas-Fermi approximation for the quadrupoles and the anisotropic VAQs. The bifurcation point (22), at which the anisotropic-VAQ branch splits off from its isotropic counterpart, is marked by the red dot.

V. VORTEX-ANTIVORTEX QUADRUPLETS (VAQS)
Similar to the 2D dipoles, which coexist with VAPs, real stationary 2D quadrupoles coexist with a branch of stable complex solutions for vortex-antivortex bound states in the form of the (isotropic) VAQ, see a typical example of the latter in Fig. 6. However, Fig. 7 demonstrates the difference from the situation for the coexistence of the dipoles and VAPs: the isotropic-VAQ branch emerges at µ = 0, rather than branching off from the quadrupole one at a finite value of µ, cf. Fig. 5. Recall that, unlike the dipoles, quadrupole solutions are, strictly speaking, always unstable, hence, indeed, a new stable branch cannot bifurcate from them. As shown in Fig. 7 , the isotropic-VAQ family is stable at 0 < µ < 5.74 (which corresponds to 0 < N < 11.06). At µ > 5.74 (N > 11.06), these modes become unstable against oscillatory perturbations, eventually evolving into the isotropic ground state (not shown here in detail). Further A typical example of an AVAQ mode is displayed in Fig. 8. The analysis demonstrates that this solution is always unstable [this fact explains why the bifurcation occurring at point (22) does not destabilize the branch of the isotropic  VAQs]. Specifically, from the bifurcation point (22) up to µ = 12.54 (N = 40.77), a perturbed AVAQ solution starts spinning motion, with an angular velocity determined by the initial perturbation, see an example in Fig. 9 (in this case, the conservation of the angular momentum is provided by the recoil effect of small-amplitude waves emitted by the perturbed AVAQ, which are not visible in Fig. 9). Eventually, the unstable spinning AVAQ transforms not into the isotropic ground state, which is typical for other species of unstable modes which were mentioned above, but, instead, into a stably rotating VAP. In several cases that were examined, this transition occurred before the AVAQ would complete half a cycle of its rotation. Finally, at µ > 12.54 (N > 40.77), the AVAQ immediately transforms into a rotating VAP, without going through the intermediate spinning stage (not shown here in detail).

VI. CONCLUSION
The recently introduced class of optical and matter-wave models, with the strength of the self-repulsive cubic nonlinearity growing from the center to periphery, gives rise to a large variety of completely stable complex 2D and 3D self-trapped states (solitons), which do not exist or are strongly unstable in other physical settings. Here we have demonstrated that, together with previously found states, the same models support several other species of 2D and 3D solitons which are not available (in a stable form) in other systems either. These are 2D dipoles and quadrupoles and 3D octupoles, as well as VAPs (vortex-antivortex pairs) and isotropic or anisotropic VAQs (vortex-antivortex quadruplets). The modes found here are stable or weakly unstable, surviving for a long time (over a long propagation distance), which makes them physically relevant objects. In addition to quiescent states, persistently spinning dipoles have been found too. The results, which were obtained by means of numerical methods and analytically, using the TFA (Thomas-Fermi approximation), demonstrate the potential of the settings with the spatially growing strength of the self-repulsion for the creation of a great variety of complex stable multidimensional modes.
Natural extensions of the present analysis may be additional analysis of the 3D setting (in particular, for 3D dipoles and quadrupoles), and its development for two-component systems, in which the multitude of nontrivial self-trapped states may be further expanded, as suggested, e.g., by results for the two-component system in 1D [48].