Dynamics of Plane Waves in the Fractional Nonlinear Schrödinger Equation with Long-Range Dispersion

: We analytically and numerically investigate the stability and dynamics of the plane wave solutions of the fractional nonlinear Schrödinger (NLS) equation, where the long-range dispersion is described by the fractional Laplacian ( − ∆ ) α /2 . The linear stability analysis shows that plane wave solutions in the defocusing NLS are always stable if the power α ∈ [ 1,2 ] but unstable for α ∈ ( 0,1 ) . In the focusing case, they can be linearly unstable for any α ∈ ( 0,2 ] . We then apply the split-step Fourier spectral (SSFS) method to simulate the nonlinear stage of the plane waves dynamics. In agreement with earlier studies of solitary wave solutions of the fractional focusing NLS, we ﬁnd that as α ∈ ( 1,2 ] decreases, the solution evolves towards an increasingly localized pulse existing on the background of a “sea” of small-amplitude dispersive waves. Such a highly localized pulse has a broad spectrum, most of whose modes are excited in the nonlinear stage of the pulse evolution and are not predicted by the linear stability analysis. For α ≤ 1, we always ﬁnd the solution to undergo collapse. We also show, for the ﬁrst time to our knowledge, that for initial conditions with nonzero group velocities (traveling plane waves), an onset of collapse is delayed compared to that for a standing plane wave initial condition. For defocusing fractional NLS, even though we ﬁnd traveling plane waves to be linearly unstable for α < 1, we have never observed collapse. As a by-product of our numerical studies, we derive a stability condition on the time step of the SSFS to guarantee that this method is free from numerical instabilities. they are unstable for α ∈ ( 0, 1 ) , and the instability is substantially affected by the wave number λ k . It shows that unstable modes appear in the low- µ region, and their number generally increases with | λ k | . If


Introduction
The fractional Schrödinger equation was first introduced in [1,2] as the path integral of Lévy trajectories. Later, it was also derived from a microscopic lattice system with long-range interactions when passing to the continuum limit [3]. In contrast to its classical counterpart, the fractional Schrödinger equation can model nonlocal interactions, which enables it to describe complex phenomena in various fields; see in [1,2,[4][5][6][7][8][9] and references therein. In this paper, we refer to these nonlocal interactions as "long-range dispersion", because in the classical limit (i.e., for α = 2 in (1) below), the corresponding term describes group velocity dispersion of linear waves. However, the nonlocality of the fractional Schrödinger equation also introduces considerable challenges in understanding the properties of its solutions [10][11][12][13]. In this paper, we analytically and numerically study the dynamics of the simplest-plane wave-solutions of the fractional nonlinear Schrödinger (NLS) equation and reveal effects of the long-range dispersion on the stability and dynamics of plane wave solutions.
We consider the one-dimensional fractional NLS of the following form [2,3,5,6,11,14]: where i = √ −1 is the imaginary unit, and u(x, t) is a complex-valued wave function. The constant γ ∈ R describes the strength of local (or short-range) nonlinear interactions, where the interactions are repulsive (defocusing) if γ > 0 or attractive (focusing) if γ < 0. The fractional Laplacian (−∆) α 2 is defined via a pseudo-differential operator [5,[15][16][17]: where F and F −1 denote, respectively, the Fourier transform and its inverse, and µ is the wave number. In a special case of α = 2, the definition in (2) reduces to the classical Laplace operator −∆, and thus Equation (1) becomes the classical (non-fractional) cubic NLS. In this paper, we will focus on the case with 0 < α ≤ 2, in which case (2) represents the action of a nonlocal operator that describes the long-range dispersion mentioned above [3,10,18,19]. Note that physical relevance of the case α < 1 has not been universally accepted. Nonetheless, this case was considered for the Majda-McLaughlin-Tabak model [7], and therefore we also include it into the scope of our analysis. The fractional NLS (1) has two important conserved quantities: the L 2 -norm, or mass of the wave function, and the total energy, where ∇ α = −(−∆) α 2 . These conserved quantities can be used as benchmarks to test the performance of numerical methods for the fractional NLS.
Plane wave solutions are the simplest nontrivial solutions of the NLS. Understanding their dynamics and stability provides insight into the behavior of more complicated solutions. In the classical (α = 2) NLS, stability of the plane wave solutions is well known to critically depend on the sign of nonlinearity: If the nonlinearity is defocusing (γ > 0), the plane waves are always stable, whereas they are subject to long-wave instability in the focusing (γ < 0) case. This instability is known as the Benjamin-Feir instability in the context of water waves or as the modulational instability in the context of plasma physics and nonlinear optics. Unstable perturbations in the classical NLS do not grow without bound. For the fractional NLS, the dynamics of even the simplest nontrivial solution-the plane wave-is still not completely understood. One contribution of this paper is to study stability of both standing and traveling plane waves with respect to small perturbations. To our knowledge, modulational instability of the plane wave solutions of the fractional NLS was studied in [20]; however, its authors considered a different-spatio-temporal, in three dimensions-model, and therefore its results cannot be compared to those presented here.
While linear stability of plane wave solutions can be studied analytically, the nonlinear stage of their dynamics has to be studied numerically. One of the most popular numerical methods in studies of the NLS has been the split-step Fourier spectral (SSFS) method [11,[21][22][23][24][25][26][27]. The SSFS method preserves exactly the mass conservation and dispersion relation of the NLS (1). However, the SSFS method could also introduce numerical instabilities which have no analytical counterparts. To interpret its results with confidence, one thus needs to ensure that the method does not introduce numerical, or spurious, instability. Such numerical instability of the SSFS method in simulating the plane wave solutions of the classical NLS has been studied in, e.g., [21,[23][24][25][26], but it is missing in the literature on the fractional NLS. A second contribution of this paper is a detailed study of the conditions that ensure numerical stability of the SSFS method when it is used to simulate plane wave solutions of the fractional NLS.
Finally, we study with the SSFS the nonlinear stage of the evolution of perturbed plane waves. Our results for standing plane waves (i.e., those not moving in the reference frame of Equation (1)) are consistent with earlier results by the authors of [14,28], which showed increased localization of the field as α decreased and, in particular, demonstrated wave collapse in finite time for α = 1. However, to our knowledge, previous studies did not consider dynamics of traveling solutions. (The only exception known to us is in [29]; however, its authors proved the existence of a localized traveling wave solution for 1 < α < 2, but did not actually find its analytical or numerical form.) Thus, a third main contribution of this paper is a study of the nonlinear dynamics of traveling (plane) waves, which turns out to be quantitatively (although not qualitatively) different from that of standing waves. This paper is organized as follows. In Section 2, the stability of plane wave solutions in the fractional NLS is studied by a linear stability analysis. In Section 3, the SSFS method is described, and its numerical stability condition is obtained. In Section 4, we numerically study the dynamics of the unstable plane wave solutions in the fractional NLS, revealing the role played by the long-range dispersion. Conclusions of this study are made in Section 5.

Plane Wave Solution and Its Stability
The fractional NLS (1) admits the plane wave solution of the form provided that the following nonlinear dispersion relation is satisfied The constants a and λ k are the amplitude and the wave number, respectively, of the plane wave. A wave packet consisting of plane waves (5) with wave numbers in a narrow interval near λ k = 0 would move with a group velocity where sgn(·) denotes the sign function. For this reason, we will refer to a plane wave with λ k = 0 as a traveling plane wave, contrasting it with a standing plane wave, which has λ k = 0. In simulations reported in subsequent sections, we will consider a finite domain of size L, i.e., [−L/2, L/2], and thus the wave number is chosen as Without loss of generality, we will assume λ k ≥ 0 in this paper. Due to the nonlinearity of NLS, the plane wave solution (5) is stable only under certain conditions. Here, we will focus on its linear stability analysis in the NLS (1) with the power α ∈ (0, 2]. Consider a perturbed solution of the form where u(x, t) is the plane wave solution in (5), and ε(x, t) is a small perturbation function satisfying |ε| 1. We assume that ε is periodic, and it can be expressed in Fourier series as with µ l = 2lπ/L, for l ∈ Z. In what follows, we set Re( ε 0 (0)) = 0 without loss of generality. Indeed, having Re( ε 0 (0)) = 0 simply rescales the amplitude a of the unperturbed wave (5) in the order O(ε).
Substituting (8) into the NLS (1) and taking the leading order terms of ε, we obtain where ε * represents the complex conjugate of ε. Taking the Fourier transform of (9) and its complex conjugate, we get the following system of ODEs: where The eigenvalues of matrix G l are where The solution of system (10) grows exponentially if at least one of the eigenvalues Λ ± l has positive real part, which occurs if In other words, if the condition (13) holds for at least one l ∈ Z\{0}, the plane wave solution (5) is unstable.
Below, we will focus on how the instability condition (13) depends on the power α and the wave number λ k . It is well known that the plane wave solutions in the classical (α = 2) NLS are always stable in the defocusing cases and unstable in the focusing cases. Moreover, stability of the plane wave solutions in the classical NLS is independent of its wave number λ k . In contrast, both these well-established facts no longer hold in the fractional (0 < α < 2) NLS: the plane wave may be unstable in the defocusing case, and its (in)stability is affected by λ k . These situations are delineated in Theorems 1 and 2 below. Theorem 1 ( Defocusing NLS with γ > 0). (i) For α ∈ (0, 2], a standing plane wave is always stable to infinitesimal perturbations. (ii) For α ∈ [1, 2], a traveling plane wave is always stable, independent of the wave number λ k . (iii) For α ∈ (0, 1), a traveling plane wave could be unstable if γ|a| 2 is sufficiently large. Moreover, the growth rate of this instability is maximized at modes with |µ l | = |λ k |, if the following condition is satisfied: When α = 2, our conclusions in Theorem 1 (i)-(ii) reduce to those in the literature [23,25] that a plane wave solution is always stable in the classical defocusing NLS, independent of its wave number. In the fractional NLS, the standing plane wave is always stable, but the stability of traveling plane waves crucially depends on the fractional order α. Figure 1 shows that for α ∈ [1,2], traveling plane waves are always stable, confirming the conclusion in Theorem 1 (ii). By contrast, they are unstable for α ∈ (0, 1), and the instability is substantially affected by the wave number λ k . It shows that unstable modes appear in the low-µ region, and their number generally increases with |λ k |. If the condition (14) is satisfied, the growth rate of the instability reaches its maximum at µ l = ±λ k (see Figure 1a,b), confirming the conclusion in Theorem 1 (iii). However, if the condition (14) is not satisfied, there exist stable bands centered at ±λ k , which are surrounded by unstable modes (these stable bands are clearly visible in Figure 1d for α ∈ (0, 0.735)). Theorem 2 (Focusing NLS with γ < 0). (i) For α ∈ (0, 2], a standing plane wave is unstable, provided that there is at least one l ∈ Z\{0} such that |µ l | α < 2|γ||a| 2 . Moreover, the number of unstable pairs is given by where · denotes the floor function. (ii) For α ∈ (0, 2], a traveling plane wave is unstable, if L and/or |γ||a| 2 are sufficiently large. Theorem 2 (i) implies that if |µ 1 | α > 2|γ||a| 2 , all standing plane waves in the focusing NLS are stable; see the cases of α ∈ (0, 0.15) in Figure 2a.
Otherwise, unstable modes appear in the low-µ region, i.e., the region of |µ l | < (2|γ||a| 2 ) 1/α . It implies that for α ∈ (0, 1), the bandwidth of instability is strongly affected by the nonlinearity |γ||a| 2 . In the special case of α 1, the bandwidth of instability is exponentially small if 2|γ||a| 2 < 1 (see Figure 2a), while if 2|γ||a| 2 > 1, this bandwidth is exponentially large (see Figure 2b, the thin light-green "band" near bottom). Figure 3a-d presents the stability situation of traveling plane waves in the focusing NLS, which is more complicated than that of standing plane waves. For α ∈ (1, 2], unstable modes always lie in a low-pass band, i.e., instability occurs for small |µ l |. By contrast, two distinctive phenomena are observed for α ∈ (0, 1]: (i) there is a "stable gap" in the low-µ region. Moreover, the bandwidth of this "stable gap" is around 1 + (2 − 2 α ) 1/α |λ k | for α ∈ (0, 1], and in the special case of α = 1 it reduces to λ k . Indeed, there is ∆ l < 0 in (12) for the wave numbers satisfying |µ l | < µ cr , where the critical wave number µ cr is defined as the positive root of |λ k + µ| α + |λ k − µ| α − 2|λ k | α . Therefore, the wave numbers satisfying |µ| ≤ µ cr are stable. (ii) The bandwidth of instability expands rapidly as α decreases. To show it, we take the Taylor expansion of (12) for |µ l | λ k and obtain Note that despite condition |µ l | λ k , the term of λ k cannot be neglected if α is sufficiently small. (c) (d) Figure 3. Contour plots of Re(Λ + l ) to show the stability of traveling plane wave solutions with a = 0.5 in the focusing NLS. For α > 1, Re(Λ + l ) > 0 (i.e., instability is present) for |µ l | 1, although their values are small. In contrast, there is a "stable gap" at low-µ region if α < 1 In fact, it is this term that is responsible for the expansion of the unstable band compared to the case of λ k = 0. Indeed, the instability condition (13) along with (16) implies the unstable modes: As the exponent 1/α > 1, it suffices for the base to exceed one in order to make |µ l | > 1, which makes the right-hand side rapidly increase as α is decreased.

Split-Step Fourier Spectral Method
In this section, we first briefly describe the split-step Fourier spectral (SSFS) method for solving the fractional NLS (1) and discuss some of its well-known properties. For our purposes, the first-order SSFS method suffices; if higher accuracy is sought in time, higher-order SSFS method are well known [11,22,27,30,31]. In Section 3.2, we investigate the numerical stability of the SSFS method when it is used to compute plane wave solutions of the NLS.

Description of the SSFS Method
Let τ > 0 be the time step, and define t n = nτ for n = 0, 1, . . . From time t = t n to t = t n+1 , one numerically solves the NLS (1) in two steps: For t ∈ (t n , t n+1 ], Equation (17) conserves the quantity |u(x, t)| ≡ |u(x, t n )|, and consequently it reduces to a linear equation for u(x, t). Integrating it in time gives the solution to (17): To obtain a numerical solution of (18), we consider a computational domain [−L/2, L/2] with periodic boundary conditions. Due to the definition of the fractional Laplacian in (2), the Fourier spectral method is the most natural method for spatial discretization of (18). Choose an even integer J > 0, and define the mesh size h = L/J and spatial grid points where the m-th Fourier mode u m (t) is defined as Substituting (19) into (18), taking the discrete Fourier transform, and integrating the resulting ordinary differential equation for u m (t), one obtains Combining (20) with (19) leads to an approximate solution of (18). Finally, denoting the numerical approximation of u(x j , t n ) by u n j , one obtains the first-order SSFS method as The SSFS method (21)- (22) has spectral accuracy in space and first-order accuracy in time; the latter can be easily improved by employing higher-order split-step schemes [11,22,27,30,31]. This method is explicit and thus easy to implement. In addition, it exactly preserves the discrete L 2 -norm (mass) of the solution, as well as the dispersion relation ω = |λ k | α + γ|a| 2 of plane wave solutions of the fractional NLS for any α ∈ (0, 2] [22]. These quantities are also preserved by higher-order SSFS method. Due to all these features, the SSFS is one of the most popular methods for solving the classical NLS; for the same reasons it can be equally effectively applied to the fractional NLS.

Stability of the SSFS Method
In this subsection, we study numerical stability of the SSFS method (21)-(22) when it is used to simulate plane wave solutions in (5). The stability condition that we will obtain is essential to enable one to simulate the fractional NLS time-efficiently. That is, it allows one to select the time step small enough so as to avoid numerical artifacts, and yet not too small so as to avoid excessively large simulation times. Our approach is substantially different from those used in [24,25], due to the facts that α = 2 and λ k = 0. Here, we assume where µ max is the maximum wave number resolved by the spectral grids. Note that this assumption does not restrict the value of λ k to be small, as one can have, e.g., λ k = 0.1µ max 1 for a sufficiently large µ max . At time t = t n , we consider a perturbed plane wave solution: where ε n j , with |ε n j | 1, is a small periodic perturbation on [−L/2, L/2] and can be expanded as: In contrast to (19), a different range of indices is used in the summation of (25). This is the main difference in the analysis for nonzero wave number λ k = 0 from the by-now standard analysis for λ k = 0 [25]. In fact, the perturbation term in (24) is u n j ε n j , which has effective wave number of (λ k + µ l ) ≡ µ k+l rather than just µ l . To ensure this wave number fits into the range of [µ min , µ max ], we have to require µ l ∈ [µ min − λ k , µ max − λ k ], resulting in asymmetric summation indices in (25). In order to emphasize this difference compared to the symmetric summation indices in Section 3.1, we use a different letter-l in (25) instead of m in (19). This labeling convention of indices will be used in the remainder of this subsection: m and l, respectively, for the Fourier transform of the solution, u m , and the multiplicative perturbation, ε l .
Setting n = n + 1 in (24) and taking its discrete Fourier transform, one obtains On the other hand, u n+1 m can be independently found from the SSFS scheme (21)- (22). To this end, we substitute u n j from (24) into (21) and take the leading order terms of ε n j to obtain Substituting (27) into (22) and taking its discrete Fourier transform, we obtain For m = k, we have used the convention Re( ε n 0 ) = 0, adopted in Section 2. Furthermore, we have split the wave number indices l ≡ (m − k) of the perturbation terms ε n m−k into the "main" and "aliased" sets defined as To understand this choice, note that when (m − k) ≡ l is within one of these sets, then the indices of both ε-and ε * -terms in each line in (28) are within the same set. For S main , this is evident, while for S aliased this follows from the fact that discrete Fourier transform "wraps" (i.e., aliases) indices modulo J when they fall outside the range [−J/2, J/2 − 1]. Note that the aliased region S aliased was considered in the SSFS stability study for the classical NLS in [24], where stability on the background of traveling waves was studied, but did not need to be considered in the SSFS stability study on the background of a standing plane wave [25], because in that case S aliased is an empty set.
Carrying on and comparing (28) with (26) and their complex conjugates, we obtain the following equations: where we let l = m − k, and denote s as The matrix with The solution of system (31) will exponentially grow if either |Λ l,+ | > 1 or |Λ l,− | > 1. From (32), we then deduce the instability condition summarized below. Lemma 1 ( General instability condition). The unstable modes of the SSFS method in simulating a plane wave solution in (5) satisfy cos θ l − γ|a| 2 τ sin θ l > 1, (34) with θ l being defined in (33).
We remark that Lemma 1 predicts both analytically unstable modes as described in (13) and numerically unstable modes introduced by the SSFS method. In what follows, we assume that the time step τ is small enough such that which is necessary to ensure accuracy of the SSFS method. We then divide our discussion of the instability condition (34) into two parts. On one hand, with assumption (35) the condition cos θ l − γ|a| 2 τ sin θ l > 1 is satisfied only when |θ l | 1, implying that unstable modes in this case occur only in the main spectral region. Furthermore, the Taylor expansion of cos θ l − γ|a| 2 τ sin θ l > 1 up to terms O(τ 2 ) yields i.e., the same condition as in (13). Therefore, the condition cos θ l − γ|a| 2 τ sin θ l > 1 in (34) predicts the analytically (i.e., truly) unstable modes of the plane wave solutions.
On the other hand, the numerical instability condition of the SSFS method is given by which occurs when θ l ≈ π + 2nπ (for n ∈ Z). Then, we have the following results: Theorem 3 (Numerical stability condition). Assume that the wave number λ k µ max and the nonlinearity |γ||a| 2 τ 1. Then, the SSFS method remains numerically stable, provided that the time step satisfies where the ratio r = λ k /µ max , and the critical value α cr (r) < 1 is defined as satisfying it is plotted in Figure A1b. Note that the condition in the second line of (37), when it applies, sets a smaller stability threshold than that in the first line.
A detailed derivation of this Theorem can be found in Appendix A. The condition (37) implies that the larger the speed of the simulated plane wave, the larger the time threshold of numerical stability. In simulating the standing plane wave solutions, the stability condition for the SSFS method reduces to τ < π 1−α h α for any α ∈ (0, 2]. Corollary 1 (Classical NLS with α = 2). (i) For the classical NLS (α = 2), condition (37) reduces to in agreement with the conclusions in [25] for λ k = 0 and [24] for λ k = 0. (ii) For the practical case of λ k /µ max ≡ r 1, the first line of condition (37) reduces to .

Remark 1.
The condition (37) suggests that for α 1, the stability threshold of time step may become O(1). However, such a large time step cannot be used solely from the accuracy perspective. Therefore, in this case a reasonable requirement to impose on the time step is (35).
The analysis of this section can be generalized to higher-order split-step methods. We have numerically verified that the stability condition (37) remains the same for a fourthorder split-step method, although the growth rate of the numerical instability has been observed to be different. This sameness of the stability condition and the difference of the instability growth rates for the SSFS applied to the fractional NLS parallels that for the case of the classical NLS [32].

Dynamics of Plane Wave Solutions
Here, we study the dynamics of unstable plane wave solutions by directly simulating the NLS (1) with the first-order SSFS method. In selected cases, we have verified that using a fourth-order SSFS method does not affect the conclusions. Our main goal here is to highlight effects of the fractional, versus classical, Laplacian, on the dynamics of the plane wave solutions beyond the linear stage of the analytical instability development.
The initial condition is chosen as a perturbed plane wave solution of the following form: where a and ε 0 are constants, and |ε 0 | 1. Note that the shape of this perturbation has a maximum at the middle of the computational window and minima at its boundaries. Therefore, for the standing plane waves, this perturbation is expected to favor the emergence of a pulse in the middle of the computational domain. Perturbations that have maxima at different locations will favor formation of pulses there. This all pertains to the initial stage of the evolution; at later times, the pulses begin to move due to their interaction with many different Fourier harmonics created from the background noise, which thus generically have different magnitudes. Nonetheless, our simulations show that various initial perturbations of a similar size may do not affect our conclusions about differences between the fractional and classical NLS.
In the following, we choose ε 0 = 10 −5 and a = 1 2 in (39). Unless stated otherwise, we set domain length L = 5π with number of grid points J = 2 13 . In all cases considered, we selected the time step τ according to condition (37) so as to avoid numerical instability.

Standing Plane Waves
In the defocusing NLS, the analysis in Theorem 1 (i) shows that the standing plane waves are always stable for any α ∈ (0, 2], and our numerical simulations confirm this. Therefore, below we report results for the standing plane wave in the focusing NLS. We chose λ k = 0 in the initial condition (39) and γ = −1 in the NLS (1).
For the parameters stated in the preamble for this Section, formula (15) predicts that for any α ∈ (0.76, 2], there exists one pair of unstable modes with l = ±1 on the spectral grid µ l , and thus the standing plane wave is linearly unstable. Figure 4 presents the time evolution of |u(x, t)| for different α. It shows that the solution remains a plane wave for a short time, and then around t = 50 (which, of course, depends on the parameters chosen in (39)) one hump appears due to the instability. For α 2, the hump retains its shape only for a short period of time, after which it dissolves into a delocalized state (a plane wave for α = 2), then reassembles itself again, and this process goes on approximately periodically.
As α decreases, the time interval between the hump reappearances decreases, and both it and the amplitude of the re-emerging pulse become more and more irregular. As α → 1 + , the solution tends to form a narrow pulse (see Figure 4 bottom row). For α ∈ (0.76, 1], the plane wave was found to collapse. (For α < 0.76, there are no unstable modes according to (15), and thus the plane wave remains stable. The more analytically unstable modes exist in the computational domain (which can be controlled either by the domain's length L or by the wave's amplitude |a|), the earlier the collapse occurs. Furthermore, the modulational instability growth rate, found from Equations (11)- (13), has an evident effect on the time of the collapse occurrence.) The above observations about collapse are consistent with findings in [14,28], which also demonstrated collapse for α < 1, albeit in their case the initial condition was a localized solitary wave. It is also similar to what has been reported about solutions of the Majda-McLaughlin-Tabak model, a model with nonlocality in both dispersion and nonlinearity, which has been proposed to mimic wave turbulence [4,7,8].   In the NLS with α equal or close to 2, the spectrum expands (almost) periodically as the unstable modes evolve. The smaller the fractional order α, the more frequently the spectrum expands and shrinks, and also the broader the expanded spectrum becomes. When α is sufficiently close to 1, the spectrum expands quickly after the analytical instability is triggered, but does not shrink back (see Figure 5 for α = 1.1).
Let us note that for a given power α, a smaller nonlinearity (i.e., the value of |γ||a| 2 ) results in a less broad spectrum of the emerging pulse (cf. Figure 6 for γ = −0.8). Additionally, the time required for the formation of a narrow pulse increases as the nonlinearity decreases. The Fourier spectrum shown in Figure 5 illustrates the role of long-range dispersion that is described by the fractional Laplacian (−∆) α 2 . Due to an interplay of nonlinearity with long-range dispersion, the instability at low wave numbers triggers creation of new harmonics in the high wave number region. The smaller the fractional power α, the stronger the long-range interactions, and thus the stronger this effect. While the linear analysis predicts this instability, it could not predict that the instability would lead to a dramatic spectral expansion. Numerical simulations were required to reveal this nonlinear stage of the instability development. Our observations are consistent with those reported in the literature [19], where a discrete (with finite inter-site spacing) version of the NLS with long-range interactions was studied. As we noted in Section 1, the continuum limit of such a discrete equation leads to the fractional NLS. In fact, it is stated in [19] that reproducing their results within the framework of a continuous fractional NLS was then an open problem; we have addressed this problem here.
To conclude this section, we present, in Figure 7a,b, snapshots of the field and its Fourier spectrum at time instances shortly before the collapse occurs (the last snapshot is taken within 0.1 time unit from the instance when the Fourier spectrum becomes essentially flat). A feature to note is the cusp in the spectrum at µ = 0. This signals a slow decay of the field u with x, which is due to the long-range interactions enabled by the Lévy index α = 0.9.

Effects of a Larger Domain
In this section, we focus on the dynamics of standing plane waves in the focusing NLS when the domain length L increases. Here, we use γ = −1 and L = 10π, instead of L = 5π used in the previous subsection. By formula (15), more unstable modes appear as the domain size increases. Specifically, in this case we find that there exist three unstable pairs (i.e., l = ±1, ±2, ±3) for any α ∈ [1.36, 2]. Their interaction with one another can result in an emergence of more than one pulse; this expectation is indeed borne out by our numerical simulations, reported in Figure 8. The field dynamics in the classical NLS (α = 2) appears to be regular. This is likely the consequence of this equation's integrability and the presence of only a small number of unstable modes, as the computational domain is still relatively short. As α decreases, the evolution of the initially emerged "humps" in the field (around t = 100) becomes more an more chaotic. This may be explained by the fact that in nonintegrable systems, pulse interaction is usually accompanied by energy transfer from a smaller to a larger pulse (see, e.g., in [33]). Eventually, such repeated interactions lead to formation of a single dominant pulse propagating in a "sea" of small-amplitude dispersive waves. Another feature, most evident from the lower panels of Figure 8, is that a spatially symmetric initial condition results in the formation of a moving pulse (i.e., an asymmetric one). This can be explained by the fact that since the linearly unstable modes, which grow initially, arise from a background noise due to the computer round-off error and thus have different amplitudes (i.e., are asymmetric).
For α ≤ 1, as in the case of the shorter domain, the solution eventually collapses as long as linearly unstable modes exist.

Effects of Nonzero Group Velocity (7) of the Plane Wave
The analysis in Theorem 1 shows that the traveling (λ k = 0) plane waves in the defocusing NLS are always stable if the fractional power α ∈ [1, 2], which is confirmed by our numerical simulations. Therefore, we will focus our following study on the traveling plane wave in the focusing NLS with α ∈ (0, 2] or in the defocusing NLS with α ∈ (0, 1). Figure 9 shows the time evolution of |u(x, t)| for three representative values of α in the focusing NLS when α > 1, where λ k = 2π/L (i.e., k = 1) in the initial condition (39). Equation (13) predicts that there is one pair of unstable modes at l = ±1 for α = 2 and 1.9. Correspondingly, the instability develops into a single soliton-like pulse that moves over a nonzero background. In the classical NLS, the reappearance of plane wave solution is observed, and the frequency of this recurrence is almost the same as that observed in Figure 4 for standing plane waves. Similar to the observations in Figure 4, a small decrease of α from the classical value α = 2 leads to a more frequent appearance of the pulse out of the quiescent state. While for α = 1.3, two unstable pairs occur at l = ±1, ±2, and the behavior of the solution becomes more chaotic due to a nonlinear interaction of these modes (initially) and others excited by them (eventually).
Our simulations show that traveling plane waves in the focusing NLS with α ∈ (0, 1] also collapse, but the occurrence of collapse is progressively delayed by the greater values of λ k (with all other parameters being fixed). Table 1 shows the collapse time in the focusing cases, where the collapse time t col is defined as the instance where the magnitude of the Fourier component at the edge of the computational spectral window reaches 10 −5 . While this number is, of course, chosen arbitrarily, replacing it with, say, 10 −3 or 10 −7 would cause only an insignificant change in the numbers reported in Table 1 and would not change the trends observed. This is because, as earlier for the case of the standing plane wave, collapse in the x-space corresponds to the continuing broadening of the solution's Fourier spectrum. Once this spectrum increases by several orders of magnitude at the edges of the computational spectral window, the solution can no longer be considered valid. Of course, we verified that doubling the width of the computational spectral window would not prevent the solution's spectrum from continuing to broaden; therefore, doubling that width only insignificantly affected the collapse times. We find that for a given λ k , collapse of solution appears earlier for a smaller α. This is the opposite of the situation with the standing wave's collapse and is explained by the magnitude of the modulational instability growth rates in the two cases. Indeed, for the standing wave, this rate decreases with decreasing α (see Figure 2) and disappears for some small enough α (see the discussion in Section 4.1), it retains its value for a traveling wave as α decreases (see Figure 3). On the other hand, for a given α, the collapse is delayed by increasing λ k . We hypothesize that this may occur because for a larger λ k , unstable modes µ l appear at larger |l| (see Figure 3); a possible causal link is discussed at the end of next paragraph.
Snapshots of the solution and its spectrum evolving towards collapse are shown in Figure 10a,b.  Figure 10. Same as Figure 7, but for λ k = 5 · (2π/L). At all simulated times before the collapse, the magnitude of the Fourier spectrum at the edge of the computational spectral window is on the order of 10 −12 . The fact that in panel (a) the "hump" at t = 286 and t = 289 appears at approximately the same location is a coincidence; this pulse moves around the computational window.
The last snapshot is taken within 0.2 time unit before the collapse occurs. A few observations follow from this figure. First, similar to Figure 7b, the spectrum has a cusp at the top, indicating some long-range "order" present in the solution, in agreement with the presence of long-range interactions due to the fractional Laplacian with a low value of α. Second, the following scenario of collapse can be deduced. The initial small corrugation grows to a finite size and develops a so-called dispersive wave region, seen as faster oscillations on the right side of each of the periods of the sinusoid for t = 156. Next, due to these faster oscillations moving through the slower corrugation, there form "humps" with slightly different amplitudes (t = 278). Due to their interaction with one another, one dominant "hump" forms (t = 286). This is similar to the behavior observed in other non-integrable systems [33], where the dominant pulse "sucks in" the energy from the ambient field and smaller pulses via collisions with them. It is this dominant pulse that subsequently undergoes collapse. Finally, this scenario explains the reason behind the delay of the collapse caused by a nonzero λ k of the initial wave. Apparently, collapse requires that a certain "amount" of the field be concentrated near a given location. As long as there is a corrugation with peaks of the (approximate) sinusoid of similar amplitude, collapse cannot form. However, as soon as "enough" energy has been absorbed by the dominant pulse, collapse becomes possible. The larger λ k , the smaller the period of the corrugation, and thus the smaller the energy within each one of its periods, and the longer it takes for a sufficient amount of energy to amass near a given location to begin the onset of collapse.
Two more remarks are in order. First, we stress that this ability of the smaller period of the initial corrugation to delay the onset of collapse was stated for the condition that α is fixed. Indeed, as one takes α so small that the results of Section 2 (see Figure 2) predict a wide band of modulationally unstable modes (e.g., for α = 0.15), the corrugation at the initial stage of the evolution is quite irregular (due to many modes growing at a similar rate). Yet, the fact that α is smaller, and hence the dispersion term (−∆) α/2 has a significantly longer range, acts to overcome the effect of this irregularity and promotes collapse; the latter is found to occur significantly earlier than for α = 0.9.
Second, while the smaller period of the initial corrugation can only delay, but not arrest, collapse, a recent study [34] has shown that a trapping potential can arrest it. This is consistent with the fact that collapse is caused by the interaction of the field over a sufficiently large range, and a suitable trapping potential is able to essentially truncate this range by localizing the field.
To conclude this subsection, we briefly comment on the defocusing NLS with α ∈ (0, 1). In this case, although the traveling plane waves are also unstable, their dynamics are substantially different from those in the focusing case. No collapse of solution is observed, and the evolution of the solution is qualitatively similar to that in the focusing NLS with α > 1, shown in Figure 9.

Conclusions
In this paper, we investigated the stability and dynamics of the plane wave solutions of the fractional NLS from both analytical and numerical perspectives. The comparison between the classical and fractional NLS were made to understand the role of long-range dispersion in the fractional NLS.
First, in Section 2, we showed that modulational stability of the plane wave solutions to small perturbations can be significantly different in the fractional NLS (1) compared to its classical counterpart when the Lévy index α ∈ (0, 1). In the defocusing (γ > 0) NLS, the standing plane wave solution (i.e., λ k = 0) is always stable for any α ∈ (0, 2]. In contrast, the traveling plane wave solution (i.e., λ k = 0) is stable only for α ∈ [1,2]. In the focusing (γ < 0) NLS, the plane wave solution can be linearly unstable for any α and λ k . Furthermore, unstable modes occur at low-frequency region for α ∈ [1, 2], but they appear at relatively high-frequency region for α ∈ (0, 1), which implies that the dynamics of the plane wave solution (before they collapse) for α ∈ (0, 1) could be more chaotic.
Second, in Section 3, we studied stability of the split-step Fourier spectral (SSFS) method in simulating the dynamics of the plane wave solution of the fractional NLS. This is essential in order to guarantee that the numerical solution does not suffer a blow-up due to purely numerical artifacts; then, one can ascertain that, if the solution undergoes collapse, it reflects the true behavior of the analytical solution. We obtained an expression for the time step threshold below which numerical instability does not occur. This threshold depends on both the parameter α, of the NLS, and λ k , of the simulated plane wave. In particular, for smaller α, this threshold can be quite large, and then the time step should be chosen to satisfy the requirement of accuracy, rather than stability, of the numerical method.
Third, in Section 4, we studied the nonlinear stage of the plane wave instability by directly simulating the NLS. For the initial condition being a perturbed standing plane wave in the focusing NLS, we found that the perturbation generically develops into a localized pulse on a nonzero background; the lower the Lévy index α > 1, the more localized the pulse. In the Fourier space this corresponds to of a significant widening of the solution spectrum, which expands well beyond the band of modulationally unstable modes found by the linear stability analysis. For α > 1, the pulse undergoes "breathing", whereby it narrows and broadens intermittently. For α ≤ 1, the pulse localization (in x-space) eventually leads to its collapse, in agreement with findings in [14,28].
We also found that the increase of the group velocity of the plane wave (so that the latter is traveling as opposed to standing) in the focusing NLS delays the formation of a localized pulse (or pulses), for a given fractional index < 1α ≤ 2, and delays the onset of collapse, for α < 1. To our knowledge, this has not been reported before. We also presented a hypothesis as to why this delay could occur. In the defocusing case, even though traveling plane waves are found to be linearly unstable for α < 1, that instability does not not lead to collapse of the solution.

Conflicts of Interest:
The authors declare no conflict of interest.
Due to the assumption (35), one has φ = −γ|a| 2 τ + O (γ|a| 2 τ) 3 . (A1) Then we can rewrite the numerical instability condition in (36) as which can be satisfied at θ l ≈ nπ, for integer |n| ≥ 1. Here, we only consider two possibilities: θ l ≈ π or −π, as the others (i.e., θ l ≈ nπ with |n| > 1) will yield τ that is larger than the stability threshold. The derivation of stability condition is slightly different for focusing and defocusing NLS. Below, we will present a detailed derivation for the focusing NLS and will briefly comment on the situation with the defocusing case at the end.
The two possibilities: θ l ≈ π and θ l ≈ −π, need to be considered separately. (i) For θ l ≈ π, the condition (A2) implies (recall that γ < 0 and thus φ > 0): For α > 1, the situation is illustrated by the solid curve in Figure A1a, which shows that the maximum value θ l,max occurs at the edges of the spectral regions (30). Taking then l at the edge of S main , so that |µ l | = µ max − λ k (see Figure A1a), yields θ l,max = θ l (µ max − λ k ). Using (33) and substituting the result into the left inequality in (A3), one obtains the first line of the stability condition (37) with the first argument of the "max" operator. The second argument of "max" comes into play when one considers the case α < 1, where, as illustrated by the dashed line in Figure A1a, θ l,max occurs in the middle of the aliased region: θ l,max = θ l (µ min + λ k ). (ii) For θ l ≈ −π, the condition (A2) implies − π < θ l < −π + 2φ. (A4) Note that the condition (A4) can be satisfied only when α < 1 and λ k = 0. In other words, the stability condition of the SSFS method in simulating the standing plane wave is always given by (A3). For α < 1 and λ k = 0, the minimum value θ l,min is attained for µ l = ±λ k rather than for µ l = 0 as for α > 1; see the dashed line in Figure A1a. It will be the condition (A4) rather than (A3) that will set the stability threshold when |θ l, min | ≥ |θ l, max |, or equivalently, |θ l (±λ k )| ≥ |θ l (µ max − λ k )|. Using (33), this yields The critical value α cr defined in Theorem 3 follows from (A5) with the equal sign and is plotted in Figure A1b. For α < α cr , substituting θ l,min into the right inequality of (A4) leads to the stability condition By assumption (35), the value of 2φ is negligible compared to π, and thus one obtains the second line of stability condition (37). Note that if condition (A6) is violated, numerical instability will occur at |µ l | ≈ λ k , or, equivalently, for µ k ≈ 0, 2λ k . Therefore, it will more readily contaminate the solution than an instability occurring for |µ l | ≈ µ max , considered previously. However, in practice such an instability can occur only for rather large time steps, which should be avoided based on the requirement of accuracy rather than stability of the numerical method.  (33) with τ exceeding the instability threshold (37) by 10%. Note that µ m ≡ µ l + λ k , as noted in Section 3.2. For smaller ratios λ k /µ max the relative size of the aliased region is smaller. (b) Relative difference |θ l 1 =k − θ l 1 =0 |/(|θ l 1 =k + θ l 1 =0 |/2) in the aliased region.
In the defocusing (γ > 0) case, we will only comment the difference in the analysis for θ l ≈ π. Following a similar analysis that led to condition (A3) in the main spectral region, we obtain the following condition of numerical instability: where we recall that φ, given by (A1), is now negative. The relative effect caused by the φ-term in (A7) is on the order O(|γ||a| 2 ) and thus, by assumption (35), can be neglected. In this approximation, the numerical stability condition of the SSFS method is still given by (37).