## Abstract

A direct perturbation method for approximating dark soliton solutions of the nonlinear Schrödinger (NLS) equation under the influence of perturbations is presented. The problem is broken into an inner region, where the core of the soliton resides, and an outer region, which evolves independently of the soliton. It is shown that a shelf develops around the soliton that propagates with speed determined by the background intensity. Integral relations obtained from the conservation laws of the NLS equation are used to determine the properties of the shelf. The analysis is developed for both constant and slowly evolving backgrounds. A number of problems are investigated, including linear and nonlinear damping type perturbations.

Perturbation theory as applied to solitons that decay at infinity, i.e. so-called bright solitons, has been developed over many years (cf. Karpman & Maslov 1977; Kodama & Ablowitz 1981; Herman 1990). The analytical work employs a diverse set of methods including perturbations of the inverse scattering transform (IST), multi-scale perturbation analysis, perturbations of conserved quantities, etc.; the analysis applies to a wide range of physical problems. In optics, a central equation that describes the envelope of a quasi-monochromatic wave train is the nonlinear Schrödinger (NLS) equation, which in normalized form reads
where *D*,*n* are constants. In this paper, we consider the NLS equation in a typical nonlinear optics context, where *D* corresponds to the group-velocity dispersion (GVD), *n*>0 is related to the nonlinear index of refraction, *z* is the direction of propagation and *t* corresponds to the retarded time. In this form, the sign of *D* determines whether the light focuses or defocuses. In the anomalous GVD (or self-focusing nonlinearity) regime, the NLS equation exhibits so-called ‘bright’ solitons which are pulses that decay rapidly at infinity. In this case, the solitons are formed because of a balance between dispersion and self-focusing cubic nonlinearity.

On the other hand, in the normal GVD (or self-defocusing nonlinearity) regime, decaying pulses broaden and bright solitons of the NLS equation do not exist. Instead, solitons can be found as localized dips in intensity that decay off of a continuous-wave (CW) background. These dark solitons, which are termed black when the intensity of the dip goes to zero and grey otherwise, are also associated with a rapid change in phase across the pulse. The experimental observations of dark solitons in both fibre optics (Emplit *et al.* 1987) and planar waveguides (Swartzlander *et al.* 1991) sparked significant interest in the asymptotic analysis of their propagation dating back two decades.

The propagation of bright solitons under perturbation is described by the adiabatic evolution of the soliton parameters; i.e. the soliton's height, velocity, position shift and phase shift. However, the non-vanishing boundary of dark solitons introduces serious complications when applying the perturbative methods developed for bright solitons. In early work, the particular case of linear loss was studied both numerically (Zhao & Bourkoff 1989) and analytically (Giannini & Joseph 1990). The analysis was specifically for black solitons and solved explicitly for higher order correction terms. These results were re-derived (Lisak *et al.* 1991) by a more straightforward method. The method was extended to grey solitons and general perturbations but only for two of the four main soliton parameters; background height and soliton depth were determined. The evolution of the background was shown to be independent of the soliton by Kivshar & Yang (1994) and the asymptotic behaviour at infinity was used to separate the propagation of the background magnitude from the rest of the soliton. The amplitude/width of the soliton ‘core’ was determined via a perturbed Hamiltonian. Of the methods proposed, many have employed perturbation theory based on IST theory. In Konotop & Vekslerchik (1994), orthogonality conditions are derived from a set of squared Jost functions (eigenfunctions of the linearized NLS operator; Kaup 1976); from these conditions, the soliton parameters are, in principle, determined. Over the years, various corrections/modifications have been made to the details (Chen *et al.* 1998; Lashkin 2004; Ao & Yan 2005). The implementation of this method, however, has fundamental flaws, which we discuss at the end of §8.

In this paper, we address a central issue systemic through all these methods. For dark solitons, finding the adiabatic evolution of the soliton parameters (background height, soliton depth, position shift and phase shift) alone is insufficient to fully characterize the evolution of a dark soliton. We find both analytically and numerically the existence of a shelf that develops around a dark soliton under perturbation. The tendency for shelves to generate around dark solitons under external perturbation (Burtsev & Camassa 1997) was used to explain discrepancies in the perturbed conservation laws. But, the analytical calculation of the core soliton parameters was not obtained. Subsequently, shelf contribution has been ignored. However, the shelf is critical in developing the perturbation theory and has a non-trivial contribution to the integrals employed to find expressions for the soliton parameters. In this paper, we use perturbed conservation laws since they can be easily derived directly from the NLS equation and do not require the associated subtleties inherent in the IST method and extend to non-integrable problems. To carry out the procedure, we employ suitable asymptotic information about the higher order perturbation terms (beyond the soliton); other than for the ‘soliton centre’, we do not need the exact higher order solution to solve the leading-order problem for the key parameters. In order to determine the soliton centre, we find and employ the first-order perturbation solution. We note that shelves in soliton perturbation theory have been found earlier in a different, but as it turns out much easier, class of problems. They were needed to effectively understand the KdV equation under perturbation (cf. Knickerbocker & Newell 1980; Ablowitz & Segur 1981). In the KdV equation, there is a small shelf produced in the wake of the soliton. The height/speed of the soliton, shelf and the additional soliton parameter that determines the centre of the soliton are all determined by perturbation theory (Kodama & Ablowitz 1981).

In bright soliton perturbation theory, Fredholm alternatives are used to derive equations for the slowly varying soliton parameters. These conditions are identical to those derived from the perturbed conservation laws for energy and momentum (Ablowitz *et al.* 2009). A similar connection can be drawn for dark solitons. Even here, however, the existence of non-decaying eigenfunctions means that normal Fredholm alternative methods are insufficient to generate the related equations and a generalized method must be employed (Nixon 2011).

The outline of this paper is as follows. In §1, we pose the problem and illustrate how the background evolves under perturbation independent of any localized solitary wave disturbances. Sections 2–5 set up the basic analysis and a prototypical problem is discussed that helps describe the ideas. The method of multiple scales is employed to find the first-order approximation for a black soliton under the action of a dissipative perturbation that decays to zero well away from the soliton core. The concept of a moving boundary layer is used to bridge the differences between the inner soliton solution and the outer background. This discrepancy between the approximate soliton solution and the background manifests itself as a shelf developing on either side of the soliton. Perturbed conservation laws are used to find the growth of the shelf in both magnitude and phase. The analytical results are shown to be in agreement with numerical simulations of the perturbed NLS equation. In §§6–8, the method is extended to grey solitons under general perturbations. Asymptotic information about the shelf is obtained from the linear first-order perturbation equation; to determine the asymptotic states and ‘centre of phase’, the complete solution of the linear problem is not required. However, to find the soliton centre we find the first-order correction term. In §§9–11, the perturbation method is applied to some physically relevant perturbations: dissipation and two photon absorption (TPA). We find that the spatial frequency of the soliton differs from that of the background on which it resides. All of the adiabatically varying core soliton parameters and the shelf have not been obtained in the many previous studies of perturbed dark solitons.

## 1. The boundary at infinity

Let us consider the dimensionless NLS equation with normal dispersion: *D*=−1, *n*=1 (we can always rescale NLS to get these normalized values) and with an additional small forcing perturbation
1.1where |*ϵ*|≪1. Further, we will assume a non-vanishing boundary value at infinity; i.e. as . The effect the perturbation has on the behaviour of the solution at infinity is independent of any local phenomena such as pulses that do not decay at infinity; i.e. dark solitons. In the case of a CW background, which is relevant to perturbation problems with dark solitons as well as in applications to lasers, we have as , and the evolution of the background at either end is given by the equation
1.2

We write *U*^{±}(*z*)=*u*^{±}(*z*)e^{iϕ±(z)}, where *u*^{±}(*z*)>0 and *ϕ*^{±}(*z*) are both real functions of *z*. Then, the imaginary/real parts of the above equations are
1.3aand
1.3b

The above equations completely describe the adiabatic evolution of the background under the influence of the perturbation *F*[*U*]. Although this is true for all choices of perturbation, we will further restrict ourselves to perturbations that maintain the phase symmetry of equation (1.1); i.e. *F*[*U*(*z*,*t*)e^{iθ}]=*F*[*U*(*z*,*t*)]e^{iθ}. As we show next, this is a sufficient condition to keep the magnitude of the background equal on either side and a property of most commonly considered perturbations. We assume that at *z*=0,*u*^{+}(0)=*u*^{−}(0), then, since *u*^{±}(*z*) satisfy the same equation, the evolution is the same for all *z*. Hence, . While this restriction is convenient, the essentials of the method presented here apply in general. The equations for the background evolution (1.3) can now be further reduced by considering the phase difference , which is the parameter related to the depth of a dark soliton (see below); here *ϕ*^{±}(*z*) represents the phase as , respectively,
1.4Thus, while the magnitude of the background evolves adiabatically, the phase difference remains unaffected by the perturbation.

Let us now focus on the evolution of a dark soliton under perturbation. To simplify our calculations, we take out the fast evolution of the background phase , so equation (1.1) becomes
1.5The dark soliton solution to the unperturbed equation is given by
1.6where the core parameters of the soliton, *A*,*B*,*t*_{0},*σ*_{0}, are all real, the magnitude of the background is and the phase difference across the soliton is , *A*≠0. When *A*=0, equation (1.6) describes a black soliton, which has a phase difference of *π*.

## 2. The first-order correction

We write the solution in terms of the amplitude and phase: *u*=*q*e^{iϕ}, where *q* and *ϕ* are both real functions of *z* and *t*, so equation (1.5) becomes
We employ the method of multiple scales by introducing a slow scale variable *Z*=*ϵz* with the parameters *A*, *B*, *t*_{0} and *σ*_{0} being functions of *Z* and expand *q* and *ϕ* as series in *ϵ*: *q*(*Z*,*z*,*t*)=*q*_{0}+*ϵq*_{1}+*O*(*ϵ*^{2}) and *ϕ*(*Z*,*z*,*t*)=*ϕ*_{0}+*ϵϕ*_{1}+*O*(*ϵ*^{2}); we have at *O*(1)
2.1aand
2.1bwith the general dark soliton solution
2.2aand
2.2bwhere . For a black soliton, the form of solution (2.2) is taken to be
2.3aand
2.3bwhere we note that in this representation *q*_{0} is allowed to be negative. At *O*(*ϵ*), we have
2.4aand
2.4b

We begin with a somewhat simpler problem than others we deal with later; i.e. consider the linear dissipative filter perturbation
2.5and for concreteness, at leading order we assume a black pulse, equation (2.3), which satisfies the boundary conditions, equation (1.4); i.e. . This leaves us with slow evolution terms
2.6If we look for a stationary solution, *q*_{1z}=*ϕ*_{1z}=0, and note that *ϕ*_{0t}=*ϕ*_{0tt}=0, then equations (2.4) reduce to
2.7aand
2.7bwhere Im[*F*]=*γq*_{0tt} and Re[*F*]=0. Thus, the equations decouple into two linear second-order differential equations and we find exact solutions
2.8aand
2.8bwhere we have chosen the free parameters to remove exponential growth and maintain the antisymmetric property of *q* and *ϕ*, respectively. Looking at the asymptotic behaviour as , we have
2.9where the superscript ^{±} indicates the value of a function as , respectively.

## 3. Boundary layer

Notice that and as . As a result, the solution to order *ϵu*≈(*q*_{0}+*ϵq*_{1})e^{i(ϕ0+ϵϕ1)} does not match the boundary conditions at infinity. Thus, our problem is now broken into two regions: the region that matches imposed non-decaying boundary condition behaviour at infinity that is unaffected by the soliton, and the region in which the *O*(*ϵ*) correction term is valid and the solution is quasi-stationary. We introduce a boundary layer in which there is a transition from the non-zero value in the perturbation term to the boundary conditions at infinity (see also Knickerbocker & Newell 1980; Kodama & Ablowitz 1981). Note in this section we will consider the more general case when is a function of *Z*=*ϵz*. We find the behaviour of this boundary layer, where the regions are matched. For this, we return to equation (1.5) and seek a solution perturbed around the solution at infinity, say , where *w* and *θ* are real functions of *z* and *t*; the equation is automatically satisfied at *O*(1) and we have at *O*(*ϵ*)
3.1After substituting equations (1.3) and (1.4) and noting , it follows that the right-hand side is actually a higher order term and may be dropped. As a corollary, to leading order, the boundary layer is independent of perturbation. We now break equation (3.1) into real and imaginary parts
3.2Employing the compatibility conditions, *w*_{zzt}=*w*_{tzz} and *θ*_{zzt}=*θ*_{tzz}, the above equations become
3.3which is the same equation for both functions, though we will need a different solution for each. This is because of the differing boundary conditions to correctly match the inner region to the outer region. On the left side of the shelf, we match the non-zero quasi-stationary shelf to the equilibrium state at infinity (which is on the left of the soliton); the boundary conditions are
3.4On the right side of the shelf (on the right of the inner soliton), the boundary conditions are
3.5

If we let , then the ‘dispersion’ relation for equations (3.3) is found to be
3.6For long waves (*k*≪1), we approximately have or . Thus, we see that long-wave solutions (i.e. |*k*|≪1) move with instantaneous velocity .

With this in mind, we look for solutions to equations (3.3) in a moving frame of reference: , *ζ*=*z* and .
3.7

We assume that derivatives with respect to *x* are small i.e. long waves. Seeking an optimal balance in equation (3.7), we take
leaving us with
3.8

There are now two similarity solutions that we find to satisfy the boundary conditions (3.4) and (3.5) derived from matching the two regions. This is done by using the change of variable *ξ*=*x*/*ζ*^{1/3} to reduce equations (3.7) to Airy-type equations (Abramowitz & Stegun 1965). Solutions are expressed in terms of integrals of the Airy function.
3.9where *a*=−2(*V*/3)^{1/3}. Note that the sign of *a* is dependent on the sign of the velocity *V* so that the form for *w* is valid at both boundaries. For the left boundary and for the right boundary .
3.10where the lower bound on the first integral is for the left boundary and for the right boundary. For the left boundary and for the right boundary . An important point is that in the case of a black soliton, there are two boundary layers moving away from the soliton solution with speed generating a shelf. The shelf must be carefully taken into consideration when dealing with the integrals employed in soliton perturbation theory, be it in the conservation laws that we will be employing or in integral secularity conditions.

## 4. Perturbed conservation laws

We still need to solve for the slowly evolving parameters *σ*_{0}(*Z*) and *t*_{0}(*Z*) for the black soliton. This can be done by deriving equations for the growth of the shelf from the perturbed conservation laws associated with the perturbed NLS (1.5) equation from only the leading-order solution and asymptotic information about the perturbed solution. The shelf is associated with the asymptotic parameters and , which are in turn expressed in terms of *σ*_{0Z} and *t*_{0Z}. We use the Hamiltonian *H*, the energy *E*, the momentum *I* and the centre of energy *R* defined below. For grey solitons, *t*_{0} proves to be more difficult to obtain than *σ*_{0}. To find *t*_{0} we employ *u*_{1} in the expansion *u*=*u*_{0}+*ϵu*_{1}+⋯ . To find *σ*_{0}, only asymptotic information is needed.
4.1a
4.1b
4.1c
4.1dwhere *u** denotes the complex conjugate.

Note that since the standard total energy would be infinite, we define the energy of a dark pulse to be the difference of the total energy and the energy of a CW of corresponding magnitude. For the unperturbed NLS equation, the first three integrals are conserved quantities while the last can be written in terms of the momentum, i.e. d*R*/d*z*=−*I*. Evolution equations for these integrals may be easily obtained from equations (1.1) and (1.4)
4.2a
4.2b
4.2c
4.2d

## 5. The black soliton

We begin with the perturbed conservation of energy
5.1Substituting in *u*=*q*e^{iϕ}, *F*[*u*]=i*γu*_{tt}, *T*=*t*−*t*_{0} expanding *q*=*q*_{0}+*ϵq*_{1}+⋯ and taking the terms up to *O*(*ϵ*), we have
5.2aAt *O*(1) equation (5.2a) is satisfied: ; At *O*(*ϵ*), we have
5.2bThis is an equation for the change in energy caused by the propagation of the shelf. Notice that on the left-hand side of equation (5.2b), we are only integrating over , the inner region around the soliton defined by the boundary layers found in the last section. Since *q*_{0} and *q*_{1} are only functions of *T* and *Z*, we can apply the fundamental theorem of calculus to arrive at
5.2cAnd, for large *z* (although in practice, only needs to be modestly larger than the full-width half-maximum), we take and , leaving us with
5.2dBy substituting in the asymptotic approximation (2.9) found earlier for , we arrive at an expression for *σ*_{0}
5.3

Next, we consider the modified conservation of momentum
5.4Again, we let *u*=*q*e^{iϕ}, *F*[*u*]=i*γu*_{tt}, *T*=*t*−*t*_{0} and use the perturbation expansion for *u* up to *O*(*ϵ*) so that equation (5.4) becomes
5.5awhich, using *ϕ*_{0T}=0, reduces in the same way as the conservation of energy to
5.5bBy substituting in the asymptotic approximations (2.9) found earlier for , we arrive at an expression for *t*_{0}
5.6Later, we will see that the above result agrees with the more general grey soliton case.

This can now be compared with direct numerical simulations. The magnitude and phase are depicted in figure 1*a*,*b*, respectively, for *z*=30 and . Here, we see the inner region discussed earlier is *t*∈(−30,30), where the asymptotic solution agrees well with the numerical solution. The remainder of the domain constitutes the outer region where the inner asymptotic solution matches with the exterior rest state; this is discussed subsequently. The boundary layer shown in figure 2*b* compares the solutions (3.9) and (3.10) to numerics and illustrates how the inner and outer solutions are connected. The propagation of this boundary layer can be seen in figure 2*a* where the contour plot illustrates the soliton down the middle and the shelf extending out from it. The speed of the boundary layer matches the speed predicted by the long-wave approximation in §3.

## 6. The grey soliton

Now we consider the general case of a grey soliton with velocity *A*(*Z*); we also recall . Let *u*=*q*e^{iϕ}, where *q*>0 and *ϕ* are real functions of *z* and *t*, and introduce the moving frame of reference and *ζ*=*z*, so that with *u*=*u*(*ζ*,*T*,*Z*), equation (1.5) becomes
6.1Then using *u*=*q*e^{iϕ}, this is now broken into imaginary and real parts, respectively,
6.2aand
6.2bWe now write equation (6.2b) in terms of the slow evolution variable *ζ*=*ϵZ* and series expansions *q*=*q*_{0}+*ϵq*_{1}+*O*(*ϵ*^{2}) and *ϕ*=*ϕ*_{0}+*ϵϕ*_{1}+*O*(*ϵ*^{2}). At *O*(1), the equations are satisfied by the soliton solution (2.2). We look for stationary solutions at *O*(*ϵ*)
6.3aand
6.3bwhere *u*_{0}=*q*_{0}e^{iϕ0} and
6.4aand
6.4b

Next we assume a shelf structure similar to the one found in our example problem will develop; this is supported by numerical computations. Consider that equation (6.3a) in the limit using and yields
6.5We assume *q*_{1} tends to a constant with respect to *t*; i.e. as . As a result . Then, *q*_{1} and *ϕ*_{1T} both tend asymptotically to constants as , which corresponds to a shelf developing around the soliton. Substituting *ϕ*_{0T} into equation (6.3b), in the limit , we get
6.6We define as the phase change across the core soliton. This is consistent with the soliton parameters *A* and *B* being expressed in terms of background magnitude, , and phase change, Δ*ϕ*_{0},
6.7Using from §1 and substituting equation (6.7) into equation (6.6) we find
6.8

## 7. Conservation laws for the grey soliton

Next we use the modified conservation equations (4.2d) to solve for the shelf parameters and as well as the slow evolution variables *A*, *σ*_{0Z}. More work is required in order to find *t*_{0}. Note that if we find *A*, then . The edge of the shelf still propagates with velocity ; however, the speed may now vary in *z*. In terms of the moving frame of reference, the boundaries of the shelf are
7.1where *S*_{L} and *S*_{R} give the position in *T* of the left and right boundaries of the shelf, respectively, at *ζ*. Note that for all *Z*; thus, the soliton cannot overtake the shelf. In figure 3, we illustrate the general structure of a perturbed dark soliton with the moving shelf. The inner region consists of the core soliton and the shelf expanding around it, while the outer region consists of the infinite boundary conditions characterized by equations (1.4). The boundaries between these regions are delineated by dotted red lines at *t*=*S*_{L} and *t*=*S*_{R}.

We begin with the evolution equation for the Hamiltonian (4.2):
7.2Substituting in *u*=(*q*_{0}+*ϵq*_{1})e^{i(ϕ0+ϵϕ1)} and changing variables to the moving frame of reference, we have up to *O*(*ϵ*)
7.3where both here and later on *u*_{0}=*q*_{0}e^{iϕ0}. The Hamiltonian is unique among the evolution equations (4.2d) in that the contribution of the shelf appears only at *O*(*ϵ*^{2}) or higher and that to *O*(*ϵ*) may be ignored. We now put in the soliton form (2.2) to get
7.4Taking a derivative with respect to *Z* of the equation , we get , which can be used to consolidate equation (7.4) down to
7.5

The evolution equations for energy (5.1) and momentum (5.4) both remain the same after transforming to the moving frame of reference
7.6and
7.7The inner region over which *q*_{1} and *ϕ*_{1} are relevant is *T*∈[*S*_{L}(*ζ*),*S*_{R}(*ζ*)], and outside this region *q*_{1}=*ϕ*_{1T}=0. At *O*(1), the equations are satisfied and at *O*(*ϵ*) we have
7.8aand
7.8bSince the integrands on the left-hand side are not functions of *ζ*, we can apply the fundamental theorem of calculus to arrive at
7.9aand
7.9bWe are left now with the evolution of the centre of energy
7.10which after transforming to the moving frame of reference is now
After rearranging some terms, we have
7.11a
7.11b
7.11c
7.11dEquation (7.11) line (*a*) yields
7.12The terms on equation (7.11) line (*b*) reproduce the energy equation (7.6) and cancel out. The terms on equation (7.11) line (*c*) are calculated up to *O*(*ϵ*) using the previous results by integrating the energy and momentum equation (7.9b)
Note that d/d*ζ*=*ϵ*d/d*Z* and *S*_{R} and *S*_{L} are *O*(1/*ϵ*) in terms of *Z*.

Putting everything together in terms of *Z*=*ϵζ* yields
7.13After some cancellations, this breaks into *O*(1) terms
7.14and *O*(*ϵ*) terms that include *t*_{0Z} and higher order energy and momentum terms, which have not been determined. The six equations, equations (6.8), (7.5), (7.9a), (7.9b) and (7.14), can now be used to solve for the set of six parameters , , *A* and *σ*_{0}. To find *t*_{0Z}, we need to employ more information—see the next section.
7.15a
7.15b
7.15c
7.15d
7.15e
7.15f
7.15g
7.15h
7.15i

We have added equations (7.15h) and (7.15i) to the list since it is often easier to use these formulations for *B*_{Z} and Δ*ϕ*_{0Z} rather than working out *B* and Δ*ϕ*_{0} explicitly and then taking derivatives.

By combining equations (7.5) and (7.9b), we arrive at
7.16which may be rewritten as
7.17If we define Δ*ϕ*_{1} as
7.18then *ϵϕ*_{1} is the phase change across the shelf. Substituting this definition along with (7.15i) into equation (7.17), we arrive at
7.19Thus, the total phase change across the inner region remains constant, which is consistent with our earlier result that (the phase change from to ) remains constant for all perturbations. As an example, figure 8 shows that the entire phase change remains consistent with the given boundary condition.

## 8. *t*_{0Z} and higher order terms

To find the final parameter *t*_{0}, we employ the first-order correction term. We look for a series solution to equation (1.1) of the form *u*=*u*_{0}+*ϵu*_{1}+*O*(*ϵ*^{2}), and at *O*(*ϵ*) we have
8.1After changing variables to the moving frame of reference , *z*=*ζ* we have
8.2Here,
8.3

If we look for stationary solutions (∂/∂*ζ*=0), this can be written as a system of coupled second-order differential equations
8.4awhere
8.4band
8.4c

In the limit , this system decouples into two second-order differential equations that are not difficult to solve and give two strictly real solutions and two strictly imaginary solutions. For each solution we found for *A*=0, we assume there exists a solution for *A*≠0 that differs only in the perpendicular direction; e.g. if satisfies equation (8.4c) with *A*=0, then there exists *u*_{I} such that satisfies equation (8.4c) with *A*≠0; namely only the second component changes and hence the system reduces to the first-order equation that is consistent with the remaining equations. Under this assumption, we find a complete set of homogeneous solutions
8.5a
8.5b
8.5cand using variation of parameters, we can obtain a particular solution, **U**_{1p}, for the forcing *G*[*u*_{0}].

After combining real and imaginary parts, the full solution to equation (8.2) is given by
8.6where *c*_{1}, *c*_{2}, *c*_{3} and *c*_{4} are functions of *Z* and *u*_{1p} is dependent on the yet to be determined *t*_{0}. We take *c*_{4}=0 to remove the exponential growth in *u*_{14} and separate out the contribution of *t*_{0Z} that appears linearly in the particular solution *u*_{1p}
8.7where
8.8so that has no unknowns left in it.

To put *u*_{1} in terms of the magnitude and phase functions *q*_{0}, *q*_{1}, *ϕ*_{0} and *ϕ*_{1}, we expand our previous approximation for *u*
8.9so that
8.10
8.11
8.12

By looking at the asymptotic behaviour of the solution *u*_{1} as , we find the equation
8.13Since *u*_{11T}, *u*_{12T} and *u*_{1pT} all go to zero in the limit , the above equation can be used to find *c*_{3}.

With this, we are now able to find a second-order differential equation for *t*_{0} from the Hamiltonian at *O*(*ϵ*^{2})
8.14awhere
8.15

On the left-hand side, we have the slow evolution of the *O*(*ϵ*) terms and the fast evolution of the *O*(*ϵ*^{2}) terms. *H*_{1} is dependent on *u*_{0} and *u*_{1} and is given by
8.16*H*_{2} is dependent on *u*_{0}, *u*_{1} and the order *ϵ*^{2} correction *u*_{2}. However, as before, we assume a stationary (in the moving frame of reference) solution *u*_{2} (as was done for *u*_{0} and *u*_{1}), then the derivative of *H*_{2} with respect to the fast evolution variable *ζ* only depends on the asymptotic behaviour of *u*_{0} and *u*_{1} and is given by
8.17Though it is not immediately obvious, we find that *c*_{1} and *c*_{2} do not contribute to the Hamiltonian in equation (8.17), so *t*_{0} is the only unknown. We take *t*_{0}(0)=0, and to find a suitable initial condition *t*_{0Z}(0), we require that the Hamiltonian be accounted for by *H*_{0} at *z*=0; i.e. the higher order terms are initially zero
8.18

Our prediction for *t*_{0} differs greatly from that given by methods based on a ‘so-called’ complete set of squared Jost function. This discrepancy may be partially explained by the assumption that the squared Jost function forms a basis for the solution space of equation (8.2). The eigenfunctions are found in the inverse scattering theory for the defocusing NLS equation with non-vanishing boundary values (Zakharov & Shabat 1973), and as a direct result, the acquired basis functions associated with the soliton are localized and bounded. However, we have solved explicitly for the first correction term, and we find that the solution has an expanding shelf. From this, we deduce that the squared Jost functions associated with the soliton are an insufficient basis.

## 9. Dissipative perturbation

Let us return to the perturbation *F*[*u*]=i*γu*_{tt}, however, we now consider the evolution of a general dark soliton with . As was the case for black solitons, the background height is found to be constant from equations (1.4); i.e. . In figure 4*a*, we see that the velocity of the soliton does not affect the velocity of the shelf, which still moves with velocity . Using the equations derived in §§7 and 8, we can now solve for all relevant parameters
9.1a
9.1b
9.1cWe also note that these results agree with the black soliton when *A*=0. In the limit , we have ; however, this has a discontinuity in its derivative, so instead we used for our black soliton calculations. As a result, there is a sign difference in from equation (2.9).

Unlike the speed of the shelf at the edges, the magnitude and phase, , depend on the soliton's velocity (which is in turn related to the soliton's depth and the phase across the soliton). As illustrated in figure 4*b*, the shelf is shallower behind the soliton for larger speeds (or smaller phase change Δ*ϕ*_{0}). The extra phase induced by the perturbation means that the spatial frequency of the soliton is different from the frequency of the CW background on which it lies. Though *σ*_{0} evolves adiabatically, the soliton eventually becomes noticeably out of phase from the background as shown in figure 5*a*. Here, the background phase (*ϕ*^{+} and *ϕ*^{−}) is constant since the fast evolution of the background phase was taken out in equation (1.1). Finally, in figure 5*b*, we show the improvement finding *t*_{0}(*Z*) makes on predicting the centre of the soliton over just using the velocity *A*.

## 10. Linear damping

We now apply our results to the case of linear damping 10.1which was both the first (Giannini & Joseph 1990) and a commonly used example used in the study of perturbed dark solitons.

In this example, we now have a moving background found by solving equation (1.4) 10.2

From equations (7.15), we can determine the slowly varying soliton parameters 10.3a 10.3b 10.3c

Figure 6*a* shows that the existence of a raised shelf and dynamics of the shelf edge are well predicted by the asymptotic theory. The background height and trough of the soliton (*A*(*Z*)) are accurately approximated by our method (figure 7*a*); this agrees with previously found approximations (Kivshar & Yang 1994). Our results for *t*_{0} and *σ*_{0} are plotted in figure 7*b*. As mentioned earlier, previous attempts using IST are not adequate (Chen *et al.* 1998; Lashkin 2004; Ao & Yan 2005). *t*_{0} was not obtained in Kivshar & Yang (1994); furthermore, no previous work has considered the evolution of the parameter *σ*_{0}.

## 11. Two-photon absorption

Dark solitons have been proposed in the development of optical switching devices (Luther-Davies *et al.* 1992). Here, materials with high nonlinearities are used to reduce the power for soliton formation and the switching threshold; however, an enhanced TPA coefficient also accompanies these materials. An example of TPA with strong defocusing nonlinearity is the semiconductor ZnSe (Skinner *et al.* 1991). This is represented by the perturbation term
11.1

From equations (7.15), we can find all parameters other than *t*_{0}, which are given below. This yields the evolution of the background height and trough height (figure 8*a*); we also find that the phase change across the core soliton does not remain constant (as had been the case for our previous examples).
11.2a
11.2b
11.2c
11.2d

As in the linear damping example, a raised shelf develops around the core soliton as seen in figure 6*b*. Again we see remarkable correlation with our predicted shelf velocity. Figure 8*a* shows strong agreement between numerics and asymptotic analysis for the evolution of both and *A*(*Z*). Figure 8 shows that while the phase change across the core soliton decreases by half, this is compensated for by the change in phase over the shelf, and the total phase change from negative infinity to positive infinity remains constant, which agrees with equation (7.19) (figure 8*b*).

## 12. Conclusion

In conclusion, we develop a novel approach to dark soliton perturbation theory that breaks the problem into an inner region around the soliton and a shelf that matches to the boundary at infinity. Under perturbation, a dark soliton develops a shelf the edge of which propagates out at a speed equal to the magnitude of the CW background. In analytical terms, the shelf arises owing to properties of the perturbation, which serve to drive mean contributions in the amplitude and growing terms in the phase. It is also found that the soliton can have a different frequency from the CW background. The method can be applied to general perturbations that can have both moving and constant backgrounds. For typical perturbations, the asymptotic approximations were calculated and were compared favourably with direct numerical results. These comparisons confirmed the existence of the analytically predicted shelf and support our results. The non-vanishing background and soliton are treated separately from the core soliton; in this way, we obtain a consistent perturbation theory for dark solitons.

## Acknowledgements

This research was partially supported by the US Air Force Office of Scientific Research, under grant FA9550-09-1-0250 and the NSF under grant DMS-0905779.

- Received December 20, 2010.
- Accepted March 7, 2011.

- This journal is © 2011 The Royal Society