\documentclass[11pt]{article}
%\documentclass[amsmath,amssymb]{revtex4}
\usepackage{fullpage}
\usepackage{psfrag}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{epsfig} 
%\usepackage{draftcopy}
\usepackage{comment} 
\usepackage{float}     % with epsfig
\restylefloat{figure}  % with epsfig

\input{/home/roger1/TeX/mymacros}

\title{Instabilities of one-dimensional stationary solutions of the cubic nonlinear Schr\"{o}dinger equation}
\date{Draft - \today}
%\author{John D. Carter}
%\email{carterj1@seattleu.edu}
% \affiliation{Seattle University\\Mathematics Department\\Seattle, WA 98122}
%\author{Bernard Deconinck}
%\email{bernard@amath.washington.edu}
%\affiliation{University of Washington\\Department of Applied Mathematics\\Seattle, WA 98195}
%\author{Harvey Segur}
%\email{segur@colorado.edu}
%\affiliation{University of Colorado at Boulder\\Department of Applied Mathematics\\Boulder, CO 80309}
%\author{Roger Thelwell}
%\email{thelwell@amath.washington.edu}
%\affiliation{University of Washington\\Department of Applied Mathematics\\Seattle, WA 98195}

\begin{document}

\maketitle
%\tableofcontents

\begin{abstract}
We consider the two-dimensional cubic nonlinear Schr\"{o}dinger (NLS) equation,
which admits a large family of one-dimensional bounded traveling-wave 
solutions.  All such stationary solutions may be written in terms of an 
amplitude and a phase.  Solutions with piecewise constant phase have been 
well studied~\cite{zr,rr2,MYS,infeldz,JCBDtp}, and no stable stationary 
solutions exist.  Here we 
consider stability of the larger class of solutions whose phase is 
dependent on the spatial dimension of the one-dimensional wave form.
These solutions are said to have nontrivial-phase (NTP).  
We study the spectral stability of such NTP solutions numerically, 
using (a spectrally based) Hill's method. We present evidence which 
suggests that all such NTP solutions are 
unstable with respect to both one-dimensional and transverse perturbations.  
Instability occurs in all cases: for both the elliptic and hyperbolic NLS equations, 
and in the focusing and defocusing case.  
\end{abstract}

\section{Introduction}
\label{S:Intro}

The cubic nonlinear Schr\"{o}dinger (NLS) equation in two spatial dimensions is
given by
\beq
\label{E:NLS}
i\psi_t + \alpha \psi_{xx} + \beta \psi_{yy} + |\psi|^2 \psi = 0.
\eeq
The NLS equation is said be {\em focusing} or
{\em attractive} in the $x$-dimension if $\alpha >0$.  If $\alpha < 0$,
NLS is said to be {\em defocusing} or {\em repulsive} in the 
$x$-dimension.  Similarly, the sign of $\beta$ leads to focusing 
or defocusing in the $y$-dimension~\cite{sulemsulem}.
The NLS equation is called {\em hyperbolic} if 
$\alpha \beta <0$ and {\em elliptic} if $\alpha \beta > 0$. 
Without loss of generality, both $\alpha$ and $\beta$ are chosen to be $\pm 1$. 

Equation \eqref{E:NLS} admits a large family of one-dimensional bounded traveling-wave solutions. All such solutions, modulo Lie group symmetries, may be written in the form~\cite{sulemsulem}
\beq
\label{E:soln_form}
\psi(x,t)=\phi(x) e^{i\theta(x) + i \lambda t},
\eeq
where $\phi(x)$ and $\theta(x)$ are real-valued functions, and $\lambda$ is a
real constant.  Solutions of the form \eqref{E:soln_form} are possible if
\begin{subequations}
\label{E:defs}
\begin{eqnarray}
\phi^2(x)&=& \disp{\alpha \left(-2 k^2 \sn^2(x,k)+B \right),}\label{E:phi_def} \\
\theta(x)&=& \disp{c\int_0^x\phi^{-2}(\xi)d\xi,}  \label{E:theta_def} \\
\lambda&=& \disp{\frac{1}{2}\alpha (3B-2(1+k^2)),} \label{E:omega_def} \\
c^2&=& \disp{-\frac{\alpha^2}{2} B(B-2k^2)(B-2),} \label{E:c_def}
\end{eqnarray}
\end{subequations}
where $c$ is a real constant.
Here $k \in [0,1]$ is the elliptic modulus of the 
Jacobi elliptic sine function, $\sn(x,k)$.   
The function $\sn(x,k)$ is periodic if $k\in[0,1)$, with 
period given by $L = 4K$, where $K=K(k)$ is defined by 
\begin{equation}
K(k)=\int_{0}^{\pi/2}\left(1-k^2\sin^2 x\right)^{-1/2} \, dx,
\end{equation}
which is the complete elliptic integral of the first kind.
When $k=0$, $\sn(x,0) = \sin(x)$ with $L = 2 \pi$.  As 
$k$ approaches 1, $\sn(x,k)$ limits to $\tanh(x)$ and $L$ approaches 
infinity~\cite{byrd}.

The solution $\psi$ is said to have 
{\em trivial-phase} (TP) if $\theta(x)$ is (piecewise) constant and  
{\em nontrivial-phase} (NTP) if $\theta(x)$ is not constant. Equivalently, the solution $\psi$ has TP if $c=0$, and has NTP if $c \ne 0.$
For every choice of $\alpha$ and $\beta$,~\eqref{E:defs} specifies a 
two-parameter family of solutions in the free parameters $k$ and $B$.  
Without loss of generality, both $\alpha$ and $\beta$ are chosen to be $\pm 1$.   
The phase contribution $\theta(x)$ of \eqref{E:theta_def} 
implicitly depends on $\alpha$ and $B$ in both \eqref{E:phi_def} 
and \eqref{E:c_def}. In order for $\phi$ and $\theta$ to be real-valued 
functions, we need $B \in  [2 k^2 , 2]$ if $\alpha =1$ or
$B \le 0$ if $\alpha = -1$. 
\begin{figure}
 \psfrag{Stokes}{Stokes}
 \psfrag{sn}{sn}
 \psfrag{cn}{cn}
 \psfrag{dn}{dn}
 \psfrag{soliton}{soliton}
 \psfrag{grey}{gray}
 \psfrag{dark}{dark}
 \psfrag{k}{k}
 \psfrag{B}{B}
 \psfrag{0}{0}
 \psfrag{-1}{-1}
 \psfrag{2}{2}
 \psfrag{1}{1}
 \bc
  \begin{tabular}{cc}
   \epsfig{file=./pics/focus_full.eps, width = 7 cm} &
   \epsfig{file=./pics/defocus_full.eps, width = 8 cm} \\
   (a) &  (b)
  \end{tabular}
  \caption{Admissible parameter space  for solutions of the form given in \eqref{E:soln_form} for (a) focusing ($\alpha = 1$) and 
(b) defocusing ($\alpha = -1$) regimes.  The gray interior region corresponds to NTP solutions.\label{F:parameter_space}}
 \ec
\end{figure}
Figure \ref{F:parameter_space} represents the $(k,B)$-parameter 
space corresponding to NTP solutions of the NLS equation.  
As $c \to 0$ 
, the phase $\theta$ approaches a (piecewise) constant, and the NTP solutions reduce to 
one of five simpler types of TP solutions: 
({\em i}) a Stokes' plane wave, ({\em ii}) a cn-type solution,
({\em iii}) a dn-type solution, ({\em iv}) an sn-type solution,   
({\em v}) a soliton-type solution.
The limiting solutions correspond to the boundaries of the regions in Fig.~
\ref{F:parameter_space}.  
Table~\ref{T:tp_limit} provides the values of $k$ and $B$ that 
cause~\eqref{E:soln_form} to limit to the simpler solutions, and also 
gives the explicit expression of $\psi$ for all but one.  A cumbersome expression in the case of the gray soliton has been omitted.
An overview of these NLS solutions is given in~\cite{carr1}.  
Details of the Jacobi elliptic functions sn, cn and dn may be found 
in~\cite{byrd}. 

\begin{table}[t]
\bc
\begin{tabular}{|c|ll|c|l|}
\hline
i$\alpha = 1$     & $k = 0$ &  $B\in[2k^2,2]$ & Stokes' plane wave & $\sqrt{B} \, e^{i(qx + \lambda t)}$\\
%\cline{2-5}
     &  $k \in(0,1) $ &  $B=2k^2$ & cn-type & $\sqrt{2} \, \cn(x,k) \, e^{i(2k^2 - 1)t}$ \\
     &  $k \in(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} \, \, \mbox{sech}(x) e^{it}$ \\
\hline
\hline
$\alpha = -1$ & $k = 0$ &  $B\le 0$ & Stokes' plane wave & $\sqrt{B} \, e^{i(qx + \lambda t)}$\\
     & $k \in(0,1) $ &  $B=0$ & sn-type & $\sqrt{2} \, \sn(x,k) e^{i(1+k^2)t}$\\
& $k = 1$ &  $B = 0$ & dark soliton & $\sqrt{2} \, \tanh(x) \, e^{i2t}$ \\
     & $k = 1$ &  $B < 0 $ & gray soliton & not incl.  \\
\hline
\end{tabular}
\caption{ Parameter values $(k,B)$ which reduce the NTP solution form \eqref{E:soln_form} to simpler TP solutions.  Here $q$ is an arbitrary wavenumber.\label{T:tp_limit}}
\ec
\end{table}

While both TP and NTP solutions are of interest,  
stability of TP solutions is well understood:
for example, see 
see~\cite{zr,KuzRubZak,kivpel,FilFraYun,rr2,MYS,JCBDsoliton,JCBDtp}. 
\begin{comment}
,although much of this attention has focused on the stability of 
solitary wave solutions. i
\end{comment}
All TP solutions are known to be unstable, see~/cite{JCBDtp} and the references therein
We know of only Infeld and 
Ziemkiewicz's paper~\cite{infeldz} which considers the stability of NTP 
solutions of the NLS equation.
The $(k,B)$-parameter space of TP solutions
is essentially one-dimensional; it forms the boundary of the NTP parameter 
spaces shown in  Fig.~\ref{F:parameter_space}.  The dimension of the NTP 
solution space is
fully two dimensional. We had hoped that this much larger parameter 
space would yield NTP solutions that would be stable under one-dimensional and/or transverse perturbations. 
The more difficult linear stability problem and the large parameter space that needs to 
be explored makes the NTP setting more complex than the TP setting.
In this paper, we investigate the spectral stability
of all NTP solutions \eqref{E:soln_form}, for all choices of $\alpha=\pm1$ and $\beta=\pm1$.  That is, we compute the eigenvalues of an approximate spectral problem in order to identify possible growing modes in the linear stability problem. A similar treatment of all TP solutions is found in~\cite{JCBDtp}.
No spectrally stable NTP solutions were found, but the search did provide important stability information.

\section{The linear stability problem}
\label{S:Linearized_NTP}

In order to study the linear stability of NTP solutions
of the NLS equation, we consider perturbations of the form
\begin{equation}
\psi_{_p}(x,y,t)=(\phi(x)+\epsilon u(x,y,t)+i \epsilon
v(x,y,t) + \mathcal{O}(\eps^2))\mbox{e}^{i\theta(x) + i\lambda t},
\label{E:npert}
\end{equation}
where $u(x,y,t)$ and $v(x,y,t)$ are real-valued functions, $\epsilon$ is a
small real parameter and $\phi(x)\mbox{e}^{i\theta(x)+i\lambda t}$ is
a NTP solution of NLS.  
Substituting (\ref{E:npert})
in (\ref{E:NLS}), linearizing and separating real and imaginary
parts leads to
\begin{subequations}
\begin{equation}
\lambda u-3\gamma\phi^2u-\beta
u_{yy}+ \alpha c^2 \frac{1}{\phi^4}u-2\alpha c\frac{1}{\phi^{3}}\phi_x v+2\alpha
c\frac{1}{\phi^{2}}v_x-\alpha u_{xx}=-v_t,\label{nreale}
\end{equation}
\begin{equation}
\lambda v-\gamma\phi^2v-\beta v_{yy}+\alpha c^2\frac{1}{\phi^{4}}v+2\alpha
c\frac{1}{\phi^{3}}\phi_x u-2\alpha c\frac{1}{\phi^{2}}u_x-\alpha v_{xx}=u_t.\label{nimage}
\end{equation}
\label{nevaluepde}
\end{subequations}
Since \eqref{nevaluepde} does not 
depend on $y$ or $t$ explicitly, we may assume that $u(x,y,t)$ and 
$v(x,y,t)$ have the form 
\begin{subequations}
\begin{equation}
u(x,y,t)=U(x,\rho,\Omega)\mbox{e}^{i\rho y+\Omega t}+c.c.,
\end{equation}
\begin{equation}
v(x,y,t)=V(x,\rho,\Omega)\mbox{e}^{i\rho y+\Omega t}+c.c.,
\end{equation}
\label{uvnontrivial}
\end{subequations}
where $\rho$ is a real constant, $U(x)$ and $V(x)$ are complex-valued
functions,  $\Omega$ is a complex constant and $c.c.$ denotes complex
conjugate.  Notice that $\rho$ is the wavenumber of the transverse 
perturbation and $\Omega$ is the exponential growth rate associated 
with $\rho$. 
If bounded $U,V$ exist such that $\Omega$ has a positive real part, 
then the amplitudes of the perturbations grow
exponentially in time and the unperturbed solution is said to be
unstable.

Upon substitution, (\ref{nevaluepde}) gives
\begin{subequations}
\label{E:neval}
\begin{equation}
\lambda U-3\gamma\phi^2U+\beta\rho^2
U+\alpha c^2\frac{1}{\phi^{4}}U-2\alpha c\frac{1}{\phi^{3}}\phi_xV+2\alpha
c\frac{1}{\phi^{2}}V_x-\alpha U_{xx}=-\Omega V,
\end{equation}
\begin{equation}
\lambda V-\gamma\phi^2V+\beta\rho^2 V+\alpha c^2\frac{1}{\phi^{4}}V+2\alpha
c\frac{1}{\phi^{3}}\phi_x U-2\alpha c\frac{1}{\phi^{2}}U_x-\alpha
V_{xx}=\Omega U.
\end{equation}
\end{subequations}
If $c=0$, then \eqref{E:neval} reduces 
to the stability analysis of TP solutions. 
This case is
examined in~\cite{JCHStp,JCBDtp,infeldz,alesh,kart,MYS} and others. 
Using the linear system \eqref{E:neval}, we are now able to 
investigate the stability of the perturbed NTP solution numerically.  

\section{Numerical investigation of spectral stability}
\label{S:Numerics}

The main difficulty for the numerical investigation of \eqref{E:neval}
is the size of the parameter space involved.  For every choice of the 
equations parameters $\alpha,
\beta$ and solution parameter pairs $(k,B)$, the spectrum of \eqref{E:neval}
needs to computed for a range of $\rho$ values to determine stability or to
analyze any instabilities.  An efficient numerical 
method is necessary.  Hill's method~\cite{BDNKhill}, which exhibits 
exponential convergence, allows for the systematic exploration of the 
large phase space encountered here.

\subsection{Hill's method}
\label{SS:method}

To apply Hill's method, Fourier expansions are needed for all coefficients 
of \eqref{E:neval}.  Using the complex Fourier form, we have
\begin{equation}
\label{E:phihat}
\begin{array}{rr}
\displaystyle{\phi^2(x) = \sum_{n = -\infty}^{\infty} Q_n e^{i 2 n \pi x / L},} &  \qquad
\displaystyle{\phi^{-2}(x) = \sum_{n = -\infty}^{\infty} R_n e^{i 2 n \pi x / L},} \\
\displaystyle{\phi^{-4}(x) = \sum_{n = -\infty}^{\infty} S_n e^{i 2 n \pi x / L},} & \qquad 
\displaystyle{\phi^{-3}(x) \phi'(x) = \sum_{n = -\infty}^{\infty} T_n e^{i 2 n \pi x / L}}, 
\end{array}
\end{equation}
where $Q_n,R_n,S_n$ and $T_n$ are the Fourier coefficients. 
We note that $\phi^2(x)$ has period $L/2$ and that $\phi(x)$ is never zero except in the TP limit cases. 

The periodicity of the coefficients in \eqref{E:neval} allows us 
to decompose the eigenfunction components $U$ and $V$ of the spectral problem in a 
Fourier-Floquet form
\beq
\label{E:UVhat}
U(x) = e^{i \mu x} \sum_{n=-\infty}^{\infty} U_n e^{-i n \pi x/ L}   
\mbox{\quad and \quad}
V(x) = e^{i \mu x} \sum_{n=-\infty}^{\infty} V_n e^{-i n \pi x/ L}.   
\eeq
The form of $U$ and $V$ in \eqref{E:UVhat} follows 
from Floquet's theorem and the observation that 
we seek eigenfunctions, which are by definition bounded.  
This decomposition has the benefit of 
admitting both periodic and anti-periodic eigenfunctions when $\mu=0$. 
Allowing the Floquet exponent $\mu$ to be different from 0 
gives rise to solutions that are either quasiperiodic or have period 
greater than $L$.  
Again, see~\cite{BDNKhill} for a more complete explanation and discussion.

Substitution of \eqref{E:phihat} and \eqref{E:UVhat} in \eqref{E:neval} 
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 

\begin{subequations}
\label{E:nevald}
\begin{multline}
\label{E:ndL1}
-\left(\lambda + \beta \rho^2 - \alpha\left(i\mu+ \frac{i n \pi}{L}\right)^2 \right) 
     U_n
+ 3 \gamma 
         \sum_{m = -\infty}^{\infty} Q_{\frac{n-m}{2}} U_m
- \alpha c^2  
         \sum_{m = -\infty}^{\infty} S_{\frac{n-m}{2}} U_m \\
 + 2 \alpha c 
         \sum_{m = -\infty}^{\infty} T_{\frac{n-m}{2}} V_m
- 2 \alpha c \left(i\mu + \frac{i n \pi}{L}\right) 
          \sum_{m = -\infty}^{\infty} R_{\frac{n-m}{2}} V_m 
    = 
        \Omega V_n 
\end{multline}
\begin{multline}
\label{E:ndL0}
\left(\lambda + \beta \rho^2 - \alpha\left(i\mu+ \frac{i n \pi}{L}\right)^2 \right) 
        V_n
- \gamma
         \sum_{m = -\infty}^{\infty} Q_{\frac{n-m}{2}} V_m
+ \alpha c^2  
         \sum_{m = -\infty}^{\infty} S_{\frac{n-m}{2}} V_m \\
+ 2 \alpha c 
         \sum_{m = -\infty}^{\infty} T_{\frac{n-m}{2}} U_m
- 2 \alpha c \left(i\mu + \frac{i n \pi}{L}\right) 
          \sum_{m = -\infty}^{\infty} R_{\frac{n-m}{2}} U_m 
    = 
        \Omega U_n. 
\end{multline}
\end{subequations}
The system of equations \eqref{E:ndL1} and \eqref{E:ndL0} 
is {\em equivalent} to the original system \eqref{E:neval}, 
and holds for all integers $n$. Here  
$\mu \in [\frac{-\pi}{K},\frac{\pi}{K})$ and
$Q_\frac{n-m}{2} = 0$ if $n-m \not\in 2\mathbb{Z}$, with 
$R_{(\cdot)},S_{(\cdot)}$ and $T_{(\cdot)}$ similarly defined.  

In practice, a 
pre-multiplication of the linear system by $\phi^4$ 
allows for an exact cosine series expansion of $\phi^2, \phi^4$ and 
$\phi^6$ to be used.  This follows from the differential equations for 
$\sn(x,k)$ and Jacobi's~\cite{Jacobi} series expansion of $\sn^2(x,k)$.
This pre-multiplication transforms the original 
eigenvalue problem into a generalized eigenvalue problem. Golub and Van Loan~\cite{GolubVanLoan} provide a brief discussion of generalized eigenvalue problems. 

\subsection{Numerical Experiments} 
\label{SS:Experiments}

By choosing a finite number of Fourier modes, thus truncating the 
exact bi-infinite system \eqref{E:nevald}, we explicitly 
construct and compute approximations to the spectral elements of 
\eqref{E:neval}.  We consider all four cases: 
{\rm (I)} 
focusing in both $x$ and $y$
($\alpha = \beta =1$), 
{\rm (II)} 
focusing in $x$ and defocusing in $y$
($\alpha = -\beta =1$),
{\rm (III)}  
defocusing in $x$ and focusing in $y$ 
($-\alpha = \beta = 1$)  and finally,
{\rm (IV)} 
defocusing in both $x$ and $y$ 
($-\alpha = -\beta =1$). 

In each case, a large number of parameter values in the two-dimensional
parameter space shown in Fig.~\ref{F:parameter_space} 
was explored numerically.  
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$ 
Fourier modes reduces the bi-infinite exact system \eqref{E:nevald} to 
an approximate $(4N+2)$-dimensional 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 \ref{T:p_vals}.
In the table, $k$ is the elliptic modulus, $B$ is the offset parameter and 
may be interpreted as a measure of the nontrivial-phase quantity $\theta$, 
$(4N+2)$ is the matrix dimension used to 
approximate the full operator, $\rho$ is the wavenumber of the perturbation in 
the $y$-dimension, and $\mu$ is the Floquet exponent. Lastly,
$\mathtt{linspace(a,b,m)}$ is a linearly spaced vector from $a$ to 
$b$ of length $m$, 
$\mathtt{logspace(a,b,m)}$ is a logarithmically spaced vector from $10^a$ to 
$10^b$ of length $m$ 
and function {\tt ceil(x)} is the smallest integer not less than $x$.

\begin{table}[t]
\bc
\begin{tabular}{|c|c|l|}
\hline
Parameter & Description & \qquad Value \\
\hline
\hline
$k$ & Elliptic modulus & \, \tt{linspace(0,1,65)} \\ 
\hline
$B$ & Shift & \begin{tabular}{l}
              For $\alpha = -1:  -\mathtt{logspace(-8,0,65)} $ \\ 
              For $\alpha = \phantom{-}1: (2k^2 + \mathtt{logspace(-8,0,65)}) \cap (2k^2 , 2)$
              \end{tabular} \\ 
 
\hline
$N$ & Fourier cutoff & \begin{tabular}{l}
                        For $\alpha=-1: 15 + \mathtt{ceil(5 k^{5})}$ \\ 
                        For $\alpha= \phantom{-}1: 10 + \mathtt{ceil(25 k^{10})}$ 
                        \end{tabular} \\
\hline
$\rho$ & Perturbation wavenumber & \, \tt{linspace(0,4,65)} \\
\hline
$\mu$ & Floquet parameter & \, $\mathtt{linspace(-\frac{\pi}{K},\frac{\pi}{K},21)}$  \\
\hline
\end{tabular}
\ec
\caption{Parameter values and ranges used in numerical experiments.{\label{T:p_vals}}}
\end{table}

\subsection{Observations}
\label{SS:obs}

First and foremost, it should be stated that {\em none} of the solutions 
considered here were found to be spectrally stable under one-dimensional and transverse perturbations.  This establishes, at 
least numerically, that {\em all} one-dimensional traveling-wave solutions of NLS of the form given by~\eqref{E:soln_form} are spectrally
unstable with respect to both longitudinal and transverse perturbations.  At this point, it remains to investigate the nature of the 
instabilities, 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 transverse
perturbations with wavenumber denoted by $\rho \in [0,4]$. Note the $\rho=0$ 
corresponds to one-dimensional, {\em i.e.} longitudinal, perturbations.
For each NLS equation ({\em i.e.} for each choice of $\alpha,\beta = \pm 1$), 
and for each parameter pair $(k,B)$, and for each perturbation of wave 
number $\rho$, a sequence of a sequence of Floquet parameters 
$\mu$ was chosen from the interval $[-\frac{\pi}{K},\frac{\pi}{K}]$. 
\begin{comment}
For each parameter triplet 
$(k,B,\rho)$, for $\alpha,\beta = \pm 1$, 
a sequence of Floquet parameters 
$\mu$ was chosen from the interval $[\frac{-\pi}{K},\frac{\pi}{K}]$. 
\end{comment}
The generalized  eigenvalues and eigenvectors were computed from the resulting 
matrix.  The eigenvalues are approximations of spectral elements 
of \eqref{E:neval}, and an approximation of the corresponding 
eigenfunctions may be reconstructed from the generalized eigenvectors.

Since a single eigenvalue with positive real part leads to 
instability of the system, the eigenvalue with largest real part over 
all choices of $\mu$ was recorded for each $(k,B,\rho)$ triplet. That is,
we compute
\begin{equation}
\Omega_{growth}(k,B,\rho) = \max_{\mu \in [-\pi/K,\pi/K]} \, Re \, \, \Omega(k,B,\rho,\mu),
\end{equation}
which we call the growth rate for a given NTP solution with solution parameters $(k,B)$ and perturbed with transverse wavenumber $\rho$. 
We reduce the dimension still further by computing the  
largest such growth rate over all sampled perturbation wavenumbers $\rho$.
This quantity, 
\begin{equation}
\Omega_{\max}(k,B) = \max_{\rho \in [0,4]}\Omega_{growth}(k,B,\rho),
\end{equation}
the maximal growth rate over all $\rho$, is plotted in the first column of 
Figs.~\ref{F:focus_growth} and~\ref{F:defocus_growth}.  This represents the maximal exponential growth 
rate that a solution with parameters $(k,B)$ can undergo in the range examined, and allows us to determine the solution which is spectrally the most unstable. We also recorded 
the minimum growth rate over all $\rho$, 
\begin{equation}
\Omega_{\min}(k,B) = \min_{\rho \in [0,4]}\Omega_{growth}(k,B,\rho),
\end{equation}
to verify that all solutions are unstable with respect to every sampled perturbation.

Every point plotted in Figs.~\ref{F:focus_growth} and \ref{F:defocus_growth} corresponds 
to a solution, and the boundaries in the figure are the boundaries of 
the regions represented in Fig.~\ref{F:parameter_space} and correspond to 
limiting TP solutions.
Fig.~\ref{F:focus_growth} corresponds to the $x$-focusing ($\alpha = 1$)
parameter range $(k,B) = (0,1) \times (2k^2,2)$ in the $\alpha = 1$ case 
of Fig.~\ref{F:parameter_space}(a). 
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]$.
Fig.~\ref{F:defocus_growth} corresponds to the $x$-defocusing ($\alpha = -1$)
parameter range of $(k,B) = (0,1) \times (-1,0)$ shown in 
Fig.~\ref{F:parameter_space}(b).  The 
transform $T_d(B) = -B$ is used in Fig.~\ref{F:defocus_growth}.
A log$_{10}$ scale is used in the vertical dimension of 
Figs.~\ref{F:focus_growt} and~\ref{F:defocus_growth}.
This causes the panels of Fig.~\ref{F:focus_growth} to become increasing sparse in their lower right corners.  
The right-hand panels of Figs. \ref{F:focus_growth} and \ref{F:defocus_growth} indicates the wavenumber $\rho$ which leads to 
maximal growth shown in the left-hand panels.  
Recall that our computations were truncated at $\rho=4$, and so in all cases 
maximums and minimums may occur at values of $\rho>4$.  

\subsubsection{Case {\rm I}: Elliptic setting with $\alpha = \beta = 1$ (Optics, BEC)}
Panels {\rm I}a and {\rm I}b of Fig.~\ref{F:focus_growth} summarize 
some properties of the 
computed instabilities in the case of focusing in both the $x$- and 
$y$-dimensions.   
\begin{comment}
Recall that the original parameter space corresponds to 
$(k,B) = (0,1) \times (2K^2,2)$,
and that the map $T_f$ is used to transform this to $(0,1) \times (0,1)$. 
This transform, coupled with the log$_{10}$ scaling, causes the plots to 
become increasingly sparse towards the bottom right corner.
\end{comment}
The lower boundary of the plot corresponds to $B = 2k^2 + (10^{-8})$, and so is 
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 $\Omega_{\max}$ in panel {\rm I}a of Fig.~\ref{F:focus_growth}.  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) as $k$ increases. It reaches a global maximum of $R_{\max} = 5.666$ around $k = 0.96$. 
Moving away from the cn boundary results in the rapid increase of $\Omega_{\max}$.  Movement away from the dn boundary results in a much slower increase in the value of $\Omega_{\max}$, as does moving away from the Stokes' wave boundary for $B$ larger than approximately 0.001.
The maximum value of $\Omega_{\max}$ over the sampled $(k,B)$ space, given by 
$R_{\max}=5.666$, is reached near $(k,B) \approx (0.99,1.98)$.  The minimum ($R_{\min} = 0.015693$) occurs near to $(k,B) \approx (0.01,0.01)$ for $\rho =4.$   

In {\rm I}b, the wavelength corresponding to the maximal growth of  
{\rm I}a is given.  In this case, the maximum instability occurs 
for the shortest wavelength samples, $\rho = 4$.
This indicates that there is a strong short-wavelength 
(large $\rho$, possible infinite) instability.

\subsubsection{Case {\rm II}: Hyperbolic setting with $\alpha = -\beta = 1$ (Deep water waves)}
Panels {\rm II}a and {\rm II}b of Fig. \ref{F:focus_growth} 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 so is 
just slightly away (in parameter space of $B$) from a cn-type solution.  The upper boundary is close to 
dn-type solution, 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, while $k=0.99$ on the left boundary. The bright soliton solution is the limiting case result when $k \to 1$.   

As in Case {\rm} I, a ridge of large instability is noticeable in 
the growth plot {\rm II}a.  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) as $k$ increases.  This ridge has a local minimum near $k=0.7$ and increases again to a global (over all admissible $(k,B)$-parameter space) maximum near $k = 0.96$.
As in the setting above, moving away from the cn-type boundary results in a rapid increase of $\Omega_{\max}$.  Moving away from the dn boundary results in a much slower increase in the value of $\Omega_{\max}$. For $B>0.001$, moving away from the boundary result in a similar slow increase in $\Omega_{\max}$.
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.~\ref{F:focus_growth}. These are centered near $(k,B) = (0.01,0.01)$ and
$(k,B) = (0.96,1.98)$, respectively.

In {\rm II}b, the wavelength corresponding to the maximal growth $R_{\max}$ of  
{\rm II}a are given.  In this case, the maximum instability of 
$R_{max}$ occurs for $\rho = 3.375.$  The surface represented by 
{\rm II}b appears to be much smoother than the surface of {\rm I}b. 

\begin{figure}[t]
\psfrag{ab 1 -1}{Ia) $\Omega_{\max}$ with $(\alpha,\beta) =  (1,-1)$}
\psfrag{ab 1 1}{IIa) $\Omega_{\max}$ with $(\alpha,\beta) =  (1,1)$}
\psfrag{flog}{$\log_{10} \,T_f(B)$}
\psfrag{k}{k}
\psfrag{Wave number 1 -1}{Ib) Associated $\rho$, $(\alpha,\beta) =  (1,-1)$  }
\psfrag{Wave number 1 1}{IIb) Associated $\rho$, $(\alpha,\beta) =  (1,1)$  }
\psfrag{rho}{$\rho$}
\begin{tabular}{cc}
 \epsfig{file=./pics/focus_plus_g.eps, width = 6cm} &
 \epsfig{file=./pics/focus_plus_r.eps, width = 6cm}  \\
 \epsfig{file=./pics/focus_minus_g.eps, width = 6cm}  &
 \epsfig{file=./pics/focus_minus_r.eps, width = 6cm}  
\end{tabular} 
\caption{Focusing in $x$-dimension plots representing a) $\Omega_{\max}$ 
and b) the wavenumber $\rho$ associated to $\Omega_{\max}$. 
R$_{\max} = \max_{k,B} \Omega_{\max}$  and 
R$_{\min} = \min_{k,B} \Omega_{\min}$. 
The vertical direction indicates the amount of 
nontrivial-phase in the underlying solution, and is plotted on a $\log_{10}$ scale. White space corresponds to ($k,B$)-parameter space that was {\em not} sampled.{\label{F:focus_growth}}}
\end{figure} 

\subsubsection{Case {\rm III}: Hyperbolic setting with $-\alpha = \beta = 1$ (Deep water waves)}
Panels {\rm III}a and {\rm III}b of Fig. \ref{F:defocus_growth} summarize 
some properties of 
the computed instabilities in the case of defocusing in the $x$-dimension 
and focusing $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 to Stokes wave solutions, while $k=0.99$ on the left boundary.   In the limit as $k \to 1$, the 
solution becomes a gray soliton.  

A distinct ridge of large instability is noticeable in the growth plot 
{\rm III}a.  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) as $k$ increases. It quickly reaches 
the global maximum (over all admissible $(k,B)$-parameter space) of $R_{\max} = 7.6375$ near $k = 0.02$ and $B=-0.0001$.
 The ridge then appears to 
decrease in amplitude as $k$ increases towards 1.
Moving away from the sn-type boundary results in a rapid increase of $\Omega_{\max}$.  Moving away from the dn-type boundary results in a much slower increase in the value of $\Omega_{\max}$. Similarly, there is much slower increase 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) \approx (0.02,0.00001)$.  The minimum exponential growth, $R_{\min}=0.015578$ is found near $(k,B) \approx (0.01,0.9)$.  Both the maximum and minimum are located near to the 
Stokes' wave boundary.   

The plot {\rm III}a indicates short wave perturbations lead to large values of 
$\Omega_{\max}$. The largest growth occurs for a 
perturbation with wavenumber of $\rho = 3.625$. 

\subsubsection{Case {\rm IV}: Elliptic setting with $-\alpha = -\beta = 1$ (Optics, BEC)}

Panels {\rm IV}a and {\rm IV}b of Fig.~\ref{F:defocus_growth} 
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$ solution, as measured in the parameter $B$. 
The left boundary of the plots, 
where $k = 0.01$, represents a region in parameter space near to the Stokes wave solutions, while $k=0.99$ on the left boundary. In the limit as $k \to 1$, the solutions become gray solitons.   


A distinct ridge of large instability is noticeable in the growth plot 
{\rm IV}a.  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) as $k$ increases, to reach 
a global maximum near $k = 0.02$ and $B=-0.0001$. The ridge then appears to 
decrease in amplitude as $k$ increases towards 1.
As in the setting above, moving away from the sn-type boundary results in a rapid increase of $\Omega_{\max}$.  Moving away from the dn-type boundary results in a much slower increase in the value of $\Omega_{\max}$. The same is true when moving away from the Stokes boundary, when $B$ is larger 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 $\Omega_{\max}$ in the panel above it. The maximum and minimum values are obtained near $(k,B) \approx (0.01,0.00009)$ and $(k,B) \approx (0.01,-1)$, respectively.  Both are located near to the Stokes' wave boundary.   

In the plot {\rm IV}b, wavenumbers corresponding to 
$\Omega_{\max}$ of {\rm IV}a are given.  It appears that a majority of 
the large values of $\Omega_{\max}$ are  
attributable to small $\rho$ (long wave) perturbations. In fact,
the largest growth occurs for $\rho = 0$, the one-dimensional perturbation. 
In contrast, short wavelength two-dimensional perturbations with wavenumber $\rho>3$ were associated with many of the large $\Omega_{\max}$ values in the previous plots.

\begin{comment}

\subsubsection{Comparisons}

Perhaps most notable in a comparison of the growth stability plots of 
Fig.~\ref{F:all_growth} is the observation that the stability 
surfaces are qualitatively very similar in all cases, and 
that the growth factor does not grow appear to grow monotonically for 
any fixed choice of $k \in (0,1)$.  Instead, the maximal growth rate 
is achieved at some point in the interior of the admissible 
$(k,B)$-parameter space.  There are other evident relationships.  
The overall qualitative structure of the growth surface is 
quite similar in both the focusing case {\rm I-II}a and in the 
defocusing case {\rm III-IV}a, suggesting that $\alpha$ may be used to distinguish these surfaces.  The position of both the 
maxima and minima appear to be influenced only slightly by the sign of the 
perturbation, both for $\alpha = 1$ and when $\alpha = -1$.  When 
$\alpha = 1$, case {\rm I-II}, the largest 
instability is achieved for relatively large values of $k$, in contrast 
the location of the maximal growth of $\alpha = -1$ in {\rm III-IV}, which 
occurs for relatively small values of $k$, and is located near the stokes 
trivial-phase limiting case.  Also, the distinct local nature of the 
maxima in the $\alpha = 1$ of cases {\rm I-II}a is different from the 
more global nature of the maxima in the $\alpha = -1$ {\rm III-IV}a 
setting.   
We also note that in all cases the minimum growth rate is found nearest 
to the stokes limit of $k=0$, and that the stability values are slightly larger in the hyperbolic 
setting {\rm II-III}a when compared to the elliptic setting {\rm I,IV}a. 

\end{comment}

\begin{figure}
\psfrag{ab -1 -1}{$IIIa) \Omega_{\max}$ with $(\alpha,\beta) =  (-1,-1)$}
\psfrag{ab -1 1}{$IVa) \Omega_{\max}$ with $(\alpha,\beta) =  (-1,1)$}
\psfrag{dlog}{$\log_{10}(-B)$}
\psfrag{k}{k}
\psfrag{Wave number -1 -1}{IIIb Associated $\rho$, $(\alpha,\beta) =  (-1,-1)$  }
\psfrag{Wave number -1 1}{IVb Associated $\rho$, $(\alpha,\beta) =  (-1,1)$  }
\psfrag{rho}{$\rho$}
\begin{tabular}{cc}
\epsfig{file=./pics/defocus_plus_g.eps, width = 6cm}  &
\epsfig{file=./pics/defocus_plus_r.eps, width = 6cm} \\
\epsfig{file=./pics/defocus_minus_g.eps, width =  6cm}  &
 \epsfig{file=./pics/defocus_minus_r.eps, width = 6cm}
\end{tabular} 
\caption{Defocusing in the $x$-dimension plots where a) is $\Omega_{\max}$ 
and b) the wavenumber $\rho$ associated with $\Omega_{\max}$. 
R$_{\max} = \max_{k,B} \Omega_{\max}$  and 
R$_{\min} = \min_{k,B} \Omega_{\min}$. 
The vertical direction indicates the amount of 
nontrivial-phase in the underlying solution, and is plotted on a $\log_{10}$ scale.  {\label{F:defocus_growth}}}
\end{figure} 

\section{Summary}

In this paper, we have considered the spectral instability of one-dimensional 
traveling-wave NTP solutions of the cubic nonlinear Schr\"{o}dinger equation.
Such solutions can be 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 part of the resulting eigenvalues was used to
determine that there are {\em no} spectrally stable NTP solutions.
The numerics indicate a well-defined ridge of maximal 
instability which is located in the $(k,B)$-parameter region associated 
with fully nontrivial-phase solutions. In addition, the numerical 
evidence indicates that exponential growth rates of dn-type solutions, 
and to a lesser extent the Stokes' wave solutions, are insensitive to 
transverse perturbation.  The growth of the cn-type and sn-type 
solutions appear to be quite sensitive to this perturbation.
In summary, numerical evidence suggests that bounded, 
nontrivial-phase one-dimensional traveling-wave solutions to the 
cubic NLS equation are unstable with respect to both one-dimensional and 
transverse perturbation.

%\bibliographystyle{plain}
\bibliographystyle{abbrv}
\bibliography{ntppaper}

\end{document}

