Nonlinear Schrödinger (NLS) equations, see  for general references, of the form
where , are omnipresent as asymptotic models for modulation of waves in applications such as hydrodynamics, nonlinear optics and Bose-Einstein condensates. Here we concentrate on the focusing cubic NLS () in dimensions. Despite the importance of the NLS equation in applications, it is notoriously difficult to treat numerically. This is due to the fact that solutions can have zones of rapid modulated oscillations, called dispersive shock waves (DSWs), and that the focusing version (1) is critical. This means that the equation can have an blow-up in finite time for initial data with an norm larger than the one of the ground state. The reader is referred to [21, 24] for an asymptotic description of the blow-up. For numerical approaches for DSWs in this context see for instance  and references therein, for blow-up in NLS see  and more recently .
These problems are only aggravated if one wants to study solutions such as the Peregrine solution  to the cubic 1D NLS,
see  for generalizations thereof, which do not vanish asymptotically, but which tend algebraically to a constant in modulus. Such solutions are often discussed as a model for rogue waves: It is claimed that the Peregrine solution has been observed experimentally in rogue wave experiments in hydrodynamics [7, 8, 9], in plasma physics  and in nonlinear optics . However, recently it has been shown numerically in [15, 4, 18] that the Peregrine solution is unstable, and its orbital instability in 1D NLS has been rigorously proven by Munõz . Interesting remarks on the Peregrine solution can be also found in . For alternative approaches via NLS equations to rogue wave phenomena see for instance  and [6, 5].
For rapidly decreasing or periodic initial data the standard numerical approach is to use Fourier methods, see  and references. Initial data with compact support allow the application of perfectly matched layers [3, 32] or transparent boundary conditions  to avoid unwanted reflections at the boundaries of the computational domains. The latter techniques, with some modification, can be successfully used in the case of algebraically decreasing data as in (2). For example, in  one of the authors presents a multi-domain spectral approach with a compactification of the infinite interval, i.e., where infinity is a regular point on the grid. The use of spectral methods for Schrödinger equations is recommended since the latter are purely dispersive equations. This means that the introduction of numerical dissipation should be strongly limited if dispersive effects like rapid oscillations are of interest.
In order to study the transverse stability of the Peregrine solution, in this paper we extend the approach of  to the 2D NLS for perturbations which are periodic or rapidly decreasing in . To this end the 1D approach in  is changed in two ways: firstly, we use a Fourier spectral method in together with the original approach in 1D from  which gives a spectral method both in and , adapted to solutions with an algebraic decay in to a finite value at infinity, and rapidly decreasing or periodic in . Secondly, we replace the implicit fourth order Runge-Kutta (IRK4) method for the time integration in  with a fourth order splitting where the IRK4 method will be only applied to the linear part of the equation. This leads to a fully explicit method the accuracy of which is studied at the example of the Peregrine solution in 1D.
The Peregrine solution can be seen as an exact -independent, and thus not localized in , solution to the 2D NLS equation (1). With the code presented in this paper, we are able to study transverse perturbations of the Peregrine solution in 2D. We consider various types of such perturbations, both localised and non-localised in . It is shown that the Peregrine solution is unstable against all studied perturbations, and that in certain cases even a blow-up can be observed.
The paper is organized as follows: in Section 2 we present the numerical approach for the time integration for the 1D code and test its convergence for the example of the Peregrine solution; the code is then extended to 2D. In Section 3 we consider perturbations of the Peregrine solution localized in . Perturbations not localized in are studied in Section 4. We add some concluding remarks in section 5.
2. Numerical approach
In this section we summarise the numerical approaches applied in this paper. Firstly, we recall the main features of the 1D code , then we present a fully explicit variant for the time integration scheme of the code and test it at the example of the Peregrine solution. The resulting code is generalised to the 2D NLS (1) where the dependence on the transverse variable is treated with Fourier methods.
2.1. Spatial discretisation
Spectral methods, see for instance 
and references therein, are an attractive choice for the numerical solution of partial differential equations (PDEs), in particular in higher dimensions, because of their excellent approximation properties for smooth functions. If the function to be approximated is analytic on the considered interval, the numerical error whilst applying spectral methods is known to decrease exponentially with the resolution, a feature which is called spectral convergence. An interesting additional advantage of spectral methods is that they minimize the introduction of numerical dissipation in the approximated PDEs. This is especially attractive in the present context of dispersive PDEs where dissipation could suppress the dispersive effects one wants to study as DSWs.
The basic idea of spectral methods is to approximate a function on a given interval by a set of global functions on this interval. The best known such method works for a periodic setting and is the truncated Fourier series, i.e., a trigonometric polynomial. It is well known that the numerical error in approximating an analytic periodic function via a truncated Fourier series decreases exponentially with the number of terms in this sum. In a finite precision approach the same is true for rapidly decreasing functions if the period is chosen large enough that the function and its first derivatives vanish with numerical precision at the domain boundaries. This is the approach we will use for the variable .
Nevertheless, Fourier techniques are also known to be much less efficient in the case of discontinuous periodic functions. In this case the numerical error in approximating the function via a truncated Fourier series is known to decrease only linearly with the number of terms in the sum, and there is a Gibbs phenomenon at the discontinuity in this case. Thus, it is considerably better to approximate functions which are only smooth on intervals by a sum of polynomials there. Slowly decreasing functions on the whole real line can be spectrally treated through a compactification of the external domain, namely
where is some constant. This means we cover the real line with two domains which are both mapped to the interval .
Note that in such a multi-domain method, in principle an arbitrary number of intervals can be used. This allows to allocate resolution were needed, whereas a compactification of the real line to a simple domain, for instance as in  or via the well known map , imposes the distribution of the numerical resolution for a given discretisation. If infinity is not a regular point as here, one can also use two infinite intervals as in , and . In the context of finite element methods, such intervals are called infinite elements. Compactification techniques are standard for elliptic equations, for instance in astrophysics, see  for a similar approach as here.
For the case of the 1D NLS, this means that we solve in interval I () the equation
and in domain II ()
Put together domain I and II cover the whole real line, and the solution has to be at least and thus in particular to be differentiable at the domain boundaries. This implies the conditions
On each of the two domains we introduce standard Chebyshev collocation points, , , and , ,
and approximate the considered function by the Lagrange interpolation polynomial sampled at the collocation points. The derivatives of the functionare approximated by the derivatives of the Lagrange polynomial which leads to the action of a Chebyshev differentiation matrix on the function, see the discussion in [25, 28]. This means that the PDEs (4) and (5
) are approximated by the finite dimensional system of ordinary differential equations (ODEs)
where in an abuse of notation we have used the same symbols for the functions
and the vector with the components, , and similarly for . Note that we always use an even number of points in domain II in order to make sure that infinity is not a grid point. This simplifies the treatment of the singularity in the equation for .
This approach can be generalised to two spatial dimensions in the following way, where we assume that is either periodic or rapidly decreasing in : in this case, after the discretisation in , the components of vector depend on (the case for is completely analogous). The -dependence is then approximated via a truncated Fourier series for ,
This means we discretise also , and use a Fast Fourier transform
(FFT) to compute the discrete Fourier transform of. Derivatives in are then approximated in standard way by multiplication of the Fourier coefficients in the sum by a factor . This yields a spectral method in both and , i.e., a combined Chebyshev-Fourier approach for the spatial coordinates. Since for a function both analytic in and , the Chebyshev and the Fourier coefficients have to decrease exponentially, we can control the spatial resolution of the solution via the behavior of the spectral coefficients. Note that the Chebyshev coefficients can be computed with a Fast Cosine Transform which is close to the FFT, see .
2.2. Explicit time integration scheme
The spatial discretisation of the previous subsection leaves us with a finite dimensional ODE system of the form
where is an dimensional matrix, is a linear operator formed essentially by a discretised version of the Laplace operator in (1), and where is a nonlinear operator free of derivatives. Since there are second derivatives in the linear part of NLS equations, the resulting system is stiff which means that explicit time integration schemes will be inefficient since stability conditions would impose severe restrictions on the size of the time steps. This is especially true since the operator has a conditioning of order , see for instance the discussion in , and similarly for .
We therefore chose in  an implicit time integration scheme to address these stability issues. Concretely, we apply a fourth order Runge-Kutta (IRK4) scheme, the Hammer-Hollingsworth method, a 2-stage Gauss scheme. The general formulation of an -stage Runge–Kutta method for the initial value problem is as follows:
where are real numbers and . For the IRK4 method used here, one has , , , , and .
The matching conditions (6) are implemented via a -method, see . The resulting implicit system for and requires an iterative solution. For this we used a simplified Newton method to invert the operators of the left hand side of (11) and get the new and ). The iteration converges rapidly in simple cases, but can take more than 100 iterations to reach high precision in some of the examples discussed in the following.
In the linear case , the solution to (11) can be written in the form
An inversion of the operator containing is not recommended since the conditioning of the matrix would be of the order . Furthermore, in order to get a unique solution, one would have to impose a condition at the domain boundaries which is in itself problematic. Therefore we used in  even in the linear case the iterative approach based on (11). In this paper however, we are solving (12) by inverting first and then in (12), each of them with a condition. Thus we only have to deal with the conditioning of the operator , but in contrast to the iterative approach we have to do that twice. In practice this implies, as will be discussed in more detail below for the Peregrine example, that the maximally achievable accuracy is of the order of instead of in the iterative solution. Nevertheless this is more than sufficient for our purposes, and it is significantly more efficient as no iterations are needed.
Using this approach for the nonlinear system (11) can be done through splitting, i.e., to split the equation into the linear part (12) and the nonlinear part which can be integrated explicitly. The linear part will then be integrated with IRK4 in the form (12). The motivation for splitting methods comes from the Trotter-Kato formula [27, 14]
where and are certain unbounded linear operators, for details see .
The idea of these methods for an equation of the form is to write the solution as
where and are sets of real numbers that represent fractional time steps. Yoshida  gave an approach which produces split step methods of any even order. We apply here a fourth order splitting. Combined with IRK4 for the linear step, this gives a fourth order splitting scheme for NLS with an inexact solution of the linear step via an implicit fourth order method. However, the latter can be done explicitly because of the linearity of the equation. The nonlinear step can be integrated as usual for NLS in explicit form since for the equation , is constant in time. Thus we get in total a fully explicit scheme which, especially in two dimensions, is much more efficient than the IRK4 scheme for the full equation with an iterative solution.
The numerical accuracy of the time integration can be controlled via the conserved quantities of the NLS equation of which there are 3 in two dimensions (the 1D cubic NLS is of course completely integrable and thus has an infinite number of conserved quantities). Since we want to treat solutions as the Peregrine solution which do not tend to zero at infinity, but satisfy , we consider a combination of the energy and the norm which is defined in such cases and which is a conserved quantity of NLS,
Due to unavoidable numerical errors, this quantity will not be exactly conserved in actual computations. However the relative (with respect to the initial value) conserved gives an indication of the numerical error (as discussed in , it typically overestimates the numerical accuracy by one to two orders of magnitude). Note that (14) is not defined for the Peregrine solution on . But it is for , the situation we study in this paper.
As a test of the above code we use the Peregrine solution (2) which can be seen as a independent solution of the 2D NLS equation. Since this solution does not test the dependence in the code, the results are the same as with the 1D code, and we concentrate on the latter to study the dependence of the numerical error on the time step. Concretely we impose the initial data and solve the NLS equation for these initial data. The numerical and the exact solution are compared for , as the numerical error we take the norm of the difference between both solutions. We always use and collocation points in the respective domains, and vary the number of time steps between and .
The Chebyshev coefficients when representing the Peregrine solution in the form where , are the Chebyshev polynomials can be seen on the left of Fig. 1, the coefficients at the final time on the right of the same figure. This shows that the solution is always well resolved in space since the Chebyshev coefficients decrease in both domains to the order of the rounding error which is of the order of .
In Fig. 2 we show the norm of the difference between numerical and exact solution for in dependence of the number of time steps . It can be seen that this is a fourth order method by comparing it to a straight line with slope in the same figure. Below an accuracy of the order of the convergence becomes slower than this since the conditioning of the Chebyshev differentiation matrices becomes important. This shows that numerical errors below can be easily reached, but that for higher precisions the code from  based on an iterative approach has to be used. On the right of Fig. 2 we show the numerically computed quantity (14) in dependence of . It shows a very similar behavior as the error and is even slightly larger. Note that we cannot consider a relative error here since (14) vanishes for the Peregrine solution.
3. Localized perturbations
In this section we consider localized perturbations of the Peregrine solution. Concretely we study the time evolution of the 2D NLS equation for the initial data
with being a constant. We consider the 4 cases and . All examples are computed with and polynomials for the -dependence, Fourier modes for and time steps.
We first consider the case and in Fig. 3, on the left the initial condition, on the right the solution for . The solution is clearly unstable in the sense that the initial perturbations grow.
To show that the solution does not stay close to the Peregrine solution, we show in Fig. 4 the solutions for and for (minimal modulus of for and ) and (one of the maxima of for and ) together with the Peregrine solution.
The situation is similar in the case and as can be seen in Fig. 5: on the left the modulus of the solution for is shown, on the right the solution for together with the Peregrine solution (dotted). Again the perturbed solution does not stay close to the Peregrine solution.
The situation changes somewhat in the case and . In this case the solution appears to blow up for (the relative conservation drops in this case below and the code is stopped). The solution for can be seen on the left of Fig. 8. The norm of the solution in dependence of time on the right of the same figure also indicates a blow-up. This would imply that the Peregrine solution is strongly unstable as a solution to the 2D focusing NLS equation.
4. Nonlocalized perturbations
In this section we study perturbations nonlocalized in , but localized in of the Peregrine solution of the form , where . The same numerical parameters as in the previous section are applied.
In Fig. 10 we show the solution on the right of Fig. 9 for and (the maximum of the modulus of ) together with the dotted Peregrine solution. It can be seen that the perturbed NLS solution does not stay close to the latter.
The solution for the case can be seen at different times in Fig. 11. The initial perturbation obviously grows, and the formed structures appear even to blow up in this case for .
A blow-up is suggested also by the norm of the solution which seemingly diverges as can be seen on the left of Fig. 12. The spectral coefficients on the right of the same figure also indicate a loss of resolution because of said blow-up.
A blow-up once more suggests that the Peregrine solution is strongly unstable as a solution to the 2D focusing NLS equation. Whereas this is numerically difficult to decide, the solution is clearly unstable in the sense that it does not stay close to the exact solution at the same time as can be seen in Fig. 13, where the solution for two values of corresponding to the shown minimum and the maximum of in blue and the corresponding Peregrine solution in red is given.
In this paper we have presented a numerical approach to the 2D focusing cubic NLS equation based on a multi-domain spectral approach with a compactified exterior domain in the -coordinate (which allows to treat infinity as a finite grid point) and a Fourier spectral method in the coordinate . This provides a spectral approach in both spatial coordinates which allows to achieve a numerical error exponentially decreasing with the resolution for functions periodic or rapidly decreasing in and rational and bounded at infinity in . For the time integration a fully explicit fourth order method was presented based on a 4th order splitting scheme and an IRK4 method for the linear step.
This code allows to treat perturbations of the Peregrine solution with high precision and efficiency. In the examples studied the known instability of the Peregrine solution in 1D is shown to be amplified in 2D. All considered perturbations grow in time. In addition, numerical evidence for strong instability of the Peregrine solution in this case are presented: the solution appears to blow up in finite time if the initial data have a (14) (which is finite on the considered torus in ) is smaller than the one of the Peregrine solution.
for which an efficient computation of matrix exponentials will be needed. In particular it will be necessary to filter the unphysical eigenvalues of the Chebyshev differentiation matrices, see the discussion in.
-  M.J. Ablowitz and T.P. Horikis, Interacting nonlinear wave envelopes and rogue wave formation in deep water, Physics of Fluids 27, 012107 (2015); https://doi.org/10.1063/1.4906770
-  H. Bailung, S. K. Sharma, and Y. Nakamura, Observation of Peregrine solitons in a multicomponent plasma with negative ions, Phys. Rev. Lett. 107, 255005 (2011).
-  J. Bérenger. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 114, 185-200 (1994).
-  M. Birem and C. Klein, Multidomain spectral method for Schrödinger equations, Adv. Comp. Math. DOI 10.1007/s10444-015-9429-9 (2015)
-  J. Bona and J.-C. Saut, Dispersive Blow-Up II. Schrödinger-Type Equations, Optical and Oceanic Rogue Waves, Chinese Annals of Math. Series B, 31, (6), 793-810 (2010).
-  J.L. Bona, G. Ponce, J.-C. Saut, C. Sparber, Dispersive blow-up for nonlinear Schr dinger equations revisited, JMPA, 102, 782–811 (2014).
-  A. Chabchoub, N. P. Hoffmann, and N. Akhmediev, Rogue wave observation in a water wave tank, Phys. Rev. Lett., 106, 204502 (2011).
-  A. Chabchoub, N. Hoffmann, M. Onorato, and N. Akhmediev, Super rogue waves: observation of a higher-order breather in water waves, Phys. Rev. X 2, 011015 (2012).
-  A. Chabchoub, N. Hoffmann, H. Branger, C. Kharif, and N. Akhmediev, Experiments on wind-perturbed rogue wave hydrodynamics using the Peregrine breather model, Physics of Fluids, 25, DOI: 10.1063/1.4824706 (2013).
-  P. Dubard, P. Gaillard, C. Klein, and V.B. Matveev, On multi-rogue wave solutions of the NLS equation and positon solutions of the KdV equation, Eur. Phys. J. Special Topics, 185, 247–258 (2010).
-  Grosch, C.E., Orszag, S.A.: Numerical solution of problems in unbounded regions: coordinate transforms. J. Comput. Phys. 25, 273 296 (1977)
-  M. Hochbruck and A. Ostermann, Exponential integrators, Acta Numerica (2010), pp. 209 286.
-  A.-K. Kassam and L. Trefethen, Fourth-Order Time-Stepping for stiff PDEs, SIAM J. Sci. Comput., 26 (2005), pp. 1214 1233.
-  T. Kato, Trotter’s product formula for an arbitrary pair of self-adjoint contraction semigroups, vol. 3, Academic Press, Boston, 1978, pp. 185 195.
-  U. Al Khawaja, H. Bahlouli, M. Asad-uz-zaman, S.M. Al-Marzoug, Modulational instability analysis of the Peregrine soliton, Commun Nonlinear Sci Numer Simulat 19, 2706 2714 (2014).
-  B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, G. Genty, N. Akhmediev, and J. M. Dudley, The Peregrine soliton in nonlinear fibre optics, Nat. Phys. 6, 790–795 (2010).
-  C. Klein and N. Stoilov, A numerical study of blow-up mechanisms for Davey-Stewartson II systems, Stud. Appl. Math., DOI : 10.1111/sapm.12214 (2018)
-  C. Klein and M. Haragus, Numerical study of the stability of the Peregrine breather, Annals of Mathematical Sciences and Applications, 2(2), 217-239 (2017).
-  C. Klein, Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equation, ETNA, 29, 116–135 (2008).
-  Y.C. Li, On the so-called rogue waves in the nonlinear Schrödinger equation, preprint, arXiv:1511.00620v1 (2015)
-  F. Merle and P. Raphael. On universality of blow-up profile for l2 critical nonlinear Schrödinger equation. Inventiones Mathematicae, 156:565 672, 2004.
-  C. Muñoz, Instability in nonlinear Schrödinger breathers, Proyecciones Vol. 36, no. 4 (2017) 653 683.
-  D.H. Peregrine, Water waves, nonlinear Schrödinger equations and their solutions, J. Austral. Math. Soc. B, 25, 16–43, DOI:10.1017/S0334270000003891 (1983).
-  C. Sulem and P.L. Sulem. The nonlinear Schrödinger equation. Springer, 1999.
-  L. N. Trefethen, Spectral Methods in Matlab, SIAM, Philadelphia, PA, 2000.
-  Trefethen, L.N., 2013. Approximation theory and approximation practice (Vol. 128). Siam.
-  H. Trotter, On the Product of Semi-Groups of Operators, Proceedings of the American Mathematical Society, 10 (1959), pp. 545 551.
-  Weideman, J.A.C. and Reddy, S.C., A Matlab differentiation matrix suite, ACM TOMS, 26 (2000), 465–519.
-  www.lorene.obspm.fr
-  H. Yoshida, Construction of higher Order symplectic Integrators, Physics Letters A, 150 (1990), pp. 262–268.
-  V.E. Zakharov and A.B. Shabat, Exact theory of two-dimensional self-focusing and one-dimensional selfmodulation of waves in nonlinear media. Sov. Phys. JETP, 34(1), 62–69 (1972); translated from Zh. Eksp. Teor. Fiz. 1, 118–134 (1971).
-  C. Zheng. A perfectly matched layer approach to the nonlinear Schrödinger wave equations. J. Comput. Phys. 227, 537-556 (2007).
-  C. Zheng, Exact nonreflecting boundary conditions for one-dimensional cubic nonlinear Schrödinger equations, J. Comput. Phys. 215, 552-565 (2006).