Open AccessArticle

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

Siwei Duo

¹,

Taras I. Lakoba

and

Yanzhi Zhang

^3,*

Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA

Department of Mathematics and Statistics, University of Vermont, Burlington, VT 05405, USA

Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409, USA

Author to whom correspondence should be addressed.

Symmetry 2021, 13(8), 1394; https://doi.org/10.3390/sym13081394

Submission received: 21 June 2021 / Revised: 20 July 2021 / Accepted: 22 July 2021 / Published: 1 August 2021

(This article belongs to the Special Issue Recent Progress in Studies of Stability of Numerical Schemes)

Download

Browse Figures

Versions Notes

Abstract

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

α \in [1, 2]

but unstable for

α \in (0, 1)

. In the focusing case, they can be linearly unstable for any

α \in (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 find that as

α \in (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

α \leq 1

, we always find the solution to undergo collapse. We also show, for the first 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 find 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.

Keywords:

fractional nonlinear Schrödinger equation; fractional Laplacian; plane wave solution; modulation instability; split-step method; numerical stability

1. 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]:

\begin{matrix} i \partial_{t} u (x, t) = {(- Δ)}^{\frac{α}{2}} u + γ {| u |}^{2} u, x \in R, t > 0, \end{matrix}

(1)

where

i = \sqrt{- 1}

is the imaginary unit, and

u (x, t)

is a complex-valued wave function. The constant

γ \in 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

{(- Δ)}^{\frac{α}{2}}

is defined via a pseudo-differential operator [5,15,16,17]:

\begin{matrix} {(- Δ)}^{\frac{α}{2}} u (x, t) : = F^{- 1} [{| μ |}^{α} F (u)], f o r α > 0, \end{matrix}

(2)

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 < α \leq 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,

\begin{matrix} N (t) : = \int_{R} {| u (x, t) |}^{2} d x = \int_{R} {| u (x, 0) |}^{2} d x = N (0), t \geq 0, \end{matrix}

(3)

and the total energy,

\begin{matrix} E (t) : = \int_{R} [| \nabla^{\frac{α}{2}} u |^{2} + \frac{γ}{2} {| u |}^{4}] d x = E (0), t \geq 0, \end{matrix}

(4)

where

\nabla^{α} = - {(- Δ)}^{\frac{α}{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.

2. Plane Wave Solution and Its Stability

The fractional NLS (1) admits the plane wave solution of the form

\begin{matrix} u (x, t) = a exp (i (λ_{k} x - ω t)), \end{matrix}

(5)

provided that the following nonlinear dispersion relation is satisfied

\begin{matrix} = | λ_{k} |^{α} + γ {| a |}^{2} . \end{matrix}

(6)

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} \neq 0

would move with a group velocity

\begin{matrix} v_{g} = \frac{\partial}{\partial λ_{k}} = α sgn (λ_{k}) {| λ_{k} |}^{α - 1}, for λ_{k} \neq 0, \end{matrix}

(7)

where

sgn (\cdot)

denotes the sign function. For this reason, we will refer to a plane wave with

λ_{k} \neq 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

\begin{matrix} λ_{k} = 2 k π / L, k \in Z . \end{matrix}

Without loss of generality, we will assume

λ_{k} \geq 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

α \in (0, 2]

. Consider a perturbed solution of the form

\begin{matrix} \tilde{u} (x, t) = u (x, t) (1 + ε (x, t)), \end{matrix}

(8)

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

\begin{matrix} ε (x, t) = \sum_{l \in Z} {\hat{ε}}_{l} (t) exp (i μ_{l} x) \end{matrix}

with

μ_{l} = 2 l π / L

, for

l \in Z

. In what follows, we set

Re ({\hat{ε}}_{0} (0)) = 0

without loss of generality. Indeed, having

Re ({\hat{ε}}_{0} (0)) \neq 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

\partial_{t} ε (x, t) = i (| λ_{k} |^{α} ε - \sum_{l \in Z} {\hat{ε}}_{l} (t) | λ_{k} + μ_{l} |^{α} exp (i μ_{l} x) - γ {| a |}^{2} (ε^{*} + ε)),

(9)

where

ε^{*}

represents the complex conjugate of

ε

. Taking the Fourier transform of (9) and its complex conjugate, we get the following system of ODEs:

\frac{d}{d t} (\begin{matrix} {\hat{ε}}_{l} \\ {\hat{ε}}_{- l}^{*} \end{matrix}) = G_{l} (\begin{matrix} {\hat{ε}}_{l} \\ {\hat{ε}}_{- l}^{*} \end{matrix}), l \in Z ∖ {0},

(10)

where

\begin{matrix} G_{l} = i (\begin{matrix} | λ_{k} |^{α} - {γ | a |}^{2} - {| λ_{k} + μ_{l} |}^{α} & - {γ | a |}^{2} \\ {γ | a |}^{2} & - | λ_{k} |^{α} + {γ | a |}^{2} + {| λ_{k} - μ_{l} |}^{α} \end{matrix}) . \end{matrix}

The eigenvalues of matrix

G_{l}

are

\begin{matrix} Λ_{l}^{\pm} = - \frac{i}{2} (| λ_{k} + μ_{l} |^{α} - {| λ_{k} - μ_{l} |}^{α}) \pm \frac{1}{2} \sqrt{Δ_{l}}, \end{matrix}

(11)

where

Δ_{l} = - [(| λ_{k} + μ_{l} |^{α} + {| λ_{k} - μ_{l} |}^{α}) - 2 {| λ_{k} |}^{α}] [(| λ_{k} + μ_{l} |^{α} + {| λ_{k} - μ_{l} |}^{α}) - 2 (| λ_{k} |^{α} - 2 γ {| a |}^{2})] .

(12)

The solution of system (10) grows exponentially if at least one of the eigenvalues

Λ_{l}^{\pm}

has positive real part, which occurs if

\begin{matrix} Δ_{l} > 0 . \end{matrix}

(13)

In other words, if the condition (13) holds for at least one

l \in 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 $α \in (0, 2]$ , a standing plane wave is always stable to infinitesimal perturbations.
(ii): For $α \in [1, 2]$ , a traveling plane wave is always stable, independent of the wave number $λ_{k}$ .
(iii): For $α \in (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:

$| λ_{k} |^{α} < \frac{{4 γ | a |}^{2}}{2 - 2^{α}}, α \in (0, 1) .$

(14)

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

α \in [1, 2]

, traveling plane waves are always stable, confirming the conclusion in Theorem 1 (ii). By contrast, they are unstable for

α \in (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} = \pm λ_{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

\pm λ_{k}

, which are surrounded by unstable modes (these stable bands are clearly visible in Figure 1d for

α \in (0, 0.735)

Theorem 2

(Focusing NLS with

γ < 0

(i): For $α \in (0, 2]$ , a standing plane wave is unstable, provided that there is at least one $l \in Z ∖ {0}$ such that $| μ_{l} |^{α} {< 2 | γ | | a |}^{2}$ . Moreover, the number of unstable pairs is given by

$n_{p} = ⌊\frac{{L (2 | γ | | a |}^{2})^{1 / α}}{2 π}⌋, f o r α \in (0, 2],$

(15)

where $⌊ \cdot ⌋$ denotes the floor function.
(ii): For $α \in (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

α \in (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

α \in (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

α \in (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

α \in (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

α \in (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

| μ | \leq μ_{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

Δ_{l} \approx - 4 (| μ_{l} |^{α} - λ_{k}^{α}) (| μ_{l} |^{α} - λ_{k}^{α} + 2 γ {| a |}^{2}) .

(16)

Note that despite condition

| μ_{l} | ≫ λ_{k}

, the term of

λ_{k}

cannot be neglected if

α

is sufficiently small.

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:

\begin{matrix} | μ_{l} | < {(λ_{k}^{α} + {2 | γ | | a |}^{2})}^{1 / α} . \end{matrix}

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.

3. 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.

3.1. Description of the SSFS Method

Let

τ > 0

be the time step, and define

t_{n} = n τ

for

n = 0, 1, \dots

From time

t = t_{n}

t = t_{n + 1}

, one numerically solves the NLS (1) in two steps:

\begin{matrix} i u_{t} (x, t) = γ {| u |}^{2} u, \end{matrix}

(17)

\begin{matrix} i u_{t} (x, t) = {(- Δ)}^{\frac{α}{2}} u . \end{matrix}

(18)

For

t \in (t_{n}, t_{n + 1}]

, Equation (17) conserves the quantity

| u (x, t) | \equiv | 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):

\begin{matrix} u (x, t) = u (x, t_{n}) e^{- i γ | u (x, t_{n}) |^{2} (t - t_{n})}, t \in [t_{n}, t_{n + 1}] . \end{matrix}

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

x_{j} = - L / 2 + j h

, for

0 \leq j \leq J

. We take

\begin{matrix} u (x_{j}, t) = \sum_{m = - J / 2}^{J / 2 - 1} {\hat{u}}_{m} (t) e^{i μ_{m} x_{j}}, \end{matrix}

(19)

where the m-th Fourier mode

{\hat{u}}_{m} (t)

is defined as

\begin{matrix} {\hat{u}}_{m} (t) = \frac{1}{J} \sum_{j = 0}^{J - 1} u (x_{j}, t) e^{- i μ_{m} x_{j}}, - \frac{J}{2} \leq m \leq \frac{J}{2} - 1 . \end{matrix}

Substituting (19) into (18), taking the discrete Fourier transform, and integrating the resulting ordinary differential equation for

{\hat{u}}_{m} (t)

, one obtains

\begin{matrix} {\hat{u}}_{m} (t) = {\hat{u}}_{m} (t_{n}) e^{- i | μ_{m} |^{α} (t - t_{n})}, t \in [t_{n}, t_{n + 1}], - \frac{J}{2} \leq m \leq \frac{J}{2} - 1 . \end{matrix}

(20)

Combining (20) with (19) leads to an approximate solution of (18).

Finally, denoting the numerical approximation of

u (x_{j}, t_{n})

u_{j}^{n}

, one obtains the first-order SSFS method as

\begin{matrix} u_{j}^{n + 1, -} = u_{j}^{n} e^{- i γ | u_{j}^{n} |^{2} τ}, \end{matrix}

(21)

\begin{matrix} u_{j}^{n + 1} = \sum_{m = - J / 2}^{J / 2 - 1} {\hat{u}}_{m}^{n + 1, -} e^{- i | μ_{m} |^{α} τ} e^{i μ_{m} x_{j}}, 0 \leq j \leq J - 1, n = 0, 1, \dots \end{matrix}

(22)

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

α \in (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.

3.2. 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

α \neq 2

and

λ_{k} \neq 0

. Here, we assume

λ_{k} ≪ μ_{max},

(23)

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:

u_{j}^{n} = {\tilde{u}}_{j}^{n} (1 + ε_{j}^{n}), with {\tilde{u}}_{j}^{n} = a e^{i (λ_{k} x_{j} - ω t_{n})},

(24)

where

ε_{j}^{n}

, with

| ε_{j}^{n} | ≪ 1

, is a small periodic perturbation on

[- L / 2, L / 2]

and can be expanded as:

ε_{j}^{n} = \sum_{l = - J / 2 - k}^{J / 2 - 1 - k} {\hat{ε}}_{l}^{n} e^{i μ_{l} x_{j}}, 0 \leq j \leq J - 1 .

(25)

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} \neq 0

from the by-now standard analysis for

λ_{k} = 0

[25]. In fact, the perturbation term in (24) is

{\tilde{u}}_{j}^{n} ε_{j}^{n}

, which has effective wave number of

(λ_{k} + μ_{l}) \equiv μ_{k + l}

rather than just

μ_{l}

. To ensure this wave number fits into the range of

[μ_{min}, μ_{max}]

, we have to require

μ_{l} \in [μ_{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,

{\hat{u}}_{m}

, and the multiplicative perturbation,

{\hat{ε}}_{l}

Setting

n = n + 1

in (24) and taking its discrete Fourier transform, one obtains

\begin{matrix} {\hat{u}}_{m}^{n + 1} = a e^{- i t_{n + 1}} \{\begin{matrix} 1 + {\hat{ε}}_{0}^{n + 1}, & m = k, \\ {\hat{ε}}_{m - k}^{n + 1}, & m \neq k, \end{matrix} - \frac{J}{2} \leq m \leq \frac{J}{2} - 1 . \end{matrix}

(26)

On the other hand,

{\hat{u}}_{m}^{n + 1}

can be independently found from the SSFS scheme (21)–(22). To this end, we substitute

u_{j}^{n}

from (24) into (21) and take the leading order terms of

ε_{j}^{n}

to obtain

u_{j}^{n + 1, -} = {\tilde{u}}_{j}^{n} e^{- {i γ | a |}^{2} τ} (1 + ε_{j}^{n} - i γ {| a |}^{2} τ (ε_{j}^{n} + {(ε_{j}^{n})}^{*})), 0 \leq j \leq J - 1 .

(27)

Substituting (27) into (22) and taking its discrete Fourier transform, we obtain

\begin{matrix} {\hat{u}}_{m}^{n + 1} = a e^{- i (t_{n} + γ | a |^{2} τ + | μ_{m} |^{α} τ)} \{\begin{matrix} 1 + {\hat{ε}}_{0}^{n}, & m = k; \\ {\hat{ε}}_{m - k}^{n} - i γ {| a |}^{2} τ ({\hat{ε}}_{m - k}^{n} + {({\hat{ε}}_{k - m}^{n})}^{*}), & (m - k) \in 𝒮_{main}; \\ {\hat{ε}}_{m - k}^{n} - i γ {| a |}^{2} τ ({\hat{ε}}_{m - k}^{n} + {({\hat{ε}}_{k - m - J}^{n})}^{*}), & (m - k) \in 𝒮_{aliased} . \end{matrix} \end{matrix}

(28)

For

m = k

, we have used the convention Re

({\hat{ε}}_{0}^{n}) = 0

, adopted in Section 2. Furthermore, we have split the wave number indices

l \equiv (m - k)

of the perturbation terms

{\hat{ε}}_{m - k}^{n}

into the “main” and “aliased” sets defined as

\begin{matrix} 𝒮_{main} = \{l | - (J / 2 - 1 - k) \leq l \leq J / 2 - 1 - k\}, \end{matrix}

(29)

\begin{matrix} 𝒮_{aliased} = \{l | - J / 2 - k \leq l \leq - J / 2 + k\} . \end{matrix}

(30)

To understand this choice, note that when

(m - k) \equiv l

is within one of these sets, then the indices of both

\hat{ε}

- and

{\hat{ε}}^{*}

-terms in each line in (28) are within the same set. For

𝒮_{main}

, this is evident, while for

𝒮_{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

𝒮_{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

𝒮_{aliased}

is an empty set.

Carrying on and comparing (28) with (26) and their complex conjugates, we obtain the following equations:

(\begin{matrix} {\hat{ε}}_{l}^{n + 1} \\ {({\hat{ε}}_{s}^{n + 1})}^{*} \end{matrix}) = A_{l} (\begin{matrix} {\hat{ε}}_{l}^{n} \\ {({\hat{ε}}_{s}^{n})}^{*} \end{matrix}), for - \frac{J}{2} - k \leq l \leq \frac{J}{2} - 1 - k,

(31)

where we let

l = m - k

, and denote s as

\begin{matrix} s = \{\begin{matrix} - l, & if l \in 𝒮_{main}; \\ - (l + J), & if l \in 𝒮_{aliased} . \end{matrix} \end{matrix}

The matrix

\begin{matrix} A_{l} & = & (\begin{matrix} {(1 - i γ | a |}^{2} τ) e^{- i (| μ_{k + l} |^{α} - | λ_{k} |^{α}) τ} & - {i γ | a |}^{2} τ e^{- i (| μ_{k + l} |^{α} - | λ_{k} |^{α}) τ} \\ {i γ | a |}^{2} τ e^{i (| μ_{k + s} |^{α} - | λ_{k} |^{α}) τ} & {(1 + i γ | a |}^{2} τ) e^{i (| μ_{k + s} |^{α} - | λ_{k} |^{α}) τ} \end{matrix}) \end{matrix}

has eigenvalues

Λ_{l, \pm} = (cos θ_{l} - γ {| a |}^{2} τ sin θ_{l}) \pm \sqrt{{(cos θ_{l} - γ {| a |}^{2} τ sin θ_{l})}^{2} - 1}

(32)

with

θ_{l} = \frac{τ}{2} (| μ_{l} + λ_{k} |^{α} + | μ_{s} + λ_{k} |^{α} - 2 {| λ_{k} |}^{α}) .

(33)

The solution of system (31) will exponentially grow if either

| Λ_{l, +} | > 1

| Λ_{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

{| γ | | a |}^{2} τ ≪ 1,

(35)

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

\begin{matrix} (| μ_{l} + λ_{k} |^{α} + | μ_{l} - λ_{k} |^{α} - 2 {| λ_{k} |}^{α}) (| μ_{l} + λ_{k} |^{α} + | μ_{l} - λ_{k} |^{α} - 2 (| λ_{k} |^{α} - {2 γ | a |}^{2})) < 0, \end{matrix}

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

(cos θ_{l} - γ {| a |}^{2} τ sin θ_{l}) < - 1,

(36)

which occurs when

θ_{l} \approx π + 2 n π

(for

n \in 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

τ < \{\begin{matrix} \frac{π^{1 - α} h^{α}}{max \{(1 + {(1 - 2 r)}^{α}) / 2, {(1 - r)}^{α}\} - r^{α}}, & i f r = 0, o r r \neq 0 a n d α > α_{cr} (r) \\ \frac{π^{1 - α} h^{α}}{(1 - 2^{α - 1}) r^{α}}, & i f r \neq 0 a n d α \leq α_{cr} (r) \end{matrix}

(37)

where the ratio

r = λ_{k} / μ_{max}

, and the critical value

α_{cr} (r) < 1

is defined as satisfying

2 - 2^{α - 1} = {((1 - r) / r)}^{α};

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

α \in (0, 2]

Corollary 1

(Classical NLS with

α = 2

(i): For the classical NLS ( $α = 2$ ), condition (37) reduces to

$\begin{matrix} τ < \frac{h^{2}}{π {(1 - λ_{k} / μ_{max})}^{2}}, \end{matrix}$

(38)

in agreement with the conclusions in [25] for $λ_{k} = 0$ and [24] for $λ_{k} \neq 0$ .
(ii): For the practical case of $λ_{k} / μ_{max} \equiv r ≪ 1$ , the first line of condition (37) reduces to

$τ < \frac{π^{1 - α} h^{α}}{(1 - α r) - r^{α} + O (r^{2})} .$

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 fourth-order 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].

4. 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:

u (x, 0) = a e^{i λ_{k} x} (1 + ε_{0} cos \frac{2 π x}{L}), - \frac{L}{2} \leq x \leq \frac{L}{2},

(39)

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 = \frac{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.

4.1. Standing Plane Waves

In the defocusing NLS, the analysis in Theorem 1 (i) shows that the standing plane waves are always stable for any

α \in (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

α \in (0.76, 2]

, there exists one pair of unstable modes with

l = \pm 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

α \to 1^{+}

, the solution tends to form a narrow pulse (see Figure 4 bottom row). For

α \in (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].

Figure 5 displays the time evolution of the spectrum of the solution

u (x, t)

, with the corresponding density dynamics shown in Figure 4.

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

{(- Δ)}^{\frac{α}{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

4.2. 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 = \pm 1, \pm 2, \pm 3

) for any

α \in [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

α \leq 1

, as in the case of the shorter domain, the solution eventually collapses as long as linearly unstable modes exist.

4.3. Effects of Nonzero Group Velocity (7) of the Plane Wave

The analysis in Theorem 1 shows that the traveling (

λ_{k} \neq 0

) plane waves in the defocusing NLS are always stable if the fractional power

α \in [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

α \in (0, 2]

or in the defocusing NLS with

α \in (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 = \pm 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 = \pm 1, \pm 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

α \in (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}

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.

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

α \in (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.

5. 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

α \in (0, 1)

. In the defocusing (

γ > 0

) NLS, the standing plane wave solution (i.e.,

λ_{k} = 0

) is always stable for any

α \in (0, 2]

. In contrast, the traveling plane wave solution (i.e.,

λ_{k} \neq 0

) is stable only for

α \in [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

α \in [1, 2]

, but they appear at relatively high-frequency region for

α \in (0, 1)

, which implies that the dynamics of the plane wave solution (before they collapse) for

α \in (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

α \leq 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 α \leq 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.

Author Contributions

Investigation, T.I.L., S.D. and Y.Z. All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

Y. Zhang was supported by the NSF grants DMS-1217000 and DMS-1620465, and T. I. Lakoba was supported by the NSF grant DMS-1217006.

Institutional Review Board Statement

Not Applicable.

Informed Consent Statement

Not Applicable.

Data Availability Statement

Not Applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Derivation of Theorem 3

In the following, we will provide a detailed derivation for the numerical stability condition (37) in Theorem 3. First, we define the quantity

ϕ

such that

\begin{matrix} cos ϕ = \frac{1}{\sqrt{1 + γ^{2} {| a |}^{4} τ^{2}}}, sin ϕ = \frac{- {γ | a |}^{2} τ}{\sqrt{1 + γ^{2} {| a |}^{4} τ^{2}}} . \end{matrix}

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

cos (θ_{l} - ϕ) < - cos ϕ,

(A2)

which can be satisfied at

θ_{l} \approx n π

, for integer

| n | \geq 1

. Here, we only consider two possibilities:

θ_{l} \approx π

- π

, as the others (i.e.,

θ_{l} \approx 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} \approx π

and

θ_{l} \approx - π

, need to be considered separately.

(i): For $θ_{l} \approx π$ , the condition (A2) implies (recall that $γ < 0$ and thus $ϕ > 0$ ):

$π < θ_{l} < π + 2 ϕ .$

(A3)

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 $𝒮_{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} \approx - π$ , the condition (A2) implies

$- π < θ_{l} < - π + 2 ϕ .$

(A4)

Note that the condition (A4) can be satisfied only when $α < 1$ and $λ_{k} \neq 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} \neq 0$ , the minimum value $θ_{l, min}$ is attained for $μ_{l} = \pm λ_{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} | \geq | θ_{l, max} |$ , or equivalently, $| θ_{l} (\pm λ_{k}) | \geq | θ_{l} (μ_{max} - λ_{k}) |$ . Using (33), this yields

$2 - 2^{α - 1} \geq {((1 - r) / r)}^{α} .$

(A5)

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

$τ < \frac{2 (π - 2 ϕ)}{2 λ_{k}^{α} - {(2 λ_{k})}^{α}} .$

(A6)

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} | \approx λ_{k}$ , or, equivalently, for $μ_{k} \approx 0, 2 λ_{k}$ . Therefore, it will more readily contaminate the solution than an instability occurring for $| μ_{l} | \approx μ_{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.

Figure A1. (a) Illustration of the main and aliased regions for

γ = - 1

λ_{k} / μ_{max} = 0.25

, where

θ_{l}

is defined in (33) with

τ

exceeding the instability threshold (37) by 10%. Note that

μ_{m} \equiv μ_{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.

Figure A1. (a) Illustration of the main and aliased regions for

γ = - 1

λ_{k} / μ_{max} = 0.25

, where

θ_{l}

is defined in (33) with

τ

exceeding the instability threshold (37) by 10%. Note that

μ_{m} \equiv μ_{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} \approx π

. Following a similar analysis that led to condition (A3) in the main spectral region, we obtain the following condition of numerical instability:

π + 2 ϕ < θ_{l} < π,

(A7)

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).

References

Laskin, N. Fractional quantum mechanics and Lévy path integrals. Phys. Lett. A 2000, 268, 298–305. [Google Scholar] [CrossRef] [Green Version]
Laskin, N. Fractional Schrödinger equation. Phys. Rev. E 2002, 66, 056108. [Google Scholar] [CrossRef] [Green Version]
Kirkpatrick, K.; Lenzmann, E.; Staffilani, G. On the continuum limit for discrete NLS with long-range lattice interactions. Comm. Math. Phys. 2013, 317, 563–591. [Google Scholar] [CrossRef] [Green Version]
Cai, D.; Majda, A.J.; McLaughlin, D.W.; Tabak, E.G. Dispersive wave turbulence in one dimension. Physica D 2001, 152/153, 551–572. [Google Scholar] [CrossRef]
Herrmann, R. Fractional Calculus: An Introduction for Physicists; World Scientific Publishing Company: Singapore, 2011. [Google Scholar]
Longhi, S. Fractional Schrödinger equation in optics. Opt. Lett. 2015, 40, 1117–1120. [Google Scholar] [CrossRef]
Majda, A.J.; McLaughlin, D.W.; Tabak, E.G. A one-dimensional model for dispersive wave turbulence. J. Nonlinear Sci. 1997, 7, 9–44. [Google Scholar] [CrossRef]
Zakharov, V.; Dias, F.; Pushkarev, A. One-dimensional wave turbulence. Phys. Rep. 2004, 398, 1–65. [Google Scholar] [CrossRef]
Zaslavsky, G.M.; Edelman, M.; Tarasov, V.E. Dynamics of the chain of forced oscillators with long-range interaction: From synchronization to chaos. Chaos 2007, 17, 043124. [Google Scholar] [CrossRef] [PubMed]
Jeng, M.; Xu, S.-L.-Y.; Hawkins, E.; Schwarz, J.M. On the nonlocality of the fractional Schrödinger equation. J. Math. Phys. 2010, 51, 062102. [Google Scholar] [CrossRef] [Green Version]
Kirkpatrick, K.; Zhang, Y. Fractional Schrödinger dynamics and decoherence. Physica D 2016, 332, 41–54. [Google Scholar] [CrossRef] [Green Version]
Luchko, Y. Fractional Schrödinger equation for a particle moving in a potential well. J. Math. Phys. 2013, 54, 012111. [Google Scholar] [CrossRef]
Duo, S.; Zhang, Y. Computing the ground and first excited states of the fractional Schrödinger equation in an infinite potential well. Commun. Comput. Phys. 2015, 18, 321–350. [Google Scholar] [CrossRef] [Green Version]
Klein, C.; Sparber, C.; Markowich, P. Numerical study of fractional nonlinear schrödinger equations. Proc. Math. Phys. Eng. Sci. 2014, 470, 20140364. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Jacob, N. Pseudo-Differential Operators and Markov Processes; Imperial College Press: London, UK, 2011; Volume I. [Google Scholar]
Oldham, K.B.; Spanier, J. The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order; Dover: Mineola, NY, USA, 2006; Volume I. [Google Scholar]
Samko, S.G.; Kilbas, A.A.; Marichev, O.I. Fractional Integrals and Derivatives; Gordon and Breach Science: Yverdon-les-Bains, Switzerland, 1993. [Google Scholar]
Duo, S.; Wang, H.; Zhang, Y. A comparative study on nonlocal diffusion operators related to the fractional Laplacian. Discret. Cont. Dyn. B 2019, 24, 231–256. [Google Scholar] [CrossRef] [Green Version]
Korabel, N.; Zaslavsky, G.M. Transition to chaos in discrete nonlinear Schrödinger equation with long-range interaction. Physica A 2007, 378, 223–237. [Google Scholar] [CrossRef] [Green Version]
Zhang, L.; He, Z.; Conti, C.; Wang, Z.; Hu, Y.; Lei, D.; Li, Y.; Fan, D. Modulational instability in fractional nonlinear Schrödinger equation. Commun. Nonlinear. Sci. Numer. Simulat. 2017, 48, 531–540. [Google Scholar] [CrossRef]
Dahlby, M.; Owren, B. Plane wave stability of some conservative schemes for the cubic Schrödinger equation. ESAIM Math. Model. Numer. Anal. 2009, 43, 677–687. [Google Scholar] [CrossRef] [Green Version]
Duo, S.; Zhang, Y. Mass conservative method for solving the fractional nonlinear Schrödinger equation. Comput. Math. Appl. 2016, 71, 2257–2271. [Google Scholar] [CrossRef] [Green Version]
Faou, E.; Gauckler, L.; Lubich, C. Plane wave stability of the split-step Fourier method for the nonlinear Schrödinger equation. Forum Math. Sigma 2014, 2, e5. [Google Scholar] [CrossRef] [Green Version]
Lakoba, T.I. Instability of the split-step method for a signal with nonzero central frequency. J. Opt. Soc. Am. B 2013, 30, 3260–3271. [Google Scholar] [CrossRef]
Weideman, J.A.C.; Herbst, B.M. Split-step methods for the solution of the nonlinear Schrödinger equation. SIAM J. Numer. Anal. 1986, 23, 485–507. [Google Scholar] [CrossRef]
Cano, B.; González-Pachón, A. Plane waves numerical stability of some explicit exponential methods for cubic Schrödinger equation. J. Comput. Math. 2016, 34, 385–406. [Google Scholar]
Bao, W.; Zhang, Y. Dynamics of the ground state and central vortex states in Bose–Einstein condensation. Math. Model. Methods Appl. Sci. 2005, 15, 1863–1896. [Google Scholar] [CrossRef]
Chen, M.; Zeng, S.; Lu, D.; Hu, W.; Guo, Q. Optical solitons, self-focusing, and wave collapse in a space-fractional Schrödinger equation with a Kerr-type nonlinearity. Phys. Rev. E 2018, 98, 022211. [Google Scholar] [CrossRef] [PubMed]
Hong, Y.; Sire, Y. A new class of traveling solitons for cubic fractional nonlinear Schrödinger equations. Nonlinearity 2017, 30, 1262–1286. [Google Scholar] [CrossRef] [Green Version]
Glassner, M.; Yevick, D.; Hermansson, B. High-order generalized propagation techniques. J. Opt. Soc. Am. B 1991, 8, 413–415. [Google Scholar] [CrossRef]
Yoshida, H. Construction of higher order symplectic integrators. Phys. Lett. A 1990, 150, 262–268. [Google Scholar] [CrossRef]
Lakoba, T.I. Instability analysis of the split-step Fourier method on the background of a soliton of the nonlinear Schrödinger equation. Numer. Meth. Part. Differ. Equ. 2012, 28, 641–669. [Google Scholar] [CrossRef]
Jordan, R.; Josserand, C. Self-organization in nonlinear wave turbulence. Phys. Rev. E 2000, 61, 1527–1539. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Qiu, Y.; Malomed, B.A.; Mihalache, D.; Zhu, X.; Peng, X.; He, Y. Stabilization of single- and multi-peak solitons in the fractional nonlinear Schrödinger equation with a trapping potential. Chaos Soliton. Fract. 2020, 140, 110222. [Google Scholar] [CrossRef]

Figure 1. Contour plots of

Re (Λ_{l}^{+})

to show the stability of traveling plane wave solutions with

a = 0.5

in the defocusing NLS. Note that the instability growth rate, indicated by the numbers on the right of the colorbars, is smaller for the smaller values of

γ

(in panels (c) and (d)).

Figure 1. Contour plots of

Re (Λ_{l}^{+})

to show the stability of traveling plane wave solutions with

a = 0.5

in the defocusing NLS. Note that the instability growth rate, indicated by the numbers on the right of the colorbars, is smaller for the smaller values of

γ

(in panels (c) and (d)).

Figure 2. Contour plots of

Re (Λ_{l}^{+})

to show the stability of standing plane wave solutions with

a = 0.5

in the focusing NLS. Note that when

α > 1

, one

Re (Λ_{l}^{+}) > 0

(i.e., instability is present) for

| μ_{l} | ≪ 1

, although their values are small.

Figure 2. Contour plots of

Re (Λ_{l}^{+})

to show the stability of standing plane wave solutions with

a = 0.5

in the focusing NLS. Note that when

α > 1

, one

Re (Λ_{l}^{+}) > 0

(i.e., instability is present) for

| μ_{l} | ≪ 1

, although their values are small.

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

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

Figure 4. Dynamics of unstable standing plane waves in the focusing NLS with

γ = - 1

and

L = 5 π

. The recurrence of plane wave solution is observed in the classical NLS.

Figure 4. Dynamics of unstable standing plane waves in the focusing NLS with

γ = - 1

and

L = 5 π

. The recurrence of plane wave solution is observed in the classical NLS.

Figure 5. Time evolution of

{log}_{10} | F [u (x, t)] |

for unstable standing plane waves in the focusing NLS with

γ = - 1

and

L = 5 π

. We show only the range

μ > 0

because the solution’s spectrum is symmetric about

μ = 0

. Also, the one-sided width of the computational spectral domain that we used was twice as large as shown in the figure.

Figure 5. Time evolution of

{log}_{10} | F [u (x, t)] |

for unstable standing plane waves in the focusing NLS with

γ = - 1

and

L = 5 π

. We show only the range

μ > 0

because the solution’s spectrum is symmetric about

μ = 0

. Also, the one-sided width of the computational spectral domain that we used was twice as large as shown in the figure.

Figure 6. Dynamics of unstable standing plane waves in the focusing NLS with

γ = - 0.8

and

L = 5 π

Figure 6. Dynamics of unstable standing plane waves in the focusing NLS with

γ = - 0.8

and

L = 5 π

Figure 7. The field and its Fourier spectrum for

L = 5 π

γ = - 1

α = 0.9

λ_{k} = 0

a = 0.5

, taken at time instances indicated in the legend. Only part of the resolved Fourier spectrum is shown in panel (b). 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}

Figure 7. The field and its Fourier spectrum for

L = 5 π

γ = - 1

α = 0.9

λ_{k} = 0

a = 0.5

10^{- 12}

Figure 8. Dynamics of unstable standing plane waves in the focusing NLS with

γ = - 1

and

L = 10 π

Figure 8. Dynamics of unstable standing plane waves in the focusing NLS with

γ = - 1

and

L = 10 π

Figure 9. Dynamics of the unstable traveling plane wave solution in the focusing NLS with

γ = - 1

and

L = 5 π

, where

λ_{k} = 2 π / L

in the initial state (39).

Figure 9. Dynamics of the unstable traveling plane wave solution in the focusing NLS with

γ = - 1

and

L = 5 π

, where

λ_{k} = 2 π / L

in the initial state (39).

Figure 10. Same as Figure 7, but for

λ_{k} = 5 \cdot (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.

Figure 10. Same as Figure 7, but for

λ_{k} = 5 \cdot (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.

Table 1. Collapse time

t_{col}

of traveling plane waves in the focusing NLS with

γ = - 1

and

L = 5 π

, rounded to the nearest 5 units.

Table 1. Collapse time

t_{col}

of traveling plane waves in the focusing NLS with

γ = - 1

and

L = 5 π

, rounded to the nearest 5 units.

$α$	0.9	0.8	0.7	0.6	0.5
$k = 1$	125	120	110	105	80
$k = 5$	290	255	245	205	85
$k = 20$	>500	355	285	240	215

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Duo, S.; Lakoba, T.I.; Zhang, Y. Dynamics of Plane Waves in the Fractional Nonlinear Schrödinger Equation with Long-Range Dispersion. Symmetry 2021, 13, 1394. https://doi.org/10.3390/sym13081394

AMA Style

Duo S, Lakoba TI, Zhang Y. Dynamics of Plane Waves in the Fractional Nonlinear Schrödinger Equation with Long-Range Dispersion. Symmetry. 2021; 13(8):1394. https://doi.org/10.3390/sym13081394

Chicago/Turabian Style

Duo, Siwei, Taras I. Lakoba, and Yanzhi Zhang. 2021. "Dynamics of Plane Waves in the Fractional Nonlinear Schrödinger Equation with Long-Range Dispersion" Symmetry 13, no. 8: 1394. https://doi.org/10.3390/sym13081394

APA Style

Duo, S., Lakoba, T. I., & Zhang, Y. (2021). Dynamics of Plane Waves in the Fractional Nonlinear Schrödinger Equation with Long-Range Dispersion. Symmetry, 13(8), 1394. https://doi.org/10.3390/sym13081394

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. Plane Wave Solution and Its Stability

3. Split-Step Fourier Spectral Method

3.1. Description of the SSFS Method

3.2. Stability of the SSFS Method

4. Dynamics of Plane Wave Solutions

4.1. Standing Plane Waves

4.2. Effects of a Larger Domain

4.3. Effects of Nonzero Group Velocity (7) of the Plane Wave

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Appendix A. Derivation of Theorem 3

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI