Flat-floor bubbles, dark solitons, and vortices stabilized by inhomogeneous nonlinear media

We consider one- and two-dimensional (1D and 2D) optical or matter-wave media with a maximum of the local self-repulsion strength at the center, and a minimum at periphery. If the central area is broad enough, it supports ground states in the form of flat-floor \textquotedblleft bubbles", and topological excitations, in the form of dark solitons in 1D and vortices with winding number $m$ in 2D. Unlike bright solitons, delocalized bubbles and dark modes were not previously considered in this setting. The ground and excited states are accurately approximated by the Thomas-Fermi expressions. The 1D and 2D bubbles, as well as vortices with $m=1$, are completely stable, while the dark solitons and vortices with $m=2$ have nontrivial stability boundaries in their existence areas. Unstable dark solitons are expelled to the periphery, while unstable double vortices split in rotating pairs of unitary ones. Displaced stable vortices precess around the central point.


I. INTRODUCTION
The consideration of models based on the nonlinear Schrödinger (NLS) equations with a spatial variation of the local coefficient in front of the cubic term makes it possible to essentially expand the variety of stable localized and extended states supported by the NLS equations. In many cases, the variation is represented by nonlinear lattices, i.e., spatially periodic modulations of the local nonlinearity strength [1]. These arrangements help to stabilize multidimensional and spatiotemporal self-trapped states in optics [2,3], photonics [4,5], and various setups for Bose-Einstein condensates (BECs) in ultracold atomic gases [6][7][8]. The interplay of such structures with parity-time symmetry [9] and gauge fields [10,11] was considered too. In this connection, it is relevant to mention that linear lattices, i.e., spatially periodic linear potentials, are commonly used for the creation and stabilization of solitons. In particular, linear lattices in combination with self-repulsive cubic nonlinearity [12] give rise to various families of gap solitons, including fundamental [13,14], subwavelength [15,16], parity-time-symmetric [17,18], subfundamental [19], composite [20,21], surface [22], moving [23,24], dark [25], quasi-discrete [26,27], and multipole [28] ones, as well as clusters built of them [29]. Further, gap solitons were predicted in moiré lattices [30,31] and in systems with quintic and cubic-quintic nonlinearities [32]. Twodimensional (2D) gap solitons [33] and vortex solitons of the gap type [34] were predicted too.
Nonlinear lattices [1] can also support many types of soliton families, including those of fundamental, dipole, and component hopfions [59] and hybrid vortices [55].
The objective of this work is to extend the potential use of the spatially modulated self-repulsive nonlinearity. There are physically relevant species of nonlinear modes that have not been studied, in this context, in previous works. One of them represents "bubbles", i.e., modes with a suppressed density in the central region, which, unlike spatially odd dark solitons, are spatially even states that do not cross zero. The bubbles were originally introduced in usual models with the uniform nonlinearity [65]. Collisions between bubbles [66] and a possibility of collapse dynamics in them [67] were studied too. Models with spatially inhomogeneous self-defocusing, which is strong near the center and weak at the periphery, are very natural settings for the study of bubbles.
Another important class of the modes is one with flat-top shapes. They have been studied, theoretically and experimentally, in many settings with competing nonlinearities, such as cubic-quintic optical media [68] and "quantum droplets", predicted in BEC with the attractive cubic nonlinearity balanced by the effective quartic self-repulsion induced by quantum fluctuations around mean-field states [69,70]. The flat-top shape is explained by the fact that the competition imposes an upper limit on the density of the wave field, hence increase of the total norm of the field leads to the expansion of the mode, keeping its amplitude constant. In particular, the upper limit imposed on the density makes the BEC in the "quantum droplets" an incompressible liquid, which explains the name of these states. The droplets have been created in the experiment [71][72][73], and their theoretical investigation was extended for states with embedded vorticity [74,75].
Unlike the flat-top-shaped states, concepts of nonlinear modes whose shape would demonstrate features such as a "flat floor" or "flat waist" were not elaborated yet, to the best of our knowledge. In this paper, we explore such possibilities and demonstrate that 1D and 2D flat-floor bubbles, 1D flat-waist dark solitons and 2D flat-waist vortex modes can be readily supported by inhomogeneous self-defocusing cubic nonlinearity, without the necessity to use composite competing nonlinearities. In addition to their interest to fundamental studies, these previously unexplored states may find applications, such as all-optical beam guiding [76,77] (using the inner "hole" in dark solitons, bubbles, or vortices as conduits steering embedded signal beams) and optical tweezers. In these contexts, the size of the hole is a crucially important parameter, which can be adjusted by means of the waist shaping.
The rest of the paper is organized as follows. The 1D and 2D models are formulated in Section II, which also includes some analytical results, such as ones based on the Thomas-Fermi approximation (TFA) and the limit case of the model with the delta-functional spatial modulation of the local nonlinearity strength. Numerical results for bubbles, dark solitons (in 1D), and vortices (in 2D) are reported in Section III. It also includes a qualitative analytical explanation for stability and instability of the dark solitons. The paper is concluded by Section IV. The nonlinear Schrödinger (NLS) equation, which describes the propagation of laser beams in self-defocusing nonlinear media, or effectively 2D matter waves in BECs (in the latter context, it is called the Gross-Pitaevskii (GP) equation [78][79][80]), is taken in the scaled form: Here E and z denote the field amplitude and propagation distance, respectively, (x, y) is the set of the transverse coordinates, with r ≡ x 2 + y 2 , and ∇ 2 = ∂ 2 x + ∂ 2 y is the paraxialdiffraction operator, in terms of optics, or the kinetic-energy operator for matter waves. In the latter case, E and z are replaced by wave function ψ and time t, respectively. In the 1D version of Eq. (1), with single coordinate x, G(r) is replaced by G(x).
The NLS equation (1) conserves the total power (alias norm), the angular momentum (in 2D, where it may be conveniently written in polar coordinates, r and θ), and the Hamiltonian, In 1D, expressions (2) and (4) are replaced by their straightforward counterparts. Numerical simulations performed in this work maintain the conservation of P , M , and H.
To construct flat-floor bubbles and flat-waist vortices and dark solitons, we define the profile of the self-defocusing nonlinearity with G(r) > 0 as where g 0 , r 0 , and g > g 0 are positive constants. In 1D, r and r 0 should be replaced by |x| and x 0 > 0. By means of rescaling, we set g 0 ≡ 1 in what follows below, while g > 1 represents the largest local strength of the self-defocusing nonlinearity, at r = 0, and r 0 is the width of the flat-floor section of the nonlinearity profile. Stationary solutions with real propagation constant k < 0 and embedded integer vorticity, m = 0, 1, 2, ... , are sought for, in polar coordinates (r, θ), as where the real stationary field amplitude U satisfies the equation The mode with m ≥ 1 may be considered as a vortex with winding number m embedded in the fundamental state (the 2D "bubble"), which corresponds to m = 0. The 1D stationary states are looked for as with real function U satisfying the equation The 1D counterpart of vortices is the dark soliton, which can be constructed by embedding the usual dark-soliton solution of the NLS equation with G = const [81,82], in the flat-top segment of the 1D fundamental state. In the uniform space, the dark soliton is characterized by its power defect, i.e., the difference of the integral power and its background value, corresponding to the constant asymptotic density, U 2 backgr = −k/G: Another characteristic of the dark soliton is the difference of its Hamiltonian and the background value, For the dark soliton in the uniform space, it is relevant to express ∆H in terms of ∆P , using Eqs. (11) and (12): Similarly to Eq. (11), one can define the power defect for the bubble, vortex, and dark soliton states in the 2D and 1D versions of the present model, with the nonlinear-modulation profile defined as per Eq. (5). The defect is a sum of two terms, which represent the difference between the actual solution and background values in the regions of r < r 0 (or |x| < x 0 , in 1D) and r > r 0 (or |x| > x 0 ): 2. The delta-functional profile in 1D The limit case of the 1D model, opposite to the flat-floor configuration, is one with the central region that is very narrow in comparison with the dark-soliton's width. In this limit, Eq. (9) is replaced by where δ(x) is the Dirac's delta-function, and the coefficient from the 1D version of Eq. (5). The interaction of field U with the delta-functional defect is accounted for by the following term in the Hamiltonian (the interaction potential), cf. Eq. (12): Note that Eq. (16) gives rise to an exact bubble solution, with offset ξ determined by the cubic equation for tanh √ −kξ : It is easy to see that Eq. (19) has a single solution for all positive values of γ. The power defect of solution (18), defined as in Eq. (11), is Note that, as it follows from Eqs. (19) and (20), the power defect satisfies condition d∆P/dk > 0, which is tantamount to the anti-Vakhitov-Kolokolov criterion. It is a necessary condition for stability of states supported by self-repulsive nonlinearities [40]. This fact suggests that the exact bubbles supported by the delta-functional nonlinearity-modulation profile may be stable (the Vakhitov-Kolokolov criterion proper, with the opposite sign, is a well-known necessary condition for the stability of localized states supported by self-attractive nonlinearities [83][84][85]). Indeed, direct simulations of perturbed evolution of these bubbles, performed in the framework of Eq.

B. Linearized equations for small perturbations
To explore the stability of the stationary solutions against small perturbations, we introduce the perturbed 1D solution as (21) where p(x) and q * (x) are small perturbations associated with an instability growth rate, λ. Substituting this ansatz in Eq.
(1) leads to the eigenvalue problem for λ, In the 2D case, the perturbed wave function with embedded integer vorticity, m, and an integer azimuthal index, n, of the perturbation is taken as which leads to the eigenvalue problem in the following form: Stationary perturbed solutions are stable if real parts of all the corresponding eigenvalues are zero, Re(λ) = 0.
C. The Thomas-Fermi approximation (TFA) TFA, which neglects the derivatives in Eq. (7), is a commonly known method for constructing stationary solutions [80]. In the framework of TFA, the approximate solution of Eq. (7) is  Im( ) Examples of the eigenvalue spectra at k = −15, g = 5, and values of x0 marked by symbols S1-S3 in panel (a): Obviously, Eq. (25) pertains to the bubble, while the 2D expression (26) predicts a 2D bubble for both m = 0 and vortices with m ≥ 1 (in the former case, the "hole" at r 2 < −m 2 /k is absent in the TFA solution). Note also that, for the radial-modulation profile (5), the inner flat area, r < r 0 , is completely covered by the zero part of TFA, see the bottom row in Eq. (26), for large values of the vorticity, The predictions produced by TFA for the 1D and 2D bubbles and vortices are compared below with numerical solutions. In 1D, the TFA does not directly produce dark solitons, but they can be introduced by means of a straightforward approximation combining Eqs. (10) and (25), which is relevant under condition √ −kx 0 ≫ 1.
The accuracy of TFA can be naturally characterized by calculating the power defects, defined by Eqs. (14) and (15), and comparing them to counterparts produced by the numerical solution of Eqs. (7) and (9). In particular, it is obvious from Eqs. (25) and (26) (with m = 0) that TFA predicts the power defect of the 1D and 2D bubbles to be proportional to |k|. For the 2D state with m = 0, it can be calculated analytically in the limit of r 0 = 0: For the 1D bubbles, TFA produces a simple result in the limit of g ≫ 1, The comparison of the TFA-predicted and numerically found ∆P 1D,2D is displayed below in Figs. 1(d) and 5(c1), respectively.

III. RESULTS AND DISCUSSION
The numerical solution of Eqs. (7) and (9)  finding algorithms for f (x) = 0, solved by means of the iteration process, x n+1 = x n −f (x n )/f ′ (x n ) where n denotes the iteration's number and f ′ stands for the first derivative [86]. In the case of the 2D equation (7), the method was applied in the polar coordinates, removing the artificial singularity at r → 0 by means of the appropriate boundary conditions, viz., dU/dr| r=0 = 0 for zero-vorticity modes, or U (r) ∼ r m ones with vorticity m ≥ 1. The input for these modes was taken as and setting m = 1 for the 1D dark solitons.
The stability of the stationary solutions was subsequently explored by solving eigenvalue equations (22) and (24), and the results were verified by direct simulations of Eq. (1) for the perturbed propagation, employing the finite-difference method for marching in time [87]. It was checked that the numerical mesh with ∆x = ∆r = 0.1 and ∆z = 0.002 was sufficient for producing fully reliable results (i.e., rerunning the calculations with essentially smaller ∆x, ∆r and ∆z does not change the findings). It is relevant to mention that numerical methods for producing stationary and dynamical solutions to NLS/GP equations with a spatially varying nonlinearity coefficients were developed and used in many previous works [1,29,[40][41][42][43], [49]- [63].

A. One-dimensional (1D) bubbles
Typical profiles of the flat-floor bubbles in 1D are displayed in Fig. 1. The transition of the ordinary bubbles into ones featuring the flat-floor profile, following the increase of width x 0 of the flat section in the nonlinearity profile, is displayed in Fig. 1(a). Naturally, the width of the flat segment in the bubble is determined by x 0 . The width is not sensitive to the strength of the nonlinearity in the flat section, g, while the minimum value of the field in the bubble mode depends on g, as seen in Fig. 1(b). Profiles of the flat-floor bubbles with different values of the propagation constant, k, are compared in Fig. 1(c).
The power defect versus propagation constant k of 1D flatfloor bubbles, as obtained from the numerical results and as predicted by TFA, are shown in Fig. 1(d). Note that the TFA and numerical findings are virtually identical, with ∆P 1D growing linearly with |k|, in exact agreement with TFA. Furthermore, for g = 20 and g 0 = 1, which corresponds to Fig.  1(d), Eq. (30) predicts the slope d (∆P 1D ) /dk = 2 √ ln 20 ≈ 3.46, which is very close to the value 3.38 produced by Fig.  1(d).
The calculation of stability eigenvalues for the flat-floor bubbles has demonstrated that they are completely stable (at least, up to k = −24). This conclusion is corroborated by direct simulations of the modes marked by A1-A3 in Fig. 1(d), as shown in Figs. 2(a-c), respectively. Fig. 3 are profiles and the stability area of 1D dark solitons, constructed with the help of the input given by Eq. (28). Profiles of the solitons with different values of width x 0 of the flat section in the nonlinearity profile are shown in Fig. 3(a), where one can see that the ordinary dark solitons transform into flat-waist ones as x 0 increases. We present the profiles of the 1D dark solitons with different values of nonlinearity strength g in Fig. 3(b), in which the amplitude of the flat-waist segment decreases with the increase of g. Then, the profiles of dark solitons with different propagation constants k are shown in Fig. 3(c). Note that the variation of the darksoliton profiles following the change of x 0 , g, and k, observed in Fig. 3, are quite similar to those for the 1D bubbles shown in Fig. 1.

Depicted in
Properties of the dark-soliton families are summarized in Fig. 3(d) by means of lines showing the relation between the soliton's power defect and propagation constant. It is seen that the TFA, based on Eqs. (28), produces the ∆P (k) dependence for the dark solitons which is virtually identical to its numerically computed counterpart.
An essential result is reported in Fig. 3(e): the stability domain for the 1D dark solitons in the (x 0 ,k) plane, produced by the solution of the eigenvalue problem based on Eq. (22). It is seen that the decrease of the width of the flat-floor area, x 0 , or the decrease of |k| lead to destabilization of the dark soliton. In either case, it tends to lose its stability when it becomes relatively wide, in comparison with width 2x 0 of the central region. Further, Fig. 4 demonstrates that the destabilization, which occurs with the decrease of x 0 at a fixed value of k 0 , is accounted for by a simple bifurcation of the center-saddle type [88], with two imaginary eigenvalues colliding at the critical point and carrying over into a pair of real ones with opposite signs.
The stability of sufficiently narrow dark solitons can be understood following the principle that the stable configuration is one minimizing the system's Hamiltonian. To implement this approach, it is necessary to compare values of ∆H, given by Eq. (12), for the dark soliton placed in the central region, with G = g, and in the peripheral area (|x| > x 0 ), with G = 1, taking into regard the constraint that this should be done for equal values of the total power. Indeed, the values of the power defect, given by Eq. (11), are different for the same k but different local values of G(x), the defect being larger for smaller G = 1 than for G = g: (31) To maintain the balance of the total power, one may assume that the imbalance, δP , is distributed in a very broad flat background, with large spatial width L, so that the constant background amplitude, U backr = √ −k, for G = 1 (see Eq. (25)), is replaced by U backr = √ −k + δU , with a small change, which absorbs the imbalance of the total power. This shift, in turn, gives rise to a small change in the quartic term of Hamiltonian (4) in the background area: where expression (32) is substituted for δU . Finally, the total difference of values of the Hamiltonian for the dark soliton with the center placed at x = 0 and at |x| > x 0 is The positiveness of δH total in Eq. (34) suggests that the configuration with a relatively narrow dark soliton, placed in the broad top-floor region, with a larger value of G = g, is stable, in comparison with its counterpart in which the dark soliton is placed in the asymptotic area, where the local nonlinearity coefficient, G = 1, is smaller. Note a crucially important contribution of term (33), which represents the change of the Hamiltonian of the indefinitely broad background due to the small change of its amplitude: without this term, δH total would have a wrong sign.
On the other hand, for smaller |k| or smaller values of x 0 , when the dark soliton does not remain very narrow in comparison with size x 0 of the flat-floor region, Fig. 3(e) demonstrates that a relatively wide dark soliton, with the center placed at x = 0, is unstable. This fact may be qualitatively explained, considering the limit case of the very narrow modulation profile approximated by Eq. (16). Indeed, treating the dark soliton as a quasi-particle in the framework of the perturbation theory [89], the Hamiltonian perturbation term (17) acts on the particle as an effective potential, obtained by the substitution of the unperturbed dark soliton (10), with its center placed at x = ξ: By itself, this potential is attractive, with a minimum at ξ = 0.
However, the effective dynamical mass of the dark soliton, treated as a particle, is negative [81,82,90], therefore the action of potential (35) Fig. 2, it is spontaneously expelled from the flat-floor area, as might be expected.

C. Two-dimensional (2D) modes: bubbles and vortices
We now turn to the consideration of 2D flat-floor bubbles (which represent the ground state) and flat-waist vortices. At first, we dwell on results for the bubbles. The profiles and intensity distributions for them are displayed in Fig. 5. In particular, Fig. 5(a1) depicts the profiles of 2D flat-floor bubbles for different values of r 0 in the spatial-modulation pattern given by Eq. (5), whose intensity distributions are presented in Figs. 5(a2-a4). It is seen that the bubbles form the flat-floor shape as r 0 increases. Similarly, Fig. 5(b1) depicts the profiles of 2D flat-floor bubbles for different values of g, whose 2D intensity distributions are shown in Figs. 5(b2-b4). In Figs. 5(b1-b4) one sees that the minimum value of the bubble's field naturally decreases with the increase of g, in accordance with Eq. (26). Generally, the profiles of these 2D flat-floor bubbles are similar to their 1D counterparts, cf. Fig. 1. Further, accord- ing to the results of the linear stability analysis and numerical simulations, the flat-floor bubbles are all stable (as one may expect for the ground states, which are actually represented by the bubbles), at least up to |k| = 24. Dependencies of the power defect, ∆P 2D , on k for the 2D flat-floor bubbles, as produced by the numerical results and TFA (see Eq. (26)), are reported in Fig. 5(c1), and the intensity distributions with different k are depicted in Figs. 5(c2-c4). Note that the numerically found and TFA-predicted ∆P 2D (k) dependencies are virtually identical, both being proportional to |k|, as said above (see, in particular, Eq. (29)).
Next we present the results for flat-waist vortices. Typical examples of such modes with winding number m = 1 are displayed in Fig. 6 for r 0 = 2 and r 0 = 5. The power defect versus k for the 2D flat-waist vortices, as produced by the numerical solution and predicted by and TFA (see Eq. (26) for different values of r 0 , is displayed in Figs. 6(a1,b1). It is seen that, unlike the 2D bubbles (cf. Fig. 5(c1)), there is a small discrepancy between the numerical findings and TFA prediction, and the power defect is not strictly proportional to |k|. Typical intensity distributions, contour plots, crosssection profiles, and phase patterns of the vortices are presented in Figs. 6(a2-a5) and 6(b2-b5). These modes naturally feature flat waists instead of the flat floor, as the vorticity stipulates vanishing of the amplitude at r = 0. The width of the flat waist of the vortices increases with r 0 , similar to the 2D bubble (ground-state) modes presented in Fig. 5. The phase patterns are usual, directly corresponding to m = 1.
The results for double flat-waist vortices, with the winding number m = 2, are depicted in Fig. 7. The dependence of the power defect on k, intensity distributions, contour plots, and cross-section profiles of such double vortices are similar to those with m = 1, while the phase patterns in panels 7(c4) and (d4) obviously represent the double vorticity. Note a visible discrepancy between the numerical findings and TFA prediction for the power defect, as a function of k, in panel (a1), similar to what is shown for the 2D bubbles in Fig. 6(b1).
An essential issue is the stability of the vortices. According to our results, obtained from the linear-stability analysis and corroborated by direct simulations, the vortices with m = 1 are completely stable, at least up to k = −24. On the other hand, the double vortices are stable only if |k| is not too large. A stability chart for them in the plane of (r 0 , k) is shown in Fig. 7(b).
Simulations of the perturbed evolution of vortices from Figs. 6 and 7 are displayed in Fig. 8. First, panels (a) and (c) of this figure demonstrate the stable propagation of the unitary and double vortices labeled C2 in Fig. 6 and D1 in Fig.  7, respectively.
Further, Fig. 8(b) displays precession of the vortex whose pivot (center), with coordinates (x c , y c ) = (4.7, 0), is initially shifted off the origin. In Fig. 8(b), the pivot performs a circular motion in the periphery of the flat-waist area of the vortex. The trajectory of the motion is shown in Fig. 9(b). In addition, Fig. 9(a) shows the precession trajectory of the unitary vortex with the initial position (x c , y c ) = (1.5, 0). In the latter case, the radius of the precession motion decreases in the course of the evolution. A general conclusion is that, if the original vortex is stable, its precessing motion is robust as well. These results resemble experimentally observed precession of vortices in BEC [91].
Unstable evolution of the vortex labeled D2 in Fig. 7 is displayed in Figs. 8(d) and 9(c), where the double vortex (with m = 2) splits into a persistently rotating pair of unitary ones, with the rotation period z = 7.39. The character of the instability corresponds to the spectrum of eigenvalues for small perturbations around the double-vortex solution, which is displayed in Fig. 7(d5). It is seen, in this figure, that the instability is accounted for by a quartet of eigenvalues with nonzero imaginary parts, (with mutually independent pairs of ±), on the contrary to the pair of purely real eigenvalues that represent the instability of the 1D dark solitons in Figs. 4(a,b). The imaginary part of the unstable eigenvalues, i.e., Im(λ) in Eq. (36), determines the angular velocity of the rotation of the emerging pair of split unitary vortices in Figs. 8(d) and 9(c). The same mechanism, dominated by the quartet of critical eigenvalues, accounts for the instability of the double vortices in the entire yellow (unstable) area in Fig. 7(b). Actually, the emergence of a quartet of complex eigenvalues is known as a generic type of the destabilizing bifurcation. In particular, it typically accounts for the onset of instability of binary vortices trapped in an external potential [92] and discrete vortex solitons [93]. Actually, the splitting instability of double vortices (without setting the emerging pair in rotation) is a commonly known property of the 2D NLS equation with the spatially uniform cubic self-repulsion, see review [94] and original works [95- 101]. In the present case, the difference is that the spatial modulation of the local nonlinearity in Eq. (5) induces an effective trap, which prevents the separation of the split pair, keeping it in the rotating state. The same trap maintains the stability of the double vortices in the blue area of Fig. 7(b) (the stabilization of multiple vortices by an external potential is known in other models [102]).
Lastly, the precession of vortices with m = 2, which are initially placed off the center, is displayed in Figs. 9(d,e). In panels (d) and (e), the double vortices, with (r 0 = 5, k = −5) and (r 0 = 3, k = −2) are, respectively, initially unstable and stable against splitting, according to Fig. 7(b). In both cases, the shifted double vortices split in rotating pairs with a small separation between the unitary vortices, which perform robust precession. Further, the separation between the unitary vortices in the pair remains constant if the original double vortex is stable (Fig. 9(e)), while the size of the pair produced by splitting of an unstable double vortex slowly increases in the course of the evolution, as seen in Figs. 9(c,d). The consideration of vortices with m ≥ 3 is left beyond the scope of the present work.

IV. CONCLUSION
We have considered the model of the optical medium or BEC with intrinsic self-repulsion. The model includes the spatial modulation of the local nonlinearity strength, with a maximum in the central area and a minimum in the periphery. The analysis has revealed new modes, viz., the 1D and 2D flat-floor bubbles, 1D flat-waist dark solitons, and 2D flatwaist vortices. In a broad parameter region, the structure of the modes is accurately predicted by the TFA (Thomas-Fermi approximation). Through the computation of eigenvalues for small perturbations, it is predicted that the 1D and 2D flatfloor bubbles, as well as 2D vortices with winding number m = 1, are completely stable, which is corroborated by direct simulations. Nontrivial findings are instability boundaries, in their existence area, for dark solitons and vortices with m = 2. Unstable dark solitons spontaneously escape into the peripheral area with a smaller local self-repulsion coefficient, which is explained by means of the analytical approximation. Unstable double vortices split in rotating pairs of unitary vortices. Lastly, originally displaced stable vortices feature robust precession around the origin.
The consideration of vortices with m ≥ 3 may be interesting too, especially the precession in the presence of the flatwaist structure.