Instability of the finite-difference split-step method on the background of a soliton of the nonlinear Schrödinger equation

# Instability of the finite-difference split-step method on the background of a soliton of the nonlinear Schrödinger equation

## Abstract

We consider the implementation of the split-step method where the linear part of the nonlinear Schrödinger equation is solved using a finite-difference discretization of the spatial derivative. The von Neumann analysis predicts that this method is unconditionally stable on the background of a constant-amplitude plane wave. However, simulations show that the method can become unstable on the background of a soliton. We present an analysis explaining this instability. Both this analysis and the instability itself are substantially different from those of the Fourier split-step method, which computes the spatial derivative by spectral discretization. We also found that the modes responsible for the numerical instability are supported by the sides of the soliton, in contrast to unstable modes of linearized nonlinear wave equations, which (the modes) are supported by the soliton’s core.

Keywords:  Finite-difference methods, Numerical instability, Nonlinear evolution equations.

PACS number(s):  02.60.-x, 02.60.Cb, 02.60.Lj, 02.70.Bf, 42.65.Tg.

## 1 Introduction

The split-step method (SSM), also known as the operator-splitting method, is widely used in numerical simulations of evolutionary equations that arise in diverse areas of science: nonlinear waves, including nonlinear optics and Bose–Einstein condensation [1][5], atomic physics [6, 7], studies of advection-reaction-diffusion equations [8][11], and relativistic quantum mechanics [12]. In this paper we focus on the SSM applied to the nonlinear Schrödinger equation (NLS)

 iut−βuxx+γu|u|2=0,u(x,0)=u0(x). (1)

Although the real-valued constants and in (1) can be scaled out of the equation, we will keep them in order to distinguish the contributions of the dispersive () and nonlinear () terms. Without loss of generality we will consider ; then solitons of (1) exist for .

The idea of the SSM is that (1) can be easily solved analytically when either the dispersive or the nonlinear term is set to zero. This alows one to seek an approximate numerical solution of (1) as a sequence of steps which alternatively account for dispersion and nonlinearity:

 for n from 1 to nmax do:¯u(x)=un(x)exp(iγ|un(x)|2Δt)(nonlinear step)un+1(x)={solution of \ iut=βuxx at t=Δtwith initial condition u(x,0)=¯u(x)(% dispersive step) (2)

where the implementation of the dispersive step will be discussed below. In (2), is the time step of the numerical integration and . Scheme (2) can yield a numerical solution of (1) whose accuracy is . Higher-order schemes, yielding more accurate solutions (e.g., with accuracy , , etc.), are known [13, 14, 6], but here we will restrict our attention to the lowest-order scheme (2); see also the paragraph after Eq. (33) below.

The implementation of the dispersive step in (2) depends on the numerical method by which the spatial derivative is computed. In most applications, it is computed by the Fourier spectral method:

 un+1(x)=F−1[exp(iβk2Δt)F[¯u(x)]]. (3)

Here and are the discrete Fourier transform and its inverse, is the discrete wavenumber:

 −π/Δx≤k≤π/Δx, (4)

and is the mesh size in . However, the spatial derivative in (2) can also be computed by a finite-difference (as opposed to spectral) method [15, 16, 17, 18]. For example, using the central-difference discretization of and the Crank–Nicolson method, the dispersion step yields:

 iumn+1−¯umΔt=β2(um+1n+1−2umn+1+um−1n+1Δx2+¯um+1−2¯um+¯um−1Δx2), (5)

where , is a point in the discretized spatial domain:  , and is the the length of the domain. Recently, solving the dispersive step of (2) by a finite-difference method has found an application in the electronic post-processing of the optical signal in fiber telecommunications [19]. That post-processing must be done in real time, and hence the speed of its implementation becomes a key factor. Ref. [20] showed that finite-difference implementations (referred to there as the Finite and Infinite Impulse Response filters; see its Eqs. (13) and (22)) of the dispersive step of (2) provide considerable speed improvement over the Fourier-based solution (3); see also Sec. 2.4.2 in [3]. Finally, let us note that the version of the NLS where the second derivative is replaced by its finite-difference approximation, as on the right-hand side (r.h.s.) of (5), is of considerable interest also in its own right [21].

In what follows we assume periodic boundary conditions:

 u(−L/2,t)=u(L/2,t); (6)

the case of other types of boundary conditions is considered in Appendix A. We will refer to the SSM with spectral (3) and finite-difference (5) implementations of the dispersive step in (2) as s-SSM and fd-SSM, respectively. Our focus in this paper will on the fd-SSM.

Weideman and Herbst [15] used the von Neumann analysis to show that both versions, s- and fd-, of the SSM can become unstable when the background solution of the NLS is a plane wave:

 upw=(A/√γ)eiωpwt,A=const,ωpw=|A|2. (7)

Specifically, they linearized the SSM equations on the background of (7):

 un=upw+~un,|~un|≪|un| (8)

and sought the numerical error in the form

 ~un=~Aeλtn−ikx,~A=const. (9)

The SSM is said to be unstable when for a certain wavenumber one has:  (i)  Re in (9), but  (ii)  the corresponding Fourier mode in the original equation (1) is linearly stable. Weideman and Herbst found that the s- and fd-SSMs on the background (7) become unstable when the step size exceeds:

 Δtthr,s≈Δx2/(π|β|),for s-SSM% (???) \& (???) (10)

and

 Δtthr,fd=Δx/√2|β||A|2% only for β>0––––––––––––––––––––,for fd-SSM (???) \& (???) (11)

respectively. Note that for , the fd-SSM simulating a solution close to the plane wave (7) is unconditionally stable. Typical dependences of the instability growth rate, Re, on the wavenumber is shown in Fig. 1. Let us emphasize that the SSM is unstable for even though both its constituent steps, (2) and either (3) or (5), are numerically stable for any .

Solutions of the NLS (and of other evolution equations) that are of practical interest are considerably more complicated than a constant-amplitude solution (7). To analyze stability of a numerical method that is being used to simulate a spatially varying solution, one often employs the so-called “principle of frozen coefficients” [22] (see also, e.g., [23, 12]). According to that principle, one assumes some constant value for the solution and then linearizes the equations of the numerical method to determine the evolution of the numerical error (see (8) and (9)). However, as we show below, this principle applied to the SSM fails to predict, even qualitatively, important features of the numerical instability.

As a first step towards developing a general framework of stability analysis on the background of a spatially varying solution, we analyzed [24] the instability of the s-SSM on the background of a soliton of the NLS:

 usol(x,t)=A√2/γsech(Ax/√−β)eiωsolt≡Usol(x)eiωsolt,ωsol=A2. (12)

First, we demonstrated numerically that the instability growth rate in this case is very sensitive to the time step and the length of the spatial domain; also, its dependence on the wavenumber is quite different from that shown in Fig. 1(a). Moreover, the instability on the background of, say, two well-separated (and hence non-interacting) solitons can be completely different from that on the background of one of these solitons. To our knowledge, such features of the numerical instability had not been reported for other numerical methods. Moreover, they could not be predicted based on the principle of frozen coefficients. We then demonstrated that all those features could be explained by analyzing an equation satisfied by the numerical error of the s-SSM:

 i~vt−ωsol~v−β(~vxx+k2π~v)+γ|usol|2(2~v+~v∗)=0, (13)

where is proportional to the continuous counterpart of defined in (8), and is defined in the caption to Fig. 1. Note that (13) is similar, but not equivalent, to the NLS linearized about the soliton:

 i~ut−ωsol~u−β~uxx+γ|usol|2(2~u+~u∗)=0. (14)

The extra -term in (13) indicates that the potentially unstable numerical error of the s-SSM has a wavenumber close to .

In this paper we present an analysis of the numerical instability of the fd- (as opposed to s-) SSM on the soliton background. This analysis is based on an equation for the numerical error which, as (13), is a modified form of the linearized NLS. However, both that equation and its analysis are qualitatively different from those [24] for the s-SSM. Moreover, the modes that render the s- and fd-SSMs unstable are also qualitatively different. Namely, for the s-SSM, the numerically unstable modes contain just a few Fourier harmonics and hence are not spatially localized. On the contrary, the modes making the fd-SSM unstable are localized and are supported by the sides (i.e., tails) of the soliton. To our knowledge, such “tail-supported” localized modes, as opposed to those supported by the soliton’s core, have not been reported before.

The main part of this manuscript is organized as follows. In Sec. II we present simulation results showing the development of numerical instability of the fd-SSM applied to a soliton. In Sec. III we derive an equation (a counterpart of (13)) governing the evolution of the numerical error, and in Sec. IV obtain its localized solutions that grow exponentially in time. In Sec. V we summarize the results. In Appendix A we show how our analysis can be modified for boundary conditions other than periodic. In Appendix B we describe the numerical method used to obtain the localized solutions in Sec. IV. In Appendix C we discuss how the instability sets in.

## 2 Numerics of fd-SSM with soliton background

We numerically simulated Eq. (1) with , , and the periodic boundary conditions (6) via the fd-SSM algorith (2) & (5). The initial condition was the soliton (12) with :

 u0(x)=sech(x)+ξ(x); (15)

the noise component with zero mean and the standard deviation was added in order to reveal the unstable Fourier components sooner than if they had developed from a round-off error.

Below we report results for two values of the spatial mesh size , where is the number of grid points:   and . We verified that, for a fixed , the results are insensitive to the domain’s length (unlike they are for the s-SSM [24]) as long as is sufficiently large. Also, at least within the range of values considered, the results depend on monotonically (again, unlike for the s-SSM).

First, let us remind the reader that the analysis of [15] on a constant-amplitude background (7) for predicted that the fd-SSM should be stable for any .1 For the soliton initial condition (15) and the parameters stated above, our simulations showed that the numerical solution becomes unstable for . For future use we introduce a notation:

 C=(Δt/Δx)2. (16)

In Fig. 2(a) we show the Fourier spectrum of the numerical solution of (1), (15) obtained by the fd-SSM with (i.e., slightly above the instability threshold) at . The numerically unstable modes are seen around the end points of the spectral axis. At , these modes are still small enough so as not to cause visible damage to the soliton: see the solid curve in Fig. 2(b). However, at a later time, the soliton begins to drift: see the dashed line in Fig. 2(b), that shows the numerical solution at . Such a drift may persist over a long time: e.g., for , the soliton still keeps on moving at . However, eventually it gets overcome by noise and loses its identity.

We observed the same scenario for several different values of , , and (for ). The direction of the soliton’s drift appears to be determined by the initial noise; this direction is not affected by the placement of the initial soliton inside the spatial domain. The time when the drift’s onset becomes visible decreases, and the drift’s velocity increases, as increases.

The soliton’s drift is a nonlinear stage of the development of the numerical instability. It will be explained in Sec. V. In the linear stage, the numerically unstable modes are still small enough so that they do not visibly affect the soliton or one another. To describe this stage, we computed a numeric approximation to the instability growth rate  Re defined in (9):

 (17)

where (see (4)). The so computed values of the instability growth rate are shown in Fig. 3 along with the results of a semi-analytical calculation presented in Sec. IV.

The above numerical results motivate the following three questions:  (i) explain the observed instability threshold (see the sentence before (16));  (ii) identify the modes responsible for the numerical instability; and  (iii) calculate the instability growth rate (see Fig. 3). In Sec. IV we will give an approximate analytical answer to question (i). However, answers to questions (ii) and (iii) will be obtained only semi-analytically, i.e., via numerical solution of a certain eigenvalue problem.

## 3 Derivation of equation for numerical error for fd-SSM

Here we will derive a modified linearized NLS that is a counterpart of (13), which was obtained for the s-SSM in [24]. In view of periodic boundary conditions (6), the finite-difference implementation (5) of the dispersive step in (2) can be written as

 un+1(x)=F−1[eiP(k)F[¯u(x)]], (18)
 eiP(k)≡1+2iβrsin2(kΔx/2)1−2iβrsin2(kΔx/2)=exp[2iarctan(2βrsin2(kΔx/2))],r=ΔtΔx2, (19)

where , were defined after (3). For , the exponent in (19) equals that in (3). However, for , they differ substantially: see Fig. 4. It is this difference that leads to the instabilities of the s- and fd-SSMs being qualitatively different.

Using Eqs. (2) and (18), one can write, similarly to Eq. (3.1) in [24], a linear equation satisfied by a small numerical error of the fd-SSM:

 F[~un+1]=eiP(k)F[eiγ|ub|2Δt(~un+iγΔt(u2b~u∗n+|ub|2~un))]. (20)

Here is defined similarly to (8), with being either or , depending on the background solution. The exponential growth of can occur only if there is sufficiently strong coupling between and in (20). This coupling is the strongest when the temporal rate of change of the realtive phase between these two terms is the smallest. In [24] we showed that this rate can be small only for those where the exponent is close to a multiple of . Using (19) (see also Fig. 4), we see that this can occur only for sufficiently high where rather than . Then:

 P(k) = π−1βrsin2(kΔx/2)+O(1r3) (21) = π−1βr−(k−kmax)2Δx24βr+O⎛⎝1r3+((k−kmax)Δx)4r⎞⎠,

where . We have also used that

 r=Δt/Δx2=C/Δt≫1, (22)

given that the numerical instability was observed in Sec. II for .

In order for the truncation of the Taylor expansions used in (21) to be self-consistent, one needs to assume that

 O(1/r2)

note that . We arbitrarily declare to be in the middle of that range: , which implies that the range of wavenumbers for which the instability develops satisfies:

 |k−kmax|=O(1/√Δx)≪kmax=π/Δx. (23)

This turns out to be consistent with the result shown in Fig. 2(a).

Substituting the first three terms on the r.h.s. of (21) into (20), using (22), and introducing a new variable

 ~vn=(eiπ)n~un=(−1)n~un, (24)

one obtains:

 F[~vn+1] = exp(−iΔtCβ{1+(k−kmax)2Δx24})× (25) F[eiγ|ub|2Δt{~vn+iγΔt(u2b~v∗n+|ub|2~vn)}].

Note that (25) describes a small change of occurring over the step , because for , the r.h.s. of that equation reduces to . Therefore we can approximate the difference equation (25) by a differential equation. To that end, first recall from (23) that the wavenumbers of are on the order of ; hence we seek

 ~vn(x)=eikmaxx~wn(x). (26)

The effective wavenumber of is then , and according to (23) varies slowly over the scale . Introducing the scaled space and wavenumber by

 χ=x/ϵ,κ=(k−kmax)ϵ,ϵ=Δx/2, (27)

one rewrites (25) as:

 Fsc[~wn+1]=exp(−iΔtCβ{1+κ2})Fsc[eiγ|ub|2Δt{~wn+iγΔt(u2b~w∗n+|ub|2~wn)}], (28)

where now is the Fourier transform with respect to the scaled variables (27). In handling the term in (25), we have used the fact that on the spatial grid , one has:

 ~v∗n(xm)=e−ikmaxxm~w∗n(xm)=e−iπm~w∗n(xm)=eiπm~w∗n(xm)=eikmaxxm~w∗n(xm).

Next, the s-SSM (2), (3) can be written as

 F[un+1]=eiβk2ΔtF[eiγ|u|2Δtu]. (29)

When and , this is equivalent to the NLS (1) plus a term proportional to

 Δt[β∂2x,γ|u|2]−u+O(Δt2), (30)

where denotes a commutator (see, e.g., Sec. 2.4 in [3]). Equation (28) has the form of a linearized Eq. (29) with a different coefficient in the dispersion term and with an extra phase. Therefore, (28) must be equivalent to a modified linearized NLS, with the modification affecting only the corresponding terms:

 i~wt+(~wχχ−~w)/(Cβ)+γ(u2b~w∗+2|ub|2~w)=0, (31)

plus a term proportional to the linearized form of the commutator (30). Neglecting that latter term as small (of order ) compared to the rest of the expression and denoting  , we rewrite (31) as:

 iψt+δψ+ψχχ/(Cβ)+γU2b(ϵχ)(2ψ+ψ∗)=0, (32)

where

 δ=−ωb−1/(Cβ). (33)

Here is either or , and is either constant or , depending on whether the background solution is a plane wave (7) or a soliton (12). The modified linearized NLS (32) for the fd-SSM is the counterpart of Eq. (13) that was derived for the s-SSM.

Our subsequent analysis of the instability of the first-order accurate fd-SSM (2) & (5) will be based on Eq. (32). The instability of the second-order accurate version of this method, where the order of the nonlinear and dispersive steps is alternated in any two consecutive full time steps [13], is the same as that of the first-order version. The instability of higher-order versions (e.g., -accurate) can be studied similarly to how that was done in Ref. [24] for the s-SSM.

The boundary conditions satisfied by are still periodic:

 ψ(−L/(2ϵ),t)=ψ(L/(2ϵ),t). (34)

This follows from the fact that satisfies the periodic boundary conditions (6) and from (26), given that for and with some integer ,

 e−ikmaxL/2=e−iMπ=eiMπ=eikmaxL/2.

There are three differences between Eq. (32) and the linearized NLS (14). Most importantly, (32) has the opposite sign of the dispersion term. This is explained by the shape of the curve for the fd-SSM in Fig. 4 at high wavenumbers, where the curvature is opposite to that at . Secondly, unlike the -term in (14), the -term in (32) with can be either positive or negative, depending on the value of . Thirdly, the “potential” (when ) is a slow function of the scaled variable . That is, solutions of (32) that vary on the scale “see” the soliton as being very wide. This should also be contrasted with the situation for the s-SSM, where the modes described by Eq. (13) “see” the soliton as being very narrow [24].

Before proceeding to find unstable modes of Eq. (32) with , let us note that (32) with confirms the result of Ref. [15] regarding the instability of the fd-SSM on the plane-wave background. Namely, for , Eq. (32) with describes the evolution of a small perturbation to the plane wave in the modulationally stable case (see, e.g., Sec. 5.1 in [3]). That is, for , there is no numerical instability, in agreement with [15]. On the other hand, for , Eq. (32) describes the evolution of a small perturbation in the modulationally unstable case, and hence the plane wave of the NLS (1) can become numerically unstable. The corresponding instability growth rate found from (32) and Eq. (5.1.8) of [3] can be shown to agree with the one that can be obtained from Eq. (37) and the next two unnumbered relations in [15]. An example of this growth rate is shown in Fig. 1(b). Also, using our (32) and Eq. (5.1.8) of Ref. [3], the threshold value of can be shown to be given by (11), in agreement with [15].

## 4 Unstable modes of modified linearized NLS for fd-SSM

In this section we focus on the case where the background solution is a soliton; hence and (see (12)). Substituting into (32) and its complex conjugate the standard ansatz [25]    and using yet another rescaling:

 X=A√−βχ≡2A√−βxΔx,D=−Cβ2A2δ≡β2(1βA2+C),Λ=Cβ2A2λ,V(y)=2Cβ2sech2(y), (35)

one obtains:

 (∂2X+D−V(ϵX)(2112))→ϕ=iΛσ3→ϕ, (36)

where is a Pauli matrix, , and stands for a transpose. If is an eigenpair of (36), then so are ,  , and , where

 σ1=(0110)

is another Pauli matrix. Note also that is defined in the same way as in (9); hence Re indicates an instability. Below we will use shorthand notations and .

We begin analysis of (36) with two remarks. First, this equation is qualitatively different from an analogous equation that arises in studies of stability of both bright [25] and dark [26] NLS solitons in that the relative sign of the first and third terms of (36) is opposite of that in [25, 26]. This fact is the main reason why the unstable modes supported by (36) are qualitatively different from unstable modes of linearized NLS-type equations, as we will see below. While the latter modes are supported by the soliton’s core (see, e.g., Fig. 3 in [27]), the unstable modes of (36) are supported by the soliton’s “tails”.

Second, from (36) and (35) one can easily establish the minimum value of parameter where an instability (i.e., ) can occur. The matrix operators on both sides of (36) are Hermitian; the operator on the r.h.s. is not sign definite. Then the eigenvalues are guaranteed to be purely imaginary when the operator on the left-hand side (l.h.s.) is sign definite [28]; otherwise they may be complex. The third term on the l.h.s. of (36) is negative definite, and so is the first term in view of (34). The second term, , is negative when

 C<1/(|β|A2). (37)

Thus, (16) and (37) yield the stability condition of the fd-SSM on the background of a soliton. We will show later that an unstable mode indeed first arises when just slightly exceeds the r.h.s. of (37).

Since the potential term in (36) is a slow function of , it may seem natural to employ the Wentzel–Kramers-Brillouin (WKB) method to analyze it. Below we show that, unfortunately, the WKB method fails to yield an analytic form of unstable modes of (36). However, it still indicates where they can exist. Away from “turning points” (see below) the WKB-type solution of (36) is:

 →ϕ=(a+eθ+/ϵ+b+e−θ+/ϵ)→φ++(a−eθ−/ϵ+b−e−θ−/ϵ)→φ−, (38)

where are some constants, and

 (θ′±)2=−D+2V±√V2−Λ2,V≡V(ϵX),θ′≡dθ/d(ϵX), (39a) →φ±=1[(θ′±)2(V2−Λ2)]1/4⎛⎜⎝√Λ±√Λ2−V2−i√Λ∓√Λ2−V2⎞⎟⎠. (39b)

At a turning point, say, , the solution (38), (39) breaks down, which can occur because the denominator in (39b) vanishes. In such a case, one needs to obtain a solution of (36) in a transition region around the turning point by expanding the potential:  , and then solving the resulting approximate equation. For a single Schrödinger equation, a well-known solution of this type is given by the Airy function. This solution is used to “connect” the so far arbitrary constants in (38) on both sides of the turning point.

Now a turning point of (36) is where: either (i) or ,  or (ii)  . The former case can be shown (see, e.g., [29]) to reduce to a single Schrödinger equation case, where the solution in the transition region is given by the Airy function. However, at present, no such transitional solution is analytically available in case (ii)2 [30, 31]. Therefore, the solutions (38) canot be “connected” by an analytic formula across such a turning point, and hence one cannot find the eigenvalues in this case.

Based on the past experience with unstable linear modes of nonlinear waves, it is reasonable to assume that that unstable solutions of (36) must be localized. This allows one to predict the location of these solutions, which will then be verified numerically. Indeed, expressions (38) and (39a) show that for and (i.e., inside the soliton’s core). The solution (38) cannot exponentially grow from, say, the left side of the soliton to the right. Indeed, on the scale of Eq. (36), the soliton is very wide and hence, if the solution (38) is of order one on the left, it would have become exponentially large on the right side of the soliton. Therefore, inside the soliton, the solution (38) must be exponentially decaying. From (39a) it follows that it is also possible for a solution with to decay outside the soliton (i.e., where ) when . Thus, one can expect that a localized mode of (36) is to be lumped somewhere at the soliton’s side and decay in both directions away from that location.

In Fig. 5 we show the first (i.e., corresponding to the greatest ) such a mode for , points (hence ), , , . For these parameters, the threshold given by the r.h.s. of (37) is , and parameter in (36) is related to by:

 D=C−1. (40)

The numerical method of finding these modes is described in Appendix B, and the modes found by this method are shown in Fig. 5(a) for different values of . In Fig. 5(b) we show the same modes obtained from the numerical solution of the NLS (1) by the fd-SSM. These modes were extracted from the numerical solution by a high-pass filter, and then the highest-frequency harmonic was factored out as per (26). The agreement between Figs. 5(a) and 5(b) is seen to be good. In Fig. 5(d) we show the location of the peak of the first unstable mode, computed both from (36) and from the numerical solution of (1), versus parameter . The corresponding values of the instability growth rate were shown earlier in Fig. 3. Let us stress that for the localized modes of (36) was found to be purely real up to the computer’s round-off error (). There also exist unstable modes with complex , but such modes are not localized and have smaller growth rates than the localized modes.

As increases from the critical value given by (37), the localized unstable mode becomes narrower and also moves toward the center of the soliton. Moreover, higher-order localized modes of (36) arise. Typical profiles of the second and third modes are shown in Fig. 6, along with the parameter for which such modes first become localized within the spatial domain. In Appendix C we demonstrate that the process of “birth” of an eigenmode that eventually (i.e., with the increase of ) becomes localized, is rather complicated. In particular, it is difficult to pinpoint the exact value of parameter where such a mode appears. Therefore, the values shown in Fig. 6 are accurate only up to the second decimal place.

## 5 Conclusions

The main contribution of this work is the extension the (in)stability analysis of the fd-SSM beyond the von Neumann analysis in order to include spatially-varying background solutions, such as the soliton (12). We showed that, as previously for the s-SSM [24], this is done via a modified equation, (32), derived for the Fourier modes that approximately satisfy the resonance condition

 |β|k2Δt=π. (41)

Note that such modified equations — (13) for the s-SSM and (32) for the fd-SSM — are different. Their analyses are also qualitatively different, and so are the modes that are found to cause the instability of these two numerical methods. For the s-SSM, these modes are almost monochromatic (i.e., non-localized) waves that “pass” through the soliton very quickly. It is this scattering of those waves on the soliton that was shown [24] to lead to their instability. In contrast, for the fd-SSM considered in this work, the dominant unstable modes are stationary relative to the soliton. Moreover, they are localized at the sides, as opposed to the core, of the soliton. To our knowledge, such localized modes were not reported before in studies of (in)stability of nonlinear waves. We consider this finding another significant contribution of this work.

Let us stress that, as previously for the s-SSM, the principle of frozen coefficients fails to even qualitatively describe the instability of the fd-SSM on the soliton background. Indeed, that principle assumes that the spatially varying coefficient in (32) can be approximated by a constant. However, such an approach predicts [15] that any solution of the NLS with , including the soliton, should be stable (see (11) and the end of Sec. III), contrary to the findings of this work. In fact, the unstable modes of (36) localized at the soliton’s sides owe their existence to the spatial variation of the soliton’s profile.

It was straightforward to obtain an approximate threshold, (37) (where is given by (16)), beyond which the numerical instability may occur. Our simulations showed that the instability does indeed occur just slightly above that threshold. However, finding the exact threshold, , remains an open problem, and so is the calculation of the growth rate of the localized unstable modes; see Appendix C.

Let us note that the instability of the fd-SSM was studied in a recent paper [18]. However, the focus of that paper was on a rigorous mathematical proof — by a completely different method than we used here — of the stability of the numerical scheme under a certain condition, rather than on the details of the development of the instability, as in our work. That condition — their Eq. (2.9) — is a constraint on the parameter , whereas our condition (37) is a constraint on the parameter . Clearly, as , our constraint allows for a greater than that of [18], and hence is sharper. Also, it is consistent with the numerical experiments reported in that paper.

Finally, let us show how our results can qualitatively explain the observed dynamics of the numerically unstable soliton — see the text after Eq. (16). Let be the field of the unstable modes on the soliton’s sides. At an early stage of the instabilty, , the amplitude of the soliton. Also, its characteristic wavenumbers are much greater than those of the soliton: see Fig. 2(a) and (23). Accordingly, substitute into the NLS (1) and discard all the high-wavenumber terms to obtain:

 i(usol)t−β(usol)xx+γusol|usol|2=−2γusol|uunst|2. (42)

This is the equation for a perturbed soliton with the perturbation being, in general, not symmetric about the soliton’s center (see Fig. 5(c)), because the modes on the left and right sides of the soliton do not “see” each other and hence can have different amplitudes. Such a perturbation is known (see, e.g., [3], Sec. 5.4.1) to cause the soliton to move, which is precisely the effect we observed.

## Acknowledgement

I thank Jake Williams for help with numerical simulations at an early stage of this work, and Eduard Kirr for a useful discussion. This research was supported in part by the NSF grant ECCS-0925706.

## Appendix A: Modified linearized NLS for fd-SSM with non-periodic boundary conditions

We consider the Dirichlet boundary conditions (b.c.) and without loss of generality (for the purposes of this derivation) set them equal to zero. Neumann or mixed b.c. can be treated similarly, and lead to similar results.

The equation for the dispersive step of the SSM is still given by (5). However, now instead of (6) we assume: , , where and are the end points of the spatial grid. Then (5) can be rewritten as [32]

 (I+(iβr/2)A)un+1=(I−(iβr/2)A)¯u, (43)

where: , similarly for , is an identity matrix, and is an tridiagonal matrix with on the main diagonal and on the sub- and super-diagonals.

The starting point of our derivation in Sec. III, Eq. (18), has exacly the same form for the case of the Dirichlet b.c., except that is replaced with  — an expansion over the complete set of the eigenvectors of ; similarly, is replaced by . The exponential in (19) that acts on the th eigenvector is replaced by

 eiPj=1−iβrλj/21+iβrλj/2, (44)

where is the corresponding eigenvalue [32]:

 λj=−4sin2(πj/(2M)). (45)

Equations (44), (45) and the middle expression in (19) coincide provided that we identify:

 k=jπ/(MΔx)=jπ/L. (46)

However, we are still a step away from proving that the modified linearized NLS for the Dirichlet b.c. case is the same as that equation for periodic b.c.. This is because , which is the Fourier symbol of the second derivative, is not the symbol of the second derivative under the transformations and . Under those transformation, the required symbol is given by (45). We will now use this observation to supply the last step and show that the modified linearized NLS for the case of Dirichlet b.c. is indeed the same as (31). This follows from (44)–(46) and a calculation that is similar to (21):

 eiP(k) ≈ −(1+1iβr(1−sin2((k−kmax)Δx/2))) (47) ≈ −(1+1iβr+sin2((k−kmax)Δx/2))iβr),

where we have used that and that for highly oscillatory eigenvectors of , one has . The last term on the r.h.s. of (47) is the desired symbol of the second derivative, and then the rest of the derivation is the same as that leading to (28). From it one obtains the same modified linearized NLS as (31). Our numerical simulations of the NLS using the fd-SSM with zero Dirichlet b.c. confirm this conclusion.

## Appendix B: Numerical solution of eigenproblem (36)

We work with (36) written in an equivalent form:

 σ3(∂2X+D−V(ϵX)(2112)−iΛ0σ3)→ϕ=i(Λ−Λ0)σ3→ϕ, (48)

where the reason to include a constant will be explained later. We discretize (48) using Numerov’s method. It approximates the equation    by a finite-difference scheme

 Φm+1−2Φm+Φm−1=ΔX212(Fm+1+10Fm+Fm−1) (49)

with accuracy ; here , , etc., and . Then, for the discretized solution () one obtains:

 (−1)k−1([1ΔX2Aper+Nper{DI−2V−iΛ0(−1)k−1I}]fk−NperVf3−k)=i(Λ−Λ0)Nperfk, (50)

where is defined after (43);   is as in (43) except that its th and th entries equal 1 (to account for the periodic boundary conditions);   has a similar structure as :   (see (49)), , , and the rest of its entries are zero;  and . Next, defining the combined vector and matrices:

 ^V=(2VVV2V),^σ3=(IOO−I),

where is the zero matrix, one rewrites (50) as:

 ^σ3[1ΔX2^Aper+D^Nper−^Nper^V−iΛ0^σ3^Nper]^f=i(Λ−Λ0)^Nper^f. (51)

This equation has the form of the generalized eigenvalue problem where is a positive definite matrix. This problem can be solved by Matlab’s command eigs. As its options, we specified that we were looking for 24 smallest-magnitude eigenvalues and the corresponding eigenmodes. Beyond the instability threshold there are several modes with complex . We visually inspected them and found that the most unstable mode was also the most localized and also had a real eigenvalue.

We verified that the eigenvalues did not change to five significant figures whether we used or ; so we used . Finally, the need to shift the eigenvalues by a occurred for relatively large : . There, the localized modes were no longer among those with the 24 smallest-magnitude eigenvalues (see above), and significantly increasing the number of modes sought caused eigs to fail. Therefore, we had to shift the eigenvalues by in order to make the mode of interest appear among the 24 smallest-eigenvalue ones.

## Appendix C: “Birth” of localized unstable mode

The instability growth rates plotted in Fig. 3 are monotonic functions of the the parameter . This, however, occurs only when is sufficiently beyond a critical value, , where the dominant (i.e., with the greatest-) unstable mode is created. Near , which is slightly above the threshold value given by the r.g.s. of (37), the evolution of the greatest- eigenvalue is quite irregular.

Below we present results about this evolution for the dominant unstable mode (shown in Fig. 5) for , (i.e., ) and the rest of the parameters being the same as listed in Sec. II, i.e.: , , and . Then, the control parameter is related to the free term in (36) by (40). While we have been unable to rigorously establish an analytical expression for , we will present a hypothesis as to what it may be. The main message that we intend to convey is that the “birth” of a localized eigenmode of Eq. (36) occurs via a complex sequence of bifurcations, in contrast to a single bifurcation that typically takes place when an unstable mode of a nonlinear wave is “born” (see, e.g., [33]).

We numerically observed that eigenmodes of (36) with emerge not only from the origin (i.e. when imaginary eigenvalues “collide” at and give rise to two real ones), but also from “collision” of eigenvalues with whereby two3 complex eigenvalues are “born”. However, the very first (i.e., for the smallest ) unstable mode does emerge out of the origin. We will now show that this mode is essentially non-localized and, moreover, it is not the mode that eventually becomes the dominant unstable mode, whose growth rate is plotted in Fig. 3. The reason why we still chose to discuss the former mode while being primarily interested in the latter one, will become clear as we proceed.

For a mode with , Eq. (36) can be split into two uncoupled Schrödinger equations:

 (∂2X+D−ν±V(ϵX))ϕ±=0,ϕ±=ϕ1±ϕ2, (52a) where ν−=1 and ν+=3. Note that ϕ± satisfy the periodic boundary conditions: ϕ±(−L/(2ϵ))=ϕ±(L/(2ϵ)). (52b)

In Fig. 7(a) we show an example of a nontrivial solution of (52). In view of the periodic boundary conditions, this figure is equivalent to Fig. 7(b). Recall from Sec. IV that the eigenmode is exponentially small inside the soliton. Then the solution shown in Fig. 7(b) can be thought of as being localized inside the valley bounded by the two “halves” of the potential. Using this observation, one can estimate the isolated values of for which one of the equations (52a), along with (52b), has a nontrivial solution, by the WKB method. The condition for the existence of a mode localized inside the valley of Fig. 7(b) is given by the Bohr–Sommerfeld formula:

 (∫Xleft−L/(2ϵ)+∫L/(2ϵ)Xright)√D−νV(ϵX)dX=π(n+12), (53)

where is either or , is an integer, and are the turning points (see Sec. IV), where

 D−νV(ϵXleft,right)=0. (54)

The number of full oscillation periods of the mode inside the valley equals ; for example, in Fig. 7, .

When , the potential in (54) can be approximated by an exponential:  . Then, using (35) and (54), we reduce (53) to

 √D∫L/(2ϵ)Xright√1−exp[−2ϵ(X−Xright)]dX=π2(n+12), (55)

with and

 Xright=12ϵln8νCβ2D. (56)

Neglecting the exponentially small terms of the order , one obtains from (55):

 √D(L−ln8νCβ2D−2(1−ln2))=ϵπ(n+12). (57)

Note that the WKB condition (53), and hence (57), is valid when is sufficiently large. In particular, it is not supposed to accurately predict the “birth” of the first unstable mode, where . Indeed, Eq. (57) predicts that such a mode (for ) emerges at