1 Introduction

The scalar nonlinear Schrödinger (NLS) equation is one of the most celebrated and well-studied nonlinear partial differential equations (PDEs) in physics and mathematics. One of the reasons for its prominence is that almost any dispersive, energy-preserving system reduces to it, in some appropriate limit. In fact, the NLS equation provides a canonical description for the envelope dynamics of a quasi-monochromatic plane wave propagating in a weakly nonlinear dispersive medium whenever dissipation can be neglected.

Historically, the origin of the NLS equation can be traced back to the 1950’s, in the work of Ginzburg and Landau [94] and Ginzburg [93] on the macroscopic theory of superconductivity, and also of Ginzburg and Pitaevskii [95] on the theory of superfluidity. Nonetheless, it was only in the 1960’s that Chiao et al [56] and Talanov [166, 167] revealed the wider physical importance of the NLS equation in connection with the phenomenon of self-focusing of an electromagnetic beam in nonlinear media. Shortly thereafter, Zakharov derived the NLS equation for small-amplitude waves in water of infinite depth [187], and Benney and Roskes generalized the derivation to the case of finite depth [28]. Essentially, in the water-wave context the (multidimensional) NLS equation is obtained from the Euler–Bernoulli equations for the dynamics of an ideal (i.e., incompressible, irrotational and inviscid) fluid under the assumption of a small-amplitude, quasi-monochromatic wave expansion. The NLS equation has since been proposed as a model in such diverse fields as plasma physics [188], nonlinear optics [100, 101], magnetic spin waves [106, 196], low-temperature physics and Bose-Einstein condensation (BEC) [145, 146], etc, thus fully justifying its relevance.

In the \(1+1\) dimensional case (one temporal and one spatial dimension, or two spatial dimensions with one being the direction of propagation, and one the transverse dimension), the scalar NLS equation can be written in terms of dimensionless variables as

$$\begin{aligned} iq_t + q_{xx} - 2\sigma |q|^2 q =0, \qquad \sigma =\pm 1 \end{aligned}$$
(1.1)

(\(q=q(x,t)\) is a complex function, and subscripts x and t denote partial differentiation throughout) where \(\sigma\) distinguishes between two inequivalent versions, characterizing, respectively, the focusing or anomalous (\(\sigma =-1\)) and the defocusing or normal (\(\sigma =1\)) dispersion regimes. The focusing NLS equation admits the usual, bell-shaped soliton:

$$\begin{aligned} q(x,t)=2i\eta e^{-2i(\xi x+2(\xi ^2-\eta ^2)t)+i \theta }\mathop {\textrm{sech}}\nolimits [2\eta (x+4\xi t-x_o)], \end{aligned}$$
(1.2)

where \(\xi ,\eta ,x_o,\theta\) are arbitrary real parameters, which are usually referred to as “bright” solitons, while the defocusing NLS only admits soliton solutions with nontrivial boundary conditions, (cf. Eq. (2.1)). These solitons with nonzero boundary conditions (NZBC) are the so-called “dark” solitons, and appear as localized dips of intensity on a nonzero background field (see Sect. 2 for details). One of the distinguishing properties of both bright and dark (scalar) solitons is that they interact trivially with each other, in the sense that they preserve asymptotically their shape, amplitude and velocity, and the only effect of the interaction is an overall shift in their position and phase.

Mathematically, the NLS Eq. (1.1) attains broad significance since, in addition to single and multisoliton solutions, it possesses an infinite number of conserved quantities, and, most notably, it is integrable via the Inverse Scattering Transform (IST). The IST is a method that allows one to linearize a class on nonlinear evolution equations (which are thereby referred to as “integrable systems”), and to reduce the solution of the initial-value problem for the nonlinear PDE to a series of linear steps, much like the Fourier transform for linear PDEs (and in fact in an engineering context the IST is often referred to as nonlinear Fourier Transform, or NFT). Notably, the IST for the focusing NLS equation was developed in the pioneer works by Zakharov & Shabat and Ablowitz, Kaup, Newell & Segur under the assumption that the potential q(xt) is rapidly decaying as \(|x|\rightarrow \pm \infty\) [7, 192].

An essential pre-requisite of the IST method is the association of the nonlinear integrable PDE with a pair of linear problems—a “Lax pair”—such that the given nonlinear PDE results as a compatibility condition between them [126]. Specifically, the matrix pair \({\textbf{U}},{\textbf{V}}\) is a Lax pair for the nonlinear PDE

$$\begin{aligned} q_t=F[x,t,q,q_x,q_{xx},\dots ] \end{aligned}$$

if the compatibility condition of the overdetermined linear system of differential equations

$$\begin{aligned} \varphi _x={\textbf{U}} \varphi , \qquad \varphi _t={\textbf{V}} \varphi \end{aligned}$$

[i.e., the equality of the mixed derivatives \(\varphi _{xt}=\varphi _{tx}\), or, equivalently, the so-called zero curvature condition \({\textbf{U}}_t-{\textbf{V}}_x={\textbf{U}}{\textbf{V}} - {\textbf{V}}{\textbf{U}}\)] is identically satisfied whenever q(xt) is a solution of the nonlinear PDE. Here \({\textbf{U}}, {\textbf{V}}\) are in general matrix functions of q(xt) and its spatial derivatives, and of a complex parameter k which is assumed to be time-independent. The solution of the initial-value problem by IST proceeds in three steps, as follows:

  1. 1.

    the direct problem—the transformation of the initial datum q(x, 0) into the transformed “scattering data” S(k, 0);

  2. 2.

    time dependence—the evolution of the scattering data often according to simple, explicitly solvable evolution equations (i.e., finding S(kt));

  3. 3.

    the inverse problem—the recovery of the evolved solution q(xt) in terms of the evolved scattering data S(kt).

Both the direct and the inverse problems make use of the first operator in the Lax pair, the “scattering problem”, while the time evolution is determined by the second operator in the Lax pair. In the direct problem, the first step amounts to constructing eigenfunction solutions of the associated linear problem. These eigenfunctions depend on both the original spatial variables and on the spectral parameter k, and are usually defined in terms of integral equations with a prescribed plane-wave behavior as \(x\rightarrow \pm \infty\). Second, with these eigenfunctions, one determines scattering data that are independent of the original spatial variables. The scattering data typically consist of: (i) a reflection coefficient \(\rho (k,t)\), with \(\rho (k,0)\) being the analog of the Fourier transform of the initial datum in the linear case; (ii) discrete eigenvalues \(k_1,\dots ,k_n\), which are values of the spectral parameter for which the scattering problem admits \(L^2\) eigenfunctions; further, each \(k_j\in {\mathbb {C}}\) contributes a soliton to the solution, and real and imaginary parts of the discrete eigenvalue may specify the amplitude and velocity of the soliton, cf. (1.2) with \(k_1=\xi +i\eta\); (iii) the normalization constants \(C_1,\dots ,C_n\) associated to each discrete eigenvalue, which specify the position and phase of each soliton (cf. (1.2) with \(x_o\) and \(\theta\) related, respectively, to the modulus and argument of \(C_1\)). The analytiticy properties of the eigenfunctions as functions of the spectral parameter k as established in the direct problem are key to the formulation of the inverse problem. In turn, in the inverse problem, one first recovers the eigenfunctions from the evolved scattering data, either in terms of Marchenko integral equations or, in recent years, by formulating a suitable Riemann–Hilbert problem for the eigenfunctions. Finally, one reconstruct the potential q(xt) in the original variables from the evolved eigenfunctions, usually via an asymptotic expansion in k of the eigenfunctions themselves.

It is worth mentioning that with the IST machinery one can do more than just solve the initial-value problem; one can also construct special solutions by positing an elementary ansatz for the scattering data S(k, 0) and then using the equations of the inverse problem to obtain the corresponding solution q(xt). For instance, multisoliton (reflectionless) solutions can be constructed in this way. Although direct methods such as Darboux transformations or Hirota’s bilinear forms might be more efficient if one is only interested in obtaining explicit solutions, the benefit of the IST construction is that as a by-product one gains insight into the spectral characterization of the solutions, which is not achievable by other methods. Also, the IST provides an effective way to study the asymptotic (long-time) behavior of the solutions, e.g., via the nonlinear steepest descent technique [61] (see Sect. 5), as well as a way to study their stability using the so-called squared eigenfunctions [91].

Vector generalizations of the scalar NLS equations arise as relevant physical models under conditions similar to those described by NLS whenever there are suitable multiple wavetrains moving with nearly the same group velocity. Of particular relevance is the so-called Manakov system [133], which accounts for the fact that in optical fibers and waveguides the propagating electric field has two polarized components transverse to the direction of propagation, and therefore the appropriate model is a two-component coupled NLS system. Quite naturally, vector NLS systems (VNLS) with any number of components have been theoretically and experimentally investigated. Certain matrix NLS systems have also been proposed as models for one-dimensional spin-1 BECs [102, 103]. In some cases, these vector and matrix NLS systems are still completely integrable, and their initial-value problem can be effectively solved by the IST.

It is important to point out that while the IST for NLS systems with rapidly decaying potentials has been available as a tool for 50 years [133, 192], and a generalization to certain non-zero boundary conditions for the scalar defocusing NLS was developed by Zakharov and Shabat already in 1973 [193] (see also [74] for a detailed and rigorous description), other problems have remained open until very recently. For instance, the IST for the defocusing vector NLS with NZBC was only fully developed in 2006 [149] (some preliminary results can be found in [88]; see also [36] for additional details), and the IST for square matrix NLS systems with NZBC was carried out shortly afterwards [104, 175] (see also [151]). Similarly, the IST for the focusing NLS with NZBC was obtained by Biondini and Kovačič in 2014 [35] (some partial results can be found in [131]), and the solution of the problem for the focusing Manakov system followed in 2015 [121]. Importantly, unlike what happens in the case of rapidly decaying potentials or in the square matrix case, where the generalization of the IST from two to an arbitrary number of components is straightforward, the IST for the defocusing Manakov system with NZBC to more than two components is a highly nontrivial problem which was addressed in [37, 150]. We should also mention that most of the above mentioned works deal with so-called “symmetric” NZBC, i.e., boundary conditions where the modulus or the vector or matrix norm of the solution is the same asymptotically as \(x\rightarrow \pm \infty\). A number of recent papers have also explored the IST with non-symmetric NZBC both for the scalar focusing and defocusing NLS [33, 44, 66, 153], but the problem for vector and matrix NLS systems with non-symmetric NZBC is still outstanding.

With the developments of the last fifteen years, the IST has seen a marked revival, and a flourishing of papers and applications devoted in particular to the study of soliton interactions and more generally the long-time asymptotics behavior of solutions, and to the investigation of the integrable nature of modulational instability, integrable turbulence, soliton gases and connections to rogue waves. The goal of this topical review paper is to provide a survey of classical and more recent results on the IST for one-dimensional NLS systems on the line (\(-\infty<x<\infty\)) with certain symmetric non-zero boundary conditions, discuss some new developments and applications, and conclude by offering some perspectives about future directions on the subject. The plan of the paper is the following: in Sect. 2 we will review the IST for focusing and defocusing scalar NLS, and in Sects. 3 and 4 the IST for square matrix and vector NLS systems will be illustrated. In each section, a general overview, highlighting the main ideas, will be presented, while referring to the relevant existing literature for proofs and further details. Finally, Sect. 5 will be devoted to some recent theoretical developments and applications of the IST, including a succinct review of a modification of the standard IST for the focusing NLS equation with NZBC recently proposed by Bilman and Miller in [30] in order to deal with arbitrary-order poles and potentially severe spectral singularities in a rigorous way. Experimental applications of the IST will also be discussed, and possible future directions in the field highlighted.

We would like to stress that, since the focus of this review is on the IST and its applications, for economy of space the description of results achieved by the use of direct methods, and of experimental results on NLS-type systems and their soliton solutions will necessarily be limited, and the associated list of references may only touch upon some of the extensive body of valuable works on the subject. Also, in this review we will not discuss the periodic problem and the associated finite gap method, and we refer the reader to [92, 98] and references therein for the periodic NLS Cauchy problem. Lastly, even though they are outside the scope of this review, we want to briefly discuss below some other important topics related to the IST for NLS systems which have been the subject of investigation over the years.

Integrable discrete NLS systems. Both the NLS and the vector NLS equations admit integrable discretizations, i.e., finite difference approximations as differential-difference (continuous in time, discrete in space) equations which remain integrable as discrete systems. The most celebrated integrable discretizations of the scalar focusing and defocusing NLS equations are the so-called Ablowitz–Ladik (AL) equations [8, 9]. The IST for the focusing AL equation was developed in [8] for rapidly decaying potentials, and generalized to NZBC in [148, 154, 178]. The IST for the defocusing AL equation with NZBC was developed in [6, 179] under the assumption that the background \(Q_o\) be strictly less than 1, and generalized to arbitrarily large background \(Q_o\) in [141]. As shown in [139], when \(Q_o>1\) the defocusing AL becomes modulationally unstable, and admits breathers and rogue wave solutions which are the analog of those of its focusing counterpart [20, 21, 139, 148]. Integrable vector discretizations of the NLS have also been proposed [14, 89, 90, 173], and the IST for the focusing version of these coupled discrete systems has been studied in [13, 15, 173] for rapidly decaying potentials. The IST with NZBC for these systems is, however, an open problem.

Initial-boundary value problems. The application of the IST for the solution of initial-boundary value problems posed on the half-line or the segment for NLS systems was pioneered by Fokas, who developed a method based on the simultaneous spectral analysis of the Lax pair, as well as on an algebraic relation (known as the “global relation”) coupling the initial conditions with all boundary values. There is by now a vast literature on the subject, dealing with both linearizable and nonlinearizable boundary conditions. We refer the interested reader to [31, 34, 49, 76,77,78,79, 105, 127] and references therein for details. We mention, however, that there is little or no work in the literature in the case of nontrivial boundary conditions.

Nonlocal NLS systems. An integrable nonlocal version of the nonlinear Schrödinger equation, in which the nonlinear term has the form \(2\sigma q(x,t)q^*(-x,t)\) in place of \(2\sigma |q(x,t)|^2\), was introduced by Ablowitz and Musslimani in 2013 [10]. The nonlocal NLS equation is PT symmetric, i.e., such that the so-called self-induced potential is invariant under the combined action of parity and time reversal symmetry. Its Lax pair can be obtained from the NLS Lax pair by simply replacing \(\sigma q^*(x,t)\) with \(\sigma q^*(-x,t)\). The initial-value problem has been solved by IST both for decaying potentials and in the case of NZBC, and soliton solutions have been obtained, via the IST as well as by direct methods. Since then, several new reverse space-time and reverse time nonlocal nonlinear integrable NLS-type equations have been introduced, including generalizations to coupled, discrete, and multidimensional integrable nonlocal NLS systems, and, most recently, space-time nonlocal NLS-type equations with space and time shifts. We refer the reader to [11, 12] for a detailed account and additional references.

Fractional NLS equations. An integrable fractional NLS equation has been introduced by Ablowitz et al in 2022 as part of a new class of integrable fractional nonlinear evolution equations which describe dispersive transport in fractional media [4]. These equations have been constructed from nonlinear integrable equations using a mathematical process which relies on completeness of eigenfunctions, dispersion relations, and inverse scattering transform techniques. A fractional integrable discrete nonlinear evolution equation has been also proposed [5], and these works will certainly spur further investigations and applications.

2 IST for the Scalar NLS with NZBC

As mentioned above, the IST for the focusing NLS equation was developed in 1972 by Zakharov and Shabat under the assumption that the potential q(xt) is rapidly decaying as \(|x|\rightarrow \pm \infty\) [192], and this problem has been extensively studied in the literature over the years, both in the focusing and in the defocusing dispersion regimes (see, for instance, [7, 13, 16, 137] for detailed accounts of the IST in the rapidly decaying case).

The situation is quite different when one is interested in potentials that do not decay as \(|x|\rightarrow \infty\). This class of potentials is particularly important for the defocusing NLS, since it admits soliton solutions with NZBC, the so-called dark/gray solitons, which have the form:

$$\begin{aligned} q(x,t)=q_o e^{-2iq_o^2t}\left[ \cos \alpha +i \sin \alpha \tanh (q_o \sin \alpha (x+2q_o\cos \alpha t -x_o)\right] \end{aligned}$$
(2.1)

with \(q_o,\alpha ,x_o\) arbitrary real parameters. Dark soliton solutions are such that

$$\begin{aligned} q(x,t)\sim q_\pm (t)\equiv q_o e^{-2i\sigma q_o^2t+i\theta _\pm } \qquad \text {as } x\rightarrow \pm \infty , \qquad \sigma =1 \end{aligned}$$
(2.2)

and appear as localized dips of intensity \(q_o \sin \alpha\) on the background field \(q_o\). The IST for the defocusing NLS equation with NZBC was first studied by Zakharov and Shabat in 1973 [193]; the problem was subsequently clarified and generalized over the years in various works [22, 23, 87, 110, 111, 128], and a detailed study can be found in the monograph [74]. In [65] some rigorous results were presented on the functional class of nondecaying potentials where the direct and inverse scattering problems can be solved, and on the analyticity properties of eigenfunctions.

As to the focusing NLS with NZBC, although solutions over a nontrivial background (2.2) with \(\sigma =-1\) have been extensively investigated over the years [17,18,19, 124, 131, 144, 168], the only early study of the IST with NZBC was by Ma in 1979 [131], and the IST was not fully developed until 2014 [35]. The main reason why the IST for the focusing NLS with NZBC remained an outstanding problem for over 40 years is perhaps related to the fact that, unlike its defocusing counterpart, the scalar focusing NLS exhibits modulational instability—also known as Benjamin-Feir instability in the context of water waves [26, 27]—namely the instability of a uniform background to long wavelength perturbations (see [191] for a thorough review of this topic). Although modulational instability might have been considered an obstacle to the development of the IST, recent studies have shown that on the contrary the IST provides an effective tool to study the nonlinear stage of modulational instability (see Sect. 5 for some additional details).

We review below the IST for both focusing and defocusing NLS with symmetric NZBC at space infinity, mostly following [35, 65, 151] though using slightly different notations and normalizations, for the sake of clarity and self-consistency among different sections of this work. Instead of (2.2), it is convenient to take both asymptotic values \(q_\pm\) to be time-independent. This is easily achieved by performing a trivial gauge transformation to remove the constant rotation due to the background field \(q_o\) (i.e., letting \({\tilde{q}}(x,t)=e^{2i\sigma q_o^2t}q(x,t)\)), thus obtaining an NLS equation in the form:

$$\begin{aligned} iq_t+q_{xx}-2\sigma (|q|^2-q_o^2)q=0 \end{aligned}$$
(2.3)

where \(\sigma\) distinguishes, as before, between the focusing (\(\sigma =-1\)) and defocusing (\(\sigma =1\)) regimes. The corresponding boundary conditions become

$$\begin{aligned} q(x,t)\sim q_\pm =q_o e^{i\theta _\pm } \qquad x\rightarrow \pm \infty \end{aligned}$$
(2.4)

with time-independent \(\theta _\pm \in {\mathbb {R}}\).

Lax pair, Riemann surface and uniformization coordinate. The Lax pair for the NLS equation in the form (2.3) can be written as

$$\begin{aligned} \varphi _x={\textbf{U}} \varphi , \qquad \varphi _t={\textbf{V}} \varphi \end{aligned}$$
(2.5a)

where \(\varphi =\varphi (x,t,k)\) is a \(2\times 2\) matrix and

$$\begin{aligned} {\textbf{U}}= & {} -ik\sigma _3+Q, \qquad {\textbf{V}}=i\sigma _3 (-2k^2+ Q_x-Q^2+\sigma q_o^2)+2kQ, \end{aligned}$$
(2.5b)
$$\begin{aligned} \sigma _3= & {} \begin{pmatrix}1 &{} 0 \\ 0 &{} -1 \end{pmatrix}, \qquad Q(x,t)=\begin{pmatrix} 0 &{} q(x,t) \\ \sigma q^*(x,t) &{} 0 \end{pmatrix}, \end{aligned}$$
(2.5c)

and the asterisk denotes complex conjugation hereafter.

The scattering problem in (2.5) for \(\sigma =-1\) is often referred to as Zakharov–Shabat or AKNS scattering problem, while when \(\sigma =1\) the corresponding operator is also known as the Dirac operator (with nonzero rest mass). With the NZBC (2.4), the asymptotic scattering problems \(\varphi _x={\textbf{U}}_\pm \varphi\), \(\varphi _t={\textbf{V}}_\pm \varphi\), obtained by replacing Q with

$$\begin{aligned} Q_\pm = \lim _{x\rightarrow \pm \infty } Q(x, t) \end{aligned}$$
(2.6)

in (2.5), have eigenvalues \(\pm i\lambda\) where \(\lambda ^2=k^2-\sigma q_o^2\). It is then convenient to think of the variable k as belonging to a Riemann surface \({\mathbb {K}}\) consisting of a sheet \({\mathbb {K}}_I\) and a sheet \({\mathbb {K}}_{II}\) each coinciding with the complex plane cut along the semilines \((\infty ,-q_o] \cup [q_o,\infty )\) if \(\sigma =1\), and along the segment \([-iq_o,iq_o]\) on the imaginary k-axis if \(\sigma =-1\), with edges glued in such a way that \(\lambda (k)\) is continuous through the cut. The branch points are the values of k for which \(k^2-\sigma q_o^2=0\), i.e., \(k= \pm \sqrt{\sigma } q_o\), with \(\sqrt{\sigma }=i\) in the focusing case, and \(\sqrt{\sigma }=1\) in the defocusing case. Following a well-established procedure first introduced in [74], we define a uniformization variable:

$$\begin{aligned} z = k + \lambda \, \end{aligned}$$
(2.7a)

whose inverse transformation is given by

$$\begin{aligned} k = \frac{1}{2}\,(z+\sigma q_o^2/z),\qquad \lambda = \frac{1}{2}\,(z-\sigma q_o^2/z). \end{aligned}$$
(2.7b)

With these definitions, in the defocusing case the branch cut on either sheet is mapped onto the real z axis; the two sheets \({\mathbb {K}}_I\) and \({\mathbb {K}}_{II}\) of the Riemann surface are, respectively, mapped onto the upper and lower half-planes of the complex z-plane; a neighborhood of \(k=\infty\) on either sheet is mapped onto a neighborhood of \(z=\infty\) or \(z=0\) depending on the sign of \(\mathop {\textrm{Im}}\nolimits k\).

Conversely, in the focusing case the branch cut on either sheet is mapped onto the circle \(C_o\) (i.e., the circle centered at \(z=0\) and of radius \(q_o\) in the complex z-plane), \({\mathbb {K}}_I\) is mapped onto the exterior of \(C_o\); \({\mathbb {K}}_{II}\) is mapped onto the interior of \(C_o\); \(z(\infty _I) = \infty\) and \(z(\infty _{II}) = 0\).

As a consequence, in the defocusing case, \(\mathop {\textrm{Im}}\nolimits \lambda >0\) in the upper-half plane (UHP) and \(\mathop {\textrm{Im}}\nolimits \lambda <0\) in the lower-half plane (LHP) of z:

$$\begin{aligned} \sigma =1: \qquad D^+=\{z\in {\mathbb {C}}: \mathop {\textrm{Im}}\nolimits z >0\},\qquad D^-=\{z\in {\mathbb {C}}: \mathop {\textrm{Im}}\nolimits z <0\}. \end{aligned}$$
(2.8a)

On the other hand, in the focusing case \(\mathop {\textrm{Im}}\nolimits \lambda\) is not sign-definite in either half-plane; instead, one has \(\mathop {\textrm{Im}}\nolimits \lambda >0\) in \(D^+\) and \(\mathop {\textrm{Im}}\nolimits \lambda <0\) in \(D^-\), where, in this case:

$$ \sigma= -1: \quad D^+=\{z\in {\mathbb{C}}:(|z|^2-q_o^2)\mathop {\textrm{Im}}\nolimits z >0\},\quad D^-= \{z\in {\mathbb {C}}:(|z|^2-q_o^2)\mathop {\textrm{Im}}\nolimits z <0\}. $$
(2.8b)

The complex z-plane and the two domains \(D^\pm\) are shown in Fig. 1 (defocusing case on the left, focusing case on the right).

Fig. 1
figure 1

Left/Right: The complex z-plane, showing the regions \(D^\pm\) where \(\mathop {\textrm{Im}}\nolimits \lambda >0\) (grey) and \(\mathop {\textrm{Im}}\nolimits \lambda <0\) (white), respectively, in the defocusing/focusing case. Also shown in the figures are the oriented contours for the Riemann–Hilbert problem (red), and the symmetries of the discrete spectrum of the scattering problem

As will be discussed below, it is precisely the sign of \(\mathop {\textrm{Im}}\nolimits \lambda\) that determines the regions of analyticity of the eigenfunctions. Note that with some abuse of notation we will rewrite all the k dependence as dependence on z wherever appropriate.

Jost solutions and analyticity. The continuous spectrum \(\Sigma _k\) consists of all values of k (on either sheet) such that \(\lambda (k)\in {\mathbb {R}}\); i.e., \(\Sigma _k= {\mathbb {R}}\cup i[-q_o,q_o]\) in the focusing case, and \(\Sigma _k={\mathbb {R}}{\setminus } (-q_o,q_o)\) in the defocusing case. The corresponding sets in the complex z-plane are \(\Sigma _z= {\mathbb {R}}\cup C_o\) and \(\Sigma _z={\mathbb {R}}\), respectively, \(C_o\) being the circle of radius \(q_o\) centered at the origin (see Fig. 1). Hereafter we will omit the subscripts on \(\Sigma\), as the intended meaning will be clear from the context. For all \(z\in \Sigma\), we can now define the Jost eigenfunctions \(\Phi (x,t,z)\) and \(\Psi (x,t,z)\) as the simultaneous solutions of both parts of the Lax pair (2.5) such that

$$\begin{aligned} \Phi (x,t,z)& \equiv\, (\phi (x,t,z)\,\, {\bar{\phi }}(x,t,z)) \\& =\, \left[ E_-(z)+o(1)\right] \textrm{e}^{-i\theta (x,t,z)\sigma _3} \qquad \textrm{as}~x\rightarrow -\infty , \end{aligned}$$
(2.9a)
$$\begin{aligned} \Psi (x,t,z)& \equiv\, ({\bar{\psi }} (x,t,z)\,\, \psi (x,t,z)) \\& = \,\left[ E_+(z)+o(1)\right] \textrm{e}^{-i\theta (x,t,z)\sigma _3} \qquad \textrm{as}~x\rightarrow +\infty , \end{aligned}$$
(2.9b)

where

$$\begin{aligned} \theta (x,t,z) = \lambda (z)(x+2k(z)t), \qquad E_\pm (z)=I_2-\frac{i}{z}\sigma _3 Q_\pm , \end{aligned}$$
(2.10)

\(E_\pm (z)\) are the (simultaneous) eigenvectors of the asymptotic Lax matrices \({\textbf{U}}_\pm =\lim _{x\rightarrow \pm \infty }{\textbf{U}}(x,t)\) and \({\textbf{V}}_\pm =\lim _{x\rightarrow \pm \infty }{\textbf{V}}(x,t)\), and \(\phi (x,t,z)\) and \({\bar{\phi }}(x,t,z)\) (resp. \({\bar{\psi }}(x,t,z)\) and \(\psi (x,t,z)\)) are the column vectors of the \(2\times 2\) matrix solutions \(\Phi (x,t,z)\) (resp. \(\Psi (x,t,z)\)). It is convenient to introduce modified eigenfunctions:

$$\begin{aligned} \left( M(x,t,z)\,\, {\bar{M}}(x,t,z)\right)& = \,\Phi (x,t,z)e ^{i\theta (x,t,z)\sigma _3} \\ \left( {\bar{N}}(x,t,z)\,\, N(x,t,z)\right)& = \,{} \Psi (x,t,z)e ^{i\theta (x,t,z)\sigma _3} , \end{aligned}$$
(2.11)

which approach \(E_\pm (z)\) as \(x\rightarrow \pm \infty\).

Following an established procedure, one can then write down integral Volterra equations for the modified eigenfunctions and prove that, under mild integrability conditions on the potential, the modified eigenfunctions M(xtz) and N(xtz) can be analytically extended in the complex z-plane where \(\mathop {\textrm{Im}}\nolimits \lambda (z)> 0\), and \({\bar{M}}(x,t,z)\) and \({\bar{N}}(x,t,z)\) can be analytically extended in the complex z-plane where \(\mathop {\textrm{Im}}\nolimits \lambda (z)< 0\). Let us denote by \(L^{1,s} ({\mathbb {R}})\) the complex Banach space of all measurable functions f(x) for which \((1 + |x|)^s f (x)\in L^{1}({\mathbb {R}})\) for \(s=0,1\). Then the following theorem holds.

Theorem 2.1

If \(q(x,t)-q_+\in L^1([x_+,+\infty ))\) and \(q(x,t)-q_-\in L^1((-\infty ,x_-])\) as functions of x for all \(t\geqslant 0\), for some \(x_\pm \in {\mathbb {R}}\), and \(D^\pm\) are defined in (2.8) for the defocusing/focusing cases (\(\sigma =\pm 1\), respectively), then for the modified eigenfunctions the following holds: M(xtz) and N(xtz) are analytic functions of z for \(z\in D^+\), and they are continuous up to \(\partial D^+ {\setminus }\left\{ \pm \sqrt{\sigma } q_o,0\right\}\); \({\bar{M}}(x,t,z)\) and \({\bar{N}}(x,t,z)\) are analytic functions of z for \(z\in D^-\), and they are continuous up to \(\partial D^- {\setminus }{\left\{ \pm \sqrt{\sigma } q_o,0\right\} }\). Furthermore, if \(q(x,t)-q_+\in L^{1,1}([x_+,+\infty ))\) and \(q(x,t)-q_-\in L^{1,1}((-\infty ,x_-])\) as functions of x for all \(t\geqslant 0\), for some \(x_\pm \in {\mathbb {R}}\), then the modified eigenfunctions are such that: M(xtz) and N(xtz) are analytic functions of z for \(z\in D^+\), and they are continuous up to \(\partial D^+{\setminus } \left\{ 0\right\}\); \({\bar{M}}(x,t,z)\) and \({\bar{N}}(x,t,z)\) are analytic functions of z for \(z\in D^-\), and they are continuous up to \(\partial D^-\setminus \left\{ 0\right\}\). Namely, in this case the eigenfunctions are also continuous up to the branch points \(\pm \sqrt{\sigma } q_o\).

The results are proven using standard Neumann series representations for the solutions of the integral equations for the eigenfunctions [35, 65, 74]. Similar analyticity properties for the corresponding columns of \(\Phi (x,t,z)\) and \(\Psi (x,t,z)\) follow trivially from the above theorem. Note that the Jost eigenfunctions in general have a pole at \(z=0\) (cf. (2.10)). The behavior of eigenfunctions (and scattering coefficients, see below) at \(z=0\) and as \(z\rightarrow \infty\) is analyzed separately, as it is crucial to properly formulate the inverse problem. Finally, we should clarify that we refer to the points \(z=\pm \sqrt{\sigma } q_o\) as branch points, even though there is no branching in the complex z-plane, simply as a reminiscence of the fact that they are the images of the actual branch points \(k=\pm \sqrt{\sigma } q_o\) for \(\lambda =\sqrt{k^2-\sigma q_o^2}\).

Scattering coefficients. As a consequence of Jacobi’s formula, any matrix solution \(\varphi (x,t,z)\) of (2.5) satisfies \(\partial _x(\det \varphi ) = \partial _t(\det \varphi ) = 0\), given that both \({\textbf{U}}\) and \({\textbf{V}}\) in (2.5) are traceless. Thus, since for all \(z\in \Sigma\) one has \(\lim _{x\rightarrow -\infty }\Phi (x,t,z)\,\textrm{e}^{i\theta \sigma _3}=E_-\) and \(\lim _{x\rightarrow +\infty }\Psi (x,t,z)\,\textrm{e}^{i\theta \sigma _3}= E_+\), it follows that

$$\begin{aligned} \det \Phi (x,t,z)= \det \Psi (x,t,z)= \det E_\pm (z)= \gamma (z)\,\qquad x,t\in {\mathbb {R}},\qquad z\in \Sigma , \end{aligned}$$
(2.12)

where

$$\begin{aligned} \gamma (z)=1- \sigma q_o^2/z^2. \end{aligned}$$
(2.13)

Letting \(\Sigma _o\) denote the continuous spectrum minus the branch points, i.e., \(\Sigma _o = \Sigma \setminus \{\pm \sqrt{\sigma }q_o\}\), we then have that \(\forall z\in \Sigma _o\) both \(\Phi\) and \(\Psi\) are two fundamental matrix solutions of the scattering problem. Hence, there exists a constant \(2\times 2\) scattering matrix S(z) such that

$$\begin{aligned} \Phi (x,t,z) = \Psi (x,t,z)\,S(z), \qquad S(z)= \begin{pmatrix} a(z) &{} {\bar{b}}(z) \\ b(z) &{} {\bar{a}}(z) \end{pmatrix}, \qquad x,t\in {\mathbb {R}},\qquad z\in \Sigma _o. \end{aligned}$$
(2.14)

Note that since \(\Phi\) and \(\Psi\) are simultaneous solutions of both parts of the Lax pair, the entries of S(z) are independent of time, and the same is true for the norming constants introduced below. Moreover, (2.12) and (2.14) also imply that \(\det S(z) = 1\). One can express the scattering coefficients in terms of the columns of the Jost eigenfuctions as:

$$a(z) = \,\textrm{Wr}(\phi ,\psi )/\gamma (z), \qquad {\bar{a}}(z)=-\textrm{Wr}({\bar{\phi }},{\bar{\psi }})/\gamma (z),$$
(2.15a)
$$b(z)= -\textrm{Wr}(\phi ,{\bar{\psi }})/\gamma (z), \qquad {\bar{b}}(z)=\textrm{Wr}({\bar{\phi }},\psi )/\gamma (z),$$
(2.15b)

[\(\textrm{Wr}(f, g)\) denotes the Wronskian determinant of the column vectors fg], which show that a(z) and \({\bar{a}}(z)\) can be analytically continued respectively on \(D^+\) and \(D^-\) (and may have poles at the branch points \(z= \pm \sqrt{\sigma }q_o\) where \(\gamma (z)\) vanishes, cf. (2.13)). On the other hand, b(z) and \({\bar{b}}(z)\) are continuous for \(z\in \Sigma _o\), generically with poles at \(z=\pm \sqrt{\sigma }q_o\), but do no admit analytic continuation off \(\Sigma\).

Finally, for \(z\in \Sigma _o\) we can use (2.11) and (2.14) to obtain

$$\begin{aligned} M(x,t,z)/a(z)= \,& {} {\bar{N}}(x,t,z)+e^{2i\theta (x,t,z)}N(x,t,z)\rho (z), \end{aligned}$$
(2.16a)
$$\begin{aligned} {\bar{M}}(x,t,z)/{\bar{a}}(z)= \,& {} N(x,t,z)+e^{-2i\theta (x,t,z)}{\bar{N}}(x,t,z){\bar{\rho }}(z), \end{aligned}$$
(2.16b)

where M(xtz)/a(z) and \({\bar{M}}(x,t,z)/{\bar{a}}(z)\) are meromorphic in \(D^+\) and \(D^-\), and we introduced the reflection coefficients

$$\begin{aligned} \rho (z)=b(z)/a(z), \qquad {\bar{\rho }}(z)={\bar{b}}(z)/{\bar{a}}(z), \qquad z\in \Sigma _o. \end{aligned}$$
(2.17)

We omit for brevity the details of the behavior of eigenfunctions, scattering coefficients and reflection coefficients at \(z=\pm \sqrt{\sigma }q_o\), they can be found in [35, 65, 74].

Symmetries. In solving an initial-value problem by IST, one has to take into account that symmetries in the potentials in the Lax pair induce symmetries in the eigenfunctions, which in turn induce corresponding symmetries in the scattering data. The symmetries for the IST with NZBC are complicated by the fact that \(\lambda (k)\) changes sign from one sheet of the Riemann surface to the other. In terms of the uniformization variable z, one needs to consider the following involutions: 1. \(z\mapsto z^*\) (UHP/LHP), implying  \((k,\lambda )\mapsto (k^*,\lambda ^*)\) (same sheet); 2. \(z\mapsto \sigma q_o^2/z\) (outside/inside \(C_o\)), implying  \((k,\lambda )\mapsto (k,-\lambda )\) (opposite sheets). Both these transformations correspond to symmetries of the scattering problem; the first one is related to the usual conjugation symmetry in the potential, \(Q^\dagger =\sigma Q\), while the second one is a direct consequence of the branching in the plane of the scattering parameter k, or, equivalently, the fact that the uniformization variable z provides a double covering of the k plane.

In the defocusing case, one can easily check that the Jost eigenfunctions are such that

$$\begin{aligned} \Phi (x,t,z)=\sigma _1 \Psi ^* (x,t,z^*) \sigma _1, \qquad \sigma _1=\begin{pmatrix} 0 &{} 1 \\ 1 &{} 0 \end{pmatrix}, \end{aligned}$$
(2.18)

for all z in \(D^+\) for the first column, and z in \(D^-\) for the second column. Correspondingly, from (2.14) it follows

$$\begin{aligned} a(z)={\bar{a}}^*(z^*) \quad z\in D^+, \qquad {\bar{b}}(z)=b^*(z) \quad z\in \Sigma . \end{aligned}$$
(2.19)

Moreover, \(\det S(z)=1\) yields

$$\begin{aligned} |a(z)|^2=1+|b(z)|^2\geqslant 1 \qquad z\in \Sigma _o. \end{aligned}$$
(2.20)

In the focusing case, the corresponding symmetries for eigenfunctions and scattering coefficients read:

$$\begin{aligned} \Phi (x,t,z)=\sigma _2 \Psi ^* (x,t,z^*) \sigma _2, \qquad \sigma _2=\begin{pmatrix} 0 &{} -i \\ i &{} 0 \end{pmatrix}, \end{aligned}$$
(2.21)

again with z in \(D^+\) for the first column, and z in \(D^-\) for the second column, and

$$\begin{aligned} a(z)={\bar{a}}^*(z^*) \quad z\in D^+, \qquad {\bar{b}}(z)=-b^*(z^*) \quad z\in \Sigma . \end{aligned}$$
(2.22)

[Note that in the focusing case the continuous spectrum \(\Sigma\) is a superset of the real axis, hence \(z^*\) is necessary in the above symmetry of b(z) when \(z\notin {\mathbb {R}}\).]

The second involution, \(z\mapsto \sigma q_o^2/z\), corresponds to a reflection with respect to the circle \(C_o\). The symmetries for the eigenfunctions are given by

$$\begin{aligned} \Phi (x,t,z)= & {} -\frac{1}{iz} \Phi (x,t,\sigma q_o^2/z)\sigma _3 Q_-, \\ \Psi (x,t,z)= \,& {} \frac{1}{iz} \Psi (x,t,\sigma q_o^2/z)\sigma _3 Q_+, \end{aligned}$$
(2.23)

[as before, z in \(D^+\) for the first column, and z in \(D^-\) for the second column] and in turn this implies

$$\begin{aligned} a(\sigma q_o^2/z)= \,& {} \frac{q_+q_-^*}{q_o^2} {\bar{a}}(z) \quad \text {for }z\in D^-,\qquad \\ {\bar{b}}(\sigma q_o^2/z)= & {} -\sigma \frac{q_+q_-}{q_o^2} b(z)\,\quad \quad \text {for }z\in \Sigma . \end{aligned}$$
(2.24)

Finally, the above relations yield the corresponding symmetries for the reflection coefficients:

$$\begin{aligned} {\bar{\rho }}(z)=\rho ^*(z^*), \qquad {\bar{\rho }}(\sigma q_o^2/z) = -\sigma \frac{q_+^2}{q_o^2} \rho (z) \qquad \forall z\in \Sigma . \end{aligned}$$
(2.25)

Discrete spectrum. The discrete spectrum is the set of all \(z\in {\mathbb {C}}{\setminus } \Sigma\) such that \(a(z)=0\) or \({\bar{a}}(z)=0\), values for which the Jost eigenfunctions are in \(L^2({\mathbb {R}})\).

In the defocusing case, (2.20) shows that a(z) and \({\bar{a}}(z)\) cannot vanish for \(z\in \Sigma _o\). Moreover, the zeros of a(k) and \({\bar{a}}(k)\) as functions of the original spectral variable k are real (a consequence of the formal self-adjointness of the scattering problem when \(\sigma =1\)), and they are also simple (see [74] for details). In [65] it is also shown that if \(q-q_\pm \in L^{1,4}({\mathbb {R}}^\pm )\), then a(k) and \({\bar{a}}(k)\) have a finite number of zeros, all of which belong to the spectral gap \(k\in (-q_o,q_o)\). In terms of the uniformization variable, this means that one can generically assume a finite number of pairs of discrete eigenvalues on the circle \(C_o\setminus \left\{ \pm q_o\right\}\), i.e.

$$\begin{aligned} {\mathcal {Z}}=\left\{ \zeta _j,\zeta _j^*\right\} _{j=1}^n, \qquad \zeta _j=q_o e^{i\alpha _j}, \quad 0<\alpha _j<\pi . \end{aligned}$$

[The second symmetry (2.24) does not produce additional discrete eigenvalues in this case, since for \(\zeta \in C_o\) one has \(q_o^2/\zeta =\zeta ^*\).] For all \(j=1,\dots ,n\), the Wronskian representations (2.15a) yield

$$\begin{aligned} \phi (x,t,\zeta _j)=b_j\psi (x,t,\zeta _j), \qquad {\bar{\phi }}(x,t,\zeta _j^*)={\bar{b}}_j{\bar{\psi }}(x,t,\zeta _j^*), \end{aligned}$$
(2.26)

for some constants \(b_j,{\bar{b}}_j\in {\mathbb {C}}\), with \(b_j^*={\bar{b}}_j\) due to the symmetries (2.18). The second symmetry (2.23) also implies

$$\begin{aligned} b_j^*=-\frac{q_-}{q_+^*}b_j, \qquad {\bar{b}}_j=-\frac{q_+^*}{q_-}{\bar{b}}_j. \end{aligned}$$

In turn, these provide the residue relations which will be needed in the inverse problem:

$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z=\zeta _j}\left[ M/a\right]= \,& {} C_je^{-2i\Omega (x,t,\zeta _j)}N(x,t,\zeta _j), \qquad \\ \mathop {\textrm{Res}}\limits _{z=\zeta _j^*}\left[ {\bar{M}}/{\bar{a}}\right]= \,& {} {\bar{C}}_je^{-2i\Omega (x,t,\zeta _j^*)}{\bar{N}}(x,t,\zeta _j^*), \end{aligned}$$
(2.27)

where the norming constants \(C_j=b_j/a^\prime (\zeta _j)\) and \({\bar{C}}_j={\bar{b}}_j/{\bar{a}}^\prime (\zeta _j^*)\) satisfy the symmetry relations

$$\begin{aligned} {\bar{C}}_j=C_j^*, \qquad C_j^*=\frac{q_+}{q_+^*}C_j. \end{aligned}$$
(2.28)

Note that under the assumptions of Theorem 2.1 the Jost solutions are continuous at the branch points \(\pm q_o\), while the scattering coefficients generically have simple poles when \(z=\pm q_o\), unless the columns of \(\Phi (x, t, z)\) and \(\Psi (x,t,z)\) become linearly dependent at either \(z=\pm q_o\) or both. In this case, a(z) and \({\bar{a}}(z)\) are nonsingular near the corresponding branch point, and the point \(z=q_o\) (or \(z=-q_o\)) is called a virtual level [74].

The discrete spectrum of the focusing NLS is more involved: discrete eigenvalues can be located anywhere in the complex z-plane, they are not necessarily simple, and one cannot a priori exclude the existence of zeros of a(z) and \({\bar{a}}(z)\) for \(z\in \Sigma\) (the latter are often referred to as spectral singularities [195], or embedded eigenvalues). For decaying potentials, there exists a characterization of the location of (real) spectral singularities, as well as sufficient conditions on q(xt) that guarantee their absence (for instance, spectral singularities are absent in the case of single lobe potentials, and certain double and multiple lobe potentials [117, 119, 120]). On the other hand, there are potentials in the Schwartz class for which discrete eigenvalues can accumulate to spectral singularities, and spectral singularities themselves can accumulate on the continuous spectrum [195]. The characterization of spectral singularities for the focusing NLS with NZBC is an open problem, to the best of our knowledge. To date, what is known in the scalar case is that there are arbitrarily small perturbations of the constant background for which discrete eigenvalues are present, both in the focusing and in the defocusing case, although there is a class of perturbations of the constant background for which no discrete eigenvalues are present [32, 43] (see also [67] for some results on the location of the discrete eigenvalues of the defocusing Zakharov–Shabat system). In the following, we will assume a finite number of simple, discrete eigenvalues and no spectral singularities, i.e., in light of the symmetries discussed above, a discrete spectrum consisting of quartets of eigenvalues:

$$\begin{aligned} {\mathcal {Z}}=\left\{ z_j,-q_o^2/z_j^*,z_j^*,-q_o^2/z_j \right\} _{j=1}^{\mathcal {N}}, \end{aligned}$$

where without loss of generality we assume \(z_j\in D^+\) with \(\mathop {\textrm{Im}}\nolimits z_j>0\) for all \(j=1,\dots ,{\mathcal {N}}\) (cf. Fig. 1). As before, if \(a(z)=0\) at \(z=z_j\in D^+\), then \({\bar{a}}(z_j^*)=0\), and the eigenfunctions must be proportional, namely:

$$\begin{aligned} \phi (x,t,z_j)=b_j\psi (x,t,z_j), \qquad {\bar{\phi }}(x,t,z_j^*)={\bar{b}}_j {\bar{\psi }}(x,t,z_j^*). \end{aligned}$$
(2.29)

Consequently, one has

$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z=z_j}\left[ M/a\right]= \,& {} C_je^{-2i\Omega (x,t,z_j)}N(x,t,z_j), \qquad C_j=b_j/a^\prime (z_j), \end{aligned}$$
(2.30a)
$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z=z_j^*}\left[ {\bar{M}}/{\bar{a}}\right]= \,& {} {\bar{C}}_je^{-2i\Omega (x,t,z_j^*)}{\bar{N}}(x,t,z_j^*), \qquad {\bar{C}}_j={\bar{b}}_j/{\bar{a}}^\prime (z_j^*), \end{aligned}$$
(2.30b)

and the norming constants are such that \({\bar{C}}_j=-C_j\), due to the symmetries (2.21) and (2.22). Moreover, because of the second symmetry, from (2.23) and (2.24) it follows

$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z=-q_o^2/z_j^*}\left[ M/a\right]= \,& {} C_{{\mathcal {N}}+j}e^{-2i\Omega (x,t,-q_o^2/z_j^*)}N(x,t,z-q_o^2/z_j^*), \end{aligned}$$
(2.30c)
$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z=-q_o^2/z_j}\left[ {\bar{M}}/{\bar{a}}\right]= \,& {} {\bar{C}}_{{\mathcal {N}}+j}e^{-2i\Omega (x,t,-q_o^2/z_j)}{\bar{N}}(x,t,-q_o^2/z_j), \end{aligned}$$
(2.30d)

with

$$\begin{aligned} C_{{\mathcal {N}}+j}=(q_o/z_j^*)^2(q_-^*/q_-){\bar{C}}_j, \qquad {\bar{C}}_{{\mathcal {N}}+j}=(q_o/z_j)^2(q_-/q_-^*)C_j. \end{aligned}$$
(2.31)

Asymptotics of eigenfunctions and scattering data. The asymptotic behavior of the modified Jost eigenfunctions for \(z\rightarrow \infty\) and \(z\rightarrow 0\) can be obtained via standard WKB expansions, or using integration by parts on the corresponding integral equations (cf. [35, 65, 151] for details). The asymptotic expansion as \(z\rightarrow \infty\) is given by

$$\begin{aligned} M(x,t,z)= & {} \begin{pmatrix} 1 +\displaystyle {\frac{i\sigma }{z}}\mathop {\int }\limits _{-\infty }^x\left[ |q(y,t)|^2-q_o^2 \right] dy\\ i\sigma q^*(x,t)/z \end{pmatrix} \\{} & {} +O(1/z^2) \qquad z\rightarrow \infty , z\in D^+ \end{aligned}$$
(2.32a)
$$\begin{aligned} {\bar{M}}(x,t,z)= & {} \begin{pmatrix} -iq(x,t)/z \\ 1-\displaystyle {\frac{i\sigma }{z}}\mathop {\int }\limits _{-\infty }^x\left[ |q(y,t)|^2-q_o^2 \right] dy \end{pmatrix} \\{} & {} +O(1/z^2) \qquad z\rightarrow \infty , z\in D^- \end{aligned}$$
(2.32b)

and

$$\begin{aligned} {\bar{N}}(x,t,z)= & {} \begin{pmatrix} 1 +\displaystyle {\frac{i\sigma }{z}}\mathop {\int }\limits _x^{\infty }\left[ |q(y,t)|^2-q_o^2 \right] dy\\ i\sigma q^*(x,t)/z \end{pmatrix} \\{} & {} +O(1/z^2) \qquad z\rightarrow \infty , z\in D^- \end{aligned}$$
(2.32c)
$$\begin{aligned} N(x,t,z)= & {} \begin{pmatrix} -iq(x,t)/z \\ 1-\displaystyle {\frac{i\sigma }{z}}\mathop {\int }\limits ^{\infty }_x\left[ |q(y,t)|^2-q_o^2 \right] dy \end{pmatrix} \\{} & {} +O(1/z^2) \qquad z\rightarrow \infty , z\in D^+. \end{aligned}$$
(2.32d)

Similarly, asymptotics as \(z\rightarrow 0\) in the appropriate region \(D^\pm\) yields

$$\begin{aligned} M(x,t,z)= & {} \begin{pmatrix} q(x,t)q_-^*/q_o^2+O(z)\\ i\sigma q_-^*/z+O(1) \end{pmatrix}, \qquad \\ {\bar{M}}(x,t,z)= & {} \begin{pmatrix} -iq_-/z+O(1)\\ q^*(x,t)q_-/q_o^2+O(z) \\ \end{pmatrix}, \end{aligned}$$
(2.33a)
$$\begin{aligned} {\bar{N}}(x,t,z)= & {} \begin{pmatrix} q(x,t)q_+^*/q_o^2+O(z)\\ i\sigma q_+^*(x,t)/z+O(1) \end{pmatrix}, \qquad \\ N(x,t,z)= & {} \begin{pmatrix} -iq_+/z+O(1)\\ q^*(x,t)q_+/q_o^2+O(z) \end{pmatrix}. \end{aligned}$$
(2.33b)

The above expansions can be used to reconstruct the scattering potential q(xt) from the solution of the inverse problem for the eigenfunctions.

Finally, inserting the above asymptotic expansions for the Jost eigenfunctions into (2.14) one shows that, as \(z\rightarrow \infty\) in the appropriate regions of the complex z-plane,

$$\begin{aligned} S(z)= I_{2} + O(1/z). \end{aligned}$$
(2.34)

The above asymptotics holds with \(\mathop {\textrm{Im}}\nolimits z\geqslant 0\) and \(\mathop {\textrm{Im}}\nolimits z\leqslant 0\) for a(z) and \({\bar{a}}(z)\), respectively, and with \(z\in \Sigma\) for b(z) and \({\bar{b}}(z)\). Similarly, one can show that as \(z\rightarrow 0\)

$$\begin{aligned} S(z)= \begin{pmatrix} q_+/q_-^* &{} 0 \\ 0 &{} q_+^*/q_-^*\end{pmatrix} + O(z), \end{aligned}$$
(2.35)

where again the asymptotics for the diagonal entries of S(z) can be extended to \(D^+\) for a(z) and to \(D^-\) for \({\bar{a}}(z)\), while the asymptotics for the off-diagonal entries hold for \(z\in \Sigma\).

Inverse problem. The inverse problem consists of constructing a map from the scattering data, that is: (i) the reflection coefficient \(\rho (z)\) on the continuous spectrum \(\Sigma\), (ii) the eigenvalues \({\mathcal {Z}} = \{z_j,z_j^*,\sigma q_o^2/z_j,\sigma q_o^2/z_j^*\}_{j=1}^n\), and (iii) the corresponding associated norming constants \(\{C_j,{\bar{C}}_j,{\hat{C}}_j,\hat{{\bar{C}}}_j \}_{j=1}^n\); back to the potential q(xt). In the IST approach, one first uses the scattering data to recover the modified eigenfunctions, then one recovers the potential q(xt) in terms of large asymptotic behavior in the spectral parameter of said eigenfunctions. In the direct problem, conditions on q(xt) are established such that the modified eigenfunctions N(xtz) and \({\bar{N}} (x,t,z)\) exist and are analytic as functions of the scattering parameter z in the regions \(D^\pm\), respectively. Similarly, under the same conditions on the potentials, the modified eigenfunctions M(xtz)/a(z) and \({\bar{M}}(x,t,z)/{\bar{a}}(z)\) are meromorphic functions of z in the regions \(D^\pm\), respectively. Therefore, in the inverse problem one assumes that the (unknown) modified eigenfunctions are sectionally meromorphic. With this assumption, the equations that relate the eigenfunctions on the continuous spectrum \(\Sigma\) can be considered to be the jump conditions of a Riemann–Hilbert problem across the contour \(\Sigma\). Indeed, from (2.16) one can write:

$$\begin{aligned} \mu ^-(x,t,z)=\mu ^+(x,t,z)(I_2-G(x,t,z)), \qquad \forall z\in \Sigma \end{aligned}$$
(2.36)

for the sectionally meromorphic functions

$$\begin{aligned} \mu ^+(x,t,z)=\left( M/a,N \right) , \qquad \mu ^-(x,t,z)=\left( {\bar{N}},{\bar{M}}/{\bar{a}} \right) , \end{aligned}$$
(2.37)

with a jump matrix

$$\begin{aligned} G(x,t,z)=\begin{pmatrix} 0 &{} \sigma e^{-2i\theta (x,t,z)}\rho ^*(z^*) \\ -e^{2i\theta (x,t,z)}\rho (z) &{} \sigma \rho (z)\rho ^*(z^*) \end{pmatrix}, \end{aligned}$$
(2.38)

which only depends on the reflection coefficient \(\rho (z) = b(z)/a(z)\). The above equations define a matrix, multiplicative, homogeneous Riemann–Hilbert problem (RHP). In order to complete the formulation of the RHP one needs a normalization condition, which in this case is the asymptotic behavior of \(\mu ^\pm\) as \(z\rightarrow \infty\). Recalling the asymptotic behavior of the Jost eigenfunctions and scattering coefficients (2.33) and (2.35), it is easy to check that

$$\begin{aligned} \mu ^\pm = I_{2} + O(1/z),\qquad z\rightarrow \infty . \end{aligned}$$
(2.39a)

On the other hand,

$$\begin{aligned} \mu ^\pm = -(i/z)\sigma _3\,Q_+ + O(1),\qquad z\rightarrow 0. \end{aligned}$$
(2.39b)

After regularization, to account for the poles at \(z = 0\) and at the discrete eigenvalues using the residue conditions (2.27) and (2.30), one can use Cauchy projectors to obtain a system of coupled integral equations for the columns of \(\mu ^\pm\). Specifically, for the analytic columns one can write

$$\begin{aligned} N(x,t,z)= & {} \begin{pmatrix} - iq_+/z \\ 1 \end{pmatrix} + \sum _{j=1}^{J}\frac{\textrm{e}^{-2i\theta (x,t,\zeta _j^*)}}{z-\zeta _j^*}\,{\bar{N}}(x,t,\zeta _j^*)\, {\bar{C}}_j \\{} & {} - \frac{1}{2\pi i}\mathop {\int }\limits _\Sigma \frac{e^{-2i\theta (x,t,\zeta )}{\bar{N}}(x,t,\zeta ){\bar{\rho }}(\zeta )}{\zeta -(z+i0)}\,\textrm{d}\zeta , \end{aligned}$$
(2.40a)
$$\begin{aligned} {\bar{N}}(x,t,z)= & {} \begin{pmatrix} 1 \\ i\sigma q_+^\dagger /z \end{pmatrix} + \sum _{j=1}^{J}\frac{\textrm{e}^{2i\theta (x,t,\zeta _j)}}{z-\zeta _j}\,N(x,t,\zeta _j)\, C_j \\{} & {} + \frac{1}{2\pi i}\mathop {\int }\limits _\Sigma \frac{e^{2i\theta (x,t,\zeta )}N(x,t,\zeta )\rho (\zeta )}{\zeta -(z-i0)}\,\textrm{d}\zeta , \end{aligned}$$
(2.40b)

where \(\sigma =1\), \(J={\mathcal {N}}\) and \(\Sigma ={\mathbb {R}}\) in the defocusing case, while in the focusing case \(\sigma =-1\), \(J=2{\mathcal {N}}\), \(\Sigma ={\mathbb {R}}\cup [-iq_o,q_o]\) and we have labelled discrete eigenvalues and associated norming constants so that: \(\zeta _j=z_j\), \(\zeta _{n+j}=-q_o^2/z_j^*\), \(C_{n+j}=(q_o/z_j)^2(q_-/q_-^*)C_j\) for \(j=1,\dots ,n\) (cf. (2.31)). In turn, comparing the behavior as \(z\rightarrow \infty\) of the second component of N(xtz) in (2.40) with (2.32) yields the reconstruction formula

$$\begin{aligned} q(x,t)= \,& {} q_+ +\sigma \sum _{j=1}^J C_j^* N_2^*(x,t,\zeta ^*_j)e^{-2i\theta (x,t,\zeta _i)} \\{} & {} +\frac{\sigma }{2\pi }\mathop {\int }\limits _\Sigma \rho ^*(\zeta )N_2^*(x,t,\zeta )e^{-2i\theta (x,t,\zeta )}\textrm{d}\zeta , \end{aligned}$$
(2.41)

where \(N_2(x,t,z)\) denotes the second component of N(xtz).

Finally, 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 D^+\):

$$\begin{aligned} a(z)=\prod _{j=1}^J\frac{z-\zeta _j}{z-\zeta _j^*}\exp {\left[ -\frac{1}{2\pi i} \mathop {\int }\limits _\Sigma \frac{\log (1-\sigma \rho (z)\rho ^*(\zeta ^*))}{\zeta -z}d\zeta \right] }, \end{aligned}$$
(2.42)

and recalling that \(a(z)\rightarrow q_+/q_-\) as \(z\rightarrow 0\), we conclude that the potential satisfies

$$\begin{aligned} \frac{q_+}{q_-}=\prod _{j=1}^J\frac{\zeta _j}{\zeta _j^*} \, \exp {\left[ -\frac{1}{2\pi i} \mathop {\int }\limits _\Sigma \frac{\log (1-\sigma \rho (z)\rho ^*(\zeta ^*))}{z}d\zeta \right] }, \end{aligned}$$
(2.43)

which is often referred to as the \(\theta\)-condition [74]. As before, in the above equations \(\sigma =1\), \(J={\mathcal {N}}\) and \(\Sigma ={\mathbb {R}}\) in the defocusing case, while in the focusing case \(\sigma =-1\), \(\Sigma ={\mathbb {R}}\cup [-iq_o,q_o]\), \(J=2{\mathcal {N}}\) and \(\zeta _j=z_j\), \(\zeta _{{\mathcal {N}}+j}=-q_o^2/z_j^*\).

We mention for completeness that the inverse problem can be equivalently formulated in terms of the so-called Marchenko (or Gelfand–Levitan–Marchenko, GLM) integral equations. First, one assumes triangular representations for the Jost solutions of the form:

$$\begin{aligned} {\bar{\psi }}(x,t,z)= \,& {} \left[ e^{i\theta (x,t,z)}I_2+\mathop {\int }\limits _x^{\infty }{\textbf{K}}(x,s;t)e^{i\theta (s,t,z)}ds\right] E_{+,1}(z), \end{aligned}$$
(2.44a)
$$\begin{aligned} \psi (x,t,z)= \,& {} \left[ e^{-i\theta (x,t,z)}I_2+\mathop {\int }\limits _x^{\infty }{\textbf{K}}(x,s;t)e^{-i\theta (s,t,z)}ds\right] E_{+,2}(z), \end{aligned}$$
(2.44b)

where \(E_{+,j}(z)\) are the columns of the eigenvector matrix \(E_+(z)\) in (2.10), and \({\textbf{K}}(x,s;t)\) is a \(2\times 2\) matrix often referred to as Marchenko kernel. The Marchenko kernel is then shown to satisfy the integral equation:

$$\begin{aligned} {\textbf{K}}(x,s;t)+{\textbf{G}}(x+y;t)+\mathop {\int }\limits _x^\infty {\textbf{K}}(x,s;t){\textbf{G}}(s+y;t)ds=0 \end{aligned}$$
(2.45a)

where \({\textbf{G}}(x+y;t)\) is completely determined by the scattering data (see [65, 74, 193] for additional details, keeping in mind that a different normalization is used for the eigenfunctions; also, [193] does not make use of the uniformization variable z). Both \({\textbf{K}}(x,s;t)\) and \({\textbf{G}}(x+y;t)\) satisfy the same symmetries (2.18) or (2.21) as the Jost eigenfunctions. The potential Q(xt) can then be reconstructed from the Marchenko kernel via

$$\begin{aligned} Q(x,t)=Q_+{\textbf{K}}(x,x;t)-\sigma _3 {\textbf{K}}(x,x;t)\sigma _3. \end{aligned}$$
(2.45b)

Solitons. As mentioned in the introduction, as a byproduct of the IST one can also obtain multisoliton solutions by setting \(\rho (z)\equiv 0\) for all \(z\in \Sigma\) and assigning an appropriate set of discrete eigenvalues and associated norming constants. With this ansatz, the inverse problem for the eigenfunctions reduces to an algebraic system, and simple, determinantal expressions for the multisoliton solutions can be obtained. For the defocusing NLS with NZBC, a one-soliton solution has the form (2.1), where \(\zeta _1=q_o e^{i\alpha }\), and the norming constant \(C_1\) is related to the soliton center \(x_o\) so that \(x_o q_o \sin \alpha =\log {\sqrt{|C_1|/(2q_o \sin \alpha )}}\). Dark solitons of the defocusing NLS have been experimentally observed early on as temporal pulses in optical fibers [73, 185], in thin magnetic films [55], in BECs [47], more recently in water waves [51], etc.

The focusing NLS, on the other hand, has a much richer soliton structure. Depending on the location of the discrete eigenvalue in the complex z-plane, one can have: (i) the Tajiri-Watanabe breather solution [168] for a discrete eigenvalue in generic position in \(D^+\); (ii) the Kuznetsov-Ma solution (which is periodic in t and homoclinic in z) when \(\zeta _1\) is located on the imaginary axis [124, 131]; (iii) the Akhmediev breather (periodic in x and homoclinic in t) when \(\zeta _1=q_oe^{i\alpha }\) (i.e., for \(\zeta _1\) on the circle \(C_o\)) [17]; (iv) finally, by taking the limit of a stationary (Kuznetsov-Ma) soliton solution as the discrete eigenvalue \(\zeta _1\) approaches the branch point \(iq_o\) with a suitable rescaling of the norming constant, one can obtain the celebrated Peregrine breather:

$$\begin{aligned} q_\P (x,t) = 1 - \frac{16 i t + 4}{4 x^2+16 t^2+1}. \end{aligned}$$
(2.46a)

It should be noted that some of these solutions have recently been observed experimentally in optical fibers [114, 115, 182] and in water waves [50], in spite of the difficulties introduced by the onset of modulational instability. We refer the reader to [35] for the explicit expressions and plots of the above solutions of the focusing NLS on a nontrivial background, and to [42, 186] for a generalization of the IST formalism to include discrete eigenvalues with higher multiplicity (see also [194] about the regularity of the multiple higher-order poles solitons of the focusing NLS equation with NZBC).

3 IST for Matrix NLS with NZBC

In reviewing the IST for multicomponent NLS systems, we will start with the square matrix NLS case because for this system the generalization of the IST is more or less straightforward, provided one suitably accounts for the multiplicity of the asymptotic eigenvalues of the Lax pair. The main complication is related to the fact that the norming constants are complex matrices, which requires novel tools and, at the same time, produces a richer soliton sector.

One of the main motivations for the renewed interest in matrix NLS systems is the increasing focus in recent years on the study of multicomponent BECs, and particularly spinor condensates. The equations describing spinor condensates exhibit nonlinear terms reflecting the SU(2) symmetry of the spin exchange interactions, and have no analog in nonlinear optics, where Kerr-type nonlinearities only depend on the squared modulii of the components. Spinor BECs formed by atoms with spin F are described by a macroscopic wave function with \(2F +1\) components, and give rise to various phenomena that are not present in single-component BECs, such as formation of spin domains, spin textures, and topological states. We refer the interested reader to the review articles [112, 113] and references therein for further details on the theory and experimental realizations of spinor BECs.

A number of theoretical works have been dealing with multicomponent vector solitons in \(F=1\) spinor BECs, which are described by a 3-component macroscopic wave function. In this review we will be concerned with the matrix NLS equation:

$$\begin{aligned} iQ_t+Q_{xx}-2\sigma QQ^\dagger Q=0,\qquad \sigma =\pm 1, \end{aligned}$$
(3.1)

where Q(xt) is a \(2\times 2\) matrix-valued function, and \(Q^\dagger\) is the Hermitian conjugate of Q. The system (3.1) with \(2\times 2\) symmetric matrix potential

$$\begin{aligned} Q(x,t)= \begin{pmatrix} q_1(x,t) &{} q_0(x,t) \\ q_0(x,t) &{} q_{-1}(x,t) \end{pmatrix} , \end{aligned}$$

has been proposed as a model to describe hyperfine spin \(F=1\) spinor Bose-Einstein condensates with either repulsive interatomic interactions and anti-ferromagnetic spin-exchange interactions (\(\sigma = +1\)), or attractive interatomic interactions and ferromagnetic spin-exchange interactions (\(\sigma =-1\)); \(q_1,q_0,q_{-1}\) are related to the vacuum expectation values of the three components of the quantum field operator in the three possible spin configurations \(1,0,-1\) [102, 103, 175]. Specifically, the matrix NLS (3.1), with a symmetric \(2\times 2\) matrix potential Q(xt) as above, is equivalent to the following three-component coupled system of equations:

$$\begin{aligned}&i\partial _t q_1+\partial ^2_{x}q_1-2\sigma q_1\left[ |q_1|^2+2 |q_o|^2 \right] -2\sigma q_o^2q_{-1}^*=0 \\&i\partial _t q_{-1}+\partial ^2_{x}q_{-1}-2\sigma q_{-1}\left[ |q_{-1}|^2+2|q_o|^2 \right] -2\sigma q_o^2q_{1}^*=0 \\&i\partial _t q_o+\partial ^2_{x}q_o-2\sigma q_o\left[ |q_1|^2+|q_o|^2 +|q_{-1}|^2\right] -2\sigma q_1q_o^*q_{-1}=0 \end{aligned}$$

for the wave function \(q_j(x,t)\) of atoms with magnetic spin quantum number \(j=0,\pm 1\). More generally, Q(xt) can be an \(m\times m\) symmetric matrix and describe higher order spin condensates. In the case of a rapidly decaying potential Q(xt), the above system has also been studied in [86], and its generalization to the \(F=2\) case has also been considered in [84, 85, 176]. The IST for the system (3.1) with NZBC was first presented in [104], following the approach proposed by Kawata and Inoue for the scalar defocusing NLS [110, 111], and some soliton solutions were obtained in [123, 175, 177]. The problem was recently revisited in [151] using a different approach, and some novel soliton solutions were presented.

In this section, we follow [151] and the notation introduced in Sect. 1 to review the IST with NZBC for the system (3.1). Like in the scalar case, the gauge transformation \(Q(x,t)={\tilde{Q}}(x,t)\,e^{-2i\sigma q_o^2t}\) allows one to reduce the system to

$$\begin{aligned} iQ_t+Q_{xx}+2\sigma \left( q_o^2I_m-QQ^\dagger \right) Q=0 \end{aligned}$$
(3.2)

where \(I_m\) denotes the \(m\times m\) identity matrix and where we have dropped the \(\tilde{}\) in the dependent variable for shortness. In this form, Eq. (3.2) admits time-independent NZBC:

$$\begin{aligned} Q(x,t)\rightarrow Q_\pm \qquad \text {as } \quad x\rightarrow \pm \infty , \end{aligned}$$
(3.3)

for which we assume the following constraints:

$$\begin{aligned} Q_{\pm }^\dagger Q_\pm =Q_\pm Q_\pm ^\dagger =q_o^2I_m, \end{aligned}$$
(3.4)

representing the analog of the symmetric boundary conditions considered in Sect. 1 for the scalar case.

It is worth mentioning that the matrix NLS equation remains integrable also in the form

$$\begin{aligned} iQ_t+Q_{xx}+2 Q\Sigma Q^\dagger \Omega =0 \end{aligned}$$
(3.5)

where \(\Sigma , \Omega\) are constant Hermitian matrices [172]. Without loss of generality, one can assume \(\Sigma , \Omega\) to be diagonal with entries \(0,\pm 1\), and therefore the only inequivalent novel reductions correspond to choosing \(\Sigma =\Omega =\sigma _3\), or \(\Sigma =-\Omega =\sigma _3\). The IST for these “mixed sign” matrix NLS equations was developed both in the case of decaying potentials and with NZBC, and soliton solutions were derived in various works [97, 142, 152]. In this review we will only be concerned with the focusing/defocusing square matrix NLS equations (3.1).

Lax pair, Riemann surface and uniformization coordinate. The Lax pair for the matrix NLS equation in the form (3.2) can be written as

$$\begin{aligned} \varphi _x = {\textbf{U}}\,\varphi \,, \qquad \varphi _t = {\textbf{V}}\,\varphi \,, \end{aligned}$$
(3.6a)

with

$$\begin{aligned} {\textbf{U}}(x,t,z)= & {} -ik\underline{\sigma }_3+\underline{Q},\qquad \\ {\textbf{V}}(x,t,z)= & {} -2ik^2\underline{\sigma }_3 +2k\underline{Q}+i\underline{\sigma }_3\left[ \underline{Q}_x +\sigma q_o^2I_{2m} -\sigma \underline{Q}\,\underline{Q}^\dagger \right] , \end{aligned}$$
(3.7a)
$$\begin{aligned} \underline{\sigma }_3= & {} \begin{pmatrix} I_m &{}0_m \\ 0_m &{}-I_m \end{pmatrix},\qquad \underline{Q}= \begin{pmatrix}0_m &{} Q\\ \sigma Q^\dagger &{}0_m\end{pmatrix}. \end{aligned}$$
(3.7b)

[Here and in the following \(0_m\) is used to denote the \(m\times m\) zero matrix.] Asymptotically as \(x\rightarrow \pm \infty\) the scattering problem becomes:

$$\begin{aligned} \varphi _x = U_\pm \,\varphi , \qquad U_\pm = -ik \underline{\sigma }_3 + \underline{Q}_\pm , \qquad \text {with } \quad \underline{Q}_\pm = \begin{pmatrix}0_m &{} Q_\pm \\ \sigma Q_\pm ^\dagger &{}0_m\end{pmatrix}, \end{aligned}$$
(3.8)

and the constraint (3.4) for the boundary conditions, which is equivalent to

$$\begin{aligned} \underline{Q}_\pm \underline{Q}_\pm ^\dagger =\underline{Q}_\pm ^\dagger \underline{Q}_\pm =q_o^{2} I_{2m}, \end{aligned}$$
(3.9)

implies that the eigenvalues of \(U_\pm\) are \(\pm i\sqrt{k^2-\sigma q_o^2}\), in this case each with multiplicity \(2\,m\). [Note that generalizing the IST for the matrix NLS system to NZBC that do not satisfy the constraint (3.9) is still an open problem.] To properly account for the branching of the eigenvalues, one can introduce the two-sheeted Riemann surface defined by \(\lambda ^2 = k^2 -\sigma q_o^2\) so that \(\lambda (k)\) is a single-valued function on this surface, and the same uniformization variable \(z=k+\lambda\) as in the scalar case, cf. (2.7), yielding the domains \(D^\pm\) in (2.8) where \(\mathop {\textrm{Im}}\nolimits \lambda\) is sign-definite (cf. Fig. 1). Like in the scalar case, \(\Sigma =(-\infty ,-q_o]\cup [q_o,\infty )\) in the defocusing case, and \(\Sigma ={\mathbb {R}}\cup [-iq_o,iq_o]\) in the focusing case, and \(\Sigma _o=\Sigma \setminus \left\{ \pm \sqrt{\sigma } q_o\right\}\).

Jost solutions and analyticity. The Jost solutions are defined in terms of the asymptotic eigenvectors of the equation defining the scattering problem, namely the first of (3.6). Taking into account (3.8), on either sheet of the Riemann surface one can write the asymptotic eigenvector matrix as

$$\begin{aligned} E_\pm (z) = I_{2m}-\frac{i}{z}\underline{\sigma }_3\underline{Q}_\pm , \end{aligned}$$
(3.10)

such that \(U_\pm E_\pm = -i\lambda \, E_\pm \,\underline{\sigma }_3\). Note that

$$\begin{aligned} \det E_\pm (z) = \left( \frac{2\lambda }{\lambda +k}\right) ^m = \gamma ^m(z), \qquad E^{-1}_\pm (z) = \frac{1}{\gamma (z)}\,[I_{2m}+\frac{i}{z}\underline{\sigma }_3 \underline{Q}_\pm ], \end{aligned}$$
(3.11a)

where \(\gamma (z)\) is as in (2.13), \(E_\pm ^{-1}\) are defined for all values of z such that \(\gamma (z)\ne 0\), i.e., away from the branch points: \(z\ne \pm iq_o\) is the focusing case (\(\sigma =-1\)), and \(z\ne \pm q_o\) in the defocusing case (\(\sigma =1\)).

As \(x\rightarrow \pm \infty\), the time evolution of the solutions of (3.6) satisfy \(\varphi _t = {\textbf{V}}_\pm \,\varphi\), with \({\textbf{V}}_\pm = -2ik^2 \underline{\sigma }_3 +2k \underline{Q}_\pm\), where we have used the boundary conditions (3.3), (3.9) and hence assumed \(\underline{Q}_x\rightarrow 0_{2\,m}\). Note that \({\textbf{V}}_\pm = 2k\,{\textbf{U}}_\pm\), hence \({\textbf{U}}_\pm\) and \({\textbf{V}}_\pm\) share the same eigenvectors. We can then define the Jost eigenfunctions \(\Phi (x,t,z)\) and \(\Psi (x,t,z)\) as the simultaneous solutions of both parts of the Lax pair (3.6) such that

$$\begin{aligned} \Phi (x,t,z)\equiv \,& {} (\phi (x,t,z)\,\, {\bar{\phi }}(x,t,z)) \\= \,& {} \left[ E_-(z)+o(1)\right] \textrm{e}^{-i\theta (x,t,z)\underline{\sigma }_3} \qquad \textrm{as}~x\rightarrow -\infty , \end{aligned}$$
(3.12a)
$$\begin{aligned} \Psi (x,t,z)\equiv \,& {} ({\bar{\psi }} (x,t,z)\,\, \psi (x,t,z)) \\= \,& {} \left[ E_+(z)+o(1)\right] \textrm{e}^{-i\theta (x,t,z)\underline{\sigma }_3} \qquad \textrm{as}~x\rightarrow +\infty , \end{aligned}$$
(3.12b)

where \(\theta (x,t,z)\) is given by (2.10) like in the scalar case, and \(E_\pm\) are given in (3.10). Modified eigenfunctions, with constant boundary conditions as \(x\rightarrow \pm \infty\), are defined like in the scalar case:

$$\begin{aligned} \left( M(x,t,z)\,\, {\bar{M}}(x,t,z)\right)= \,& {} \Phi (x,t,z)e ^{i\theta (x,t,z)\underline{\sigma }_3} \qquad \\ \left( {\bar{N}}(x,t,z)\,\, N(x,t,z)\right)= \,& {} \Psi (x,t,z)e ^{i\theta (x,t,z)\underline{\sigma }_3} . \end{aligned}$$
(3.13)

Then, the analog of Theorem 2.1 holds. In particular, assume the \(m\times m\) potential matrix Q(xt) in the matrix NLS equation is such that \(Q(x,t)-Q_+\in L^{1,1}([x_+,+\infty ))\) and \(Q(x,t)-Q_-\in L^{1,1}((-\infty ,x_-])\) as matrix functions of x for all \(t\geqslant 0\), for some \(x_\pm \in {\mathbb {R}}\), and that the boundary conditions \(Q_\pm\) satisfy the constraint (3.9). The modified eigenfunctions of the scattering problem (3.6) defined by (3.12) and (3.13) are such that M(xtz) and N(xtz) are analytic functions of z for \(z\in D^+\), and they are continuous up to \(\partial D^+\setminus \left\{ 0\right\}\); \({\bar{M}}(x,t,z)\) and \({\bar{N}}(x,t,z)\) are analytic functions of z for \(z\in D^-\), and they are continuous up to \(\partial D^-{\setminus } \left\{ 0\right\}\).

Scattering coefficients. Like in the scalar case, \(\Phi\) and \(\Psi\) are two fundamental matrix solutions of the scattering problem. Hence there exists a constant \(2m\times 2m\) scattering matrix S(z) such that

$$\begin{aligned} \Phi (x,t,z) = \Psi (x,t,z)\,S(z), \qquad S(z)= \begin{pmatrix} a(z) &{} {\bar{b}}(z) \\ b(z) &{} {\bar{a}}(z) \end{pmatrix}, \qquad x,t\in {\mathbb {R}},\qquad z\in \Sigma _o, \end{aligned}$$
(3.14)

where now \(a(z),b(z){\bar{a}}(z),{\bar{b}}(z)\) are \(m\times m\) blocks. Using (3.14) one can easily verify that:

$$\begin{aligned} \det a(z)= \,& {} \textrm{Wr}(\phi ,\psi )/\textrm{Wr}({\bar{\psi }}, \psi )\equiv \det (\phi \,\, \psi )/ \gamma ^m, \end{aligned}$$
(3.15a)
$$\begin{aligned} \det {\bar{a}}(z)= \,& {} \textrm{Wr}({\bar{\psi }},{\bar{\phi }})/\textrm{Wr}({\bar{\psi }}, \psi )\equiv \det ({\bar{\psi }}\,\, {\bar{\phi }})/ \gamma ^m, \end{aligned}$$
(3.15b)

where \(\textrm{Wr}(f, g)\) denotes the Wronskian determinant of the \(2\,m\times m\) matrices f and g, i.e., the determinant of the \(2m\times 2m\) matrix whose columns are given by the m columns of f and the m columns of g. While in the scalar case the analogs of such Wronskian representations allow one to establish the analyticity properties of the diagonal entries of the scattering matrix (cf. Sect. 2), in this case they only provide a proof of analyticity for the determinants of the diagonal blocks a(z) and \({\bar{a}}(z)\). However, one can follow the approach in [151] and obtain integral representations for the scattering matrix in terms of analytic eigenfunctions to prove that the upper diagonal block a(z) is analytic in \(D^+\) and continuous up to \(\Sigma _o\equiv \partial D^+{\setminus } \left\{ \pm \sqrt{\sigma }q_o\right\}\), and the lower diagonal block \({\bar{a}}(z)\) is analytic in \(D^-\) and continuous up to \(\Sigma _o\equiv \partial D_-{\setminus }\left\{ \pm \sqrt{\sigma }q_o\right\}\) (cf. (2.8)). The off-diagonal blocks of the scattering matrix S(z), namely, b(z) and \({\bar{b}}(z)\), are only defined for \(z\in \Sigma _o\), where they are continuous, but in general they do not admit analytic continuation off \(\Sigma _o\).

Finally, for \(z\in \Sigma _o\) we can use (3.13) and (3.14) to write

$$\begin{aligned} M(x,t,z)a^{-1}(z)= \,& {} {\bar{N}}(x,t,z)+e^{2i\theta (x,t,z)}N(x,t,z)\rho (z), \end{aligned}$$
(3.16a)
$$\begin{aligned} {\bar{M}}(x,t,z){\bar{a}}^{-1}(z)= \,& {} N(x,t,z)+e^{-2i\theta (x,t,z)}{\bar{N}}(x,t,z){\bar{\rho }}(z), \end{aligned}$$
(3.16b)

where \(M(x,t,z)a^{-1}(z)\) and \({\bar{M}}(x,t,z){\bar{a}}^{-1}(z)\) are meromorphic in \(D^+\) and \(D^-\), and we introduced (matrix) reflection coefficients

$$\begin{aligned} \rho (z)=b(z)a^{-1}(z), \qquad {\bar{\rho }}(z)={\bar{b}}(z){\bar{a}}^{-1}(z), \qquad z\in \Sigma _o. \end{aligned}$$
(3.17)

Symmetries. In analogy to the scalar case, we will need to consider the following involutions: 1. \(z\mapsto z^*\) (UHP/LHP), implying  \((k,\lambda )\mapsto (k^*,\lambda ^*)\) (same sheet); 2. \(z\mapsto \sigma q_o^2/z\) (outside/inside \(C_o\)), implying  \((k,\lambda )\mapsto (k,-\lambda )\) (opposite sheets). Both these transformations correspond to symmetries of the scattering problem; the first one is related to the usual conjugation symmetry in the potential, \(\underline{Q}^\dagger =\sigma \underline{Q}\), while the second one is a direct consequence of the branching in the plane of the scattering parameter k, or, equivalently, the fact that the uniformization variable z provides a double covering of the k plane. In addition, we will also have to consider a third symmetry that corresponds to assuming \(Q^T=Q\), which in terms of \(\underline{Q}\) can be written as

$$\begin{aligned} \underline{Q}=-\underline{\sigma }_2 \underline{Q}^T \underline{\sigma }_2 \qquad \text {where} \qquad \underline{\sigma }_2= \begin{pmatrix} 0_m &{} i I_m \\ -i I_m &{} 0_m \end{pmatrix}, \end{aligned}$$
(3.18)

\(\underline{\sigma }_2\) being a \(2m\times 2m\) generalization of the \(2\times 2\) Pauli matrix \(\sigma _2\). Determining the corresponding symmetries for the eigenfunctions in the matrix case is more complicated than in the scalar case, since for the first and third symmetry it involves considering appropriate bilinear combinations of eigenfunctions. We refer the reader to [151] for details, and we limit ourselves to summarizing below the results for the symmetries of the scattering data, since these are the only symmetries needed in the inverse problem.

Proposition 3.1

Let the \(m\times m\) potential matrix Q(xt) in the matrix NLS equation be such that \(Q(x,t)-Q_+\in L^1([x_+,+\infty ))\) and \(Q(x,t)-Q_-\in L^1((-\infty ,x_-])\) as matrix functions of x for all \(t\geqslant 0\), for some \(x_\pm \in {\mathbb {R}}\). For all \(z\in \Sigma _o\) the matrix reflection coefficients \(\rho (z)=b(z)a^{-1}(z)\) and \({\bar{\rho }}(z)={\bar{b}}(z){\bar{a}}^{-1}(z)\) defined in terms of the blocks of the scattering matrix S(z) in analogy to (2.17) satisfy the following symmetries:

$$\begin{aligned} \rho (z)=\sigma {\bar{\rho }}(z^*), \qquad \rho (\sigma q_o^2/z)=-\frac{\sigma }{q_o^2}Q_+^\dagger {\bar{\rho }}(z)Q_+^\dagger . \end{aligned}$$
(3.19)

Moreover, the diagonal blocks of the scattering matrix a(z) and \({\bar{a}}(z)\) satisfy the following symmetry properties for \(z\in D^-\cup \Sigma _o\):

$$\begin{aligned} \det {\bar{a}}(z)= \,& {} (\det a (z^*))^* \end{aligned}$$
(3.20a)
$$\begin{aligned} a(\sigma q_o^2/z)= \,& {} \frac{1}{q_o^2}Q_+{\bar{a}}(z)Q_-^\dagger \qquad \Rightarrow \qquad \\ \det {\bar{a}}(z)= \,& {} \frac{q_o^{2m}}{\det Q_+(\det Q_-)^*}\det a(\sigma q_o^2/z). \end{aligned}$$
(3.20b)

If, in addition, Q(xt) is a symmetric matrix (i.e., \(Q(x,t)=Q^T(x,t)\) for all \(x,t\in {\mathbb {R}}\)), then

$$\begin{aligned} \rho ^T(z)=\rho (z), \qquad {\bar{\rho }}^T(z)={\bar{\rho }}(z) \qquad z\in \Sigma _o \end{aligned}$$
(3.20c)

and

$$\begin{aligned} {\bar{a}}(z)=a^*(z^*) \qquad z\in D^-\cup \Sigma _o. \end{aligned}$$
(3.20d)

If \(Q(x,t)-Q_+\in L^{1,1}([x_+,+\infty ))\) and \(Q(x,t)-Q_-\in L^{1,1}((-\infty ,x_-])\) as matrix functions of x for all \(t\geqslant 0\), for some \(x_\pm \in {\mathbb {R}}\), all the above symmetries also extend to include the branch points, and hence are valid for \(z\in \Sigma\), and \(z\in D^-\cup \Sigma\), respectively.

Discrete eigenvalues and norming constants. The discrete spectrum of the scattering problem is the set of all values \(z\in {\mathbb {C}}\setminus \Sigma\) such that the scattering problem admits eigenfunctions in \(L^2({\mathbb {R}})\), and in this case they coincide with the zeros of \(\det a(z)\) in \(D^+\) and those of \(\det {\bar{a}}(z)\) in \(D^-\) where, according to (3.15a), the Jost eigenfunctions become linearly dependent. In the defocusing case, it is easy to show that for any \(z\in \Sigma\) and for any \(\xi \in {\mathbb {C}}^m\) the symmetries of the scattering data imply:

$$\begin{aligned} \Vert a(z)\xi \Vert ^2=\Vert \xi \Vert ^2+\Vert b(z)\xi \Vert ^2, \qquad \Vert {\bar{a}}(z)\xi \Vert ^2=\Vert \xi \Vert ^2+\Vert {\bar{b}}(z)\xi \Vert ^2. \end{aligned}$$

As a consequence, \(a(z)\xi =0\) implies \(\xi =0\), and therefore \(\det a(z)\ne 0\) for all \(z\in \Sigma\). The same of course holds for \({\bar{a}}(z)\). In the focusing case, however, one cannot exclude the existence of zeros of \(\det a(z)\) and \(\det {\bar{a}}(z)\) for \(z\in \Sigma\), i.e., of spectral singularities. No results are currently available in the literature regarding the location of eigenvalues and spectral singularities (or sufficient constraints on the potential for their absence) for the general square matrix ZS/AKNS scattering problem with NZBC, except for the result in [118] where it is shown that if the \(L^1\)-norm of the matrix potential is such that \(||Q||_1<\pi /2\) the scattering problem has no discrete eigenvalues or spectral singularities, and if \(||Q||_1=\pi /2\) there are no discrete eigenvalues, but there could be a spectral singularity. In fact, as discussed in Sect. 2, to the best of our knowledge the characterization of spectral singularities is an open problem even for the \(2\times 2\) ZS/AKNS scattering problem (corresponding to the scalar focusing NLS) with nonzero boundary conditions. In the following, when dealing with the focusing ZS/AKNS \(2m\times 2m\) scattering problem with nonzero boundary conditions, we will assume that there are no spectral singularities, namely that \(\det a(z)\ne 0\) for all \(z\in \Sigma\). This also implies \(\det {\bar{a}}(z)\ne 0\) for all \(z\in \Sigma\), because of (3.20a).

The symmetries of a(z) and \({\bar{a}}(z)\) imply that the discrete spectrum of the scattering problem (3.6) is the same as in the scalar case, given by:

$$\begin{aligned} \sigma= & {} -1 \quad \text {(focusing case):} \qquad {\mathcal {Z}} = \{z_j,-q_o^2/z_j^*,z_j^*,-q_o^2/z_j,\}_{j=1}^{\mathcal {N}} \\ \sigma= & {} 1 \quad \text {(defocusing case):} \qquad {\mathcal {Z}}=\{ \zeta _j, \zeta _j^*\}_{j=1}^{\mathcal {N}} \end{aligned}$$

where in the focusing case each first pair is in \(D^+\) (and we assume without loss of generality \(\mathop {\textrm{Im}}\nolimits z_n>0\)) and each second pair is in \(D^-\); in the defocusing case the eigenvalues are on the circle \(C_o\) (cf. Fig. 1).

Corresponding to each pair or quartet of simple discrete eigenvalues (namely, discrete eigenvalues that correspond to simple poles of \(a^{-1}(z)\) in \(D^+\) and \({\bar{a}}^{-1}(z)\) in \(D^-\)), one can introduce a pair or a quartet of \(m\times m\) norming constants \(\left\{ C_j,{\bar{C}}_j\right\}\) and \(\left\{ C_j,{\bar{C}}_j,{\hat{C}}_j,\hat{{\overline{C}}}_j\right\}\) such that

$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z=z_j}\big [ M(x,t,z)a^{-1}(z) \big ]= \,& {} \textrm{e}^{2i\theta (x,t,z_j)}N(x,t,z_j)\, C_j, \end{aligned}$$
(3.21a)
$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z=z_j^*}\big [{\bar{M}}(x,t,z){\bar{a}}^{-1}(z) \big ]= \,& {} \textrm{e}^{-2i\theta (x,t,z_j^*)}{\bar{N}}(x,t,z_j^*){\bar{C}}_j, \end{aligned}$$
(3.21b)
$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z={\hat{z}}_j\equiv \sigma q_o^2/z_j^*}\big [M(x,t,z)\, a^{-1}(z) \big ]= \,& {} \textrm{e}^{2i\theta (x,t,{\hat{z}}_j)}N(x,t,{\hat{z}}_j)\, {\hat{C}}_j, \qquad \\ {\hat{C}}_j= \,& {} \frac{1}{(z_j^*)^2}Q_+^\dagger {\bar{C}}_j Q_+^\dagger , \end{aligned}$$
(3.21c)
$$\begin{aligned} \mathop {\textrm{Res}}\limits _{z={\hat{z}}_j^*\equiv \sigma q_o^2/z_j}\big [ {\bar{M}}(x,t,z)\, {\bar{a}}^{-1}(z) \big ]= \,& {} \textrm{e}^{-2i\theta (x,t,{\hat{z}}_j^*)}{\bar{N}}(x,t,{\hat{z}}_j^*)\, \hat{{\overline{C}}}_j, \qquad \\ \hat{{\overline{C}}}_j= \,& {} \frac{1}{z_j^2}Q_+ C_j Q_+ . \end{aligned}$$
(3.21d)

Of course, in the defocusing case only the first two equations are considered.

Generically, when the discrete eigenvalues are simple zeros of \(\det a(z)\) and \(\det {\bar{a}}(z)\), the norming constants have zero determinant (and, in particular, in the case \(m=2\) this means that they are rank 1 matrices). On the other hand, it is possible to have discrete eigenvalues which are zeros of order \(p>1\) of \(\det a(z)\) and \(\det {\bar{a}}(z)\) but still poles of order 1 of \(a^{-1}(z)\) and \({\bar{a}}^{-1}(z)\). For instance, in the case \(m=2\), \(\det a(z)\) can have a zero of order 2 at a discrete eigenvalue \(z_j\), but \(a(z_j)=0_2\) and consequently \(Ma^{-1}\) still has a first order pole at \(z_j\) (the same of course holds at each of the symmetric points, because of the symmetries). In this case, the corresponding norming constants are full rank (i.e., their determinant is nonzero). This has no counterpart in the scalar NLS, and it has the consequence that the nature of the solitons depends on whether the associated norming constants are rank-one matrices (giving rise to “ferromagnetic” solitons), or full rank (corresponding to “polar” solitons).

Finally, one needs to determine the symmetries of the norming constants. The second involution requires \({\hat{C}}_j\) and \(\hat{{\overline{C}}}_j\) to be related to \(C_j\) and \({\bar{C}}_j\) as specified in (3.21). In [151] it is shown that, as a consequence of the first symmetry, the norming constants \(C_j,{\bar{C}}_j\) satisfy the same symmetry as in the matrix NLS with zero boundary conditions [13], namely:

$$\begin{aligned} {\bar{C}}_j = \sigma C_j^\dagger . \end{aligned}$$
(3.22a)

Moreover, the third symmetry also requires that \(C_j\) and \({\bar{C}}_j\) be symmetric matrices:

$$\begin{aligned} C_j=C_j^T, \qquad {\bar{C}}_j={\bar{C}}_j^T. \end{aligned}$$
(3.22b)

It is important to point out that the above symmetries have been obtained in [151] indirectly, using the reconstruction formula in the inverse problem to link the symmetries of the norming constants to the symmetries of the potential Q(xt), and in the case of decaying potentials the corresponding symmetry was obtained by assuming that the scattering relation (3.16) can be extended off the continuous spectrum and up to the location of the discrete eigenvalue, which is a fairly strong and not always justified assumption. To our knowledge, a rigorous derivation of all the symmetries of the norming constants within the framework of the direct scattering problem is still outstanding even in case of decaying potentials, and would represent a nice mathematical advance in the IST for matrix NLS equations.

Asymptotics of eigenfunctions and scattering data. The asymptotic properties of the eigenfunctions and the scattering matrix are needed to properly define the inverse problem. Moreover, the next-to-leading-order behavior of the eigenfunctions allows one to reconstruct the potential from the solution of the RHP. The asymptotic behavior of the modified Jost eigenfunctions for \(z\rightarrow \infty\) and \(z\rightarrow 0\) can be obtained via standard WKB expansions, or using integration by parts on the corresponding integral equations (cf. [151] for details). The asymptotic expansion as \(z\rightarrow \infty\) is given by

$$\begin{aligned} M(x,t,z)= & {} \begin{pmatrix} I_m+ \displaystyle {\frac{i\sigma }{z}} \mathop {\int }\limits _{-\infty }^x[Q(x^\prime ,t)Q^\dagger (x^\prime ,t) -q_o^2I_m]dx^\prime + O(1/z^2) \\ \displaystyle {\frac{i\sigma }{z}} Q^\dagger (x,t)+O(1/z^2) \end{pmatrix} \qquad \\{} & {} z\rightarrow \infty , z\in D^+, \end{aligned}$$
(3.23a)
$$\begin{aligned} {\bar{M}}(x,t,z)= & {} \begin{pmatrix} -\displaystyle {\frac{i}{z}}Q(x,t)+O(1/z^2) \\ I_m- \displaystyle {\frac{i\sigma }{z}} \mathop {\int }\limits _{-\infty }^x[Q^\dagger (x^\prime ,t)Q(x^\prime ,t) -q_o^2I_m]dx^\prime + O(1/z^2) \end{pmatrix} \qquad \\{} & {} z\rightarrow \infty , z\in D^-, \end{aligned}$$
(3.23b)
$$\begin{aligned} {\bar{N}}(x,t,z)= & {} \begin{pmatrix} I_m+ \displaystyle {\frac{i\sigma }{z}} \mathop {\int }\limits _x^{\infty }[Q(x^\prime ,t)Q^\dagger (x^\prime ,t) -q_o^2I_m]dx^\prime + O(1/z^2) \\ \displaystyle {\frac{i\sigma }{z}} Q^\dagger (x,t)+O(1/z^2) \end{pmatrix} \qquad \\{} & {} z\rightarrow \infty , z\in D^-, \end{aligned}$$
(3.23c)
$$\begin{aligned} N(x,t,z)= & {} \begin{pmatrix} -\displaystyle {\frac{i}{z}} Q(x,t) +O(1/z^2) \\ I_m- \displaystyle {\frac{i\sigma }{z}} \mathop {\int }\limits _{-\infty }^x[Q^\dagger (x^\prime ,t)Q(x^\prime ,t) -q_o^2I_m]dx^\prime + O(1/z^2) \end{pmatrix} \qquad \\{} & {} z\rightarrow \infty , z\in D^+. \end{aligned}$$
(3.23d)

Similarly, asymptotics as \(z\rightarrow 0\) in the proper region \(D^\pm\) yields:

$$\begin{aligned} M(x,t,z)= & {} \begin{pmatrix} QQ_-^\dagger /q_o^2+O(z) \\ i\sigma Q_-^\dagger /z +O(1) \end{pmatrix}, \qquad {\bar{M}}(x,t,z) \begin{pmatrix} -iQ_-/z+O(1) \\ Q^\dagger Q_-/q_o^2+O(z) \end{pmatrix}, \end{aligned}$$
(3.24a)
$$\begin{aligned} {\bar{N}}(x,t,z)= & {} \begin{pmatrix} QQ_+^\dagger /q_o^2+O(z) \\ i\sigma Q_+^\dagger /z +O(1) \end{pmatrix}, \qquad N(x,t,z) \begin{pmatrix} -iQ_+/z+O(1) \\ Q^\dagger Q_+/q_o^2+O(z) \end{pmatrix}. \end{aligned}$$
(3.24b)

Finally, inserting the above asymptotic expansions for the Jost eigenfunctions into (3.14) one shows that, as \(z\rightarrow \infty\) in the appropriate regions of the complex z-plane,

$$\begin{aligned} S(z)= I_{2m} + O(1/z). \end{aligned}$$
(3.25)

The above asymptotics holds with \(\mathop {\textrm{Im}}\nolimits z\geqslant 0\) and \(\mathop {\textrm{Im}}\nolimits z\leqslant 0\) for a(z) and \({\bar{a}}(z)\), respectively, and with \(z\in \Sigma\) for b(z) and \({\bar{b}}(z)\). Similarly, one can show that as \(z\rightarrow 0\)

$$\begin{aligned} S(z)= \frac{1}{q_o^2}\begin{pmatrix} Q_+Q_-^\dagger &{} 0_m \\ 0_m &{} Q_+^\dagger Q_-\end{pmatrix} + O(z), \end{aligned}$$
(3.26)

where again the asymptotics for the block diagonal entries of S(z) can be extended to \(D^+\) for a(z) and to \(D^-\) for \({\bar{a}}(z)\), while the asymptotics for the off-diagonal blocks hold for \(z\in \Sigma\).

Inverse problem. The starting point for the formulation of the inverse problem for the eigenfunctions is (3.16), which we now regard as a relation between eigenfunctions analytic in \(D^+\) and those analytic in \(D^-\). We introduce the sectionally meromorphic matrices

$$\begin{aligned} \mu ^+(x,t,z) = (M\,a^{-1}\,\, \,N),\qquad \mu ^-(x,t,z) = ({\bar{N}}\,\,\, {\bar{M}}\,{\bar{a}}^{-1}), \end{aligned}$$
(3.27)

where superscripts ± distinguish between analyticity in \(D^+\) and \(D^-\), respectively. From (3.16) we then obtain the jump condition

$$\begin{aligned} \mu ^-(x,t,z)= \mu ^+(x,t,z)\,(I_{2m}-G(x,t,z)) \qquad z\in \Sigma , \end{aligned}$$
(3.28)

where the jump matrix is

$$\begin{aligned} G(x,t,z) = \begin{pmatrix} 0_m &{} -\textrm{e}^{-2i\theta (x,t,z)} {\bar{\rho }}(z) \\ \textrm{e}^{2i\theta (x,t,z)}\rho (z) &{} \rho (z){\bar{\rho }}(z) \end{pmatrix}. \end{aligned}$$
(3.29)

Recalling the asymptotic behavior of the Jost eigenfunctions and scattering coefficients, it is easy to check that

$$\begin{aligned} \mu ^\pm= \,& {} I_{2m} + O(1/z),\qquad z\rightarrow \infty , \end{aligned}$$
(3.30a)
$$\begin{aligned} \mu ^\pm= & {} -(i/z)\underline{\sigma }_3\,\underline{Q}_+ + O(1),\qquad z\rightarrow 0. \end{aligned}$$
(3.30b)

Subtracting the asymptotic behaviors as \(z\rightarrow \infty\) and \(z\rightarrow 0\), as well as the pole contributions, one can write a formal solution to the above RHP as:

$$\begin{aligned} \mu (x,t,z)= \,& {} I_{2m} - (i/z)\underline{\sigma }_3\,\underline{Q}_+ + \sum _{j=1}^{J}\frac{\mathop {\textrm{Res}}\limits _{\zeta _j}\mu ^+}{z-\zeta _j} + \sum _{j=1}^{J}\frac{\mathop {\textrm{Res}}\limits _{\zeta _j^*}\mu ^-}{z-\zeta _j^*} \\{} & {} + \frac{1}{2\pi i}\mathop {\int }\limits _\Sigma \frac{\mu ^+(x,t,\zeta )}{\zeta -z}\,G(x,t,\zeta )\,\textrm{d}\zeta , \qquad z\in {\mathbb {C}}\setminus \Sigma , \end{aligned}$$
(3.31)

where \(J={\mathcal {N}}\) in the defocusing case, while in the focusing case \(J=2{\mathcal {N}}\) and \(\zeta _j= z_j\) and \(\zeta _{j+{\mathcal {N}}} = \sigma q_o^2/z_j^*\) for \(j=1,\dots ,{\mathcal {N}}\). Explicitly, taking into account the residue conditions (3.21), one can obtain the following expressions for the analytic columns of \(\mu ^\pm\):

$$\begin{aligned} N(x,t,z)= & {} \begin{pmatrix} - iQ_+/z \\ I_m \end{pmatrix} + \sum _{j=1}^{J}\frac{\textrm{e}^{-2i\theta (x,t,\zeta _j^*)}}{z-\zeta _j^*}\,{\bar{N}}(x,t,\zeta _j^*)\, {\bar{C}}_j \\{} & {} - \frac{1}{2\pi i}\mathop {\int }\limits _\Sigma \frac{e^{-2i\theta (x,t,\zeta )}{\bar{N}}(x,t,\zeta ){\bar{\rho }}(\zeta )}{\zeta -(z+i0)}\,\textrm{d}\zeta , \end{aligned}$$
(3.32a)
$$\begin{aligned} {\bar{N}}(x,t,z)= & {} \begin{pmatrix} I_m \\ i\sigma Q_+^\dagger /z \end{pmatrix} + \sum _{J}^{2{\mathcal {N}}}\frac{\textrm{e}^{2i\theta (x,t,\zeta _j)}}{z-\zeta _j}\,N(x,t,\zeta _j)\, C_j \\{} & {} + \frac{1}{2\pi i}\mathop {\int }\limits _\Sigma \frac{e^{2i\theta (x,t,\zeta )}N(x,t,\zeta )\rho (\zeta )}{\zeta -(z-i0)}\,\textrm{d}\zeta , \end{aligned}$$
(3.32b)

where in the focusing case \(C_{j+{\mathcal {N}}}={\hat{C}}_j\) and \({\bar{C}}_{j+{\mathcal {N}}}=\hat{{\bar{C}}}_{h}\) for \(j=1,\dots ,{\mathcal {N}}\) (cf. (3.21)).

Finally, the last step is to reconstruct the potential from the solution of the RHP. From (3.31), one obtains the asymptotic behavior of \(\mu ^\pm (x,t,z)\) as \(z\rightarrow \infty\), and comparing it with (3.23d) we then find:

$$\begin{aligned} Q(x,t)= \,& {} Q_+ + i \sum _{j=1}^{J}\textrm{e}^{-2i\theta (x,t,\zeta _j^*)}\,{\bar{N}}^{\textrm{up}}(x,t,\zeta _j^*)\, {\bar{C}}_j \\{} & {} + \frac{1}{2\pi }\mathop {\int }\limits _\Sigma e^{-2i\theta (x,t,\zeta )}{\bar{N}}^{\textrm{up}}(x,t,\zeta ){\bar{\rho }}(\zeta )\,\textrm{d}\zeta , \end{aligned}$$
(3.33)

where the superscript \(^{\textrm{up}}\) denotes the upper \(m\times m\) block of the modified eigenfunction \({\bar{N}}\).

We mention that the analog of the trace formula (2.42) for \(\det a(z)\) can be obtained, and this in turn can provide a relationship between \(\det Q_+/Q_-\) which plays the role of a \(\theta\)-condition. On the other hand, finding an explicit expression for the analytic matrix a(z) in terms of a minimal set of scattering data (i.e., reflection coefficients, discrete eigenvalues and norming constants) is an open problem even in the case of decaying potentials, for which a closed-form solution is only known in the pure soliton case (see [52] for a derivation in the context of coupled Maxwell–Bloch equations).

Finally, we point out that the inverse problem can also be formulated in terms of Marchenko integral equations, and the analog of Eqs. (2.44) and (2.45) can be obtained.

Solitons. The reconstruction formula (3.33) can be used to obtain soliton solutions of the matrix NLS with NZBC. Some solutions are matrix generalizations of the scalar solitons of the focusing and defocusing NLS described in Sect. 2, but in the matrix case each solution has a polar and a ferromagnetic counterpart, depending on whether the associated norming constant has zero determinant or not. Physically, the signature of ferromagnetic solitons is that they have nonzero total spin, and exhibit a kink-like behavior in the form of domain walls between the \(\pm 1\) components and the 0-component. Importantly, while in the decaying case all matrix one-soliton solutions are reducible via unitary transformations to a combination of oppositely-polarized scalar solitons, when a nonzero background is present, not all matrix one-soliton solutions are reducible to a simple combination of scalar solutions (see [129] for the focusing case, and [3] for the defocusing case). Finally, in the focusing case, by taking suitable limits of all the solutions on a nonzero background as the discrete eigenvalue approaches the branch point, one obtains families of rogue wave (i.e., rational) solutions, which generalize the Peregrine solution (2.46). Some examples of soliton solutions for the focusing and defocusing matrix NLS (3.2) are shown in Figs. 2, 3, 4, 5.

Fig. 2
figure 2

Ferromagnetic 1-soliton solution of the focusing matrix NLS (\(q_1\), \(q_{-1}\) and \(q_0\) from left to right) with \(Q_+=I_2\), \(\zeta _1=1+2i\) and the entries of \(C_1\) are chosen to be \(C_{11}=C_{22}=C_{12}=1+i\) (\(\det C_1=0\)). From [151]

Fig. 3
figure 3

Polar 1-soliton solution of the focusing matrix NLS (\(q_1\), \(q_{-1}\) and \(q_0\) from left to right) with \(Q_+=I_2\), \(\zeta _1=2i\) and the entries of \(C_1\) chosen to be \(C_{11}=C_{22}=0\) and \(C_{12}=C_{21}=1\) (\(\det C_1\ne 0\)). Note the solution is homoclinic in x and periodic in t, like its scalar counterpart, the Kuznetsov-Ma breather of the focusing NLS equation. From [151]

We refer the reader to [3, 123, 129, 151, 175, 177] for the explicit expressions and additional plots of the soliton solutions of matrix NLS (3.2) with NZBC. To the best of our knowledge, a generalization of the IST formalism to include discrete eigenvalues with higher multiplicity has not been developed for the square matrix NLS with NZBC.

4 IST for Vector NLS with NZBC

In many applications, vector NLS (VNLS) systems are the key governing equations. Physically, the VNLS equation arises under conditions similar to those described by NLS whenever there are multiple wavetrains moving with nearly the same group velocity. Moreover, VNLS also models systems where the electromagnetic field has more than one component, e.g., in optical fibers and waveguides, where the propagating electric field has two polarized components transverse to the direction of propagation. The dimensionless system

$$\begin{aligned} i{\textbf{q}}_t+{\textbf{q}}_{xx}-2\sigma \Vert {\textbf{q}}\Vert ^2{\textbf{q}}=0 \qquad \sigma =\pm 1, \end{aligned}$$
(4.1)

where \({\textbf{q}}(x,t)\) is a two-component, complex vector function, was originally introduced by Manakov [133] for \(\sigma =-1\) as an asymptotic model governing the propagation of the electric field in a waveguide. Subsequently, Menyuk showed in [134] that in optical fibers with constant birefringence, assuming certain nonlinear (four-wave mixing) terms are neglected, the two polarization components of the complex electromagnetic field envelope orthogonal to direction of propagation along a fiber satisfy asymptotically the following nondimensional equations:

$$\begin{aligned} i(u_t+\delta u_x)+ d u_{xx}+(|u|^2+\alpha |v|^2)= \,& {} 0, \\ i(v_t+\delta v_x)+ d v_{xx}+(\alpha |u|^2+|v|^2)= \,& {} 0, \end{aligned}$$

where \(\delta\) represents the group velocity “mismatch” between the components, 2d is the group velocity dispersion, and \(\alpha\) is a constant depending on the polarization properties of the fiber. The physical phenomenon of birefringence implies that the phase and group velocities of the electromagnetic wave are different for each polarization component. When \(\alpha \ne 1\) the above system is not integrable. However, averaging over the fast birefringence fluctuations that are normally observed in a communications environment [135] yields \(\delta =0\) and \(\alpha =1\)—that is, the system reduces to the VNLS (4.1), which therefore attains broader relevance.

The focusing (\(\sigma =-1\)) VNLS admits soliton solutions which are the analog of the bright soliton (1.2), multiplied by a constant, norm-1 polarization vector \({\textbf{p}}\). Manakov was able to integrate the VNLS system (4.1) by the IST method in the case of a rapidly decaying potential. He also showed that, unlike scalar soliton, vector solitons generically interact nontrivially: while the soliton velocities are preserved and the interaction is still elastic in the sense that the total energy is conserved, there can be a redistribution of energy between the components expressed by a change in the asymptotic form of the polarization vectors—a “polarization shift”.

Fig. 4
figure 4

Rogue wave solution of the focusing matrix NLS (\(q_1\), \(q_{-1}\) and \(q_0\) from left to right). \(Q_+=I_2\), \(\zeta _1=i\) and the norming constant is chosen as in Fig. 10 of [129]

Fig. 5
figure 5

1-soliton profiles for ferromagnetic states for the defocusing matrix NLS with \(Q_+=I_2\), generated by a discrete eigenvalue \(\zeta _1 = e^{0.927i}\) and a norming constant \(C_1\) with, from left to right: \(|C_{11}|/|C_{22}| = 4\),  \(|C_{11}|/|C_{22}| = 1\) and \(|C_{11}|/|C_{22}| = 4/9\). In each plot, the three components \(q_1\) (black solid line), \(q_o\) (red dot-dashed line) and \(q_{-1}\) (blue dotted line) are shown. From [3]

The square matrix NLS discussed in Sect. 3 can be considered more generally for an \(n\times m\) potential matrix Q(xt) (obviously, without the requirement that it be symmetric), and the latter can still be expressed as the compatibility condition of the Lax pair (3.6), with \(Q,Q^\dagger\) rectangular matrices of appropriate sizes. Then one can obtain a vector reduction by choosing, for instance, \(n=1\) and arbitrary m or vice-versa, and the square matrix case corresponds, of course, to \(m=n\). In the case of rapidly decaying potentials, this framework has been exploited to generalize the IST and Manakov’s results for soliton interactions to a VNLS with any number of components (see [13] for details). When NZBC are considered, however, the vector case (2 or, more generally, m-component VNLS) and the square matrix case are completely different as far as the development of the IST is concerned. While in the scalar and square matrix NLS the scattering (Jost) eigenfunctions yield a complete set of analytic eigenfunctions, this is not the case for the VNLS, and therefore one needs to develop a procedure to supplement the missing analytic eigenfunctions. This was done for the defocusing and focusing Manakov system with NZBC in [121, 149] respectively, by generalizing an idea originally introduced by Kaup in the context of the IST for the three-wave interaction equation [109], which used cross products of eigenfunctions of the adjoint scattering problem to obtain additional analytic eigenfunctions. While in the focusing case this procedure can be generalized in a more or less straightforward way to arbitrary m, in the defocusing case \(m=2\) is significantly different from \(m>2\), and a completely different approach needed to be developed, which will be discussed in Sect. 4.2.

4.1 Manakov System

Here we review the essential elements of the IST for the defocusing and focusing Manakov system with NZBC. We omit the calculations, referring the reader to [1, 36, 121, 149] for further details.

Lax pair and uniformization coordinate. As usual, we rewrite the Manakov system (4.1) in the form

$$\begin{aligned} i{\textbf{q}}_t+{\textbf{q}}_{xx}-2\sigma (\Vert {\textbf{q}}\Vert ^2-q_o^2)\,{\textbf{q}}=0 \qquad \sigma =\pm 1, \end{aligned}$$
(4.2)

with \({\textbf{q}}=(q_1(x,t),q_2(x,t))^T\) satisfying the following (constant) NZBC at infinity:

$$\begin{aligned} \lim _{x\rightarrow \pm \infty }{\textbf{q}}(x,t) ={\textbf{q}}_{\pm } = {\textbf{q}}_o e^{\pm i\theta }, \end{aligned}$$
(4.3)

where \(q_o = \Vert {\textbf{q}}_o\Vert > 0\) is the amplitude of the (symmetric) nontrivial background. Note that (4.3) restricts the boundary conditions \({\textbf{q}}_\pm\) to be parallel to each other. When \({\textbf{q}}_+\) and \({\textbf{q}}_-\) are not parallel, a different approach to the IST is required (see [2]), but it will not be considered here. Moreover, the asymptotic phases were chosen as \(\theta _+=-\theta _-=\theta\) without loss of generality thanks to the phase invariance of Eq. (4.2).

The Lax pair for the system (4.2) is given by

$$\begin{aligned} \varphi _x = {\textbf{U}} \,\varphi ,\qquad \varphi _t = {\textbf{V}}\,\varphi , \end{aligned}$$
(4.4a)

where

$$\begin{aligned} {\textbf{U}}(x,t,k)= & {} -ik J + Q,\qquad \\ {\textbf{V}}(x,t,k)= & {} -2ik^{2} J -2k Q+ i J (Q_x - \sigma Q^2 +\sigma q_o^2 ), \end{aligned}$$
(4.4b)
$$\begin{aligned} J= \,& {} \textrm{diag}(1,-1,-1),\qquad Q(x,t) = \begin{pmatrix} 0&{} {\textbf{q}}^T \\ \sigma {\textbf{q}}^* &{} 0_{2\times 2} \end{pmatrix}. \end{aligned}$$
(4.4c)

Like in the previous sections, in order to remove the branching of the eigenfunctions in the spectral parameter k due to their dependence on \(\lambda\) given by \(\lambda ^2=k^2-\sigma q_o^2\), we introduce the uniformization variable z defined via the conformal mapping (2.7), and identify the domains \(D^\pm\) in (2.8) where \(\mathop {\textrm{Im}}\nolimits \lambda\) is sign definite (cf. Fig. 1).

Jost eigenfunctions and scattering data. As usual, the Jost eigenfunctions are defined as

$$\begin{aligned} \Phi (x,t,z)\equiv\, & {} [\phi _1(x,t,z),\phi _2(x,t,z),\phi _3(x,t,z)] \\=\, & {} \left[ E_-(z) + o(1)\right] e^{i\Theta (x,t,z)}, \qquad x\rightarrow - \infty , \end{aligned}$$
(4.5)
$$\begin{aligned} \Psi (x,t,z)\equiv & {} [\psi _1(x,t,z),\psi _2(x,t,z),\psi _3(x,t,z)] \\= & {} \left[ E_+(z) + o(1)\right] e^{i\Theta (x,t,z)}, \qquad x\rightarrow + \infty , \end{aligned}$$
(4.6)

where

$$\begin{aligned} \Theta= \,& {} {\varvec{\Lambda }} x-{\varvec{\Omega }} t = \textrm{diag}(\theta _1, \theta _2,-\theta _1), \end{aligned}$$
(4.7a)
$$\begin{aligned} {\varvec{\Lambda }}= & {} \textrm{diag}(-\lambda ,k,\lambda ), \qquad {\varvec{\Omega }}=\textrm{diag}(2k\lambda ,-(k^2+\lambda ^2),-2k\lambda ), \end{aligned}$$
(4.7b)

and \(k,\lambda\) are thought of as functions of the uniformization variable z via (2.7). In (4.5), \(E_\pm (z)\) are the matrices of asymptotic eigenvectors of the Lax pair, which can be chosen as follows:

$$\begin{aligned} E_\pm (z)=\begin{pmatrix} z &{} 0 &{} -\sigma q_o^2/z \\ i\sigma {\textbf{q}}_\pm ^* &{} -i{\textbf{q}}_\pm ^\perp &{} -i\sigma {\textbf{q}}_\pm ^* \end{pmatrix}, \end{aligned}$$
(4.8)

where for any two-component vector \({\textbf{p}}=(p_1,p_2)^T\) we take \({\textbf{p}}^\perp =(p_2,-p_1)^T\).

If \({\textbf{q}}(x,t)\) approaches \({\textbf{q}}_\pm\) sufficiently rapidly as \(x\rightarrow \pm \infty\) for all t, the Jost eigenfunctions are well-defined and continuous for all \(x,t,z\in \Sigma {\setminus }\left\{ 0\right\}\), where, as before, \(\Sigma ={\mathbb {R}}\) in the defocusing case, and \(\Sigma ={\mathbb {R}}\cup C_o\) in the focusing case. Furthermore, since \(\det E_\pm (z)=-2\sigma q_o^2\lambda (z)\ne 0\) for any \(z\in \Sigma _o\equiv \Sigma {\setminus } \{\pm \sqrt{\sigma } q_o\}\), \(\Phi (x,t,z)\) and \(\Psi (x,t,z)\) are two fundamental matrix solutions of the Lax pair (4.4a), and hence there exists a \(3 \times 3\) invertible matrix A(z) such that

$$\begin{aligned} \Phi (x,t,z) = \Psi (x,t,z) A(z),\qquad z \in \Sigma _o. \end{aligned}$$
(4.9)

As usual, \(A(z) = (a_{ij}(z))\) is referred to as the scattering matrix. For future convenience, we also introduce the inverse matrix \(B(z):= A^{-1}(z) = (b_{ij}(z))\). As before, the eigenfunctions have a pole at \(z=0\), and the behavior of eigenfunctions and scattering coefficients as \(z\rightarrow 0\) and as \(z\rightarrow \infty\) is analyzed separately.

When it comes to the analyticity of the Jost eigenfunctions, the vector problem is different than the scalar and square matrix ones discussed previously, and the difference can be traced to the presence of the additional eigenvalue ik in the \(3\times 3\) scattering problem (cf. (4.7)), which implies that the columns \(\phi _2,\psi _2\) of the Jost eigenfunctions cannot be analytically continued in either \(D^+\) or \(D^-\) where \(\mathop {\textrm{Im}}\nolimits \lambda\) is sign definite, since their analyticity depends on the sign of \(\mathop {\textrm{Im}}\nolimits (\lambda \pm k)\). Specifically, for the Manakov system the following theorem holds (cf. [149] and [121]).

Theorem 4.1

Let \({\textbf{q}}(x,t)-{\textbf{q}}_+\in L^1([x_+,+\infty ))\) and \({\textbf{q}}(x,t)-{\textbf{q}}_-\in L^1((-\infty ,x_-])\) as functions of x for all \(t\geqslant 0\), for some \(x_\pm \in {\mathbb {R}}\), and \(D^\pm\) are defined in (2.8) for the defocusing/focusing cases (\(\sigma =\pm 1\), respectively).

For the defocusing Manakov system, \(\phi _1(x,t,z)\) and \(\psi _3(x,t,z)\) are analytic functions of z for \(z\in D^+\), and they are continuous up to \(\partial D^+ {\setminus }\left\{ \pm \sqrt{\sigma } q_o,0\right\}\); \(\phi _3(x,t,z)\) and \(\psi _1(x,t,z)\) are analytic functions of z for \(z\in D^-\), and they are continuous up to \(\partial D^- {\setminus }\left\{ \pm \sqrt{\sigma } q_o,0\right\}\). On the other hand, \(\phi _2(x,t,z)\) and \(\psi _2(x,t,z)\) are continuous for \(z\in \Sigma _o\) but cannot be extended off it.

For the focusing Manakov system, the following columns of \(\Phi (x,t,z)\) or, correspondingly, \(\Psi (x,t,z)\) can be analytically extended onto the corresponding regions of the complex z-plane:

$$\begin{aligned}{} & {} \phi _1: \, \, \, z \in D_{1}, \qquad \phi _2: \, \, \, \mathop {\textrm{Im}}\nolimits z<0, \qquad \phi _3: \, \, \, z \in D_{4}, \end{aligned}$$
(4.10a)
$$\begin{aligned}{} & {} \psi _1: \, \, \, z \in D_2, \qquad \psi _2: \, \, \, \mathop {\textrm{Im}}\nolimits z >0, \qquad \psi _3: \, \, \, z \in D_{3}, \end{aligned}$$
(4.10b)

where the domains of analyticity \(D_1,\dots ,D_4\) are

$$\begin{aligned} D_{1}= & {} \{z: \mathop {\textrm{Im}}\nolimits z>0 \wedge |z|>q_o \}, \qquad D_2 = \{z: \mathop {\textrm{Im}}\nolimits z<0 \wedge |z|>q_o \}, \end{aligned}$$
(4.11a)
$$\begin{aligned} D_{3}= & {} \{z: \mathop {\textrm{Im}}\nolimits z<0 \wedge |z|<q_o \}, \qquad D_{4} = \{z: \mathop {\textrm{Im}}\nolimits z>0 \wedge |z|<q_o \}. \end{aligned}$$
(4.11b)

Furthermore, if \({\textbf{q}}(x,t)-{\textbf{q}}_+\in L^{1,1}([x_+,+\infty ))\) and \({\textbf{q}}(x,t)-{\textbf{q}}_-\in L^{1,1}((-\infty ,x_-])\) as functions of x for all \(t\geqslant 0\), for some \(x_\pm \in {\mathbb {R}}\), then the eigenfunctions \(\phi _j,\psi _j\) for \(j=1,3\) are also continuous up to the branch points \(\pm \sqrt{\sigma } q_o\).

In both cases, there is no complete set of analytic eigenfunctions in any given domain, and therefore one needs to introduce auxiliary eigenfunctions to obtain a complete set. As mentioned at the beginning of Sect. 4, for the \(3\times 3\) scattering problem considered here this can be achieved by generalizing the approach Kaup introduced for the three-wave interaction equation [109], based on the observation that the cross products of adjoint scattering eigenfunctions can be used to define eigenfunctions of the scattering problem. Specifically, one considers the formal “adjoint” of the Lax pair (4.4a):

$$\begin{aligned} {\tilde{\varphi }}_x = {\textbf{U}}^*{\tilde{\varphi }}, \qquad {\tilde{\varphi }}_t = {\textbf{V}}^*{\tilde{\varphi }}, \end{aligned}$$
(4.12)

and defines the adjoint Jost solutions as the simultaneous solutions \({\tilde{\Phi }}(x,t,z)\), \({\tilde{\Psi }}(x,t,z)\) of (4.12) such that

$$\begin{aligned} {\tilde{\Phi }}(x,t,z)\equiv\, & {} [{\tilde{\phi }}_1, {\tilde{\phi }}_2, {\tilde{\phi }}_3]= [{\tilde{E}}_-(z)+o(1)]e^{-i\Theta (x,t,z)}, \qquad x\rightarrow - \infty , \end{aligned}$$
(4.13)
$$\begin{aligned} {\tilde{\Psi }}(x,t,z)\equiv\, & {} [{\tilde{\psi }}_1, {\tilde{\psi }}_2, {\tilde{\psi }}_3]= [{\tilde{E}}_+(z)+o(1)]e^{-i\Theta (x,t,z)}, \qquad x\rightarrow + \infty , \end{aligned}$$
(4.14)

with \({\tilde{E}}_\pm (z) = E_\pm ^*(z^*)\) for \(z\in \Sigma\). Under the same assumptions as in Theorem 4.1: \({\tilde{\phi }}_1\) and \({\tilde{\psi }}_3\) are analytic in \(D^-\), \({\tilde{\phi }}_3\) and \({\tilde{\psi }}_3\) are analytic in \(D^+\), and \({\tilde{\phi }}_2\) and \({\tilde{\psi }}_2\) do not admit analytic continuation off \(\Sigma\) in the defocusing case, and they are analytic for \(\mathop {\textrm{Im}}\nolimits z>0\) and \(\mathop {\textrm{Im}}\nolimits z<0\), respectively, in the focusing case.

How to proceed to define the auxiliary analytic eigenfunctions (and, consequently, how to set up the inverse problem) depends on whether the focusing or defocusing cases are considered, and a fully unified treatment along the lines of Sects. 2 and 3 is not as feasible. Therefore, from this point on, we will discuss the two problems separately.

4.1.1 Defocusing Manakov System

Auxiliary eigenfunctions. In the defocusing case, two additional solutions of the original Lax pair (4.4a) can then be obtained from the adjoint Jost solutions (4.12) as:

$$\begin{aligned} \chi (x,t,z)= & {} - e^{i \theta _2(x,t,z)} J [{\tilde{\phi }}_{3}(x,t,z) \times {\tilde{\psi }}_{1}(x,t,z)], \end{aligned}$$
(4.15a)
$$\begin{aligned} {\bar{\chi }}(x,t,z)= & {} - e^{i \theta _2(x,t,z)} J [{\tilde{\phi }}_{1}(x,t,z) \times {\tilde{\psi }}_{3}(x,t,z)]. \end{aligned}$$
(4.15b)

By construction, \(\chi (x,t,z)\) is analytic for \(z\in D^+\), while \({\bar{\chi }}(x,t,z)\) is analytic for \(z\in D^-\). Moreover, one can show that, for \(z \in \Sigma _o\):

$$\begin{aligned} \chi (z)= \,& {} 2\lambda (z)[b_{33}(z) \psi _{2}(z) - b_{23}(z) \psi _{3}(z)] \\= \,& {} 2\lambda (z)[a_{11}(z) \phi _{2}(z) - a_{21}(z) \phi _{1}(z)], \end{aligned}$$
(4.16a)
$$\begin{aligned} {\bar{\chi }}(z)= \,& {} 2\lambda (z)[b_{21}(z) \psi _{1}(z) - b_{11}(z) \psi _{2}(z)] \\= \,& {} 2\lambda (z)[a_{23}(z) \phi _{3}(z) - a_{33}(z) \phi _{2}(z)], \end{aligned}$$
(4.16b)

where the (xt)-dependence in the eigenfunctions is omitted for simplicity. One can now define two complete sets of eigenfunctions for the system (4.4a) as:

$$\begin{aligned} \Phi ^+(x,t,z)= \,& {} (\phi _{1}(x,t,z),\,\chi (x,t,z),\,\psi _{3}(x,t,z)),\quad z\in D^+\cup \Sigma \setminus \left\{ 0\right\} , \end{aligned}$$
(4.17a)
$$\begin{aligned} \Phi ^-(x,t,z)= \,& {} (\psi _{1}(x,t,z),\,{\bar{\chi }}(x,t,z),\,\phi _{3}(x,t,z)),\quad z\in D^-\cup \Sigma \setminus \left\{ 0\right\} , \end{aligned}$$
(4.17b)

with \(\Phi ^\pm (x,t,z)\) analytic for \(z\in D^\pm\), respectively. Furthermore, one can show that

$$\begin{aligned} \det \Phi ^+(x,t,z)= & {} -4q_0^2\, \lambda ^2(z)a_{11}(z)b_{33}(z)e^{i \theta _2(x,t,z)},\qquad D^+\cup \Sigma , \end{aligned}$$
(4.18a)
$$\begin{aligned} \det \Phi ^-(x,t,z)= \,& {} 4q_0^2\, \lambda ^2(z) a_{33}(z)b_{11}(z)e^{i \theta _2(x,t,z)},\qquad D^-\cup \Sigma , \end{aligned}$$
(4.18b)

where \(a_{11}(z)\) and \(b_{33}(z)\) are analytic in \(D^+\), while \(b_{11}(z)\) and \(a_{33}(z)\) are analytic in \(D^-\).

Symmetries. The Lax pair admits the usual two symmetries. The first symmetry corresponds to the mapping \(z \mapsto z^*\), under which in the defocusing case one has

$$\begin{aligned} A^\dag (z) = \Gamma (z)B(z)\Gamma ^{-1}(z),\qquad \Gamma (z) = \textrm{diag}(-q_o^2/z, 2\lambda (z), z), \qquad z \in {\mathbb {R}}. \end{aligned}$$
(4.19)

The symmetries for the analytic entries of the scattering matrix are then extended off the real z-axis using Schwarz reflection principle, and yield

$$\begin{aligned} b_{11}^*(z^*)=a_{11}(z), \quad z\in D^+, \qquad b_{33}^*(z^*)=a_{33}(z), \quad z\in D^-. \end{aligned}$$
(4.20)

The second symmetry corresponds to the mapping \(z \mapsto {\hat{z}}^*=q_0^2/z\), and it leads to the following relations for the scattering coefficients:

$$\begin{aligned} A({\hat{z}}^*) = \Pi \,A(z)\,\Pi ^{-1},\qquad \Pi = \begin{pmatrix}0 &{} 0 &{} -1\\ 0 &{} 1 &{} 0\\ -1 &{} 0 &{} 0\\ \end{pmatrix}, \quad z\in {\mathbb {R}}, \end{aligned}$$
(4.21)

giving

$$\begin{aligned} a_{11}({\hat{z}}^*)=a_{33}(z) \quad z\in D^-, \qquad b_{11}({\hat{z}}^*)=b_{33}(z) \quad z\in D^+. \end{aligned}$$
(4.22)

The reflection coefficients that will be needed in the inverse problem are:

$$\begin{aligned} \rho _1(z) = \frac{b_{31}(z)}{b_{11}(z)},\quad \rho _2(z) = \frac{a_{12}(z)}{a_{11}(z)},\quad {\bar{\rho }}_2(z) = \frac{b_{21}(z)}{b_{11}(z)}. \end{aligned}$$
(4.23a)

Note that only two of the above three coefficients are independent, since using the first symmetry one can show that \({\bar{\rho }}_2^*(z^*) = q_o^2\rho _2(z)/(q_0^2 - z^2)\).

Discrete spectrum. According to Eqs. (4.18), the columns of \(\Phi ^\pm (x,t,z)\) become linearly dependent at the zeros of \(a_{11}(z)\) and \(b_{33}(z)\), and those of \(a_{33}(z)\) and \(b_{11}(z)\), respectively, and these zeros are the discrete eigenvalues of the scattering problem. For simplicity, in the following we assume that none of the above scattering coefficients vanishes on the real z-axis (i.e., we exclude spectral singularities), and that all zeros are simple (the case of multiple discrete eigenvalues has been considered in [36]). Let \(D_o=\{z \in {\mathbb {C}}: |z| < q_o \}\) denote the disk of radius \(q_o\), so that \(C_o=\partial D_o\). Importantly, unlike the scalar and square matrix defocusing NLS, the zeros are not confined to the circle \(C_o\). Zeros on \(C_o\) correspond to proper discrete eigenvalues, which give rise to \(L^2\) eigenfunctions of the scattering problem, while zeros off \(C_o\) are generalized eigenvalues, which do not correspond to \(L^2\) eigenfunctions (cf. [149] for additional details). Discrete eigenvalues on the circle \(C_o\) correspond to dark-dark (DD) solitons (i.e., solitons that appear as localized dips of intensity on the background field \(q_o\) in both components), while discrete eigenvalues off \(C_o\) give rise to dark-bright (DB) solitons (i.e., solitons that appear as dark solitons in one component and bright in the other component), which have no analog in the scalar defocusing NLS. Then we label \(\zeta _n\) for \(n=1,\dots ,{\mathcal {N}}_1\) the zeros of \(a_{11}(z)\) on \(C_o\cap {\mathbb {C}}^+\), and \(z_n\) for \(n=1,\dots ,{\mathcal {N}}_2\) the zeros of \(a_{11}(z)\) in \(D_o\cap {\mathbb {C}}^+\). For each pair \(\zeta _n,\zeta _n^*\) of discrete eigenvalues on the circle, one then has \(a_{11}(\zeta _n) = a_{33}(\zeta _n^*) = 0\), and

$$\begin{aligned} \phi _{1}(x,t,\zeta _n)= & {} c_n \psi _{3}(x,t,\zeta _n),\quad \\ \phi _{3}(x,t,\zeta _n^*)= & {} {\bar{c}}_n \psi _{1}(x,t,\zeta _n^*), \qquad n = 1,\dots ,{\mathcal {N}}_1, \end{aligned}$$
(4.24a)

with

$$\begin{aligned} {\bar{c}}_n = c_n,\quad c_n^*=\frac{\zeta _n^*}{\zeta _n}\dfrac{ b_{11}'(\zeta _n^*)}{ a_{33}'(\zeta _n^*)}{\bar{c}}_n. \end{aligned}$$
(4.24b)

Similarly, for each quartet of eigenvalues \(z_n,z_n^*, {\hat{z}}_n, {\hat{z}}_n^*\) off the circle, one has \(a_{11}(z_n) = b_{11}(z_n^*) = a_{33}({\hat{z}}_n^*) = b_{33}({\hat{z}}_n) = 0\), and

$$\begin{aligned} \phi _{1}(x,t,z_n)= \,& {} d_n \chi (x,t,z_n) \equiv d_n {\bar{\chi }}(x,t,{\hat{z}}_n^*),\quad \\ {\bar{\chi }}(x,t,z_n^*)= \,& {} {\bar{d}}_n \phi _{+,1}(x,t,z_n^*) \equiv -{\bar{d}}_n \psi _{3}(x,t,{\hat{z}}_n), \end{aligned}$$
(4.25a)

for all \(n = 1,\dots ,{\mathcal {N}}_2\), with

$$\begin{aligned} {\bar{d}}_n^* = 2\,{\hat{z}}_n^*\,\lambda ({\hat{z}}_n^*)b_{11}({\hat{z}}_n^*)\,d_n. \end{aligned}$$
(4.25b)

All details can be found in [36, 149] (keeping in mind that here we use for the eigenfunctions the normalizations of [149], not those of [36]).

Asymptotics of eigenfunctions and scattering data. The asymptotic behavior of the Jost eigenfunctions as \(z\rightarrow 0\) and as \(z\rightarrow \infty\) in the proper domains of analyticity can be obtained via a WKB expansion on the scattering problem:

$$\begin{aligned} \phi _1(x,z,t)e^{-i\theta _1}\sim & {} \begin{pmatrix} z \\ i{\textbf{q}}^*(x,t)\end{pmatrix}, \qquad \\ \psi _3(x,z,t)e^{i\theta _1}\sim & {} \begin{pmatrix} {\textbf{q}}^T(x,t){\textbf{q}}_+^*/z \\ i{\textbf{q}}^*_+ \end{pmatrix} \qquad z\rightarrow \infty , \quad z\in D^+ \end{aligned}$$
(4.26a)
$$\begin{aligned} \phi _1(x,z,t)e^{-i\theta _1}\sim & {} \begin{pmatrix} z{\textbf{q}}^T(x,t){\textbf{q}}_-^*/q_o^2 \\ i{\textbf{q}}_-^*\end{pmatrix}, \qquad \\ \psi _3(x,z,t)e^{i\theta _1}\sim & {} \begin{pmatrix} q_o^2/z \\ i{\textbf{q}}^*(x,t)\end{pmatrix} \qquad z\rightarrow 0, \quad z\in D^+ \end{aligned}$$
(4.26b)
$$\begin{aligned} \phi _3(x,z,t)e^{i\theta _1}\sim & {} -\begin{pmatrix} {\textbf{q}}^T(x,t){\textbf{q}}_-^*/z \\ i{\textbf{q}}^*_-\end{pmatrix}, \qquad \\ \psi _1(x,z,t)e^{-i\theta _1}\sim & {} \begin{pmatrix} z \\ i{\textbf{q}}^*(x,t)\end{pmatrix} \qquad z\rightarrow \infty , \quad z\in D^- \end{aligned}$$
(4.26c)
$$\begin{aligned} \phi _3(x,z,t)e^{i\theta _1}\sim & {} -\begin{pmatrix} q_o^2/z \\ i{\textbf{q}}^*(x,t)\end{pmatrix}, \qquad \\ \psi _1(x,z,t)e^{-i\theta _1}\sim & {} \begin{pmatrix} z{\textbf{q}}^T(x,t){\textbf{q}}_+^*/q_o^2 \\ i{\textbf{q}}^*_+ \end{pmatrix} \qquad z\rightarrow 0, \quad z\in D^- \end{aligned}$$
(4.26d)

as well as

$$\begin{aligned} \chi (x,t,z)e^{-i\theta _2}\sim & {} -\begin{pmatrix} {\textbf{q}}^T(x,t){\textbf{q}}_-^\perp \\ iz{\textbf{q}}_-^\perp \end{pmatrix}, \quad \\ {\bar{\chi }}(x,t,z)e^{i\theta _2}\sim & {} \begin{pmatrix} {\textbf{q}}^T(x,t){\textbf{q}}_+^\perp \\ iz{\textbf{q}}_+^\perp \end{pmatrix}\qquad z\rightarrow \infty , \end{aligned}$$
(4.26e)
$$\begin{aligned} \chi (x,t,z)e^{-i\theta _2}\sim & {} \begin{pmatrix} {\textbf{q}}^T(x,t){\textbf{q}}_+^\perp \\ iq_o^2 {\textbf{q}}_+^\perp \end{pmatrix}, \quad \\ {\bar{\chi }}(x,t,z)e^{-i\theta _2}\sim & {} -\begin{pmatrix} {\textbf{q}}^T(x,t){\textbf{q}}_-^\perp \\ iq_o^2 {\textbf{q}}_-^\perp \end{pmatrix} \quad z\rightarrow 0, \end{aligned}$$
(4.26f)

from which one can then obtain the asymptotics of the analytic scattering coefficients in the appropriate domains:

$$\begin{aligned} a_{11}(z)\sim \left\{ \begin{array}{ll} 1 &{} \quad z\rightarrow \infty \\ {\textbf{q}}_+^T{\textbf{q}}_-^*/q_o^2 &{} \quad z\rightarrow 0 \end{array} \right. , \qquad b_{33}(z)\sim \left\{ \begin{array}{ll} {\textbf{q}}_-^T{\textbf{q}}_+^*/q_o^2 &{} \quad z\rightarrow \infty \\ 1 &{} \quad z\rightarrow 0 \end{array} \right. , \end{aligned}$$
(4.27a)

(the behavior of \(a_{33}(z)\) and \(b_{11}(z)\) follows from the symmetries (4.20)). Finally, as mentioned above, under mild decay conditions on the potential one can ensure the Jost eigenfunctions are continuous at \(z=\pm q_o\), while the scattering coefficients generically have poles at the branch points but such that the reflection coefficients are finite and unimodular. We refer to [36, 149] for details.

Inverse problem. The starting point for reconstructing the solution \({\textbf{q}}(x,t)\) from the scattering data is the scattering relation (4.9), which yields a jump condition for a matrix Riemann–Hilbert problem. Specifically, we introduce the sectionally meromorphic matrices

$$\begin{aligned} \mu ^+(x,t,z)= & {} \left[ \frac{\phi _1}{a_{11}}e^{-i\theta _1}, \, \frac{\chi }{2\lambda b_{33}}e^{-i\theta _2},\, \psi _3e^{i\theta _1}\right] , \qquad \\ \mu ^-(x,t,z)= & {} \left[ \psi _1 e^{-i\theta _1}, \, -\frac{{\bar{\chi }}}{2\lambda b_{11}}e^{-i\theta _2},\, \frac{\phi _3}{a_{33}}e^{i\theta _1}\right] , \end{aligned}$$
(4.28)

where superscripts \(^\pm\) distinguish between analyticity in \(D^+\) and \(D^-\), respectively. Then, using  (4.9) and (4.16a), as well as the symmetries of the scattering coefficients, then obtain the jump condition

$$\begin{aligned} \mu ^+(z)=\mu ^-(z)\left( I_3-G(z)\right) \qquad z\in \Sigma , \end{aligned}$$
(4.29)

where the jump matrix G(xtz) is given by

$$\begin{aligned} G(z)=\begin{pmatrix} \rho _1(z)\rho _1({\hat{z}}^*) &{} -[{\bar{\rho }}_2(z)+\rho _1(z){\bar{\rho }}_2({\hat{z}}^*)]e^{i(\theta _1-\theta _2)} &{} -\rho _1(z)e^{2i\theta _1} \\ [\rho _2(z)+\rho _1({\hat{z}}^*)\rho _2({\hat{z}}^*)]e^{i(\theta _2-\theta _1)} &{} \rho _2({\hat{z}}^*){\bar{\rho }}_2({\hat{z}}^*)e^{i\theta _1} &{} \rho _2({\hat{z}}^*)e^{i(\theta _1+\theta _2)} \\ \rho _1({\hat{z}}^*)e^{-2i\theta _1} &{} {\bar{\rho }}_2({\hat{z}}^*)e^{-i\theta _2} &{} 0\end{pmatrix}. \end{aligned}$$
(4.30)

[the (xt)-dependence in the eigenfunctions and in \(\theta _j\) has been omitted for brevity]. Subtracting the asymptotic behavior as \(z\rightarrow \infty\), the pole at \(z=0\) and at the discrete eigenvalues, one can then reduce the solution of the inverse problem to a set of integral equations for the following three modified vector eigenfunctions:

$$\begin{aligned} M_1(z)= \,& {} \psi _{3}(z) \,\textrm{e}^{i \theta _1(z)},\quad z\in {\mathbb {C}}^+, \end{aligned}$$
(4.31a)
$$\begin{aligned} M_2(z)= \,& {} \psi _{1}(z) \,\textrm{e}^{-i \theta _1(z)}, \quad M_3(z) = \dfrac{{\bar{\chi }}(z)}{2\lambda (z)\,b_{11}(z)}\,\textrm{e}^{-i\theta _2(z)},\quad z\in {\mathbb {C}}^-. \end{aligned}$$
(4.31b)

Note that \({\bar{\chi }}(z)/{(2\lambda (z)\,b_{11}(z))}\) remains analytic also at the zeros of \(b_{11}(z)\). Moreover, these eigenfunctions satisfy the integral equations

$$\begin{aligned} M_1(z)= & {} - \begin{pmatrix} {\hat{z}}^*\\ iq_+^* \end{pmatrix} + \sum _{n=1}^{{\mathcal {N}}_1} \frac{{\bar{\gamma }}_n\,e^{-2i \theta (\zeta _n)}}{z-\zeta _n^*} M_2(\zeta _n^*) - \sum _{n=1}^{{\mathcal {N}}_2} \frac{{\bar{{\delta }}}^*_n \, e^{i(\theta _2(z_n)- \theta _1(z_n))}}{z-{\hat{z}}_n^*} M_3({\hat{z}}_n^*) \\{} & {} +\frac{1}{2 \pi i} \mathop {\int }\limits _{-\infty }^{\infty } \frac{\textrm{d}\zeta }{\zeta - z}\Big [ \rho _1(\zeta ) M_2(\zeta )\,e^{2i\theta _1(\zeta )} \\{} & {} - \rho _2({\hat{\zeta }}^*) M_3(\zeta )\,e^{i(\theta _2(\zeta )+\theta _1(\zeta ))} \Big ], \quad z\in {\mathbb {C}}^+, \end{aligned}$$
(4.32a)
$$\begin{aligned} M_2(z)= & {} \begin{pmatrix} z\\ iq_+^* \end{pmatrix} + \sum _{n=1}^{{\mathcal {N}}_1} \frac{\gamma _n z\,e^{-2i \theta _1(\zeta _n)}}{z-\zeta _n} M_1(\zeta _n) - \sum _{n=1}^{{\mathcal {N}}_2} \frac{{\bar{{\delta }}}^*_n z\,e^{i(\theta _2(z_n)- \theta _1(z_n))}}{{\hat{z}}_n^*(z-z_n)} M_3({\hat{z}}^*_n) \\{} & {} -\frac{z}{2 \pi i} \mathop {\int }\limits _{-\infty }^{\infty } \frac{\textrm{d}\zeta }{\zeta (\zeta - z)} \Big [ \rho _1({\hat{\zeta }}^*)\,e^{-2i \theta _1(\zeta )}M_1(\zeta ) \\{} & {} + \rho _2(\zeta ) M_3({\hat{\zeta }}^*) e^{i (\theta _2(\zeta )-\theta _1(\zeta ))} \Big ], \quad z\in {\mathbb {C}}^-, \end{aligned}$$
(4.32b)
$$\begin{aligned} M_3(z)= & {} \begin{pmatrix} 0\\ iq_+^\perp \end{pmatrix} - \sum _{n=1}^{{\mathcal {N}}_2} \frac{{\bar{\delta }}_n z}{(z-{\hat{z}}_n)(z-z_n^*)} M_2(z_n^*) e^{i(\theta _1(z_n^*)-\theta _2(z_n^*))} \\{} & {} -\frac{1}{2 \pi i} \mathop {\int }\limits _{-\infty }^{\infty } \frac{\textrm{d}\zeta }{\zeta - z} [{\bar{\rho }}_2(\zeta )\,e^{i(\theta _1(\zeta )-\theta _2(\zeta ))} M_2(\zeta ) \\{} & {} + {\bar{\rho }}_2({\hat{\zeta }}^*)M_1(\zeta )e^{-i(\theta _1(\zeta )+\theta _2(\zeta ))}] ,\quad z \in {\mathbb {C}}^-, \end{aligned}$$
(4.32c)

where, as before, the (xt)-dependence was omitted for brevity. In the above system: (i) \(\gamma _n, {\bar{\gamma }}_n\) are the norming constants associated to the discrete eigenvalues \(\zeta _n, \zeta _n^*\) as given in (4.33) below [the first symmetry can then be used to express \(\gamma _n\) in terms of \({\bar{\gamma }}_n\), see (4.34)]; (ii) \(\delta _n, {\bar{\delta }}_n\) are the norming constants associated to the discrete eigenvalues \(z_n, z_n^*\) given in (4.33) [note that the second symmetry has been used to eliminate from the system the norming constants associated to the other 2 eigenvalues in the quartet; moreover, the first symmetry allows to express \({\bar{\delta }}_n\) in terms of \(\delta _n\), see (4.34) below]; (iii) reflection coefficients \(\rho _1(z)\), \(\rho _2(z)\) given by (4.23a). The norming constants are defined in terms of the coefficients in (4.25) and (4.24) as follows:

$$\begin{aligned} \gamma _n= & {} \frac{c_n}{\zeta _n a'_{11}(\zeta _n)},\qquad {\bar{\gamma }}_n = \frac{{\bar{c}}_n}{a'_{33}(\zeta _n^*)}, \end{aligned}$$
(4.33a)
$$\begin{aligned} \delta _n= & {} \frac{d_n}{z_n a'_{11}(z_n)},\qquad {\bar{\delta }}_n = - \frac{{\bar{d}}_n}{z_n^*b'_{11}(z_n^*)}, \end{aligned}$$
(4.33b)

and they can be shown to satisfy the symmetries

$$\begin{aligned} \gamma _n =-{\bar{\gamma }}_n/\zeta _n^*, \quad {\bar{\gamma }}_n^*=-{\bar{\gamma }}_n, \quad {\bar{\delta }}_n^*=2\,{\hat{z}}_n^*\lambda (z_n)b_{11}({\hat{z}}_n^*)\delta _n. \end{aligned}$$
(4.34)

Recall that \({\hat{z}}=q_o^2/z^*\) and note, in particular, that \({\bar{\gamma }}_n\) must be purely imaginary.

Similarly to the scalar and square matrix case, one can obtain a trace formula for the coefficient \(a_{11}(z)\) in terms of scattering data, and also obtain as a byproduct the \(\theta\)-condition relating the scattering data to the asymptotic phase difference of the potential, respectively given by:

$$\begin{aligned} \displaystyle a_{11}(z)= & {} \prod _{n=1}^{{\mathcal {N}}_1} \frac{z - \zeta _n}{z - \zeta _n^*}\,\prod _{n=1}^{{\mathcal {N}}_2} \frac{z - z_n}{z - z_n^*}\,\exp \Big \{ -\frac{1}{2 \pi i} \mathop {\int }\limits _{-\infty }^{\infty }\dfrac{\log [1- R(\zeta )]}{\zeta -z}d \zeta \Big \}, \end{aligned}$$
(4.35a)
$$\begin{aligned} \displaystyle e^{i\Delta \,\theta }= & {} \prod _{n=1}^{{\mathcal {N}}_1} \frac{\zeta _n}{\zeta _n^*}\,\prod _{n=1}^{{\mathcal {N}}_2} \frac{z_n}{z_n^*}\,\exp \Big \{ -\frac{1}{2 \pi i} \mathop {\int }\limits _{-\infty }^{\infty }\dfrac{\log [1- R(\zeta )]}{\zeta }d\zeta \Big \}, \end{aligned}$$
(4.35b)

with \(\rho _1(z),\,\rho _2(z)\) as in (4.23a) and

$$\begin{aligned} R(z) = \frac{z^2}{q_0^2}|\rho _1(z)|^2 + \frac{q_0^2}{(z^2 - q_0^2)} |\rho _2(z)|^2,\qquad z \in {\mathbb {R}}. \end{aligned}$$
(4.36)

The trace formula (4.35a) and the symmetry \(b_{11}(z)=a_{11}^*(z^*)\) (cf. (4.19)) then allow to express \(\delta _n\) in terms of \({\bar{\delta }}_n\) (cf. (4.34)) in the linear system (4.32).

Finally, one can reconstruct the potential \({\textbf{q}}(x,t)\) of the Manakov system from \(M_1,M_2,M_3\) by comparing their asymptotics to that of the eigenfunctions obtained from the direct scattering problem, obtaining

$$\begin{aligned} {\textbf{q}}(x,t)= \,& {} {\textbf{q}}_+ +i \sum _{n=1}^{{\mathcal {N}}_1} \frac{{\bar{\gamma }}_n}{\zeta _n^*} \left[ M_{1}^{\textrm{dn}}(x,t,\zeta _n)\right] ^* e^{-2i \theta _1(x,t,\zeta _n)} \\{} & {} - i \sum _{n=1}^{{\mathcal {N}}_2} \frac{{\bar{{\delta }}}_n}{{\hat{z}}_n} \left[ M_3^{\textrm{dn}}(x,t,{\hat{z}}_n^*)\right] ^* e^{i (\theta _1(x,t,z_n^*)-\theta _2(x,t,z_n^*))} \\{} & {} -\frac{1}{2 \pi } \mathop {\int }\limits _{-\infty }^{\infty } \Big \{\rho _1^*({\hat{\zeta }}) \left[ M^{\textrm{dn}}_{1}(x,t,\zeta )\right] ^* e^{2i \theta _1(x,t,\zeta )} \\{} & {} -\rho _2^*(\zeta ) \left[ M_3^{\textrm{dn}}(x,t,{\hat{\zeta }})\right] ^* e^{i (\theta _1(x,t,\zeta )-\theta _2(x,t,\zeta ))} \Big \} \frac{\textrm{d}\zeta }{\zeta }, \end{aligned}$$
(4.37)

where the superscript \(^{\textrm{dn}}\) denotes the lower two components of the vector eigenfunctions.

To the best of our knowledge, obtaining a Gelfand–Levitan–Marchenko formulation of the inverse problem for the defocusing Manakov system with NZBC is an open problem. Difficulties arise due to the lack of analyticity of the Jost eigenfunctions \(\phi _2,\psi _2\). The defect of analyticity is circumvented by introducing the auxiliary eigenfunctions \(\chi ,{\bar{\chi }}\), but the latter are not scattering eigenfunctions, i.e., they do not have a simple, plane-wave behavior as \(x\rightarrow \pm \infty\), and hence they likely do not admit triangular representations analogous to (2.44).

Solitons. In addition to dark solitons, which are the solutions (2.1) of the scalar defocusing NLS multiplied by an arbitrary constant unit vector \({\textbf{p}}\), the defocusing Manakov system admits dark-bright (DB) solitons:

$$\begin{aligned} q_1(x,t)= \,& {} q_oe^{i(\theta +\alpha )}\left[ \cos \alpha -i\sin \alpha \tanh (\nu (x-2\xi t-x_o))\right] \end{aligned}$$
(4.38a)
$$\begin{aligned} q_2(x,t)= & {} -i\sin \alpha \sqrt{q_o^2-|z_1|^2}e^{i(\xi x-(\xi ^2-\nu ^2)t+\phi )}\textrm{sech}(\nu (x-2\xi t-x_o)) \end{aligned}$$
(4.38b)

which have no analog in the scalar case. The DB solution is expressed in terms of the discrete eigenvalue \(z_1=|z_1|e^{i\alpha }\equiv \xi +i \nu\) with \(-q_o<\xi <q_o\) and \(0<\nu <q_o\), so \(|z_1|<q_o\). Note that a DD soliton is obtained when the discrete eigenvalue is chosen as \(\zeta _1=q_o e^{e^{i\alpha }}\).

DD and DB soliton solutions for the defocusing Manakov system had been discovered by direct methods [116, 156, 165], and further investigated in various works [75, 136, 138, 143], but the IST developed in [149] actually provided their spectral characterization. DB soliton states were first observed in pioneering experiments in optics [53, 54] and more recently in binary BECs [25, 48, 99] (see [112] for an excellent review about solitons in coupled nonlinear Schrödinger models). Note that, unlike vector bright solitons of the focusing Manakov system, DD and DB defocusing solitons interact trivially, they do not exhibit any polarization shift [58]. We also mention [161], where theoretical and numerical evidence is provided for the potential realization of the Peregrine soliton in repulsive two-component BECs, by reducing a two-component defocusing (i.e., repulsive interaction) setting to an effective focusing single-component one for a minority component in the condensate.

4.1.2 Focusing Manakov System

Unlike the defocusing case, in the focusing Manakov system all of the columns of \(\Phi (x,t,z)\) and \(\Psi (x,t,z)\) are analytic in some portion of the complex z-plane (cf. Theorem 4.1). Below, we describe how to complete the formulation of the direct problem and adapt the inverse problem for the focusing Manakov system. We follow [121], but using the normalizations of [149] for the eigenfunctions, so as to provide a more unified description.

Auxiliary eigenfunctions. One can use again the columns of the adjoint eigenfunctions (4.13) to define four new solutions (auxiliary eigenfunctions) of the original Lax pair (4.4a):

$$\begin{aligned} \chi _1(x,t,z)= \,& {} \textrm{e}^{i\theta _2(x,t,z)} [{\tilde{\psi }}_{1} \times {\tilde{\phi }}_{2}](x,t,z), \qquad \\ \chi _2(x,t,z)= \,& {} \textrm{e}^{i\theta _2(x,t,z)} [{\tilde{\phi }}_{1} \times {\tilde{\psi }}_{2}](x,t,z), \end{aligned}$$
(4.39a)
$$\begin{aligned} \chi _3(x,t,z)= \,& {} \textrm{e}^{i\theta _2(x,t,z)} [{\tilde{\psi }}_{2} \times {\tilde{\phi }}_{3}](x,t,z), \qquad \\ \chi _4(x,t,z)= \,& {} \textrm{e}^{i\theta _2(x,t,z)} [{\tilde{\phi }}_{2} \times {\tilde{\psi }}_{3}](x,t,z). \end{aligned}$$
(4.39b)

Note that four different auxiliary eigenfunctions are introduced, in contrast to the defocusing case, where only two auxiliary eigenfunctions were required. By construction, each auxiliary eigenfunction \(\chi _j(x,t,z)\) is analytic for \(z \in D_j\). Then one can define \(3 \times 3\) matrices

$$\begin{aligned} \Phi _1(x,t,z)= \,& {} (\phi _{1}(x,t,z),\psi _{2}(x,t,z),\chi _1(x,t,z)), \qquad z \in D_1, \end{aligned}$$
(4.40a)
$$\begin{aligned} \Phi _2(x,t,z)= \,& {} (\psi _{1}(x,t,z),\phi _{2}(x,t,z),\chi _2(x,t,z)), \qquad z \in D_2, \end{aligned}$$
(4.40b)
$$\begin{aligned} \Phi _3(x,t,z)= \,& {} (\chi _3(x,t,z),\phi _{2}(x,t,z),\psi _{3}(x,t,z)), \qquad z \in D_3, \end{aligned}$$
(4.40c)
$$\begin{aligned} \Phi _4(x,t,z)= \,& {} (\chi _4(x,t,z),\psi _{2}(x,t,z),\phi _{3}(x,t,z)), \qquad z \in D_4, \end{aligned}$$
(4.40d)

each of which is analytic in one of the four fundamental domains. Also, under the same assumptions as in Theorem 4.1, one can show that the following scattering coefficients can be analytically extended off of \(\Sigma\) in the following regions:

$$\begin{aligned}{} & {} a_{11}: \, \, z \in D_1, \qquad a_{22}: \, \, \mathop {\textrm{Im}}\nolimits z <0, \qquad a_{33}: \, \, z \in D_4, \end{aligned}$$
(4.41a)
$$\begin{aligned}{} & {} b_{11}: \, \, z \in D_2, \qquad b_{22}: \, \, \mathop {\textrm{Im}}\nolimits z >0, \qquad b_{33}: \, \, z \in D_3. \end{aligned}$$
(4.41b)

Symmetries. Like in all previous cases, we need to consider two symmetries: the transformation \(z\mapsto z^*\) (mapping the upper-half plane into the lower-half plane and viceversa), and \(z\mapsto -q_o^2/z\) (mapping the exterior of the circle \(C_o\) into the interior, and viceversa). For brevity, we will omit the symmetries for the eigenfunctions (see [121] for details), and only provide the induced symmetries in the scattering data which are needed in the inverse problem. Adapting the results in [121] to the normalizations used here, one finds that as a result of the first symmetry the scattering matrix and its inverse satisfy the symmetry relation:

$$\begin{aligned} A^\dagger (z^*) = C(z) B(z) C^{-1}(z), \qquad C(z) = \textrm{diag}(z^2\gamma (z),q_o^2,q_o^2\gamma (z)) \qquad z \in \Sigma , \end{aligned}$$
(4.42)

where \(\gamma (z)=1+q_o^2/z^2\). In particular, the analytic scattering coefficients are such that

$$\begin{aligned} b_{11}(z)= \,& {} a_{11}^{*}(z^{*}) \quad z \in D_2, \qquad b_{22}(z)=a_{22}^{*}(z^{*}) \qquad \mathop {\textrm{Im}}\nolimits z>0, \qquad b_{33}(z) = a_{33}^{*}(z^{*}) \qquad z \in D_3. \end{aligned}$$
(4.43)

Furthermore, the second symmetry implies

$$\begin{aligned} A(-q_o^2/z) = {\tilde{\Pi }}(z) A(z) {\tilde{\Pi }}^{-1}(z), \qquad {\tilde{\Pi }}(z) = \begin{pmatrix} 0 &{} 0 &{} -1 \\ 0 &{} 1 &{} 0 \\ 1 &{} 0 &{} 0 \end{pmatrix} \qquad z \in \Sigma . \end{aligned}$$
(4.44)

The analyticity of the scattering matrix then allows to conclude that

$$\begin{aligned} a_{11}(z)= \,& {} a_{33}(-q_o^2/z), \quad z \in D_1, \qquad b_{11}(z) = b_{33}(-q_o^2/z), \quad z \in D_2, \end{aligned}$$
(4.45a)
$$\begin{aligned} b_{22}(z)= \,& {} b_{22}(-q_o^2/z), \quad \mathop {\textrm{Im}}\nolimits z \ge 0, \qquad a_{22}(z) = a_{22}(-q_o^2/z), \quad \mathop {\textrm{Im}}\nolimits z \le 0. \end{aligned}$$
(4.45b)

Finally, the following reflection coefficients will appear in the inverse problem:

$$\begin{aligned} \rho _1(z) = \frac{a_{21}(z)}{a_{11}(z)}, \quad \rho _2(z) = \frac{a_{31}(z)}{a_{11}(z)}, \quad \rho _3(z) = \frac{a_{32}(z)}{a_{22}(z)}. \end{aligned}$$
(4.46)

Using the symmetries one can express \(\rho _3(z)\) in terms of \(\rho _1(z),\rho _2(z)\) showing that, as expected, only two of the reflection coefficients are independent.

Discrete spectrum. Recalling (4.40), one can show that

$$\begin{aligned}{} & {} \textrm{Wr}\,\Phi _1(x,t,z) \propto a_{11}(z)b_{22}(z) \,\textrm{e}^{i\theta _2(x,t,z)}, \quad z \in D_1, \end{aligned}$$
(4.47a)
$$\begin{aligned}{} & {} \textrm{Wr}\,\Phi _2(x,t,z) \propto a_{22}(z)b_{11}(z) \,\textrm{e}^{i\theta _2(x,t,z)}, \quad z \in D_2, \end{aligned}$$
(4.47b)
$$\begin{aligned}{} & {} \textrm{Wr}\,\Phi _3(x,t,z) \propto a_{22}(z)b_{33}(z) \,\textrm{e}^{i\theta _2(x,t,z)}, \quad z \in D_3, \end{aligned}$$
(4.47c)
$$\begin{aligned}{} & {} \textrm{Wr}\,\Phi _4(x,t,z) \propto a_{33}(z)b_{22}(z) \,\textrm{e}^{i\theta _2(x,t,z)}, \quad z \in D_4, \end{aligned}$$
(4.47d)

so that in each of the domains the columns of \(\Phi _j(x,t,z)\) become linearly dependent at the zeros of \(a_{11}(z)\), \(b_{22}(z)\), etc. Based on the symmetries (4.43) and (4.45), the above relationships imply that the discrete eigenvalues appear in symmetric quartets \(\left\{ z_j, z_j^*, -q_o^2/z_j, -q_o^2/z_j^*\right\}\). In [121], discrete eigenvalues are classified as belonging to the first, second or third type depending on whether they are zeros of \(a_{11}(z)\), \(b_{22}(z)\) or both, and the following theorem is proven.

Theorem 4.2

Let \(z_o\in D_1\) be a discrete eigenvalue of the scattering problem. That is, \(a_{11}(z_o)b_{22}(z_o) = 0\). Then the following are true:

  1. (i)

    If \(z_o\) is an eigenvalue of the first type (i.e., such that \(a_{11}(z_o)=0\) but \(b_{22}(z_o)\ne 0\)), there exist constants \(c_o\), \({\hat{c}}_o\), \(\check{c}_o\), and \({\overline{c}}_o\) such that

    $$\begin{aligned} \phi _{1}(x,t,z_o)&= c_o \chi _1(x,t,z_o)/b_{22}(z_o), \quad \chi _2(x,t,z_o^*) = {\hat{c}}_o \psi _{1}(x,t,z_o^*), \\ \chi _3(x,t,-q_o^2/z_o^*)&= \check{c}_o \psi _{3}(x,t,-q_o^2/z_o^*), \quad \phi _{3}(x,t,-q_o^2/z_o) = {\overline{c}}_o \chi _4(x,t,-q_o^2/z_o). \end{aligned}$$
  2. (ii)

    If \(z_o\) is an eigenvalue of the second type (i.e., such that \(b_{11}(z_o)=0\) but \(a_{22}(z_o)\ne 0\)), there exist constants \(d_o\), \({\hat{d}}_o\), \(\check{d}_o\), and \({\overline{d}}_o\) such that

    $$\begin{aligned} \chi _1(x,t,z_o)&= d_o\psi _{2}(x,t,z_o), \quad \phi _{2}(x,t,z_o^*) = {\hat{d}}_o\chi _2(x,t,z_o^*), \\ \phi _{2}(x,t,-q_o^2/z_o^*) &= \check{d}_o\chi _3(x,t,-q_o^2/z_o^*), \quad \chi _4(x,t,-q_o^2/z_o) = {\overline{d}}_o\psi _{2}(x,t,-q_o^2/z_o). \end{aligned}$$
  3. (iii)

    If \(z_o\) is an eigenvalue of the third type (i.e., such that \(a_{11}(z_o)=b_{22}(z_o)= 0\)), then \(\chi _1(x,t,z_o)=\chi _2(x,t,z_o^*)= {\textbf{0}}\), and there exist constants \(f_o\), \({\hat{f}}_o\), \(\check{f}_o\), and \({\overline{f}}_o\) such that

    $$\begin{aligned} \phi _{1}(x,t,z_o)= & {} f_o \psi _{2}(x,t,z_o), \quad \phi _{2}(x,t,z_o^*) = {\hat{f}}_o \psi _{1}(x,t,z_o^*), \\ \phi _{2}(x,t,-q_o^2/z_o^*)= & {} \check{f}_o \psi _{3}(x,t,-q_o^2/z_o^*), \quad \phi _{3}(x,t,-q_o^2/z_o) = {\overline{f}}_o \psi _{2}(x,t,-q_o^2/z_o). \end{aligned}$$

Note that, in contrast to the defocusing case, where zeros of the analytic scattering coefficients off \(C_o\) do not lead to bound states, in the focusing Manakov system all discrete eigenvalues correspond to \(L^2\) eigenfunctions.

Inverse problem. In order to formulate the inverse problem as a RHP, one defines the piecewise meromorphic function \(\mu (x,t,z)\) as \(\mu (x,t,z) = \mu _j(x,t,z)\) for \(z \in D_j\) (\(j=1,\ldots ,4\)), where

$$\begin{aligned} \mu _1(x,t,z)= \,& {} \Phi _1 \textrm{e}^{- i \Theta } [\mathop {\textrm{diag}}\nolimits (a_{11},1,b_{22})]^{-1} \quad z\in D_1, \end{aligned}$$
(4.48a)
$$\begin{aligned} \mu _2(x,t,z)= \,& {} \Phi _2 \textrm{e}^{- i \Theta } [\mathop {\textrm{diag}}\nolimits (1,a_{22},b_{11})]^{-1} \quad z\in D_2, \end{aligned}$$
(4.48b)
$$\begin{aligned} \mu _3(x,t,z)= \,& {} \Phi _3 \textrm{e}^{- i \Theta } [\mathop {\textrm{diag}}\nolimits (b_{33},a_{22},1)]^{-1} \quad z\in D_3, \end{aligned}$$
(4.48c)
$$\begin{aligned} \mu _4(x,t,z)= \,& {} \Phi _4 \textrm{e}^{- i \Theta } [\mathop {\textrm{diag}}\nolimits (b_{22},1,a_{33})]^{-1} \quad z\in D_4. \end{aligned}$$
(4.48d)

Then \(\mu _j(x,t,z)\) satisfy the jump conditions

$$\begin{aligned} \mu ^+(x,t,z) = \mu ^-(x,t,z) [I_3 - \textrm{e}^{i \Theta (x,t,z)}{\textbf{L}}(z) \textrm{e}^{- i \Theta (x,t,z)}], \quad z \in \Sigma , \end{aligned}$$
(4.49)

where \(\mu = \mu ^+\) for \(z \in D^+ = D_1 \cup D_3\) and \(\mu = \mu ^-\) for \(z \in D^- = D_2 \cup D_4\) (namely, \(\mu ^+ = \mu _1\) for \(z \in D_1\), \(\mu ^+ = \mu _3\) for \(z \in D_3\), \(\mu ^- = \mu _2\) for \(z \in D_2\), and \(\mu ^- = \mu _4\) for \(z \in D_4\)) and where the superscripts ± denote, respectively, projections from the left and the right of the appropriate contour in the complex z-plane. Here \(\Sigma = \Sigma _1\cup \Sigma _2\cup \Sigma _3\cup \Sigma _4\), where \(\Sigma _j\) is the boundary of \({\overline{D}}_j \cap {\overline{D}}_{j+1\!\mod 4}\) (oriented so that \(D^+\) is always to its left), and the matrix \({\textbf{L}}(x,t,z)\) is specified on each portion of the contour in terms of the reflection coefficients. The expression of the jump matrix \({\textbf{L}}(x,t,z)\) is quite lengthy, we refer to [121] for details.

The above RHP can be formally solved by converting it into a mixed system of algebraic-integral equations by subtracting the asymptotic behavior at infinity, by regularizing (i.e., subtracting the pole at \(z=0\) and any pole contributions from the discrete spectrum), and then applying Cauchy projectors. Finally, the potential \({\textbf{q}}(x,t)\) is recovered in terms of formal solution of the RHP for the eigenfunctions by taking into account that

$$\begin{aligned} \psi _1(x,t,z)e^{-i\theta _1(x,t,z)}\sim \begin{pmatrix} 1 \\ -i{\textbf{q}}^*(x,t)\end{pmatrix} \qquad z\rightarrow \infty , z\in D_1. \end{aligned}$$
(4.50)

Again, we omit the explicit expressions for brevity, referring to [121] for details.

Like in the defocusing case, a Gelfand–Levitan–Marchenko formulation of the inverse problem for the focusing Manakov system with NZBC is currently not available, to the best of our knowledge. Since in the focusing case all columns of the Jost eigenfunctions \(\Phi ,\Psi\) are analytic in some portion of the complex z-plane, it is more likely that a suitable generalization of the triangular representations (2.44) exist.

Solitons. Each discrete eigenvalue of type I,II,III produces a soliton solution of the focusing Manakov system. Solutions of type I correspond to soliton solutions of the scalar focusing NLS equation with NZBC described in Sect. 2 multiplied by a constant vector \({\textbf{p}}\). Solutions of type II and type III are the analogue of the dark-bright soliton solutions of the defocusing Manakov system with NZBC in (4.38), and they differ from each other only for the maximum of the bright component (in fact, solutions of type III can be obtained formally by taking the analytic continuation of solutions of type II when the eigenvalue is taken inside the circle \(C_o\) as opposed to outside, upon proper redefinition of the norming constants). In [121] only simple discrete eigenvalues were considered. Multiple pole solutions were recently constructed by using the dressing method to study resonant interactions of breathers [158]. To the best of our knowledge, solitons of the focusing Manakov system on a nontrivial background have not been experimentally observed so far.

4.2 3-component Defocusing NLS with NZBC

If one considers the Lax pair (4.4a) with \({\textbf{q}}(x,t)\) an m-component vector and \(m>2\), the eigenvalue ik in the scattering problem becomes a multiple eigenvalue with multiplicity \(m-1\). In the focusing case, one still has \(m-1\) Jost eigenfunctions and \(m-1\) adjoint Jost eigenfunctions which are analytic for \(\mathop {\textrm{Im}}\nolimits z<0\), and \(m-1\) which are analytic for \(\mathop {\textrm{Im}}\nolimits z>0\). Therefore, there is no conceptual difficulty in generalizing the construction in the preceding section, namely using cross products of each of the \(2(m-1)\) adjoint eigenfunctions which are analytic for \(z\in {\mathbb {C}}^\pm\) with \({\tilde{\phi }}_j\) and \({\tilde{\psi }}\) for \(j=1,3\) like in (4.39), to obtain 2m independent eigenfunctions \(\chi _j(x,t,z)\) which are analytic in one of the domains \(D_j\) for \(j=1,\dots ,4\), and then combine them like in (4.40) to have \((m+1)\times (m+1)\) eigenfunctions \(\Phi _j\), \(j=1,\dots ,4\) each of which is analytic in one of the four fundamental domains. This problem has been recently addressed in [130]. On the other hand, in the m-component defocusing VNLS the number of non-analytic eigenfunctions increases with m (specifically, \(2(m-1)\) out of the \(2(m+1)\) scattering eigenfunctions are not analytic, and the cross product construction only allows one to build 2 auxiliary analytic eigenfunctions. Hence, a fundamentally new approach is needed to generalize the IST with NZBC for the defocusing VNLS with \(m>2\) components. This problem was solved in [37, 150], where the tensor approach used by Beals, Deift and Tomei for m-th order scattering problems [24] was generalized to deal with nonzero boundary conditions and degenerate eigenvalues.

Broadly speaking, this approach is consistent with the usual development of the IST. Namely, for the direct problem: (i) Find complete sets of sectionally meromorphic eigenfunctions for the scattering operator that are characterized by their asymptotic behavior. (ii) Identify a minimal set of data that describes the relations among these eigenfunctions, and which therefore defines the scattering data. For the inverse problem: reconstruct the scattering operator (and in particular the potential) from its scattering data. Specific features of this approach, including those related to the extension of Beals, Deift, and Tomei’s work, are the following ones: (a) Typically, eigenfunctions and scattering data are defined for values of the scattering parameter in the continuous spectrum (e.g., the real axis in the case of decaying potentials, and \(\Sigma\) for NZBC), and are then extended to the complex plane. In this approach, the eigenfunctions are first defined away from the continuous spectrum, and the appropriate limits as the scattering parameter approaches the continuous spectrum are then evaluated. (b) The use of forms (tensors constructed by wedge products of columns of the matrix eigenfunctions), which simplifies the investigation of the analyticity properties of the eigenfunctions by reducing it to the study of Volterra equations. (c) Departure from \(L^2\)-theory: as already pointed out in Sect. 4.1, bounded eigenfunctions are insufficient to characterize the discrete spectrum in the defocusing case when the order of the scattering operator exceeds two.

We give below a succinct review of the IST of the 3-component defocusing VNLS equation:

$$\begin{aligned} i {\textbf{q}}_t + {\textbf{q}}_{xx} - 2 ( \Vert {\textbf{q}}\Vert ^2 - q_o^2 ) {\textbf{q}} = {\textbf{0}}, \end{aligned}$$
(4.51)

where \({\textbf{q}}(x,t) = (q_1,q_2,q_3)^T\), and with NZBC \(\lim _{x \rightarrow \pm \infty } {\textbf{q}}(x,t) = {\textbf{q}}_\pm\), \(\Vert {\textbf{q}}_\pm \Vert =q_o>0\). The generalization to more than 3 components can be obtained following the general ideas of [37, 150] and, although the details can be cumbersome, we do not expect significant conceptual differences.

Lax pair and uniformization coordinate. The VNLS Eq. (4.51) admits the following \(4\times 4\) Lax pair:

$$\begin{aligned} \varphi _x = {\textbf{U}}\,\varphi , \qquad \varphi _t = {\textbf{V}}\,\varphi , \end{aligned}$$
(4.52a)

where

$$\begin{aligned} {\textbf{U}}(x,t,k)= & {} -ik J + Q, \qquad \\ {\textbf{V}}(x,t,k)= \,& {} 2ik^2 J - i J (Q_x - Q^2 +q_o^2 ) - 2k Q, \end{aligned}$$
(4.52b)
$$\begin{aligned} J= & {} \begin{pmatrix} 1&{}{\textbf{0}}^T \\ {\textbf{0}}&{}-I_3 \end{pmatrix}\!,\qquad Q(x,t) = \begin{pmatrix} 0&{} {\textbf{q}}^\dagger \\ {\textbf{q}} &{} O \end{pmatrix}\!. \end{aligned}$$
(4.52c)

We assume that the boundary conditions \({\textbf{q}}_\pm\) at \(x\rightarrow \pm \infty\) are parallel, the generalization to non-parallel and/or asymmetric (\(\Vert {\textbf{q}}_+\Vert \ne \Vert {\textbf{q}}_-\Vert\)) boundary conditions is a completely open problem. In this, case, thanks to the U(N) invariance of Eq. (4.51), without loss of generality they can be chosen in the form \({\textbf{q}}_\pm = (0,0,q_\pm )^T\), with \(q_\pm = q_o\,\textrm{e}^{i \theta _\pm }\) and \(\theta _\pm\) arbitrary real constants.

As for the scalar defocusing NLS and Manakov system, the scattering problem is formally self-adjoint, but the continuous spectrum \({\mathbb {R}}\setminus (-q_o,q_o)\) exhibits a gap, and the Jost solutions are expressed in terms of \(\lambda = (k^2-q_o^2)^{1/2}\). As before, to deal effectively with the branching of \(\lambda\), we use the uniformization variable \(z= k + \lambda\) defined in (2.7) (cf. Fig. 1), with the original variables given by \(k= {\textstyle \frac{1}{2}}(z+q_o^2/z)\) and \(\lambda = {\textstyle \frac{1}{2}}(z-q_o^2/z)\). Expressing all functional dependence on k and \(\lambda\) in the IST in terms of z then eliminates the branching.

Jost solutions and analyticity and auxiliary eigenfunctions. For all \(z \in {\mathbb {R}}\) one defines the Jost solutions \(\Phi (x,t,z)\) and \(\Psi (x,t,z)\) as the simultaneous solutions of both parts of the Lax pair with the free-particle asymptotic behavior

$$\begin{aligned} \Phi (x,t,z)= & {} \left[ E_-(z) + o(1) \right] \textrm{e}^{i \Theta (x,t,z)}, \qquad x \rightarrow - \infty , \qquad \\ \Psi (x,t,z)= & {} \left[ E_+(z) + o(1) \right] \textrm{e}^{i \Theta (x,t,z)}, \qquad x \rightarrow + \infty , \end{aligned}$$
(4.53)

with \(\Theta (x,t,z) = \mathop {\textrm{diag}}\nolimits (\theta _1,\dots ,\theta _4) = \Lambda x - \Omega t\), where

$$\begin{aligned} i\Lambda (z)= \,& {} i\mathop {\textrm{diag}}\nolimits (-\lambda ,k,k,\lambda ), \\ -i\Omega (z)= & {} -i\mathop {\textrm{diag}}\nolimits (-2k \lambda ,k^2 {+} \lambda ^2, k^2{+}\lambda ^2, 2k \lambda ) \end{aligned}$$

are, respectively, the eigenvalue matrices of \(U_\pm = \lim _{x\rightarrow \pm \infty }U\) and \(V_\pm =\lim _{x\rightarrow \pm \infty }V\), and

$$\begin{aligned} E_\pm (z) = I + J\,Q_\pm /(iz), \end{aligned}$$
(4.54)

are the eigenvector matrices. As discussed previously, only the first and last columns of \(\Phi ,\Psi\) admit analytic continuation onto the complex z-plane. To circumvent this problem, in [150] a fundamental set of meromorphic eigenfunctions \(\Xi ^\pm (x,t,z)\) in each half plane was constructed using an appropriate extension of the scattering problem to higher-dimensional tensors. [Hereafter, subscripts ± denote normalization as \(x\rightarrow \pm \infty\), whereas superscripts ± denote analyticity or meromorphicity in the upper/lower half of the complex z-plane.]

For \(z\in {\mathbb {R}}\), the meromorphic eigenfunctions can be written in terms of the Jost eigenfunctions as follows:

$$\begin{aligned} \Xi ^+(x,t,z)= \,& {} \Phi (x,t,z) \alpha (z) = \Psi (x,t,z) \beta (z), \end{aligned}$$
(4.55a)
$$\begin{aligned} \Xi ^-(x,t,z)= \,& {} \Phi (x,t,z) {\tilde{\beta }}(z) = \Psi (x,t,z) {\tilde{\alpha }}(z), \end{aligned}$$
(4.55b)

where \(\alpha\) and \({\tilde{\alpha }}\) are upper triangular matrices, while \(\beta\) and \({\tilde{\beta }}\) are lower triangular ones. (In particular, the diagonal entries of \(\alpha\) and \({\tilde{\beta }}\) are all unity.) We refer the reader to [37, 150] for details.

Scattering coefficients. Since \(\det \Phi =\det \Psi = \gamma (z)\,\textrm{e}^{2i[kx - (k^2 + \lambda ^2)t]}\), where \(\gamma (z) = \det E_\pm = 1 - q_o^2/z^2\), for \(z\in {\mathbb {R}}\setminus \{\pm q_o\}\) one can introduce the scattering matrix \(A(z) = (a_{ij})\) via

$$\begin{aligned} \Phi (x,t,z) = \Psi (x,t,z) A(z), \end{aligned}$$
(4.56)

with \(\det A(z) = 1\). Also, since the Jost solutions solve the t-part of the Lax pair as well, all spectral data are time-independent.

Using Eq. (4.56) one can obtain triangular decompositions of the scattering matrix \(A(z) = \beta (z) \alpha ^{-1}(z) = {\tilde{\alpha }}(z) {\tilde{\beta }}^{-1}(z)\), similarly to the N-wave interactions [137]. In turn, these decompositions allow one to express the entries of \(\alpha ,\dots ,{\tilde{\beta }}\) in terms of the minors of A, denoted as

$$\begin{aligned} A_{{\big (\begin{array}{c} i_1,\dots ,i_p\\ k_1,\dots ,k_p \end{array}\big )}} = \det \begin{pmatrix} a_{i_1k_1} &{} \dots &{} a_{i_1k_p} \\ \vdots &{} \ddots &{} \vdots \\ a_{i_pk_1} &{} \cdots &{} a_{i_pk_p} \end{pmatrix}, \end{aligned}$$

where \(1 \leqslant i_1< \cdots < i_p \leqslant 4\) and similarly for \(k_1,\dots ,k_p\). In particular, the upper and lower principal minors of A are, respectively, determinants of the form

$$\begin{aligned} A_{[1,\dots ,p]} = A_{{\big (\begin{array}{c} 1,\dots ,p\\ 1,\dots ,p \end{array}\big )}},~ A_{[p,\dots ,N]} = A_{{\big (\begin{array}{c} p,\dots ,4\\ p,\dots ,4 \end{array}\big )}}, \quad 1 \leqslant p \leqslant 4. \end{aligned}$$

Importantly, Eqs. (4.55) allow one to obtain the analyticity properties of the minors of A. The upper principal minors \(A_{[1]},A_{[1,2]},A_{[1,2,3]}\) are analytic in the upper half plane (UHP), while the lower principal minors \(A_{[2,3,4]},A_{[3,4]},A_{[4]}\) are analytic in the lower half plane (LHP). In addition, the following minors are also analytic:

$$\begin{aligned} A_{{\left (\begin{array}{c} 1,2\\ 1,3 \end{array}\right )}}, A_{{\left (\begin{array}{c} 1,3\\ 1,2 \end{array}\right )}}:~ \mathop {\textrm{Im}}\nolimits z > 0; \quad A_{{\left (\begin{array}{c} 2,4\\ 3,4 \end{array}\right )}}, A_{{\left (\begin{array}{c} 3,4\\ 2,4 \end{array}\right )}}:~ \mathop {\textrm{Im}}\nolimits z < 0. \end{aligned}$$

Using these results, one can write a fundamental set of analytic eigenfunctions in either half plane as

$$\begin{aligned} \chi ^\pm (x,t,z) = \Xi ^\pm (x,t,z) \Delta ^\pm (z), \end{aligned}$$
(4.57)

where

$$\begin{aligned} \Delta ^+= & {} \mathop {\textrm{diag}}\nolimits (1,A_{[1]},A_{[1,2]},A_{[1,2,3]}),\\ \Delta ^-= & {} \mathop {\textrm{diag}}\nolimits (A_{[2,3,4]},A_{[3,4]},A_{[4]},1). \end{aligned}$$

Note \(\phi _{1} = \chi ^+_1\) and \(\psi _{4} = \chi ^+_4\) are analytic in the UHP, while \(\psi _{1} = \chi ^-_1\) and \(\phi _{4} = \chi ^-_4\) are in the LHP.

Symmetries. The richness of the 3-component VNLS equation compared to the Manakov system comes from the discrete spectrum and symmetries, as we discuss next. Like in all previous cases, the Lax pair admits two involutions: \(z\mapsto z^*\), mapping the UHP into the LHP and viceversa, and \(z\mapsto q_o^2/z\) mapping the exterior of the circle \(C_o\) into the interior and viceversa. The behavior of the analytic eigenfunctions under these symmetries is different, however, and is obtained by first noting that, for \(z\in {\mathbb {R}}\),

$$\begin{aligned} \Phi (x,t,z)= \,& {} J[\Phi ^\dagger (x,t,z)]^{-1}C = \Phi (x,t,q_o^2/z)\Pi _-, \quad \\ \Psi (x,t,z)= \,& {} J[\Psi ^\dagger (x,t,z)]^{-1}C = \Psi (x,t,q_o^2/z)\Pi _+ \end{aligned}$$
(4.58)

with \(\Pi _\pm (z) = \mathop {\textrm{diag}}\nolimits (0,1,1,0) - iJQ_\pm ^*/z\) and \(C(z) = \mathop {\textrm{diag}}\nolimits (\gamma ,-1,-1,-\gamma )\). One then expresses the non-analytic Jost eigenfunctions in Eqs. (4.58) in terms of the columns of \(\chi ^\pm\) via Eqs. (4.55) and (4.57), and uses the Schwarz reflection principle to lift the resulting relations off the real axis. In particular, wherever the eigenfunctions are analytic, one has

$$\begin{aligned} \psi _{1}^*(x,t,z^*)= & {} - \frac{\textrm{e}^{-2 i \theta _2}}{A_{[1,2]}A_{[1,2,3]}} \, J L[\chi _2^+,\chi _3^+,\psi _{4}], \end{aligned}$$
(4.59a)
$$\begin{aligned} \psi _{1}(x,t,q_o^2/z)= \,& {} (iz/q_\pm )\,\psi _{4}, \end{aligned}$$
(4.59b)
$$\begin{aligned}{ [\chi _2^+(x,t,z*)]^*}= - \frac{\textrm{e}^{- 2 i \theta _2}}{A_{[4]}\gamma } \, J L [\psi _{1},\chi _3^-,\phi _{4}], \end{aligned}$$
(4.59c)
$$\begin{aligned} \chi _2^+(x,t,q_o^2/z)= & {} \frac{\textrm{e}^{i \Delta \theta }}{A_{[3,4]}} \left ( A_{[4]} \chi _2^- + A_{{\left (\begin{array}{c} 3,4\\ 2,4 \end{array}\right )}} \chi _3^- \right ), \end{aligned}$$
(4.59d)

and similarly for the other analytic eigenfunctions. Here,

$$\begin{aligned} L[{\textbf{u}},{\textbf{v}},{\textbf{w}}] = \det \begin{pmatrix} u_1&{}u_2&{}u_3&{}u_4 \\ v_1&{}v_2&{}v_3&{}v_4 \\ w_1&{}w_2&{}w_3&{}w_4 \\ {\textbf{e}}_1&{}{\textbf{e}}_2&{}{\textbf{e}}_3&{}{\textbf{e}}_4 \end{pmatrix} \end{aligned}$$

is the multilinear and totally antisymmetric operator which generalizes the familiar cross-product to four dimensions, and \(\Delta \theta = \theta _+ - \theta _-\). Corresponding symmetries exist for the scattering data. In particular, wherever the minors are analytic, one has

$$\begin{aligned} A_{[1]}(z)= & {} A_{[2,3,4]}^*(z^*) = \textrm{e}^{i \Delta \theta } A_{[4]}(q_o^2/z), \end{aligned}$$
(4.60a)
$$\begin{aligned} A_{[1,2]}(z)= \,& {} A_{[3,4]}^*(z^*), \quad A_{[1,2,3]}(z) = A_{[4]}^*(z^*), \end{aligned}$$
(4.60b)
$$\begin{aligned} A_{{\left (\begin{array}{c} 1,2\\ 1,3 \end{array}\right )}}(z)= \,& {} - A_{{\left (\begin{array}{c} 3,4\\ 2,4 \end{array}\right )}}^*(z^*) = \textrm{e}^{i \Delta \theta } A_{{\left (\begin{array}{c} 2,4\\ 3,4 \end{array}\right )}}(q_o^2/z), \end{aligned}$$
(4.60c)
$$\begin{aligned} A_{{\left (\begin{array}{c} 1,3\\ 1,2 \end{array}\right )}}(z)= \,& {} - A_{{\left (\begin{array}{c} 2,4\\ 3,4 \end{array}\right )}}^*(z^*), \end{aligned}$$
(4.60d)

as well as

$$\begin{aligned} \textrm{e}^{i\Delta \theta } A_{[1,2]}(z) A_{[1,2]}^*(q_o^2/z^*) = A_{[1]}(z) A_{[1,2,3]}(z) + A_{{\big (\begin{array}{c} 1,2\\ 1,3 \end{array}\big )}}(z) A_{{\big (\begin{array}{c} 1,3\\ 1,2 \end{array}\big )}}(z). \end{aligned}$$
(4.60e)

These symmetries play a crucial role in the characterization of the discrete spectrum. Indeed, it is the presence of \(L[\cdot ]\) and the non-principal analytic minors in Eqs. (4.59), (4.60e) that make the discrete spectrum and the corresponding soliton solutions of the 3-component case much richer than those in the Manakov system.

Discrete spectrum. The discrete spectrum is comprised of the values of \(z\in {\mathbb {C}}\) for which the columns of \(\chi ^\pm\) are linearly dependent. Equations (4.55)–(4.57)) yield

$$\begin{aligned} \det \chi ^+= A_{[1]} A_{[1,2]} A_{[1,2,3]} \gamma \textrm{e}^{2 i \theta _2}, \end{aligned}$$
(4.61a)
$$\begin{aligned} \det \chi ^-= A_{[4]} A_{[3,4]} A_{[2,3,4]} \gamma \textrm{e}^{2 i \theta _2}. \end{aligned}$$
(4.61b)

The scattering problem is self-adjoint, so bound states can only occur for \(k\in {\mathbb {R}}\), i.e., \(|z|=q_o\), and these give rise to dark solitons. On the other hand, the analytic principal minors in Eq. (4.61) can have zeros for \(|z|\ne q_o\). As in the Manakov system, such zeros yield dark-bright solitons.

The symmetries imply that discrete eigenvalues appear in symmetric quartets \(Z_n=\{z_n,z_n^*,\hat{z}_n,\hat{z}_n^*\}\), where \(\hat{z} = q_o^2/z^*\). When the non-principal analytic minors are identically zero, each quartet yields a dark-bright soliton in which the bright component is aligned exclusively with either the first or the second component of \({\textbf{q}}(x,t)\), while the dark part is along the third component of \({\textbf{q}}(x,t)\). On the other hand, some of the extra analytic minors might be nonzero at one of the discrete eigenvalues. Specifically, for each quartet \(Z_n\) we consider the following configuration of simple zeros for the analytic minors: \(A_{[1]}(z_n)=0\) and \(A_{[1,2]}(z_n)A_{[1,2,3]}(z_n)\ne 0\), as well as

$$\begin{aligned} A_{{\left (\begin{array}{c} 1,2\\ 1,3 \end{array}\right )}}(z_n)=0,\quad A_{{\left (\begin{array}{c} 1,3\\ 1,2 \end{array}\right )}}(z_n)\ne 0. \end{aligned}$$
(4.62)

The eigenfunctions can then be shown to be related as follows:

$$\begin{aligned} \phi _{1}(x,t,z_n)= \,& {} b_n \chi _2^+(x,t,z_n), \end{aligned}$$
(4.63a)
$$\begin{aligned} \chi _2^-(x,t,z_n^*)= \,& {} d_n\psi _{1}(x,t,z_n^*), \end{aligned}$$
(4.63b)
$$\begin{aligned} \chi _3^+(x,t,\hat{z}_n)= \,& {} e_n \psi _{4}(x,t,\hat{z}_n), \end{aligned}$$
(4.63c)
$$\begin{aligned} \phi _{4}(x,t,\hat{z}_n^*)= \,& {} f_n \chi _2^-(x,t,{\hat{z}}_n^*) + g_n{\tilde{\chi }}_3^-(x,t,{\hat{z}}_n^*), \end{aligned}$$
(4.63d)

where \(d_n\), \(e_n\) and \(f_n\) are determined in terms of \(b_n\) via the symmetries, \({\tilde{\chi }}_3^-(x,t,z) = \chi _3^-(x,t,z)/A_{[3,4]}(z)\) is finite at \({\hat{z}}_o^*\), and \(g_n\) is proportional to the nonzero analytic non-principal minor.

Inverse problem. The inverse problem is formulated in terms of a Riemann–Hilbert problem (RHP) for the sectionally meromorphic matrix \(\mu (x,t,z) = \mu ^\pm\) for \(\mathop {\textrm{Im}}\nolimits z\gtrless 0\), with

$$\begin{aligned} \mu ^+= & {} \bigg (\frac{\phi _{-,1}}{A_{[1]}},\frac{\chi _2^+}{A_{[1,2]}} - \frac{A_{{\left (\begin{array}{c} 1,3\\ 1,2 \end{array}\right )}}\chi _3^+}{A_{[1,2]}A_{[1,2,3]}}, \frac{\chi _3^+}{A_{[1,2,3]}},\phi _{+,4}\bigg )\textrm{e}^{-i\Theta }, \\ \mu ^-= & {} \bigg (\phi _{+,1},\frac{\chi _2^-}{A_{[2,3,4]}}, \frac{\chi _3^-}{A_{[3,4]}} - \frac{A_{{\big (\begin{array}{c} 2,4\\ 3,4 \end{array}\big )}}\chi _2^-}{A_{[3,4]}A_{[2,3,4]}},\frac{\phi _{-,4}}{A_{[4]}}\bigg )\,\textrm{e}^{-i\Theta }. \end{aligned}$$

Indeed, manipulating the scattering relation (4.56), one obtains the jump condition

$$\begin{aligned} \mu ^+ = \mu ^-(I - \textrm{e}^{-iK\Theta }L\textrm{e}^{iK\Theta }), \quad z \in {\mathbb {R}}, \end{aligned}$$
(4.64)

where \(K = \mathop {\textrm{diag}}\nolimits (-1,1,1,-1)\) and L(z) is explicitly determined in terms of the reflection coefficients of the problem: \(\rho _1(z) = {a_{21}}/{a_{11}}\), \(\rho _2(z) = {a_{31}}/{a_{11}}\) and \(\rho _3(z) = {a_{41}}/{a_{11}}\).

Since \(\mu ^\pm (x,t,z)\rightarrow I\) as \(z\rightarrow \infty\), one can use Cauchy projectors to reduce the RHP to a system of linear integral equations. In addition, if a nontrivial discrete spectrum is present, as usual one must supplement the system with appropriate algebraic equations, obtained by computing the residues of \(M^\pm (x,t,z)\) at the discrete eigenvalues. The residues of \(\mu ^\pm (x,t,z)\) at each point of \(Z_n\) can then be computed via Eqs. (4.63) and expressed in terms of a single, arbitrary complex vector norming constant \({\textbf{c}}_n=(c_{1,n},c_{2,n})^T\).

Finally, computing the asymptotic behavior of the solution of the RHP as \(z\rightarrow \infty\) and comparing with the asymptotic behavior obtained from the direct problem allows one to write down a reconstruction formula for the solution of the VNLS equation (4.51):

$$\begin{aligned} q_j(x,t) = -i\lim _{z\rightarrow \infty }z\, \mu _{j+1,1}(x,t,z), \qquad j=1,2,3. \end{aligned}$$
(4.65)

As usual, in the reflectionless case [\(\rho _j(z)\equiv 0\)] the RHP reduces to a linear algebraic system and one obtains the pure soliton solutions.

Solitons. As mentioned above, the 3-component VNLS admits DBB soliton solutions, with two bright and one dark components. When only one quartet of discrete eigenvalues is present, letting \(z_o=\nu _o+i\eta _o = |z_o|\textrm{e}^{i\alpha _o}\), with \(|z_o|<q_o\) and a vector norming constant \({\textbf{c}}_o\), the DBB soliton has the form

$$\begin{aligned} q_j(x,t)= & {} -ip_{j,o} \,w_o\sin \alpha _o\,\textrm{e}^{i(\nu _o x-(\nu _o^2-\eta _o^2)t)}\mathop {\textrm{sech}}\nolimits [\eta _o(x-2v_ot - x_o)], \quad j=1,2, \end{aligned}$$
(4.66a)
$$\begin{aligned} q_3(x,t)= \,& {} q_+ \textrm{e}^{i\alpha _o} ( \cos \alpha _o -i\sin \alpha _o \tanh [\eta _o(x-2v_ot - x_o)] ), \end{aligned}$$
(4.66b)

where \(w_o = \smash {\sqrt{q_o^2-|z_o|^2}}\), and the unit-norm polarization vector for the bright components is simply \({\textbf{p}}_o = (p_{1,o},p_{2,o})^T = {\textbf{c}}_o/\Vert {\textbf{c}}_o\Vert\).

Unlike DB soliton solutions of the Manakov system, DBB solitons exhibit nontrivial polarization interactions (see Fig. 6), i.e., energy (and phase) exchanges between the bright components of the interacting solitons, similarly to those arising in the focusing Manakov system with decaying potentials. Importantly, in recent experimental works [29], DBB solitons were observed in a BEC using a method based on local spin rotations which simultaneously imprint suitable phase and density distributions. This subsequently enabled the observation of the striking collisional properties of the emerging multi-component solitons [125], and the results showed a remarkable quantitative agreement with our analytical predictions in [38, 155].

Fig. 6
figure 6

A 2-soliton solution of the defocusing 3-component VNLS equation exhibiting a polarization shift, obtained for \(q_o=1\) with \(z_1=i/2\) (stationary soliton), \(z_2=(1+i)/4\) (moving soliton) and norming vectors \({\textbf{c}}_1 = (1,0)^T\) and \({\textbf{c}}_2=(1,1+i/2)^T\). Note how the bright component of soliton 1 is initially aligned exclusively with \(q_1\), but acquires a component along \(q_2\) as a result of the interaction. From [38]

5 Applications, Open Problems and Future Directions

A robust IST for the focusing NLS equation with NZBC. One of the important recent developments of the IST with NZBC was proposed by Bilman and Miller in the context of the scalar, focusing NLS equation with the purpose of dealing with arbitrary order poles and potentially severe spectral singularities [30]. In the previous sections, we have only considered simple discrete eigenvalues (i.e., poles of order 1 in the RHP), and we have routinely assumed initial conditions for which the scattering problem has no spectral singularities. The first assumption was only for convenience of presentation, and, as discussed above, higher order poles in the RHP can and have been considered in various works for both scalar and coupled NLS systems, both for decaying potentials and with NZBC. On the other hand, including spectral singularities in the IST is problematic, and, when NZBC are considered, there is currently no explicit constraint on the initial condition that allows one to exclude their presence in the spectrum. Furthermore, the Peregrine solution (2.46) described in Sect. 2 corresponds to a particular type of spectral singularity at the edges \(\pm i q_o\) of the continuous spectrum \(\Sigma\). Although it can be obtained from other solutions as a limit, it had not been obtained directly as the solution of a Cauchy problem via IST. The work [30] proposes a rigorous way to generalize the IST to capture rogue waves—in this context, rational (or semi-rational) solutions like (2.46)—and it allows one to also handle other types of spectral singularities, as well as poles of arbitrary order in unified way.

The IST in [30] is formulated in the k-plane, and not in terms of the uniformization variable z, and the background is chosen as \(q_\pm =1\), but the general ideas can be adapted to arbitrary \(q_\pm\). The main problem with the traditional IST formulation concerns the nature of the boundary values taken by the matrix \(\mu (x,t,\cdot )\) in the RHP on the continuous spectrum \(\Sigma _k\). The direct problem shows that these boundary values are continuous functions of k except at the points \(\pm iq_o\), where \(\mu\) is singular. When the initial condition coincides with the background potential, with the normalization used in [30] the eigenfunctions in the RHP have \((k\mp q_o)^{-1/4}\) singularities at the branch points, and a growth condition \(O(k\mp q_o)^{-1/4}\) at the branch points is expected, and in fact it is required to ensure uniqueness of solution for the inverse problem.

As discussed in Sect. 2, spectral singularities are zeros of the analytic scattering coefficients a(k) and \({\bar{a}}(k)\) that fall on the continuous spectrum \(\Sigma _k\). At these points, the sectionally meromorphic matrix \(\mu (x,t,\cdot )\) in the RHP fails to have a well-defined nontangential limit from one side or the other of the contour. Spectral singularities can occur in the interior of \(\Sigma _k\), or at its edges, i.e., the branch points \(\pm i q_o\). In the rapidly decaying case, it was shown by Zhou [195] that there exist Schwartz class initial conditions for which spectral singularities, possibly even an infinite number of them, are accumulation points of zeros of a(k) and \({\bar{a}}(k)\). A priori, a similar situation could occur in the case of NZBC. On the other hand, dealing with possible spectral singularities at the branch points is unavoidable with a nonzero background if one wants to include the Peregrine solution in the class of admissible potentials. Moreover, a careful analysis shows that the Peregrine solution is indistinguishable from the background at the level of the scattering data, and the associated eigenfunctions only differ by the rate at which they blow up as k approaches the branch points.

The key idea of the “robust” IST developed in [30] is to formulate the RHP for the eigenfunctions on a contour that completely avoids the branch cut \([-iq_o,iq_o]\). Specifically, in a suitable, proper subset \(\Sigma ^\prime _k={\mathbb {R}}\setminus [-r,r]\) of the continuous spectrum \(\Sigma _k={\mathbb {R}}\cup [-iq_o,iq_o]\), the reflection coefficient is defined as in the traditional IST, and an associated jump matrix is introduced. The radius \(r>q_o\) that bounds \(\Sigma _k^\prime\) away from the origin depends on the initial condition, and a (non-sharp) estimate for the lower bound of r is explicitly provided in [30]. For \(|k|=r\), the Jost solutions and the scattering coefficients a(k), \({\bar{a}}(k)\) can be computed in the standard way, and an appropriate jump matrix is defined on the circle of radius \(r>q_o\). Essentially, the function \(\mu (x,t,k)\) in the robust IST coincides with the one in the traditional IST for \(|k|>r\), and outside this disk \(\mu (x,t,k)\) is analytic except for \(k\in {\mathbb {R}}\setminus [-r,r]\), where it has a jump which, again, coincides with the standard jump. On the other hand, \(\mu (x,t,k)\) has a different definition inside the disk \(|k|<r\), and it is analytic except on the branch cut. \(\mu (x,t,k)\) is continuous on the upper and lower semicircles of radius r, and the boundary values from each side are related by a suitable jump condition. Note that the ensuing RHP does not require to build in poles from the discrete eigenvalues (regardless of their number and order), or to specify the growth rate of the eigenfunctions, as these pieces of information are encoded in the jump across the circle \(|k|=r\). Also, one does not need to make assumptions on the type or severity of spectral singularities for \(|k|<r\). This method can be applied to other scalar integrable PDEs, and possibly extended to matrix and vector NLS systems.

Long-time asymptotics. One of the notable applications of the IST is related to the study of the long-time asymptotic behavior of solutions of integrable equations. In the rapidly decaying case, the long-time behavior of solutions of the defocusing NLS equation (and of the focusing NLS equation without solitons) was studied early on, the absence of a discrete spectrum being a significant simplification. The solution q(xt) is found to decay like \(t^{-1/2}\) as \(t\rightarrow \infty\), the leading term (without any error estimate) obtained by Zakharov and Manakov [190] and independently by Segur and Ablowitz [163, 164] already in 1976 from the study of the GLM integral equations in terms of which the inverse problem was originally formulated. As shown later on using the RHP and the steepest descent technique, estimates on the size of the error depend on smoothness and decay assumptions on the initial condition. Specifically, in [60, 62, 63] it was proved that if the initial condition has a high degree of smoothness and decay, then the error is \(O(t^{-1}\log t)\). This estimate was subsequently improved, under the much weaker assumption that the initial condition belong to the weighted Sobolev space \(H^{1,1}({\mathbb {R}})\) [64]. In [68] the asymptotic analysis of the RHP based on \({\bar{\partial }}\)-problems, rather than the asymptotic analysis of singular integrals on contours, was employed to obtain sharp asymptotics for the solution of the defocusing NLS, leading to an error estimate \(O(t^{-3/4})\) with the same regularity assumption that the initial condition be in \(H^{1,1}({\mathbb {R}})\). The \({\bar{\partial }}\)-steepest descent method was then utilized in [57] to obtain the leading order approximation to the solution of the defocusing NLS with NZBC for large times in the solitonic region of space-time, with error bound which decays as \(t\rightarrow \infty\) for initial data in a suitable class. The problem was recently revisited in [184], where a different asymptotic expansion is obtained in the solitonless region, with an \(O(t^{-1/2})\) contribution from the continuous spectrum, and an error term \(O(t^{-3/4})\).

Following the earlier results obtained in the absence of a discrete spectrum, the long-time asymptotic analysis of the initial-value problem for the focusing NLS with rapidly decaying data was fully developed in [45] using the \({\bar{\partial }}\)-generalization of the nonlinear steepest descent method. For an initial datum in \(H^{1,1}({\mathbb {R}})\) with a simple discrete spectrum and no spectral singularities, it is shown in [45] that in any fixed space-time cone the leading order term of the solution as \(t\rightarrow \infty\) is a multi-soliton solution whose parameters are modulated by soliton-soliton and soliton-radiation interactions as one moves through the cone, with an error estimate \(O(t^{-3/4})\). Special cases of initial conditions producing a single spectral singularity, and an infinite number of solitons were considered in [107, 108]. In the case of the focusing NLS with NZBC, some earlier results were obtained in [46] for the long-time dynamics of solutions with step-like initial conditions. Notably, in [40] the IST-based analysis of the long-time asymptotics of the focusing NLS with NZBC was used to characterize the nonlinear stage of modulational instability. Asymptotically in time, the spatial domain splits into three regions: two far-field regions, in which the solution is approximated by its initial value up to a phase shift, on either side of a central region in which the solution exhibits oscillatory behavior described by slow modulations of the periodic traveling wave solutions of the focusing NLS. The full analysis of the associated RHP in the long-time limit was developed in [41] under the assumption that no discrete spectrum is present, and generalized in [39] to initial conditions that allow for the presence of discrete spectrum, which unveiled different possible interaction regimes of soliton-radiation interaction (soliton transmission, soliton trapping, and a mixed regime in which the soliton transmission or trapping is accompanied by the formation of a nondispersive soliton-generated wake) depending of the location of the discrete eigenvalue in the spectral plane.

Some recent results on the focusing Manakov system with rapidly decaying potentials can be found in [83, 183], but the investigation of the long-time asymptotics of the initial-value problem for the focusing or defocusing Manakov systems with NZBC is a completely open problem.

Applications of the IST. The last few years have seen a blossoming of exciting new applications of the IST to several experimental settings. In [157], the IST is used to tackle the problem of identifying exact solutions of the focusing NLS on a nontrivial background that exhibit rogue wave events of statistical significance. A new approach to the classification of rogue waves is proposed, in which appropriate, locally coherent structures are isolated from a globally incoherent wave train, and then analyzed by means of a numerical IST procedure. In [82], a theoretical model of the asymptotic stage of the noise-induced modulational instability (MI) is developed within the IST formalism. Specifically, the model relies on ensembles of N-soliton bound states of the focusing NLS with a nontrivial background, with a semiclassical distribution for the discrete eigenvalues, and random phases for the norming constants. The numerical study revealed a remarkable agreement between the spectral (Fourier) and statistical properties of the long-time evolution of the MI, and those of the multisoliton, random-phase bound states, and the approach can be likely extended to a large class of nonlinear integrable turbulence problems.

Other recent applications of the IST involved analysing focusing and defocusing NLS systems with various box-type (piecewise constant) initial configurations. For instance, in [81] the IST was used to investigate the solutions of the focusing NLS equation with an initial condition corresponding to a spatially extended box-shaped wave field, with a plane wave (a condensate) in the middle, and zero at the edges of the box. The scattering data consist of the continuous spectrum, and the soliton eigenvalues and associated norming constants, with the number of solitons, N, proportional to the box width. The continuous spectrum is then removed from the scattering data, and analytic expressions for the corrections to the norming constants that result from the removal of the continuous spectrum are obtained. The corrected soliton parameters produce a spatially symmetric N-soliton solution, which converges asymptotically as \(N\rightarrow \infty\) to the condensate, thus providing an effective solitonic model for the condensate. In [159] the exact analytical tools provided by the IST for the defocusing NLS with NZBC were exploited to achieve an on-demand generation of dark soliton arrays. The relevant physical model model is that of a 1-dimensional, harmonically trapped scalar BEC composed of repulsively interacting atoms, and the response of such a system to box-type initial configurations was analyzed. The physical parameters of each of the enucleated solitons can be obtained by computing the discrete eigenvalues of the scattering problem. As expected, the size of the box directly affects the number, size and velocity of the solitons, while the initial phase determines the parity of the solutions. The analytical predictions for the soliton amplitudes and velocities showed remarkable agreement with direct numerical simulations, both in the presence of and without harmonic confinement. These results were generalized in [160] to the defocusing Manakov system with box-type initial conditions, leading to the controlled creation of dark-bright soliton trains in repulsive two-component BECs.

We should also mention that the numerical implementation of the IST has attracted a good deal of attention recently. The presence of oscillatory dispersive tail, and the fact that when solitons are present, they may never escape the dispersion, limit the efficiency of traditional numerical methods at capturing the solution of dispersive PDEs for large t, as the computational cost of time-stepping methods grows rapidly in time. These limitations make the numerical implementation of the IST very attractive, since the computational cost to compute the solution for any give xt is seen to be bounded. In [170], Trogdon and Olver solved the initial-value problem for the focusing and defocusing NLS equations in the rapidly decaying case numerically by implementing the IST. Their numerical IST scheme has two major components: the first is the use of a Chebyshev collocation method for solving RHPs [140, 171], and the second is the deformation of contours in the spirit of the method of nonlinear steepest descent. The computation of the scattering data and of the NLS solution are shown to both be spectrally convergent sufficiently smooth and rapidly decaying initial conditions.

In certain niches engineering communities, the IST as a signal processing technique is known as the non-linear Fourier transform (NFT), and several numerical methods have also been implemented to detect and characterize solitons within wave profiles [69, 180, 181]. In [147], for instance, a new algorithm aimed at computing the phase shift of solitons in processes governed by the Korteweg-de Vries (KdV) equation was proposed. In numerical examples, the new algorithm is found to perform reliably even in cases where existing algorithms break down. In [174] an IST-based approach to digital signal processing in optical communications was explored, showing through numerical modeling how an initial wave profile can be recovered from the scattering data of the received optical signal, without direct backward propagation. In a similar context, [80] exploited the mathematical properties of the NLS equation and its multisoliton solutions to mitigate nonlinear signal distortions in optical fibers. Choosing N discrete eigenvalues with the same imaginary parts, a Soliton Orthogonal Frequency Division Multiplexing (SOFDM) method is proposed as an alternative to the traditional OFDM scheme, which likewise allows for the use of efficient fast Fourier transform algorithms to recover the data. In the proposed approach, signal modulation is applied to the kernel of the GLM equations that define the inverse problem for the eigenfunctions. The approach can also be used to control signal parameters in the presence of a continuous spectrum.

Future directions. In addition to the various open problems highlighted in each of the previous sections, future directions in the field will likely involve attempts to generalize the IST to other types of physically relevant boundary conditions, such as periodic boundary conditions, or more general x-dependent boundary conditions. Even more ambitiously, we think the frontier is currently represented by a generalization of the IST to slowly decaying potentials, and in turn to potentials with slow decay to certain constant boundary conditions, which are related to the emergent theory of soliton gases.

Generally speaking, the existence of a Lax pair does not guarantee per se the solution of the initial-value problem, and the question as to whether an initial-value problem for an integrable PDE could be solved by a suitable generalization of the IST outside the classes of rapidly decaying or periodic potentials was raised by various distinguished researchers over the years (see, for instance, [59, 122, 132]). In reference to the KdV equation, in [70] it was recently pointed out that the IST is not yet developed to a level which would satisfy a pragmatic physicist, in that little or nothing is currently known when the initial datum is neither (sufficiently rapidly) decaying at infinity nor periodic, and the same is true for potentials that do not decay sufficiently rapidly to a constant background. The question of extending the IST to initial data that are simply bounded or slowly decaying is obviously of great practical importance, as is the case of a random statistically uniform field, which is associated to integrable turbulence [189]. In this regard, a new approach to the IST based upon Hankel operators was proposed in [162], where the initial-value problem for the KdV equation was solved for an arbitrary initial datum bounded from below, decaying sufficiently rapidly as \(x\rightarrow +\infty\), but otherwise unrestricted and, in particular, without imposing any boundary condition as \(x\rightarrow -\infty\). The work was motivated by the questions raised in [70], and it offers at least a proof of principle that an affirmative answer to the question is possible, although the result could rely on the unidirectional nature of the KdV equation, and hence not be generalizable to other integrable PDEs. These questions are currently completely open, and naturally lead to questions about the spectrum of a scattering operator when the potential is simply a bounded function. In [70] the reverse question was addressed, namely, assuming the scattering problem has two allowed bands and a suitable RHP is defined by a pair of positive Hölder continuous functions \(r_1(k),r_2(k)\) on the allowed bands, the associated potentials (called primitive potentials) are shown to be bounded, but are neither periodic nor quasi-periodic, and do not decay as \(x\rightarrow \pm \infty\). In the special case in which one of the two functions—say, \(r_2(k)\)—is identically zero, the RHP that characterizes the primitive potentials was subsequently shown to arise as the limit as \(N\rightarrow \infty\) of a gas of N-solitons [96]. The primitive potentials approach can be generalized to a wide class of spectral problems, allowing to construct non-decaying, non-periodic solutions to the various associated nonlinear integrable systems, such as the focusing and defocusing NLS equations. The soliton condensates and breather gas of the focusing NLS equation studied in [71, 72, 169] are likely related to the notion of primitive potentials, as in the KdV case.

With the developments of the last fifteen years, the IST, which burst into the scene in the early 70’s with ground-breaking research, has seen a marked revival, a flourishing of papers and new, exciting applications which, in our opinion, bode well for the future.