Reduced dynamics for one and two dark soliton stripes in the defocusing nonlinear Schr\"odinger equation: a variational approach

We study the dynamics and pairwise interactions of dark soliton stripes in the two-dimensional defocusing nonlinear Schr\"odinger equation. By employing a variational approach we reduce the dynamics for dark soliton stripes to a set of coupled one-dimensional"filament"equations of motion for the position and velocity of the stripe. The method yields good qualitative agreement with the numerical results as regards the transverse instability of the stripes. We propose a phenomenological amendment that also significantly improves the quantitative agreement of the method with the computations. Subsequently, the method is extended for a pair of symmetric dark soliton stripes that include the mutual interactions between the filaments. The reduced equations of motion are compared with a recently proposed adiabatic invariant method and its corresponding findings and are found to provide a more adequate representation of the original full dynamics for a wide range of cases encompassing perturbations with long and short wavelengths, and combinations thereof.


I. INTRODUCTION
In the past two decades, the study of coherent structures in the form of dark solitons has been a theme pervading a wide range of areas within Physics. Early examples were more focused on classical physics (including mechanical [1] and electrical [2] lattices), nonlinear optical [3], as well as magnetic systems [4]. More recently, however, there has been a host of additional systems including notably a wide variety of experiments in atomic Bose-Einstein condensates summarized, e.g., in Refs. [5,6], but also realizations in electromagnetic [7], hydrodynamic [8], acoustic [9], plasma [10] and excitonpolariton [11] systems among others.
A major thrust within the subject has been offered by the extensive experimental accessibility to such nonlinear waves offered in the context of atomic (but also exciton-polariton) condensates. Here, there has been a variety of techniques of formation of the structures, including wave interference [12,13], phase imprinting [14] and rapid dragging of laser beams through a trapped BEC [15]. Additionally, the structures have offered a variety of intriguing insights in the dynamics through their potential dynamical instabilities (and their avoidance through suitable manipulation of length scales [16]) which may lead to a variety of vortical (in 2D) and vortex line/ring structures (in 3D), as has also been documented experimentally [17][18][19][20]. More recently extensions to multi-component settings have been pursued in the form of dark-bright and dark-dark solitonic states [21] and even spinorial realizations of such structures have recently been identified [22]. At the same time, there has * URL: http://nlds.sdsu.edu been a growing interest in applications including possibilities of atomic matter-wave interferometers [23] and proposals for the use of dark solitons as qubits in BECs [24].
The study of transverse ("snaking") dynamical instabilities of dark solitons in higher dimensional settings has been a topic of extensive interest since early on [25], with much of the early work summarized in the review [26]. Recently, there has been a surge of further interest in the subject [27][28][29] fueled by an adiabatic invariant (AI) based insight enabling the derivation of effective equations for the dark soliton stripes (in 2D) and planes (in 3D), but also the ability of this methodology to tackle ring (in 2D, but also in 3D in the form of spherical) solitons [30]. This methodology was not only seen to have the right long wavelength limit (a fundamental prerequisite for such a theory). Additionally in many cases, including those of ring and spherical solitons, it provided with unprecedented accuracy and simplicity analytical approximations for the frequencies of vibration and destabilization of the coherent structures that were tested in both linearization (spectral) computations, as well as in the fully nonlinear dynamics of the system.
Nevertheless, there is reason to believe that the theory can be further improved. On the one hand, while the above AI methodology captures the correct long wavelength limit, it does not a priori capture the restabilization of perturbations of large wavenumbers (above a certain k c ). Moreover, as it stands, the theory is developed solely for the evolutionary dynamics of the center of the dark solitonic stripes (or planes), but does not arise as a coupled theory for the evolution of the center and the width (or velocity) of the structure. For these reasons, but also in order to obtain a theory with a definitive Lagrangian framework (avoiding the issue of identification of canonical variables and) enabling the systematic derivation of Euler-Lagrange equations of motion, we develop herein an alternative variational approach.
At the heart of our analysis lies a two-dimensional generalization of the classic work of Ref. [31] considering the variational characterization of the dynamics of a dark soliton in 1D. Here, we endow the dark soliton stripe (DSS) with a center, a width and a speed that are transversely dependent (as in the recent AI work of Refs. [27][28][29]), yet we substitute the relevant ansatz of one or two stripes in the Lagrangian of the model. The subsequent derivation of the Euler-Lagrange equations of motion for two dynamical variables (e.g. position and velocity) constitutes the basis for our further analysis of the theory and its improved, as we will show, agreement with the full numerical results. We discuss some of the limitations of the method, such as its qualitative but not quantitative capture of the restabilization of large wavenumbers and offer relevant amendments that provide the optimal characterization, to our knowledge, of the motion of single and multiple dark solitonic stripes available to date. The manuscript is structured as follows. Section II describes the VA methodology and puts forwards the reduced equations of motion for a single DSS and for two (symmetrically displaced) interacting DSSs. Section IV is devoted to testing the validity of the VA approach and to compare its predictions against a previous filament technique based on adiabatic invariance (AI) [27][28][29]. Finally, in Sec. V we summarize our results and give possible avenues for future research.

II. VARIATIONAL APPROACH FOR ONE STRIPE
Our starting point will be the prototypical 2D NLS model [3,6]: where u 0 is the magnitude of the background and α = ∓1 accounts for the elliptic and hyperbolic NLS cases. In our numerical results we will focus on the elliptic NLS case (i.e., α = −1). Nonetheless, for genericity's sake, we keep α in our analysis herein. In the 1D case (α = 0), the NLS admits the following exact traveling dark soliton solution: where We recall that k accounts for the velocity of the background while v denotes the velocity of the dark soliton itself. To variationally follow the dark soliton as a quasi-1D stripe (or filament) in the 2D NLS (1), we consider the 2D extension u 2D = u 2D (x, y, t) by allowing the position X(y, t) of the dark soliton to be a function of y and t (in the spirit also of earlier works such as Refs. [27][28][29]), while keeping a stationary background (k = 0), as follows: where we consider the inverse width of the DSS D to be, in general, decoupled (in the dark soliton core) from the background level B. Nonetheless, we enforce A 2 + B 2 = u 2 0 =constant to keep the DSS from affecting the tails at x = ±∞. Note that when B = D and v = X t = A, the DSS ansatz (3) reduces to the exact dark soliton (2) mounted on a stationary background (k = 0).
The NLS (1) is derived from the renormalized (in the x-direction) Lagrangian: where the averaged (i.e., integrated) Lagrangian along the x-direction may be written as We note that the term 1 − u 2 0 /|u| 2 is introduced to renormalize the momentum while the term |u| 2 − u 2 0 renormalizes the power. This renormalization in introduced so that the (1D) averaged Lagrangian L y converges when evaluated over the 1D dark soliton (2) [30,31].
Let us now evaluate the 2D Lagrangian (4) over the ansatz (3) by assuming that the dark soliton moves locally and thus it does not affect the tails (background) in A and B. Namely, we will use the following approximation for the y-derivative of the 2D ansatz (3): Alternatively, the above approximation can be thought of assuming a generic traveling profile u ≈ f [D(t, y)(x − X(t, y))] that yields, using the chain rule, u y = u x Dy D (x − X) − X y . Finally, after evaluating the expression (6) we enforce that, due to the balance outside the dark soliton core region, D = B and A 2 + B 2 = u 2 0 which leads to the simplified average Lagrangian: As per the Euler-Lagrange prescription, we now take variations in the variables X and B which yield, respectively: These coupled partial differential equations (PDEs) represent the equations of motion for the 2D ansatz (3) in terms of the variables A, B, and X. Using the relation A 2 + B 2 = u 2 0 , we can decouple and rewrite these coupled PDEs in terms of the variables A and X: where B 2 = u 2 0 − A 2 . Let us now study the (linear) stability for the reduced VA PDE system of Eqs. (10) and (11). To this end, we linearize the system around the equilibrium configuration A = η and X = 0 by considering A = η+a and X = ηt+ξ for a, ξ ≪ 1 to find: Linear stability of plane waves [a(y, t), ξ(y, t)] = a 0 e i(ky−ωt) , b 0 e i(ky−ωt) the dispersion relation: Depicted is the eigenfrequency ω as a function of the transverse wavenumber k for µ = 1. Same notation and layout as in Fig. 1 where now the thick green dashed line corresponds to the improved VA prediction of Eq. (22). For completeness we also show the real part of the eigenfrequency for the full NLS system (thin blue dashed curves) and the improved VA (green dotted line).
As depicted in Fig. 1 (see green dashed curves), Eq. (14) describes a dispersion relation ω = ω(k) that has qualitatively the correct shape when compared to the stability of the full (original) NLS model (1) (see solid curves), obtained numerically. However, we note that, despite having the correct trend, the dispersion curves fail to accurately capture quantitatively the actual growth rates and wavenumber cutoffs for stability (i.e., wavenumbers k c such that ω(k c ) = 0). This qualitative shape improves on the previous reduced description for the dynamics of DSSs based on the AI assumption [27][28][29]. The reduced AI PDE for DSSs yields, for η = 0, the following (linear) dispersion relation [29]: which simply predicts a linear growth (see thin black line in Fig. 1) of the instability rates as the wavenumber increases. Therefore, inspired by the correct qualitative shape of the dispersion curves from the VA reduced model, we now propose a modification of the VA model so as to also quantitatively match more adequately the dispersion curves. In order to modify the VA approach to appropriately capture the correct instability growth rates, we choose to modify the reduced VA PDEs so to match the correct wavenumber cutoffs k c . We thus amend our Lagrangian (7) by considering the phenomenological variant: where we have introduced the function so as to match the wavenumber cutoffs k c obtained from an asymptotic approximation to the linear stability problem [32]. With this modification, the improved equations of motion (for the elliptic case α = −1) for a DSS yields: The corresponding linearization for this new reduced PDE model yields: which in turn yields the linear dispersion By construction, Eq. (22) captures, as expected, the critical wave numbers k 2 c = h( u 2 0 − η 2 ) where instability for k < k c changes to stability for k > k c for all values of the propagation velocity A = η, see Fig. 2. Furthermore, as it can be observed from Fig. 2, not only the critical cutoffs k c match, but the maximum growth rates are also well approximated with the improved VA equations. Therefore, we expect that the modified Lagrangian (16), leading to the VA equations (18)- (19), is able to give a good description for the DSS dynamics for all wavenumbers. We will return to this point in Sec. IV, where we will present our numerical comparisons. This is in contrast with the AI approach [27][28][29] which, by construction, is valid for small wavenumbers and should thus be expected to be more adequate there (as will be again seen in Sec. IV).

III. VARIATIONAL APPROACH FOR TWO STRIPES
In this section we develop a VA approach to follow the interaction of two DSSs akin to what was obtained using the AI approach in Ref. [29]. Consider again the defocusing (α = −1) 2D NLS Eq.(1), which, as mentioned before, may be derived from the renormalized (in the x-direction) Lagrangian (4)- (5), which we now split as follows: where In order to variationally follow the 1D two-dark soliton extension of Eq. (3) as two DSSs in 2D, we consider the following 2D symmetric extension by allowing each dark soliton to have its own position and keep a stationary background (k = 0): where we use the following short-hand definitions: where X 1 (y, t) and X 2 (y, t) are the spatio-temporal locations of the two DSSs with X 1 < X 2 . The two-DSS ansatz (27) has an overall background level B and velocity A (A > 0 representing the two DSSs moving outwards and A < 0 inwards) where, as for the one DSS case, the relation A 2 + B 2 = u 2 0 =constant remains valid. It is important to mention here that a more elaborate two DSS ansatz with independent velocities for each dark soliton does not allow for a tractable VA approach (i.e., integrals that can be obtained in closed form in the Lagrangian). Therefore, we restrict our attention to the ansatz (27) that has the drawback of having a common (symmetric) velocity for both dark solitons so that a solution that is not symmetric, namelyẊ 2 = −Ẋ 1 will eventually tend to a symmetric one (Ẋ 2 = −Ẋ 1 ) as time evolves (see Sec. IV for more details).
Let us now perform the VA approach using the two DSS ansatz (27). Letting ∆ ≡ X 2 − X 1 > 0 (i.e., the two dark solitons do not overlap, nor change relative positions) and assuming B∆ ≫ 1 (i.e., the two dark solitons are relatively well separated), we have the following useful approximations: where the precise form of the functions L, F , G and H is not particularly important at this stage as they all contribute to the same order (e −2B∆ ) and will all be combined appropriately in the final, explicit results below [cf. Eq. (38) in what follows]. Using this information, we compute: Using the above we can simplify the integrand of L 1 as: where the approximation, has been used with z 1t = Dt D z 1 − DX 1t and z 2t = Dt D z 2 − DX 2t . As it was the case for the functions L, F , G, and H above, the precise form of the function R is not relevant at this stage as it will be combined in the explicit Lagrangian given in Eq. (38) below. We thus get, upon integration, where On the other hand for the L 2 integral, using exact integrations yields: where g (z) ≡ −16e 2z (e 2z −1) 5 1 + (9 + 12z) e 2z −(9 − 12z) e 4z −e 6z .
Finally, considering the transverse dependence of u, from Eq. (28), we rewrite (27) 1 (y, t)) + 1]+iA. As in the single dark stripe case, assuming that B and A do not depend on y, directly but through the relations D = B and A 2 + B 2 = u 2 0 for D = D(y, t), yields: Assuming small and slow displacements, we can safely neglect the cross terms D y X 1y and D y X 2y from z 2 1y , z 2 2y and z 1y z 2y and thus we obtain: where Now, recalling the introduction of the factor h(B) in Eq. (17) to improve the agreement of the growth rates with the numerically observed ones, we finally obtain the effective Lagrangian: According to the Euler-Lagrange prescription, taking variations over X 1 , X 2 , and B and recalling that ∆ ≡ X 2 − X 1 and A 2 + B 2 = u 2 0 , yields: Equations (39)-(41) represent one of the main results of this work where the VA methodology has been employed to reduce the full dynamics of two interacting DSSs in two spatial dimensions to these three (1+1)D coupled PDEs. It is interesting to note that, if one starts with a symmetric initial configuration X 2 (y, t = 0) = −X 1 (y, t = 0), given the symmetry of Eqs. (40) and (41), the dynamics preserves this symmetry at all times [i.e., X 2 (y, t) = −X 1 (y, t) is an invariant manifold of the dynamics during the evolution] and hence the coupled PDEs can be reduced to solely the equations Eq. (40) for X 1 and Eq. (39) for A.

IV. NUMERICAL RESULTS
In this section we corroborate that the dynamical reduction, for a single DSS and for two interacting DSSs, obtained through the VA approach is indeed valid under a wide range of initial conditions. Moreover, we compare the VA results with the AI methodology put forward in Refs. [27][28][29] and showcase where the two display similar results, as well as where the former represents a significant improvement over the latter. It is important to stress that, from now on, we use the improved VA model that includes the h(B) term introduced in Eq. (17). Let us start by considering a single DSS. The DSS is always unstable (snaking instability) to transverse wavenumbers k such that 0 ≤ k ≤ k c (see previous section). Therefore, if one considers a spatial domain (x, y) ∈ [−L x , L x ] × [−L y , L y ], with a DSS aligned in the y-direction (i.e., using the same notation as in the previous section), there will always be snaking instability provided that L y > π/k c . Conversely, if the domain is too small, unstable transverse modes do not "fit" inside the domain and thus the DSS is rendered stable (recall for instance how this property can be used to arrest the instability of DSSs in trapped BECs in Ref. [16]).
In order to test the validity of the improved VA approach, we performed a series of simulations for the dynamical evolution of a single DSS under different perturbations using second order finite differencing in space [33] with fourth-order Runge-Kutta with periodic boundary condition along the y-direction and mod-squared Dirichlet (MSD) boundary conditions [34] along the x-direction in order to avoid any undesired effects from the boundaries. The simulations depicted in Figs. 3-5 are aimed at controllably testing the dynamics of perturbations with different wavenumbers in the domain L x = L y = 20. In particular, Fig. 3 depicts the dynamics ensuing from an initially stationary [η = A(y, t = 0) = 0] DSS that has been perturbed with the longest possible wavelength (satisfying the periodic boundary conditions in the ydirection). In addition to testing the VA reduced equations of motion, we also implement the AI methodology put forward in Refs. [27][28][29]. As it can be seen from the figure, both the improved VA [see the thick green (gray) curves] and the AI [see the dark blue (black) dotted curves] methodologies give a reliable description of the snaking dynamics up to the point where the DSS loses  its transverse dark-soliton-like profile as it nucleates vortices of alternate signs at the nodes of the perturbation mode. A few remarks are in order at this point. Firstly, neither the VA nor the AI methods were designed, by construction, to follow the stripe after losing its darksoliton-like stripe shape. Thus, as DSSs tend to decay into vortex patterns, there will always be a point in time where the VA (or AI) basic assumptions will be violated (most notably the ansatz of a dark solitonic stripe being an accurate descriptor of the full 2D field) and thus the  Fig. 3 but for a perturbation of the initial position of the dark soliton given by X(y, t = 0) = 10 j=1 εj sin[kπ(y + ϕj )/Ly] with εj = 0.01, ϕj = (j − 1)Lyπ/5, and k = j. See Supplemental Material movie-snake-1-10 for an animation depicting the corresponding dynamics [35]. dynamical reduction will no longer be valid. On the other hand, as it can be noticed from Fig. 2, for low wavenumbers pertaining the case of Fig. 3, both VA and AI are able to appropriately predict the correct growth rate of perturbations. However, as one notices from Fig. 2, both VA and AI tend to slightly overestimate the growth rates. This is precisely what it is observed from Fig. Fig. 3 where both VA and AI reductions tend to slightly "run faster" in the dynamical destabilization. Note however, that the VA's overestimation is slightly smaller than the one obtained through the AI. As we see next, this issue with the AI overestimation will become more acute for larger perturbation wavenumbers. Figure 4 depicts a similar case as the one depicted in Fig. 3 but for the mode with the second largest wavelength. Again, both the improved VA and AI tend to slightly overestimate the growth rate of the perturbation with the VA approximation being better than the AI one. In order to test a more complex scenario where the essence of the improvement of the theory proposed in this work is most dramatically evident, we perturbed the stationary DSS with a combination of the first 10 modes. The results are depicted in Fig. 5. As it can be observed from the figure, the VA reduction does an excellent job at following the dynamical stabilization of the DSS. On the other hand, as we are now introducing larger wavenumbers, the AI is clearly less accurate (as we expected from its spectral predictions) and tends to considerably overestimate the growth rates. In fact, it can be noticed that the VA approximation is even able to track the full nonlinear dynamics of the DSS up to the point (and even, arguably, slightly beyond, considering e.g. the snapshot at t = 18.5) where it has broken into a pattern including vortex-antivortex pairs. Furthermore, as the configuration contains larger wavenumbers, the instability sets in faster and, thus, the slight growth rate overestimation provided by the VA is now downplayed over the span of time before the DSS breaks into vortices. This highlights that in a typical scenario where several unstable modes are present, the improved VA will be an excellent reduced description of the full dynamics of single DSSs.
We now proceed to test the improved VA prediction for the interaction of two DSSs. Similar to the one DSS case, let us probe both the long wavelength scenario and multiple, mixed mode case. Figure 6 depicts the evolution for two stationary (side-by-side) DSSs symmetrically perturbed by the longest possible wavelength in the provided domain. As the picture shows, the VA is able to track extremely well the DSS dynamics up to the point where the DSSs touch, reconnect, and split into patterns involving vortical structures. It is important to note that we have chosen a configuration where the interactions between the two stripes are quite not trivial. This can be noticed by the fact that the top portions of the DSSs (for y ≈ 8) initially get closer to each other (due to the individual growth rates of the snaking instability for each DSS) and then repel each other (t > 20) once the DSS proximity is such that the mutual dark soliton repulsion dominates the dynamics. A more compelling case can be made by testing the VA prediction for the two DSSs by starting with an initial condition containing a non-trivial combination of the first ten modes; see Fig. 7. As for the one DSS case presented in Fig. 5, the two DSS VA reduction is also able to adequately track the two DSSs even well after they break up into individual vortices; see in particular the snapshots between t = 18 and t = 26. It is relevant to mention that, as in single DSS case, two DSSs case may also be approximated by the AI methodology [29]. However, one needs to keep in mind that the AI approach is valid for cases containing long wavelengths. Thus, for cases containing shorter wavelengths (larger wavenumbers), as it is particularly the case depicted in Fig. 7, the VA is a very good approximation while the AI will fail in a similar manner akin to the results presented in Fig. 5 for the single DSS.
Finally, for completeness, we briefly study the case where the perturbation on each of the two DSSs is not symmetric. In this case, due to the chosen ansatz (containing a shared velocity term between the two DSSs), the VA dynamics necessary leads, by construction, to a symmetric configuration. Therefore, in principle, one would not expect the VA to give a meaningful prediction for the two DSSs when perturbed asymmetrically. Nonetheless, we tested that for perturbations with small asymmetries the VA does indeed reasonably well at describing the evolution of the two DSSs, despite its obvious shortcoming in that the evolution of the full dynamics is no longer symmetric. As an example, we depict in Fig. 8 a couple of cases where the left DSS is left unperturbed  where the right soliton position is perturbed by a linear combination of the first ten modes with (weak) strengths of 0.01 (top series of panels) and 0.001 (bottom series of panels). As the figure shows, after some transient, the original NLS (1) dynamics evolves such that the left DSS develops undulations (exerted by the right DSS) that are in fact approximately symmetric with respect to the undulations of the right DSS. Therefore, it is not completely surprising that this behavior is, to some extent, followed by the VA as the results in Fig. 8 show.

V. CONCLUSIONS
In the present work we deployed the VA methodology in order to provide an improved description of the transversely unstable dynamics of single and multiple (two symmetrically perturbed) DSSs in the defocusing two-dimensional nonlinear Schrödinger equation. The method consists of applying the VA to suitable ansätze that consist of one or two DSSs whose transverse position and velocity are functions of time and the transverse spatial variable. In this manner, we are able to reduce the original (2+1)D NLS into two coupled PDEs in (1+1)D for the dark soliton's position and velocity. It is also important to highlight here that the standard VA is qualitatively adequate, yet quantitatively misses FIG. 8: (Color online) Snaking dynamics of two DSSs for µ = 1 perturbed by a non-symmetric perturbation where only the right DSS is perturbed. Similar as in Fig. 7 but for an unperturbed left DSS given by X1(y, t = 0) = −x0 and a right DSS perturbed by X2(y, t = 0) = x0 + 10 j=1 εj sin[kπ(y + ϕj)/Ly] with x0 = 2, ϕj = (j − 1)Lyπ/5, and k = j (i.e., same as in Fig. 7), where the perturbation strengths for the right DSS are: εj = 0.01 (top series of panels) and εj = 0.001 (bottom series of panels). Note that although the dynamics for both cases is very similar, the case with the weakest perturbation (see bottom panels) displays a perturbation growth at a slower time scale (contrast the difference in times between the two cases). Importantly also note the deviation induced by the asymmetry between the (by construction symmetric) VA and the full dynamics. See Supplemental Material movie-2snakes-0-0-1-10-A and movie-2snakes-0-0-1-10-B for animations depicting the corresponding dynamics.
the spectral instability features of the single DSS. In view of that, we have proposed an amended variant of the VA which captures the critical wavenumbers k c of transverse restabilization (for different speeds) and consequently performs quantitatively better over the entire range of wavenumbers. Our numerical results indicate that for short wavenumbers the VA does a slightly better job at predicting the growth rates of perturbations compared to the adiabatic invariant technique developed in Refs. [27][28][29]. However, it is for perturbations containing larger wavenumbers (including a non-trivial mix for short and long wavenumbers) that the advantage of the VA methodology is more apparent. Specifically, the reduced equations of motion obtained through the VA are able to closely track the dynamics for one and two interacting DSSs. In fact, the reduced VA methodology is not only successful in tracking the DSSs dynamics up to the point where the DSSs break up into vortices; surprisingly, it continues to adequately trace the position of the remnants of the former stripe even for times slightly beyond its destabilization and breakup into vortical structures. Our numerical results suggest that the reduced VA equations constitute a viable methodology for theoretically reducing [to a (1+1)D setting] and numerically closely tracking the dynamics of DSSs over a wide range of dynamical scenarios (i.e., containing a wide range of perturbations from short to long wavelengths ones).
It is worth mentioning that the VA methodology seems not to be tractable for a general two dark soliton ansatz. Therefore, in our presentation we focus on a more specific ansatz where the two dark solitons share velocities and are thus pushing to stay symmetric (each one being the mirror image of the other one). It would be interesting to generalize the results presented herein by a suitable, and tractable, VA (or other) methodology that would be able to successfully track two DSSs for arbitrary positions and velocities. This would allow to not only treat the general two DSSs case, but the general case of N interacting stripes. Of course, there are other directions that are relevant to pursue as well, based on the program that has been recently developed at the level of the adiabatic invariant methodology. Some of them concern ring dark soliton structures, multi-component (e.g. darkbright) solitonic patterns, as well as three-dimensional states such as planar or spherical solitons. From a theoretical perspective, understanding better how to justify an amendment like the one proposed herein from first principles (that significantly improves the quantitative tracking of the transversely unstable modes) is also an important challenge for future studies.