The black solitons of one-dimensional NLS equations

Laurent DI MENZA and Clément GALLO

Analyse Numérique et EDP, Université Paris-Sud Orsay

91405 ORSAY Cedex, FRANCE

Abstract. In this paper, we prove a criterion determining if a black soliton

solution to a one-dimensional nonlinear Schrödinger equation is linearly stable
or not. This criterion handles with the sign of the limit at 0 of the Vakhitov-
Kolokolov function. For some nonlinearities, we compute numerically the black
soliton and the Vakhitov-Kolokolov function. We then show that linearly un-
stable black solitons are also orbitally unstable. In the Gross-Pitaevskii case,
we rigorously prove the linear stability of the black soliton. Finally, we nu-
merically investigate the dynamical stability of these solutions solving both
linearized and fully nonlinear equation with a finite differences algorithm.

1 Introduction
We consider here the Nonlinear Schrödinger equation (NLS) in the one-dimensional spatial
case:

i∂t u + ∂x2 u + f (|u|2)u = 0, (t, x) ∈ R2 ,
(1)
u(0) = u0 ,

where f is a real valued function and |u0(x)|2 → ρ0 as x → ∞. In all that follows, we

will assume that f (ρ0 ) = 0 and f ′ (ρ0 ) < 0. In nonlinear optics, this equation is a model
for the evolution of the complex amplitude envelope of a electric field in optical fibers, or
for the self-guided beams in planar waveguides. We refer to [KLD] for a list of physically
relevant nonlinearities f . NLS equation also appears in the description of defectons or in
the theory of one-dimensional ferromagnetic or molecular chains (see [BGMP], [BP], [B]).
Equation (1) has two invariants which are at least formally conserved by the flow: the
charge

1
 Z
 
I(u) = |u|2 − ρ0 dx
R
and the energy
Z
 
E(u) = |∇u|2 + V (|u|2) dx,
R
R ρ0
where V (r) = r f (s)ds is the potential.
In this paper, we focus on stationary solutions u(t, x) = ϕ(x) to (1). Thus we are faced
with the second-order differential equation

ϕ′′ + f (|ϕ|2)ϕ = 0. (2)

In [dB] (see also [BGMP], [BP]), A. de Bouard studies one special kind of these solutions:
the stationary bubbles, which may exist under some conditions on the nonlinearity f .
Namely, if V vanishes in the interval (0, ρ0 ) and if the largest zero of V in this interval is
of multiplicity one, there does exist a stationary bubble solving (1), that is a real, strictly
1/2
positive, even function, growing on R+ , tending to ρ0 at infinity.
Now, if V is positive on [0, ρ0 ), (1) admits an other kind of stationary solution: a black
1/2
soliton, that is a real, odd function tending to ±ρ0 at ±∞. As an example, for the Gross-
Pitaevskii equation, which is (1) with f (r) = √ 1 − r, ρ0 = 1, V (r) = (1 − r)2 /2, the black
soliton is explicitly known as ϕ(x) = tanh(x/ 2). The term “black solitons” comes from
optics, where they can be experimentally detected as black spots on bright backgrounds
(see [KLD], [KK], [MS]).
In [dB], the stationary bubbles solving (1) are shown to be all unstable. The proof is
based on a subtle analysis of the spectrum of the operator obtained by linearization of (1)
near the bubble. On the contrary, it was first thought by physicists that the black solitons
were all stable (see [KLD], [MS]). However, Y.S. Kivshar and W. Krolikowski were the first
who observed numerically unstable black solitons ([KK]). In the same paper, they justify
by variational arguments a criterion determining if a dark soliton (and in particular, a
black soliton) is stable or not. We recall that a dark soliton solving (1) is a solution of
the form u(t, x) = uv (x − vt), where v ∈ R is the speed of the soliton, and |uv (x)|2 → ρ0
as x → ∞ (see [G1]). This criterion reduces to the analysis of the sign of the derivative
dPv /dv, where Pv is the momentum of the dark soliton of speed v, given by
Z +∞  
ρ0
Pv = Im ∂x uv uv 1 − dx.
−∞ |uv |2

Namely, uv is stable if dPv /dv > 0 and unstable if dPv /dv < 0. This criterion has been
rigorously proved by Z. Lin ([L]). However, the proof is based on the hydrodynamical
form of Equation (1). A solution to (1) is written as u = (ρ0 − r)1/2 eiθ , and Z. Lin makes
the analysis on the system satisfied by (r, θx ), which has a Hamiltonian structure. This
analysis is valid for non-zero speeds (or for stationary bubbles), because the dark soliton is

2
in that case a bubble, which in particular means that |uv (x)| > 0 for x ∈ R. This approach
breaks down for the black solitons, which vanish at one point.
The goal of this paper is to fill as far as we can this gap concerning the study of the
stability of the black solitons. We make a spectral analysis as it was done by A. de Bouard
on the stationary bubbles. Denoting by ϕ a black soliton solving (1), we look for solutions
to (1) of the form ϕ + u1 + iu2 , where u1 and u2 are real valued. The linearization of (1)
near (u1 , u2) = (0, 0) writes
   
∂ u1 u1
=A , (3)
∂t u2 u2

where  
0 L1
A= ,
−L2 0
L1 = −∂x2 − f (ϕ2 ) = −∂x2 + q1 ,
L2 = −∂x2 − f (ϕ2 ) − 2ϕ2 f ′ (ϕ2 ) = −∂x2 + q1 + q2 .
We will see that f (ϕ2 ) and ϕ2 f ′ (ϕ2 ) converge exponentially at infinity. It follows as in
[dB] that the essential spectrum σe (A) of A is contained in iR.
Under this framework, we establish the following criterion of linear stability of the black
soliton1 .

Theorem 1.1 We assume that f ∈ C 3 (R+ )2 . L1 is a self-adjoint operator on L2 with

domain H 2 , its essential spectrum is [0, +∞), and L1 has a unique negative eigenvalue λ0
which is simple.
The Vakhitov-Kolokolov function

g(λ) := < (L1 − λ)−1 ϕ′ , ϕ′ >, (4)

defined for λ ∈ (λ0 , 0), is of class C 1 , increases on (λ0 , 0), and g(λ) → −∞ as λ ↓ λ0 . If
the limit of g(λ) as λ ↑ 0 is strictly positive, then A has a positive eigenvalue, and ϕ is
linearly unstable, while if it is non-positive, the spectrum of A lies in iR, and ϕ is linearly
stable.

As for the stationary bubbles studied in [dB], we prove that the black solitons which are
linearly unstable are also orbitally unstable.

Theorem 1.2 Let f ∈ C 3 (R+ ) be such that there exists a black soliton ϕ solving (1). We
assume that ϕ is linearly unstable. Then it is also orbitally unstable, in the sense that there
exists ε > 0 such that for every δ > 0, there exists v0 ∈ H 1 with kv0 kH 1 6 δ such that,
1
The notion of linear stability that we consider has to be understood as spectral stability: ϕ is said to
be linearly stable if the spectrum σ(A) is a subset of {λ, Reλ 6 0}, and linearly unstable in the other case.
2
Only f ∈ C 1 (R+ ) is needed for the sufficient condition of linear stability.

3
if u(t) ∈ ϕ + C((−T∗ , T ∗ ), H 1 ) denotes the solution to (1)3 with initial data ϕ + v0 , there
exists t ∈ (0, T ∗ ) such that inf s∈R ku(t) − τs ϕkL∞ > ε, where τs f (x) = f (x − s).

Note that we used in Theorem 1.2 the notion of orbital instability associated to the
L∞ distance, and not H 1 distance as it was done in [dB]. The reason is that we observe
numerically that the black solitons which are linearly stable, are also L∞ -orbitally stable,
but H 1 -orbitally unstable. More precisely, we observe that the solution of (1) we obtain
by perturbation of a linearly stable black soliton differs from the black soliton from a small
constant, on a spreading domain. A slightly modified version of our proof shows that such
a kind of instability also holds for the stationary bubbles, at least in dimension one.
It has also to be mentioned that since ϕ does not vanish at infinity, the instability result
above implies that

ε
∃ε > 0, ∀δ > 0, ∀u0 ∈ ϕ + H 1 with ku0 − ϕkH 1 < δ, ∃ t, inf inf ku(t) − eiθ τs ϕkL∞ > .
θ∈R s∈R 2
That is the reason why we call in the sequel orbital instability what seems to be at first
glance only instability modulo translations. Indeed, if inf s∈R ku(t) − τs ϕkL∞ > ε and θ ∈ R,
either |eiθ − 1|kϕkL∞ 6 ε/2, and thus for every s ∈ R,

ku(t) − eiθ τs ϕkL∞ > ku(t) − τs ϕkL∞ − k(1 − eiθ )τs ϕkL∞ > ε − ε/2 = ε/2,
or |eiθ − 1|kϕkL∞ > ε/2, which directly implies
√
ku(t) − eiθ τs ϕkL∞ > lim |u(t, x) − eiθ ϕ(x − s)| = ρ0 |1 − eiθ | > ε/2,
x→∞
√
since u(t) ∈ ϕ + H 1 and ϕ(x) → ρ0 = kϕkL∞ as x tends to infinity.
As it will be seen, Theorem 1.1 is a numerically efficient criterion in order to check
linear stability or linear instability. On the contrary, given a nonlinearity f such that there
exists a black soliton solving (1), it is not easy to prove rigorously via this criterion the
linear stability or the linear instability of the black soliton. We can however do it for the
Gross-Pitaevskii equation, thanks to the explicit expression of the black soliton ϕ and the
spectral properties of L1 .
√
Theorem 1.3 The black soliton ϕ(x) = tanh(x/ 2) of the Gross-Pitaevskii equation is
linearly stable.

This paper is organized as follows. In section 2, we give a definition of the black solitons
and a necessary and sufficient condition on the nonlinearity f ensuring that black solitons
exist. Section 3 is devoted to the proof of Theorem 1.1. Theorem 1.2 is proved in section
3
It is not difficult to see that the Cauchy problem (1) is locally well-posed in ϕ + H 1 , since i∆ generates
a group on H 1 and v 7→ f (|ϕ + v|2 )(ϕ + v) is locally Lipschitz continuous on H 1 . For a more detailed proof
as well as conditions under which the problem is globally well-posed, see [G2]. (−T∗ , T ∗ ) denotes here the
maximal interval of existence of the solution in ϕ + H 1 .

4
4. In section 5, we show that Theorem 1.1 is numerically efficient: we have checked for
several nonlinearities if the condition in Theorem 1.1 is satisfied or not. Theorem 1.3 is
proved in section 6. Finally, we numerically study in section 7 the dynamical linear and
nonlinear stability of the black solitons, for the pure power defocusing equation as well as
for the saturated equation. Each time, a finite differences scheme both in time and space
is used. We prove in the appendix some technical lemmas which are used in section 4.

Notations. Throughout this paper, C denotes a harmless positive constant which can
change from line to line.
If A is an unbounded operator on a Banach space, we denote by σ(A) its spectrum, by
σe (A) its essential spectrum and by ρ(A) its resolvent set.
For s ∈ R, p > 1, we define the spaces Hs := H s × H s and Lp := Lp × Lp , endowed with
their natural norms. Finally, < ·, · > denotes the L2 scalar product.

2 Existence of a black soliton

Definition 2.1 Let ρ0 > 0. A black soliton of (1) is a solution of (1) of the form

u(t, x) = ϕv (x − vt) = av eiθv (x − vt)

where v ∈ R, ϕv ∈ C 1 (R), av (0) = 0, av (x) > 0 for x ∈ R⋆ , θv is of class C 1 on R⋆ , and

1/2
piecewise C 1 on R, with θv (0+ ) = π + θv (0− ), and av → ρ0 , ∂x av → 0, ∂x θv → 0 as
x → ±∞.

Even if this definition does not a priori exclude the case where black solitons are travel-
ling at a non-zero speed v, we will see thereafter that black solitons may only be stationary
solutions (i.e. v = 0).
Next, we give a necessary and sufficient condition on the nonlinearity f which ensures
the existence of a black soliton to (1).

Theorem 2.1 Under the assumption f ∈ C(R+ ) and ρ0 > 0, there exists a black soliton ϕ
1/2
of (1), of speed v ∈ R, with amplitude tending to ρ0 at infinity if and only if the following
conditions are satisfied:
(i) f (ρ0 ) = 0
(ii) v = 0 R
ρ
(iii) V (r) := r 0 f (s)ds > 0 for r ∈ [0, ρ0 ).
Moreover, if these conditions are satisfied, up to the multiplication by some eiθ , ϕ is unique,
odd, and may be assumed to be real, positive for x > 0.
Under the further assumption f ∈ C 1 (R+ )4 and f ′ (ρ0 ) < 0 (f ′ (ρ0 ) 6 0 is a consequence of
(i) and (iii)),
√
ρ0 − ϕ(x) = e−cx(1+o(1)) ,
4
It is sufficient to assume only that f is differentiable at ρ0 .

5
 √
ϕ′ (x) ∼ c( ρ0 − ϕ(x)),
√
ϕ′′ (x) ∼ −c2 ( ρ0 − ϕ(x))
p
as x → +∞, where c := −2ρ0 f ′ (ρ0 ) is the sound speed.

Proof. Let us assume the existence of a black soliton ϕv = av eiθv = aeiθ of (1), with
1/2
amplitude tending to ρ0 at infinity.
ϕ′v = (a′ + iaθ′ )eiθ tends to 0 as x → ±∞, and ϕv solves

−ivϕ′v + ϕ′′v + f (|ϕv |2 )ϕv = 0, (5)

1/2
thus |ϕ′′v (x)| → |f (ρ0 )|ρ0 =: l0 as x → ±∞, and therefore f (ρ0 ) = 0. Indeed, let
us assume by contradiction that l0 > 0. There exists A > 0 such that x > A implies
3l0 /4 6 |ϕ′′v (x)| 6 5l0 /4. ϕv solves (5), thus, since ϕ′v and ϕ′′v are bounded, ϕv and ϕ′v are
uniformly continuous, and so is ϕ′′v . It follows that there exists η > 0 such that |x1 −x2 | 6 η
implies |ϕ′′v (x1 ) − ϕ′′v (x2 )| 6 l0 /4. If x1 , x2 > A and |x1 − x2 | = η,

|ϕ′v (x1 )| + |ϕ′v (x2 )| > |ϕ′v (x1 ) − ϕ′v (x2 )|

Z x2 Z x2
′′
= ϕv (x)dx > (|ϕ′′v (x1 )| − l0 /4)dx > l0 η/2. (6)
x1 x1

Let ε ∈ (0, l0 η/4) and B > A such that x > B implies |ϕ′v (x)| 6 ε. For x1 , x2 > B and
|x1 − x2 | = η, (6) yields a contradiction.
For x ∈ R⋆ , ϕ′v (x) = (a′ (x) + ia(x)θ′ (x))eiθ(x) . Letting x tend to 0± in this equality
and using the regularity properties of ϕ and θ, it follows that ϕ′ (0± ) = a′ (0± )eiθ(0± ) . Since
θ(0+ ) = π + θ(0− ), we get a′ (0− ) = −a′ (0+ ).
Moreover, we deduce from (5) that ϕv ∈ C 2 (R), and a straightforward computation yields
for x ∈ R⋆ ,
 
1 ρ0
θx (x) = v 1 − (7)
2 a(x)2
−axx (x) = fv (a(x)2 )a(x), (8)
2 ρ2 Rρ v2 (r−ρ0 )2
where fv (r) := f (r) + v4 (1 − r20 ). We also define Vv (r) := r 0 fv (s)ds = V (r) − 4r
.
Multiplying (8) by ax and integrating between x1 and x2 (x1 , x2 ∈ R⋆± ) leads to

−ax (x1 )2 + ax (x2 )2 = −Vv (a(x1 )2 ) + Vv (a(x2 )2 ).

Letting x2 tend to ±∞ (depending on the sign of x1 ), we get

−a′ (x1 )2 = −Vv (a(x1 )2 ), x1 ∈ R⋆ .

Therefore Vv (r) > 0 for r ∈ [0, ρ0 ). In particular, v = 0.

Since a > 0 and a(0) = 0, we have a′ (0+ ) > 0. Now, if a′ (0+ ) = 0, a ≡ 0 and ϕ cannot

6
be a black soliton. Thus, in a neighborhood of 0 in R+ , a is the solution of the Cauchy
Problem
 ′
a (x) = V (a(x)2 )1/2 ,
(9)
a(0) = 0.

Therefore a is strictly increasing as long as V (a2 ) > 0. Let us assume by contradiction that
there exists r ∈ (0, ρ0 ) such that V (r) = 0. Let r0 := min{r ∈ (0, ρ0 ), V (r) = 0} ∈ (0, ρ0 ).
Since V > 0 on (0, ρ0 ), we must have f (r0 ) = −V ′ (r0 ) = 0. The intermediate value theorem
1/2 1/2
ensures that there exists x0 ∈ R⋆+ such that a(x0 ) = r0 , and then a ≡ r0 , which is a
contradiction to the fact that ϕ is a black soliton.
On the other side, if the assumptions of the statement are satisfied, we easily construct
1/2
a black soliton ϕ = aeiθ with amplitude tending to ρ0 at infinity, in the following way :
we show that the solution a to (9) is global. Indeed, let us denote by [0, X ∗ ) the maximal
1/2
interval of existence of the solution to (9). If there exists x0 > 0 such that a(x0 ) = ρ0 ,
1/2
a′ (x0 ) = 0 and the uniqueness in the Cauchy-Lipschitz Theorem yields a ≡ ρ0 , which
1/2
is a contradiction with a(0) = 0. Thus 0 6 a(x) < ρ0 for x ∈ [0, X ∗ ), X ∗ = +∞, and
1/2
a(x) ↑ ρ0 as x → +∞. Next, we extend a to R in such a way that a is even, and we
choose θ ≡ θ+ ∈ R on R⋆+ , θ ≡ θ− = π + θ+ on R⋆− .
From now on, we assume f ∈ C 1 (R+ ) and f ′ (ρ0 ) < 0. It follows from (9) that
√
a′ (x) ∼ c( ρ0 − a(x)), (10)
x→∞

which may be integrated as

√
ρ0 − a(x) = e−cx(1+o(1)) . (11)
x→∞

The last asymptotics follows from the identity a′′ (x) = −f (a(x)2 )a(x). 

3 A linear stability criterion

In this section, we prove Theorem 1.1. Let us denote by ϕ the soliton obtained under
the conditions (i)-(ii)-(iii) of Theorem 2.1, with the extra assumption f ∈ C 3 (R+ ) and
f ′ (ρ0 ) < 0. In our study of the linear stability of ϕ, we use the same notations as in [dB].
The only difference is that here, ϕ denotes a black soliton, whereas it denotes a stationary
bubble in [dB]. As we mentioned it in the introduction, the linearization of (1) near the
black soliton ϕ yields to the system (3).
We first focus on the spectra of L1 and L2 which are self-adjoint operators on L2 . As
for the stationary bubbles studied in [dB], since f (ϕ(x)2 ) → 0 and −2ϕ(x)2 f ′ (ϕ(x)2 ) → c2
as x → ∞, σe (L1 ) = [0, +∞) and σe (L2 ) = [c2 , +∞). The fact that ϕ solves (1) writes
L1 ϕ = 0. The differentiation of this equality yields L2 ϕ′ = 0. ϕ admits exactly one zero,
whereas ϕ′ does not vanish. Therefore Sturm’s theory (see [BS]) ensures that 0 is the

7
lowest eigenvalue of L2 and that it is simple. It also suggests that L1 has a unique strictly
negative eigenvalue which is simple. That is what we prove in the next lemma5 .

Lemma 3.1 L1 has an unique strictly negative eigenvalue λ0 , which is simple. The asso-
ciated eigenvector u0 can be chosen positive.

Proof of Lemma 3.1 If v ∈ H 2 ,

Z
 
< L1 v, v >= |∂x v|2 − f (ϕ2 )|v|2 dx =: Q(v, v),
R
where Q is defined on H . Since σe (L1 ) = [0, +∞), the existence of v ∈ H 1 such that
1

Q(v, v) < 0 will imply the existence of a negative eigenvalue λ0 of L1 . Let χ ∈ C0∞ (R)
supported in [−2, 2], with χ(0) = 1. For δ > 0, we define the function ϕδ on R by

|ϕ(x)| if |x| > δ,
ϕδ (x) :=
ϕ(δ) if − δ 6 x 6 δ,
and for n, δ > 0,

ϕδ,n (x) := ϕδ (x)χ(x/n).

It is clear that ϕδ,n ∈ H 1 . Since ϕ′ ∈ L1 , thanks to the asymptotics given in Theorem 2.1,
so does ϕ′δ . Thus
Z 
2
Q(ϕδ,n , ϕδ,n ) = ϕ′δ (x)2 χ(x/n)2 + ϕδ (x)ϕ′δ (x)χ(x/n)χ′ (x/n)
R n

1 2 ′ 2 2 2 2
+ 2 ϕδ (x) χ (x/n) − f (ϕ(x) )ϕδ (x) χ(x/n) dx
n
Z
 ′ 
→ ϕδ (x)2 − f (ϕ(x)2 )ϕδ (x)2 dx.
n→∞ R

In spite of the fact that ϕδ 6∈ H 1 , we will denote this very last quantity by Q(ϕδ , ϕδ ).
Z Z
′ 2 2 2
Q(ϕδ , ϕδ ) = (ϕ (x) − f (ϕ(x) )ϕ(x) )dx + −f (ϕ(x)2 )ϕ(δ)2 dx
Z|x|>δ |x|6δ

= (ϕ′ (x)2 − f (ϕ(x)2 )ϕ(x)2 )dx

|R {z }
=0
Z δ
 ′ 2 
+ −ϕ (x) + f (ϕ(x)2 )(ϕ(x)2 − ϕ(δ)2 ) dx
−δ
′ 2
= −2δϕ (0) + o(δ),
δ→0
5
Note that in the case of the stationary bubbles studied in [dB], ϕ was positive while ϕ′ vanished exactly
once on R. Thus L1 was positive, whereas L2 had an unique negative eigenvalue. This explains why in
the sequel, the roles of L1 and L2 are exchanged in comparison with [dB].

8
where we have used the fact that L1 ϕ = 0. Thus we can choose δ > 0 small enough
such that Q(ϕδ , ϕδ ) < 0, and then n large enough in such a way that Q(ϕδ,n , ϕδ,n ) < 0.
Therefore L1 has a smallest eigenvalue λ0 < 0. It is a classical result (see for instance
[BS]) that λ0 is simple, and that there exists a positive associated eigenvector u0. Let
us assume by contradiction that L1 has another eigenvalue λ1 ∈ (λ0 , 0), and let u1 be an
associated eigenvector. Since u1 ⊥u0 and u0 > 0, u1 must vanish at some x0 ∈ R and the
Sturm’s theory implies that ϕ vanishes in both intervals (−∞, x0 ) and (x0 , +∞), which is
a contradiction, because ϕ vanishes only at 0. 

As we already mentioned in the introduction, the essential spectrum of A is purely

imaginary, as it may be shown in a similar way that for the stationary bubbles (see [dB]).
Moreover, if λ is an eigenvalue of A with associated eigenvector (u1, u2 ), −λ is also an
eigenvalue, with associated eigenvector (−u1 , u2 ). Thus the linear instability of ϕ is equiv-
alent to the existence of an eigenvalue of A with non-zero real part. If λ (6= 0) is such an
eigenvalue and if (u1 , u2 ) is an associated eigenvector, we have

L1 u2 = λu1 ,
−L2 u1 = λu2 .
−1
In particular, since L2 is self-adjoint, u2 ∈ (ker L2 )⊥ = (ϕ′ )⊥ and u1 ∈ −λL˜2 u2 +Span(ϕ′ ),
where we have denoted by L̃2 the restriction of L2 to (ϕ′ )⊥ (it is clear that L̃2 is invertible).
Therefore, taking the scalar product with u2 (which is orthogonal to ϕ′ ), and using the
fact that L̃2 is a self-adjoint coercive operator on (ϕ′ )⊥ , we obtain

< L1 u2 , u2 >
−λ2 = −1 .
< L̃2 u2 , u2 >

We deduce that −λ2 ∈ R, thus λ ∈ R ∪ iR. Therefore, if λ is an eigenvalue of A with

non-zero real part, λ ∈ R∗ , and the quantity

< L1 v, v > Q(v, v)

inf −1 = inf −1
v ∈ (ϕ′ )⊥ ∩ H 2 < L˜2 v, v > v ∈ (ϕ′ )⊥ ∩ H 1 < L̃2 v, v >
v 6= 0 v 6= 0
−1
is strictly negative. Since L˜2 is positive, we also have

Q(v, v)
α := inf 2
< 0.
v ∈ (ϕ′ )⊥ ∩ H 1 kvkL2
v 6= 0

We next prove that if α < 0, the infimum α is reached. Indeed, let (vn )n be a sequence
of (ϕ′ )⊥ ∩ H 1 such that kvn kL2 = 1 and Q(vn , vn ) → α.
Z Z
|∂x vn | dx = Q(vn , vn ) + f (ϕ2 )|vn |2 dx 6 supQ(vn , vn ) + kf kL∞ (0,ρ0 ) ,
2
R R n

9
thus (vn )n is bounded in H 1 . We can extract a subsequence (also denoted by vn for
convenience), such that vn ⇀ v∗ in H 1 . We have < v∗ , ϕ′ >= lim < vn , ϕ′ >= 0, and the
R n R
fast decrease of f (ϕR2 ) at infinity ensuresRby compactness that f (ϕ2 )|vn |2 → f (ϕ2 )|v∗ |2 .
On the other side, |∂x v∗ |2 dx 6 lim inf |∂x vn |2 dx, thus
Z Z Z Z
Q(v∗ , v∗ ) = |∂x v∗ | dx − f (ϕ )|v∗ | dx 6 lim inf |∂x vn | dx − lim f (ϕ2 )|vn |2 dx = α.
2 2 2 2
R R R R

Moreover, kv∗ kL2 6 lim inf kvn kL2 = 1. Assume by contradiction that kv∗ kL2 < 1. Then
Q(v∗ , v∗ ) α
2
6 < α < 0,
kv∗ kL2 kv∗ k2L2
Q(v∗ ,v∗ )
which is a contradiction with the definition of α. Therefore kv∗ kL2 = 1, kv∗ k2 2
= α, and
L
−1
the infimum α is reached at v∗ . Hence v∗ satisfies in H the Euler-Lagrange equation

L1 v∗ = αv∗ + βϕ′, (12)

where β ∈ R. We deduce that v∗ ∈ H 2 = D(L1 ). It can be easily seen that α > λ0 .

Indeed, it is clear that α > λ0 , and if α = λ0 , (12) implies

0 =< v∗ , (L1 − λ0 )u0 >=< (L1 − λ0 )v∗ , u0 >= β < ϕ′ , u0 >,

where u0 denotes a positive eigenvector of L1 associated with the eigenvalue λ0 . Since ϕ′

and u0 are both positive, it follows that β = 0. Eq. (12) with β = 0 implies that v∗ is
colinear to u0 , because λ0 is the only negative eigenvalue of L1 and is simple. This is a
contradiction, because u0 6∈ (ϕ′ )⊥ .
Since σ(L1 ) ∩ (λ0 , 0) = ∅, we can define on the interval (λ0 , 0) the Vakhitov-Kolokolov
fonction g by

g(λ) := < (L1 − λ)−1 ϕ′ , ϕ′ > . (13)

We have just seen that if ϕ is linearly unstable, the infimum α ∈ (λ0 , 0) of the ratio
Q(v, v)/kvk2L2 when v describes (ϕ′ )⊥ ∩H 1 is reached at some v∗ ∈ (ϕ′ )⊥ which satisfies the
Euler-Lagrange equation (12). This last equation can be rewritten as v∗ = β(L1 − α)−1 ϕ′ ,
and thus g(α) = 0.
Now, g is of class C 1 on (λ0 , 0) and g ′(λ) = k(L1 − λ)−1 ϕ′ k2L2 > 0, thus g is increasing
on (λ0 , 0). Moreover, if we denote by Π0 the orthogonal projection on u⊥ 0,

1 < u0 , . >
(L1 − λ)−1 = 2
u0 + Π0 (L1 − λ)−1 Π0 ,
λ0 − λ ku0 kL2
and thus
1 < u 0 , ϕ ′ >2
g(λ) = + < (L1 − λ)−1 Π0 ϕ′ , Π0 ϕ′ > . (14)
λ0 − λ ku0 k2L2

10
 Letting λ decrease to λ0 , the first term of the right-hand side tends to −∞, whereas the
second one remains bounded. It follows that g(λ) → −∞ as λ ↓ λ0 . Therefore g vanishes
on (λ0 , 0) if and only if the limit of g(λ) as λ ↑ 0 is strictly positive.
At this stage, we have proved that if the limit of g(λ) as λ ↑ 0 is equal to or less than
0, then the black soliton is linearly stable.
Conversely, let us assume that limλ↑0 g(λ) > 0. Let α ∈ (λ0 , 0) be such that g(α) = 0.
Denoting v := (L1 − α)−1 ϕ′ , v ∈ H 5 because ϕ′ ∈ H 4 and f (ϕ2 ) ∈ Cb3 , and we have
0 = g(α) =< v, ϕ′ >=< v, (L1 − α)v > .
Thus < L1 v, v >= αkvk2L2 < 0.
Since L̃2 is a strictly positive self adjoint operator on (ϕ′ )⊥ , one can define the contin-
−1/2 −1/2
uous operator L˜2 : L2 ∩ (ϕ′ )⊥ 7→ H 1 ∩ (ϕ′ )⊥ . Moreover, L̃2 is continuous from
−1/2
H ∩(ϕ ) into H , because f ∈ C (R+ ). Since v ∈ (ϕ ) ∩H , w0 := L˜2
3 ′ ⊥ 4 3 ′ ⊥ 3
v ∈ (ϕ′ )⊥ ∩H 4 .
We have
1/2 1/2
< L˜2 L̃1 L̃2 w0 , w0 >=< L1 v, v >< 0,
where L̃1 := ΠL1 Π is the restriction of L1 to (ϕ′ )⊥ and Π denotes the orthogonal projection
on (ϕ′ )⊥ . In particular,
1/2 1/2
< L˜2 L˜1 L̃2 w, w >
γ := inf < 0.
w ∈ (ϕ′ )⊥ ∩ H 4 kwk2L2
w 6= 0
Next, we use the two following lemmas.
Lemma 3.2 γ > −∞.
1/2 1/2
Lemma 3.3 The essential spectrum of the self adjoint operator on (ϕ′ )⊥ , L̃2 L˜1 L˜2 ,
with domain H 4 ∩ (ϕ′ )⊥ , is included in R+ .
Let us temporarily admit these lemmas. Let (wn )n∈N be a sequence in (ϕ′ )⊥ ∩ H 4 such
1/2 1/2
that for every n, kwn kL2 = 1 and < L˜2 L˜1 L˜2 wn , wn >→ γ. Thus
1/2 1/2
< (L˜2 L˜1 L̃2 − γ)wn , wn >→ 0
1/2 1/2
as n goes to infinity, which implies that 0 belongs to the spectrum of L˜2 L˜1 L̃2 − γId .
1/2 1/2
Next, thanks to Lemma 3.3, σe (L˜2 L˜1 L˜2 − γId ) ⊂ [−γ, ∞). Therefore 0 < −γ is not
1/2 1/2
in the essential spectrum of L˜2 L˜1 L̃2 − γ, which means that it is an eigenvalue. Let us
1/2
choose an associated eigenvector z0 ∈ H 4 ∩ (ϕ′ )⊥ , and define u2 := L̃2 z0 ∈ H 3 ∩ (ϕ′ )⊥ as
√
well as u1 := L1 u2 / −γ. We may easily see that (u1 , u2) is an eigenvector of A associated
√
to the eigenvalue −γ > 0, which implies that ϕ is linearly unstable. 

In order to complete the proof of Theorem 1.1, it remains to prove Lemma 3.2 and
Lemma 3.3.

11
Proof of Lemma 3.2. A similar lemma is proved in [dB]. For every w ∈ (ϕ′ )⊥ ∩ H 4 ,
we have
1/2 1/2 1/2 1/2
< L˜2 L˜1 L̃2 w, w > = < L̃1 L̃2 w, L̃2 w >
1/2 1/2
= < (L2 − q2 )L˜2 w, L̃2 w >
> kL2 wk2 2 − kq2 kL∞ < L˜2 w, w >
L
kq2 k2L∞
> kL2 wk2L2 − kq2 kL∞ kL2 wkL2 kwkL2 > − kwk2L2 .
4
We have used a Young inequality for the last estimate. It follows that γ > −kq2 k2L∞ /4 >
−∞. 

1/2 1/2
Proof of Lemma 3.3. It is quite clear that L˜2 L̃1 L̃2 , with domain H 4 ∩(ϕ′ )⊥ defines
a self-adjoint operator on L2 ∩ (ϕ′ )⊥ . Next, we write
1/2 1/2 1/2 1/2
L̃2 L̃1 L˜2 = L˜2 Π(L2 − q2 )ΠL˜2
1/2 1/2 1/2 1/2
= L˜2 (L̃2 − c2 )L̃2 + L̃2 Π(c2 − q2 )ΠL˜2
2 1/2 1/2
= L˜2 − c2 L˜2 + L˜2 Π(c2 − q2 )ΠL˜2 ,
1/2 1/2 2
and we will show that L˜2 Π(c2 − q2 )ΠL˜2 is a relatively compact perturbation of L̃2 −
1/2 1/2
c2 L˜2 . Indeed, it suffices to show that the image by L̃2 Π(c2 − q2 )ΠL˜2 of the ball
2 2
B := {f ∈ D(L̃2 − c2 L̃2 ), kf k2L2 6 1, k(L̃2 − c2 L̃2 )f k2L2 6 1}

is a compact set in L2 . If f ∈ B, then

2 2 2
kL̃2 f k2L2 =< L˜2 f, f >6 kL̃2 f kL2 kf kL2 6 k(L̃2 − c2 L˜2 )f kL2 kf kL2 + c2 kL̃2 f kL2 kf kL2 .

Thus
kL̃2 f k2L2 6 1 + c2 kL̃2 f kL2
and there exists a constant C0 > 0 such that for every f ∈ B,

kL˜2 f kL2 6 C0 .

Thus, for every f ∈ B,

2 2
kL̃2 f kL2 6 k(L̃2 − c2 L˜2 )f kL2 + c2 kL̃2 f k 6 1 + c2 C0 =: C1 ,
2 1/2
which means that B is a bounded set in D(L̃2 ) = H 4 ∩ (ϕ′ )⊥ . Next, L˜2 is continuous
4 ′ ⊥ 3 ′ ⊥ α 2
from H ∩ (ϕ ) into H ∩ (ϕ ) . Moreover, since ∂ (q2 (x) − c ) → 0 for α 6 2 as x goes to
infinity, the multiplication by c2 − q2 is compact from H 3 ∩ (ϕ′ )⊥ into H 2 . Π is continuous
1/2
from H 2 into itself because ϕ′ ∈ H 2 . Finally, L̃2 continuously maps H 2 ∩ (ϕ′ )⊥ into

12
 1/2 1/2
L2 ∩ (ϕ′ )⊥ . Therefore L̃2 Π(c2 − q2 )ΠL̃2 (B) is compact in L2 . We then deduce from
the Weyl’s criterion that
1/2 1/2 2
σe (L̃2 L˜1 L˜2 ) = σe (L˜2 − c2 L̃2 ).
2
Next, since L̃2 − c2 L˜2 is the restriction of L22 − c2 L2 to the orthogonal (ϕ′ )⊥ of the kernel
of L2 ,
2
σe (L̃2 − c2 L˜2 ) = σe (L22 − c2 L2 ).
The operator L22 − c2 L2 is a differential operator which may be explicitly computed:

L22 − c2 L2 = T + S,

where
T = ∂x4 − c2 ∂x2 ,
S = 2(c2 − r(x))∂x2 − 2r ′(x)∂x − r ′′ (x) + r(x)(r(x) − c2 )
and r(x) = q1 (x) + q2 (x). It is clear that D(T ) = H 4 ⊂ D(S). Turning into Fourier
variables, we easily compute σe (T ) = σ(T ) = R+ . Let U = C\R+ . Since the operator
S(−1 − T )−1 is compact on L2 , Theorem A1 in [H] ensures that one of the two following
properties occurs: either U ⊂ C\σe (T + S), or U ⊂ σ(T + S). The second possibility is
excluded since T + S is self adjoint, and thus σ(T + S) ⊂ R. Therefore σe (T + S) ⊂ R+ ,
and the lemma follows. 

4 Linear instability implies orbital instability

As it is mentioned by A. de Bouard in [dB], the linear instability of the stationary bubbles
implies their orbital instability. Our aim here is to prove Theorem 1.2, which states that if
a black soliton is linearly unstable, it is also orbitally unstable. As in the previous section,
we assume that f ∈ C 3 (R+ ) and f ′ (ρ0 ) < 0.
Throughout this section, we denote as before by A the operator on L2 × L2 , defined by

∀(u1 , u2 ) ∈ D(A) = H2 ,
A(u1 , u2 ) = (L1 u2 , −L2 u1 ),
and we assume that σ(A) 6⊂ iR. The analysis developed in section 3 ensures then that σ(A)
is the union of iR and of some pairs of opposite real eigenvalues with finite multiplicity.
We denote by λm > 0 the greatest eigenvalue of A. We also define an operator A1 on H1
by

∀(u1 , u2 ) ∈ D(A1 ) = H3 ,
A1 (u1 , u2 ) = A(u1 , u2 ).
Our proof of Theorem 1.2 follows the lines of that of Theorem 4.1 in [dB] (see also
Theorem 2 in [HPW]). However, in order to obtain a result of L∞ orbital instability

13
(and not only H 1 orbital instability), more accuracy is required for the description of
the spectrum of the group generated by A1 . This will be done thanks to the Spectral
Mapping Theorem of J. Prüss [P], using a method introduced by F. Gesztesy, C.K.R.T.
Jones, Y. Latushkin and M. Stanislavova in [GJLS] for localized solitary waves of Nonlinear
Schrödinger Equations. Once we have proved that σ(eA1 t ) = eσ(A)t , orbital instability is
obtained by perturbation of the black soliton in the direction of the eigenfunction associated
to the maximal real eigenvalue of A. Since the proof may be adapted to the case of a one
dimensional stationary bubble, it answers in this case to the question asked in [dB], Remark
4.1.
First, let us remark that σ(A1 ) = σ(A). Indeed, similarly to A, A1 is a relatively
compact perturbation of a differential operator with constant coefficients, which has iR
as essential spectrum. Thus σe (A1 ) = iR. Next, if (u1 , u2 ) ∈ H2 is an eigenvector of A
associated to the eigenvalue λ, ∂x2 u1 , ∂x2 u2 ∈ H 1 and thus (u1 , u2 ) ∈ H3 = D(A1 ). Therefore
λ ∈ σ(A1 ).
Next, if (u1 , u2 ) ∈ D(A1 ), we have

Re hA1 (u1 , u2 ), (u1, u2 )iH1

 
= Re −∂x2 u2 + q1 u2 , u1 H 1 − −∂x2 u1 + q1 u1 + q2 u1 , u2 H1
= Re [hq1 u2 , u1 iH 1 − h(q1 + q2 )u1 , u2 iH 1 ] .

In particular,

|Re hA1 (u1 , u2 ), (u1, u2 )iH1 | 6 (2kq1 kH 1 + kq2 kH 1 )ku1 kH 1 ku2kH 1

6 Ck(u1 , u2)k2H1 .

It follows that A1 is the generator of a C0 -group on H1 , which will be denoted by eA1 t .

We are going to prove that the spectrum of eA1 t equals the set eσ(A)t . This will be done
using the following Spectral Mapping Theorem, which is due to J. Prüss.

Theorem 4.1 Let X be a Hilbert space, A an unbounded operator on X and S(t) a C0

semigroup on X, with generator A. Let µ 6= 0. Then µ ∈ ρ(S(t)) if and only if E := {λ ∈
C, eλt = µ} ⊂ ρ(A) and sup{k(λ − A)−1 kL(X) , λ ∈ E} < ∞.

In order to apply Theorem 4.1 to the operator A1 , we will first prove uniform estimates
on the resolvent of A for λ ∈ Γa,τ0 := {a + iτ, |τ | > τ0 }, where a ∈ R∗ and τ0 > 0 is large.
We will then easily deduce uniform estimates on the resolvent of A1 , for λ ∈ Γa,τ0 .
Let us now prove the uniform estimates on the resolvent of A. Our method is similar to
the one used by F. Gesztesy, C.K.R.T. Jones, Y. Latushkin and M. Stanislavova in [GJLS]
for localized solitary waves. For λ = a + iτ , where a, τ ∈ R∗ , we have

λ − A = Bλ + MQ ,
where

14
  
λ ∂x2
Bλ =
−∂x2 + c2 λ
and
   
0 Q1 0 f (ϕ2 )
MQ = = .
Q2 0 −f (ϕ2 ) − 2ϕ2 f ′ (ϕ2 ) − c2 0
 
0 −∂x2
Since the spectrum of is included in iR and λ 6∈ iR, Bλ is
−(−∂x2 + c2 ) 0
invertible, and λ − A may be rewritten as
 
λ − A = Bλ Id + Bλ−1 MQ .
Next, we have the following lemma, which is proved in the appendix.

Lemma 4.1 Let a 6= 0. Then there exists Ca > 0 such that for every λ = a + iτ , |τ | > 1,

kBλ−1 kL(L2 ) 6 Ca .

Let Tλ be the continuous operator on L2 defined by

Tλ := Bλ−1 MQ
 
−∂x2 (λ2 − ∂x2 (−∂x2 + c2 ))−1 Q2 λ(λ2 − ∂x2 (−∂x2 + c2 ))−1 Q1
= . (15)
λ(λ2 − ∂x2 (−∂x2 + c2 ))−1 Q2 −(−∂x2 + c2 )(λ2 − ∂x2 (−∂x2 + c2 ))−1 Q1

We will next prove that each of the four block operators on L2 in the right hand side of
(15) tends to zero in L(L2 ). This will imply that for τ0 large enough, kTλ kL(L2 ) 6 1/2 as
soon as λ ∈ Γa,τ0 . We deduce that Id + Tλ is invertible in this case and thanks to Lemma
4.1,
k(λ − A)−1 kL(L2 ) 6 2Ca .
Let us define now

Kλ1 = Q(x)(−∂x2 )(λ2 − ∂x2 (−∂x2 + c2 ))−1 ,

Kλ2 = Q(x)λ(λ2 − ∂x2 (−∂x2 + c2 ))−1

and

Kλ3 = Q(x)(−∂x2 + c2 )(λ2 − ∂x2 (−∂x2 + c2 ))−1 ,

where Q = Q1 or Q2 , in such a way that Kλj ∈ L(L2 ), j = 1, 2, 3 are the adjoints of
∗
the block operators in the right-hand side of (15). Since kKλj kL(L2 ) = kKλj kL(L2 ) , we are

15
reduced to prove that kKλj kL(L2 ) tends to zero as τ goes to infinity. As in [GJLS], this will
be done as follows.
We first remark that for j = 1, 2, 3, Kλj writes Qgλj (i∂x ), where

x2
gλ1 (x) = ,
λ2 + x2 (x2 + c2 )
λ
gλ2 (x) =
λ2 + x2 (x2 + c2 )
and

x2 + c2
gλ3 (x) = .
λ2 + x2 (x2 + c2 )
Thus Kλj is an integral operator with kernel (x, y) 7→ Q(x)gˇλj (x − y), where ǧ denotes
the inverse Fourier transform of g.
Note that Q converges exponentially to zero, therefore it clearly belongs to L2 . Thus,
since |gλ1 | 6 |gλ3 |, it suffices to prove that for j = 2, 3, gλj ∈ L2 (R) and kgλj kL2 → 0 as τ
tends to infinity, which will imply that the Kλj are Hilbert-Schmidt operators which tend
to 0 in L(L2 ) as τ goes to infinity. More generally, we prove the following lemma.

Lemma 4.2 Let λ = a + iτ where a 6= 0. For j = 2, 3, for every q > 1, gλj belongs to
Lq (R) and kgλj kLq → 0 as τ tends to infinity.

Next, we show that the uniform estimate on k(λ − A)−1 kL(L2 ) for λ ∈ Γa,τ0 that we just
obtained imply a similar uniform estimate for k(λ − A)−1 kL(H1 ) . Indeed, let a ∈ R∗ and
τ0 > 0, C > 0 such that for every λ = a + iτ ∈ Γa,τ0 ,

k(λ − A)−1 kL(L2 ) 6 C.

Then if λ ∈ Γa,τ0 , F ∈ H1 and V = (λ − A)−1 F , then

kVkL2 6 CkF kL2 . (16)

We also get that V ∈ H3 and

(λ − A)∂x V = ∂x F − MQ′ V,
 
0 Q′1
where MQ′ = . Thus
Q′2 0

k∂x VkL2 6 C(k∂x F kL2 + (kq1′ kL∞ + kq2′ kL∞ )kVkL2 )

6 CkF kH1 . (17)

16
 It follows from (16) and (17) that (λ − A1 )−1 : H1 7→ H1 is continuous, with a norm
uniformly bounded as λ stays in Γa,τ0 .

At this stage, we conclude by Theorem 4.1 that

σ(eA1 t ) = eσ(A1 )t = eσ(A)t . (18)

In particular, the spectral radius of eA1 t is eλm t . Once (18) has been established, the proof
of Theorem 1.2 follows the broad lines of Theorem 2 in [HPW] and Theorem 4.1 in [dB].

It can easily be seen (see for instance [G2]) that there exists a neighborhood U of 0 in
H such that for every (v1 , v2 ) ∈ U, there is a unique solution u ∈ ϕ + C([0, 1], H 1 + iH 1 ) to
1

(1), such that u(0) = ϕ + v1 + iv2 . For t ∈ [0, 1] and (v1 , v2 ) ∈ U, we define the nonlinear
flow

T (t)(v1 , v2 ) := (Reu(t) − ϕ, Imu(t)).

Using these notations, Equation (1) equivalently expresses as
Z t
A1 t
U(t) = T (t)(v1 , v2 ) = e (v1 , v2 ) + eA1 (t−τ ) g(T (τ )(v1 , v2 ))dτ, (19)
0

where g : H1 7→ H1 is defined by

 
g(v1 , v2 ) = f (ϕ2 ) − f ((ϕ + v1 )2 + v22 ) v2 ,


f ((ϕ + v1 )2 + v22 )(ϕ + v1 ) − f (ϕ2 )(ϕ + v1 ) − 2ϕ2 f ′ (ϕ2 )v1 .

It is quite clear that there exists a0 , c0 > 0 such that U contains the ball of radius a0 in
H and kg(v1 , v2 )kH1 6 c0 k(v1 , v2 )k2H1 as soon as k(v1 , v2 )kH1 6 a0 . If k(v1 , v2 )kH1 < a0 , let
1

t∗ = sup{t ∈ [0, 1], kU(t)kH1 < a0 } > 0. Since the spectral radius of eA1 t is eλm t , for every
η > 0, there exists M > 0 such that keA1 t kL(H1 ) 6 Me(λm +η)t . Thus for every t ∈ [0, t∗ ],
Z t
(λm +η)t
kU(t)kH1 6 Me kU0 kH1 + Ma0 e(λm +η)(t−τ ) c0 kU(τ )kH1 dτ,
0
−(λm +η)t
where U0 = U(0) = (v1 , v2 ). Multipliying by e and using the Gronwall lemma, we
get

kU(t)kH1 6 MkU0 kH1 e(λm +η+M c0 a0 )t , t ∈ [0, t∗ ]. (20)

Thus, one can choose a1 > 0 small enough (for instance a1 = e−(λm +η+M c0 a0 ) a0 /M), in
such a way that kU0 kH1 < a1 implies kU(t)kH1 6 a0 for every t ∈ [0, 1]. For such a choice
of U0 , for every t ∈ [0, 1], it follows from (19) and (20) that

17
 Z t
A1 t 3 (λm +η)t
kU(t) − e U0 kH1 6 c0 M e kU0 k2H1 e(2M c0 a0 +λm +η)τ dτ
0
6 bkU0 k2H1 , (21)

where b = c0 M 3 e(λm +η) (e(2M c0 a0 +λm +η) − 1)/(2Mc0 a0 + λm + η). From now on, we fix
η ∈ (0, λm ). Let Cs be the norm of the continuous Sobolev embedding H1 ⊂ L∞ . Let Um
be an eigenvector of A1 associated to the eigenvalue λm , normalized in such a way that
kUm kH1 = 1. We fix a2 ∈ (0, a1 ) and ∆ = Mba2 e−(λm +η) /(eλm −η − 1), where a2 is small
enough, such that ∆ < 1/2 and ∆Cs /(1 − ∆) < kUm kL∞ /2. Let R = 1/(1 − ∆) > 1, let
N ∈ N, and ε = a2 e−λm N /R. Let us finally choose U0 = εUm .
We first show by induction on n that, denoting by T ∗ (U0 ) the maximal time of existence
in H1 of the solution to (19) with initial data U0 , we have for every n ∈ {0, . . . , N},
T ∗ (U0 ) > n and

kU(n)kH1 6 εReλm n . (22)

For n = 0, kU(0)kH1 = ε < εR by choice of U(0) and R. Let us now take n ∈ {1, . . . , N},
and assume that (22) holds at orders k ∈ {0, . . . , n − 1}. Then

kU(n − 1)kH1 6 εReλm (n−1) 6 εReλm N = a2 < a1 ,

therefore the time of existence of U is at least n, and we have

kU(n)kH1 6 kU(n) − enA1 U0 kH1 + kenA1 U0 kH1

n−1
X  
6 e(n−1−k)A1 U(k + 1) − eA1 U(k) + kenA1 U0 kH1
k=0 H1
n−1
X
6 Me(λm +η)(n−1−k) bkU(k)k2H1 + eλm n kU0 kH1 ,
k=0

because of (21), and since U0 is an eigenvector of enA1 associated to the eigenvalue eλm n .
The induction assumption then ensures that

 
λm n −λm η(n−1) e(λm −η)n − 1
kU(n)kH1 6 εe Mbe e εR2 λm −η +1
e −1
6 εeλm n (R∆ + 1) = εeλm n R,

which is the announced result. Note that the computation we just made also proves that
for every n ∈ {0, . . . , N},

18
 kU(n) − enA1 U0 kH1 6 εeλm n R∆. (23)

Thus, for n ∈ {0, . . . , N},

kU(n)kL∞ > kenA1 U0 kL∞ − kU(n) − enA1 U0 kL∞

εeλm n
> εeλm n (kUm kL∞ − Cs ∆R) > kUm kL∞ . (24)
2
For n = N, we get

εeλm N a2
kU(N)kL∞ > kUm kL∞ = kUm kL∞
2 2R
a2 a2
= (1 − ∆)kUm kL∞ > kUm kL∞ . (25)
2 4

Finally, we show that U(N) stays far from the curve Γ = {v = (ϕ(· − h) − ϕ, 0), h ∈ R}
in L∞ , which will complete the proof of Theorem 1.2.

Let us first define, for ρ > 0, the cones

n o
Σρ = w ∈ L∞ , w − kwkL∞ Uf
m 6 ρkwkL∞
L∞
and
±
n o
f
Σρ = w ∈ L∞ , w ∓ kwk L∞ ϕe′ 6 ρkwk L∞ ,
L∞

where Uf e′ ′ ′
m = Um /kUm kL∞ and ϕ = (ϕ , 0)/k(ϕ , 0)kL∞ .
We claim that for n ∈ {0, . . . , N}, we have U(n) ∈ Σρ1 , where ρ1 = 4∆R/kUm kL∞ .
Indeed, since enA1 U0 = eλm n U0 = εeλm n Um and thanks to (23) and (24),

kU(n)kL∞
U(n) − kU(n)kL∞ Uf
m 6 kU(n) − enA1 U0 kL∞ + εeλm n Um − Um
L∞ kUm kL∞ L∞
kU(n)kL ∞
6 kU(n) − enA1 U0 kL∞ + εeλm n − kUm kL∞
kUm kL∞
6 kU(n) − enA1 U0 kL∞ + kenA1 U0 kL∞ − kU(n)kL∞
4∆R
6 2∆Rεenλm 6 kU(n)kL∞ .
kUm kL∞
±
g
Next, we show that if v ∈ Γ is close to 0, then v ∈ Σρ2 for some ρ2 > 0. Let v ∈ Γ
√
with kvkL∞ < ρ0 /2. There exists h ∈ R such that v = (v1 , 0) with v1 = ϕ(· + h) − ϕ. We

19
 √ √

assume for instance h > 0. Since ϕ is a diffeomorphism from − ϕ−1 ( ρ0 /2), ϕ−1( ρ0 /2)
√ √ √
onto (− ρ0 /2, ρ0 /2), and since v1 (0) = ϕ(h) ∈ (0, ρ0 /2), the Mean Value Theorem
induces √ √
ρ0 −1 ′ ρ0
h< sup (ϕ ) (s) 6 =: α.
2 |s|6 ρ0 /2
√ 2 inf
−1 √
|ϕ′ (s)|
|s|6ϕ ( ρ0 /2)

It follows from the triangle inequality and a Taylor formula that

v1 v1 ϕ′
v − kvkL∞ ϕe′ = h −
L∞ h h L∞ kϕ′ kL∞ L∞
 
v1 kv1 /hkL∞
6 h − ϕ′ ′
+ kϕ kL∞ 1−
h L∞ kϕ′ kL∞
v1
6 2h − ϕ′ 6 2h2 kϕ′′ kL∞ . (26)
h L∞

Moreover, for h 6 α, since ϕ′ is decreasing on R+ , the Taylor expansion at order one at

x = 0 yields

ϕ(h) kv1 kL∞ kvkL∞

h = R1 6 ′
= ′ . (27)
ϕ′ (sh)ds ϕ (α) ϕ (α)
0

+
g
We deduce from (26) and (27) that v ∈ Σρ2 with ρ2 = 2kvkL∞ kϕ′′ kL∞ /ϕ′ (α)2 . If h 6 0,
−
g
similar computations induce v ∈ Σρ . 2
±
g
Finally, if v ∈ Σρ2 , the triangle inequality gives

kU(N) − vkL∞ > kU(N)kL∞ Uf e′

m ∓ kU(N)kL∞ ϕ
L∞

− U(N) − kU(N)kL∞ Uf
m − (kU(N)kL∞ − kvkL∞ )ϕe′ − kvkL∞ ϕe′ ∓ v ,
L∞ L∞ L∞

thus, since U(N) ∈ Σρ1 ,

(1 + kϕe′ kL∞ )kU(N) − vkL∞ > kU(N)kL∞ Uf e′

m ∓ϕ − ρ1 kU(N)kL∞ − ρ2 kvkL∞ . (28)
L∞

Up to a change of a2 , one may assume that ρ1 = 4∆R/kUm kL∞ < kUf e′

m ∓ ϕ kL∞ /2,
√
2Cs a2 6 ρ0 /2 and 8Cs a2 kϕ′′ kL∞ /ϕ′ (α)2 < kUf e′
m ∓ ϕ kL∞ /4. Under these conditions, if
v ∈ Γ, two cases may occur: either kvkL∞ > 2kU(N)kL∞ , which implies by (25)

a2
kU(N) − vkL∞ > kU(N)kL∞ > kUm kL∞ , (29)
4

20
or kvkL∞ 6 2kU(N)kL∞ . In this last case, it follows from (22) that
√
kvkL∞ 6 2Cs a2 6 ρ0 /2.
Thus v = (ϕ(· + h) − ϕ, 0) with h 6 α, and (27) holds, as well as (28). Moreover, using
once again (22)

8kU(N)k2L∞ kϕ′′ kL∞

ρ2 kvkL∞ 6
ϕ′ (α)2
8Cs εeλm N Rkϕ′′ kL∞
6 kU(N)kL∞
ϕ′ (α)2
8Cs a2 kϕ′′ kL∞ 1 f
= ′ 2
kU(N)kL∞ 6 Um ∓ ϕe′ kU(N)kL∞ . (30)
ϕ (α) 4 L∞

Using (28), (29) and (30), we finally obtain, for every v ∈ Γ,

!
a2 a2
kU(N) − vkL∞ > min kUm kL∞ , Uf e′
m ∓ϕ kUm kL∞ ,
4 16(1 + kϕe′ kL∞ ) L∞

and the proof of Theorem 1.2 is now complete.

5 Numerical check of linear stability

In this section, we numerically investigate the linear stability of the black solitons for
different kinds of nonlinearities, using the criterion given in Theorem 1.1. First, it is
possible to find an approximate value Φ of the black soliton ϕ using a differential solver,
since ϕ(0) = 0 and the value ϕ′ (0) is explicitly known. In fact, if we multiply (2) by ϕ′
and integrate between 0 and infinity, it can be shown that
p
ϕ′ (0) = V (0).
For instance, for the pure power defocusing equation

∂u ∂ 2 u
i + 2 + (1 − |u|2σ )u = 0, (31)
∂t ∂x
which corresponds to f (r) = 1 − r σ , we have
r
′ σ
ϕ (0) = .
σ+1
It is sufficient to compute the approximate solution of (2) on a sufficiently large domain
in such a way that the derivative almost vanishes at the boundary. For some prescribed λ,
the Vakhitov-Kolokolov function is given by

21
 g(λ) =< (L1 − λId )−1 ϕ′ , ϕ′ >:=< ψ, ϕ′ >,
where ψ ∈ L2 (R) satisfies the differential equation (L1 − λId )ψ = ϕ′ . This equation can be
solved with use of finite differences, seeking the approximate values ψj of the solution at
gridpoints xj = jh where j = −J, · · · , J. In this case, the discretized linear system can be
written as MΨ = Φ′ , where Ψ = (ψ−J , ψ−J+1 , · · · , ψJ−1 , ψJ )T ∈ R2J+1 , Φ′ ∈ R2J+1 stands
for the approximate value of ϕ′ computed from Φ = (ϕ−J , ϕ−J+1, · · · , ϕJ−1 , ϕJ )T and the
matrix M is given by
 
2 −1 0 ··· 0
 . .. 
 −1 2 −1 . . .   
1  . . . 

M= 2 0 . . . . . . 0  + diag q1 (xj ) − λ, −J ≤ j ≤ J .
h  . 
.
 .. . . −1 2 −1 
0 ··· 0 −1 2
The approximation G of the Vakhitov-Kolokolov function g is then given by G(λ) =
< Ψ, Φ >. In Figures 1 and 2 are respectively plotted the profiles of black solitons obtained
for (31) considered with different values of σ (σ = 0.5, 1, 2 and 3) and for each case the
Vakhitov-Kolokolov function with respect to log (−λ).
1 −1

0.8 −1.1

0.6 −1.2

0.4 sigma = 0.5 −1.3

sigma = 1
sigma = 2
0.2 sigma = 3 −1.4

0 −1.5

−0.2 −1.6

−0.4 −1.7 sigma = 0.5

sigma = 1
sigma = 2
−0.6 −1.8 sigma = 3

−0.8 −1.9

−1 −2
−6 −4 −2 0 2 4 6 −20 −18 −16 −14 −12 −10 −8 −6
x log(−lambda)

Figure 1: Black solitons for different Figure 2: G with respect to log(−λ) for
values of σ. different values of σ.

It can be observed that the function G is always strictly negative when λ tends to
zero and admits at each time a strictly negative limit at zero. This numerically checks
the linear stability in this case. It can be also observed that this limit seems to be an
increasing function of σ. Nevertheless, when σ becomes large, it has been found that the
limit value l(σ) when λ tends to zero is always strictly negative. As indicated in Figure 3,

22
 −1.1

−1.2

−1.3

−1.4
0 5 10 15 20 25 30 35 40

Figure 3: Plot of l(σ) versus σ.

it seems that l := lim l(σ) exists and remains strictly negative (we have l ≃ −1.13). This
σ→∞
result suggests the linear stability of black solitons for any nonnegative value of σ.
We then focus on the case of the saturated equation

∂u ∂ 2 u  1 1 
i + 2+ − u = 0, (32)
∂t ∂x (1 + a|u|2 )2 (1 + a)2
where a is a prescribed positive parameter. Here, the nonlinear term still vanishes for
|u|2 = 1 and remains positive between the two values 0 and 1. The value at the origin of
the black soliton derivative is given by
√
′ a
ϕ (0) = .
1+a
It is still possible here to compute the soliton and the Vakhitov-Kolokolov function. In
Figures 4 are shown the profiles of black solitons for different values of a. It can be noticed
that contrary to the first equation investigated, the black soliton derivative at zero is no
more a monotonous function of the parameter (in fact, the maximum value obtained is
reached for a∗ = 1). In figure 5 is plotted the value G(λ0 ) for λ0 = −10−5 for different
values of a. We observe that for sufficiently large values of a (say, a > aVK ), the Vakhitov-
Kolokolov function becomes positive if λ is close to zero, which ensures that the black
soliton is linearly stable when a < aVK ≃ 7.47 and linearly unstable when a > aVK . This
kind of instability threshold has been already observed in the literature (see [KK]).

23
 1
0.8
0.8

a = 0.2
0.6 a = 0.5 0.6
a=1
a=2
0.4
0.4
0.2

0 0.2

−0.2
0
−0.4

−0.6 −0.2

−0.8
−0.4
−1
−6 −4 −2 0 2 4 6 6 6.5 7 7.5 8 8.5 9 9.5 10
x a

Figure 4: Black solitons for different a. Figure 5: G(λ0 ) versus a.

6 The linear stability of the black soliton of the 1D

Gross-Pitaevskii equation
In the cubic defocusing case f (r) = 1 − r, the black soliton ϕ is explicitly known:
x
ϕ(x) = tanh √ .
2
Thus
1
ϕ′ (x) = √ ,
2 cosh2 √x
2

1
f (ϕ(x)2 ) = ,
cosh2 √x
2

and
d2 1
L := L1 = − − .
dx 2
cosh2 √x
2

We also know that the lowest eigenvalue of L is λ0 = −1/2, and an associated eigenvector
is
1
u0(x) = .
cosh √x2
The explicit form of this data enables us to prove the linear stability of the black soliton
of the Gross-Pitaevskii equation stated in Theorem 1.3.

Proof of Theorem 1.3. Let λ ∈ (−1/2, 0), and

v := (L − λ)−1 ϕ′ ∈ H 2 ⊂ L∞ .

24
Then it is clear that v is even, v(x) → 0 as x → ∞, and v solves

v(x) 1
v ′′ (x) = − 2 √x −
√ − λv(x). (33)
cosh 2 2 cosh2 √x
2

For every integer n > 1, let us define

Z +∞
v(x)
un := dx
−∞ cosh2n √x2

and Z +∞
1
an := dx.
−∞ cosh2n √x
2
√
Note that u1 = 2g(λ). The strategy of the proof is as follows. First, multiplying (33)
by 1/ cosh2n √x2 and integrating over R, we obtain an induction relation between un and
un+1 . Summing this relation, we express the first term in the asymptotic expansion of un
in function of u1 . Next, using the definition of un , we express the same quantity in function
of v(0). The comparison of these two expressions yields a link between u1 and v(0). In
particular, this relation implies that v(0) > 0 as soon as u1 > 0. Finally, we show that the
condition v(0) > 0 is inconsistent with the fact that v(x) → 0 as x → ∞.

Let us as announced multiply (33) by 1/ cosh2n √x2 and sum the obtained equality over
R. We get
Z +∞
v ′′ (x) 1
2n √x dx = −un+1 −
√ an+1 − λun . (34)
−∞ cosh 2
2

Using integration by parts, the left hand side of (34) can be expressed as
Z +∞
v ′′ (x) 2
2n √x dx = 2n un − (2n + 1)nun+1 . (35)
−∞ cosh 2

From (34) and (35) we infer the following induction relation:

1 + 2nλ2 1 an+1
un+1 = 1 1 un +
√ 2
. (36)
1 + 2n − 2n2 2 2n + n − 1

(36) yields a relation between un+1 and u1 :

n
! n n
!
Y 1 + 2lλ2 1 X Y 1 + 2lλ2 ak+1
un+1 = u1 + √ . (37)
l=1
1 + 2l1 − 2l12 2 k=1 l=k+1 1 + 2l1 − 2l12 2k 2 + k − 1

25
 We will next compute asymptotic expansions of each term
Q∞ in (37).
We begin with the
λ
term containing u1 . First remark that the product Pλ := l=1 1 + 2l2 converges. Thus
n 
Y 
λ
1 + 2 = Pλ + o(1),
2l
l=1

with Pλ ∈ (0, 1) (because λ ∈ (−1/2, 0)). On the other side, there exists a constant P > 0
such that n  
Y 1 1 √
1 + − 2 = P n(1 + o(1)).
2l 2l
l=1

Therefore
n
!
Y 1 + 2lλ2 Pλ
= √ (1 + o(1)). (38)
l=1
1 + 2l1 − 2l12 P n

We begin the computation of the asymptotic expansions of the two others terms in (37)
with the following lemma:

Lemma 6.1 Let ϕ ∈ C ∞ (R), bounded, even, with ϕ(0) 6= 0. Then

Z +∞ Z
ϕ(x) ϕ(0) ∞ e−y
= √
2n dx n→∞ √ dy(1 + o(1)).
−∞ cosh x n 0 y

Proof. We write
Z +∞ Z ∞
ϕ(x)
dx = 2 ϕ(x)e−2n ln cosh x dx,
−∞ cosh2n x 0

and we make the change of variables y = 2n ln cosh x, such that

Z +∞ Z ∞ y
ϕ(x) ϕ(arccoshe 2n ) −y y dy
2n dx = 2
p y e e 2n
−∞ cosh x 0 e −1
n 2n
Z ∞ −y
1 y e
= √ ϕ̃( ) √ dy,
n 0 n y

where for y > 0, y√

e2 y y
ϕ̃(y) := ϕ(arccoshe ) √ y 2
e −1
and ϕ̃(0) = ϕ(0). ϕ̃ is continuous on R+ , and the dominated convergence theorem yields
Z ∞ Z ∞
y e−y e−y
ϕ̃( ) √ dy −→ ϕ(0) √ dy,
0 n y n→∞ 0 y

and the lemma follows. 

26
 √
We apply this lemma on the one side with ϕ(x) = v( 2x), on the other side with
ϕ(x) = 1. We obtain respectively
√ √
√ Z ∞ v( 2x) Z
2v(0) ∞ e−y
un = 2 2n dx =
√ √ dy(1 + o(1)) (39)
−∞ cosh x n 0 y

and
√ Z ∞ −y
√ Z ∞
dx 2 e
an = 2 2n =√ √ dy(1 + o(1)). (40)
−∞ cosh x n 0 y

We next give an asymptotic expansion of the last term in (37). We rewrite it under the
form
n n
!
1 X Y 1 + 2lλ2 ak+1
√ 1 1
2 k=1 l=k+1 1 + 2l − 2l2 2k 2 + k − 1
n
! n
1 Y 1 + 2lλ2 X ak+1
= √ 1 1
  .
2 l=1 1 + 2l − 2l2 k=1 Qk 1+ 2λ
2
l=1
2l
1+ 1 − 1
(2k + k − 1)
2l 2l2

The series in the right hand side converges and we will denotes its sum by Sλ . Indeed, it
follows from (38) and (40) that
√ √ Z ∞ −y
ak+1 2P k e 1 1
  = √ √ dy (1 + o(1)) = O( ).
Qk 1+ λ2
2 2 kPλ 0 y k 2 k2
l=1
2l
1+ 1 − 1
(2k + k − 1)
2l 2l2

Therefore, using once again (38), we get

n n
!
1 X Y 1 + 2lλ2 ak+1 P λ Sλ
√ 1 1 2
= √ (1 + o(1)). (41)
2 k=1 l=k+1 1 + 2l − 2l2 2k + k − 1 P 2n

We inject (39), (38) and (41) into (37), and we obtain the relation
√ Z ∞ −y
e Pλ Sλ
2v(0) √ = (u1 + √ ). (42)
0 y P 2

For λ ∈ (−1/2, 0), it is easy to√see that Sλ√> S0 ∈ (0, +∞). Therefore, if we assume by
contradiction that u1 > −S0 /2 2 > −Sλ /2 2, (42) implies v(0) > 0.

the sequel, λ ∈ (−1/2, 0) is fixed, and we assume by contradiction that u1 >

In √
−S0 /2 2, and thus that v(0) > 0. We will show that this implies that v(x) → −∞
as x → +∞, which is a contradiction with v ∈ L2 .

27
 Let us first define, for x ∈ R,
x
w(x) := v(x) cosh √ . (43)
2
Then for x ∈ R, we have
w(x)
v(x) = ,
cosh √x2
x
′ w ′ (x) 1 tanh √2
v (x) = −√ x w(x),
cosh √x2 2 cosh √2
x
′′ w ′′ (x) 2 tanh √2 ′ w(x) 1 w(x)
v (x) = x −√ x w (x) − + .
cosh √2 2 cosh √2 3 x
cosh √2 2 cosh √x2
Since v satisfies (33), it follows that w solves
√ x 1 1
w ′′(x) = 2 tanh √ w ′(x) − √ x
− ( + λ)w(x). (44)
2 2 cosh √2 2

Remark that w(0) = v(0) > 0, w is even and w ′ (0) = 0. Evaluating (44) at x = 0 and
using the fact that λ > −1/2, we thus have
1 1
w ′′ (0) = − √ − ( + λ)w(0) < 0.
2 2

Therefore there exists η0 > 0 such that if x ∈ [0, η0 ], w ′′ (x) < 0 and w(x) > 0. Thus we
can define

x1 := sup{x > 0, w(y) > 0 and w ′′ (y) < 0 for every y ∈ [0, x]}.
If x1 = +∞, we have w ′′(y) < 0 on [0, ∞) which implies that w ′ strictly decreases on R+ ;
thus for x > 1, w ′ (x) 6 w ′(1) < w ′ (0) = 0. Therefore for x > 1, w(x) 6 w ′ (1)(x−1)+w(1).
In particular, w(x) < 0 for x large enough, which is a contradiction with the assumption
x1 = +∞. Thus x1 < ∞.
Next, we show that w(x1 ) = 0 or w ′′ (x1 ) = 0. Indeed, if it was not the case, thanks
to the continuity of w and w ′′ , there would exist a neighbourhood of x1 in R+ on which
w > 0 and w ′′ < 0, which would yield a contradiction with the definition of x1 . We next
prove that w ′′ (x1 ) < 0, and thus w(x1 ) = 0. Indeed,
√ x1 1 1
w ′′ (x1 ) = 2 tanh √ w ′(x1 ) − √ x
− ( + λ)w(x1 )
2 2 cosh √12 2
1
6 −√ < 0,
2 cosh √x12

28
because w ′ is strictly decreasing on [0, x1 ] (thus w ′(x1 ) < w ′(0) = 0), 1/2 + λ > 0 and
w ′ (x1 )
w(x1 ) > 0. Therefore v(x1 ) = w(x1 ) = 0, and v ′ (x1 ) = cosh x
√1
< 0. It follows that there
√ 2

exists η1 > 0 such that −1/ 2 < v(x) < 0 for x ∈ (x1 , x1 + η1 ]. Thus we can define
√
x2 := sup{x > x1 , −1/ 2 < v(y) 6 0 for y ∈ [x1 , x)} > x1 + η1 > x1 .
√
For x ∈ [x1 , x2 ), since −λ > 0, v(x) 6 0 and 1/ 2 + v(x) > 0, we have
 
′′ 1 1
v (x) = − √ + v(x) − λv(x) < 0.
cosh2 √x2 2

Therefore v ′ is decreasing on [x1 , x2 ), and√

v(x) 6 v ′ (x1 )(x − x1 ) if x ∈ [x1 , x2 ). Now, since
v ′ (x1 ) < 0, we have x2 < ∞, v(x2 ) = −1/ 2 and v ′ (x2 ) < v ′ (x1 ) < 0.
From now on, we distinguish the two cases:

√ 1
1. x2 > xλ := 2arccosh √−λ ,
2. x2 < xλ .
1
In the first case, for x > x2 , cosh2 x
√
+ λ 6 0. v ′ (x2 ) < 0, therefore there exists η2 > 0,
2
v ′ (x) < 0 for x ∈ [x2 , x2 + η2 ], and the definition of x3 > x2 as

x3 := sup{x > x2 , v ′ (y) < 0 for y ∈ (x2 , x)} > x2

√
makes sense. v is decreasing on  v(x) 6 −1/ 2 for every x ∈ [x2 , x3 ). It
 [x2 , x3 ), thus
follows that for x ∈ [x2 , x3 ), − cosh12 √x + λ v(x) 6 0, and
2

!
1 1 1 1 1
v ′′ (x) = − √ − + λ v(x) 6 − √ < 0.
2 cosh2 √x
2
cosh2 √x
2
2 cosh2 √x
2

Therefore v ′ is decreasing on [x2 , x3 ). This implies that x3 = +∞, and for x > x2 ,
v(x) 6 v(x2 ) + v ′ (x2 )(x − x2 ), thus v(x) → −∞ as x → +∞.
We are now concerned with the second case x2 < xλ . We introduce

x4 := sup{x ∈ (x2 , xλ ), v ′(y) < 0 for y ∈ [x2 , x)}

and the open set Ω := {x ∈ (x2 , x4 ), v ′′(x) > 0}. Ω is the disjointed union of its connex
parts ωi = (ai , bi ) for i ∈ I, where I is a finite or countable set. For x ∈ ωi , v ′′ (x) > 0, and
thus integrating twice, for x ∈ ωi ,

v ′ (x) > v ′ (ai ),

v(x) − v(ai ) > v ′ (ai )(x − ai ).

29
Using (33) and the fact that x < xλ , it follows that
!
1 1 1
v ′′ (x) 6 ′
2 √x + λ (−v(ai ) − v (ai )(x − ai )) −
√ ,
cosh 2 2 cosh2 √x
2

which may be integrated between ai and x. After an integration by parts, we obtain

 
′ ′
√ x ai
v (x) 6 v (ai ) − ( 2v(ai ) + 1) tanh √ − tanh √ − λv(ai )(x − ai )
2 2
√ Z x√
v ′ (ai ) 2 ′ x ′ y
−λ (x − ai ) − v (ai ) 2 tanh √ (x − ai ) + v (ai ) 2 tanh √ dy. (45)
2 2 ai 2
√
For x ∈ (x√2 , x4 ), v ′ (x) < 0, thus v(x) < v(x2 ) = −1/ 2. ′
√ In particular,′′ v (ai ) < 0 and
′′
v(ai ) + 1/ 2 < 0 (indeed, ai > x2 because v (x2 ) = λ/ 2 < 0 while v (ai ) = 0). For
x ∈ ωi , it follows from (45), the mean value theorem and the inequality tanh √x2 6 1 that
" #
1 1
v ′ (x) 6 v ′ (ai ) − (v(ai ) + √ ) + λv(ai ) (x − ai )
2 cosh2 ai
√
2
| {z }
=−v′′ (ai )=0.

v ′ (ai ) 2 ′
√ Z x y
−λ (x − ai ) − v (ai ) 2 (1 − tanh √ )dy
2 ai 2
 Z x 
′
√ y λ(x − ai ) 2
= v (ai ) 1 − 2 (1 − tanh √ )dy −
ai 2 2
 Z 
√ ∞
y √
6 v ′ (ai ) 1 − 2 (1 − tanh √ )dy = v ′ (ai )(1 − 2 ln(1 + e− 2x2 )). (46)
x2 2
√
We provisionally admit the fact that (1 − 2 ln(1 + e− 2x2 )) > 0. Let us assume by
contradiction that v ′ (x4 ) = 0.
If x4 6∈ Ω, there exists η > 0 such that v ′′ (x) 6 0 for x ∈ [x4 − η, x4 ], and thus
0 = v ′ (x4 ) 6 v ′ (x4 − η) < 0, which is a contradiction.
Thus x4 ∈ Ω. Let (cn )n be a sequence in Ω such that cn → x4 . For each n, cn ∈ ωi(n) ,
and cn < bi(n) 6 x4 . Thus bi(n) → x4 . If bi(n) was constant (equal to x4 = bi0 ) for n > n0 ,
we would have by (46)
√
0 = v ′ (x4 ) = v ′ (bi0 ) 6 v ′ (ai0 )(1 − 2 ln(1 + e− 2x2
)) < 0,

which is a contradiction. Thus bi(n) takes infinitely many values. We next remark that if
y ∈ (0, xλ ) satisfies v ′′ (y) = 0, it follows from (33) that
1
v(y) = − √ =: f (y).
2(1 + λ cosh2 √y )
2

30
Here, we have 0 = v ′′ (bi(n) ) → v ′′ (x4 ), thus v ′′ (x4 ) = 0, and v(x4 ) = f (x4 ). Next,

′
sinh √x42 cosh √x42 ′
f (x4 ) = λ  2 < 0 = v (x4 ).
2 √
1 + λ cosh x42

This implies that there exists η4 > 0 such that v(x) < f (x) for x ∈ [x4 −η4 , x4 ). For n large
enough, bi(n) ∈ [x4 − η4 , x4 ), v ′′ (bi(n) ) = 0, thus v(bi(n) ) = f (bi(n) ). This is a contradiction.
If x4 < xλ , the definition of x4 implies that we must have v ′ (x4 ) = 0, which we have
just seen to be impossible. Thus x4 = xλ and v ′ (x4 ) < 0. Then, we conclude as in the case
x2 > xλ , with x2 replaced by x4 .

To complete the proof, it remains to prove the estimate

 √ 
2 ln 1 + e− 2x2 < 1 (47)

in the case x2 < xλ , which will be done by bounding x2 from below. For x ∈ [0, x2 ],
λ + cosh12 √x > 0, and we have seen that v is decreasing on [0, x2 ], thus v(x) 6 v0 = v(0) for
2
x ∈ [0, x2 ], and
!
1 1
v ′′ (x) = − λ + v(x) − √
cosh2 √x
2
2 cosh2 √x
2
!
1 1
> − λ+ 2 √x v0 − √
cosh 2 2 cosh2 √x
2
 
1 1
> − v0 + √ ,
2 cosh2 √x2

because −λv0 > 0. Integrating this inequality between 0 and x ∈ [0, x2 ], we get
√ x √
|v ′(x)| = −v ′ (x) 6 (1 + v0 2) tanh √ 6 (1 + v0 2).
2
Using the mean value Theorem, it follows that
1 √
√ + v0 = |v(x2 ) − v(0)| 6 (1 + v0 2)x2 .
2
√
Therefore x2 > 1/ 2, and
√
2 ln(1 + e− 2x2
) 6 2 ln(1 + e−1 ) < 1.

This completes the proof. 

31
7 Dynamical stability
In this section, we numerically investigate the dynamical stability of the black solitons.
For this purpose, we use finite differences in both time and space.

7.1 Description of the numerical method

Space and time steps h and δt being given, all the computations will be made in a rect-
angular grid at spatial points xj = jh and at discrete times tn = nδt. Let us define unj as
the approximate value of u(xj , tn ). For any given smooth function w, we also define the
operators Lx as

1   ∂2w 
(Lx w)j = 2 wj+1 − 2wj + wj−1 = 2
+ O(h2 )
h ∂x j

that can be seen as the discretization of the differential operator ∂ 2 /∂x2 with use of Taylor
formulas. We then consider the symmetric Crank-Nicolson discretization of (1): setting
un+1/2 = (un + un+1 )/2, the system writes
i n+1 n+1/2 2 n+1/2
(uj − unj ) + (Lx un+1/2 )j + f (|uj | )uj =0
δt
for n > 0 and |j| 6 J, with J given. This is a nonlinear algebraic discrete system that is
solved at each time step using a fixed point method. This way of discretizing the system
enables us to have the conservation of the discrete first invariant on the whole space
X 
In = |unj |2 − 1 = I0 , ∀ n ≥ 0.
j∈Z

7.2 The boundary problem

When dealing with asymptotic behaviour of time evolution problems, numerics may become
delicate since for obvious CPU reasons, computations have to be made in a bounded
interval of R, often referred to as the computational domain. In this case, an inappropriate
boundary condition can lead to artificial reflection and may cause a wrong approximation
of the exact solution. In all the cases that will be investigated, the solutions under study
tend to different limit values at ±∞. The most natural way to proceed is to set these
values for the numerical approximation at each extremity of the computational domain.
This is justified because given a black soliton ϕ; the Cauchy problem for (1) is well-posed
in ϕ + H 1 (see [G2]). We show here that even in the linear case

∂u ∂ 2 u
i + 2 =0 (48)
∂t ∂x
with boundary conditions −1 and 1 at −∞ and ∞, this treatment can seriously affect the
solution asymptotics. In order to avoid such problems, a numerical trick consists in adding

32
a diffusion term in (48), that will come into play near the boundary. We thus replace (48)
by

∂u ∂ 2 u
(i − α(x)) + 2 = 0, (49)
∂t ∂x
where the function α vanishes at the interior of the numerical domain. This means that the
original problem is turned into a diffusive-like equation where the diffusion term only plays
a role close to the boundary. Consequently, waves that could reflect on the boundary are
absorbed. The term α is referred to in the literature as a sponge factor. Such a function
is plotted in Figure 6. Remark that since we deal with a modified equation, conservation
of invariants for initial problem no more holds.

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
−50 −40 −30 −20 −10 0 10 20 30 40 50
x

Figure 6: Plot of α.

We now show the results obtained for the numerical resolution of (48) considered with the
initial data u0 (x) = tanh x, for the following choice of grid parameters: h = 0.1, J = 200,
δt = 0.1 and a sponge factor given by

α(x) = exp(−(|x| − L)2 /4),

where L = Jh. In Figures 7 and 8, we view the solution computed with the sponge
factor compared to both solution computed with nonhomogeneous Dirichlet conditions
u(t, −L) = −1, u(t, L) = 1 and solution considered in the larger domain ] − 3L, 3L[,
for which the solution is not affected by the boundary on the small domain. It can be
observed that the first solution is a much better approximation of the reference one than
the solution obtained with Dirichlet conditions. In all what follows, such sponge terms will

33
be used for our tests. Note that in this case, there is no need in using a numerical scheme
which preserves the charge on the numerical domain since the modified equation is no more
conservative, even if the lack of conservation is only relevant close to the boundary.
re(u(t,x)) at time t = 15 re(u(t,x)) at time t = 15
1.5

Dirichlet conditions 1.4 Dirichlet conditions

Sponge factor Sponge factor
Solution on a larger domain Solution on a larger domain
1 1.3

1.2

0.5
1.1

1
0

0.9

−0.5 0.8

0.7
−1
0.6

0.5
−1.5
−20 −15 −10 −5 0 5 10 15 20 4 6 8 10 12 14 16 18 20
x x

Figure 7: Profile of the real part of Figure 8: Profile of the real part of
the solutions of (48) with Dirichlet the solutions of (48) with Dirichlet
boundary conditions, with sponge fac- boundary conditions, with sponge fac-
tor and computed on the larger domain tor and computed on the larger domain
] − 60, 60[, plot on the small domain ] − 60, 60[, T = 15, zoom near the right
] − 20, 20[, T = 15. boundary.

7.3 The pure power defocusing equation

Since the stability condition stated in Theorem 1.1 deals with linear stability, we first
perform computations for the linearized equation

∂u ∂ 2 u
i + 2 + f (ϕ2 )u + 2ϕ2 f ′ (ϕ2 )Reu = 0, (50)
∂t ∂x
with f (r) = 1 − r σ (which is an other way to write system (3)) in order to check if a
stability threshold occurs. We start from the initial condition
2
ε(x) = e−0.01x , (51)
that can be seen as a perturbation of zero. Note that due to the linearity of (50), this
perturbation is not needed to be small. We numerically solve equation (50) using the same
semi-implicit discretization of the linear term as in the nonlinear case and plot the evolution
of the maximal amplitude in order to look for the generation of an unstable mode. For the
sake of clarity, we have preferred to view the profile of k(t) := log(1 + ku(t)kL∞ )/t: if k

34
tends to zero as t is large, then no exponential growth is obtained and the black soliton can
be considered as linearly stable. If k tends to some positive constant λ, then instability
occurs with the corresponding growth rate given by λ.

Our first simulations showed successive drops of k, caused by periodic changes of the
sign of the real part. It has to be noticed that this phenomenon still occurs if the sponge
term is removed. Surprisingly, dealing with well-adapted grid parameters enables us to
delay it. In most of the tests that will be discussed, we have decided to deal with large
numerical domains considered with small space step h, for the following reasons: firstly, if
we look for an unstable mode for which k tends to a positive limit value, the space domain
has to be large enough in such a way that the eigenvector is well-localized. Furthermore,
the black soliton ϕ computed using a differential solver (see Section 5) has to be accurately
described around the origin where the variations of ϕ are most significant: this requires to
work with reasonably small values of h.

We have computed the evolution of k for different values of σ (σ = 1, σ = 2 and σ = 5)

with parameters J = 50000, h = 0.002, δt = 0.2 and N = 10000. Figure 9 indicates that k
decreases to 0: there is no positive eigenvalue for the linearized operator, which illustrates
linear stability for any value of σ.

−3
x 10

sigma = 1
18 sigma = 2
sigma = 5
16

14
log(1+||u(t)||∞)/t

12

10

0
200 400 600 800 1000 1200 1400 1600 1800
t

Figure 9: Plot of k versus time, σ = 1, σ = 2 and σ = 5.

We then numerically solved problem (31), injecting as initial data a small perturbation
of the black soliton numerically calculated: we take u0 = ϕ + ε, with
2
ε(x) = q(1 + i)e−0.01x cos x. (52)

35
 We took the following parameters : h = 0.01, J = 20000, δt = 1, N = 1000 and
q = 3.10−4. For any value of σ, the numerical solution seems to travel with a velocity
depending on σ. This may suggest that the black soliton is not stable. However, we
translate the solution in such a way that the real part changes its sign at the origin and we
study the L∞ error with respect to the black soliton. We plot in Figure 10 the evolution of
the L∞ norm of the error between the black soliton and the translated numerical solution
obtained for σ = 1, with respect to time. It can be observed that this difference is bounded
through time, which is a grant of stability. Other tests performed for different values of
σ have shown a similar behaviour. Note that it can be observed small oscillations of the
error caused by the translation of the solution from an integer number of grid cells and
the large value of the derivative of the black soliton at x = 0. Thus, the period observed
in Figure 10 corresponds to the time the solution takes to go through a grid cell.

−3
x 10
4

3.5

2.5

1.5

0.5

0
0 100 200 300 400 500 600 700 800 900 1000
t

Figure 10: Plot of the L∞ error between the translated solution and the exact black soliton
versus time, σ = 1.

2
We then perturbed the initial data with an oscillatory function ε(x) = 0.5e−x sin 5xϕ(x).
In our test, we set J = 800, h = δt = 0.05 and for simulations performed until final time
T = 5. In Figures 11 is plotted the solution profile at successive times t = 1, t = 3 and
t = 5. It can be observed a propagation to the right of the wavetrain, with a similar effect
to the left side for negative x. Since its length seems quite constant through time, the
global L2 error between the numerical solution and the black soliton will remain bounded.
Since it is worth investigating the behaviour at the left of the perturbed wave, we also
define the local L2 error, that is computed around the origin (in our simulations, we chose
the domain D = [−10, 10]). In Figure 12 are plotted the two errors and it can be noticed
that once the perturbation has left D, then the approximate solution mimics the soliton
profile and the local L2 error becomes small, while the global error stays bounded since

36
the travelling perturbation keeps the same amplitude. This effect could be referred as local
asymptotic stability, meaning that in compact domains of the real line, the black soliton
attracts the solutions for large times.
1.1 0.15
t=1 Whole domain
t=3 Small domain
t=5

0.1

0.05

0.9 0
0 5 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x time

Figure 11: Snapshots of the perturbed Figure 12: Evolution of global and local
solution at different times (t = 1, t = L2 error between the perturbed solu-
3 and t = 5) around the black soliton tion and the black soliton versus time.
profile, σ = 1.

These two simulations have shown different behaviours for the solution of (1): the first
perturbation leads us to a travelling wave that remains close to the black soliton in L∞
norm, whereas the second perturbation generates local oscillations going away to infinity
and a stationary global profile (that is with no translation effect). This could suggest that
the perturbed solution mimics travelling solitons of (1) with velocity depending on the
initial perturbation, the velocity being equal to zero in the last experiment. Such states
are referred as travelling bubbles and are explictly given in the case σ = 1 by
r r
v2  v 2 x − vt  v
ϕv (x − vt) = 1 − tanh 1− √ + i√ , (53)
2 2 2 2
see for instance [G0]. Note that for v 6= 0, the phase shift of the bubble differs from the
black soliton one. In particular, the numerical boundary conditions u = ±1 at x = ±xmax ,
(xmax = Jh) are not well-adapted to describe these bubbles. However, it is possible to
compare the numerical solution with such a bubble on a subdomain D ⊂] − xmax , xmax [:
first, up to a multiplication by a constant phase term exp iθ, the imaginary part of exp iθ unj
is observed to be constant, say equal to γ (that does not depend √ on time), on D. The
numerical value of the velocity is then deduced from (53) by v = 2γ and can be compared
to the velocity of the solution deduced from Figure 10 by the relation v = ∆d/∆t, where
∆d stands for the distance covered by the bubble between two local extrema of the L∞
profile (taking into account that the distance h is covered during a period) and ∆t is

37
the corresponding√ time. For instance, the computations made for perturbation (51) with
q = 3 10 give 2γ ≃ 2.086 10 and show a good agreement with ∆d/∆t ≃ 2.083 10−4
−4 −4

(the values are chosen in such a way that ∆t ∼ T in the plot). For σ > 1, the explicit form
of travelling bubbles is not known but it is still possible to compute the velocity v of the
bubble from the minimum intensity η0 of the solution: indeed, setting

(r − 1)2
Vv (r) := V (r) − v 2 ,
4r
η0 is the unique zero of Vv on the interval (0, 1) (see [G1]); consequently, v expresses as
p
V (η0 )η0
v=2 , (54)
1 − η0
the value η0 being determined from the numerical solution. Starting from the same pertur-
bation as in the previous test for the resolution of (31) with σ = 2, (54) gives v ≃ 2.77 10−4
which is again very close to the observed velocity that has been found close to 2.74 10−4 .
Other experiments have shown that this still holds for other values of σ and other ini-
tial perturbations. Furthermore, the velocity of the numerical travelling bubble essentially
depends on the perturbation profile at the origin.

7.4 The saturated equation

As in the previous section, we first investigate the linearized equation (50), starting from
(51) as initial data. We used the following parameters : J = 3 105 , h = 1/3 10−2 , δt = 10
and N = 1000. The choice of such a δt is justified by the smallness of the nonlinearity
which prescribes large characteristic time scales. We also plot in Figures 13, 14, 15 and
16 snapshots of the real and imaginary parts of the solution of (50) at different times in
the two cases a = 7.4 and a = 7.6 (note that 7.4 < aVK < 7.6). It can be observed that
in the case a = 7.4, the imaginary part spreads out in the numerical domain but keeps
the same order of magnitude through time. The solution we obtain here closely looks like
the derivative of the black soliton, which is not that surprising since L2 ϕ′ = 0, meaning
that ϕ′ is a stationary solution to (50). On the contrary, the profile of the imaginary part
obtained for a = 7.6 is amplified without any alteration of its shape, after a transient
regime at moderate times. This suggests that in the latter case, an unstable mode has
been captured from initial data (51). It points out the existence of a threshold value al for
which we observe both dynamical linear instability if a > al and stability if a < al . Note
that al ≃ aVK .
In order to study more precisely the amplification rate of the solution, k has been
plotted in the cases a = 7, a = 8 and a = 9 (see Figure 17). It gives another illustration of
linear instability of the black soliton above al . In figure 18, we have plotted k∞ := lim k(t)
t→∞
versus a for a > al . In each case, this limit has been estimated by extrapolation from
computations until t = 104 . It has to be pointed out that computing the value of the
instability growth rate is delicate when a becomes too close to the critical threshold: the

38
 25
t = 1000 t = 1000
t = 2000
0.1 t = 2000
t = 3000
t = 3000
20
0.05

15 0

−0.05
10

−0.1

5
−0.15

−0.2
0
5 10 15 20 25 30 35 100 200 300 400 500 600 700 800 900 1000
x
x

Figure 13: Plot of the real part at dif- Figure 14: Plot of the imaginary part
ferent times, a = 7.4. at different times, a = 7.4.

700
t = 1000 t = 1000
t = 2000 t = 2000
t = 3000 t = 3000
600 6

500 4

400
2

300
0

200
−2

100
−4

0
0 5 10 15 20 25 30 35 0 100 200 300 400 500 600 700
x x

Figure 15: Plot of the real part at dif- Figure 16: Plot of the imaginary part
ferent times, a = 7.6. at different times, a = 7.6.

code only gives an approximate solution that shares the same properties of the exact one
up to a numerical tolerance. Moreover, the nonlinearity is always very small for values of
a close to al , which requires our experiments to be made for sufficiently large times in such
a way that the nonlinear contribution can be observed.
We also solved Equation (32) starting from the Gaussian modulated initial perturbation
(52). In our computations, we chose δt = 1, J = 20000, h = 0.01, q = 3.10−4 in the cases
a = 5 and a = 9 (see Figures 19 and 20). For a = 5, it can be noticed that the profile
of the real part cannot be distinguished from the black soliton. However, the solution has
been translated to the right: the real part now vanishes at x = 0.12. In the case a = 9,
the situation completely changes: the translation to the right is now much faster and the
profile of the translated real part differs from the black soliton.

39
 −3
x 10
12 0.02
a=7
a=8 0.018
a=9
10
0.016

0.014
8

0.012
k(t)

k∞
6 0.01

0.008

4
0.006

0.004
2

0.002

0
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12
t a

Figure 17: Plot of k versus time, a = 7, Figure 18: Plot of k∞ as a function of

a = 8 and a = 9. a.

1.5 1.5
Solution at final time Solution at final time
Black soliton Black soliton

1
1

0.5

0.5

−0.5

−0.5
−1

−1.5 −1
−200 −150 −100 −50 0 50 100 150 200 −200 −150 −100 −50 0 50 100 150 200
x x

Figure 19: Profile of the real part of the Figure 20: Profile of the real part of the
solution at final time T = 1000 com- solution at final time T = 1000 com-
pared with the black soliton, a = 5. pared with the black soliton, a = 9.

In order to see more precisely the influence of parameter a, we view in Figure 21 on the
same plot the errors between the black soliton and the translated solution of (32) obtained
for different values of a: a = 6, a = 7, a = 8 and a = 9, with the same values of both grid
parameters and initial perturbation. As in the test performed in the linearized case, this
figure suggests the occurrence of a threshold value anl ∈ (7, 8): taking a < anl leads us to a
bounded variation of the error, whereas for a > anl the error seems to increase with time.
Note that the value of anl is in agreement with the threshold values obtained on the one
hand with the Vakhitov-Kolokolov function (aVK ) and on the other hand with the study
of the linearized system (al ). In the cases where the translated profiles do not differ much
from the black soliton (that is for a = 6 and a = 7), the oscillations that have been already

40
observed in the pure power case (see Figure 10) remain visible for large times.

0.2

a=6
0.18
a=7
a=8
0.16
a=9

0.14

0.12

0.1

0.08

0.06

0.04

0.02

0
100 200 300 400 500 600 700 800 900 1000

Figure 21: Plot of the error between the black soliton and the translated solution of (32)
a = 6, a = 7, a = 8 and a = 9.

In the stable case, the moving soliton can be identified as a travelling bubble, as for
the pure power nonlinearity. Following the same strategy as before, we have measured the
velocity v = 9.86 10−5 of the solution for a = 2, which gives a good agreement with respect
to the velocity v = 9.73 10−5 obtained by (54).

7.5 Conclusion
We have investigated the stability of black solitons for different nonlinearities using different
points of view. First, the Vakhitov-Kolokolov function g has been computed numerically.
Thanks to Theorem 1.1, linear stability holds if and only if g(λ) has a negative limit
when λ tends to zero. Secondly, the linearized equation has been solved starting from
a perturbation of the black soliton. Linear instability manifests itself by the emergence
of an exponentially increasing unstable mode. Finally, the fully nonlinear evolution has
been computed and has shown that even in the linearly stable cases, the black soliton
may be translated with a velocity depending on the initial perturbation and behaves like a
travelling bubble on a spreading domain. The evaluation of the error between the computed
solution and the initial data does not seem well-adapted for the stability analysis. The
stability has to be studied modulo translations: this is exactly orbital stability.

41
 Our numerical simulations bring to the fore that in the pure power case, the black
soliton is linearly stable and is also orbitally stable. In the saturated equation, these three
approaches tend to exhibit the same threshold value aVK ≃ al ≃ anl for the parameter
a. Let us denote by ac this threshold value. If a < ac , the black soliton is stable in
all the above mentioned meanings, whereas instability occurs when a > ac . The need to
investigate orbital stability instead of stability in usual sense may be surprising, taking
into account previous results in the literature concerning travelling bubbles for the ψ 3 − ψ 5
nonlinear Schrödinger equation (see [BP], [dB]).

8 Appendix
We prove here some technical lemmas which have been used in section 4. All these proofs
are slightly modified versions of proofs given in [GJLS].

Proof of Lemma 4.1. Note first that since λ 6∈ iR, −λ2 6∈ R+ = σ(−∂x2 (−∂x2 + c2 )).
Thus the operator λ2 − ∂x2 (−∂x2 + c2 ) is invertible, and Bλ−1 may be expressed as

 
λ(λ2 − ∂x2 (−∂x2 + c2 ))−1 −∂x2 (λ2 − ∂x2 (−∂x2 + c2 ))−1
Bλ−1 = . (55)
−(−∂x2 + c2 )(λ2 − ∂x2 (−∂x2 + c2 ))−1 λ(λ2 − ∂x2 (−∂x2 + c2 ))−1

Thus the lemma reduces to prove uniform bounds on each four operators on L2 in the
right-hand side of (55). Passing into Fourier variables, we get

|λ|
kλ(λ2 − ∂x2 (−∂x2 + c2 ))−1 kL(L2 ) 6 sup 2 2 2 2
,
ξ∈R |λ + ξ (ξ + c )|

ξ2
k − ∂x2 (λ2 − ∂x2 (−∂x2 + c2 ))−1 kL(L2 ) 6 sup 2 2 2 2
ξ∈R |λ + ξ (ξ + c )|

and

ξ 2 + c2
k(−∂x2 + c2 )(λ2 − ∂x2 (−∂x2 + c2 ))−1 kL(L2 ) 6 sup 2 2 2 2
.
ξ∈R |λ + ξ (ξ + c )|

Next,

|λ| |λ| (a2 + τ 2 )1/2 (a2 + 1)1/2

sup 2 2 2 2
6 = 6 .
ξ∈R |λ + ξ (ξ + c )| |Im(λ2 )| 2|a||τ | 2|a|
Similarly,

c2 c2 c2
6 6 .
|λ2 + ξ 2 (ξ 2 + c2 )| 2|a||τ | 2|a|
Finally, if ξ 2 (ξ 2 + c2 ) > 2τ 2 , we also have ξ 2 (ξ 2 + c2 ) − τ 2 > ξ 2 (ξ 2 + c2 )/2, and

42
 ξ2 ξ2
6 6 C,
((a2 + ξ 2 (ξ 2 + c2 ) − τ 2 )2 + 4a2 τ 2 )1/2 a2 + ξ 2 (ξ 2 + c2 )/2
while if ξ 2 (ξ 2 + c2 ) 6 2τ 2 ,

 1/2 √
ξ2 ξ2 1 ξ 2 (ξ 2 + c2 ) 2
2 2 2 2 2 2 2 2 1/2
6 6 6 ,
((a + ξ (ξ + c ) − τ ) + 4a τ ) 2|a||τ | 2|a| τ2 2|a|

which completes the proof of the lemma. 

Proof of Lemma 4.2. We first prove the result for j = 3. We write

Z ∞
(x2 + c2 )q
kgλ3 kqLq = 2 dx
[(a2 + x2 (x2 + c2 ) − τ 2 )2 + 4a2 τ 2 ]q/2
0
= 2(I1 + I2 ),

where

Z √
c/ 2
(x2 + c2 )q
I1 :=
0 [(a2 + x2 (x2 + c2 ) − τ 2 )2 + 4a2 τ 2 ]q/2
Z √
c/ 2
(3c2 /2)q
6 dx −→ 0,
0 (2|a||τ |)q |τ |→∞

and

Z ∞
1 (x2 + c2 )q
I2 := √   2 2 2 2 2 2 q/2 dx.
c/ 2 (2|a||τ |)q
1 + a +x 2|a||τ
(x +c )−τ
|

x2 (x2 +c2 )+a2 −τ 2 τ 2 +c4 /4−a2

The change of variable y = 2|a||τ |
yields, setting τ ′ = 2|a||τ |
, if τ 2 > a2 −c4 /4+
1,

 q
Z ∞ (τ 2
+ c 4
/4 − a 2 q/2−3/4
) ( y
+ 1) 1/2
+ c2
1 1 τ ′ 2 4
2(τ +c /4−a ) 2 1/2
I2 = q−1  1/2 dy.
4 c4
−τ ′ (2|a||τ |) 2 q/2 y 1/2 c2 y 1/2
2|a||τ | (1 + y ) ( τ ′ + 1) − 2(τ 2 +c4 /4−a2 )1/2 ( τ ′ + 1)

c4 c2
Next, y > 2|a||τ |
− τ ′ implies ( τy′ + 1)1/2 > (τ 2 +c4 /4−a2 )1/2
, thus for q > 1,

43
 
1/2 q
Z ∞ 3 y
+1
1 (τ 2 + c4 /4 − a2 )q/2−3/4 2 τ′
I2 6 
1/2 1/2
1/2 dy
4 (2|a||τ |)q−1 −τ ′ 1 y y
(1 + y 2)q/2 2 τ′
+1 τ′
+1
Z ∞ y

q/2−3/4
−1/2 +1 τ′
6 C|τ | dy
−τ ′ (1 + y 2 )q/2
Z τ′ Z ∞ !
q/2−3/4
dy y
6 C|τ |−1/2 2 q/2
+ ′q/2−3/4 (1 + y 2 )q/2
dy
−τ ′ (1 + y ) τ′ τ
 Z ∞ 
−1/2 ′ 1−q −q/2−3/4
6 C|τ | 1 + |τ | y dy → 0.
1 |τ |→∞

We now perform similar computations for j = 2: if τ > τ0 where τ0 > |a2 − c4 /4|1/2
large enough, in such a way that for every τ > τ0 , c4 /(8|a||τ |) − τ ′ < −τ ′ /2 < 0 and
is √
1/ 2 > c2 /(τ 2 + c4 /4 − a2 )1/2 , we compute

Z ∞
(a2 + τ 2 )q/2
kgλ2 kqLq = 2 dx
[(a2 + x2 (x2 + c2 ) − τ 2 )2 + 4a2 τ 2 ]q/2
0
Z
(a2 + τ 2 )q/2 1 ∞
= 2 (τ 2 + c4 /4 − a2 )−3/4 (1 + y 2)−q/2
(2|a||τ |)q−1 4 8|a||τ
c4
|
−τ ′

 1/2 −1/2  −1/2

y c2 y
× + 1 − + 1 dy
τ′ 2(τ 2 + c4 /4 − a2 )1/2 τ′
(a2 + τ 2 )q/2
= (J1 + J2 ),
2(2|a||τ |)q−1(τ 2 + c4 /4 − a2 )3/4
where
Z ∞
dy
J2 := 
1/2 1/2
1/2
−τ ′ /2 y c2 y
(1 + y 2 )q/2 τ′
+1 − 2(τ 2 +c4 /4−a2 )1/2 τ′
+1
Z ∞
dy
6 C
−∞ (1 + y 2)q/2
and

Z −τ ′ /2
dy
J1 := 
1/2 1/2
1/2
c4
−τ ′ y c2 y
8|a||τ | (1 + y 2 )q/2 τ′
+1 − 2(τ 2 +c4 /4−a2 )1/2 τ′
+1
Z −τ ′ /2
1 dy
6 
q/2 
1/2 1/2
1/2
c4
−τ ′ τ′ 2 y c2 y
8|a||τ | 1+ 2 τ′
+1 − 2(τ 2 +c4 /4−a2 )1/2 τ′
+1

44
 y

1/2 c2
We make the change of variables z = τ′
+1 − 2(τ 2 +c4 /4−a2 )1/2
, in such a way that for
|τ | > τ0 ,
Z √
1/ 2
−q |τ ′ |dz
J1 6 C|τ | 1/2
6 C|τ |1−q ,
0 z
which completes the proof of Lemma 4.2. 

References
[B] I.V. Barashenkov, Stability criterion for dark solitons, Phys. Rev. Lett., 77 (1996),
1193-1197.

[BGMP] I.V. Barashenkov, A.D. Gocheva, V.G. Makhankov, I.V. Puzynin, Sta-
bility of the soliton-like “bubbles”, Physica D, 34 (1989), 240-254.

[BP] I.V. Barashenkov, E.Y. Panova, Stability and evolution of the quiescent and trav-
eling solitonic bubbles, Physica D, 69 (1993), 114-134.

[BS] F.A. Berezin, M.A. Shubin, The Schrödinger equation, Kluwer Academic Publish-
ers (1991).

[dB] A. de Bouard, Instability of stationnary bubbles, SIAM J. Math. Anal. 26 (3) (1995)
566-582.

[G0] C. Gallo, Propriétés qualitatives d’ondes solitaires solutions d’EDP non linéaires
dispersives, Ph. D. Thesis, Paris-Sud Orsay University (2005).

[G1] C. Gallo, The dark solitons of the one-dimensional nonlinear Schrödinger equation,
submitted.

[G2] C. Gallo, The Cauchy Problem for defocusing Nonlinear Schrödinger equations with
non-vanishing initial data at infinity, submitted.

[GJLS] F. Gesztesy, C.K.R.T. Jones, Y. Latushkin, M. Stanislavova, A Spectral

Mapping Theorem and Invariant Manifolds for Nonlinear Schrödinger Equations, In-
diana University Math. Journal 49 (2000) 221-243.

[H] D. Henry, Geometric theory of semilinear parabolic equations, Lecture Notes in Math-
ematics 840, Springer, Berlin (1981), p. 136.

[HPW] D.B. Henry, J.F. Perez, W.F. Wreszinski, Stability Theory for Solitary-Wave
Solutions of Scalar Field Equations, Commun. Math. Phys. 85 (1982), 351-361.

[KK] Y.S. Kivshar, W. Krolikowski, Instabilities of dark solitons, Opt. Lett. 20 (1995),
1527-1529.

45
[KLD] Y.S. Kivshar, B. Luther-Davies, Dark solitons: physics and applications,
Physics Reports 298 (1998), 81-197.

[MS] D.J. Mitchell, A.W. Snyder, Stability of nonlinear guided waves, J. Opt. Soc.
Am. B 10 (1993), 1572-1580.

[P] J. Prüss, On the spectrum of C0 -semigroups, Trans. Amer. Math. Soc. 284, 2 (1984),
847-857.

[L] Z. Lin, Stability and instability of traveling solitonic bubbles, Adv. Diff. Eq., 7 (2002),
897-918.

46