Modulation stability analysis of exact multidimensional solutions to the generalized nonlinear Schrödinger equation and the Gross-Pitaevskii equation using a variational approach

Nikola Z. Petrović; Najdan B. Aleksić; Milivoj Belić

doi:10.1364/OE.23.010616

1. Introduction

Modulation instability (MI) is a general phenomenon occurring in solutions of the nonlinear wave equation wherein small perturbations to the solution grow exponentially over time, often producing singularities [1]. Determining whether MI occurs is of prime interest in the field od non-linear optics, where many different forms of wave equation naturally occur [2,3]. In particular, the generalized nonlinear Schrödinger equation (NLSE) is an important generic model that is of great use in NL optics [4–7]. The MI in the NLSE with cubic nonlinearity has been studied in [8], as well as in various related systems: the discrete NLSE [9], the NLSE with loss and a derivative term [10], the cubic-quintic NLSE [11], and others [12, 13]. Experimental confirmation of MI in the NLSE with varying coefficients was given in [14]. MI has also been extensively studied in metamaterials, especially in materials with negative index of refraction [15–23]. In particular, it was determined that metamaterials have different stability properties from ordinary materials [15]. The discovery of metamaterials has opened up the possibility of forming stable spatial solutions through the management of the sign of refraction.

Recently, a number of exact solutions have been found for the various forms of NLSE [24–26] and for the Gross-Pitaevskii (GP) equation in Bose-Einstein condensates [27,28], using the F-expansion technique and the principle of harmonic balance [29, 30]. Two types of solutions were found: travelling wave solutions and solitary wave (SW) solutions. The solutions developed in [25] and [26] were for the (3+1)-dimensional ((3+1)D) NLSE with anomalous and normal dispersion, respectively. The solutions developed in [24] were for the (2+1)D NLSE. Meanwhile, the solutions in [28] describe stable spatiotemporal SWs under the influence of a sinusoidal diffraction/dispersion parameter, while the solutions in [27] required external gain in order to maintain the amplitude. Reference [28] briefly addressed the stability of obtained solutions, while [25] and [26] left the stability analysis to be performed in future work. It is the goal of this paper to provide a complete stability analysis of the obtained solutions in all of these papers. While this may look like a restricted goal aimed at specific solutions, it is nonetheless useful because it might be applicable to multidimensional solutions found in other equations that are of questionable stability. To the best of our knowledge, this is the first time MI analysis of this form has been applied to these systems.

The stability of exact soliton solutions to the NLSE in various forms is an important question that requires careful and thoughtful answer [31]. In (1+1)D, bright and dark soliton solutions of the NLSE in Kerr medium with cubic nonlinearity are unconditionally stable for the, respectively, self-focusing and self-defocusing nonlinearity [4]. However, in homogenous bulk media with a self-focusing cubic Kerr nonlinearity one cannot have unconditionally stable solutions of the NLSE in two and three dimensions. Nevertheless, great interest has been generated when it was suggested that (2+1)D generalized NLSE with varying coefficients may lead to stable 2D solitons [32]. The stabilizing mechanism has been the sign-alternating Kerr nonlinearity in a layered medium. The solutions obtained in [24–26, 28] all resulted in the alternating sign of Kerr nonlinearity and therefore it is worth investigating whether those solutions are also stable. In addition, there is strong indication that the SW solutions can be combined into multiple solutions using the self-similar method [33, 34] and that the individual components can interact with each other without affecting each other’s form [34], which is a defining characteristic of solitons. Still, we will only use the term “solitary wave solution” to describe these solutions throughout our paper.

We will use a variational approach described in [35,36] to explore the modulation instability of the solitary wave solutions obtained in [24–26, 28]. Localized two- and three-dimensional solutions of the cubic nonlinear Schrödinger equation in [24–26] are extended one dimensional solutions. Namely, the intensity of solutions is homogenous in two out of three spatial dimensions for solutions obtained in [25, 26]. In the case of solutions obtained in [24] it is homogeneous in one out of two dimensions. It is in these homogenous directions, due to the nonlinearity, that modulation instability can develop. For this reason, of particular interest is the analysis of modulation stability of solutions in the direction of homogeneity.

2. Generalized nonlinear Schrödinger equation

We confine our analysis to the (2+1) and (3+1)-dimensional NLSEs considered in [24–26], and use the notation introduced there. We consider the generalized (3+1)D NLSE with varying coefficients and Kerr nonlinearity, developed in [25, 26]

i \frac{\partial u}{\partial z} + \frac{β (z)}{2} (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + s \frac{\partial^{2} u}{\partial t^{2}}) + χ (z) {| u |}^{2} u = i δ (z) u .

The functions β, χ, and δ stand for the diffraction, nonlinearity, and the gain or loss coefficients, respectively. Our goal is to verify whether the solutions to Eq. (1) developed in [25, 26] are modulationally stable (MS) or modulationally unstable (MU). In Eq. (1), s = −1 for the normal dispersion and s = 1 for the anomalous dispersion. For s = 0 we have the two-dimensional time independent NSLE studied in [24], which will also be covered in our stability analysis.

The solution to Eq. (1) described in [26] is given as:

\begin{array}{l} u & = & {(α)}^{3 / 2} f_{0} e^{\int_{0}^{z} δ d z} (F (θ) + ε \sqrt{\frac{c_{0}}{c_{4}}} \frac{1}{F (θ)}) \\ exp (i (a (z) (x^{2} + y^{2} + s t^{2}) + b (z) (x + y + t) + e (z))), \end{array}

where F is a Jacobi elliptic function (JEF), α is the chirp function:

α = \frac{1}{1 + 2 a_{0} \int_{0}^{z} β d z}

and:

θ = k (z) x + l (z) y + m (z) t + ω (z),

for nonlinearity χ(z) given as:

χ (z) = - c_{4} \frac{β (z)}{α f_{0}^{2}} χ_{0} exp (- 2 \int_{0}^{z} γ d z),

where we define:

χ_{0} = (k_{0}^{2} + l_{0}^{2} + s m_{0}^{2}) .

The parameters a, b, k, l, m, ω and e are:

a = α a_{0}, b = α b_{0},

k = α k_{0}, l = α l_{0}, m = α m_{0},

ω = ω_{0} - α (k_{0} + l_{0} + s m_{0}) b_{0} \int_{0}^{z} β d z,

e = e_{0} + (α / 2) (c χ_{0} - (2 + s) b_{0}^{2}) \int_{0}^{z} β d z,

where

c = c_{2} - 6 ε \sqrt{c_{0} c_{4}}

and c₀, c₂ and c₄ are parameters related to the Jacobi elliptic function parameter M [24]. We will limit our attention to the cases where ε = 0, hence c = c₂. Throughout this paper, we also take ω₀ = e₀ = 0. For convenience, we define the parameter ē as follows:

\bar{e} = e - \frac{c}{2} \int_{0}^{z} β α^{2} d z .

We now make the following gauge transformation on the solution and the coordinates:

\begin{array}{l} u \to G & = & u exp (- \int_{0}^{z} δ z) \cdot \\ exp (- i (a (x^{2} + y^{2} + s t^{2}) + b (x + y + t) + \bar{e})) / (f_{0} α^{3 / 2} {| χ_{0} c_{4} |}^{1 / 2}), \end{array}

x \to x^{'} = α (x - ς),

y \to y^{'} = α (y - ς),

t \to t^{'} = α (t - s ς),

z \to z^{'} = \int_{0}^{z} α^{2} β d z,

where

ς (z) = b_{0} \int_{0}^{z} β d z

is the solution to the equation:

\frac{\partial ς (z)}{\partial z} = β (z) (2 a (z) ς (z) + b (z)),

to obtain the following equation:

i \frac{\partial G}{\partial z^{'}} + \frac{1}{2} (\frac{\partial^{2} G}{\partial x^{' 2}} + \frac{\partial^{2} G}{\partial y^{' 2}} + s \frac{\partial^{2} G}{\partial t^{' 2}}) + σ {| G |}^{2} G = 0,

where σ = sgn(−c₄χ₀). Equation (18) is much more suitable for stability analysis than Eq. (1) because all the coefficients next to G and its derivatives are constant. Also, the wave propagation now necessarily happens along a straight line, unlike in solutions shown in [25, 26].

The stationary solutions $F = G exp (- i c_{2} \int_{0}^{z} β α^{2} d z / 2)$ , where G is the solution of Eq. (18), contain the whole range of solutions from [25] and [26], for different values of c₂ and c₄, in the direction (k₀, l₀, m₀). Since ε = 0, F will be equal to some JEF. Adjusting parameters (k₀, l₀, m₀) to these values corresponds to the rotation and re-scaling of the coordinate system (x′, y′, t′). Without loss of generality, we can put f₀ = 1 and also k₀ = 0, l₀ = 0, m₀ = 1, for temporal solutions, and k₀ = 0, l₀ = 1, m₀ = 0, for spatial solutions. As mentioned, the amplitude |G| is homogenous in two of the three spatial/temporal dimensions (i.e. in the plane perpendicular to the direction (k₀, l₀, m₀) of inhomogeneity.) It is in this plane that, owing to nonlinearity, the modulation instability can develop. Hence, it is of interest to analyze MI of perturbations in the plane of homogeneity of |G|. We chose x′ as the perpendicular direction of modulation. In this case, χ₀ = 1 for spatial SWs, whereas χ₀ = s for temporal SWs. If the SW is spatial, perturbations along the temporal axis t′ are also possible. In that case, we will consider the x′ and t′ to be swapped, so that x′ always denotes the axis of perturbation.

3. Variational approach to the modulation stability

We now use the variational approach to examine MI, following [8]. We analyze two periodic traveling wave planar solutions: F = sn and F = cn, which reduce in the limit of M = 1 to two solitary wave solutions: tanh (the dark SW) and sech (the bright SW), respectively. We study spatial and temporal periodic solutions for both cases: the normal and the anomalous dispersion.

The idea of the variational approach is to introduce a perturbation in the solutions and analyze the behavior of the perturbation. To this end, we assume:

G = G_{0} (1 + U (z) cos (K x^{'})),

where G₀ is the unperturbed solution of Eq. (18), U(z) = U_r(z) + iU_i(z) is the total amplitude perturbation, K is the wave number of the perturbation and x is the direction in which the perturbation occurs, as described in the previous section. One then constructs, according to standard procedure, the corresponding Lagrangian to Eq. (18):

L = \frac{i}{2} (G \frac{\partial G^{*}}{\partial z^{'}} - G^{*} \frac{\partial G}{\partial z^{'}}) + \frac{1}{2} {| \nabla G |}^{2} - σ {| G |}^{4},

where G^* is the complex conjugate of G, and |∇G|² = |∂G/∂x′|² + |∂G/∂y′|² + s|∂G/∂t′|². One next performs an averaging of the Lagrangian over all three transverse coordinates, to obtain:

〈 L 〉 = \int (\frac{i}{2} (G \frac{\partial G^{*}}{\partial z^{'}} - G^{*} \frac{\partial G}{\partial z^{'}}) + \frac{1}{2} {| \nabla G |}^{2} - σ {| G |}^{4}) d x^{'} d y^{'} d t^{'} .

Note that the Lagrangian is averaged over one period of perturbation in the direction of perturbation and in the direction of the SW it has been averaged from −2K(M) to 2K(M) for F = cn and from −K(M) to 3K(M) for F = sn (these boundaries converge to −∞ and ∞ for M = 1, i.e. to solitary waves). Here

K (M) = \int_{0}^{π / 2} {(1 - M {sin}^{2} t)}^{- 1 / 2} d t

is the complete elliptic integral of the first kind and M is the parameter of JEFs. The total action is now defined as:

Λ = \int_{- \infty}^{\infty} 〈 L 〉 d z .

It remains invariant in the transformation of coordinates from u to G.

Substituting the new formula for G in Eq. (19) into the effective Lagrangian given in Eq. (22), we vary U_r and U_i in the standard procedure [8], to obtain Euler-Lagrange differential equations for U_r = U_r(z) and U_i = U_i(z), as follows:

\frac{\partial}{\partial z} U_{r} = \frac{1}{2} K^{2} α^{2} β (κ U_{i}),

\frac{\partial}{\partial z} (κ U_{i}) = - \frac{1}{2} (K^{2} - κ σ d) α^{2} β U_{r},

where the parameter κ is defined as follows: κ = 1 for the perturbations along the spatial coordinates (spatial perturbations) and κ = s for the temporal perturbation of spatial SWs. The parameter d is defined as:

d = d_{cn}^{(M)} = \frac{8}{3} \frac{(2 M - 1) (E (M) - E (am (5 K (M) | M) | M)) - 2 (2 - 5 M + 3 M^{2}) K (M)}{E (M) - E (am (5 K (M) | M) | M)) - 4 (M - 1) K (M))}

for F = cn(·|M),

d = d_{dn}^{(M)} = \frac{8}{3} (2 - M + (1 - M) \frac{K (M)}{2 E (M)})

for F = dn(·|M) and

d = d_{sn}^{(M)} = \frac{8}{3} \frac{(M + 1) (E (am (4 K (M) | M) | M)) - 2 (2 + M) K (M)}{(E (am (4 K (M) | M) | M) - 4 K (M))}

for F = sn(·|M). Here K(M) = F(π/2|M) and E(M) = E(π/2|M) are the complete elliptic integrals of the first and second kind, respectively;

F (u | M) = \int_{0}^{u} {(1 - M {sin}^{2} t)}^{- 1 / 2} d t

and

E (u | M) = \int_{0}^{u} {(1 - M {sin}^{2} t)}^{1 / 2} d t

are the incomplete elliptic integrals of the first and second kind, respectively, and am(u, M) = F⁻¹(u, M) is the amplitude of the Jacobi elliptic functions. The dependence of coefficients

d_{cn}^{(M)}

,

d_{dn}^{(M)}

and

d_{sn}^{(M)}

on the elliptic parameter M (0 < M < 1) is shown in Fig. 1(a). For bright SWs

d_{cn}^{(M = 1)} = d_{dn}^{(M = 1)} = 8 / 3

and for dark SWs

d_{sn}^{(M = 1)} = 4

.

Fig. 1 (a) Nonlinearity parameter d for solutions cn, sn and dn. (b) The growth rate parameter γ for dark an bright SWs, as a function of K for the case κσ = 1. Modulational instability occurs for values of K depicted in the respective graphs. The solid lines represent the theoretical calculation of K using Eq. (30), and the square and circle dots are values of γ measured using numerical simulations, in which the dark and bright SWs, respectively, were perturbed by a small wave of the given wave number K.

Download Full Size | PDF

The solution to Eqs. (23) and (24) can now be written as:

U_{r} = U_{0} cosh (\frac{γ ξ}{1 + 2 a_{0} ξ}),

U_{i} = U_{0} \frac{2 γ}{K^{2}} sinh (\frac{γ ξ}{1 + 2 a_{0} ξ}),

where:

γ = K \sqrt{(σ κ d - K^{2})} / 2

and

ξ = \int_{0}^{z} β d z .

For small values of ξ, i.e. ξ << 1, the modulus of the perturbation amplitude can be approximated to within second order of ξ to be:

| U | = U_{0} (1 + σ κ \frac{d γ^{2}}{2 K^{2}} (1 - 4 a_{0} ξ) ξ^{2}) .

3.1. Case without chirp

In the case without chirp, i.e. for a₀ = 0, the solutions given in Eqs. (28) and (29) to Eqs. (23) and (24) become:

U_{r} = U_{0} cosh (γ ξ),

U_{i} = U_{0} \frac{2 γ}{K^{2}} sinh (γ ξ) .

The dynamics of the overall evolution of the total perturbation U are determined by the growth rate parameter γ given in Eq. (30), also known as the modulation-instability growth rate [4], and the function β(z). For κσ = −1, the growth rate parameter

γ = i \bar{γ} = i K \sqrt{K^{2} + d} / 2

is imaginary for all values of K, and, consequently, the solutions G are modultionally unconditionally stable for any function β(z). This case occurs for temporal cn-SWs with normal dispersion (s = −1) and temporal sn-SWs with anomalous dispersion (s = 1) in the self-defocusing media.

In the opposite case, κσ = 1, which holds for the temporal or spatial cn-solutions for anomalous dispersion s = 1 or temporal sn-solutions if s = −1 in the focusing media, more interesting dynamics of perturbations occur. For the analysis of these dynamics, let us assume β to be of the following form: β(z) = β₀ + β₁ sin(2πz/Z), where Z = 2π/Ω is the wavelength of β. For the spatial perturbation with wavenumber $K > \sqrt{d}$ the growth rate is zero (since γ = iγ̄ is imaginary) and the perturbation amplitude has an oscillatory solution in the following range:

\sqrt{1 - \frac{d}{K^{2}} {sin}^{2} (\bar{γ} β_{1} Z / π)} \leq | U / U_{0} | \leq 1 .

In Fig. 2 we see that the mode corresponding to

K = 2 > \sqrt{8 / 3}

is stable for both a constant and a sinusoidal form of β (plots (a) and (b), respectively).

Fig. 2 Perturbation amplitude growth for κσ = 1, a₀ = 0 and d = 8/3 as a function of propagation distance z. Black curves are numerical results, while the red curves are analytical results: (a) β₀ = 1, β₁ = 0, top to bottom: numerical results for $K = \sqrt{d / 2}$ , analytic results for $K = \sqrt{d / 2}$ , analytical results for K = 2, numerical results for K = 2, (b) β₀ = 0, β₁ = 1, Z = 1, top to bottom: numerical results for $K = \sqrt{d / 2}$ , analytical results for $K = \sqrt{d / 2}$ , analytical results for K = 2, numerical results for K = 2.

Download Full Size | PDF

If $K < \sqrt{d}$ , then |U| grows exponentially at a rate of γβ₀, as seen in Fig. 2(a). Consequently, if β₀ ≠ 0 the solution is unstable. In Fig. 1(b) we plot γ as a function of wavenumber K for bright (d = 8/3) and dark (d = 4) SWs in the case without dispersion management: β(z) = β₀ = 1. In both cases, γ has a maximum γ = Γ = d/4 for modes $K = \sqrt{d / 2}$ . Numerical simulations of Eq. (18) confirm the analytical prediction for the growth rate γ. It also follows that the amplitude of |U| can be made to be stable if the mean value of the management function is zero β₀ = 0, and the period of oscillations of β is small (i.e. Z << 1). This is due to the finite limit imposed on ξ under these conditions. The variation of the perturbation amplitude in this case is:

1 \leq | U / U_{0} | \leq \sqrt{1 + \frac{d}{K^{2}} {sinh}^{2} (\bar{γ} β_{1} Z / π)} .

3.2. Case with chirp

In the case where a₀ > 0 we analyze the evolution of perturbation modes with maximal growth rates as functions of z. Expanding the solution for large ξ to within the first order of 1/ξ we find:

| U | = U_{0} \sqrt{1 + σ κ \frac{d}{K^{2}} {sinh}^{2} (\frac{γ}{2 a_{0}})} (1 - \frac{C}{a_{0} ξ}),

where:

C = \frac{γ}{8 a_{0}} \frac{σ κ \frac{d}{K^{2}} sinh (\frac{γ}{a_{0}})}{1 + σ κ \frac{d}{K^{2}} {sinh}^{2} (\frac{γ}{2 a_{0}})} .

This implies saturation after linear growth for constant β, i.e. β₀ ≠ 0 and β₁ = 0. Indeed, some sort of saturation of the perturbation amplitude is to be expected in this case, since the solution dissipates due to the form of the chirp function [25]. The analytical results agree well with the numerical results, as can be seen in Fig. 3(a).

Fig. 3 Perturbation amplitude growth for κσ = 1, d = 8/3 and $K = \sqrt{d / 2}$ as a function of propagation distance z for systems with chirp. Black curves are numerical results, while the red curves are analytical results: (a) β₀ = 1, β₁ = 0, a₀ = 0.1, top: numerical results, bottom: analytical results, (b) β₀ = 0, β₁ = 1, Z = 1, dashed lines represent plots for a₀ = 0.1, top to bottom: analytical results for a₀ = 0.1, numerical results for a₀ = 0.1, analytical results for a₀ = 0.3, numerical results for a₀ = 0.3.

Download Full Size | PDF

The overall behavior of MIs is otherwise similar to the case without chirp. Typical behavior of |U| is presented in Fig. 3, for different values of parameters. In the case of dispersion management β₀ = 0, the perturbation amplitude has oscillatory behavior with a maximum variation that depends on the period Z of the management function β. Stabilization can be achieved by reducing the period of the management function. The dependence of the maximum perturbation amplitude on z, d and a₀ is given in Fig. 4. While an increase in parameter d in general increases the amplitude of |U|, as is seen in Fig. 4(a), an increase in a₀ reduces the amplitude of |U|, as seen in Fig. 4(b). This result is in agreement with results obtained in Fig. 3(b) where values a₀ = 0.3 and a₀ = 0.1 were compared.

Fig. 4 Maximum amplitude of perturbation for $K = \sqrt{d / 2}$ plotted against z = β₁Z/π and: (a) d for a₀ = 0.05, (b) a₀ for d = 8/3.

Download Full Size | PDF

The entire stability analysis is presented in Table 1. We see that depending on the choice of s, σ and whether the SW is spatial or temporal we have eight distinct cases for examining the stability of our solutions. In the case of spatial SWs, perturbations can occur in both the spatial and temporal directions, whereas in the case of temporal SWs they can only occur in spatial directions. For a SW to be stable we must have κσ = −1 in all directions of perturbation, otherwise it is conditionally stable, i.e. only for β₀ = 0. A SW is dark if the direction of the SW and the nonlinear term, i.e. σ are of the opposite sign, otherwise the SW is bright. In the 2D time-independent case we no longer need to consider temporal SWs and distortions.

Table 1. Stability cases

View Table

4. Numerical simulations

We now use computer simulations to simulate the behavior of our solutions when a small perturbation is introduced. Observing the rate of change of the amplitude of G in our simulations, we can then measure the value of γ and compare it with the theoretical expectation given in Eq. (30). The split-step FFT simulations produced the points on the plot in Fig. 1(b), which agree well with theoretical expectations denoted by the continuous line. Most importantly, the solutions cease to exponentially increase at precisely the values predicted by the theory of MI.

We see in Figs. 5 – 7 the main results of our simulations in scenarios involving instability. Starting from the initial form of the solution (Figs. 5–6(a)), we can see that the perturbation rapidly increases (Figs. 5–6(b)) and ultimately completely abandons the original form of the solution (Figs. 5–6(c)). Figure 5 depicts the time evolution of a bright SW, while Fig. 6 depicts the time evolution of a dark SW. Owing to the difficulty that arises from the boundary conditions for the dark SW, due to a change in the sign of u, we have instead run a simulation of two dark SWs with periodic boundary conditions, as is the standard practice. Periodic solutions show a similar pattern of instability formation as seen in Fig. 7.

Fig. 5 Development of modulation instability for the bright SW for three different values of z. Here, x is the direction of perturbation, y is the direction of the SW and t is the remaining transverse direction. Bright colors, i.e. towards the color red (the center in Fig. 5(a)), indicate a higher value of |u|².

Download Full Size | PDF

Fig. 6 Development of modulational instability for the dark SW for three different values of z. Here, x is the direction of perturbation, y is the direction of the SW and t is the remaining transverse direction. Red color (away from the center in Fig. 6(a)), indicates a higher value of |u|².

Download Full Size | PDF

Fig. 7 Development of modulational instability for the dark traveling wave (F = sn) for three different values of z. Here, x is the direction of perturbation, y is the direction of the traveling wave and t is the remaining transverse direction. Parameters are M = 0.5 and $K = \sqrt{d / 2}$ . Blue color (at the top, bottom and the three central stripes in Fig. 7(a)), indicates a lower value of |u|².

Download Full Size | PDF

5. Analysis of the stability of the Gross-Pitaevskii equation

In this section we apply the results and methods of Sec. 2 to provide a stability analysis of our solutions to the Gross-Pitaevskii equation (GPE), obtained in [28]. This section expands upon the results briefly summarized in [28] and provides additional data. We now examine the GPE, i.e.

i \partial_{t} u + \frac{β (t)}{2} (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + \frac{\partial^{2} u}{\partial z^{2}}) + χ (t) {| u |}^{2} u + η (t) (x^{2} + y^{2} + z^{2}) u = i δ (t) u .

The coefficient η refers to the strength of the quadratic potential. The main part of our analysis will be to transform the starting Eq. (39) into a form more amenable to stability analysis [38]. The propagation variable is now t instead of z.

The key differences between Eq. (39) and Eq. (1), apart from the addition of the quadratic potential, is the change in the longitudinal direction from z to t. Hence, there is no longer a distinction between normal and anomalous dispersion. This greatly simplifies the stability analysis. We restrict out attention to SWs found in [28], as the solutions found in [27] do not have a stable amplitude when they are not artificially maintained with a nonzero gain. The solutions in [28] are found under the condition that β and η are proportional trigonometric functions of the same sign. We will call their amplitudes β₀ and η₀. The parameters b, k, l and m are:

b = p b_{0},

k = p k_{0}, l = p l_{0}, m = p m_{0},

ω = ω_{0} - q (k_{0} + l_{0} + m_{0}) b_{0} .

The parameters p and q, as well as the chirp function a, are defined to be:

p = \sqrt{\frac{η_{0}}{η_{0} - 2 a_{0}^{2} β_{0}}} sech (τ (t) + τ_{0}),

q = \frac{\sqrt{η_{0} β_{0}}}{\sqrt{2} (η_{0} - 2 a_{0}^{2} β_{0})} (tanh (τ (t) + τ_{0}) - tanh τ_{0}),

a = \sqrt{\frac{η_{0}}{2 β_{0}}} tanh (τ (t) + τ_{0}),

where:

τ_{0} = arctanh (a_{0} \sqrt{\frac{2 β_{0}}{η_{0}}})

and:

τ (t) = \sqrt{\frac{2 α_{0}}{β_{0}}} \int_{0}^{t} β (t) d t .

Formula (5) now becomes:

χ (t) = - c_{4} \frac{β (t)}{p f_{0}^{2}} χ_{0} exp (- 2 \int_{0}^{t} δ d t),

where

χ_{0}^{2} = k_{0}^{2} + l_{0}^{2} + m_{0}^{2}

and p is defined in Eq. (43). For convenience we set χ₀ = 1, f₀ = 1 and β₀ = 1. We also define:

\bar{e} = e - \frac{q c}{2},

where q is given in Eq. (44) and

c = c_{2} - 6 ε \sqrt{c_{0} c_{4}}

. Again, for ε = 0 we have c = c₂.

We now define:

G = \frac{u}{p^{3 / 2} f_{0} \sqrt{| c_{4} χ_{0} |}} exp (- \int_{0}^{t} δ d t) exp (- i a (x^{2} + y^{2} + z^{2}) - i b (x + y + z) - i \bar{e})

and employ similar changes of coordinates as in Sec. 2:

x \to x^{'} = p (x - ζ)

y \to y^{'} = p (y - ζ)

z \to z^{'} = p (z - ζ)

t \to t^{'} = \int_{0}^{t} p^{2} β d t,

where p is defined in Eq. (43) and

\frac{\partial ζ (t)}{\partial t} = β (t) (2 a (t) ζ (t) + b_{0} p (t))

. The form of the function ζ(t) in the transformation is of no immediate interest, other than the fact that it depends on p(t) and τ(t), which are given in Eq. (43) and Eq. (47), respectively. The transformation gives us the following equation for G:

i \frac{\partial G}{\partial t^{'}} + (\frac{\partial^{2} G}{\partial x^{' 2}} + \frac{\partial^{2} G}{\partial y^{' 2}} + \frac{\partial^{2} G}{\partial z^{' 2}}) - σ {| G |}^{2} G = 0,

where σ = −sgn(c₄). Qualitatively, this is the same equation as Eq. (18), except for the change of the longitudinal variable from z to t. The new variable t′, which only depends on t, involves an integral over β that can change sign. This will be important in the analysis of Eq. (55). Just as in the case of NLSE, we can place without the loss of generality the z′ axis in the direction of inhomogeneity of our extended solitary solutions, i.e. assume k₀ = l₀ = 0 and m₀ = 1, and put x′ as the axis of perturbation.

Equation (55) is the usual (3+1)D nonlinear Schrödinger equation with constant coefficients, which is prone to instabilities and the wave function collapse. Instabilities in G translate into instabilities of the general solution u. This would bode disaster for the stability of exact traveling wave and solitary solutions found, were it not for the possibility of diffraction and nonlinearity management [7] in Eq. (55), thanks to the form of the primed variables. We find that, for the choice of coefficients α(t) and β(t) made in [28], the typical extended SW solutions of Eq. (55) do not collapse when perturbed, but keep oscillating in a typical breathing behavior.

We now consider the perturbation of G in this plane for the two fundamental solutions, the dark F = sn and the bright F = cn SWs, where F = Gexp(−iqc₂/2) in the form:

G = G_{0} (1 + U (t) cos (K x^{'})),

where U(t) = U_r(t) + iU_i(t) is the complex amplitude, and K is the wavenumber of the perturbation in the direction perpendicular to z′. In a standard linear stability analysis, as was already done for the NLSE in Sec. 2, the perturbation is substituted into Eq. (55) and linear first-order differential equations for U_r and U_i are obtained. Plugging in the perturbation, one obtains:

\frac{\partial U_{r}}{\partial t} = \frac{1}{2} K^{2} p^{2} β U_{i},

\frac{\partial U_{i}}{\partial t} = - \frac{1}{2} (K^{2} - σ d) p^{2} β U_{r},

where d is defined as in Eqs. (25)–(27). The solutions of Eqs. (57) and (58) determine the dynamics of the modulational instability. Equations (57) and (58) can be solved analytically to yield:

U_{r} (t) = U_{0} cosh (γ q (τ)),

U_{i} (t) = U_{0} \frac{γ}{K^{2}} sinh (γ q (τ)),

where:

γ = K \sqrt{(d - K^{2})} / 2 .

A similar analysis as in Sec. 3.1 can be performed in this case, and one obtains a much more simplified analysis with respect to Table 1 since there is no case s = 1, and there are no temporal perturbations. One obtains that the dark SWs are always stable, and the bright SWs are always conditionally stable.

We now restrict our attention to the bright SWs, i.e. σ = 1. The the modulus of the perturbation amplitude is given as:

| U | = U_{0} {(1 + \frac{σ d}{K^{2}} {sinh}^{2} (γ q (τ)))}^{1 / 2} .

A graph and a detailed analysis describing the behavior of perturbations were given in [28]. Here we focus our attention on the dependence of the maximum perturbation on parameters d, η₀ and a₀. It is worth noting that unlike in the case of NLSE the parameter q is limited assuming a reasonable choice of η₀ and a₀. Thus assuming η and β are proportional and of equal sign, finding stable solutions is likely even in the absence of dispersion management.

To achieve the maximum amplitude of |U| we take t → ∞ and $K = \sqrt{d / 2}$ . In that case, along with β₀ = 1, we obtain

γ q = - \frac{d}{4 \sqrt{2} (\sqrt{η_{0}} - \sqrt{2} a_{0})} .

Since σd/K² = 2 and sinh is an increasing function it follows that growth of |U| follows the growth of the magnitude of γq. In other words, d and a₀ contribute to an increase in |U| while η₀ contributes to a decrease in |U|. As solutions approach the singularity threshold

a_{0} = \sqrt{η_{0} / 2}

, the amplitude of perturbations also blows up.

6. Conclusion

In this paper, we have analyzed the stability of solutions of the (3+1)D NLSE with normal or anomalous dispersion and the (2+1)D time-independent NLSE. For the (2+1)D solutions, we obtained stability for dark solitary waves and conditional stability for bright solitary waves, meaning we need to apply dispersion management to keep the solitary waves stable, i.e. the diffraction/dispersion coefficient must oscillate around 0. The management function dynamically stabilizes the nonlinear structure of the transversal perturbation of the soliton if its mean value is zero. Reducing the period of the management function can be achieved by arbitrary limitations of the perturbation level. If, however, the mean value is different from zero, the amplitude modulation perturbation exponentially increases in the case of a solution without chirp and linearly with saturation in the case with chirp. In the former case we have Lyapunov instability. In the latter case the saturation is an indirect consequence of the dissipation of the solution.

For the (3+1)D case we obtain stability for the temporal bright solitons for normal dispersion and dark solitons for anomalous dispersion. All other types of solitons are conditionally stable. For the Gross-Pitaevskii equation we obtain that dark solitons are always stable and bright solitons are always conditionally stable. The obtained results are verified using computer simulations.

Acknowledgments

This publication was made possible by NPRP Grants #5-674-1-114 and #6-021-1-005 from the Qatar National Research Fund (a member of the Qatar Foundation). Work at the Institute of Physics in Belgrade was supported by the Ministry of Education and Science of the Republic of Serbia under the projects OI 171006 and III 46016.

References and links

1. G. P. Agrawal, Nonlinear Fiber Optics, 4th ed. (Academic, 2007).

2. L. Ostrovskii, “Propagation of wave packets and space-time self focusing in a nonlinear medium,” Sov. Phys. JETP 24, 797–800 (1969).

3. Y. S. Kivshar and M. Peyrard, “Modulational instabilities in discrete lattices,” Phys. Rev. A 46, 3198–3205 (1992). [CrossRef] [PubMed]

4. N. Akhmediev and A. Ankiewicz, Solitons (Chapman and Hall, 1997).

5. Y. S. Kivshar and G. Agrawal, Optical Solitons, from Fibers to Photonic Crystals (Academic, 2003).

6. A. Hasegava and M. Matsumoto, Optical Solitons in Fibers (Springer, 2003). [CrossRef]

7. B. Malomed, Soliton Management in Periodic Systems (SpringerUS, 2006).

8. Z. Rapti, P. G. Kevrekidis, A. Smerzi, and A. R. Bishop, “Variational approach to the modulational instability,” Phys. Rev. E 69, 017601 (2004). [CrossRef]

9. Y. S. Kivshar and M. Salerno, “Modulational instabilities in the discrete deformable nonlinear Schrodinger equation,” Phys. Rev. E 49, 3543–3546 (1994). [CrossRef]

10. P. Shukla and J. Rasmussen, “Modulational instability of short pulses in long optical fibers,” Opt. Lett. 11, 171–173 (1986). [CrossRef] [PubMed]

11. M. Saha and A. Sarma, “Modulation instability in nonlinear metamaterials induced by cubicquintic nonlinearities and higher order dispersive effects,” Opt. Commun. 291, 321–325 (2013). [CrossRef]

12. H. Sakaguchi and Y. Kageyama, “Modulational instability and breathing motion in the two-dimensional nonlinear Schrodinger equation with a one-dimensional harmonic potential,” Phys. Rev. E 88, 053203 (2013). [CrossRef]

13. M. Potasek, “Modulation instability in an extended nonlinear Schrdinger equation,” Opt. Lett. 12, 921–923 (1987). [CrossRef] [PubMed]

14. M. Centurion, M. A. Porter, Y. Pu, P. G. Kevrekidis, D. J. Frantzeskakis, and D. Psaltis, “Modulational instability in a layered Kerr medium: theory and experiment,” Phys. Rev. Lett. 97, 234101 (2006). [CrossRef]

15. S. Wen, Y. Wang, W. Su, Y. Xiang, X. Fu, and D. Fan, “Modulation instability in nonlinear negative-index material,” Phys. Rev. E 73, 036617 (2006). [CrossRef]

16. S. Wen, Y. Xiang, W. Su, Y. Hu, X. Fu, and D. Fan, “Role of the anomalous self-steepening effect in modulation instability in negative-index material,” Opt. Express 14, 1568–1575 (2006). [CrossRef] [PubMed]

17. Y. Xiang, S. Wen, X. Dai, Z. Tang, W. Su, and D. Fan, “Modulation instability induced by nonlinear dispersion in nonlinear metamaterials,” J. Opt. Soc. Am. B 24, 3058–3063 (2007). [CrossRef]

18. Y. Xiang, X. Dai, S. Wen, and D. Fan, “Modulation instability in metamaterials with saturable nonlinearity,” J. Opt. Soc. Am. B 28, 908–916 (2011). [CrossRef]

19. X. Dai, Y. Xiang, S. Wen, and D. Fan, “Modulation instability of copropagating light beams in nonlinear metamaterials,” J. Opt. Soc. Am. B 26, 564–571 (2009). [CrossRef]

20. A. Sarma and M. Saha, “Modulational instability of coupled nonlinear field equations for pulse propagation in a negative index material embedded into a Kerr medium,” J. Opt. Soc. Am. B 28, 944–948 (2011). [CrossRef]

21. N. Lazarides and G. P. Tsironis, “Coupled nonlinear Schrödinger field equations for electromagnetic wave propagation in nonlinear left-handed materials,” Phys. Rev. E 71, 036614 (2005). [CrossRef]

22. I. Kourakis and P. K. Shukla, “Nonlinear propagation of electromagnetic waves in negative-refraction-index composite materials,” Phys. Rev. E 72, 016626 (2005). [CrossRef]

23. I. Kourakis, N. Lazarides, and G. P. Tsironis, “Self-focusing and envelope pulse generation in nonlinear magnetic metamaterials,” Phys. Rev. E 75, 067601 (2007). [CrossRef]

24. W. P. Zhong, R. H. Xie, M. Belić, N. Petrović, G. Chen, and L. Yi, “Exact spatial soliton solutions of the two-dimensional generalized nonlinear Schrdinger equation with distributed coefficients,” Phys. Rev. A 78, 023821 (2008). [CrossRef]

25. M. Belić, N. Petrović, W. P. Zhong, R. H. Xie, and G. Chen, “Analytical light bullet solutions to the generalized (3+ 1)-dimensional nonlinear Schrdinger equation,” Phys. Rev. Lett. 101, 123904 (2008). [CrossRef]

26. N. Petrović, M. Belić, W. Zhong, R. Xie, and G. Chen, “Exact spatiotemporal wave and soliton solutions to the generalized (3+ 1)-dimensional Schrdinger equation for both normal and anomalous dispersion,” Opt. Lett. 34, 1609–1611 (2009). [CrossRef]

27. N. Z. Petrović, M. Belić, and W. P. Zhong, “Spatiotemporal wave and soliton solutions to the generalized (3+1)-dimensional Gross-Pitaevskii equation,” Phys. Rev. E 81, 016610 (2010). [CrossRef]

28. N. Z. Petrović, N. B. Aleksić, A. A. Bastami, and M. R. Belić, “Analytical traveling-wave and solitary solutions to the generalized Gross-Pitaevskii equation with sinusoidal time-varying diffraction and potential,” Phys. Rev. E 83, 036609 (2011). [CrossRef]

29. Y. B. Zhou, M. L. Wang, and Y. M. Wang, “Periodic wave solutions to a coupled KdV equations with variable coefficients,” Phys. Lett. A 308, 31–36 (2003). [CrossRef]

30. Y. B. Zhou, M. L. Wang, and T. D. Miao, “The periodic wave solutions and solitary wave solutions for a class of nonlinear partial differential equations,” Phys. Lett. A 323, 77–88 (2004). [CrossRef]

31. G. P. Agrawal, “Modulation instability induced by cross-phase modulation,” Phys. Rev. Lett. 59, 880–883 (1987). [CrossRef] [PubMed]

32. I. Towers and B. A. Malomed, “Stable (2+ 1)-dimensional solitons in a layered medium with sign-alternating Kerr nonlinearity,” J. Opt. Soc. Am. 19, 537–543 (2002). [CrossRef]

33. S. A. Ponomarenko and G. P. Agrawal, “Interactions of chirped and chirp-free similaritons in optical fiber amplifiers,” Opt. Express 15, 2963–2973 (2007). [CrossRef] [PubMed]

34. C. Dai, Y. Wang, and J. Zhang, “Analytical spatiotemporal localizations for the generalized (3+1)-dimensional nonlinear Schrödinger equation,” Opt. Lett. 35, 1437–1439 (2010). [CrossRef] [PubMed]

35. F. Ndzana II, A. Mohamadou, and T. C. Kofane, “Modulational instability in the cubicquintic nonlinear Schrödinger equation through the variational approach,” Opt. Commun. 275, 421–428 (2007). [CrossRef]

36. E. Wamba, A. Mohamadou, and T. C. Kofane, “A variational approach to the modulational instability of a Bose-Einstein condensate in a parabolic trap,” J. Phys. B: At. Mol. Opt. Phys. 41, 225403 (2008). [CrossRef]

37. N. Petrović, H. Zahreddine, and M. Belić, “Exact spatiotemporal wave and soliton solutions to the generalized (3+ 1)-dimensional nonlinear Schrödinger equation with linear potential,” Phys. Scripta 83, 065001 (2011). [CrossRef]

38. Y. Castin, Coherent Atomic Matter Waves, R. Kaiser, C. Westbrook, and F. David, eds. (Les Houches Session LXXII, Springer, 2001).

Modulation stability analysis of exact multidimensional solutions to the generalized nonlinear Schrödinger equation and the Gross-Pitaevskii equation using a variational approach

Abstract

1. Introduction

2. Generalized nonlinear Schrödinger equation

3. Variational approach to the modulation stability

3.1. Case without chirp

3.2. Case with chirp

4. Numerical simulations

5. Analysis of the stability of the Gross-Pitaevskii equation

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (7)

Tables (1)

Equations (63)

Optics Express

	s	σ	SW	pert.	κ	type	stability (3D)	stability (2D)
1	1	1	S-spatial	S T	1 1	cn	CS-Conditionally Stable	CS
2	1	1	T-temporal	S	1	cn	CS	–
3	1	−1	S	S T	1 1	sn	S-stable	S
4	1	−1	T	S	1	sn	S	–
5	−1	1	S	S T	1 −1	cn	CS	CS
6	−1	1	T	S	1	sn	CS	–
7	−1	−1	S	S T	1 −1	sn	CS	S
8	−1	−1	T	S	1	cn	S	–