Instabilities of one-dimensional stationary solutions of the cubic nonlinear Schrödinger equation

Roger J. Thelwell†, John D. Carter‡, Bernard Deconinck†

†Department of Applied Mathematics ‡Mathematics Department
University of Washington  Seattle University
Seattle, WA 98195-2420  Seattle, WA 98122

Draft - November 7, 2005

Abstract

The two-dimensional cubic nonlinear Schrödinger equation admits a large family of onedimensional bounded traveling-wave solutions. All such solutions may be written in terms of an amplitude and a phase. Solutions with piecewise constant phase have been well studied previously. Some of these solutions were found to be stable with respect to one-dimensional perturbations. No solutions are stable with respect to two-dimensional perturbations. We consider stability of the larger class of solutions whose phase is dependent on the spatial dimension of the one-dimensional wave form. We study the spectral stability of such nontrivial-phase solutions numerically, using Hill’s method. We present evidence which suggests that all such nontrivial-phase solutions are unstable with respect to both one- and two-dimensional perturbations. Instability occurs in all cases: for both the elliptic and hyperbolic nonlinear Schrödinger equations, and in the focusing and defocusing cases.

PACS numbers:02.30.Jr, 02.70.Hm, 03.75.Nt, 42.65.St, 52.35.Mw, 92.10.Hm Keywords: stability, nonlinear Schrödinger equation, NLS, spectral, phase profile, elliptic, hyperbolic, focusing, defocusing, stationary

1 Introduction

The cubic nonlinear Schrödinger (NLS) equation in two spatial dimensions is given by

iψ_{t} + αψ_{xx} + βψ_{yy} + |ψ|^{2}ψ = 0.  (1)

The NLS equation is said be focusing or attractive in the x-dimension if α > 0. If α < 0, the NLS equation is said to be defocusing or repulsive in the x-dimension. Similarly, the sign of β leads to focusing or defocusing in the y-dimension. The NLS equation is called hyperbolic if αβ < 0 and elliptic if αβ > 0. The elliptic NLS equation is used to model the propagation of pulses in optical fibers [1], Langmuir waves in plasmas [2], and in the description of Bose-Einstein condensates [3, 4]. The hyperbolic NLS equation is used to describe gravity waves in deep water [5], cyclotron waves in plasmas [6], and optical waves in planar waveguides [7].

Equation (1) admits a large family of one-dimensional bounded traveling-wave solutions. All such solutions, up to Lie group symmetries [8], may be written in the form

ψ(x, t) = ϕ(x)e^{iθ(x)+iλt},  (2)

where ϕ(x) and θ(x) are real-valued functions, and λ is a real constant [9, 10]. Solutions of the form (2) are possible if

ϕ^{2}(x) = α (−2k^{2} sn^{2}(x, k) + B) ,  (3a)
θ(x) = c ∫_{0}^{x} ϕ^{−2}(ξ)dξ,  (3b)
λ = \frac{1}{2}α(3B − 2 − 2k^{2}),  (3c)
c^{2} = −\frac{α^{2}}{2} B(B − 2k^{2})(B − 2),  (3d)

where B and c are real constants. Here k ∈ [0, 1] is the elliptic modulus of the Jacobi elliptic sine function, sn (x, k). The function sn (x, k) is periodic if k ∈ [0, 1), with period given by L = 4K, where K = K(k) is defined by

K(k) = ∫_{0}^{π/2} (1 − k^{2} sin^{2} x)^{−1/2} dx,  (4)

and is known as the complete elliptic integral of the first kind. When k = 0, sn (x, 0) = sin(x) with L = 2π. As k approaches 1, sn (x, k) approaches tanh(x) and L approaches infinity [11]. Although ϕ(x) is periodic, the solution ψ(x, t) is typically not L-periodic in the x-dimension, because the periods of e^{iθ} and ϕ are typically non-commensurate.

The solution ψ is said to have trivial-phase (TP) if θ(x) is (piecewise) constant and nontrivialphase (NTP) if θ(x) is not constant. Equivalently, the solution ψ has TP if c = 0, and has NTP if c ≠ 0. For every choice of α and β, (3) specifies a two-parameter family of NLS solutions with the free parameters k and B. Without loss of generality, we choose both α and β to be ±1. The phase contribution θ(x) given in (3b) implicitly depends on α and B in both (3a) and (3d). In order for ϕ and θ to be real-valued functions, we need B ∈ [2k^{2}, 2] if α = 1 or B ≤ 0 if α = −1. Figure 1 displays the regions of (k, B)-parameter space that correspond to NTP solutions of the focusing and defocusing NLS equations. By varying B so that c → 0, θ approaches a (piecewise) constant, and the NTP solutions reduce to one of five types of TP solutions: (i) a Stokes’ plane wave, (ii) a cn-type solution, (iii) a dn-type solution, (iv) an sn-type solution, (v) a soliton-type solution. The limiting solutions correspond to the boundaries of the regions in Fig. 1. Table 1 provides the values of k and B that cause (2) to reduce to TP solutions, and also gives the explicit expression for ψ in the TP limit. Gray solitons, all of which are NTP solutions, result when k = 1 and B < 0. An overview of stationary NLS solutions is given in [9, 10]. Properties of the Jacobi elliptic functions cn, dn and sn may be found in [11].

The stability analysis of TP solutions has been well investigated, see for example [12, 13, 14, 15, 16, 17]. While some TP solutions are stable under one-dimensional perturbations (the bright soliton [18, 19] and sn-type solutions [17]), all TP solutions are known to be unstable under twodimensional perturbations [17]. The only work with results regarding the stability of some NTP solutions that we are familiar with is the work of Infeld and Ziemkiewicz [20]. In their work, an additional damping term led to stable NTP solutions. However, without this term all NTP solutions they considered were unstable; our results agree with this conclusion. The work presented here differs from [20] in that we consider the entire parameter space of NTP solutions, so that all one-dimensional stationary solutions of the NLS equation have now been investigated.

The (k, B)-parameter space of the NTP solution is two-dimensional, whereas the (k, B)-parameter space of TP solutions is essentially one-dimensional; it forms the boundary of the NTP (k, B)- parameter spaces shown in Fig. 1. Thus the TP solutions are only a co-dimension one subset of all PSfrag replacements

Stokes’ sn cn dn

soliton gray dark

k B 0 -1 2 1 PSfrag replacements

Stokes’ sn cn dn soliton gray dark k B 0 -1 2 1

(a) α = +1  (b) α = −1

Figure 1: Admissible parameter space for solutions of the form given in (2) for (a) focusing and (b) defocusing regimes. The shaded interior region corresponds to NTP solutions. The region in the defocusing regime is unbounded below.

bounded traveling-wave solutions of (1). The aim of this work is to investigate if (1) has any onedimensional traveling-wave solutions that are stable with respect to either one- or two-dimensional perturbations. Since no solutions are known that are stable with respect to two-dimensional perturbations in the TP setting, we focus on such perturbations, although one-dimensional perturbations are considered as a special case.

We investigate the spectral stability of all NTP solutions (2), for all choices of α = ±1 and β = ±1. Although all NTP solutions are found to be unstable, our investigations do produce important information about the nature of the instabilities of these NTP solutions.

2 The linear stability problem

In order to study the linear stability of NTP solutions of the NLS equation, we consider perturbed solutions of the form

ψ_{p}(x, y, t) = (ϕ(x) + ϵu(x, y, t) + iϵv(x, y, t) + O(ϵ^{2})) e^{iθ(x)+iλt},  (5)

where u(x, y, t) and v(x, y, t) are real-valued functions, ϵ is a small real parameter and ϕ(x) e^{iθ(x)+iλt} is an NTP solution of NLS. Substituting (5) in (1), linearizing and separating real and imaginary parts leads to

λu − 3γϕ^{2}u − βu_{yy} + αc^{2} \frac{1}{ϕ^{4}} u − 2αc\frac{1}{ϕ^{3}} ϕ_{x}v + 2αc\frac{1}{ϕ^{2}} v_{x} − αu_{xx} = −v_{t},  (6a)

λv − γϕ^{2}v − βv_{yy} + αc^{2} \frac{1}{ϕ^{4}} v + 2αc\frac{1}{ϕ^{3}} ϕ_{x}u − 2αc\frac{1}{ϕ^{2}} u_{x} − αv_{xx} = u_{t}.  (6b)

k value B value  Solution type  ψ
α = 1 k = 0  B ∈ [2k^{2}, 2] Stokes’ plane wave \sqrt{B} e^{iλt}
k ∈ (0, 1) B = 2k^{2}  cn-type  \sqrt{2} k cn (x, k) e^{i(2k^{2}−1)t}
k ∈ (0, 1) B = 2  dn-type  \sqrt{2} dn (x, k) e^{i(2−k^{2})t}
k = 1  B = 2  bright soliton  \sqrt{2} sech(x)e^{it}
α = −1 k = 0  B ≤ 0  Stokes’ plane wave \sqrt{B} e^{iλt}
k ∈ (0, 1) B = 0  sn-type  \sqrt{2} k sn (x, k)e^{i(1+k^{2})t}
k = 1  B = 0  dark soliton  \sqrt{2} tanh(x) e^{i2t}

Table 1: Parameter values (k, B) which reduce NTP solutions to TP solutions.

Since (6) does not depend on y or t explicitly, we assume that u(x, y, t) and v(x, y, t) have the forms

u(x, y, t) = U(x, ρ, Ω) e^{iρy+Ωt} + c.c.,  (7a)
v(x, y, t) = V (x, ρ, Ω) e^{iρy+Ωt} + c.c.,  (7b)

where ρ is a real constant, U(x) and V (x) are complex-valued functions, Ω is a complex constant and c.c. denotes complex conjugate. Notice that ρ is the transverse wavenumber of the perturbation and Ω is the exponential growth rate associated with ρ. If bounded U, V exist such that Ω has a positive real part, then the amplitudes of the perturbations grow exponentially in time and the unperturbed solution is unstable.

Upon substitution, (6) yields the spectral problem

λU − 3γϕ^{2}U + βρ^{2}U + αc^{2} \frac{1}{ϕ^{4}} U − 2αc\frac{1}{ϕ^{3}} ϕ_{x}V + 2αc\frac{1}{ϕ^{2}} V_{x} − αU_{xx} = −ΩV,  (8a)

λV − γϕ^{2}V + βρ^{2}V + αc^{2} \frac{1}{ϕ^{4}} V + 2αc\frac{1}{ϕ^{3}} ϕ_{x}U − 2αc\frac{1}{ϕ^{2}} U_{x} − αV_{xx} = ΩU.  (8b)

If c = 0, then (8) reduces to the stability problem for TP solutions. This case is examined in [13, 20, 21, 22, 23, 17] and others. Using the linear system (8), we are able to investigate the stability of all NTP solution numerically. We only consider the stability of NTP solutions in this paper. The reader may wish to consult some of the above mentioned references for the stability analysis of their limiting special cases.

3 The numerical investigation of spectral stability: Hill’s method

The main difficulty for the numerical investigation of (8) is the size of the parameter space involved. For every choice of the equations parameters α, β and solution parameter pairs (k, B), the spectrum of (8) needs to be computed for a range of ρ values in order to determine stability or to analyze any instabilities. An efficient numerical method is necessary. Hill’s method allows for the systematic and efficient exploration of the large phase space encountered here, due to its exponential convergence [24].

To apply Hill’s method, Fourier expansions are needed for all coefficient functions of (8). Using the complex Fourier form, we have

\begin{matrix}ϕ^{2}(x) = ∑_{n=−∞}^{∞} Q_{n} e^{i2nπx/L}, & ϕ^{−2}(x) = ∑_{n=−∞}^{∞} R_{n} e^{i2nπx/L}, \\ ϕ^{−4}(x) = ∑_{n=−∞}^{∞} S_{n} e^{i2nπx/L}, & ϕ^{−3}(x)ϕ′(x) = ∑_{n=−∞}^{∞} T_{n} e^{i2nπx/L}, \end{matrix}  (9)

where Q_{n}, R_{n}, S_{n} and T_{n} are Fourier coefficients. Note that ϕ^{2}(x) has period L/2 and that ϕ(x) is never zero except in some TP limit cases.

The periodicity of the coefficient functions in (8) allows us to decompose the eigenfunction components U and V of (8) in a Fourier-Floquet form

U(x) = e^{iμx} ∑_{n=−∞}^{∞} U_{n} e^{inπx/L} and V (x) = e^{iμx} ∑_{n=−∞}^{∞} V_{n} e^{inπx/L},  (10)

where μ ∈ [\frac{−π}{K} , \frac{π}{K}). The form of U and V in (10) follows from Floquet’s theorem and the observation that eigenfunctions are bounded, by definition. This decomposition has the benefit of admitting both L-periodic and L-anti-periodic eigenfunctions when μ = 0. Recall that ψ is typically only quasiperiodic. Allowing μ to be different from 0 gives rise to solutions that are either quasiperiodic or have period greater than 2L.

Substitution of (9) and (10) into (8) and equating Fourier coefficients allows us to write equations for U_{n} and V_{n} as a coupled bi-infinite system of difference equations given by

− (λ + βρ^{2} − α (iμ + \frac{inπ}{L} )^{2}) U_{n} + 3γ ∑_{m=−∞}^{∞} Q_{\frac{n−m}{2}} U_{m} − αc^{2} ∑_{m=−∞}^{∞} S_{\frac{n−m}{2}} U_{m}

+ 2αc ∑_{m=−∞}^{∞} T_{\frac{n−m}{2}} V_{m} − 2αc (iμ + \frac{inπ}{L} ) ∑_{m=−∞}^{∞} R_{\frac{n−m}{2}} V_{m} = ΩV_{n}, (11a)

(λ + βρ^{2} − α (iμ + \frac{inπ}{L} )^{2}) V_{n} − γ ∑_{m=−∞}^{∞} Q_{\frac{n−m}{2}} V_{m} + αc^{2} ∑_{m=−∞}^{∞} S_{\frac{n−m}{2}} V_{m}

+ 2αc ∑_{m=−∞}^{∞} T_{\frac{n−m}{2}} U_{m} − 2αc (iμ + \frac{inπ}{L} ) ∑_{m=−∞}^{∞} R_{\frac{n−m}{2}} U_{m} = ΩU_{n}, (11b)

for all integers n. Here Q_{\frac{n−m}{2}} = 0 if \frac{n−m}{2} ∉ ℤ, with R_{(·)}, S_{(·)} and T_{(·)} defined similarly. The system of equations (11) is equivalent to the original system (8).

Remarks

• In practice, a pre-multiplication of the linear system by ϕ^{4} allows for the exact Fourier series expansion of ϕ^{2}, ϕ^{4} and ϕ^{6} to be used. This follows from the differential equations for sn (x, k) and Jacobi’s series expansion of sn^{2}(x, k) [25]. This pre-multiplication transforms the original eigenvalue problem into a generalized eigenvalue problem [26].

• Note that Hill’s methods enables one to compute the spectrum of a linear operator with periodic coefficients. Despite the fact that the solution (2) is typically quasiperiodic, the coefficient functions of the linear stability problem (8) are always periodic.

Parameter  Description  Value

     k  Elliptic modulus  linspace(0,1,65)

B Shift

For α = −1 : −logspace(−8, 0, 65)
For α = 1 : (2k^{2} + logspace(−8, 0, 65)) ∩ (2k^{2}, 2)

N Fourier cutoff

For α = −1 : 15 + ceil(5k^{5})
For α = 1 : 10 + ceil(25k^{10})
ρ  Perturbation wavenumber linspace(0,4,65)
μ  Floquet parameter  linspace(−\frac{π}{K} , \frac{π}{K} , 21)

Table 2: Parameter values and ranges used in numerical experiments. Only perturbations of NTP solutions are considered.

4 Numerical experiments

By choosing a finite number of Fourier modes, the exact bi-infinite system (11) is truncated. We explicitly construct and compute approximations to the spectral elements (i.e. eigenvalues or elements of the continuous spectrum) of (8) by finding the eigenvalues of the truncation of (11). We consider all four cases individually: (I) focusing in both x and y (α = β = 1), (II) focusing in x and defocusing in y (α = −β = 1), (III) defocusing in x and focusing in y (−α = β = 1) and finally, (IV) defocusing in both x and y (−α = −β = 1).

In each case, a large number of parameter values in the two-dimensional parameter space shown in Fig. 1 was explored numerically. The (k, B)-parameter values considered correspond to NTP solutions, and do not include TP solutions or gray solitons. Approximately 5.2 million generalized eigenvalue problems were considered, the size of each determined by the cutoff mode N of the underlying Fourier series. A truncation to N positive Fourier modes reduces the exact bi-infinite system (11) to an (4N + 2)-dimensional approximate problem. For several choices of k and B, a value of N = N(k, B) was chosen to ensure that the resulting eigenvalues had converged to within a measured tolerance. A simple polynomial was used to fit this data. This information, and details related to other parameter ranges used in the experiments, are included in Table 2. In the table, k is the elliptic modulus and B is the offset parameter (as in (3)), (4N + 2) is the dimension used to approximate (8), ρ is the wavenumber of the perturbation in the y-dimension, and μ is the Floquet exponent. Lastly, linspace(a, b, m) is a linearly spaced vector from a to b of length m, logspace(a, b, m) is a logarithmically spaced vector from 10^{a} to 10^{b} of length m and ceil(x) is the smallest integer greater than or equal to x.

5 Observations from the numerical investigation

First and foremost, it should be stated that none of the NTP solutions considered here were found to be spectrally stable under one-dimensional (ρ = 0) or two-dimensional (ρ ≠ 0) perturbations. This establishes, at least numerically, that all one-dimensional traveling-wave solutions of NLS of the form (2) are spectrally unstable with respect to either one- or two-dimensional perturbations. At this point, it remains to investigate the nature of the instabilities and their corresponding growth rates, so as to better understand the dynamics of this important class of solutions of the NLS equation.

Using Hill’s method we numerically considered the instabilities due to perturbations with wavenumber ρ ∈ [0, 4]. Note that ρ = 0 corresponds to one-dimensional perturbations. For each NLS equation (i.e. for each choice of α, β = ±1), and solution (i.e. each parameter pair (k, B)), and for each perturbation of wavenumber ρ, an equally-spaced sequence of Floquet parameters μ was chosen from the interval [−\frac{π}{K} , \frac{π}{K}]. The generalized eigenvalues and eigenvectors were computed for the matrix that results from a truncation of (11). The generalized eigenvalues are approximations of spectral elements of (8), and an approximation of the corresponding eigenfunctions may be reconstructed from the eigenvectors.

Since a single eigenvalue with positive real part leads to instability of the system, the approximate eigenvalue with largest real part over all choices of μ was recorded for each (k, B, ρ) triplet. That is, we compute

Ω_{growth}(k, B, ρ) = max_{μ∈[−π/K,π/K]} Re [Ω(k, B, ρ, μ)] ,  (12)

which we call the (most unstable) growth rate. The value Ω_{growth} represents the largest exponential growth rate a given NTP solution with parameters (k, B) will experience when perturbed with transverse wavenumber ρ. It also allows us to determine the perturbation to which the NTP solution is spectrally the most unstable. We reduce the dimension further by computing the largest such growth rate over all sampled perturbation wavenumbers ρ. This quantity,

Ω_{max}(k, B) = max_{ρ∈[0,4]} Ω_{growth}(k, B, ρ),  (13)

the maximal growth rate over all ρ, is plotted in the first columns of Figs. 2 and 3. The value Ω_{max} represents the maximal exponential growth rate that a solution with parameters (k, B) can undergo in the range examined. We also recorded the minimum growth rate over all ρ,

Ω_{min}(k, B) = min_{ρ∈[0,4]} Ω_{growth}(k, B, ρ),  (14)

to verify that all solutions are unstable with respect to every sampled perturbation.

Every point plotted in Figs. 2 and 3 corresponds to an NTP solution for which we considered the linear stability analysis, and the boundaries in the figures are the boundaries of the regions represented in Fig. 1, corresponding to limiting TP solutions. Fig. 2 corresponds to the x-focusing (α = 1) parameter range (k, B) = (0, 1) × (2k^{2}, 2) in the α = 1 case of Fig. 1a. The one-to-one transform T_{f}(B) = (B − 2k^{2})/(2 − 2k^{2}) is used to normalize the range of B. This maps the interval [2k^{2}, 2] to [0, 1]. Figure 3 corresponds to the x-defocusing (α = −1) parameter range of (k, B) = (0, 1) × (−1, 0) shown in Fig. 1b. The transform T_{d}(B) = −B is used in Fig. 3. A log_{10} scale is used in the vertical dimension of Figs. 2 and 3. This causes the panels of Fig. 2 to become increasing sparse in their lower right corners. The right-hand panels of Figs. 2 and 3 indicate the wavenumber ρ that leads to maximal growth shown in the left-hand panels. Recall that our computations were limited to ρ ∈ [0, 4].

5.1 Case I: Elliptic setting with α = β = 1

Panels Ia and Ib of Fig. 2 summarize some properties of the computed instabilities in the case of focusing in both the x- and y-dimensions. The lower boundary of the plot corresponds to B = 2k^{2} + (10^{−8}), and is therefore only slightly away (in the parameter space of B) from a cn-type solution. The upper boundary is close to dn-type solutions, with B = 1.99. The left boundary of the plots, where k = 0.01, represents a region in parameter space near to Stokes’ wave solutions. The entire right-hand boundary, where k = 0.99, is near to the bright soliton limit case which occurs at (k, B) = (1, 2).

A distinct ridge of large instability is noticeable in the plot of Ω_{max} in panel Ia of Fig. 2. The ridge appears to begin near the zero solution at (k, B) = (0, 0), and remains close to the cn limit boundary (within approximately .02 units, remembering the log_{10} scaling) as k increases. Moving away from the cn boundary results in the rapid increase of Ω_{max}. Movement away from the dn boundary results in a slower increase in the value of Ω_{max}, as does moving away from the Stokes’ wave boundary for B greater than approximately 0.001. The maximum value of Ω_{max} over the sampled (k, B) space, given by R_{max} = 5.666, is reached at (k, B) = (0.99, 1.98). This growth rate should be compared to the maximal growth rate of the corresponding TP case [17] which is R_{max} = 1. The minimum (R_{min} = 0.015693) occurs for (k, B) = (0.01, 0.01) for ρ = 4. Note that this is on the boundary of the computational domain and that at (k, B) = (0, 0), corresponding to the spectrally stable ψ ≡ 0 solution, R_{min} = 0.

In panel Ib, the wavelength corresponding to the maximal growth of Ia is given. In this case, the maximum instability occurs for the shortest wavelength samples, ρ = 4. This indicates that there is a strong short-wavelength instability.

5.2 Case II: Hyperbolic setting with α = −β = 1

Panels IIa and IIb of Fig. 2 summarize some properties of the computed instabilities in the case of focusing in the x-dimension and defocusing in the y-dimension. The lower boundary of the plot corresponds to B = 2k^{2} + (10^{−8}), and is therefore only slightly away (in the parameter space of B) from a cn-type solution. The upper boundary is close to dn-type solutions, with B = 1.99. The left boundary of the plots, where k = 0.01, represents a region in parameter space near to Stokes’ wave solutions. The entire right-hand boundary, where k = 0.99, is near to the bright soliton limit case which occurs at (k, B) = (1, 2).

As in Case I, a ridge of large growth rate is noticeable in the growth plot shown in panel IIa. The ridge appears to begin near the zero solution at (k, B) = (0, 0), and remains close to the cn-type limit boundary (within approximately .02 units, remembering the log_{10} scaling) as k increases. This ridge has a local minimum near k = 0.7 and increases to a global (over all admissible (k, B)-parameter space) maximum at k = 0.96. As in the setting above, moving away from the cn-type boundary results in a rapid increase of Ω_{max}. Moving away from the dn boundary results in a slower increase in the value of Ω_{max}. For B > 0.001, moving away from the boundary result in a similar slow increase in Ω_{max}. For k > 0.96, it appears that the limiting value of Ω_{max} is consistent with the bright soliton results of [18, 19]. The maximum (R_{max} = 6.1141) and minimum (R_{min} = 0.012535) growth rates span a slightly larger range than the similar values in Fig. Ia. These occur at (k, B) = (0.01, 0.01) and (k, B) = (0.96, 1.98), respectively. The maximal growth rate should be compared to the maximal growth rate of the corresponding TP case [17] which is R_{max} = 1.

In panel IIb, the wavelength corresponding to the maximal growth R_{max} of IIa is given. In this case, the maximum instability of R_{max} occurs for ρ = 3.375. The surface shown in IIb appears to be more smooth than the surface of Ib.

5.3 Case III: Hyperbolic setting with −α = β = 1

Panels IIIa and IIIb of Fig. 3 summarize some properties of the computed instabilities in the case of defocusing in the x-dimension and focusing in the y-dimension. The lower limit of the plot corresponds to B = −10^{−8}, and so is just slightly away from the sn-type solution. The left boundary of the plots, where k = 0.01, represents a region in parameter space near Stokes’ wave solutions, while k = 0.99 on the right boundary is near to the gray soliton limit.

A distinct ridge of large instability is noticeable in the growth plot displayed in panel IIIa. The ridge appears to begin near the zero solution limit at (k, B) = (0, 0), and remains close to the sn limit boundary (within approximately .02 units, remembering the log_{10} scaling) as k increases. It quickly reaches the global maximum (over all admissible (k, B)-parameter space) of R_{max} = 7.6375 at k = 0.02 and B = −0.0001. The ridge then decreases in amplitude as k increases towards 1. Moving away from the sn-type boundary results in a rapid increase of Ω_{max}. Moving away from the dn-type boundary results in a slower increase in the value of Ω_{max}. Similarly, the increase is slower when moving away from the Stokes’ wave limit for B < −0.001. The maximum exponential growth rate, R_{max} = 7.6375, occurs for (k, B) = (0.02, −0.00001). This growth rate should be compared to the maximal growth rate of the corresponding TP case [27] which is R_{max} = 1. The minimum exponential growth, R_{min} = 0.015578 is found at (k, B) = (0.01, −0.9). Both the maximum and minimum are located near the Stokes’ wave boundary. By restricting ρ = 0 and allowing B to approach zero, Ω_{max} → 0, and the one-dimensional stability result of the sn-type TP solution of [17] is recovered.

The plot IIIb indicates short-wave perturbations lead to large values of Ω_{max}. The largest growth occurs for a perturbation with wavenumber of ρ = 3.625.

5.4 Case IV: Elliptic setting with −α = −β = 1

Panels IVa and IVb of Fig. 3 summarize some properties of the computed instabilities in the case of defocusing in both the x- and y-dimensions. The lower limit of the plot corresponds to B = −10^{−8}, and so is just slightly away from the sn-type solution. The left boundary of the plots, where k = 0.01, represents a region in parameter space near to Stokes’ wave solutions, while k = 0.99 on the right boundary is near to the gray soliton limit.

A distinct ridge of large instability is noticeable in the growth plot shown in panel IVa. The ridge appears to begin near the trivial limit k = 0 and B = 0, and remains close to the sn limit boundary (within approximately .02 units, remembering the log_{10} scaling) as k increases, to reach a global maximum at k = 0.02 and B = −0.0001. The ridge then appears to decrease in amplitude as k increases towards 1. As in Case III, moving away from the sn-type boundary results in a rapid increase of Ω_{max}. Moving away from the dn-type boundary results in a slower increase in the value of Ω_{max}. The same is true when moving away from the Stokes’ boundary, when B is less than approximately -0.001. The maximum exponential growth rate, R_{max} = 7.6456, and the minimum, R_{min} = 0.0001556, span a slightly larger range of values than do the values of Ω_{max} in panel IIIa. The maximum and minimum values are obtained at (k, B) = (0.01, −0.00009) and (k, B) = (0.01, −1), respectively. Both are located near the Stokes’ wave boundary. The maximal growth rate should be compared to the maximal growth rate of the corresponding TP case [27] which is R_{max} ≈ 0.26. As in Case III, restricting ρ = 0 and allowing B to approach zero results in Ω_{max} → 0, and the one-dimensional stability result of [17] for the sn-type TP solution is recovered.

In panel IVb, wavenumbers corresponding to Ω_{max} of IVa are given. It appears that a majority of the large values of Ω_{max} are attributable to small-ρ (long-wave) perturbations. In fact, the largest growth occurs for ρ = 0, the one-dimensional perturbation. This should be contrasted with the three previous cases, where short-wavelength two-dimensional perturbations with wavenumber ρ > 3 were associated with the largest Ω_{max} values.

6 Summary

In this paper, we considered the spectral instability of one-dimensional traveling-wave nontrivialphase (NTP) solutions of the cubic nonlinear Schrödinger equation. Such solutions are expressed in terms of Jacobi elliptic functions. An exact spectral form of the linearized operator is truncated and used to construct an associated generalized eigenvalue problem. The positive real parts of the resulting eigenvalues were used to determine that there are no stable NTP solutions of the NLS equation.

Numerical results indicate a well-defined ridge of large growth rate located in the (k, B)- parameter region associated with nontrivial-phase solutions. This implies that the most unstable NTP solutions are more unstable than any TP solution, in the sense that they exhibit larger exponential growth rates. Further, for all cases the exponential growth rate Ω_{max} increases when moving away from the limiting TP solutions. This divergence is gradual in some cases, but very sharp in other cases, as discussed above.

In summary, numerical evidence suggests that all bounded, nontrivial-phase one-dimensional traveling-wave solutions of the cubic NLS equation are unstable with respect both one-dimensional and two-dimensional perturbations.

We gratefully acknowledge support from the National Science Foundation: NSF-DMS-FRG:0139847 (RJT), NSF-DMS-FRG:0139771 (JDC) and NSF-DMS-FRG:0351466 (BD).

References

[1] S. V. Manakov. On the theory of two-dimensional stationary self-focusing of electromagnetic waves. Sov. Phys. JETP, 38:248, 1974.

[2] H. L. Pecseli. Solitons and weakly nonlinear waves in plasmas. IEEE Transactions Plasma Science, 13(2):53, 1985.

[3] E. P. Gross. Hydrodynamics of a superfluid condensate. J. Math Phys., 4(2):195, 1963.

[4] L. P. Pitaevskii. Vortex lines in an imperfect Bose gas. Sov. Phys. JETP, 13(2):451, 1961.

[5] D. J. Benney and G. J. Roskes. Wave instabilities. Stud. Appl. Math, 377(48), 1969.

[6] J. R. Myra and C. S. Liu. Self-moudulation of ion Bernstein waves. Phys. Fluids, 225(23), 1980.

[7] S.-P. Gorza, N. Roig, Ph. Emplit, and M. Haelterman. Snake instability of a spatiotemporal bright soliton strip. Phys. Rev. Letters, 92(8), 2004.

[8] P. L. Sulem and C. Sulem. Nonlinear Schrödinger Equations: Self-focusing and Wave Collapse. Springer-Verlag, New York, NY, 1999.

[9] L. D. Carr, C. W. Clark, and W. P. Reinhardt. Stationary solutions of the one-dimensional nonlinear Schrödinger equation. I. Case of repulsive nonlinearity. Phys. Rev.A, 62:63610, 2000.

[10] L. D. Carr, C. W. Clark, and W. P. Reinhardt. Stationary solutions of the one-dimensional nonlinear Schrödinger equation. II. Case of attractive nonlinearity. Phys. Rev.A, 62:63611, 2000.

[11] P. F. Byrd and M. D. Friedman. Handbook of Elliptic Integrals for Engineers and Physicists. Springer-Verlag, Berlin, 1954.

[12] V. E. Zakharov and A. M. Rubenchik. Instability of waveguides and soliton in nonlinear media. Sov. Phys. JETP, 38(3):494–500, 1974.

[13] D. U. Martin, H. C. Yuen, and P. G. Saffman. Stability of plane wave solutions of the twospace-dimensional nonlinear Schrödinger equation. Wave Motion, 2:215–229, 1980.

[14] E. A. Kuznetsov, A. M. Rubenchik, and V. E. Zakharov. Soliton stability in plasmas and hydrodynamics. Phys. Rep., 142:103–165, 1986.

[15] S. E. Fil’chenkov, G. M. Fraiman, and A. D. Yunakovskii. Instability of periodic solutions of the nonlinear Schrödinger equation. Sov. J. Plasma Phys., 18(8):961–966, 1987.

[16] K. Rypdal and J. J. Rasmussen. Stability of solitary structures in the nonlinear Schrödinger equation. Phys. Scripta, 40:192, 1989.

[17] J. D. Carter and B. Deconinck. Stability of trivial-phase solutions of the two-dimensional cubic nonlinear Schrödinger equation. Submitted for publication, 2005.

[18] Y. S. Kivshar and D. E. Pelinovsky. Self-focusing and transverse instabilities of solitary waves. Phys. Rep., 331(4):118–195, 2000.

[19] B. Deconinck, D. E. Pelinovsky, and J. D. Carter. Transverse instabilities of deep-water solitary waves. Submitted for publication, 2005.

[20] E. Infeld and J. Ziemkiewicz. Stability of complex solutions of the nonlinear Schrödinger equation. Acta Phys. Pol., A59(3):255–275, 1981.

[21] V. A. Aleshkevich, A. A. Egorov, Y. V. Kartashov, V. A. Vysloukh, and A. S. Zelenina. Stability of spatiotemporal cnoidal waves in cubic nonlinear media. Phys. Rev. E, 67:066605, 2003.

[22] Y. V. Kartashov, V. A. Aleshkevich, V. A. Vysloukh, A. A. Egorov, and A. S. Zelenina. Transverse modulational instability of (2+1)-dimensional cnoidal waves in media with cubic nonlinearity. J. Opt. Soc. Am. B., 20(6):1273–1284, 2003.

[23] J. D. Carter and H. Segur. Instability in the two-dimensional cubic nonlinear Schrödinger equation. Phys. Rev. E., 68(4):045601, 2003.

[24] B. Deconinck and J. N. Kutz. Computing spectra of linear operators using Hill’s method. Submitted for publication, 2005.

[25] C. G. Jacobi. Fundamenta nova theoriae functionum ellipticarum. Königsberg, 1829.

[26] G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore, MD, 1996.

[27] A. Davey and K. Stewartson. On three-dimensional packets of surface waves. Proc. R. Soc. London Ser. A, 338:101–110, 1974.

Ia. Ω_{max} with (α, β) = (1, +1)  Ib. Corresponding ρ, (α, β) = (1, +1)

PSfrag replacements log_{10} T_{f}(B)

PSfrag replacements

log_{10} T_{f}(B)

k  k
IIa. Ω_{max} with (α, β) = (1, −1)  IIb. Corresponding ρ, (α, β) = (1, −1)

PSfrag replacements log_{10} T_{f}(B)

PSfrag replacements

log_{10} T_{f}(B)

k  k

Figure 2: Focusing in the x-dimension. The first column contains surface plots of Ω_{max}(k, B) and the second column contains the surface plots of the maximizing wavenumber ρ(k, B). R_{max} = max_{k,B} Ω_{max} and R_{min} = min_{k,B} Ω_{min}. White space corresponds to (k, B)-parameter space that was not sampled.

IIIa. Ω_{max} with (α, β) = (−1, +1)  IIIb. Corresponding ρ, (α, β) = (−1, +1)

PSfrag replacements log_{10}(−B)

PSfrag replacements

log_{10}(−B)

k  k
IVa. Ω_{max} with (α, β) = (−1, −1)  IVb. Corresponding ρ, (α, β) = (−1, −1)

PSfrag replacements log_{10}(−B)

PSfrag replacements

log_{10}(−B)

k  k

Figure 3: Defocusing in the x-dimension. The first column contains surface plots of Ω_{max}(k, B) and the second column contains surface plots of the maximizing wavenumber ρ(k, B). R_{max} = max_{k,B} Ω_{max} and R_{min} = min_{k,B} Ω_{min}.