1. Introduction
The mathematical modelling of physical phenomena often leads to a certain class of nonlinear partial differential equations (PDEs) known as integrable systems. Distinguished features of integrable systems are that: (i) they admit soliton solutions, and (ii) their initial-value problem can be effectively linearised via the Inverse Scattering Transform (IST). One of the prototypical integrable systems is the nonlinear Schrödinger (NLS) equation:
\begin{equation*}
iq_t+q_{xx}-2\sigma |q|^2q=0, \qquad \sigma=\pm 1,
\end{equation*}where 
$q(x,t)$ is a complex function of 
$x,t\in \mathbb{R}$, subscripts denote partial differentiation throughout, and 
$\sigma$ distinguishes between anomalous and normal dispersion, with 
$\sigma=-1$ corresponding to the ‘focusing’ NLS, and 
$\sigma=1$ to the ‘defocusing’ NLS. The NLS equations appear as universal models for weakly dispersive nonlinear wave trains, and have been derived in such diverse fields as deep water waves, plasma physics, nonlinear optical fibres, low-temperature physics and Bose-Einstein condensates, magneto-static spin waves and more [Reference Benney and Roskes6, Reference Hasegawa and Tappert22, Reference Hasegawa and Tappert23, Reference Kalinikos, Kovshikov and Patton26, Reference Pethick and Smith45, Reference Zakharov55, Reference Zakharov56, Reference Zvezdin and Popkov58].
Unlike its focusing counterpart, the defocusing NLS does not admit localised bright solitons, exponentially decaying as 
$|x|\to \infty$, but it possesses dark soliton solutions, which appear as localised dips of intensity over a nonzero background [Reference Kevrekidis, Frantzeskakis and Carretero-Gonzàlez33, Reference Kivshar and Luther-Davies34]. We consider the defocusing NLS equation in the form:
\begin{equation}
iq_{t}+q_{xx}-2(|q|^{2}-q_{0}^{2})q=0,
\end{equation}with the constant nonzero boundary conditions
\begin{equation}
q(x,t)\rightarrow q_{\pm}\equiv q_{0}e^{i\theta_\pm} \qquad \text{as } x\rightarrow\pm\infty,
\end{equation}where 
$q_{0} \gt 0$ is the (symmetric) amplitude of the background and 
$\theta_\pm\in \mathbb{R}$ are the asymptotic phases; the additional linear term proportional to 
$q$ in (1.1) is introduced via a gauge transformation 
$q(x,t)\mapsto e^{-2iq_0^2t}q(x,t)$ in order to allow for time-independent boundary conditions. The exact single dark soliton solution of (1.1) is given by:
\begin{equation}
q(x,t)=
e^{i\sigma_1}\left\{k_{1}+i\Lambda_{1}\tanh\left[\Lambda_{1}(x+2k_{1}t-x_{1})\right]\right\},
\end{equation}where 
$-q_0 \lt k_1 \lt q_0$ determines the soliton velocity, while 
$\Lambda_1$ fixes its amplitude (i.e., the depth of the soliton below the background 
$q_0$) and 
$x_1,\sigma_1\in \mathbb{R}$ give the soliton centre and phase, respectively. Note that the dark soliton amplitude and velocity are not independent, since they are related to each other and to the background amplitude 
$q_0$ via:
\begin{equation}
k_1^2+\Lambda_1^2=q_0^2.
\end{equation} In addition, the phases 
$\theta_\pm$ of the background are related to the soliton parameters via:
\begin{equation}
\theta_\pm =\sigma_1\pm \textrm{Arg}(k_1+i|\Lambda_1|).\end{equation}One limitation of integrable models is that, in general, the systems considered in physical experiments are non-integrable. On the other hand, the theoretical predictions for the soliton solutions in integrable cases provide an extremely valuable tool for investigating non-integrable solitary waves in regimes that are not too far from the integrable ones. As such, researchers rely on perturbation-based techniques of related integrable systems, when possible, to study how solitons and their evolution are affected by the inclusion in the mathematical description of terms that account for small dissipation, small linear/nonlinear loss, etc.
The perturbation theory for solitons that decay rapidly at infinity has been extensively explored since the late 1970s – using various different methods such as multi-scale perturbation analysis, IST-based techniques, eigenfunction-expansion techniques, perturbations of conserved quantities and direct numerical simulations [Reference Elgin15, Reference Herman24, Reference Karpman and Maslov27, Reference Kaup29–Reference Kaup and Newell31, Reference Kivshar and Malomed35, Reference Kodama and Ablowitz38, Reference Yang50, Reference Yang51]. On the other hand, the nonvanishing background characteristic of dark solitons introduces significant difficulties when one attempts to apply these perturbative approaches developed for the rapidly decaying case.
As mentioned above, for the scalar defocusing NLS equation, dark solitons are completely determined by the four parameters: 
$q_0$, 
$\Lambda_1$ (or 
$k_1)$, 
$x_1$, 
$\sigma_1$, which in general under perturbation develop an adiabatic evolution (e.g., they acquire a non-trivial dependence on a slow time variable 
$T=\varepsilon t$). Some early works investigated the perturbation of black (i.e., stationary dark) solitons in lossy fibres numerically [Reference Zhao and Bourkoff57] and later on analytically [Reference Giannini and Joseph20, Reference Lisak, Anderson and Malomed41]. In [Reference Giannini and Joseph20], a direct perturbation theory was developed based on a formal series expansion of the solution where the 
$x$ and 
$\varepsilon t$ dependencies are separated out, and the first two terms in the expansion are computed by direct integration. The method developed in [Reference Lisak, Anderson and Malomed41], based on perturbed conserved quantities, was subsequently extended to grey (i.e., non-stationary dark) solitons and to generic perturbations, but only two of the four main soliton parameters, 
$q_0$ and 
$\Lambda_1$, were determined. In [Reference Kivshar and Yang36], it was demonstrated that, under perturbation, the background evolves independently of the soliton. By separating the background amplitude from the soliton ‘core’, the authors were able to determine the soliton’s amplitude and width through a Hamiltonian method based on perturbed conservation laws. Other works investigated instability-induced dynamics of dark solitons and oscillations of dark solitons in trapped Bose-Einstein condensates [Reference Pelinovsky, Frantzeskakis and Kevrekidis43, Reference Pelinovsky, Kivshar and Afanasjev44, Reference Smirnov, Smirnova, Ostrovskaya and Kivshar48]. It is important to note that, for dark solitons, the adiabatic evolution of the soliton parameters alone does not fully characterise the perturbed solution. This is because the perturbation produces a moving shelf on both sides of the soliton. The presence of this shelf – confirmed both numerically and analytically – was in fact used in [Reference Burtsev and Camassa8] to account for discrepancies observed in the perturbed conservation laws, although the soliton core parameters were not determined analytically. Shelves emerging from the wake of a bright soliton had already been observed in soliton perturbation theories for the Korteweg-de Vries (KdV) equation [Reference Ablowitz and Segur2, Reference Knickerbocker and Newell37], the fifth-order KdV equation [Reference Yang52] and the complex modified KdV [Reference Yang53, Reference Yang54]. For this class of problems, it was possible to determine the bright soliton parameters and the shelf parameters by either direct perturbation theory [Reference Kodama and Ablowitz38] or by using the squared eigenfunctions, as detailed in [Reference Yang54].
To date, the most comprehensive analysis of dark-soliton perturbations for the scalar defocusing NLS is that of Ablowitz et al. in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1]. Using a multi-scale expansion together with perturbed conservation laws, the authors derived both the magnitude and phase of the shelf, as well as the adiabatic evolution of all soliton parameters, showing the emergence of a moving boundary layer that links the inner soliton core to the outer background. This approach was recently generalised to describe the effect of small perturbation to the dark-bright solitons of the coupled NLS equation, the so-called Manakov system [Reference Chernyavskiy, Frantzeskakis, Horikis, Koutsokostas and Prinari11] (see also [Reference Rothos47], where the perturbed evolution of the soliton parameters and the background were determined from variations of the Riemann-Hilbert problem (RHP)).
An alternative approach to soliton perturbation theory in the rapidly decaying case was introduced in [Reference Elgin15, Reference Kaup29, Reference Kaup and Newell31]. In this method, based on the IST, one calculates the variations of the scattering data of the soliton due to the perturbation, and then uses inverse scattering to reconstruct the perturbed solution. In the process, squared eigenfunctions (i.e., quadratic combinations of Jost eigenfunctions and their adjoints) play a critical role. Another soliton perturbation theory that solves the first-order perturbation equation directly was also developed in [Reference Fogel, Trullinger, Bishop and Krumhansl17, Reference Kaup30, Reference Keener and Mclaughlin32] (see the monograph [Reference Yang54] and references therein). In this method, the first-order perturbation equation is solved by expanding its solution into a set of complete eigenfunctions of the linearisation operator. This method does not explicitly use the IST, and it is often easier to apply. But its connection to the IST is still critical, since the eigenfunctions of the linearisation operator are simply the squared eigenfunctions in the IST. We will refer to this latter perturbative approach interchangeably as ‘integrable’ perturbation theory or eigenfunction-expansion-based perturbation theory.
Since the early 1990s, numerous efforts have been made to extend the integrable perturbation theory to dark solitons. In 1994, Konotop and Vekslerchik derived orthogonality conditions from a set of squared eigenfunctions for the scalar defocusing NLS equation on a constant background, from which all soliton parameters can, in principle, be obtained [Reference Konotop and Vekslerchik39]. However, this early work did not take into account the background evolution induced by the perturbation. Subsequent attempts at rigorous proofs of the completeness of the squared eigenfunctions followed: in [Reference Chen, Chen and Huang10] Chen, Chen & Huang showed completeness of the squared eigenfunctions in the case of one-soliton, and used it to investigate analytically a linear damping perturbation, and in [Reference Chen and Chen9] applied the method again to a self-steepening perturbation. In [Reference Huang, Chi and Chen25], Huang, Chi & Chen used a generalised Marchenko equation and extended the completeness proof of the previous work [Reference Chen, Chen and Huang10] to the multi-soliton case. The proof in [Reference Chen, Chen and Huang10] was then claimed to be incorrect by Ao & Yan in [Reference Ao3, Reference Ao and Yan4], based on the observation that the complete set should have two, not just one, continuous spectrum basis vector, which resulted in different predictions for the soliton velocity and the first-order correction. In a subsequent work by Ao & Yan [Reference Ao and Yan5], the results of [Reference Chen, Chen and Huang10, Reference Huang, Chi and Chen25] and [Reference Ao3, Reference Ao and Yan4] were then declared to be ‘equivalent’ under an appropriate ‘transformation between two integral variables’. As a matter of fact, the difference between the two can be traced to the fact that the earlier work used one of the symmetries of the scattering problem to reduce the number of eigenfunctions involved in the closure relation to only the ones that are linearly independent. It is also worth mentioning that in [Reference Lashkin40] squared eigenfunctions were used (though without explicitly referring to them, or to their completeness) to develop an eigenfunction-expansion-based perturbation theory for the defocusing NLS on a background, with no reference to the work in [Reference Chen, Chen and Huang10, Reference Huang, Chi and Chen25] (conversely, the papers [Reference Ao3–Reference Ao and Yan5] do not reference [Reference Lashkin40]). An integrable perturbation theory for the dark-bright solitons of the defocusing Manakov equation was also developed in [Reference Mylonas, Rothos, Kevrekidis and Frantzeskakis42].
Importantly, in all the earlier implementations of the dark-soliton perturbation theory that is based on eigenfunction expansions [Reference Ao3–Reference Ao and Yan5, Reference Chen and Chen9, Reference Chen, Chen and Huang10, Reference Huang, Chi and Chen25, Reference Konotop and Vekslerchik39, Reference Lashkin40, Reference Mylonas, Rothos, Kevrekidis and Frantzeskakis42], some of the orthogonality conditions were flawed, because the discrete eigenmodes that were used did not account for the contributions from the poles at the branch points of the integral term in the completeness relation. As we will discuss in detail in Section 5 and Section 7, this resulted in erroneous evolution equations for (at least) some of the soliton parameters. Moreover, none of the above references attempted to determine the radiation shelf that develops around the dark soliton, or presented comparisons of the theoretical predictions with numerical simulations. The existence of the shelf has long been confirmed by numerical simulations, and, as was shown in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], it is critical in developing the perturbation theory and contributes to the integrals used to determine the evolution of the soliton parameters. For instance, the adiabatic evolution of the soliton centre in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] is markedly different from the one obtained by expanding the solution in terms of squared eigenfunctions in [Reference Ao3–Reference Ao and Yan5, Reference Chen and Chen9, Reference Chen, Chen and Huang10, Reference Huang, Chi and Chen25, Reference Konotop and Vekslerchik39, Reference Lashkin40, Reference Rothos47], and it was conjectured in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] that, because of the existence of the expanding shelf, the squared eigenfunctions associated with the soliton are an insufficient basis, questioning in fact the existence of a closure relation for this problem.
The goal of this work is to revisit the perturbation theory of the scalar defocusing NLS on a nontrivial background based on the squared eigenfunctions and develop it so that it can correctly predict the slow-time evolution of all the dark soliton parameters, as well as the radiation shelves emerging on the sides of the soliton. First, we prove completeness of the squared eigenfunctions for general potentials. A crucial difference with respect to the earlier works is the need to properly account for the singularities of the scattering data at the points 
$\pm q_0$, which are branch points for the continuous spectrum. Taking contributions from such singularities out of the integral term of the closure relation and leaving the remaining integral as a principal-value integral proves to be critical in this eigenfunction-expansion-based perturbation theory, because this leads to the correct discrete eigenmodes to be used in orthogonality conditions of the perturbation theory. Then we use the 1-soliton closure relation and suitable suppression of secular growth to determine the adiabatic evolution of the soliton parameters 
$k_1(T),\Lambda_1(T)$, as well as a condition relating the evolution of 
$\sigma_1(T)$ and 
$x_1(T)$ under perturbations of order 
$\varepsilon$ as functions of a slow time variable 
$T=\varepsilon t$. The latter condition is missing in all the previous works that used eigenfunction-expansion-based perturbation theories. More importantly, the first-order correction integral to the dark soliton computed via squared eigenfunctions is shown to contain a pole due to singularities of the scattering data at the branch points. Analysis of this first-order correction integral leads to predictions for the height and velocity of the shelves. Moreover, this same analysis provides the spatial phase gradient on each side of the shelf, as well as a formula for the slow-time evolution of the core soliton’s phase, which in turn allows one to determine the dependence of the soliton centre on the slow time. All the results are corroborated by direct numerical simulations, and compared with the results of the direct perturbation theory in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], and with the earlier works using integrable perturbation theory. In particular, the estimates we obtain for all soliton parameters agree with the ones in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] to order 
$\varepsilon$, the only difference being for the soliton centre, which, unlike the other parameters, obtained from perturbed conserved quantities up to 
$\mathcal{O}(\varepsilon)$, is determined in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] in terms of differential equations obtained from the Hamiltonian at 
$\mathcal{O}(\varepsilon^2)$. We will also explain the shortcomings of the predictions for the adiabatic evolution of the soliton parameters in the earlier works within the framework of the eigenfunction-expansion-based perturbation theory.
The plan of the paper is the following. In Section 2 we give an overview of the IST in order to set the notations and present the properties of eigenfunctions and scattering data that are necessary for the following sections. In Section 3 we derive the closure relation of squared eigenfunctions for general potentials by calculating the variations in the IST generalised to the case of a nonzero background. In Section 4 we discuss the linearisation operator and obtain the 1-soliton closure relation. The multiple scale perturbation theory is developed in Section 5. Section 6 investigates the first-order correction, the shelf and evolution of the phase. The comparison with direct numerical simulations and earlier results is provided in Section 7, with explicit applications to linear and nonlinear damping, dissipation and self-steepening. Finally, Section 8 is devoted to some concluding remarks and open problems, and more technical calculations are collected in the appendices.
2. Overview of the IST
In this section, we give a succinct overview of the IST for the defocusing NLS with nonzero boundary conditions, to set the notations and review the properties of eigenfunctions and scattering data, which are required for the following sections. Additional details can be found, for instance, in [Reference Faddeev and Takhtajan16, Reference Prinari46].
We start by considering the (unperturbed) defocusing NLS in the form (1.1) with the constant symmetric nonzero boundary conditions
\begin{equation}
q(x,t)\rightarrow q_{\pm}\equiv q_{0}e^{\pm i\theta}\qquad \text{as }x\rightarrow\pm\infty, \qquad q_0 \gt 0,\ \theta\in \mathbb{R}.
\end{equation} The asymptotic phases can be chosen as 
$\pm \theta$ without loss of generality on account of the invariance of the PDE upon multiplication by an arbitrary phase. [Note that arbitrary asymptotic phases 
$\theta_\pm$ as 
$x\to \pm \infty$ can be simply accounted for by replacing 
$\theta=(\theta_+-\theta_-)/2$ throughout.] As shown in [Reference Cuccagna and Jenkins13, Reference Demontis, Prinari, van der Mee and Vitale14, Reference Gallo18, Reference Gallo19], the IST can be formulated as a bijection for an initial condition such that 
$q(x,0)-\tilde{q}(x)\in H^{1,1}(\mathbb{R})$, where
\begin{equation*}
\tilde{q}(x)=q_o\left[ \cos \theta +i \sin \theta \tanh x \right],
\end{equation*}
\begin{equation*}
H^\ell(\mathbb{R}):=\left\{
f:\mathbb{R} \to \mathbb{C} \ \text{s.t.}\
f^{(k)}\in L^2(\mathbb{R}), \ k=0,\dots,\ell
\right\}, \quad
H^{1,\ell}(\mathbb{R}):=L^{2,\ell}(\mathbb{R})\cap H^\ell(\mathbb{R}),
\end{equation*} 
$L^p(\mathbb{R})$ are the standard Lebesgue spaces, and the weighted spaces 
$L^{p,s}(\mathbb{R})$ have norms defined as
\begin{equation*}
||f||_{L^{p,s}(\mathbb{R})}:=\left(
\int_\mathbb{R} \langle x \rangle^{2s}|f(x)|^p\, dx \right)^{1/p},
\qquad
\langle x\rangle:=\sqrt{1+x^2}.
\end{equation*} Additional smoothness and decay of the reflection coefficient are established in [Reference Cuccagna and Jenkins13] under stronger assumptions on the initial condition (e.g., by considering 
$H^{1,\ell}(\mathbb{R})$ for 
$\ell=3/2,2$, see also [Reference Gkogkou, Prinari and Trogdon21] for a review).
The Lax pair of the defocusing NLS is given by
\begin{equation}
v_{x}=Uv,\qquad U=-ik\sigma_{3}+Q,\qquad \sigma_{3}=\begin{bmatrix}1&0\\0&-1\end{bmatrix},\qquad Q=\begin{bmatrix}0&q\\q^{*}&0\end{bmatrix},
\end{equation}
\begin{equation}
v_{t}=Vv,\qquad
V=i\sigma_{3}(-2k^{2}+Q_{x}-Q^{2}+q_{0}^{2})+2kQ.
\end{equation}
$\lambda$ of (2.2a) as 
$x\to \pm \infty$ are related to the spectral parameter 
$k$ by 
$\lambda^{2}=k^{2}-q_{0}^{2}$. Thus, 
$\lambda=\lambda(k)$ is a multivalued function with branch points at 
$k=\pm q_{0}$, and a two-sheeted Riemann surface cut along 
$(-\infty, -q_0)\cup (q_0,+\infty)$ can be defined such that 
$\mathop{\rm Im}\nolimits \lambda\ge 0$ on one sheet, and 
$\mathop{\rm Im}\nolimits \lambda\le 0$ on the other. A uniformisation variable 
$z$ can be introduced in the following way
\begin{equation}
z=k+\lambda\,, \qquad k=\frac{1}{2}(z+q_{0}^{2}z^{-1})\,, \qquad \lambda=\frac{1}{2}(z-q_{0}^{2}z^{-1}),
\end{equation}such that both sheets of the Riemann surface on which 
$k$ is defined are mapped into a single complex 
$z$-plane. The matrix Jost solutions that satisfy both parts of the Lax pair (2.2) are defined by the boundary conditions
\begin{equation}
\Phi(x,t,z)=\begin{bmatrix}
\phi(x,t,z)&\bar\phi(x,t,z)
\end{bmatrix} \sim E_{-}(z)e^{-i\Omega(x,t,z)\sigma_{3}} \qquad \text{as } x\rightarrow-\infty,
\end{equation}
\begin{equation}
\Psi(x,t,z)=\begin{bmatrix}
\bar\psi(x,t,z)&\psi(x,t,z)
\end{bmatrix} \sim E_{+}(z)e^{-i\Omega(x,t,z)\sigma_{3}}\qquad \text{as } x\rightarrow+\infty,
\end{equation}
\begin{equation}
E_{\pm}(z)=I-iz^{-1}\sigma_{3}Q_{\pm}\equiv \begin{bmatrix} 1 & -iq_\pm/z \\ iq^{*}_\pm/z & 1 \end{bmatrix},
\quad
Q_\pm =\begin{bmatrix}0&q_\pm \\q^{*}_\pm &0\end{bmatrix},
\end{equation}
\begin{equation}
\Omega(x,t,z)=\lambda(z)(x+2k(z)t).
\end{equation}
$I$ denotes the 
$2\times 2$ identity matrix. It can be shown that if 
$q\to q_\pm$ sufficiently rapidly as 
$x\to \pm \infty$, the eigenfunctions 
$\phi$ and 
$\psi$ are analytic for 
$\mathop{\rm Im}\nolimits z \gt 0$ and continuous for 
$\mathop{\rm Im}\nolimits z\ge 0$, while 
$\bar\phi$ and 
$\bar\psi$ are analytic for 
$\mathop{\rm Im}\nolimits z \lt 0$ and continuous for 
$\mathop{\rm Im}\nolimits z\le 0$, in both cases including at the images of the branch points 
$z=\pm q_0$. Furthermore, 
$\Phi$ and 
$\Psi$ have a singularity at 
$z=0$ (cf. (2.5a)), and
\begin{equation}
\det\Phi(x,t,z)=\det\Psi(x,t,z)=\gamma(z)\,,\qquad \gamma(z)=1-q_{0}^{2}z^{-2},
\end{equation}which is nonzero except at 
$z=\pm q_{0}$. For 
$z\in\mathbb{R}\backslash\{\pm q_{0}\}$, the Jost solutions can be related via:
\begin{equation}
\Phi(x,t,z)=\Psi(x,t,z)S(z)\,,\qquad S(z)=\begin{bmatrix}
a(z)&\bar{b}(z)\\b(z)&\bar{a}(z)
\end{bmatrix}\,,\qquad \det S(z)=1.
\end{equation} The scattering coefficients 
$a(z)$ and 
$\bar{a}(z)$ can be analytically continued into the upper and lower half planes, respectively, while 
$b(z)$ and 
$\bar{b}(z)$ are only defined on the real axis. Importantly, since the eigenfunctions are defined as simultaneous solutions of the Lax pair, the scattering coefficients are time-independent. Specifically, they can be written as:
\begin{equation}
a(z)=\frac{W(\phi(x,t,z),\psi(x,t,z))}{\gamma(z)},\qquad \bar{a}(z)=\frac{W(\bar{\psi}(x,t,z),\bar{\phi}(x,t,z))}{\gamma(z)},
\end{equation}
\begin{equation}
b(z)=\frac{W(\bar\psi(x,t,z),\phi(x,t,z))}{\gamma(z)},\qquad \bar{b}(z)=\frac{W(\bar\phi(x,t,z),\psi(x,t,z))}{\gamma(z)},
\end{equation}
$W$ denotes the Wronskian determinant, showing that generically all scattering coefficients are singular as 
$z\to \pm q_0$, unless the Jost functions become linearly dependent at either 
$\pm q_{0}$, in which case 
$\pm q_{0}$ is called a ‘virtual level’, and the scattering coefficients are 
$\mathcal O(1)$ as 
$z$ approaches 
$\pm q_0$. Importantly, for any reflectionless/pure soliton potential, both 
$\pm q_0$ are virtual levels, and all scattering coefficients are finite as 
$z\to \pm q_0$.
The eigenfunctions satisfy the following symmetries corresponding to the involution 
$z\mapsto z^{*}$:
\begin{equation}
\Phi(x,t,z)=\sigma_{1}\Phi^{*}(x,t,z^{*})\sigma_{1},\qquad
\Psi(x,t,z)=\sigma_{1}\Psi^{*}(x,t,z^{*})\sigma_{1},
\end{equation}which in turn imply the following symmetries for the scattering coefficients:
\begin{equation}
a(z)=\bar{a}^{*}(z^{*})\quad \text{for } \mathop{\rm Im}\nolimits z\geq 0, \qquad
b(z)=\bar{b}^{*}(z) \quad \text{for } z\in\mathbb{R}.
\end{equation} Due to a second involution in the Lax pair, namely 
$z\mapsto q_{0}^{2}/z$, we also have the symmetries:
\begin{equation}
\Phi(x,t,z)=-\frac{1}{iz}\Phi(x,t,q_{0}^{2}/z)\sigma_{3}Q_{-}\, ,\qquad
\Psi(x,t,z)=\frac{1}{iz}\Psi(x,t,q_{0}^{2}/z)\sigma_{3}Q_{+}\, ,
\end{equation}and the scattering coefficients satisfy:
\begin{equation}
a(q_0^2/z)=e^{2i\theta}\bar{a}(z)\quad \text{for } \mathop{\rm Im}\nolimits z\geq 0, \qquad
b(q_0^2/z)=-\bar{b}(z)\quad \text{for } z\in\mathbb{R}.
\end{equation}Eq. (2.7) can be written as:
\begin{equation}
\frac{\phi(x,t,z)}{a(z)}=\bar{\phi}(x,t,z)+\rho(z)\psi(x,t,z),\qquad
\frac{\bar \phi(x,t,z)}{\bar a(z)}={\psi}(x,t,z)+\bar\rho(z)\bar \psi(x,t,z),
\end{equation}where the reflection coefficients are given by:
\begin{equation}
\rho(z)=\frac{b(z)}{a(z)},\qquad \bar\rho(z)=\frac{\bar b(z)}{\bar a(z)}.
\end{equation} Generically, it is known that 
$a(z)$ has a finite number of simple zeros that lie on the circle 
$|z|=q_{0}$ (minus the branch points 
$\pm q_{0}$), and due to the symmetry (2.10), if 
$a(\zeta_{j})=0$ then 
$\bar{a}(\zeta_{j}^{*})=0$. So, the discrete eigenvalues come in pairs
\begin{equation}
\{\zeta_{j},\zeta_{j}^{*}\}_{j=1}^{J},\qquad \zeta_{j}=q_{0}e^{i\alpha_{j}},\qquad 0 \lt \alpha_{j} \lt \pi.
\end{equation}At each discrete eigenvalue, (2.8) implies that the Jost eigenfunctions become proportional:
\begin{equation}
\phi(x,t,\zeta_{j})=b_{j}\psi(x,t,\zeta_{j}),\qquad \bar\phi(x,t,\zeta_{j}^{*})={b}^{*}_{j}\bar\psi(x,t,\zeta_{j}^{*}),
\end{equation}for some constant 
$b_{j}\in \mathbb{C}$ satisfying the symmetry 
$b_{j}^{*}=-b_{j}$. With this in mind, the residue contributions from each pair of discrete eigenvalues are
\begin{equation}
\mathop{\rm Res}\limits_{z=\zeta_{j}}\frac{\phi(x,t,z)}{a(z)}=C_{j}\psi(x,t,\zeta_{j}),\qquad
\mathop{\rm Res}\limits_{z=\zeta_{j}^{*}}\frac{\bar \phi(x,t,z)}{\bar a(z)}=C_{j}^{*}\bar \psi(x,t,\zeta^{*}_{j}),
\end{equation}where 
$C_{j}=b_{j}/a^\prime(\zeta_{j})$ is the norming constant associated with the discrete eigenvalue (hereafter, 
$^\prime$ denotes differentiation with respect to the scattering parameter 
$z$). Moreover, the first symmetry in (2.12) gives:
\begin{equation*}
\bar{a}'(z)=-e^{-2i\theta}\frac{q_0^2}{z^2}a'(q_0^2/z),
\end{equation*}and therefore the norming constants are such that
\begin{equation}
C_{j}^{*}=
e^{2i(\theta-\alpha_{j})}
C_{j}.
\end{equation}Note that, like the scattering coefficients and the reflection coefficients, the norming constants are time-independent because the Jost eigenfunctions have been defined as simultaneous solutions of the Lax pair.
The inverse scattering problem is then solved by formulating (2.13) as an RHP across the real 
$z$ axis. After applying Cauchy projectors to (2.13) and accounting for the residues, the RHP for the Jost eigenfunctions 
$\psi(x,t,z)$ and 
$\bar{\psi}(x,t,z)$ is converted into the following linear system of algebraic-integral equations:
\begin{align}
\psi(x,t,z)e^{-i\Omega(x,t,z)}&=\begin{bmatrix}
-iq_{+}/z\\1
\end{bmatrix}+\sum_{j=1}^{J}\frac{e^{-i\Omega(x,t,\zeta_{j}^*)}\bar{\psi}(x,t,\zeta_{j}^{*})C_{j}^{*}}{z-\zeta_{j}^{*}} \nonumber\\
&\quad -\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{e^{-i\Omega(x,t,\zeta)}\bar \psi(x,t,\zeta)\rho^{*}(\zeta)}{\zeta-(z+i0)}d\zeta,
\end{align}
\begin{align}
\bar \psi(x,t,z)e^{i\Omega(x,t,z)}&=\begin{bmatrix}
1\\ iq_{+}^{*}/z
\end{bmatrix}+\sum_{j=1}^{J}\frac{e^{i\Omega(x,t,\zeta_{j})}\psi(x,t,\zeta_{j})C_{j}}{z-\zeta_{j}} \nonumber\\
&\quad +\,\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{e^{i\Omega(x,t,\zeta)}\psi(x,t,\zeta)\rho(\zeta)}{\zeta-(z-i0)}d\zeta.
\end{align}
$q(x,t)$ to (1.1)) can be recovered from the large-
$z$ asymptotic behaviour of the eigenfunctions.
3. Completeness of squared eigenfunctions for general potentials
The main purpose of this section is to prove completeness of the squared eigenfunctions for general potentials. We first show how to calculate variations of the scattering data in terms of a variation of the potential, and vice versa, which in turn will yield the adjoint squared eigenfunctions and the squared eigenfunctions. Then, we use the variations to prove the completeness result. The methodology follows [Reference Yang54], while also accounting for the nonzero boundary conditions. Throughout this section, we suppress explicit dependence on 
$t$.
3.1. Variation in the scattering data from a variation in the potential
Since the eigenfunction 
$\Phi$ defined in (2.4a) satisfies (2.2a), a variation 
$\delta\Phi$ in the eigenfunction corresponding to a variation 
$\delta Q$ in the potential must satisfy
\begin{equation}
\delta\Phi_{x}=-ik\sigma_{3}\delta\Phi+Q\delta\Phi+(\delta Q)\Phi\,,\qquad \lim_{x\rightarrow-\infty}\delta\Phi=0.
\end{equation}One can check that the solution of (3.1) is given by
\begin{equation}
\delta\Phi(x,z)=\Phi(x,z)\int_{-\infty}^{x}\Phi^{-1}(y,z)\delta Q(y)\Phi(y,z)dy.
\end{equation} Taking the limit as 
$x\rightarrow+\infty$, we have 
$\Phi=\Psi S\rightarrow E_{+}e^{-i\Omega\sigma_{3}}S$, so that (3.2) becomes
\begin{equation}
\delta S(z)=\int_{-\infty}^{\infty}\Psi^{-1}(y,z)\delta Q(y)\Phi(y,z)dy.
\end{equation}If we write the eigenfunctions as
\begin{equation}
\Phi=\begin{bmatrix}
\phi&\bar\phi
\end{bmatrix}=\begin{bmatrix}
\phi_{1}&\bar\phi_{1}\\\phi_{2}&\bar\phi_{2}
\end{bmatrix},\qquad \Psi=\begin{bmatrix}
\bar\psi&\psi
\end{bmatrix}=\begin{bmatrix}
\bar\psi_{1}&\psi_{1}\\\bar\psi_{2}&\psi_{2}
\end{bmatrix},
\end{equation}then their inverses are
\begin{equation}
\Phi^{-1}=\frac{1}{\gamma(z)}\begin{bmatrix}
\hat{\bar\phi}\\\hat\phi
\end{bmatrix}=\frac{1}{\gamma(z)}\begin{bmatrix}
\bar\phi_{2}&-\bar\phi_{1}\\
-\phi_{2}&\phi_{1}
\end{bmatrix}\,,\qquad \Psi^{-1}=\frac{1}{\gamma(z)}\begin{bmatrix}
\hat\psi\\\hat{\bar\psi}
\end{bmatrix}=\frac{1}{\gamma(z)}\begin{bmatrix}
\psi_{2}&-\psi_{1}\\
-\bar\psi_{2}&\bar\psi_{1}
\end{bmatrix},
\end{equation}where we recall that 
$\gamma(z)=\det\Phi=\det\Psi=1-q_{0}^{2}/z^{2}$ (cf (2.6)).
With these definitions, (3.3) gives entrywise:
\begin{equation}
\delta a(z)=\frac{1}{\gamma(z)}\int_{-\infty}^{\infty}\hat\psi(y,z)\delta Q(y)\phi(y,z)dy,\quad
\delta b(z)=\frac{1}{\gamma(z)}\int_{-\infty}^{\infty}\hat{\bar\psi}(y,z)\delta Q(y)\phi(y,z)dy,
\end{equation}
\begin{equation}
\delta \bar a(z)=\frac{1}{\gamma(z)}\int_{-\infty}^{\infty}\hat{\bar\psi}(y,z)\delta Q(y)\bar\phi(y,z)dy,\quad
\delta \bar b(z)=\frac{1}{\gamma(z)}\int_{-\infty}^{\infty}\hat\psi(y,z)\delta Q(y)\bar\phi(y,z)dy.
\end{equation}
$\rho=b/a$ and 
$\bar\rho=\bar{b}/\bar{a}$ (cf. (2.14)) are found to be
\begin{equation}
\delta\rho(z)=\frac{1}{\gamma(z)a(z)^{2}}\int_{-\infty}^{\infty}\hat\phi(y,z)\delta Q(y)\phi(y,z)dy,
\end{equation}
\begin{equation}
\delta\bar\rho(z)=\frac{1}{\gamma(z)\bar a(z)^{2}}\int_{-\infty}^{\infty}\hat{\bar \phi}(y,z)\delta Q(y)\bar\phi(y,z)dy,
\end{equation}
$\hat\phi=a\hat{\bar\psi}-b\hat\psi$ and 
$\hat{\bar\phi}=\bar a \hat\psi-\bar b\hat{\bar\psi}$ have been used. If we define the inner product.
\begin{equation}
\langle f,g\rangle=\int_{-\infty}^{\infty}f(y)^{T}g(y)dy,
\end{equation}then (3.7) can be written as
\begin{equation}
\delta\rho(z)=\frac{1}{\gamma(z)a(z)^{2}}\left\langle\chi,\begin{bmatrix}
\delta q\\\delta q^{*}
\end{bmatrix}\right\rangle,\qquad
\delta\bar\rho(z)=\frac{1}{\gamma(z)\bar a(z)^{2}}\left\langle\bar\chi,\begin{bmatrix}
\delta q\\\delta q^{*}
\end{bmatrix}\right\rangle,
\end{equation}where the so-called adjoint squared eigenfunctions are defined as:
\begin{equation}
\chi(x,z)=\begin{bmatrix}
-\phi_{2}(x,z)^{2}\\\phi_{1}(x,z)^{2}
\end{bmatrix},\qquad \bar\chi(x,z)=\begin{bmatrix}
\bar\phi_{2}(x,z)^{2}\\-\bar\phi_{1}(x,z)^{2}
\end{bmatrix}.
\end{equation} [Note that it is customary to define the inner product as in (3.8), without explicit complex conjugation, with the implication that the left argument is a member of an appropriate dual space. The adjoint 
$M^{A}$ of an operator 
$M$ is defined such that 
$\langle M^{A}f,g\rangle=\langle f,Mg\rangle$.]
The analyticity properties of the Jost eigenfunctions in Section 2 imply that 
$\chi(x,z)$ is analytic for 
$\mathop{\rm Im}\nolimits z \gt 0$, and 
$\bar{\chi}(x,z)$ is analytic for 
$\mathop{\rm Im}\nolimits z \lt 0$. Note also that on account of the symmetries (2.9)-(2.10), we have that the two equations in (3.9) are equivalent.
Remark 1. Since 
$\gamma(\pm q_{0})=0$, two possible scenarios can occur: (i) If the unperturbed reflection coefficient is nonzero, then 
$a(z)$ has simple poles at 
$z=\pm q_{0}$, so we have 
$\delta\rho(\pm q_{0})=0$ and the property 
$|\rho(\pm q_{0})|=1$ is preserved. For the present discussion, we assume this to be the case. (ii) If the unperturbed reflection coefficient is zero, then 
$a(\pm q_{0})$ is finite, meaning that 
$\delta\rho$ becomes singular as 
$z\rightarrow\pm q_{0}$. Notably, this happens in the case of pure solitons, to be discussed later. The implication is that for infinitesimally small variations in the potential 
$\delta q$, the resulting variation in the reflection coefficient 
$\delta\rho$ is finite. This causes the generation of a radiation shelf, a phenomenon that is also observed in the perturbation theory associated with the KdV equation, due to a similar singularity that occurs as the spectral parameter approaches zero (see, e.g., [Reference Herman24, Reference Karpman and Maslov27]).
3.2. Variation in the potential from a variation in the scattering data
In matrix form, the relations between the Jost eigenfunctions given in (2.13) can be expressed as
\begin{equation}
\bar\mu=\mu G,\qquad
G=\begin{bmatrix}
1&e^{-2i\Omega(x,z)}\bar\rho(z)\\
-e^{2i\Omega(x,z)}\rho(z)&1-\rho(z)\bar\rho(z)
\end{bmatrix}, \qquad z\in \mathbb{R},
\end{equation}
\begin{equation}
\mu=\begin{bmatrix}
e^{i\Omega(x,z)}\displaystyle{\frac{\phi(x,z)}{a(z)}}\ & e^{-i\Omega(x,z)}\psi(x,z)
\end{bmatrix},\quad \bar\mu=\begin{bmatrix}
e^{i\Omega(x,z)}\bar{\psi}(x,z)\ & e^{-i\Omega(x,z)}
\displaystyle{\frac{\bar{\phi}(x,z)}{\bar{a}(z)}}
\end{bmatrix},
\end{equation}
$z$-axis. To start with, we assume that 
$a(z)$ and 
$\bar{a}(z)$ have no zeros, i.e., 
$\mu$ and 
$\bar\mu$ are analytic functions of 
$z$ in the upper and lower half planes, respectively. A variation of (3.11) reads:
\begin{equation}
\delta\bar\mu=(\delta \mu)\, G+\mu\, \delta G\,, \qquad \delta G=\begin{bmatrix}
0&e^{-2i\Omega(x,z)} \delta\bar\rho(z) \\-e^{2i\Omega(x,z)}\delta\rho(z)\ \ &-\rho(z)\delta\bar\rho(z)-\bar\rho(z)\delta\rho(z)
\end{bmatrix}.
\end{equation} If we write 
$G=\mu^{-1}\bar\mu$, and introduce the definitions
\begin{equation}
\nu=(\delta\mu)\mu^{-1},\qquad \bar\nu=(\delta\bar\mu) \bar\mu^{-1},\qquad F=-\mu(\delta G) \bar\mu^{-1},
\end{equation}then we have the RHP:
\begin{equation}
\nu-\bar\nu=F.
\end{equation} After simplification, 
$F$ is given by
\begin{equation}
F=\Psi\begin{bmatrix}
0&-\delta\bar\rho(z)\\
\delta\rho(z)&0
\end{bmatrix}\Psi^{-1}.
\end{equation} Using the large-
$z$ asymptotic behaviour of the Jost eigenfunctions (which can be found, for example, in [Reference Prinari46]) in their respective half-planes, we have that
\begin{equation}
\nu,\bar\nu\sim\begin{bmatrix}
\mathcal O(1/z)&-i\delta q/z\\
i\delta q^{*}/z&\mathcal O(1/z)
\end{bmatrix},\qquad|z|\rightarrow\infty.
\end{equation} Since 
$\nu,\bar\nu\rightarrow0$ as 
$|z|\rightarrow\infty$, the solution on (3.14) can be written down in terms of Cauchy projectors as
\begin{equation}
\nu(x,z)=\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{F(x,\zeta)}{\zeta-(z+i0)}d\zeta,\qquad
\bar\nu(x,z)=\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{F(x,\zeta)}{\zeta-(z-i0)}d\zeta.
\end{equation} Comparing the large-
$z$ asymptotic behaviour of (3.17) with (3.16) it can be deduced that
\begin{equation}
\begin{bmatrix}
\delta q(x)\\
\delta q^{*}(x)
\end{bmatrix}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{1}{\gamma(\zeta)}\Big\{\eta(x,\zeta)\delta\rho(\zeta)+\bar\eta(x,\zeta)\delta\bar\rho(\zeta)\Big\}d\zeta,
\end{equation}where the squared eigenfunctions are defined as:
\begin{equation}
\eta(x,z)=\begin{bmatrix}
\psi_{1}(x,z)^{2}\\
\psi_{2}(x,z)^{2}
\end{bmatrix},\qquad \bar\eta(x,z)=\begin{bmatrix}
\bar\psi_{1}(x,z)^{2}\\
\bar\psi_{2}(x,z)^{2}
\end{bmatrix}.
\end{equation} Again, the analyticity properties of the Jost eigenfunctions in Section 2 imply that 
$\eta(x,z)$ is analytic for 
$\mathop{\rm Im}\nolimits z \gt 0$, and 
$\bar{\eta}(x,z)$ is analytic for 
$\mathop{\rm Im}\nolimits z \lt 0$. Moreover, due to the symmetries (2.11) and (2.12), the squared eigenfunctions and reflection coefficients satisfy
\begin{equation}
\bar\eta(x,z)=-\frac{q_{0}^{2}}{z^{2}}e^{-2i\theta}\eta(x,q_{0}^{2}/z),\qquad\delta\bar\rho(z)=-e^{2i\theta}\delta\rho(q_{0}^{2}/z),
\end{equation}which allows one to combine the two terms in (3.18) into one by making a change of variables 
$\zeta'=q_{0}^{2}/\zeta$ in the second integral:
\begin{equation}
\begin{bmatrix}
\delta q(x)\\
\delta q^{*}(x)
\end{bmatrix}=
\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{1}{\gamma(\zeta)}\eta(x,\zeta)\delta\rho(\zeta)\left(1-\frac{q_{0}^{2}}{\zeta^{2}}\right)d\zeta=
\frac{1}{2\pi}\int_{-\infty}^{\infty}\eta(x,\zeta)\delta\rho(\zeta)d\zeta.
\end{equation}This reduction is in contrast to the case of zero boundary conditions, where both squared eigenfunctions need to be retained since there is no symmetry analogous to (3.20).
3.3. Completeness relation for general potentials
As shown in the last section, a variation in the potential induced by a variation in the reflection coefficient can be expressed in terms of the squared eigenfunction 
$\eta(x,z)$ defined in (3.19) via Eq. (3.21). Conversely, a variation in the reflection coefficient induced by a variation in the potential is expressed in terms of the adjoint squared eigenfunction 
$\chi(x,z)$ in (3.10) via (3.9). Inserting (3.21) into (3.9), we have
\begin{equation}
\delta\rho(z)=\frac{1}{2\pi\gamma(z)a(z)^{2}}\int_{-\infty}^{\infty}\big\langle\chi(z),\eta(z)\big\rangle\delta\rho(\zeta)d\zeta,
\end{equation}from which we obtain the inner product between the squared eigenfunction 
$\eta$ and its adjoint 
$\chi$:
\begin{equation}
\big\langle\chi(z),\eta(\zeta)\big\rangle=2\pi\gamma(z)a(z)^{2}\delta(z-\zeta).
\end{equation}On the other hand, substituting (3.9) into (3.21) and using the definition of the inner product (3.8) gives
\begin{equation}
\begin{bmatrix}
\delta q(x)\\
\delta q^{*}(x)
\end{bmatrix}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{1}{\gamma(\zeta)a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta\begin{bmatrix}
\delta q(y)\\\delta q^{*}(y)
\end{bmatrix}dy,
\end{equation}which leads to the closure/completeness relation:
\begin{equation}
\int_{-\infty}^{\infty}\frac{1}{\gamma(\zeta)a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta=2\pi\delta(x-y)I.
\end{equation} Recall that both (3.23) and (3.25) have been derived under the assumption that 
$a(z)$ has no zeros and that the unperturbed reflection coefficient 
$\rho(z)$ is nonzero. If this is the case, the integrand has no poles, the set 
$\{\eta(x,z):z\in\mathbb{R}\}$ is complete, and (3.25) is the closure relation. Generically, 
$a(z)$ has a finite number of simple zeros 
$\{z=\zeta_{j},\;\text{Im}\,\zeta_{j} \gt 0\}_{j=1}^{J}$ that lie on the circle 
$|z|=q_{0}$. Furthermore, in any reflectionless case, poles at 
$z=\pm q_{0}$ coming from the fact that 
$\gamma(\pm q_{0})=0$ also must be accounted for. Therefore, a more general closure relation is required. The general closure relation can be obtained by simply replacing the integral in (3.25) with one taken over a contour in the upper-half 
$z$-plane from 
$-\infty+i0$ to 
$+\infty+i0$ that passes above the circle 
$|z|=q_{0}$. Moreover, when discrete eigenvalues are present (i.e., when 
$a(z)$ has zeros) the orthogonality of the continuous eigenfunctions as given in (3.23) still holds for 
$z,\zeta\in\mathbb{R}$, and will be supplemented with additional orthogonality relations between the discrete squared eigenfunctions that arise due to the residues at the complex zeros of 
$a(z)$.
The most compact form of the general closure relation is:
\begin{equation}
\int_{\Gamma}\frac{1}{\gamma(\zeta)a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta=2\pi\delta(x-y)I,
\end{equation}where 
$\Gamma$ is a contour from 
$-\infty+i0$ to 
$+\infty+i0$ that passes above the circle of radius 
$q_{0}$ in the upper half plane. A justification of the generalisation of (3.25) to the contour integral (3.26) is provided in Appendix A. Explicitly accounting for the residue contributions from the poles at the discrete eigenvalues 
$\zeta_{j}$ (given in Appendix A), the closure relation can be written as
\begin{align}
2\pi\delta(x-y)I &=\int_{C}\frac{1}{\gamma(\zeta)a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta \nonumber\\ &\quad - 2\pi i\sum_{j=1}^{J}\frac{1}{\gamma(\zeta_{j}){a}'(\zeta_{j})^{2}}\Big[\eta(x,\zeta_{j})\Theta(y,\zeta_{j})^{T}+\eta'(x,\zeta_{j})\chi(y,\zeta_{j})^{T}\Big],\end{align}where 
$C$ is a contour along the real axis indented in the upper-half plane to avoid 
$\pm q_{0}$, and we have defined
\begin{equation}
\Theta(x,\zeta_{j})=\chi'(x,\zeta_{j})+\frac{\gamma(\zeta_{j}){a}''(\zeta_{j})-\gamma'(\zeta_{j}){a}'(\zeta_{j})}{\gamma(\zeta_{j}){a}'(\zeta_{j})}\chi(x,\zeta_{j}).
\end{equation} As before, prime denotes differentiation with respect to 
$z$. Note that in the case of zero boundary conditions, the closure relation contains two continuous and four discrete squared eigenfunctions. In the present case, the number of independent eigenfunctions is halved by accounting for the symmetry 
$z\mapsto q_{0}^{2}/z$.
Finally, we can account for the possibility of simple poles on the real axis at 
$z=\pm q_{0}$ (which are due to 
$\gamma(\pm q_0)=0$, and are present whenever 
$a(\pm q_0)$ is finite, which occurs for reflectionless potentials) by replacing the integral over 
$C$ with a principal value integral over the real axis, while accounting for the residues
\begin{equation}
\mathop{\rm Res}\limits_{z=\pm q_0}\frac{\eta(x,z)\chi(y,z)^{T}}{\gamma(z)a(z)^{2}}= \frac{1}{\gamma'(\pm q_{0})a(\pm q_{0})^{2}}\eta(x,\pm q_{0})\chi(y,\pm q_{0})^{T},
\end{equation}with a factor of 
$\pi i$. Again, we stress that (3.29) is only nonzero in the reflectionless case, for which 
$a(\pm q_{0})$ is finite. In the case of nonzero reflection, where 
$a(z)$ has poles at 
$\pm q_{0}$ (see Section 2), the earlier form of the closure relation (3.27) (with integral taken over the real axis rather than 
$C$) is sufficient. With these residues incorporated, the full closure relation is:
\begin{align}2\pi\delta(x-y)I&={\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{1}{\gamma(\zeta)a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta-\pi i\sum_{\pm}\frac{1}{\gamma'(\pm q_{0})a(\pm q_{0})^{2}}\eta(x,\pm q_{0})\chi(y,\pm q_{0})^{T}\nonumber\\ &\quad - 2\pi i\sum_{j=1}^{J}\frac{1}{\gamma(\zeta_{j}){a}'(\zeta_{j})^{2}}\Big[\eta(x,\zeta_{j})\Theta(y,\zeta_{j})^{T}+\eta'(x,\zeta_{j})\chi(y,\zeta_{j})^{T}\Big].\end{align} Here and throughout the paper, 
${\int{\kern-8.5pt}-}$ denotes the Cauchy principal value integral. It is worth noting at this point that while several previous works (including [Reference Chen, Chen and Huang10, Reference Konotop and Vekslerchik39]) derived essentially equivalent completeness relations, the half residues due to the poles at 
$\pm q_{0}$ were not explicitly accounted for. The inclusion of these contributions plays an important role in the perturbation theory, on which we will elaborate in Section 6.
4. Linearisation operator and completeness relation for a single soliton
In this section, we discuss the linearisation operator of the NLS equation with nonzero boundary conditions and its relation to the squared eigenfunctions, and we derive the explicit expression of the completeness relation for a 1-soliton solution.
4.1. The linearisation operator
Upon directly taking a variation of the NLS equation (1.1), we get
\begin{equation}
i\delta q_{t}(x,t)+\delta q_{xx}(x,t)-2(2|q(x,t)|^{2}-q_{0}^{2})\delta q(x,t)-2q(x,t)^{2}\delta q^{*}(x,t)=0.
\end{equation}Combining (4.1) with its complex conjugate, a variation in the potential should satisfy
\begin{equation}
\mathcal{L}\begin{bmatrix}\delta q(x,t)\\\delta q^{*}(x,t)\end{bmatrix}=0,
\end{equation}where the linearisation operator is
\begin{equation}
\mathcal{L}=\begin{bmatrix}
i\partial_{t}+\partial_{xx}-2(2|q(x,t)|^{2}-q_{0}^{2})&-2q(x,t)^{2}\\
2q^{*}(x,t)^{2}&i\partial_{t}-\partial_{xx}+2(2|q(x,t)|^{2}-q_{0}^{2})
\end{bmatrix}.
\end{equation} Inserting (3.21) into (4.2) shows that the squared eigenfunction 
$\eta$ is an eigenfunction of the linearisation operator with zero eigenvalue, namely
\begin{equation}
\mathcal{L}\eta(x,t,z)=0.
\end{equation} Moreover, the formal adjoint of 
$\mathcal{L}$ is
\begin{equation}
\mathcal{L}^{A}=\begin{bmatrix}
-i\partial_{t}+\partial_{xx}-2(2|q(x,t)|^{2}-q_{0}^{2})&2q^{*}(x,t)^{2}\\
-2q(x,t)^{2}&-i\partial_{t}-\partial_{xx}+2(2|q(x,t)|^{2}-q_{0}^{2})
\end{bmatrix},
\end{equation}and from
\begin{equation}
\langle\chi(z),\mathcal{L}\eta(\zeta)\rangle=\langle\mathcal{L}^{A}\chi(z),\eta(\zeta)\rangle=0,
\end{equation}we deduce that the adjoint squared eigenfunction 
$\chi$ satisfies
\begin{equation}
\mathcal{L}^{A}\chi(x,t,z)=0.
\end{equation}Apart from an arbitrary phase, the exact single dark soliton solution of (1.1) is
\begin{equation}
q(x,t)=u(\xi)=k_{1}+i\Lambda_{1}\tanh\xi,\qquad \xi\equiv\Lambda_{1}(x+2k_{1}t-x_{1}),
\end{equation}where 
$\zeta_{1}=k_{1}+i\Lambda_{1}=q_{0}e^{i\theta}$. Thus, when performing a linearisation around the 1-soliton solution, it proves convenient to express the linearisation operator and its adjoint in terms of the travelling coordinate 
$\xi$ rather than 
$x$ (while allowing for potential explicit dependence on 
$t$). In this case, we have 
$\partial_{t}\rightarrow\partial_{t}+2\Lambda_{1}k_{1}\partial_{\xi}$, 
$\partial_{x}\rightarrow\Lambda_{1}\partial_{\xi}$ and
\begin{equation}
\mathcal{L}=i\partial_{t}+L,\qquad \mathcal{L}^{A}=-i\partial_{t}+L^{A},
\end{equation}where
\begin{equation}
L=\begin{bmatrix}
\Lambda_{1}^{2}\partial_{\xi}^{2}+2i\Lambda_{1}k_{1}\partial_{\xi}-2(2|u(\xi)|^{2}-q_{0}^2)&-2u(\xi)^{2}\\
2u^{*}(\xi)^{2}&-\Lambda_{1}^{2}\partial_{\xi}^{2}+2i\Lambda_{1}k_{1}\partial_{\xi}+2(2|u(\xi)|^{2}-q_{0}^{2})
\end{bmatrix},
\end{equation}
\begin{equation}
L^{A}=\begin{bmatrix}
\Lambda_{1}^{2}\partial_{\xi}^{2}-2i\Lambda_{1}k_{1}\partial_{\xi}-2(2|u(\xi)|^{2}-q_{0}^2)&2u(\xi)^{*2}\\
-2u(\xi)^{2}&-\Lambda_{1}^{2}\partial_{\xi}^{2}-2i\Lambda_{1}k_{1}\partial_{\xi}+2(2|u(\xi)|^{2}-q_{0}^{2})
\end{bmatrix}.
\end{equation}
$\psi(x,t,z)$ and 
$\bar{\psi}(x,t,z)$ with 
$J=1$, 
$\rho(z)\equiv0$, and discrete eigenvalue and norming constant given by:
\begin{equation}
\zeta_{1}=k_{1}+i\Lambda_{1}, \qquad C_{1}=-2\Lambda_{1}e^{2\Lambda_{1}x_{1}},
\end{equation}and 
$\phi(x,t,z)$ can be recovered from the above via:
\begin{equation}
\phi(x,t,z)=a(z)\bar{\psi}(x,t,z),
\end{equation} (cf. (2.13) with 
$\rho(z)\equiv 0$) as well as the explicit expression for 
$a(z)$ in the 1-soliton case, namely:
\begin{equation}
a(z)=\frac{z-\zeta_1}{z-\zeta_1^*}.
\end{equation}Eqs. (3.19) and (3.10) then give the following explicit expression for the squared eigenfunctions:
\begin{equation}
\eta(x,t,z)=\begin{bmatrix}
[iq_{+}z^{-1}(z-\zeta_{1}^{*})+\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2\\
[(z-\zeta_{1}^{*})-i\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2\end{bmatrix}\frac{e^{2i\Omega(x,t,z)}}{(z-\zeta_{1}^{*})^{2}},
\end{equation}
\begin{equation}
\chi(x,t,z)=\begin{bmatrix}
-[iq_{-}z^{-1}(z-\zeta_{1})-\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2\\
[(z-\zeta_{1})+i\Lambda_{1} e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2
\end{bmatrix}\frac{e^{-2i\Omega(x,t,z)}}{(z-\zeta_{1}^{*})^{2}}.
\end{equation}
\begin{equation}
\mathcal{Z}(\xi,z)=\eta(x,t,z)e^{-4i\lambda(k-k_{1})t-2i\lambda x_{1}}(z-\zeta_{1}^{*})^{2},
\end{equation}
\begin{equation}
\Upsilon(\xi,z)=\chi(x,t,z)e^{4i\lambda(k-k_{1})t+2i\lambda x_{1}}(z-\zeta_{1}^{*})^{2},
\end{equation}
$z$-dependent denominators have been removed to simplify later calculations, and 
$k=k(z)$ and 
$\lambda=\lambda(z)$ are as defined in (2.3). Straightforward calculations yield the explicit forms:
\begin{equation}
\mathcal{Z}(\xi,z)=\begin{bmatrix}
[iq_{+}z^{-1}(z-\zeta_{1}^{*})+\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2\\
[(z-\zeta_{1}^{*})-i\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2\end{bmatrix}e^{2i\frac{\lambda(z)}{\Lambda_{1}}\xi},
\end{equation}
\begin{equation}
\Upsilon(\xi,z)=\begin{bmatrix}
-[iq_{-}z^{-1}(z-\zeta_{1})-\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2\\
[(z-\zeta_{1})+i\Lambda_{1} e^{-\xi}\mathop{\rm sech}\nolimits\xi]^2
\end{bmatrix}e^{-2i\frac{\lambda(z)}{\Lambda_{1}}\xi}.
\end{equation}
\begin{equation}
\langle\Upsilon(z),\mathcal{Z}(\zeta)\rangle=\frac{1}{\Lambda_{1}}\int_{-\infty}^{\infty}\Upsilon(\xi',z)^{T}\mathcal{Z}(\xi',\zeta)d\xi'=2\pi\gamma(z)|z-\zeta_{1}|^{4}\delta(z-\zeta).
\end{equation} Note that the division by 
$\Lambda_{1}$ is necessary because, here and throughout the paper, the inner product is taken with respect to 
$x$ and not 
$\xi$. Substituting (4.15a) into (4.4) and (4.15b) into (4.7) with (4.9) in mind, we find that 
$\mathcal{Z}$ and 
$\Upsilon$ are eigenfunctions of 
$L$ and 
$L^{A}$, respectively:
\begin{equation}
L\mathcal{Z}(\xi,z)=4\lambda(k-k_{1})\mathcal{Z}(\xi,z),\qquad
L^{A}\Upsilon(\xi,z)=4\lambda(k-k_{1})\Upsilon(\xi,z).
\end{equation} The discrete eigenfunctions associated with the soliton eigenvalue 
$\zeta_{1}$ can be obtained by directly evaluating (4.16) at 
$z=\zeta_1$, from which we get
\begin{equation}
\mathcal{Z}_1(\xi):= \mathcal{Z}(\xi,\zeta_{1})=\Lambda_{1}^{2}\begin{bmatrix}
1\\-1
\end{bmatrix}\mathop{\rm sech}\nolimits^{2}\xi,\qquad
\Upsilon_1(\xi):=\Upsilon(\xi,\zeta_{1})=-\Lambda_{1}^{2}\begin{bmatrix}
1\\1
\end{bmatrix}\mathop{\rm sech}\nolimits^{2}\xi.
\end{equation} Furthermore, setting 
$z=\zeta_{1}$ in (4.18) shows that
\begin{equation}
L\mathcal{Z}_1(\xi)=
L^{A}\Upsilon_1(\xi)=0.
\end{equation} Differentiating (4.16) with respect to 
$z$, the generalised discrete eigenfunctions are found to be
\begin{equation}
\mathcal{Z}^\prime_1(\xi):=\left. \partial_{z}\mathcal{Z}(\xi,z)\right|_{z=\zeta_1}
=-2i\Lambda_{1} \{ \begin{bmatrix}q_{0}^{2}\zeta_{1}^{-2}\\-1\end{bmatrix}e^{-\xi}\mathop{\rm sech}\nolimits\xi-k_{1}\zeta_{1}^{-1}\begin{bmatrix}
1\\-1
\end{bmatrix}\xi\mathop{\rm sech}\nolimits^{2}\xi \},
\end{equation}
\begin{equation}
\Upsilon^\prime_1(\xi) :=\left. \partial_{z}\Upsilon(\xi,z)\right|_{z=\zeta_1}
=2i\Lambda_{1}\{\begin{bmatrix}q_{0}^{2}\zeta_{1}^{-2}\\1\end{bmatrix}e^{\xi}\mathop{\rm sech}\nolimits\xi+k_{1}\zeta_{1}^{-1}\begin{bmatrix}
1\\1
\end{bmatrix}\xi\mathop{\rm sech}\nolimits^{2}\xi\},
\end{equation}
\begin{equation}
L\mathcal{Z}^\prime_1(\xi)=-4\Lambda_{1}^{2}\zeta_{1}^{-1}\mathcal{Z}_1(\xi),\qquad
L^{A}\Upsilon^\prime_1(\xi)=-4\Lambda_{1}^{2}\zeta_{1}^{-1}\Upsilon_1(\xi).
\end{equation} It is worth noting that the identification of the discrete modes 
$\mathcal{Z}_1$ and 
${\mathcal{Z}}_1'$ explicitly resolves the algebraic multiplicity of the zero eigenvalue of the linearised operator 
$L$. The eigenfunction 
$\mathcal{Z}_1$ corresponds to the Goldstone mode from translational invariance, while the generalised eigenfunction 
${\mathcal{Z}}_1'$ corresponds to the mode from scaling invariance (see, e.g., [Reference Bilas and Pavloff7]), confirming that the zero eigenvalue is indeed of algebraic multiplicity two.
4.2. 1-soliton completeness relation
Here we specialise the completeness relation to the 1-soliton case. Returning to the generic completeness relation (3.26) and substituting (4.15) as well as (4.13) for 
$a(z)$, we have
\begin{align}
\int_{\Gamma}\frac{1}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\mathcal{Z}(\xi,\zeta)\Upsilon(\xi',\zeta)^{T}d\zeta =2\pi\delta(\xi-\xi')I.
\end{align} After calculating the residue due to 
$\zeta_{1}$, the completeness relation can be written as
\begin{align}2\pi\delta(\xi-\xi')I&=\int_{C}\frac{1}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\mathcal{Z}(\xi,\zeta)\Upsilon(\xi',\zeta)^{T}d\zeta\nonumber\\ &\quad +\,\frac{\pi\zeta_{1}}{4\Lambda_{1}^{3}} \Big[\mathcal{Z}_1(\xi)\Xi_1(\xi')^{T}+\mathcal{Z}^\prime_1(\xi)\Upsilon_1(\xi')^{T}\Big],\end{align}where, as before, 
$C$ is indented in the upper half plane to avoid 
$\pm q_{0}$, and we have defined the new generalised eigenfunction
\begin{equation}
\Xi_1(\xi)=\Upsilon_1^\prime(\xi)+\frac{2ik_{1}}{\Lambda_{1}\zeta_{1}}\Upsilon_1(\xi).
\end{equation} One can verify directly that 
$\langle\Xi_1(\xi),\mathcal{Z}^\prime_1(\xi)\rangle=0$. Now, we account for the residues at the branch points. Setting 
$z=\pm q_{0}$ in (4.18) and noting that 
$\lambda(\pm q_{0})=0$ shows that
\begin{equation}
L\mathcal{Z}^\pm_0(\xi)=L^{A}\Upsilon^\pm_0(\xi)=0, \qquad
\mathcal{Z}^\pm_0(\xi):=\mathcal{Z}(\xi,\pm q_{0}), \quad
\Upsilon_0^\pm(\xi):=\Upsilon(\xi, \pm q_0).
\end{equation} A straightforward calculation shows that the squared eigenfunctions evaluated at 
$\pm q_{0}$ can be written in terms of the dark soliton 
$u(\xi)$ as defined in (4.8) as well as the discrete eigenfunctions given in (4.19) in the following way:
\begin{equation}
\mathcal{Z}_0^\pm(\xi)=2(k_{1}\mp q_{0})\begin{bmatrix}
-u(\xi)\\u^{*}(\xi)\end{bmatrix}-\mathcal{Z}_1(\xi),\qquad
\Upsilon_0^\pm (\xi)=2(k_{1}\mp q_{0})\begin{bmatrix}
u^{*}(\xi)\\u(\xi)\end{bmatrix}-\Upsilon_1(\xi).
\end{equation} Pulling out the residue contribution to the integral in (4.24) due to 
$\pm q_{0}$, the 1-soliton completeness relation can be written as:
\begin{align}
2\pi\delta(\xi-\xi')I&={\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{1}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\mathcal{Z}(\xi,\zeta)\Upsilon(\xi',\zeta)^{T}d\zeta \nonumber\\
&\quad -\,{\pi i}\sum_{\pm}\frac{\pm1}{{8q_{0}}(k_{1}\mp q_{0})^{2}}\mathcal{Z}_0^\pm(\xi)\Upsilon_0^\pm(\xi')^{T} \nonumber \\
&\quad +\,\frac{\pi\zeta_{1}}{4\Lambda_{1}^{3}} \Big[\mathcal{Z}_1(\xi)\Xi_1(\xi')^{T}+\mathcal{Z}^\prime_1(\xi)
\Upsilon_1(\xi')^{T}\Big].
\end{align} Although it may appear that the contributions from 
$\pm q_{0}$ generate additional discrete eigenfunctions, it is actually the case that these contributions modify the discrete eigenfunctions already present in the completeness relation. Indeed, using (4.27), we find that
\begin{align}
\sum_{\pm}\frac{\pm1}{{8q_{0}}(k_{1}\mp q_{0})^{2}}\mathcal{Z}_0^\pm(\xi)\Upsilon_0^\pm(\xi')^{T}&=\frac{1}{2\Lambda_{1}^{2}}\Bigg\{\begin{bmatrix}
-u(\xi)\\u^{*}(\xi)
\end{bmatrix}\Upsilon_1(\xi')^{T} \nonumber \\
&\quad +\,\mathcal{Z}_1(\xi)\begin{bmatrix}
u^{*}(\xi')\\u(\xi')
\end{bmatrix}^{T}+\frac{k_{1}}{\Lambda_{1}^{2}}\mathcal{Z}_1(\xi)\Upsilon_1(\xi')^{T}\Bigg\}.\end{align} The last term in (4.29) can be combined with the 
$\mathcal{Z}_1\Xi_1^{T}$ term in (4.28). Specifically, we have
\begin{equation}
\frac{\pi\zeta_{1}}{4\Lambda_{1}^{3}}\mathcal{Z}_1\Xi_1^{T}-\frac{\pi k_{1}}{2\Lambda_{1}^{4}}\mathcal{Z}_1\Upsilon_1^{T}=\frac{\pi\zeta_{1}}{4\Lambda_{1}^{3}}\mathcal{Z}_1\Upsilon_1^{\prime\, T},
\end{equation}which then gives the completeness relation in the form:
\begin{align}
2\pi\delta(\xi-\xi')I&={\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{1}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\mathcal{Z}(\xi,\zeta)\Upsilon(\xi',\zeta)^{T}d\zeta \nonumber \\
&\quad -\,\frac{\pi i}{2\Lambda_{1}^{2}}\left\{\begin{bmatrix}
-u(\xi)\\u^{*}(\xi)
\end{bmatrix}\Upsilon_1(\xi')^{T} +\mathcal{Z}_1(\xi)\begin{bmatrix}
u^{*}(\xi')\\u(\xi')
\end{bmatrix}^{T}\right\}\nonumber\\ &\quad +\,\frac{\pi\zeta_{1}}{4\Lambda_{1}^{3}} \Big[\mathcal{Z}_1(\xi)\Upsilon_1^\prime(\xi')^{T}+\mathcal{Z}_1^\prime(\xi)\Upsilon_1(\xi')^{T}\Big].
\end{align} Finally, the remaining two terms in braces can be absorbed into the 
$\mathcal{Z}_1\Upsilon_1^{\prime\, T}$ and 
$\mathcal{Z}^\prime_1\Upsilon^{T}_1$ terms. In particular, if we define
\begin{equation}
\tilde{\mathcal{Z}}_1(\xi)=\mathcal{Z}_1^\prime(\xi)-\frac{2i\Lambda}{\zeta_{1}}\begin{bmatrix}
-u(\xi)\\u^{*}(\xi)
\end{bmatrix},\qquad
\tilde{\Upsilon}_1(\xi)=\Upsilon_1^\prime(\xi)-\frac{2i\Lambda_{1}}{\zeta_{1}}\begin{bmatrix}
u^{*}(\xi)\\u(\xi)
\end{bmatrix},
\end{equation}then the final expression for the 1-soliton closure relation is:
\begin{align}2\pi\delta(\xi-\xi')I&={\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{1}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\mathcal{Z}(\xi,\zeta)\Upsilon(\xi',\zeta)^{T}d\zeta\nonumber\\ &\quad +\,\frac{\pi\zeta_{1}}{4\Lambda_{1}^{3}} \Big[\mathcal{Z}_1(\xi)\tilde\Upsilon_1(\xi')^{T}+\tilde{\mathcal{Z}}_1(\xi)\Upsilon_1(\xi')^{T}\Big].\end{align} The new generalised eigenfunction 
$\tilde{\mathcal{Z}}_{1}$ is a linear combination of the Goldstone modes from scaling invariance and from phase invariance. Note that the discrete eigenfunctions satisfy the following orthogonality conditions (which can either be read off from the closure relation or verified directly):
\begin{equation} \langle\Upsilon_{1},\mathcal{Z}_{1}\rangle=\langle\tilde\Upsilon_{1},\tilde{\mathcal{Z}}_{1}\rangle=0,\qquad\langle\tilde\Upsilon_{1},\mathcal{Z}_{1}\rangle=\langle\Upsilon_{1},\tilde{\mathcal{Z}}_{1}\rangle=8\Lambda_{1}^{3}\zeta_{1}^{-1}.
\end{equation}5. Perturbation theory
Consider the defocusing NLS equation with an additional small forcing perturbation:
\begin{equation}
iq_{t}(x,t)+q_{xx}(x,t)-2(|q(x,t)|^{2}-q_{0}^{2})q(x,t)=\varepsilon F[q(x,t)], \qquad 0 \lt \varepsilon \ll 1,
\end{equation}and introduce the slow time scale 
$T=\varepsilon t$ so that 
$\partial_{t}\rightarrow\partial_{t}+\varepsilon\partial_{T}$.
Remark 2. An exact dark soliton solution of (5.1) with 
$\varepsilon=0$ has the asymptotic behaviour 
$q(x,t)\rightarrow q_{0}e^{i\theta_{\pm}}$ as 
$x\rightarrow\pm\infty$ with generic phases 
$\theta_{\pm}$. Following the technique of [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1, Reference Chernyavskiy, Frantzeskakis, Horikis, Koutsokostas and Prinari11], substituting this into (5.1) in the limit 
$x\rightarrow\pm\infty$ while allowing 
$q_{0}$ and 
$\theta_{\pm}$ to depend on 
$T$ yields
\begin{equation}
q_{0T}=\mathop{\rm Im}\nolimits\left\{e^{-i\theta_{\pm}}F[q_{0}e^{i\theta_{\pm}}]\right\},\qquad\theta_{\pm T}=\mp\frac{1}{q_{0}}\mathop{\rm Re}\nolimits\left\{e^{-i\theta_{\pm}}F[q_{0}e^{i\theta_{\pm}}]\right\}.
\end{equation} If the perturbation is such that 
$F[q_{0}e^{i\theta_{\pm}}]=F[q_{0}]e^{i\theta_{\pm}}$, then we have
\begin{equation}
q_{0T}=\mathop{\rm Im}\nolimits F[q_{0}],\qquad (\Delta\theta)_{T}=0,
\end{equation}where 
$\Delta\theta=\theta_{+}-\theta_{-}$ is the asymptotic phase difference. The first formula in (5.3) can be used to determine the evolution of the background amplitude of a perturbed dark soliton a priori. In the method outlined below, one instead determines the parameters 
$k_{1}$ and 
$\Lambda_{1}$ as in (4.8) and recovers the background through 
$q_{0}=\sqrt{k_{1}^{2}+\Lambda_{1}^{2}}$, which necessarily will give the same result as (5.3).
When 
$\varepsilon=0$, the exact dark soliton solution satisfying the boundary conditions in (2.1) is 
$q(x,t)=u(\xi)$ given by (4.8). Now, suppose that the soliton parameters depend on the slow time scale:
\begin{equation}
k_{1}=k_{1}(T),\qquad \Lambda_{1}=\Lambda_{1}(T),\qquad x_{1}=x_{1}(T).
\end{equation} In this case, we redefine the travelling coordinate 
$\xi$ as
\begin{equation}
\xi=\Lambda_{1}(T)\left(x+2\int_{0}^{t}k_{1}\,ds-x_{1}(T)\right).
\end{equation} Note that since the slow-time-dependent velocity must be integrated, it is actually the case that the dependence of 
$k_{1}$ on an even slower time scale 
$\varepsilon^2 t$ could contribute nontrivially to the soliton centre. Thus, to ensure a correct prediction for the soliton centre, one would need to include another time scale in the perturbation theory. This is discussed in the context of a specific example in Remark 3 of Section 7.3, but is beyond the scope of this work.
Furthermore, to account for slow evolution of the overall phase, we introduce another parameter
\begin{equation}
\sigma_{1}=\sigma_{1}(T),\qquad \sigma_{1}(0)=0,
\end{equation}and let
\begin{equation}
q(x,t)\approx e^{i\sigma_{1}}\big[u(\xi)+\varepsilon\tilde{q}(\xi,t)\big],
\end{equation}where 
$\tilde{q}$ is the first-order correction and explicit dependence on 
$T$ is omitted for brevity. We will later show that the slowly evolving parameter 
$\sigma_{1}(T)$ corresponds to inevitable secular growth in the first-order correction, which can be absorbed into the phase at least for small times. Substituting (5.7) into (5.1) leads at 
$\mathcal{O}(\varepsilon)$ to the following equation:
\begin{equation}
(i\partial_{t}+L)A(\xi,t)=W(\xi),\qquad A=\begin{bmatrix}
\tilde q\\\tilde{q}^{\,*}
\end{bmatrix},\qquad W=\begin{bmatrix}
w[u]\\
-w[u]^*
\end{bmatrix},
\end{equation}
\begin{equation*}
w[u]=e^{-i\sigma_{1}}F[ue^{i\sigma_{1}}]-ie^{-i\sigma_{1}}\partial_{T}(ue^{i\sigma_{1}}),
\end{equation*}
$L$ is the same as in (4.10). Note that
\begin{align}
e^{-i\sigma_{1}}\partial_{T}(ue^{i\sigma_{1}})&=k_{1T}-\sigma_{1T}\Lambda_{1}\tanh\xi+i\Lambda_{1T}(\tanh\xi+\xi \mathop{\rm sech}\nolimits^{2}\xi) -ix_{1T}\Lambda_{1}^{2}\mathop{\rm sech}\nolimits^2\xi+i\sigma_{1T}k_{1}. \end{align} Throughout, the subscript 
$T$ denotes the partial derivative with respect to 
$T=\varepsilon t$. According to the fully reduced 1-soliton closure relation (4.33), we can expand 
$A$ and 
$W$ in terms of the complete set of squared eigenfunctions 
$\{\mathcal{Z}(\xi,\zeta),\mathcal{Z}_1(\xi),\tilde{\mathcal{Z}}_1(\xi)\}$ with respective adjoint eigenfunctions 
$\{\Upsilon(\xi,\zeta),\Upsilon_1(\xi),\tilde\Upsilon_1(\xi)\}$ as follows:
\begin{equation}
W(\xi)=\frac{\zeta_{1}}{8\Lambda_{1}^{3}}\left[\langle\tilde\Upsilon_{1},W\rangle\mathcal{Z}_1(\xi)+\langle\Upsilon_1,W\rangle\tilde{\mathcal{Z}}_1(\xi)\right]+\frac{1}{2\pi}{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{\langle\Upsilon(\zeta),W\rangle}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\mathcal{Z}(\xi,\zeta)d\zeta,
\end{equation}
\begin{equation}
A(\xi,t)=\frac{\zeta_{1}}{8\Lambda_{1}^{3}}\left[\langle\tilde\Upsilon_{1},A(t)\rangle\mathcal{Z}_1(\xi)+\langle\Upsilon_1,A(t)\rangle\tilde{\mathcal{Z}}_1(\xi)\right]+\frac{1}{2\pi}{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{\langle\Upsilon(\zeta),A(t)\rangle}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\mathcal{Z}(\xi,\zeta)d\zeta.
\end{equation}
\begin{equation}
i\partial_{t}\langle\Upsilon_1,A(t)\rangle=\langle\Upsilon_1,W\rangle,\qquad
i\partial_{t}\langle\tilde\Upsilon_{1},A(t)\rangle-4\Lambda_{1}^{2}\zeta_{1}^{-1}\langle\Upsilon_1,A(t)\rangle =\langle\tilde\Upsilon_1,W\rangle.
\end{equation}To suppress secular growth in (5.11), we need to enforce the orthogonality conditions
\begin{equation}
\langle\Upsilon_1,W\rangle=\langle\tilde\Upsilon_1,W\rangle=0.
\end{equation}The first orthogonality condition above yields the evolution of the soliton velocity explicitly in terms of the perturbation:
\begin{equation}
k_{1T}=\frac{1}{2}\,\text{Im}\int_{-\infty}^{\infty}F_{0}\mathop{\rm sech}\nolimits^{2}\xi\,d\xi, \qquad F_{0}=e^{-i\sigma_{1}}F[ue^{i\sigma_{1}}].
\end{equation} As we will show with several examples in Section 7, once the perturbation is specified, the second orthogonality condition in (5.12) yields two separate conditions on the soliton parameters (one corresponding to eliminating the divergent integrals). One of these conditions can be used to fully determine the dependence of 
$\Lambda_{1}$ on the slow time 
$T$, and the other relates the evolution of the soliton phase and centre, 
$\sigma_{1}$ and 
$x_{1}$.
Moreover, equating coefficients of the continuous eigenfunctions gives
\begin{equation}
i\partial_{t}\langle\Upsilon(\zeta),A(t)\rangle+4\lambda(k-k_{1})\langle\Upsilon(\zeta),A(t)\rangle=\langle\Upsilon(\zeta),W\rangle.
\end{equation}The solution of (5.14) with zero initial condition is
\begin{equation}
\langle\Upsilon(\zeta),A(t)\rangle=\frac{\langle\Upsilon(\zeta),W\rangle}{4\lambda(k-k_{1})}\left[1-e^{4i\lambda(k-k_{1})t}\right].
\end{equation}Using the above equation in (5.10b) yields for the first-order correction:
\begin{equation}
A(\xi,t)=\frac{1}{2\pi}{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{1}{\gamma(\zeta)|\zeta-\zeta_{1}|^{4}}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{4\lambda(\zeta)\big(k(\zeta)-k_{1}\big)}\langle\Upsilon(\zeta),W\rangle\mathcal{Z}(\xi,\zeta)d\zeta.
\end{equation} Note that 
$\gamma(\zeta)=2\zeta^{-1}\lambda(\zeta)$, so this could be written as
\begin{equation}
A(\xi,t)=\frac{1}{2\pi}{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{\zeta}{|\zeta-\zeta_{1}|^{4}}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{8\lambda(\zeta)^2\big(k(\zeta)-k_{1}\big)}\langle\Upsilon(\zeta),W\rangle\mathcal{Z}(\xi,\zeta)d\zeta,
\end{equation}which, up to an inessential difference in the normalisation of the squared eigenfunctions, has the same form as the one given in [Reference Chen, Chen and Huang10], with the important caveat that the integral has to be considered in the principal value sense, due to the singularities at 
$\zeta=\pm q_0$ which arise from 
$\lambda(\zeta)$. We stress that due to the fact that the integrand has poles at 
$\pm q_0$, it is crucial for this integral to be considered in the principal value sense in order for it to be well-defined. This also explains why one has to explicitly single out the contributions at 
$\pm q_0$ from the closure relation.
It is worth highlighting that the issue of the poles at 
$\pm q_{0}$ was explicitly acknowledged in [Reference Konotop and Vekslerchik39]. Our present approach to account for these poles is to modify the adjoint discrete squared eigenfunctions appropriately (i.e., enforce orthogonality with 
$\tilde\Upsilon_{1}$ rather than with 
$\Upsilon'_{1}$) such that the first-order correction is reduced to a principal value integral. On the other hand, the approach of [Reference Konotop and Vekslerchik39] was to enforce the additional orthogonality condition 
$\langle\Upsilon_{0}^{\pm},W\rangle=0$. However, this condition is too strong. In fact, we will later show that the inner product 
$\langle\Upsilon_{0}^{\pm},W\rangle$ is generically nonzero, a fact that is crucial to understanding the generation of the shelf.
6. First-order correction, shelf, and evolution of the phase
As we will show in this section, the first-order correction can be effectively used to describe the shelf that develops on each side of the soliton, which in fact results from the singularities of the scattering data at the points 
$\pm q_0$, and to obtain estimates for the height and velocity of the shelf. Furthermore, the asymptotic behaviour of the first-order correction provides the spatial phase gradient on each side of the shelf, as well as a formula for the slow-time evolution of the core soliton phase, which in turn allows one to determine the dependence of the soliton centre on the slow time 
$T$.
Taking the first component of (5.17), on account of (5.7) and (5.8a) the first-order correction 
$\tilde{q}$ is given by:
\begin{align} \tilde{q}(\xi,t)&=\frac{1}{2\pi}{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{\zeta}{|\zeta-\zeta_{1}|^{4}}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{8\lambda(\zeta)^2\big(k(\zeta)-k_{1}\big)}\langle\Upsilon(\zeta),W\rangle \Big[i(\zeta_{1}-q_{0}^{2}\zeta^{-1})+\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big]^2 e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta. \end{align} First, we study the first-order correction on the right side of the shelf, i.e., the asymptotic limit 
$\xi\rightarrow+\infty$:
\begin{equation} \tilde{q}^{+}\sim-\frac{1}{2\pi}{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{\zeta}{|\zeta-\zeta_{1}|^{4}}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{8\lambda(\zeta)^2\big(k(\zeta)-k_{1}\big)}\langle\Upsilon(\zeta),W\rangle\big(\zeta_{1}-q_{0}^{2}\zeta^{-1}\big)^2 e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta.
\end{equation} For large 
$\xi$, the 
$\xi$-dependent exponential in (6.2) is rapidly oscillating, so the primary contributions to the principal value integral come from neighbourhoods around the simple poles 
$\zeta=\pm q_{0}$ (recall 
$\lambda(\pm q_{0})=0$). It will be convenient to define
\begin{equation}
\mathcal{F}(\zeta)=-\frac{\zeta(\zeta_{1}-q_{0}^{2}\zeta^{-1})^{2}}{16\pi|\zeta-\zeta_{1}|^{4}\big(k(\zeta)-k_{1}\big)}\langle\Upsilon(\zeta),W\rangle,
\end{equation}and split the integral into two parts
\begin{equation}
\tilde{q}^{+}\sim \mathcal{F}(-q_{0})\underbrace{{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}0}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{\lambda(\zeta)^2}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta}_{I_{-}}+\mathcal{F}(q_{0})\underbrace{{\int{\kern-11.2pt}-}_{{\kern-3.5pt}0}^{{\kern2.2pt}\infty}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{\lambda(\zeta)^2}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta}_{I_{+}}.
\end{equation}After simplification, including using (4.27), we have
\begin{equation}
\mathcal{F}(\pm q_{0})=\pm\frac{i\alpha e^{i\theta}}{16\pi(k_{1}\mp q_{0})},
\end{equation}where 
$\alpha$ is defined by
\begin{equation}
\left\langle\begin{bmatrix}u^{*}\\u\end{bmatrix},W\right\rangle=:i\alpha,\qquad\alpha\in\mathbb{R}.
\end{equation} Note that due to the form of 
$W$, the inner product in (6.6) is purely imaginary regardless of the choice of perturbation. The integrals 
$I_\pm$ in (6.4) are computed explicitly through lengthy but straightforward calculations detailed in Appendix C. This calculation also naturally reveals the right boundary of the shelf to be 
$\xi\sim2\Lambda_{1}(k_{1}+q_{0})t$. In particular, putting everything together into (6.4), we find that 
$\tilde{q}^{+}\sim0$ for 
$\xi\gg2\Lambda_{1}(k_{1}+q_{0})t$, and inside the shelf region 
$1\ll\xi\ll2\Lambda_{1}(k_{1}+q_{0})t$ we have
\begin{equation}
\tilde{q}^{+}\sim\frac{\alpha e^{i\theta}}{8q_{0}(k_{1}+q_{0})}+\frac{\alpha e^{i\theta}}{4\Lambda_{1}(k_{1}+q_{0})}i\xi+\frac{\alpha e^{i\theta}}{2}it.
\end{equation} This suggests that there is inevitable growth in 
$\xi$ and 
$t$ in the first-order correction resulting from the singularities at 
$\pm q_{0}$. Now, according to our original ansatz (5.7), the full approximate solution in the right shelf region should be expressed as
\begin{equation}
q^{+}\sim e^{i\sigma_{1}}\left[q_{0}e^{i\theta}+\varepsilon\tilde{q}^{+}\right].
\end{equation} Substituting the expression for 
$\tilde{q}^{+}$, (6.8) is equivalent up to 
$\mathcal{O}(\varepsilon)$ to
\begin{equation}
q^{+}\sim e^{i\theta+i\sigma_{1}}\left[q_{0}+\varepsilon\frac{\alpha}{8q_{0}(k_{1}+q_{0})}\right]\left[1+\frac{\alpha}{4q_{0}\Lambda_{1}(k_{1}+q_{0})}i\varepsilon\xi+\frac{\alpha}{2q_{{0}}}i\varepsilon t\right].
\end{equation} We can account for the terms that grow with 
$T=\varepsilon t$ and 
$\varepsilon\xi$ by moving them into the exponent, i.e., by approximating the above as
\begin{equation}
q^{+}\sim \left[q_{0}+\frac{\alpha}{8q_{0}(k_{1}+q_{0})}\right]\exp\left\{i\theta+i\sigma_{1}+\frac{\alpha}{4q_{0}\Lambda_{1}(k_{1}+q_{0})}i\varepsilon\xi+\frac{\alpha}{2q_{{0}}}i\varepsilon t\right\}.
\end{equation} The purpose of the introduction of the 
$T$-dependent phase parameter 
$\sigma_{1}$ is now evident. Namely, we can choose the slow-time evolution of the soliton phase to be
\begin{equation}
\sigma_{1T}=-\frac{\alpha}{2q_{0}},
\end{equation}such that the secular growth in time resulting from the inevitable time dependence of the first-order correction integral is suppressed. Therefore, we construct the full approximate solution in the right shelf region as
\begin{equation}
q^{+}\sim (q_{0}+\varepsilon h^{+})e^{i\theta+i\varepsilon\varphi^+},
\end{equation}where the height of the shelf and the spatial phase gradient are, respectively, given by
\begin{equation}
h^{+}=\frac{\alpha}{8q_{0}(k_{1}+q_{0})},\qquad\varphi_{\xi}^{+}=\frac{\alpha}{4q_{0}\Lambda_{1}(k_{1}+q_{0})}.
\end{equation} We later empirically demonstrate for several examples that this approximation of the spatial and temporal growth of the first-order correction as a complex exponential does accurately reflect the dynamics, at least for a moderate time interval. That said, we hypothesise that the breakdown of this approximation could be a reason why our predictions for the soliton phase, which coincide with the ones in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], do not necessarily retain their validity up to 
$\mathcal{O}(1/\varepsilon)$, as we will discuss in more details in Section 7 and in the conclusions.
Note that (6.11) serves as a generic formula for the evolution of the soliton phase parameter, which depends (through 
$\alpha$, see (6.6)) on the form of the perturbation. Recall that the orthogonality conditions from Section 5 provide three differential equations for four parameters. When supplemented with this formula, the soliton parameters are completely determined. Also, observe from the above that the formation of the shelf is directly connected to the evolution of 
$\sigma_{1}$. Namely, if for a given perturbation we have 
$\alpha=0$, then no shelf develops and 
$\sigma_{1}$ is constant (see, e.g., the self-steepening perturbation discussed in Section 7.4).
Repeating these calculations on the left side of the shelf, i.e., in the asymptotic limit 
$\xi\rightarrow-\infty$, we find that the left boundary of the shelf is 
$\xi\sim2\Lambda_{1}(k_{1}-q_{0})t$, with 
$\tilde{q}^{-}\sim0$ for 
$\xi\ll2\Lambda_{1}(k_{1}-q_{0})t$, and in the left shelf region 
$2\Lambda_{1}(k_{1}-q_{0})t\ll\xi\ll-1$:
\begin{equation}
\tilde{q}^{-}\sim-\frac{\alpha e^{-i\theta}}{8q_{0}(k_{1}-q_{0})}+\frac{\alpha e^{-i\theta}}{4\Lambda_{1}(k_{1}-q_{0})}i\xi+\frac{\alpha e^{-i\theta}}{2}it.
\end{equation}Similarly, we construct the full approximate solution in the shelf region as
\begin{equation}
q^{-}\sim (q_{0}+\varepsilon h^{-})e^{i\theta+i\varepsilon\varphi^-},
\end{equation}where
\begin{equation}
h^{-}=-\frac{\alpha}{8q_{0}(k_{1}-q_{0})},\qquad \varphi_{\xi}^{-}=\frac{\alpha}{4q_{0}\Lambda_{1}(k_{1}-q_{0})}.
\end{equation}From these formulae, it can be seen that, regardless of the perturbation, we have
\begin{equation}
\Lambda_{1}\varphi_{\xi}^{\pm}=\pm 2h^{\pm},
\end{equation}which is in agreement with the results in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1].
7. Applications
In this section, we apply our results to the perturbed NLS equation (5.1) with various small forcing perturbations: linear damping, nonlinear damping (also referred to as two-photon absorption, or TPA, in nonlinear optics, which is the terminology used in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1]), dissipation and self-steepening. All the results are corroborated by direct numerical simulations, and compared with the results of the direct perturbation theory, and with the earlier works using perturbation theory based on the squared eigenfunctions.
Throughout the following subsections, we compare our predictions with numerical simulations of the perturbed defocusing NLS equation on a nonzero background. The numerical computations were performed using the method proposed by Kassam and Trefethen in [Reference Kassam and Trefethen28], which is a modification of the fourth-order exponential time differencing (ETDRK4) method (see also [Reference Cox and Matthews12, Reference Trefethen49]) that employs contour integration to avoid numerical instabilities characteristic of stiff PDEs. The spatial discretisation is done using a Fourier integral representation. As such, the (non-decaying) initial condition must be multiplied by a rapidly decaying bump function to justify the use of the Fourier basis. To ensure that any disturbances resulting from this truncation remain sufficiently far from the dark soliton, we use a large computational domain with the length scale of the bump function being orders of magnitude larger than the length scale of the perturbed dark soliton and its shelf. For a typical run, the width of the bump function was taken to be on the order of 
$10^{4}$. The time step and grid spacing were both taken to be 
$\sim10^{-3}$.
7.1. Linear damping
Consider Eq. (5.1) with a linear damping perturbation:
\begin{equation}
F_{0}=-iu=-ik_{1}+\Lambda_{1}\tanh\xi.
\end{equation}In this case, (5.13) implies that
\begin{equation}
k_{1T}=-k_{1}.
\end{equation}We now demonstrate that two distinct constraints on the soliton parameters arise from the second orthogonality condition. To this end, a straightforward calculation shows that
\begin{align}\langle\Upsilon_{1}',W\rangle &=-\frac{4q_{0}}{\zeta_{1}}\int_{-\infty}^{\infty}\Big\{\left[\mathop{\rm Im}\nolimits(e^{-i\theta}F_{0})-\mathop{\rm Re}\nolimits(e^{-i\theta}D)\right]e^{\xi}\mathop{\rm sech}\nolimits\xi +\frac{k_{1}}{q_{0}}\left[\mathop{\rm Im}\nolimits(F_{0})-\mathop{\rm Re}\nolimits(D)\right]\xi\mathop{\rm sech}\nolimits^{2}\xi\Big\}d\xi, \end{align}where 
$D=e^{-i\sigma_{1}}\partial_{T}(ue^{i\sigma_{1}})$, which is given explicitly in (5.9). Note that
\begin{equation}
\mathop{\rm Re}\nolimits(D)=k_{1T}-\sigma_{1T}\Lambda_{1}\tanh\xi,
\end{equation}
\begin{align}
\mathop{\rm Re}\nolimits(e^{-i\theta}D)&=\frac{1}{q_{0}}\Big(k_{1}k_{1T}-\sigma_{1T}\Lambda_{1}k_{1}(\tanh\xi-1) +\,\Lambda_{1}\Lambda_{1T}(\tanh\xi+\xi\mathop{\rm sech}\nolimits^{2}\xi)-x_{1T}\Lambda_{1}^{3}\mathop{\rm sech}\nolimits^{2}\xi\Big),
\end{align} where we have used 
$\cos\theta=k_{1}/q_{0}$ and 
$\sin\theta=\Lambda_{1}/q_{0}$, and in the linear damping case we have
\begin{equation}
\mathop{\rm Im}\nolimits(F_{0})=-k_{1},\qquad
\mathop{\rm Im}\nolimits(e^{-i\theta}F_{0})=-\frac{1}{q_{0}}\Big(k_{1}^{2}+\Lambda_{1}^{2}\tanh\xi\Big).
\end{equation}
\begin{equation}
\langle{\Upsilon}'_{1},W\rangle = 4\Lambda_{1}\zeta_{1}^{-1}\left\{(\Lambda_{1T}+\Lambda_{1})\int_{-\infty}^{\infty}e^{\xi}\tanh\xi\mathop{\rm sech}\nolimits\xi \,d\xi+\sigma_{1T}k_{1}+\Lambda_{1T}-2x_{1T}\Lambda_{1}^{2}\right\}.
\end{equation}The integral that remains is divergent. Next, we find
\begin{equation}
\langle(u^{*},u)^{T},W\rangle
=-\frac{2i}{\Lambda_{1}}\int_{-\infty}^{\infty}\Big\{|u|^{2}+\mathop{\rm Re}\nolimits(u^{*}D)\Big\}d\xi,
\end{equation}with
\begin{equation}
|u|^{2}=k_{1}^{2}+\Lambda_{1}^{2}\tanh^{2}\xi,
\end{equation}
\begin{equation}
\mathop{\rm Re}\nolimits(u^{*}D)=k_{1}k_{1T}+\Lambda_{1}\Lambda_{1T}\tanh\xi(\tanh\xi+\xi\mathop{\rm sech}\nolimits^{2}\xi)-x_{1T}\Lambda_{1}^{3}\tanh\xi\mathop{\rm sech}\nolimits^{2}\xi.
\end{equation}
\begin{equation}
\langle(u^{*},u)^{T},W\rangle=-2i\left\{(\Lambda_{1T}+\Lambda_{1})\int_{-\infty}^{\infty}\tanh^{2}\xi \,d\xi+\Lambda_{1T}\right\}.
\end{equation}Again, the remaining integral is divergent. To eliminate the divergence from both inner products above, it suffices to enforce:
\begin{equation}
\Lambda_{1T}=-\Lambda_{1}.
\end{equation} With 
$\Lambda_{1}$ now fully determined, (7.5) and (7.8) reduce to
\begin{equation}
\langle{\Upsilon}'_{1},W\rangle=4\Lambda_{1}\zeta_{1}^{-1}(\sigma_{1T}k_{1}-\Lambda_{1}-2x_{1T}\Lambda_{1}^{2}),
\end{equation}
\begin{equation}
\langle(u^{*},u)^{T},W\rangle=2i\Lambda_{1}.
\end{equation}
$\sigma_{1T}$, was not considered; (ii) the discrete squared eigenfunctions were not properly modified to account for the contributions from 
$\pm q_{0}$, that is, the orthogonality condition being enforced was 
$\langle\Upsilon'_{1},W\rangle=0$. For these reasons, the conclusion in [Reference Chen, Chen and Huang10] is that the evolution of the soliton centre is given by 
$x_{1T}=-1/2\Lambda_{1}$. This prediction is not consistent with physical intuition (or numerical simulations), as it is independent of the direction that the soliton is travelling. Putting (7.10a) and (7.10b) together as in (4.32), we find the proper orthogonality condition to be
\begin{equation}
\langle\tilde\Upsilon_{1},W\rangle=4\Lambda_{1}\zeta_{1}^{-1}(\sigma_{1T}k_{1}-2x_{1T}\Lambda_{1}^{2})=0,
\end{equation}which provides another condition on the soliton parameters, namely,
\begin{equation}
\sigma_{1T}k_{1}=2x_{1T}\Lambda_{1}^{2}.
\end{equation} Finally, upon supplementing this condition with the formula for 
$\sigma_{1T}$ obtained in (6.11), all four soliton parameters are uniquely determined. In particular, using (7.10b), we find
\begin{equation}
\sigma_{1T}=-\frac{\Lambda_{1}}{q_{0}}.
\end{equation}The corresponding evolution of the soliton centre is
\begin{equation}
x_{1T}=-\frac{k_{1}}{2q_{0}\Lambda_{1}}.
\end{equation}Solving for the explicit time evolution of all soliton parameters:
\begin{equation}
q_{0}(T)=q_{0}(0)e^{-T},\qquad
k_{1}(T)=k_{1}(0)e^{-T},\qquad
\Lambda_{1}(T)=\Lambda_{1}(0)e^{-T},
\end{equation}
\begin{equation}
\displaystyle\sigma_{1}(T)=-\frac{\Lambda_{1}(0)}{q_{0}(0)}T+\sigma_{1}(0),\qquad
x_{1}(T)=\frac{k_{1}(0)}{2q_{0}(0)\Lambda_{1}(0)}(1-e^{T})+x_{1}(0).
\end{equation}
$\mathcal{O}(\varepsilon^{2})$. Figure 1 shows a comparison of the modulus and phase of the predicted solution with a numerical simulation for 
$\varepsilon=0.02$ with snapshots taken at moderate times.Top Row: Comparison of the modulus (left) and phase (right) of the predicted (dashed red) and numerical (solid black) solutions for a dark soliton under the influence of the linear damping perturbation 
$F[q]=-iq$ with 
$\varepsilon=0.02$ at time 
$t=10$. The initial soliton parameters are 
$q_{0}(0)=2$, 
$k_{1}(0)=-1/2$, 
$x_{1}(0)=\sigma_{1}(0)=0$. Bottom Row: The same at 
$t=20$, in which a slight deviation in our prediction for the phase from the numerics is already visible.
Note that here and in all subsequent examples, since the background evolves on the slow time scale, for the sake of comparison, we numerically solve the version of the NLS equation without the background factored out (i.e., (1.1) without the 
$q_{0}^{2}$ term) and express the approximate solution as
\begin{equation}
q(x,t)\approx\left\{k_{1}(T)+i\Lambda_{1}(T)\tanh\left[\Lambda_{1}(T)\left(x+2\int_{0}^{t}k_{1}\,ds-x_{1}(T)\right)\right]\right\}e^{i\sigma_{1}(T)-2i\int_{0}^{t}q_{0}^{2}\,ds}.
\end{equation} Figure 2 shows our predicted evolution of the soliton centre and phase with 
$t$, compared to numerical measurements of the corresponding quantities.
The predicted evolution of the parameters 
$\sigma_{1}$ (left) and 
$x_{1}$ (right) with time (dashed red), compared with numerical measurements (dotted black). The parameters are the same as in Figure 1. Note that both predictions are accurate for a moderate time interval, but start to deviate as time increases. The blue dashed line in the right panel shows the prediction of the method of [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] for the soliton centre, which is seen to be more accurate for long time. A similar deviation in the prediction for the soliton phase for larger times can be seen in Fig. 7 in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1].
From this figure it is clear that, as mentioned in the introduction, the predictions for the centre and phase begin to deviate from the numerics before 
$t=\mathcal{O}(1/\varepsilon)$. We postulate that this breakdown in the approximation could be due to the secular growth in the first-order correction integral, and that one would likely be required to compute the 
$\mathcal{O}(\varepsilon^{2})$ term in the expansion to extend the time interval of validity. Note that the prediction for the soliton centre presented in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], which is given in terms of a second-order differential equation found using information at 
$\mathcal{O}(\varepsilon^{2})$, does remain accurate for longer time intervals (though the same discrepancy in the phase is present in that work). It is worth noting that in Figure 2, and in some subsequent figures, one can perceive a slight staggering of the points corresponding to the numerical measurements of the centre and phase. This is due to the fact that the centre could lie between two grid points (the spacing of which was constrained in practice due to the necessity for a very large computational domain), and is not an error in the simulation itself. The error is on the order of the grid spacing, and does not affect the overall trend.
Lastly, the predicted height and phase gradient of the shelf in the linear damping case are given by
\begin{equation}
h^{\pm}=\pm\frac{\Lambda_{1}}{4q_{0}(k_{1}\pm q_{0})},\qquad\varphi_{\xi}^{\pm}=\frac{1}{2q_{0}(k_{1}\pm q_{0})}.
\end{equation} From these formulae, it can be seen that since 
$|k_{1}| \lt q_{0}$ the shelf is raised on both sides of the soliton, and the phase gradient is positive on the right side and negative on the left side of the soliton. In Figure 3, we show a space-time plot from which the development of the shelf around the soliton and its propagation can be seen. Figure 4 shows a snapshot of the shelf region, along with the predictions for the height and phase gradient.
A space-time plot showing the development of a raised shelf around the soliton under linear damping. The blue dashed lines denote the boundaries of the shelf region, and the red dashed line corresponds to the predicted path of the soliton core, accounting for the slow evolution of the velocity. For visualisation purposes, the evolution of the background has been artificially removed. In the right panel, a comparison of the predicted decrease in the background 
$q_{0}$ (top) and trough amplitude 
$|k_{1}|$ (bottom) with numerical measurements is shown up to 
$t=1/\varepsilon=50$. The initial parameters are the same as in Figure 1.
The modulus (left) and phase (right) of a numerical simulation of a dark soliton under the influence of the linear damping perturbation 
$F[q]=-iq$ with 
$\varepsilon=0.02$ at 
$t=10$. The blue dashed lines denote the predicted boundaries of the shelf region. The red (resp., green) dashed lines represent the predictions for the shelf height and phase gradient on the right (resp., on the left) of the soliton. The initial parameters are the same as in Figure 1.
7.2. Nonlinear damping
Consider Eq. (5.1) with a nonlinear damping perturbation
\begin{equation}
F_{0}=-i|u|^{2}u=(q_{0}^{2}-\Lambda_{1}^{2}\mathop{\rm sech}\nolimits^{2}\xi)(-ik_{1}+\Lambda_{1}\tanh\xi).
\end{equation}In this case, (5.13) implies that
\begin{equation}
k_{1T}=-\left(\frac{2}{3}k_{1}^{2}+\frac{1}{3}q_{0}^{2}\right)k_{1}.
\end{equation}Following similar steps as in Section 7.1, we substitute
\begin{equation}
\mathop{\rm Im}\nolimits(e^{-i\theta}F_{0})=-\frac{1}{q_{0}}(q_{0}^{2}-\Lambda_{1}^{2}\mathop{\rm sech}\nolimits^{2}\xi)(k_{1}^{2}+\Lambda_{1}^{2}\tanh\xi),
\end{equation}into (7.1), and after simplification we find
\begin{align}
\langle\Upsilon_{1}',W\rangle &= \frac{4\Lambda_{1}}{\zeta_{1}}\Bigg\{\left(\Lambda_{1}q_{0}^{2}+\frac{2}{3}\Lambda_{1}k_{1}^{2}+\Lambda_{1T}\right)\int_{-\infty}^{\infty}e^{\xi}\tanh\xi\mathop{\rm sech}\nolimits\xi \,d\xi +\,\Lambda_{1T}-\frac{2}{3}\Lambda_{1}q_{0}^{2}+\sigma_{1T}k_{1}-2x_{1T}\Lambda_{1}^{2}\Bigg\}. \end{align}Furthermore, one can verify that
\begin{equation}
\langle(u^{*},u)^{T},W\rangle=-2i\left\{\left(\Lambda_{1}q_{0}^{2}+\frac{2}{3}\Lambda_{1}k_{1}^{2}+\Lambda_{1T}\right)\int_{-\infty}^{\infty}1d\xi-4q_{0}^{2}\Lambda_{1}+\frac{4}{3}\Lambda_{1}^{3}-\Lambda_{1T}\right\}.
\end{equation}To eliminate the divergence in both (7.21) and (7.22), set
\begin{equation}
\Lambda_{1T}=-\left(\frac{2}{3}k_{1}^{2}+q_{0}^{2}\right)\Lambda_{1}.
\end{equation} Note that combining the evolution of 
$k_{1}$ and 
$\Lambda_{1}$, we get the evolution of the background
\begin{equation}
q_{0T}=-q_{0}^{3}.
\end{equation} With 
$\Lambda_{1}$ determined, (7.21) and (7.22) become
\begin{equation}
\langle\Upsilon_{1}',W\rangle= \frac{4\Lambda_{1}}{\zeta_{1}}\{-\left(\frac{2}{3}k_{1}^{2}+\frac{5}{3}q_{0}^{2}\right)\Lambda_{1}+\sigma_{1T}k_{1}-2x_{1T}\Lambda_{1}^{2}\},
\end{equation}
\begin{equation}
\langle(u^{*},u)^{T},W\rangle=2i\left(\frac{2}{3}k_{1}^{2}+\frac{5}{3}q_{0}^{2}\right)\Lambda_{1}.
\end{equation}
\begin{equation}
\langle\tilde\Upsilon_{1},W\rangle=4\Lambda_{1}\zeta_{1}^{-1}(\sigma_{1T}k_{1}-2x_{1T}\Lambda_{1}^{2}),
\end{equation}which is the same as in the previous example. So, we must again enforce the condition
\begin{equation}
\sigma_{1T}k_{1}=2x_{1T}\Lambda_{1}^{2}.
\end{equation}From (6.11), using (7.25b) we find the evolution of the phase to be
\begin{equation}
\sigma_{1T}=-\left(\frac{2}{3}k_{1}^{2}+\frac{5}{3}q_{0}^{2}\right)\frac{\Lambda_{1}}{q_{0}},
\end{equation}which in turn determines the evolution of the soliton centre to be
\begin{equation}
x_{1T}=-\left(\frac{2}{3}k_{1}^{2}+\frac{5}{3}q_{0}^{2}\right)\frac{k_{1}}{2q_{0}\Lambda_{1}}.
\end{equation}The evolution equations for the velocity, amplitude and phase coincide with the ones obtained from the method of [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] (note that no prediction for the soliton centre is given in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] for this example).
In this case, the system of nonlinear differential equations for the soliton parameters must be solved numerically. Figure 5 shows snapshots of a comparison between our predictions (with the evolution equations for the soliton parameters solved numerically and inserted into (7.16)) with a simulation. Qualitatively, the nonlinear damping case resembles the linear damping case, with the prediction for the phase showing a noticeable discrepancy from the numerics before 
$t=\mathcal{O}(1/\varepsilon)$. While the prediction for the centre holds up reasonably well, a deviation does occur, illustrated in Figure 6. Finally, Figure 7 shows a spacetime plot of the development of the shelf, and Figure 8 shows the shelf region and the predictions for the height and phase gradient, which in this case are given by
\begin{equation}
h^{\pm}=\pm\frac{\Lambda_{1}}{4q_{0}(k_{1}\pm q_{0})}\left(\frac{2}{3}k_{1}^{2}+\frac{5}{3}q_{0}^{2}\right),\qquad\varphi_{\xi}^{\pm}=\frac{1}{2q_{0}(k_{1}\pm q_{0})}\left(\frac{2}{3}k_{1}^{2}+\frac{5}{3}q_{0}^{2}\right).
\end{equation}Like in the linear damping case, the shelf is raised on both sides of the soliton.
Top Row: Comparison of the modulus (left) and phase (right) of the predicted (dashed red) and numerical (solid black) solutions for a dark soliton under the influence of the nonlinear damping perturbation 
$F[q]=-i|q|^{2}q$ with 
$\varepsilon=0.02$ at time 
$t=5$. The initial soliton parameters are 
$q_{0}(0)=2$, 
$k_{1}(0)=-1/2$, 
$x_{1}(0)=\sigma_{1}(0)=0$. Bottom Row: The same at 
$t=10$, from which a deviation in the phase can be seen.
The predicted evolution of the parameters 
$\sigma_{1}$ (left) and 
$x_{1}$ (right) with time (dashed red), compared with numerical measurements (dotted black). The initial parameters are the same as in Figure 5.
A space-time plot showing the development of a raised shelf around the soliton under nonlinear damping. The blue dashed lines denote the boundaries of the shelf region, and the red dashed line corresponds to the predicted path of the soliton core, accounting for the slow evolution of the velocity. For visualisation purposes, the evolution of the background has been artificially removed. In the right panel, a comparison of the predicted decrease in the background 
$q_{0}$ (top) and trough amplitude 
$|k_{1}|$ (bottom) with numerical measurements is shown up to 
$t=1/\varepsilon=50$. The initial parameters are the same as in Figure 5.
The modulus (left) and phase (right) of a numerical simulation of a dark soliton under the influence of the nonlinear damping perturbation 
$F[q]=-i|q|^{2}q$ with 
$\varepsilon=0.02$ at 
$t=5$. The blue dashed lines denote the predicted boundaries of the shelf region. The red (resp., green) dashed lines represent the predictions for the shelf height and phase gradient on the right (resp., on the left) of the soliton. The initial parameters are the same as in Figure 5.
7.3. Dissipation
Consider Eq. (5.1) with a dissipative perturbation
\begin{equation}
F_{0}=iu_{xx}=2\Lambda_{1}^{3}\mathop{\rm sech}\nolimits^{2}\xi\tanh\xi.
\end{equation} In this case, since 
$F_{0}$ is purely real, (5.13) implies that
\begin{equation}
k_{1T}=0.
\end{equation} Returning to (7.1) and putting in the new expression for 
$F_{0}$ gives
\begin{gather}\langle\Upsilon_{1}',W\rangle =\frac{4}{\zeta_{1}}\int_{-\infty}^{\infty}\Big\{2\Lambda_{1}^{4}e^{\xi}\mathop{\rm sech}\nolimits^{3}\xi\tanh\xi+\Lambda_{1}\Lambda_{1T}(e^{\xi}\tanh\xi\mathop{\rm sech}\nolimits\xi+\xi e^{\xi}\mathop{\rm sech}\nolimits^{3}\xi) \nonumber\\
-x_{1T}\Lambda_{1}^{3}e^{\xi}\mathop{\rm sech}\nolimits^{3}\xi
-\sigma_{1T}\Lambda_{1}k_{1}(\tanh\xi-1)e^{\xi}\mathop{\rm sech}\nolimits\xi-k_{1}\sigma_{1T}\Lambda_{1}\xi\tanh\xi\mathop{\rm sech}\nolimits^{2}\xi\Big\}d\xi.\end{gather}Computing the convergent integrals, this reduces to
\begin{equation}
\langle\Upsilon_{1}',W\rangle = 4\Lambda_{1}\zeta_{1}^{-1}\left\{\Lambda_{1T}\int_{-\infty}^{\infty}e^{\xi}\tanh\xi\mathop{\rm sech}\nolimits\xi \,d\xi+\frac{4}{3}\Lambda_{1}^{3}+\sigma_{1T}k_{1}+\Lambda_{1T}-2x_{1T}\Lambda_{1}^{2}\right\}.
\end{equation}Next, we find
\begin{equation}
\langle(u^{*},u)^{T},W\rangle=-2i\left\{\Lambda_{1T}\int_{-\infty}^{\infty}\tanh^{2}\xi \,d\xi+\frac{4}{3}\Lambda_{1}^{3}+\Lambda_{1T}\right\}.
\end{equation}The only way to remove the divergence from both (7.34) and (7.35) is to set
\begin{equation}
\Lambda_{1T}=0.
\end{equation}With this, we have
\begin{equation}
\langle\Upsilon_1^\prime,W\rangle= 4\Lambda_{1}\zeta_{1}^{-1}\left(\frac{4}{3}\Lambda_{1}^{3}+\sigma_{1T}k_{1}-2x_{1T}\Lambda_{1}^{2}\right),
\end{equation}
\begin{equation}
\langle(u^{*},u)^{T},W\rangle=-\frac{8}{3}i\Lambda_{1}^{3}.
\end{equation}
\begin{equation}
\langle\tilde{\Upsilon}_1,W\rangle=4\Lambda_{1}\zeta_{1}^{-1}(\sigma_{1T}k_{1}-2x_{1T}\Lambda_{1}^{2}),
\end{equation}which leads to the same condition on the centre and phase that was found in both the previous examples,
\begin{equation}
\sigma_{1T}k_{1}=2x_{1T}\Lambda_{1}^{2}.
\end{equation}Substituting (7.37b) into (6.11) gives
\begin{equation}
\sigma_{1T}=\frac{4\Lambda_{1}^{3}}{3q_{0}},
\end{equation}and the corresponding evolution of the soliton centre is
\begin{equation}
x_{1T}=\frac{2\Lambda_{1}k_{1}}{3q_{0}}.
\end{equation}Solving explicitly for all soliton parameters:
\begin{equation}
q_{0}(T)=q_{0}(0),\;\;\;k_{1}(T)=k_{1}(0),\;\;\;\Lambda_{1}(T)=\Lambda_{1}(0),
\end{equation}
\begin{equation}
\displaystyle\sigma_{1}(T)=\frac{4\Lambda_{1}(0)^{3}}{3q_{0}(0)}T+\sigma_{1}(0),\;\;\;x_{1}(T)=\frac{2\Lambda_{1}(0)k_{1}(0)}{3q_{0}(0)}T+x_{1}(0).
\end{equation}
$T$. In [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], a quadratic prediction that is valid for longer time intervals is obtained, but we again note that the method used in that work employs information at 
$\mathcal{O}(\varepsilon^{2})$ to determine the centre, suggesting that a higher order correction would be needed in the present method as well. The reason why our theory does not capture the quadratic term is explained in Remark 3 below.Top Row: Comparison of the modulus (left) and phase (right) of the predicted (dashed red) and numerical (solid black) solutions for a dark soliton under the influence of the dissipative perturbation 
$F[q]=iq_{xx}$ with 
$\varepsilon=0.02$ at time 
$t=10$. The initial soliton parameters are 
$q_{0}(0)=2$, 
$k_{1}(0)=-1/2$, 
$x_{1}(0)=\sigma_{1}(0)=0$. Bottom Row: The same at 
$t=20$, from which it can be seen that our prediction for the soliton centre deviates by an amount consistent with a comparison with [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], which suggests that the next correction to the centre would be quadratic in 
$T$. Note that the shift in the phase in the bottom right panel is due to the shift in the centre; aside from this, the phase prediction remains accurate, as can be seen from the left panel in Fig. 10.
The discrepancy is illustrated in Figure 10, where the predicted phase (which does remain accurate for long times) and centre are compared to numerical measurements. Note that our linear prediction for the centre captures the initial slope of the curve.
The predicted evolution of the parameters 
$\sigma_{1}$ (left) and 
$x_{1}$ (right) with time (dashed red), compared with numerical measurements (dotted black) under a dissipative perturbation. The parameters are the same as in Figure 9. Note that the prediction for the phase remains accurate for long times, but the linear prediction for the centre only remains valid for a short time interval. The blue dashed line in the right panel shows the quadratic prediction for the soliton centre obtained using the method of [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], which accurately captures the curvature.
Remark 3. In [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], 
$\mathcal{O}(\varepsilon^{2})$ information is used to obtain the following formula (after notational changes) for the soliton centre, which is plotted in blue in the right panel of Figure 10:
\begin{equation}
x_{1}(T)=x_{1}(0)+\frac{2\Lambda_{1}(0)k_{1}(0)}{3q_{0}(0)}T+\frac{8\Lambda_{1}(0)^{3}k_{1}(0)}{9q_{0}(0)}T^{2}.
\end{equation} A possible explanation as to why the present theory is unable to capture the quadratic term in (7.43) is that this term could arise from a dependence of the velocity parameter 
$k_{1}$ on the slower time scale 
$\hat T=\varepsilon T=\varepsilon^2 t$. In particular, since
\begin{equation*}
\frac{8\Lambda_{1}(0)^{3}k_{1}(0)}{9q_{0}(0)}T^{2}=2\int_{0}^{T}\frac{8\Lambda_{1}(0)^{3}k_{1}(0)}{9q_{0}(0)}\, \hat{s}\,d\hat{s},
\end{equation*}an identical correction to the location of the soliton centre can be obtained by taking the evolution of the velocity parameter to be
\begin{equation}
k_{1}(\hat{T})=k_{1}(0)-\frac{8\Lambda_{1}(0)^{3}k_{1}(0)}{9q_{0}(0)}\hat{T}.
\end{equation} Recall that the parameter 
$k_{1}$ controls not only the velocity but also the trough amplitude of the dark soliton. As such, if the correction found in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] for the centre is to be interpreted instead as a correction to 
$k_{1}$ as in (7.44), then there should be an observable slow decrease in the trough amplitude over longer time scales. Numerical evidence suggests that this is the case, and that (7.44) could provide an appropriate correction to the trough amplitude. In the bottom right panel of Figure 11, numerical measurements of the trough amplitude are plotted in comparison with the hypothetical 
$\varepsilon^2 t$-dependent correction to 
$k_{1}$.
A space-time plot showing the development of a depressed shelf around the soliton under a dissipative perturbation. The blue dashed lines denote the boundaries of the shelf region, and the red dashed line corresponds to the path of the soliton core. In the top right panel, the background amplitude 
$q_{0}$ is shown to remain constant up to at least 
$t=400$. Each mark on the vertical axis represents an increment of 
$10^{-9}$. In bottom right panel, the slow decrease in the trough amplitude 
$|k_{1}|$ is shown up to 
$t=100$. The numerical measurements are compared to the hypothetical 
$\varepsilon^{2}t$-dependent formula given in (7.44), which was obtained by interpreting the quadratic correction to the soliton centre found in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] instead as a next-order correction to the velocity, as described in Remark 3. The initial parameters are the same as in Figure 9.
By similar logic, on account of (7.16) we also speculate that a discrepancy in the phase parameter 
$\sigma_{1}$ (e.g., in the case of linear and nonlinear damping) could arise due to 
$\hat{T}$-dependence of the background amplitude 
$q_{0}$. In the present case of dissipation, the numerics indicate that 
$q_{0}$ does in fact remain constant over long time intervals, and no dependence on 
$\hat{T}$ is detected, as shown in the top right panel of Figure 11. This is consistent with the fact that our prediction for 
$\sigma_{1}$ does remain accurate in this case (cf. Figure 10).
The height of the shelf and the phase gradient are given by
\begin{equation}
h^{\pm}=\mp\frac{\Lambda_{1}^{3}}{3q_{0}(k_{1}\pm q_{0})},\qquad\varphi_{\xi}^{\pm}=-\frac{2\Lambda_{1}^{2}}{3q_{0}(k_{1}\pm q_{0})}.
\end{equation} In this case, we see that since 
$|k_{1}| \lt q_{0}$, the shelf is depressed on both sides of the soliton, and the phase gradient is positive on the left and negative on the right. A space-time plot showing the development and propagation of the shelf is displayed in Figure 11, and a snapshot of the shelf region including the predicted height and phase gradient is shown in Figure 12.
The modulus (left) and phase (right) of a numerical simulation of a dark soliton under the influence of the dissipative perturbation 
$F[q]=iq_{xx}$ with 
$\varepsilon=0.02$ at 
$t=10$. The blue dashed lines denote the predicted boundaries of the shelf region. The red (resp., green) dashed lines represent the predictions for the shelf height and phase gradient on the right (resp., on the left) of the soliton. The initial parameters are the same as in Figure 9.
7.4. Self-steepening
We note that while the condition 
$\sigma_{1T}k_{1}=2x_{1T}\Lambda_{1}^{2}$ has appeared in all three previous examples, it is not universal. For example, in the case of a self-steepening perturbation
\begin{equation}
F_{0}=-i(|u|^{2}u)_{x}=-\Lambda_{1}^{2}\big(2ik_{1}\Lambda_{1}\mathop{\rm sech}\nolimits^{2}\xi\tanh\xi-k_{1}^{2}\mathop{\rm sech}\nolimits^{2}\xi-3\Lambda_{1}\mathop{\rm sech}\nolimits^{2}\xi\tanh^{2}\xi\big),
\end{equation}one finds after applying our method that 
$k_{1T}=0$, 
$\Lambda_{1T}=0$, 
$\sigma_{1T}=0$ and no prominent shelf is formed. The absence of the shelf is tied to the fact that for this perturbation, we have that
\begin{equation}
\langle\Upsilon_{0}^{\pm},W\rangle=0,
\end{equation} (cf. Eq. (4.27)). Thus, the singularity in the first-order correction integral (6.1) is fully removed, and in turn there is no secular growth. For the same reason, any perturbation satisfying (7.47) will not generate a shelf. In the self-steepening case, the only soliton parameter affected is the centre, whose evolution is obtained directly from the orthogonality condition 
$\langle\tilde\Upsilon_{1},W\rangle=0$ and is given by
\begin{equation}
x_{1T}=2k_{1}^{2}+\Lambda_{1}^2.
\end{equation} So, the centre is pushed forward linearly by the perturbation, which is similar to the situation for bright solitons under a self-steepening effect (see [Reference Yang54]). We note that our predictions here are significantly different than those given in [Reference Chen and Chen9] for the same perturbation. In that work, the amplitude and velocity of the dark soliton are predicted to evolve as a result of the perturbation. However, this is not supported by the numerics, which confirm that these parameters remain constant under the perturbation. Figure 13 demonstrates that our prediction for the soliton centre (7.48) remains accurate even up to 
$t=1/\varepsilon$.
Left: The modulus of a numerical simulation of a dark soliton under the influence of the self-steepening perturbation 
$F[q]=-i(|q|^{2}q)_{x}$ with 
$\varepsilon=0.02$ at 
$t=50$. The red dashed line denotes the location of the predicted soliton centre, while the blue dashed line marks where the centre would have been in the absence of the perturbation. Right: A comparison of our prediction for the evolution of the centre with the numerically measured values. The initial soliton parameters are 
$q_{0}(0)=2$, 
$k_{1}(0)=-1/2$, 
$x_{1}(0)=\sigma_{1}(0)=0$.
8. Conclusions
In this work, we revisited the integrable perturbation theory of the scalar defocusing NLS on a nontrivial background, addressing the conjecture posed in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] that, because of the existence of an expanding shelf that develops on each side of the soliton due to the perturbation, the squared eigenfunctions associated with the soliton are an insufficient basis, and questioning the existence of a closure relation for this problem, and the ability of integrable perturbation theory to accurately describe the radiation shelf.
We provided a proof of the completeness of the squared eigenfunctions, a crucial difference with respect to the earlier works lying in properly accounting for the singularities of the scattering data at the points 
$\pm q_0$, which are branch points for the continuous spectrum. Moreover, we showed that the first-order correction can in fact be used to describe the soliton shelf, which results from the singularities of the scattering data at the points 
$\pm q_0$, and we obtained estimates for the magnitude, velocity and phase gradient on each side of the shelf. Using the 1-soliton closure relation, suitable suppression of secular growth, and the asymptotic behaviour of the first-order correction we determined the adiabatic evolution of the soliton parameters as functions of a slow time variable 
$T=\varepsilon t$.
All the results were corroborated by direct numerical simulations, and compared with the results of the direct perturbation theory in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], and with the earlier works using perturbation theory based on the squared eigenfunctions. In particular, we showed that the estimates we obtained for all soliton parameters agree with the ones in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] to order 
$\varepsilon$, the only difference being for the soliton centre, which, unlike the other parameters, obtained from perturbed conserved quantities up to 
$\mathcal{O}(\varepsilon)$, is determined in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1] in terms of differential equations obtained from the Hamiltonian at 
$\mathcal{O}(\varepsilon^2)$. We also explained the shortcomings of the predictions for the adiabatic evolution of the soliton parameters in the earlier works within the framework of the integrable perturbation theory.
For some perturbations, our predictions for the soliton centre and phase do not match the numerical data on the time scale of 
$T=\varepsilon t=\mathcal{O}(1)$. A similar problem for the phase can also be seen in some of the examples considered in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1]. The reason for this discrepancy has been discussed in Remark 3, and it lies in the fact that since the slow-time-dependent velocity must be integrated in time when passing to the co-moving frame (cf Eq. (5.5)), it is actually the case the dependence of 
$k_1$ on an even slower time scale 
$\varepsilon^2t$ could contribute nontrivially to the soliton centre. Thus, to ensure a correct prediction for the soliton centre, one would need to include another time scale in the perturbation theory. A similar argument would apply to the soliton phase 
$\sigma_1$, on account of the integration in time of the slow-varying background amplitude 
$q_0$ in Eq. (7.16). This indicates that in order to achieve correct estimates up to 
$t=\mathcal{O}(1/\varepsilon)$ for the soliton centre and phase, one would need to compute the second-order correction terms. It is important to stress that this is an intrinsic feature of the problem, and not of the perturbative approach, as both our estimates, based on the integrable perturbation theory, and the ones in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], based on direct perturbation theory, with multiple scale expansions and perturbed conservation laws, exhibit the same behaviour. In fact, one could argue that similar considerations should be applied even in the bright soliton perturbation theory, where the same transformation to the co-moving frame is routinely applied.
We believe the present work may open up a new vein of research on integrable perturbation theory for integrable systems on a nonzero background, and the methodology can also be generalised to discrete and coupled integrable systems. Another natural follow-up is to explore ways to improve the accuracy of the estimates for the adiabatic evolution of the soliton parameters for times up to 
$\mathcal{O}(1/\varepsilon)$. This could be done by introducing an additional time scale 
$\varepsilon^2t$ and a second-order correction to the potential, which will require handling quadratic contributions in the first-order correction appearing in the orthogonality conditions. Alternatively, one could generalise the approach used to determine the slow-time evolution of the soliton centre in [Reference Ablowitz, Nixon, Horikis and Frantzeskakis1], which makes use of the Hamiltonian at 
$\mathcal{O}(\varepsilon^2)$, and consider also the other perturbed conserved quantities to 
$\mathcal{O}(\varepsilon^2)$ in order to get additional corrections to the soliton phase, amplitude and velocity.
Acknowledgements
BP and NJO gratefully acknowledge partial support for this work from the NSF, under grant DMS-2406626. NJO would like to thank Dr. Sathyanarayanan Chandramouli for assistance with the numerics. BP acknowledges the Fulbright Foundation in Greece and the Fulbright program, and the Mathematics Department of the University of Ioannina, Greece, for the kind hospitality during the revision of this work. The authors would also like to thank the anonymous reviewers for their comments, which helped improve the paper.
Appendix A. General completeness relation
To justify the general completeness relation (3.26), note that as 
$|z|\rightarrow\infty$, we have that 
$a(z)\sim 1$, 
$\bar{a}(z)\sim 1$ and
\begin{equation}
\eta(x,z)\chi(y,z)^T\sim e^{2i\lambda(z)(x-y)}\,\text{diag}\,(0,1),\quad\bar\eta(x,z)\bar\chi(y,z)^T\sim e^{-2i\lambda(z)(x-y)}\,\text{diag}\,(1,0),
\end{equation}and we recall that 
$\lambda(z)=(z-q_0^2/z)/2\sim z/2$ and 
$z\to \infty$. These are obtained directly from the large-
$z$ asymptotic behaviour of the Jost eigenfunctions, which can be found, for example, in [Reference Gkogkou, Prinari and Trogdon21, Reference Prinari46] where the same notations as ours are used. From the analyticity properties discussed in Sections 2 and 3 and the above asymptotics, we then get the identities
\begin{equation}
\int_{\Gamma}\frac{1}{\gamma(\zeta)^{2}a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta=\int_{\Gamma}e^{i\zeta(x-y)}d\zeta\,\text{diag}\,(0,1),
\end{equation}
\begin{equation}
\int_{\bar\Gamma}\frac{1}{\gamma(\zeta)^{2}\bar a(\zeta)^{2}}\bar\eta(x,\zeta)\bar\chi(y,\zeta)^{T}d\zeta=\int_{\bar\Gamma}e^{-i\zeta(x-y)}d\zeta\,\text{diag}\,(1,0),
\end{equation}
$\Gamma$ is a contour from 
$-\infty+i0$ to 
$+\infty+i0$ that passes above all zeros 
$\zeta_{j}$ of 
$a(z)$, and 
$\bar\Gamma$ is a contour from 
$-\infty-i0$ to 
$+\infty-i0$ that passes below all zeros 
$\zeta_{j}^{*}$ of 
$\bar{a}(z)=a^{*}(z^{*})$. For reference, Figure 14 displays the contours 
$\Gamma$ and 
$\bar\Gamma$, as well as the additional contours mentioned later in this section. Since the integrands on the right-hand sides in (A.2) are entire functions that decay exponentially as the radius of the contour of integration increases in their respective half-planes, both contours can be deformed to the real axis, and using 
$\int_{-\infty}^{\infty}e^{\pm i\zeta(x-y)}d\zeta=2\pi\delta(x-y)$ we get
\begin{equation}
\int_{\Gamma}\frac{1}{\gamma(\zeta)^{2}a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta+
\int_{\bar\Gamma}\frac{1}{\gamma(\zeta)^{2}\bar a(\zeta)^{2}}\bar\eta(x,\zeta)\bar\chi(y,\zeta)^{T}d\zeta=2\pi\delta(x-y)I.
\end{equation}In the upper-half plane, the contour 
$\Gamma$ (blue) passes above all zeros of 
$a(z)$ and the contour 
$C$ (orange) passes below the zeros of 
$a(z)$ while being indented around the branch points 
$z=\pm q_{0}$. Their counterparts 
$\bar\Gamma$ (red) and 
$\bar{C}$ (green) are shown in the lower-half plane. The contour 
$\bar\Gamma'$ (magenta) represents the image of 
$\bar\Gamma$ after the change of variable 
$z'=q_{0}^{2}/z$.
Using (3.20) and the similar symmetry satisfied by the adjoint squared eigenfunctions
\begin{equation}
\bar\chi(x,z)=\frac{q_{0}^{2}}{z^{2}}e^{-2i\theta}\chi(x,q_{0}^{2}/z),
\end{equation}as well as 
$\bar{a}(z)=e^{-2i\theta}a(q_{0}^{2}/z)$, we rewrite the second integral to get
\begin{eqnarray}
2\pi\delta(x-y)I&=&\int_{\Gamma}\frac{1}{\gamma(\zeta)^{2}a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta\nonumber\\ &&-\,
\int_{\bar\Gamma}\frac{1}{\gamma(q_{0}^{2}/\zeta)^{2} a(q_{0}^{2}/\zeta)^{2}}\eta(x,q_{0}^{2}/\zeta)\chi(y,q_{0}^{2}/\zeta)^{T}d\zeta.\end{eqnarray} Now, make the change of variable 
$\zeta'=q_{0}^{2}/\zeta$ in the second integral. The contour 
$\bar\Gamma$ is mapped to a contour 
$\bar\Gamma'$ in the upper-half plane that begins at 
$-0+i0$, travels clockwise while remaining within the circle of radius 
$q_{0}$, and finally approaches 
$+0+i0$ (shown in Figure 14). Since all potential poles are on the circle 
$|z|=q_{0}$, integrating along 
$\bar\Gamma'$ is equivalent to integrating along 
$\Gamma$ with the opposite orientation. With this, the two integrals in (A.5) can be combined into
\begin{equation}
\int_{\Gamma}\frac{1}{\gamma(\zeta)^{2}a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}\left(1-\frac{q_{0}^{2}}{\zeta^{2}}\right)d\zeta=\int_{\Gamma}\frac{1}{\gamma(\zeta)a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta.
\end{equation} To explicitly include the contributions from the zeros of 
$a(z)$, which are double poles of the integrand above, note that
\begin{equation}
\left(\int_{C}-\int_{\Gamma}\right)\frac{1}{\gamma(\zeta)a(\zeta)^{2}}\eta(x,\zeta)\chi(y,\zeta)^{T}d\zeta=2\pi i\sum_{j=1}^{J}\mathop{\rm Res}\limits_{z=\zeta_j}\frac{\eta(x,z)\chi(y,z)^{T}}{\gamma(z)a(z)^{2}},
\end{equation}where 
$C$ is a contour from 
$-\infty+i0$ to 
$+\infty+i0$ that follows the real axis while being indented in the upper half plane to avoid 
$z=\pm q_{0}$. The residue at each 
$\zeta_{j}$ is found to be:
\begin{gather}\frac{1}{\gamma(\zeta_{j})a'(\zeta_{j})^{2}}\Bigg[\frac{\gamma(\zeta_{j})a''(\zeta_{j})-\gamma'(\zeta_{j}){a}'(\zeta_{j})}{\gamma(\zeta_{j}){a}'(\zeta_{j})}\eta(x,\zeta_{j})\chi(y,\zeta_{j})^{T} +\,\eta'(x,\zeta_{j})\chi(y,\zeta_{j})^{T}+\eta(x,\zeta_{j})\chi'(y,\zeta_{j})^{T}\Bigg]. \end{gather}Substituting the above expression into (3.26) leads to (3.27).
Appendix B. Derivation of 1-soliton Jost eigenfunctions
In the case of no reflection and a single soliton (
$\rho\equiv0$ and 
$J=1$), (2.19) reduces to the linear algebraic system
\begin{equation}
\psi(x,t,z)e^{-i\Omega(x,t,z)}=\begin{bmatrix}
-iq_{+}/z\\1
\end{bmatrix}+\frac{e^{-i\Omega(x,t,\zeta_{1}^{*})}\bar\psi(x,t,\zeta_{1}^{*})C_{1}^{*}}{z-\zeta_{1}^{*}},
\end{equation}
\begin{align}
\bar\psi(x,t,z)e^{i\Omega(x,t,z)}=\begin{bmatrix}
1 \nonumber\\ iq_{+}^{*}/z
\end{bmatrix}+\frac{e^{i\Omega(x,t,\zeta_{1})}\psi(x,t,\zeta_{1})C_{1}}{z-\zeta_{1}}.
\end{align}
$\zeta_{1}=k_{1}+i\Lambda_{1}$ with 
$k_{1}\in\mathbb{R}$, 
$\Lambda_{1} \gt 0$. We need to review here some properties of the scattering coefficient 
$a(z)$. Taking into account its analyticity properties, its zeros, and the symmetries, one can obtain the following representation (trace formula) for 
$a(z)$ for 
$z\in \mathbb{C}^+$:
\begin{equation}
a(z)=\prod_{j=1}^J\frac{z-\zeta_j}{z-\zeta_j^*}\exp{\left[ -\frac{1}{2\pi i} \int_{-\infty}^{\infty} \frac{\log(1-\rho(z)\rho^*(\zeta^*))}{\zeta -z}d\zeta \right]}\,,
\end{equation}where 
$J$ is the number of (simple) discrete eigenvalues. Recalling that 
$a(z)\rightarrow q_+/q_-$ as 
$z\rightarrow 0$, we conclude that the potential satisfies
\begin{equation}
\frac{q_+}{q_-}=\prod_{j=1}^J\frac{\zeta_j}{\zeta_j^*} \, \exp{\left[-\frac{1}{2\pi i} \int_{-\infty}^{\infty} \frac{\log(1-\rho(z)\rho^*(\zeta^*))}{z}d\zeta \right]}\,,
\end{equation}which is often referred to as the 
$\theta$-condition [Reference Faddeev and Takhtajan16]. Due to our choice of boundary conditions 
$q_{\pm}=q_{0}e^{\pm i\theta}$, the 
$\theta$-condition with 
$J=1$ and 
$\rho(z)\equiv 0$ demands that 
$\zeta_{1}=q_{+}$. Furthermore, note that 
$\Omega_{1}:=\Omega(x,t,\zeta_{1})=-\Omega(x,t,\zeta_{1}^{*})$. With this in mind, evaluating (B.1a) at 
$z=\zeta_{1}$ and (B.1b) at 
$z=\zeta_{1}^{*}$ gives
\begin{equation}
\psi(x,t,\zeta_{1})=\begin{bmatrix}
-i\\1
\end{bmatrix}e^{i\Omega_{1}}+\frac{e^{2i\Omega_{1}}\bar\psi(x,t,\zeta_{1}^{*})C_{1}^{*}}{\zeta_{1}-\zeta_{1}^{*}},
\end{equation}
\begin{equation}
\bar\psi(x,t,\zeta_{1}^{*})=\begin{bmatrix}
1\\i
\end{bmatrix}e^{i\Omega_{1}}+\frac{e^{2i\Omega_{1}}\psi(x,t,\zeta_{1})C_{1}}{\zeta_{1}^{*}-\zeta_{1}}.
\end{equation}
$\zeta_{1}-\zeta_{1}^{*}=2i\Lambda_{1}$ and taking 
$C_{1}\in\mathbb{R}$ according to (2.18) with 
$\alpha_{1}=\theta$, solving for 
$\bar\psi(x,t,\zeta_{1}^*)$ yields
\begin{equation}
\bar\psi(x,t,\zeta_{1}^*)=e^{i\Omega_{1}}\begin{bmatrix}
1\\i
\end{bmatrix}\frac{2\Lambda_{1}}{2\Lambda_{1}-C_{1}e^{2i\Omega_{1}}}.
\end{equation} Defining a new parameter (the soliton centre) 
$x_{1}$ through 
$C_{1}=-2\Lambda_{1}e^{2\Lambda_{1}x_{1}}$, this becomes
\begin{equation}
\bar\psi(x,t,\zeta_{1}^*)=e^{i\Omega_{1}}\begin{bmatrix}
1\\i
\end{bmatrix}\frac{1}{1+e^{2i\Omega_{1}+2\Lambda_{1}x_{1}}}.
\end{equation}This can be written in terms of the travelling coordinate defined in (4.8),
\begin{equation}
\xi=\Lambda_{1}(x+2k_{1}t-x_{1})=-i\Omega_{1}-\Lambda_{1}x_{1},
\end{equation}in which case we have
\begin{equation}
\bar\psi(x,t,\zeta_{1}^*)=e^{i\Omega_{1}}\begin{bmatrix}
1\\i
\end{bmatrix}\frac{1}{1+e^{-2\xi}}=e^{i\Omega_{1}}\begin{bmatrix}
1\\i
\end{bmatrix}\frac{1}{2}e^{\xi}\mathop{\rm sech}\nolimits\xi.
\end{equation}Substituting (B.8) into (B.1a) and again using (B.7), we obtain
\begin{equation}
\psi(x,t,z)e^{-i\Omega(x,t,z)}=\begin{bmatrix}
-iq_{+}/z\\1
\end{bmatrix}-\frac{\Lambda_{1}e^{-2\xi}}{z-\zeta_{1}^{*}}\begin{bmatrix}
1\\i
\end{bmatrix}\mathop{\rm e}\nolimits^{\xi}\mathop{\rm sech}\nolimits\xi.
\end{equation}By components, this is:
\begin{equation}
\psi_{1}(x,t,z)=\Big(-iq_{+}z^{-1}-\frac{1}{z-\zeta_{1}^{*}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{i\Omega(x,t,z)},
\end{equation}
\begin{equation}
\psi_{2}(x,t,z)=\Big(1-i\frac{1}{z-\zeta_{1}^{*}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{i\Omega(x,t,z)}.
\end{equation}
$\bar\psi(x,t,z)$ can then be recovered through the symmetry (2.9),
\begin{equation}
\bar{\psi}_{1}(x,t,z)=\psi_{2}^{*}(x,t,z^{*})=\Big(1+i\frac{1}{z-\zeta_{1}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{-i\Omega(x,t,z)},
\end{equation}
\begin{equation}
\bar\psi_{2}(x,t,z)=\psi_{1}^{*}(x,t,z^{*})=\Big(iq_{-}z^{-1}-\frac{1}{z-\zeta_{1}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{-i\Omega(x,t,z)}.
\end{equation}
$\phi(x,t,z)$ can be found through the relation (4.12),
\begin{equation}
\phi_{1}(x,t,z)=\Big(\frac{z-\zeta_{1}}{z-\zeta_{1}^{*}}+i\frac{1}{z-\zeta_{1}^{*}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{-i\Omega(x,t,z)},
\end{equation}
\begin{equation}
\phi_{2}(x,t,z)=\Big(iq_{-}z^{-1}\frac{z-\zeta_{1}}{z-\zeta_{1}^{*}}-\frac{1}{z-\zeta_{1}^{*}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{-i\Omega(x,t,z)},
\end{equation}
$\bar\phi(x,t,z)$ can be recovered through symmetries,
\begin{equation}
\bar\phi_{1}(x,t,z)=\phi_{2}^{*}(x,t,z^{*})=\Big(-iq_{+}z^{-1}\frac{z-\zeta_{1}^{*}}{z-\zeta_{1}}-\frac{1}{z-\zeta_{1}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{i\Omega(x,t,z)},
\end{equation}
\begin{equation}
\bar\phi_{2}(x,t,z)=\phi_{1}^{*}(x,t,z^{*})=\Big(\frac{z-\zeta_{1}^{*}}{z-\zeta_{1}}-i\frac{1}{z-\zeta_{1}}\Lambda_{1}e^{-\xi}\mathop{\rm sech}\nolimits\xi\Big)e^{i\Omega(x,t,z)}.
\end{equation}Appendix C. Evaluation of the integrals in Eq. (6.4)
In this Appendix, we provide the detailed calculations of the integrals 
$I_\pm$ in Eq. (6.4). First, consider the integral 
$I_{-}$ given by
\begin{equation}
I_{-}={\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{\lambda(\zeta)^2}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta,
\end{equation}which has a pole at 
$\zeta=-q_{0}$. First, using 
$\lambda(\zeta)=\frac{1}{2}(1-q_{0}\zeta^{-1})(\zeta+q_{0})$, this can be written as
\begin{equation}
I_{-}=4\,{\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{(\zeta+q_{0})^{2}}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{(1-q_{0}\zeta^{-1})^{2}}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta.
\end{equation}Applying integration by parts,
\begin{equation}
I_{-}=4\,{\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{\zeta+q_{0}}\partial_{\zeta}\left\{\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{(1-q_{0}\zeta^{-1})^{2}}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}\right\}d\zeta.
\end{equation}Expanding the derivative, we can split (C.3) into
\begin{equation}
I_{-}=J_{1}+i\xi J_{2}+it J_{3},
\end{equation}where the three integrals are given by
\begin{equation}
J_{1}={\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{\zeta+q_{0}}\frac{-8q_{0}\zeta^{-2}}{(1-q_{0}\zeta^{-1})^{3}}\big[1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}\big]e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta,
\end{equation}
\begin{equation}
J_{2}={\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{\zeta+q_{0}}\frac{1}{(1-q_{0}\zeta^{-1})^{2}}\frac{8\lambda'(\zeta)}{\Lambda_{1}}\big[1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}\big]e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta,
\end{equation}
\begin{equation}
J_{3}={\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{\zeta+q_{0}}(-16)\frac{\lambda'(\zeta)(k(\zeta)-k_{1})+\lambda(\zeta){k}'(\zeta)}{(1-q_{0}\zeta^{-1})^{2}}e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta.
\end{equation}
$J_{1}$, using the fact that the dominant contribution to the integral comes from the pole 
$\zeta\approx-q_{0}$, in which case 
$\lambda(\zeta)\approx\zeta+q_{0}$ and 
$k(\zeta)\approx-q_{0}$, we get
\begin{equation}
J_{1}\sim\frac{1}{q_{0}}\,{\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{\zeta+q_{0}}\big[1-e^{-4i(\zeta+q_{0})(k_{1}+q_{0})t}\big]e^{\frac{2i}{\Lambda_{1}}(\zeta+q_{0})\xi}d\zeta
\end{equation}To evaluate principal value integrals of this form, we can use the formula
\begin{equation}
{\int{\kern-11.2pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.2pt}\infty}\frac{e^{ip\vartheta}}{p}dp= i\pi\mathop{\rm sgn}\nolimits \vartheta.
\end{equation} To this end, 
$J_{1}$ can be split again into
\begin{equation}
J_{1}\sim\frac{1}{q_{0}}\,{\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{\zeta+q_{0}}e^{\frac{2i}{\Lambda_{1}}(\zeta+q_{0})\xi}d\zeta-\frac{1}{q_{0}}\,{\int{\kern-11.5pt}-}_{{\kern-4.5pt}-\infty}^{{\kern2.5pt}0}\frac{1}{\zeta+q_{0}}e^{\frac{2i}{\Lambda_{1}}(\zeta+q_{0})\left[\xi-2\Lambda_{1}(k_{1}+q_{0})t\right]}d\zeta.
\end{equation}Applying the formula (C.7) gives
\begin{equation}
J_{1}\sim\frac{\pi i}{q_{0}}\left\{1-\mathop{\rm sgn}\nolimits\big[\xi-2\Lambda_{1}(k_{1}+q_{0})t\big]\right\}.
\end{equation} If 
$\xi\gg2\Lambda_{1}(k_{1}+q_{0})t$, then 
$J_{1}\sim0$. In the right shelf region 
$1\ll\xi\ll2\Lambda_{1}(k_{1}+q_{0})t$, 
$J_{1}\sim2\pi i/q_{0}$. Following similar steps for 
$J_{2}$, we obtain
\begin{equation}
J_{2}\sim\frac{2\pi i}{\Lambda_{1}}\left\{1-\mathop{\rm sgn}\nolimits\big[\xi-2\Lambda_{1}(k_{1}+q_{0})t\big]\right\}.
\end{equation}which gives in the shelf region 
$J_{2}=4\pi i/\Lambda_{1}$. Lastly, noting that 
$\lambda'(\zeta)\approx1$ and 
$k'(\zeta)\approx0$ near the pole, we get
\begin{equation}
J_{3}\sim4\pi i(k_{1}+q_{0})\,\mathop{\rm sgn}\nolimits\big[\xi-2\Lambda_{1}(k_{1}+q_{0})t\big],
\end{equation}which shows that
\begin{equation}
J_{3}\sim\begin{cases}
4\pi i(k_{1}+q_{0}),\quad&\xi\gg2\Lambda_{1}(k_{1}+q_{0})t\\
-4\pi i(k_{1}+q_{0}),\quad&1\ll\xi\ll2\Lambda_{1}(k_{1}+q_{0})t
\end{cases}.
\end{equation} We now move on to evaluate the integral 
$I_{+}$, given by
\begin{equation}
I_{+}=4\,{\int{\kern-11.2pt}-}_{{\kern-3.5pt}0}^{{\kern2.2pt}\infty}\frac{1}{(\zeta-q_{0})^{2}}\frac{1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}}{(1+q_{0}\zeta^{-1})^{2}}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta,
\end{equation}which has a pole at 
$\zeta=q_{0}$. Similarly to the previous calculations, applying integration by parts splits the integral into
\begin{equation}
I_{+}=H_{1}+i\xi H_{2}+itH_{3},
\end{equation}where
\begin{equation}
H_{1}={\int{\kern-11.2pt}-}_{{\kern-3.5pt}0}^{{\kern2.5pt}\infty}\frac{1}{\zeta-q_{0}}\frac{8q_{0}\zeta^{-2}}{(1+q_{0}\zeta^{-1})^{3}}\big[1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}\big]e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta
\end{equation}
\begin{equation}
H_{2}={\int{\kern-11.2pt}-}_{{\kern-3.5pt}0}^{{\kern2.2pt}\infty}\frac{1}{\zeta-q_{0}}\frac{1}{(1+q_{0}\zeta^{-1})^{2}}\frac{8\lambda'(\zeta)}{\Lambda_{1}}\big[1-e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}\big]e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta
\end{equation}
\begin{equation}
H_{3}={\int{\kern-11.2pt}-}_{{\kern-3.5pt}0}^{{\kern2.2pt}\infty}\frac{1}{\zeta-q_{0}}(-16)\frac{\lambda'(\zeta)(k(\zeta)-k_{1})+\lambda(\zeta){k}'(\zeta)}{(1+q_{0}\zeta^{-1})^{2}}e^{4i\lambda(\zeta)(k(\zeta)-k_{1})t}e^{2i\frac{\lambda(\zeta)}{\Lambda_{1}}\xi}d\zeta.
\end{equation}
\begin{equation}
H_{1}\sim\frac{\pi i}{q_{0}}\{1-\mathop{\rm sgn}\nolimits\big[\xi-2\Lambda_{1}(k_{1}-q_{0})t\big]\},
\end{equation}
\begin{equation}
H_{2}\sim\frac{2\pi i}{\Lambda_{1}} \{ 1-\mathop{\rm sgn}\nolimits\big[\xi-2\Lambda_{1}(k_{1}-q_{0})t\big] \}. clscurl;,
\end{equation} Since we are studying the large positive 
$\xi$ limit, and 
$k_{1} \lt q_{0}$, we have that 
$H_{1}\sim0$ and 
$H_{2}\sim0$. Thus, the contributions to the height and phase gradient of the right side of the shelf only come from the pole 
$\zeta=-q_{0}$. The third integral is found to be
\begin{equation}
H_{3}\sim4\pi i(k_{1}-q_{0})\mathop{\rm sgn}\nolimits\big[\xi-2\Lambda_{1}(k_{1}-q_{0})t\big],
\end{equation}which is limit gives 
$H_{3}\sim4\pi i(k_{1}-q_{0})$. To summarise, putting everything into (6.4), we have that if 
$\xi\gg2\Lambda_{1}(k_{1}+q_{0})t$,
\begin{equation}
\tilde{q}^{+}\sim-\mathcal{F}(-q_{0})4\pi (k_{1}+q_{0})t-\mathcal{F}(q_{0})4\pi (k_{1}-q_{0})t=0,
\end{equation}as expected, and in the right shelf region 
$1\ll\xi\ll2\Lambda_{1}(k_{1}+q_{0})t$,
\begin{equation}
q^{+}\sim\mathcal{F}(-q_{0})\left[\frac{2\pi i}{q_{0}}-\frac{4\pi}{\Lambda_{1}}\xi+4\pi (k_{1}+q_{0})t\right]-\mathcal{F}(q_{0})4\pi (k_{1}-q_{0})t.
\end{equation}Putting in (6.5), this reduces to the expression given in Section 6. The calculation of the left side of the shelf from the large negative 
$\xi$ limit proceeds in an analogous way. In that case, the contributions to the height and phase gradient of the shelf come only from the pole 
$\zeta=q_{0}$.
Open access
