Dynamics and stabilization of bright soliton stripes in the hyperbolic-dispersion nonlinear Schr\"odinger equation

TL;DR

This paper investigates the dynamics and stability of bright soliton stripes in a two-dimensional nonlinear Schrödinger equation with hyperbolic dispersion, analyzing instabilities and stabilization methods.

Contribution

It introduces a variational approximation for the transverse dynamics and demonstrates stabilization of soliton stripes using a channel-shaped potential.

Findings

The variational approximation effectively models transverse instabilities.

Channel-shaped potentials can fully stabilize bright soliton stripes.

The study characterizes snaking and necking instabilities in the system.

Abstract

We consider the dynamics and stability of bright soliton stripes in the two-dimensional nonlinear Schr\"odinger equation with hyperbolic dispersion, under the action of transverse perturbations. We start by discussing a recently proposed adiabatic-invariant approximation for transverse instabilities and its limitations in the bright soliton case. We then focus on a variational approximation used to reduce the dynamics of the bright-soliton stripe to effective equations of motion for its transverse shift. The reduction allows us to address the stripe's snaking instability, which is inherently present in the system, and follow the ensuing spatiotemporal undulation dynamics. Further, introducing a channel-shaped potential, we show that the instabilities (not only flexural, but also those of the necking type) can be attenuated, up to the point of complete stabilization of the soliton stripe.

Equations

$$i{u}_t=-\frac{1}{2}{u}_{xx}+\frac{1}{2}{u}_{yy}-|u{|}^{2}u+V\left(x,y\right)u.$$

$$u\left(x,y,t\right)=A\mathrm{sech}\left[A\left(x-x_0-vt\right)\right]e^{i\left(vx+\mu t+{\theta }_0\right)},$$

$$\mu =\frac{A^2-v^2}{2}.$$

$$u(x,y,t)$$

$$H_{2D}=\frac{1}{2}{\iint }_{-\infty }^{\infty }\left[|u_x{|}^2-|u_y{|}^2-|u{|}^4\right]dxdy.$$

$$u(x,y,t)=A\mathrm{sech}\left[A\left(x-X(y,t)\right)\right]\times$$

$$\exp \left\{i\left[v(x-X(y,t))+\frac{1}{2}\left(A^2+v^2)\right)t\right]\right\}.$$

$$H_{2D}=\frac{1}{2}\iint \left[(1-X_y^2(y,t))|u_x{|}^2-|u{|}^4\right]dxdy.$$

$$H_{2D}=\frac{1}{2}{\int }_{-\infty }^{+\infty }\left[(1-X_y^2)\left(\frac{1}{3}{A}^3+A{X}_t^2\right)-\frac{2}{3}{A}^3\right]dy.$$

$${\int }_{-\infty }^{+\infty }{\left[(1-X_y^2)\left(\frac{1}{3}{A}^3+A{X}_t^2\right)-\frac{2}{3}{A}^3\right]}_{t}dy=0.$$

$${\int }_{-\infty }^{+\infty }[X_{yy}\left(\frac{A^3}{3}+A{X}_t^2\right)$$

$$+2X_y{X}_{yt}{X}_t-(1-X_y^2)X_{tt}]X_{t}dy=0,$$

$$(1-X_y^2)X_{tt}=X_{yy}\left(\frac{A^3}{3}+A{X}_t^2\right)+2X_y{X}_{yt}{X}_t.$$

$$X_{tt}=\frac{A^3}{3}{X}_{yy},$$

$${\omega }^2+\frac{A^2}{3}{k}^2=0.$$

$$L$$

$$u\left(x,y,t\right)=A\left(y,t\right)\mathrm{sech}\left[A\left(y,t\right)\left(x-\xi \left(y,t\right)\right)\right]\times$$

$$\exp \left\{i\left[v(y,t)(x-\xi (y,t))+\theta (y,t)\right]\right\},$$

$$L={\int }_{-\infty }^{+\infty }\mathcal{L}dy,$$

$$\mathcal{L}$$

$$V^{\mathrm{eff}}(A,\xi )={\int }_{-\infty }^{+\infty }V(x,y)h(A,x-\xi )dx,$$

$$\delta \theta :A_t-{\left[A\left({\theta }_y-v{\xi }_y\right)\right]}_y=0,$$

$$E=\iint |u(x,y){|}^{2}dxdy,$$

$${\partial }_t\iint |u(x,y,t){|}^{2}dxdy=2A\left({\theta }_y-v{\xi }_y\right){|}_{y=-\infty }^{y=+\infty },$$

$$\delta v:{\xi }_t=v-{\xi }_y\left(v{\xi }_y-{\theta }_y\right)+\frac{{\pi }^2}{12A}{\left(\frac{v_y}{A}\right)}_y,$$

$$\delta \xi :v_t=\frac{-1}{3A}{\left(A^3{\xi }_y\right)}_y-v_y\left(v{\xi }_y-{\theta }_y\right)-\frac{1}{2A}{V}_{\xi }^{\mathrm{eff}}(A,\xi ),$$

$$-\frac{12+{\pi }^2}{36}\left(\frac{A_{yy}}{A}-\frac{A_y^2}{2A^2}\right)={\theta }_t-v{\xi }_t-\frac{A^2}{2}\left({\xi }_y^2+1\right)+\frac{v^2}{2}$$

$$+\frac{{\pi }^2}{24}\frac{v_y^2}{A^2}-\frac{1}{2}{\left(v{\xi }_y-{\theta }_y\right)}^2+\frac{1}{2}{V}_A^{\mathrm{eff}}(A,\xi ),$$

$${\theta }_t=\frac{1}{2}{v}^2-\frac{1}{2}\left(v^2{\xi }_y^2-{\theta }_y^2\right)+\frac{{\pi }^{2}v}{12A}{\left(\frac{v_y}{A}\right)}_y-\frac{{\pi }^2}{24}\frac{v_y^2}{A^2}$$

$$-\frac{12+{\pi }^2}{36}\left(\frac{A_{yy}}{A}-\frac{A_y^2}{2A^2}\right)+\frac{A^2}{2}\left({\xi }_y^2+1\right)-\frac{1}{2}{V}_A^{\mathrm{eff}}(A,\xi ).$$

$$A_t$$

$${\xi }_t$$

$$v_t$$

$${\theta }_t$$

Full text

I Introduction

Nonlinear Schrödinger (NLS) equations are a class of universal models governing the nonlinear propagation of waves in dispersive and diffractive media sulem ; ablowitz . In particular, in one dimension (1D) the interplay of the self-focusing cubic nonlinearity and paraxial diffraction, or, alternatively, anomalous group-velocity dispersion (GVD) gives rise to the commonly known bright solitons Zakharov , which are replaced by dark solitons in the case of the normal GVD Agrawal ; siam . The NLS equations play an equally important role in 2D and 3D settings, where the diffraction and GVD act together, in the combination with the cubic self-focusing or defocusing review ; sulem ; JiankeBook:2010 ; Fibich ; siam . In the case of the anomalous sign of the GVD and self-focusing sign of the nonlinearity, 2D and 3D solitons (in particular, “light bullets” in optics bullets ) exist too as formal solutions to the multidimensional NLS equations, but they are strongly unstable.

An alternative possibility is to consider the interplay of the paraxial diffraction in one transverse direction (assuming that the setting is two-dimensional) with the normal GVD and self-focusing, which is usually called, as a whole, the NLS equation with the hyperbolic 2D dispersion. The latter is a model of wide interest sulem ; CoDiTri ; CoTri ; GHS1 ; GHS2 ; GHS3 , appearing in a number of different applications. These include but are not limited to deep water waves abseg ; zakh , cyclotron waves in plasmas sen ; myra and nonlinear optics trillo1 ; trillo2 , to name a few. The hyperbolicity is realized as the opposite signs of the diffraction and GVD terms in the equation, see Eq. (1) below. In this case, 2D bright solitons do not exist, but one can consider quasi-1D solitons, which are self-trapped in the transverse direction, and uniformly extended along the longitudinal coordinate. Because such bright-soliton stripes may have obviously important physical realizations, a significant problem concerns the stability of the stripes Kuznetsov:1986 . It is well known that they are actually subject to the snaking (flexural) instability.

An initial purpose of the present work is to extend the known analytical approach to the description of the long-wave snaking instability, so as to explore its manifestations for perturbations with finite wavelengths, as well as the development of the instability beyond the usual limit of the linear approximation. The main objective of the work is to elaborate a possibility to attenuate and, eventually, completely suppress the instability by trapping the stripe in an underlying channel potential. This is somewhat in the spirit of our earlier effort to stabilize dark solitons via a potential barrier in the self-defocusing NLS case dark1 , although the stabilization is trickier here, as we will explain in detail below, due to the accessibility of the solitary waves to different types of instabilities.

The rest of the paper is organized as follows. The model, its basic features and theoretical setup are given in Sec. II. Effective longitudinal equations of motion for the soliton stripes are derived in Sec. III, using the method based on the use of the adiabatic invariant, recently quite successfully used in the case of dark solitons. Subsequently, a more elaborate (and more accurate in the present setting) variational approximation is developed. In Sec. IV we numerically investigate the dynamics and stability, as produced by different effective stripe-evolution equations, and compare the results with those directly produced by the underlying NLS equation. In Sec. V we introduce an appropriately crafted external potential that is able to control (and eventually eliminate) all potential instabilities. Finally, in Sec. VI we give a brief recap of the results and outline directions for further work.

II Model and Theoretical Setup

As per the above discussion, our model of choice will be the 2D NLS equation with the hyperbolic dispersion in the presence of the external potential, $V(x,y)$:

[TABLE]

Here $u(x,y,t)$ is a complex-valued field, and, as indicated above, the dispersions along the $x$- and $y$-directions have opposite signs. This setting corresponds, for the fixed sign of the cubic nonlinearity, to self-focusing and defocusing in the former and latter directions, respectively. In terms of the realization in a planar optical waveguide in the particular context of nonlinear optics, the evolution variable $t$ is actually the propagation distance, while $x$ is the transverse coordinate, and $y$ is the reduced temporal coordinate Agrawal . In this case, terms $u_{xx}$ and $u_{yy}$, with the signs adopted in Eq. (1), represent, respectively, the paraxial diffraction and normal GVD, while the cubic term stands for the usual Kerr nonlinearity. In what follows below, we consider the settings without the external potential, $V(x,y)=0$, and with a 1D potential, $V(x)$, that will help to “guide”, and potentially stabilize the soliton stripes.

In the absence of the potential, Eq. (1) admits a solution in the form of quasi-1D bright-soliton stripes,

[TABLE]

with

[TABLE]

This is a straightforward 2D homogeneous extension of the commonly known 1D bright soliton of the focusing NLS equation. It contains four free parameters which determine its height (and inverse width) $A$, velocity $v$, initial position $x_{0}$, and initial phase $\theta_{0}$.

The study of the stability of bright soliton stripes has a time-honored history, including some prevalent misconceptions about the type and origin of some of its ensuing instabilities (see for instance Ref. JiankeBook:2010 for a comprehensive review of the subject). As bright soliton solutions are stable in 1D, it is natural to study the stability once they are uniformly extended in 2D along the $y$ direction. In the elliptic dispersion case [with a minus sign in front of $u_{yy}$ in Eq. (1), which corresponds to the anomalous GVD in optics Agrawal ], the NLS is prone to collapse if the (squared) $L^{2}$ norm (alias total energy, in terms of optics) of a localized mode exceeds the critical value, which coincides with the corresponding norm of the Townes solitons sulem ; Fibich . The norm of soliton stripes, as they extend to $y=\pm\infty$ is infinite, hence the stripes are always prone to the modulational instability, which is related to the collapse-onset trend and is called necking instability Kuznetsov:1986 . It eventually leads the soliton stripe to breakup into individual lumps that are subsequently led to collapse (individually, or after their merger with other such lumps) review .

On the other hand, in the framework of Eq. (1), which includes the hyperbolic dispersion, the $y$-direction is effectively (in tandem with the nonlinearity) a defocusing one, and, thus, it should provide a stabilizing effect against necking. To some degree, the necking is indeed reduced due to the defocusing nature of the sign combination of the $y$-dispersion (normal GVD in optics) and nonlinearity. However, as it was detailed in Ref. Bernard-Peli:2006 and is explained below, the bright soliton stripe in the hyperbolic case still suffers from (relatively weak; see below) necking instabilities. On the other hand switching from elliptic to hyperbolic dispersion introduces strong snaking transverse instability on the bright-soliton stripe Kuznetsov:1986 ; Kivshar-Peli:2000 ; fermions . This instability was originally studied by Zakharov and Rubenchik Zakharov-Rubenchik:1974 —and later addressed in other works, see for instance Refs. Peli:2001 ; Bernard-Peli:2006 — using a perturbative approach for small wavenumbers $k$ of the snaking modulation, up to the second order. The results from these studies, naturally limited to small wavenumbers, incorrectly predicted that the snaking instability was only present in a window starting at $k=0$ (mainly because of the introduction of a wavenumber cutoff), when it is indeed present at all wavenumbers. These conclusions initially led to some misconceptions like the fact that transverse instabilities were only limited (as it occurs in the elliptic case) to a finite range of wavenumbers and, thus, that short-wavelength transverse perturbations might be stable. Later, numerical and analytical studies Cohen:1976 ; Saffman-Yuen:1978 ; Anderson:1979 ; Bernard-Peli:2006 ; Alexe confirmed that the snaking and necking instabilities are indeed present beyond the window predicted in Ref. Zakharov-Rubenchik:1974 (see Ref. Peli2 for a rigorous proof of these facts). In fact, Ref. Bernard-Peli:2006 , using an Evans-function technique Kapitula-Sandstede:1998 , has clarified the origin of a secondary instability that precludes the stabilization of short-wavelength perturbations.

The, previously known, instability properties for the bright-soliton stripe in the hyperbolic-dispersion case are summarized in Fig. 1 where the instability branches are denoted by ‘S’ and ‘N’ for snaking and necking instabilities, respectively. The instability spectrum can be numerically obtained by perturbing the exact stationary solution, $u_{0}(x)$, given by Eqs. (2) and (3) with $x_{0}=v=\theta_{0}=0$. In what follows below, we select the stationary solution by fixing $\mu=1/2$ in Eq. (3). To perform the stability analysis, $u_{0}$ is perturbed as per the standard expression JiankeBook:2010 ,

[TABLE]

where $(\cdot)^{\star}$ stands for complex conjugation, and $\lambda$ yields the instability growth rate (if it has a positive real part) of the transverse modes with real wavenumber $k$. Thus, for each value of $k$, one can compute the spectrum of eigenvalues $\lambda$, substituting expression (4) in Eq. (1), linearizing it with respect to the perturbations and solving the ensuing eigenvalue problem for $\lambda$ and the corresponding eigenvector pair $(a(x),b(x))$.

Further, the stability spectrum and eigenmodes corresponding to the snaking and necking instabilities for the perturbation wavenumber $k=0.5$ are depicted in Fig. 2. As seen in the figure, the snaking perturbation mode, which is responsible for creating transverse undulations of the position of the bright-soliton stripe, naturally features an odd shape in the $x$-direction. On the other hand, the necking instability is, also quite naturally, represented by an even mode, which represents longitudinal modulations of the local density. As mentioned above, one may intuitively expect that the necking instability should not be present if the model is effectively defocusing in the $y$-direction. Nonetheless, the instability of the necking type is present, due to the fact that the eigenvalues that give rise to this instability originate from a collision between the eigenvalue at the edge of the continuous spectrum for small values of $k$ and an eigenvalue bifurcating from the origin (at $k=0$), see Refs. Bernard-Peli:2006 ; Yang-Pelin for details. The respective necking eigenmode inherits properties of its “parent” eigenmodes leading to the bifurcation. Namely, the mode at the edge of the continuum spectrum with small $k$ (the first “parent”) is extended (periodic) along the $x$-direction, while the mode bifurcating from the origin (the second “‘parent”) is an even one, which is essentially localized in the region occupied the underlying steady state. Therefore, the emerging necking mode is an even one, weakly localized (wider than the stationary solution), with oscillatory tails in the $x$-direction. These observations will be crucial when we describe the dynamics of the bright soliton-stripe using the variational approximation in Sec. III.2, where this necking mode cannot be captured. Furthermore, the mixed nature of the necking mode, together with the fact that the real part of its eigenvalue is about twice as small as the one corresponding to the snaking instability, and that necking is not present for small values of $k$, make this mode somehow elusive when performing direct numerical simulations (and, by the same token, elusive in possible experiments), which might be the reason why it was usually missed.

III Reduced effective equations for the soliton stripe

In this section we report two methods to derive effective reduced dynamics of the bright-soliton stripes. The techniques, based on the adiabatic-invariant (AI) approach AI_us:2017 ; AI_us:2018 and the variational approximation progress , rely on a suitable “projection” of the soliton stripe onto a subset of collective coordinates which account for basic features of the stripe, such as, in particular, its center position along the $x$-axis and its dependence on the transverse variable $y$.

III.1 The use of the adiabatic invariant (AI)

We adapt here the AI reduction approach, which was developed for dark-soliton stripes in Ref. AI_us:2017 ; AI_us:2018 , to the case of bright solitons. The method relies on using the conservation of the Hamiltonian to derive an equation of motion for the center of the soliton (along the $x$-direction). In the absence of the external potential ($V=0$), the 2D Hamiltonian corresponding to Eq. (1) is

[TABLE]

To approximate the dynamics of the soliton stripe, we consider the 1D soliton ansatz (2), where the position of the stripe, $X(y,t)$, is assumed to be a function of the transverse coordinate and time:

[TABLE]

Within the AI approximation, we assume that deformations of the stripe weakly affect the velocity, and do not affect the amplitude and width of the transverse soliton’s profile. We have in mind here a nearly stationary soliton scenario, and in assuming an invariant amplitude/width, we attempt to tackle the phenomenon of snaking, rather than that of necking, involving width and amplitude variation. As we are adopting a traveling wave ansatz of the general form $u=f(x-X)\exp(ig(x-X))$, it immediately follows that $|u_{y}|=|X_{y}|\cdot|u_{x}|$, and, therefore, the Hamiltonian (5) can be rewritten as

[TABLE]

The evaluation of this expression for the stripe ansatz given by Eq. (6) after identifying that $v=X_{t}$ and subsequent integration in the $x$-direction yields

[TABLE]

It is important to reiterate here that the steps above rely on the AI assumption that the velocity $v$ is only weakly-dependent, on $y$ (and $t$) through $v=X_{t}(y,t)$.

We now use the fact that the Hamiltonian is conserved, i.e., $dH/dt=0$, to obtain

[TABLE]

We apply the integration by parts in Eq. (8), assuming that the variation of the stripe vanishes at $y\rightarrow\pm\infty$,

[TABLE]

from where it follows that the reduced AI equation for the soliton’s transverse position is

[TABLE]

If we consider small perturbations from the initial straight soliton stripe, the linearization of Eq. (9) yields

[TABLE]

which recovers the well-known snaking instability of the bright soliton stripe Kuznetsov:1986 , which is described by the following dispersion relation, assuming $X\sim\exp(i(ky-\omega t))$,

[TABLE]

We stress that the main result of the AI approach is not to predict the linear instability of the stripe. The AI goes further, as Eq. (9) should in principle allow one to follow the nonlinear evolution of the snaking past its initial instability. In Sec. IV below we will compare the AI prediction, based on Eq. (9), with the full numerically simulated evolution of the original NLS equation (1).

III.2 The variational approximation (VA)

III.2.1 The full VA

As presented above, the AI reduction approach, has an intrinsic shortcoming: as we only use a single conserved quantity (the Hamiltonian), this approach can only produce a single equation of motion. In the case of the dark-soliton stripe, as put forward in Refs. AI_us:2017 ; AI_us:2018 , the single equation, governing the transverse shift of the dark-soliton’s center, is sufficient to describe its dynamics —as the width of the dark soliton is determined by the height of the constant background into which the dark soliton is embedded and the soliton speed. However, in the present case, the bright soliton, in principle, has more degrees of freedom (through the potential variation of its amplitude) and, thus, it will be relevant to develop an approach that can predict the evolution of these degrees of freedom along the stripe. In principle, a more general AI description could be obtained by using the Hamiltonian formulation together with its corresponding canonical variables. However, at this point, it remains elusive in our setting precisely what is the right choice of canonical variables that could describe the soliton filament’s dynamics, so as to use Hamilton’s equations around such a canonical formulation. Therefore, instead of using a Hamiltonian formulation, we turn to a more straightforward Lagrangian one, enabling the derivation of Euler-Lagrange equations for both the center position and the width (or equivalently amplitude) of the solitary wave. In that vein, we adapt a quasi-2D VA methodology, as proposed in Ref. Malomed92 , to develop the full description of the dynamics of the bright-soliton stripe. The VA is based on the Lagrangian of Eq. (1) (this time, it includes the external potential $V$):

[TABLE]

Based on the exact bright-soliton solution of Eq. (2), we introduce the following 2D ansatz for the respective stripe:

[TABLE]

where, in contrast with the AI ansatz (6), we now have four $t$- and $y$-dependent parameters. These four variables represent the amplitude (or inverse width) $A(y,t)$ of the soliton, its location $\xi(y,t)$ [which plays the same role as $X(y,t)$ in the AI ansatz (6)], the $x$-velocity $v(y,t)$, and the phase $\theta(y,t)$ [the solution with $A\left(y,t\right)=$constant, $v(y,t)=$constant, $\xi(y,t)=vt$ and $\theta(y,t)=\frac{1}{2}v^{2}t+\frac{1}{2}A^{2}t+\theta_{0}$ amounts to the exact bright-soliton stripe in the absence of any external potential].

Substituting the VA ansatz (12) in the Lagrangian (11), we perform its $x$-direction integration for the trial function (2), which yields

[TABLE]

with the Lagrangian density

[TABLE]

where we do not explicitly write the dependence on $(y,t)$ of the dynamical variables, and an effective (averaged over the ansatz) potential is introduced:

[TABLE]

where $h(A,r)\equiv A^{2}\operatorname{sech}^{2}\left(Ar\right)$. Note that the effective potential depends on time and the longitudinal coordinate $y$ through its dependence on $A(y,t)$ and $\xi(y,t)$.

Now, following the standard VA method, we perform variations of parameters $\theta$, $v$, $\xi$, and $A$ to produce the following Euler-Lagrange equations. First, it is

[TABLE]

which is tantamount to the continuity equation, which secures the conservation of the squared $L^{2}$ norm (alias total energy, in terms of the optical realization, which is different from the above-mentioned Hamiltonian),

[TABLE]

which diverges (linearly, as a function of the domain size) for the stripes. From Eq. (2) it follows $\int_{-\infty}^{+\infty}|u(x,y,t)|^{2}dx=2A$, hence the integration of Eq. (16) with respect to $y$ produces:

[TABLE]

which corresponds to the energy flux along the $y$-axis.

Next, the variation of the Lagrangian produced by Eqs. (13) and (14) with respect to $v$ and $\xi$ yields

[TABLE]

[TABLE]

where Eq. (16) was used to simplify Eq. (19), and the subscript in $V_{\xi}^{\mathrm{eff}}\equiv\partial_{\xi}V^{\mathrm{eff}}$ denotes the derivative with respect to $\xi$ once the explicit dependence on $A(y,t)$ and $\xi(y,t)$ in $V^{\mathrm{eff}}$ has been introduced. Finally, the variation of the Lagrangian with respect to $A$ produces the following equation,

[TABLE]

where $V_{A}^{\mathrm{eff}}\equiv\partial_{A}V^{\mathrm{eff}}$. In combination with Eq. (18), Eq. (20) yields

[TABLE]

Thus Eqs. (16), (18), (19), and (21) provide a full approximation for the spatio-temporal evolution of the bright-soliton stripe’s parameters $A$, $\xi$, $v$ and $\theta$, respectively. In Sec. IV we compare this approximation to the respective numerical solution of Eq. (1).

III.2.2 Slowly varying solutions

To cast the VA in a more explicit form, we here adopt additional assumptions, neglecting certain terms in Eqs. (16), (18), (19) and (21). Namely, we focus on the slow evolution of the soliton-stripe’s parameters, treating the first derivatives $v_{y}$, $A_{y}$, etc. as first-order small quantities, and neglecting higher-order terms, such as $v_{y}^{2}$, $v_{y}\,w_{y}$, etc., with the notable exception of those corresponding to position $\xi(y,t)$ which are associated with the snaking dynamics. In the framework of this slowly-varying assumption, the set of equations (16), (18), (19) and (21) reduce to the following nonlinear system of simplified coupled equations for the stripe’s parameters:

[TABLE]

In particular, Eq. (23) provides the first order correction, $\sim$ $v_{yy}$, to the lowest-order relation between the position and velocity, $\xi_{t}=v$. Actually, this correction is the improvement provided by the VA in comparison to the (simplest variant developed above of the) AI-based approach.

Note also that, according to the slowly-varying assumption, Eqs. (23) and (24) for the position and velocity decouple from the spatial derivatives of the other two variables, the stripe’s amplitude and phase. This observation implies that, for the slow transverse stripe’s dynamics, the snaking (associated with $\xi$ and $v$) is only weakly coupled to necking (associated with $A$ and $\theta$). In Sec. IV, we corroborate this conclusion by means of the comparison to numerically generated results.

III.2.3 The linear-stability analysis

In the presence of a localized external potential in the $x$ direction [$V(x,y)=V(x)$] in Eq. (1), with a minimum at $x=0$, we can readily expand the potential as $V^{\mathrm{eff}}=V^{\mathrm{eff}}\left(A,0\right)+\frac{1}{2}V_{\xi\xi}^{\mathrm{eff}}\left(A,0\right)\xi^{2}$ for small values of $\xi$, as the first derivative of the effective potential vanishes, $V_{\xi}^{\mathrm{eff}}\left(A,0\right)=-\int_{-\infty}^{+\infty}V\left(x\right)h_{x}\left(A,x\right)dx=0$, due to the parities of the profiles $(x)$ and $h(A,x)$. Therefore, for a fixed value of $A$, the substitution of planar waves in Eqs. (23) and (24), $\xi(y,t)=Be^{i\left(ky-\omega t\right)}$ and $v(y,t)=Ce^{i\left(ky-\omega t\right)}$ yields the dispersion relation,

[TABLE]

[the linearization of full VA equations (18) and (19) leads to the same result]. For instance, in the case of a delta-functional external potential, with $V(x)=-\epsilon\delta(x)$, $V_{\xi\xi}^{\mathrm{eff}}\left(A,0\right)=2\epsilon A^{4}$, the VA predicts stabilization of the snaking perturbations at $k^{2}>3\epsilon A$. It is important to note that the dispersion relation in Eq. (26) gives a global approximation for the instability growth rates as we did not assume any condition on the smallness of the wavenumber $k$. This is in contrast with previous results found by Zakharov and Rubenchik (ZR) using a fourth order perturbation analysis in the wavenumber, which, in the absence of external potential and for unit amplitude, i.e. $A=1$, yields Zakharov-Rubenchik:1974 ; Bernard-Peli:2006 ; Kivshar-Peli:2000

[TABLE]

Figure 3 shows the comparison of the linear-instability spectra, as predicted by the AI and VA approaches, namely Eqs. (10) and (26), together with the perturbative result of Eq. (27), against the full numerical solution of Eq. (1) in the absence of the external potential. As expected, all these spectra approximations asymptotically coincide in the $k\rightarrow 0$ limit. It is relevant to mention that Ref. Kivshar-Peli:2000 also considered this $k\rightarrow 0$ limit where effective equations for the dynamics of the bright stripe were obtained (see Eq. (5.36) in Ref. Kivshar-Peli:2000 ). However, the VA and ZR linear spectra remain approximately valid at larger $k$ as well, contrary to what is the case for the simplified AI model, which is only valid for small $k$. Also, conceived as a global method (i.e., independent of $k$), the VA seemingly gives a better approximation than the ZR perturbative approach of Eq. (27) for wavenumbers close to the critical threshold $k=k_{\rm cr}\simeq 1.08$. However, the VA fails to capture the secondary instability of the snaking branch which the numerical solutions reveal past $k=1$. This deficiency of the VA is explained by the fact that the secondary bifurcation of the snaking branch at $k>1$ involves (as it is also true for the bifurcation of the necking branch) a mode belonging to the continuum spectrum, with slower decaying oscillatory tails, that the VA is not designed to capture in the context of the ansatz of Eq. (12). In Fig. 4 we present two examples corresponding to the destabilization dynamics for the bright soliton stripe due to snaking for long and short wavelengths. As it can be noticed from panels (a) and (c), the destabilization for long wavelengths, corresponding to small values of $k$ before the bifurcation at approximately $k=k_{\rm cr}\simeq 1.08$, keeps the transverse profile of the stripe close to a unimodal (sech-like) hump. On the other hand, see panels (b) and (d), for short wavelengths, $k>k_{\rm cr}$; here, the transverse profile develops strong tails that cannot be captured by the sech-shaped ansatz profile.

In Sec. V we report results obtained in the presence of the external potential. Of particular interest is the eventual suppression of different types of the instability of the soliton stripe by appropriately designed potentials, at different wavenumbers $k$.

IV Numerical results: soliton-stripe dynamics

We now aim to verify the validity of the results predicted by the AI and (especially so) the VA method, by going beyond the prediction for the stability of the bright-soliton stripe. Indeed, as both the AI and the VA generate approximate (reduced) equations of motion for the dynamics of the stripe, the reduced equations may be valid not just for small perturbations around the steady state, but also for large (nonlinear) perturbations. However, for the reasons explained previously, these approximate methods fail to reveal the necking dynamics in the case of the hyperbolic dispersion, therefore below we focus on the snaking dynamics. It is relevant to stress that, as the snaking perturbations have a relatively large instability growth rate, in comparison to necking, it is expected that the snaking dynamics ought to dominate over necking. Below we numerically integrate the underlying equation (1) by means of a combination of the standard second-order finite-difference algorithm in space and fourth-order Runge-Kutta scheme in time. The same results have also been produced using spectral discretization in space. As concerns equations produced by the AI and VA methods, they were solved numerically by dint of standard spectral methods in space together with Runge-Kutta integration in time.

To compare the different levels of the approximation —corresponding to the AI, full VA, and slowly varying VA—, we start with the stationary soliton-stripe solution $u_{0}$ given in Eq. (2) at $t=0$ with zero initial velocity ($v=0$) and centered at $x_{0}=0$ with $A=1$ [i.e., with $\mu=1/2$, as given by Eq. (3)] which we then perturb to create the initial condition in the form of

[TABLE]

where $A_{0}(y)$ and $X_{0}(y)$ are, respectively, the amplitude and position perturbations with amplitudes $\varepsilon_{A}$ and $\varepsilon_{X}$. The perturbation may be (i) snaking applied to the initial position of the soliton stripe along the $y$-direction (i.e., only $\varepsilon_{X}\not=0$), (ii) a necking deformation involving a perturbation of the soliton-stripe’s amplitude in the $y$-direction. (i.e., only $\varepsilon_{A}\not=0$), or (iii) both. As corroboration of the VA results presented in Sec. III.2, we show in panels (a)–(d) of Fig. 5 a typical example of the numerically simulated evolution involving a combination of snaking and necking perturbations. Panels (a) and (c) correspond to a bright-soliton stripe perturbed solely by a snaking perturbation, while panels (b) and (d) correspond to a mixture of snaking and necking. All the perturbations were chosen with the same size and shape (see caption to Fig. 5 for details). The motivation for this set of simulations is to corroborate the conclusion drawn from the VA, according to which the necking only weakly affects snaking. Figure 5 indeed suggests that this is true. In this example, the snaking dynamics (see the top panels) is almost identical for cases initiated with and without necking perturbations. We have also checked that other cases, corresponding to variation of the relative amplitudes of the snaking and necking perturbations, lead to similar results supporting this conclusion. On the other hand, we have observed that if the initial perturbation is purely seeded in terms of necking (with no perturbation in snaking), snaking is not initiated, at least for the time scales considered herein. A typical example showing this is depicted in panels (e) and (f) of Fig. 5. This feature is quite useful as neither the AI or VA reduced model is able to appropriately capture the necking dynamics. Therefore, one may use the AI and VA methods to exclusively follow snaking dynamics without the need to take care of the necking perturbations. Accordingly, in what follows below we focus exclusively on the snaking dynamics. Nonetheless, it is important to mention that necking perturbations are always present (even if not introduced initially), as snaking always induces the necking behavior, cf. the example depicted in the left panels of Fig. 5.

We now compare the snaking dynamics obtained from simulations of Eq. (1) and the AI and VA reductions. An example of this comparison is depicted in Fig. 6. Panel (a) depicts the full simulations for the snaking, along with  results produced by the reductions, as indicated in the caption. All models have the same linearization for long wavelengths, therefore the dynamics for short times ($t<1$) is very similar for both reductions, closely matching results of the direct simulations. However, the AI reduction (shown by green curves) overestimates the growth rate of all perturbations (cf. Fig. 3) and thus its predictions quickly cease to match results of the simulations. In fact, the AI dynamics blows up shortly after $t=6$. On the other hand, the VA predictions closely follow the simulated snaking at longer times. The slowly varying VA (see black curves) is able to predict the snaking up to $t=14$. However, shortly after $t=16$, this approximation blows up, as was the case for the AI reduction. In contrast, the full VA method (see red curves) is able to accurately follow the numerically simulated snaking at all times up to the point were the stripe starts to break up into individual bright patches ($t>14$). It is remarkable that even when the stripe is in the process of breaking up (at $t\simeq 14$), up to the emergence of clearly separated bright patches ($t>16$) the VA reduction still very closely follows the location of the numerically simulated quite intricate pattern of snaking.

To further confirm the validity of the VA, we depict in panels (b) and (c) of Fig. 6 several cross sections at $y=$ constant, as produced by the full simulations of the NLS equation (1), and by the VA reductions. At $t=14$ [panel (a)], the simulated NLS dynamics (solid blue curves) starts to develop asymmetric tails that are more pronounced for values of $y$ where the snaking is further away from $x=0$. Despite the fact that the VA method produces, by construction, a symmetric (i.e., even in $x$) profile, the VA profiles closely follow the numerically simulated counterparts. Notice the importance of the higher order terms in controlling both the center position (which is still reasonably accurate for the VA), but also importantly the amplitude of the solution; the absence of these terms in the reduced VA can be clearly seen to lead to its less accurate representation already at $t=14$. However, at $t=16$ [panel (b)], the slowly varying VA starts to drift away from the numerically produced NLS dynamics, in terms of snaking and necking alike, the latter pertaining to the amplitude of the bright-soliton stripe. At this stage, the amplitudes predicted by the slowly varying VA are several times larger than the original one. Nonetheless, even at this stage, when the numerically simulated $x$-profiles are far away from the VA-based soliton (symmetric) ansatz (see, for instance, the numerically simulated NLS profile at $t=16$ and $y=4\pi/3$), the full VA is able to approximate well the location and amplitude of the stripe. In addition to the typical example displayed in Fig. 6, we have checked that other perturbations lead to similar results (not shown here in detail). Thus, the direct simulations corroborate the validity of the VA method especially in its full form, and the limitations of the simplified AI reduction.

V Using the external potential to control and eliminate instabilities

In this section, we aim to demonstrate how the external potential in Eq. (1) can be used to attenuate and, possibly, eliminate instabilities of the bright-soliton stripe. To suppress the snaking instability, it is natural to propose a channel-shaped potential uniform in the $y$-direction, with a transverse Gaussian profile:

[TABLE]

where $\sigma$ and $V_{0}<0$ are the width and strength of the Gaussian. $V_{0}>0$ (or, similarly, the edge of the 2D medium surface ) may be used to stabilize dark-soliton stripes dark1 ; Adhikari ; Ponz ; Boaz ; dark2 ; dark3 , while it tends to enhance the instability of the bright one (not shown here in detail).

Figure 7 depicts the instability spectra, computed and depicted as in Fig. 2, corresponding to the stability of a stationary bright-soliton stripe with $\mu=1/2$ [see Eq. (3)], in the presence of potential (29) with $\sigma=0.5$ and various values of $V_{0}$ (the insets show the respective transverse profiles of the soliton and potential). As it seen in the figure, even a very weak potential (see case of $V_{0}=-0.05$) shifts the snaking-instability branch from the origin to the right and thus indicates that long-wavelength perturbation cease being unstable. As the strength of the potential channel increases, the snaking branch shifts further to the right and, at the same time, it becomes shorter, showing not only that a smaller interval of the snaking wavenumbers remains unstable, but also that the magnitude of the instability becomes smaller. It is worthy to note that the necking-instability branch is also attenuated by the channel potential, in a similar manner to the snaking branch. When the strength of the potential attains the value $V_{0}=-1.2$ in Eq. (29), the snaking instability is completely eliminated. The stabilization becomes perfect at $V_{0}=-1.4$, when the necking instability is completely suppressed too. We have checked that the randomly perturbed bright-soliton stripe indeed develops no instability at $V_{0}\leq-1.4$ (not shown here in detail). Thus, the appropriately designed channel potential is capable to provide a full stabilization of the soliton stripes.

Some remarks are in order at this stage. Firstly, the surprising stabilization of the necking —while we were trying to suppress snaking— might be understood by recalling that the necking-instability eigenmode in the free space is not fully localized around the soliton’s core, due to the tails attached to it (in the $x$ direction). The channel potential is able to curtail these tails, by trapping the soliton in the channel, thus helping to tame this eigenmode into an “innocuous” one. On the other hand, we have tried other values of the channel’s width $\sigma$ in Eq. (29)], concluding that only values of $\sigma$ close to the soliton’s width provide complete suppression of the instabilities. In particular, narrow potentials tend to produce a secondary instability, whereby modes close to $k=0$ become unstable (not shown here in detail). A detailed analysis of the secondary instability falls outside the scope of the current work.

Finally, the prediction of the stabilization of the bright-soliton stripe against the snaking perturbations in the channel potential, as produced by the VA, has been tested too. Namely, the application of the linear-stability condition, provided by Eq. (26), to the Gaussian potential given by Eq. (29), gives rise to the stability predictions depicted by dashed red curves in Fig. 7. In the figure we only present the VA results for $-0.8\leq V_{0}\leq 0$ since for values of $V_{0}$ below $-0.5$ the VA prediction starts to deteriorate and fails to give as accurate of a prediction for the snaking instability window. This drawback is a direct consequence of the choice of reduced ansatz that we used in Eq. (2). In particular, ansatz (2) assumes that the height and inverse width of the soliton are equal (to $A$). We have checked that this assumption approximately holds for all the statics and dynamics presented herein in the absence of external potential ($V_{0}=0$). However, in the presence of the external potential ($V_{0}<0$), while the height of the bright soliton does decrease as $V_{0}$ decreases, its width does not increase and is kept approximately constant by the confining nature of the potential (see insets in Fig. 7). This shortcoming of the VA could be mended by considering an ansatz including the width of the soliton as a new, independent, variational parameter. However, this approach would also require the addition of a chirp (phase) variational parameter (as a conjugate variable to the soliton width) that would considerably complicate the VA methodology. Such an extension falls outside of the scope of the present manuscript. Nonetheless, it is important to stress that the results presented above suggest that the VA is able to approximate the tendency of the channel potential towards stabilization (i.e., towards reduction of the snaking instability growth rate) very well for moderate potential strengths, with $\left|V_{0}\right|<0.5$.

VI Conclusions & Future Challenges

We have studied the stability and dynamics of transverse instabilities, mainly of the snaking type, for bright-soliton stripes in the 2D NLS equations with the dispersion of the hyperbolic type. Such settings can be readily implemented in terms of the spatiotemporal propagation in self-focusing optical waveguides with the normal GVD (among other physical settings). The analysis reported here fully accounts for the mechanisms underlying the transverse snaking instability of the stripe, with the help of approximate reductions of the underlying NLS equation, based on the AI (adiabatic-invariant) and VA (variational approximation) methods, which are applied to suitable ansätze. While the applicability of the results provided by the simplified AI approximation considered here is limited, the VA method, applied to the soliton stripe with parameters (position, velocity, amplitude, and phase) that depend not only on time, but also on the longitudinal spatial coordinate, yields an accurate description of the snaking stability for the stripe. Furthermore, we have shown that the reduced dynamical equations of motion, which are also produced by the VA, are quite efficient in describing the full nonlinear dynamics of the unstable bright-soliton stripe, all the way to the stage at which the stripe breaks up into individual bright spots. Further, introducing an appropriately designed channel potential, we were able to fully suppress the instabilities, of both the snaking and necking types. As concerns the stabilization of the stripe against the snaking perturbations, the VA provides good estimates of the (in)stability eigenvalues in the case of small up to moderate strengths of the channel potential.

The present paper not only complements previous studies on the subject but also puts forward efficient mechanisms for the control of the instabilities. It is worth mentioning that the leading instability estimates, corresponding to small wavenumbers of the snaking perturbations, originally due to Zakharov and Rubenchik Zakharov-Rubenchik:1974 , are captured and improved by our VA method. We envisage potential applications of this methodology to other higher-dimensional problems, such as dark-soliton stripes, cf. Refs. dark1 ; Adhikari ; Ponz ; Boaz ; dark2 ; dark3 . In the framework of the present model, it may be relevant to construct stripes whose transverse structure is similar not to the ground state of the Gaussian trapping potential (29), but to one of its excited bound states, provided that the finite-depth potential supports such states (in terms of optics, the latter condition corresponds to the condition that the channel waveguide is a multimode one multimode ). In fact, the self-trapping nonlinearity may help to create bound states even if they do not exist in the linear limit Carr . A challenging possibility is to develop the analysis for bright-soliton filaments in 3D media with the hyperbolic dispersion, where the stabilization may be provided by a fiber-like trapping channel burnt in the bulk material (an experimental technique which can creates such channels is well known Jena ). Another significant challenge consists of improving the adiabatic invariant methodology, by utilizing an amplitude dependent ansatz in it and considering the resulting canonical formalism and the associated equations of motion. This would be a significant addition to the arsenal of methods of relevance to both bright and dark filamentary structures. Such studies are currently in progress and will be reported in future publications.

Acknowledgements

LACA gratefully acknowledges the permit of absence from IPN and the financial support of the Fulbright-García Robles program during his research stay at SDSU. PGK and RCG gratefully acknowledge the support by the National Science Foundation under grants PHY-1602994 and PHY-PHY-1603058.