Semiclassical limit of the focusing NLS: Whitham equations and the Riemann-Hilbert Problem approach

In this paper we consider the focusing Nonlinear Schrödinger (fNLS) equation

(1.1)

$$i\varepsilon\partial_{t}\psi+\frac{1}{2}\varepsilon^{2}\partial_{x}^{2}\psi+|\psi|^{2}\psi=0,$$

where $x\in{\mathbb{R}}$ and $t\geq 0$ are space-time variable and $\varepsilon>0$. In the semiclassical (small dispersion) limit we take $\varepsilon\rightarrow 0^{+}$. The fNLS is a basic model for self-focusing and self-modulation, for example, it governs nonlinear transmission of light in optical fibers; it can also be derived as a modulation equation for a broad class of nonlinear systems. It was first integrated (with $\varepsilon=1$) by Zakharov and Shabat [38], who produced a Lax pair for it and used the inverse scattering procedure to describe the time evolution of general decaying potentials ($\lim_{|x|\to\infty}\psi(x,0,\varepsilon)=0$) in terms of the scattering data, that is, radiation and solitons.

Integrability of the fNLS enables the existence of the finite-gap solutions, the non-decaying quasi-periodic potentials with many remarkable properties [1]. Slow modulations of these potentials are described by the Whitham modulation equations [37], a system of quasilinear equations derived by the averaging procedure applied to fNLS. The remarkable feature of the fNLS evolution is that the modulated finite-gap potentials arise as asymptotic solutions in the initial value problems for the fNLS with non-oscillating initial conditions (e.g. Gaussian or sech initial data for (1.1)). This has been understood in the framework of the Riemann-Hilbert problem (RHP) approach to the inverse scattering yielding both the finite-gap potentials and their modulations within the rigorous asymptotic description of the semiclassical fNLS evolution. However, the connections between the Whitham modulation theory and the RHP approach have not been properly traced.

The main results of the present paper:

• •

  We link the main objects of the Whitham theory (Riemann invariants, their characteristic velocities, fundamental differentials $dp,dq$) with the corresponding objects coming from the RHP approach, which are expressed in terms of the $g$-function from the RHP approach and the corresponding hyperelliptic Riemann surface $\mathcal{R}$;

• •

  We derive several equivalent forms of the modulation equations (by which we mean systems of transcendental equations for the branchpoints of $\mathcal{R}$), and use them to show that various forms of the generalized hodograph solutions of the Whitham equations are satisfied;

• •

  we use the $g$-function to establish the role of the potential for the differentials $dp,dq$ in determining the breaking (phase transition) curves and double point velocities for the fNLS.

The first (spatial) equation of the Lax pair for the fNLS is known as Zakharov - Shabat (ZS) system

(1.2)

$$i\varepsilon\frac{d}{dx}W=\begin{pmatrix}z&\psi\\ \bar{\psi}&-z\end{pmatrix}W~,$$

where $z$ is a spectral parameter and $W$ is a $2$ by $2$ matrix-function.

The scattering data, corresponding to the initial data $\psi(x,0,\varepsilon)$, consist of the reflection coefficient $r_{0}(z,\varepsilon)$, as well as of the points of the discrete spectrum, if any, together with their norming constants. The time evolution of the scattering data is simple ([38]). The corresponding evolution of a given potential is obtained through the inverse scattering transformation (IST) of the evolving scattering data. The inverse scattering transform for the fNLS equation (1.1) at the point $x,t$ can be cast as the following matrix Riemann-Hilbert Problem (RHP) in the spectral $z$-plane: find a $2\times 2$ matrix-valued function $m(z)=m(z;x,t,\varepsilon)$, which depends on the asymptotic parameter $\varepsilon$ and the external parameters $x,t$, such that: i) $m(z)$ is analytic in ${\mathbb{C}}\backslash\Gamma$, where the contour $\Gamma={\mathbb{R}}$ has the natural orientation; ii)

(1.3)

$$m_{+}=m_{-}\begin{pmatrix}1+r\bar{r}&\bar{r}\cr r&1\cr\end{pmatrix}=m_{-}V~$$

on the contour $\Gamma$, where $r(z,\varepsilon)=r_{0}(z,\varepsilon)\exp[{2i\over\varepsilon}(2z^{2}t+zx)]$ and $m_{\pm}(z)=\lim_{\delta\rightarrow 0}m(z\pm i\delta)$ with $\delta>0$ and $z\in{\mathbb{R}}$; iii) $\lim_{z\rightarrow\infty}m(z)={\bf 1},$ where ${\bf 1}$ denotes the identity matrix. In the presence of solitons the contour $\Gamma$ contains additional small circles around the eigenvalues with the corresponding jump-matrices (see, for example, [30] or [19]).

A complementary approach to the analysis of nonlinear dispersive equations is provided by the Whitham modulation theory which produces a system of quasilinear partial differential equations governing slow modulations of periodic or quasi-periodic waves. A prominent area where the Whitham method proved extremely useful is the theory of dispersive shock waves (DSWs) [16], [11]

We consider an $N$-phase quasi-periodic solution of the fNLS equation (1.1),

(1.4)

$$\psi=\psi_{N}(\varepsilon^{-1}\eta_{1},\varepsilon^{-1}\eta_{2},\dots,\varepsilon^{-1}\eta_{N};{\boldsymbol{\alpha}},\bar{\boldsymbol{\alpha}}),\quad\hbox{where}\quad\eta_{j}=k_{j}x+\omega_{j}t+\eta^{0}_{j}.$$

Here $\boldsymbol{\alpha}\equiv\{\alpha_{j}\}_{j=0}^{N}\in\mathbb{C}^{N}$ and c.c. are the branchpoints of the hyperelliptic Riemann surface of genus $N$ on which (1.4) lives [1], the (normalized by $\varepsilon$) wavenumbers $k_{j}$ and the frequencies $\omega_{j}$ are defined in terms of the branchpoints, and $\eta_{j}^{0}$ are arbitrary initial phases. Now, if the parameters $\alpha_{j}$, $\bar{\alpha}_{j}$ are allowed to vary slowly in space and time they must satisfy the Whitham modulation equations,

(1.5)

$$(\alpha_{j})_{t}=V_{j}^{(N)}({\boldsymbol{\alpha}},\boldsymbol{\bar{\alpha}})(\alpha_{j})_{x},\quad(\bar{\alpha}_{j})_{t}=\overline{V}_{j}^{(N)}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})(\bar{\alpha}_{j})_{x},\qquad j=0,\dots,N\,,$$

so that $\alpha_{j}$, $\bar{\alpha}_{j}$ are Riemann invariants. The characteristic velocities $V_{j}^{(N)}$, $\overline{V}_{j}^{(N)}$ are expressed in terms of the Riemann invariants $\alpha_{j}$ through hyperelliptic (for $N\geq 2$) or complete elliptic (for $N=1$) integrals. For $N=0$ the Whitham system (1.5) has the form

(1.6)

$$(\alpha_{0})_{t}=(\tfrac{3}{2}\alpha_{0}+\tfrac{1}{2}\bar{\alpha}_{0})(\alpha_{0})_{x},\qquad(\bar{\alpha}_{0})_{t}=(\tfrac{3}{2}\bar{\alpha}_{0}+\tfrac{1}{2}\alpha_{0})(\bar{\alpha}_{0})_{x}\,,$$

and is equivalent to the dispersionless limit of (1.1)

(1.7)

$$\rho_{t}+(\rho u)_{x}=0,\quad u_{t}+uu_{x}-\rho_{x}=0\,,$$

where the Riemann invariants and characteristic speeds in (1.6) are expressed in terms of the hydrodynamic “density” $\rho=|\psi|^{2}\geq 0$ and “velocity” $u=-i\varepsilon^{-1}(\arg\psi)_{x}$ as

(1.8)

$$\alpha_{0}=-(\frac{u}{2}+i\sqrt{\rho}),\qquad V_{0}^{(0)}=\tfrac{3}{2}\alpha_{0}+\tfrac{1}{2}\bar{\alpha}_{0}=-(u+i\sqrt{\rho})\,.$$

One can see that the characteristics of (1.7) are complex unless $\rho=0$ implying nonlinear modulational instability of the NLS equation (1.1) in the long-wave limit, and hence, ill-posedness of the initial-value problem for (1.6) for all but analytical initial data.

For $N\geq 1$ the characteristic velocities $V_{j}^{(N)}$ in (1.5) are also complex, however, vanishing of the imaginary parts for some $V_{j}$’s is possible. This “partial hyperbolicity” property makes the Whitham systems for $N\geq 1$ radically different compared to the genus zero case and has major implications in terms of stability of some solutions [10].

There are (at least) two ways of looking at the Whitham equations. Originally, they were introduced in [35] as the equations obtained by the averaging of dispersive conservation laws over the family of periodic or quasi-periodic solutions. Later, Whitham put his method in the very general variational principle framework [36]. Another approach leading to the same set of equations is the multiple-scale (nonlinear WKB) expansions method [24]. If the original dispersive equation is IST integrable, the associated Whitham system turns out to be also integrable in the sense which will be explained later on.

Using the finite-gap theory of the KdV equation [26], [21] Flaschka, Forest and McLaughlin [14] showed that the endpoints of the spectral bands of quasi-periodic finite-gap potentials of the quantum-mechanical Schrödinger operator are Riemann invariants of the Whitham modulation system associated with the KdV equation. Analogous result for the fNLS equation was obtained by Forest and Lee [15] and Pavlov [27].

The second context in which the Whitham equations appear is due to Lax, Levermore and Venakides [22], [33], [23] who derived them as the equations governing the zero dispersion limit of the KdV equation. The dispersion parameter $\varepsilon$ determines the typical scale of nonlinear oscillations in the KdV solution so the limit as $\varepsilon\to 0$ exists only in a weak sense. Venakides provided the bridge between the Flaschka-Forest-McLaughlin (wave packet averages) and Lax-Levermore (weak limits) results by developing the higher-order Lax-Levermore theory [34]. Finally, the all-encompassing approach to the small-dispersion KdV was developed by Deift, Venakides and Zhou [5], who introduced the nonlinear steepest descent method for the RHP associated with the (semi-classical) inverse scattering problem.

The Whitham equations naturally arise in the RHP construction as the equations governing the evolution of the spectral branch points $\alpha_{j}$, $\bar{\alpha}_{j}$. To be precise, the RHP theory yields the hodograph solution to the Whitham equations as the combination of the moment conditions and the Boutroux conditions. Although the connection of the Whitham equations with the RHP construction of the semi-classical limit is generally well known it has not been explored to any depth for the fNLS equation. The apparent reason for such an omission is twofold: (i) the existing rigorous analyses of the IVPs for the fNLS equation (see [19], [30] and references therein), performed within the RHP framework, already contain all the information about slow modulations of the solution so there is no need to recover it with the aid of a more restricted Whitham approach; (ii) the Whitham equations (1.5) are elliptic so their application to the problems outside the firmly established facts of the existence and convergence of the relevant solution was suspect. Nevertheless, the success of the application of the Whitham theory to dispersively modified hyperbolic conservation laws, particularly in the DSW theory for the KdV and defocusing NLS equations (see [11] and references therein), provides a strong incentive for the development of a similar theory for the elliptic, focusing case, which is also supported by the extensive numerical evidence that the key features of the small-dispersion “hyperbolic” nonlinear dynamics, such as the co-existence of smooth and rapidly oscillating regions and weak convergence, hold true for at least some cases of the semi-classical fNLS evolution (see, e.g., [3], [4], [10]).

The modulation theory approach to the description of the small dispersion limit of integrable “hyperbolic” equations, such as KdV or defocusing NLS, involves solution of a nonlinear free boundary problem for the associated Whitham modulation system via the so-called matching regularization procedure. This procedure represents an extension of the original Gurevich-Pitaevskii method of the modulation description of a DSW in the KdV equation [13] and prescribes the solution genus increase every time it undergoes a gradient catastrophe. The modulation solutions of different genera are “glued” in a spectial way along the breaking curves which are free boundaries and whose determination is part of the solution (see e.g. [12] for the detailed description of the matching regularization procedure for the KdV equation with monotone initial conditions).

The matching regularization procedure can be relatively easily implemented in the problems involving the fNLS solutions with $N=0$ and $N=1$. This was done in [7], [20], [10] for the fNLS dam-break problem. The modulation solution of [7], [20] was rigorously confirmed in [17] within the RHP analysis of the semi-classical fNLS with the square barrier initial data. The breaks involving $N\geq 2$ as in the first break for the sech potential ([19], [30]), or any higher breaks have not been considered within the Whitham theory with the only exception [10], where the modulation solution beyond the second break was used to predict the generation of rogue waves (note, however, that the second break was considered in [25] in the framework of the RHP approach).

The above discussion strongly suggests that it would be highly desirable to develop a method for solving the fNLS-Whitham equations in problems involving formation of oscillatory regions characterized by the genus $N\geq 2$. It is also clear that the RHP analysis, in particular, the nonlinear steepest descent method with the $g$-function mechanism, can provide valuable clues to the structure of the modulation solutions. Indeed, certain elements of the RHP approach offer an elegant way to circumvent the matching regularization procedure which could be quite awkward when matching modulation solutions with $N\geq 2$. Thus, a closer exploration of the interconnections between the RHP theory of the small-dispersion fNLS limit and the counterpart Whitham modulation theory seems a worthwhile task. Indeed, an appropriate combination of the Whitham theory and some key elements of the RHP which could be termed a “formal RHP analysis”, complemented by careful numerical simulations recently proved very effective for solving problems of immediate physical interest [10].

The paper is organized as follows. In Sections 2, 3 we present the general Whitham theory approach and the RHP approach respectively, as applied to semiclassical fNLS. These two approaches are compared in the case of genus zero in Section 4. In Section 5 we derive an explicit expression for the $g$-function (in the determinantal form) in higher genera regions for the box-type potentials and discuss its connections with the corresponding hyperelliptic Riemann surface $\mathcal{R}$. Solutions of the Whitham equations in terms of $g$-function are discussed in Section 6. The results of this section are based on Theorem 6.1 about three equivalent forms of modulation equations. The results of Sections 5, 6 in the case of analytic potentials were discussed in Section 7. Transitions between the regions of different genera and the characteristic velocities along breaking curves are discussed in Section 8. In particular, we show that, in the square barrier (“box”) potential case, a pair of collided branchpoints on the breaking curve always has real characteristic velocity.

The Whitham equations for the fNLS 1.1 can be represented as a single generating conservation equation [15], [27]

(2.1)

$$\partial_{t}dp=\partial_{x}dq,$$

where $dp(z,\alpha,\bar{\alpha})$ and $dq(z,\alpha,\bar{\alpha})$ are certain meromorphic differentials of the second kind (the quasimomentum and the quasienergy) on the hyperlliptic Riemann surface $\mathcal{R}$ of genus $N$ defined by the radical $R(z)=\sqrt{\prod_{j=0}^{N}(z-\alpha_{j})(z-\bar{\alpha}_{j})}$ where $z$ is the complex spectral parameter. The branch points $\boldsymbol{\alpha}=(\alpha_{0},\dots,\alpha_{N})$ and c.c. are the points of simple spectrum of the periodic Zakharov-Shabat operator (1.2).

The quasimomentum and quasienergy differentials $dp$ and $dq$ are uniquely defined by the following properties [15]:

(a) $dp$ and $dq$ have the poles of order two and three respectively at $\infty^{\pm}$ on $\mathcal{R}$ and no other poles;

(b) the expansions of $dp$ and $dq$ in the local coordinate $z=\zeta^{-1}$ near $\infty^{\pm}$ are

(2.2)

$$\begin{split}dp\sim\pm[-\frac{1}{\zeta^{2}}+\hbox{holomorphic part}]\ \hbox{near}\ \infty^{\pm},\\ dq\sim\pm[-\frac{2}{\zeta^{3}}+\hbox{holomorphic part}]\ \hbox{near}\ \infty^{\pm}\end{split}$$

(c) $dp$ and $dq$ satisfy the normalization conditions

(2.3)

$$\oint\limits_{\hat{\gamma}_{j}}dp=\oint\limits_{\hat{\gamma}_{j}}dq=0\,,\qquad j=1,\dots,N,$$

where $\hat{\gamma}_{j}$ is a clockwise loop around the branchcut connecting $\alpha_{j}$ and $\bar{\alpha}_{j}$ (the $\mathbf{A}$-cycle).

Note that $dp,dq$ are Schwarz symmetrical differentials, so that normalization conditions (2.3) are equivalent to Boutroux normalization conditions for $dp,dq$: all the cycles of $dp,dq$ on $\mathcal{R}$ are real.

The integrals over the $\mathbf{B}$-cycles, canonically conjugated to the $\mathbf{A}$-cycles, give the fundamental wavenumbers $k_{j}$ and frequencies $\omega_{j}$:

(2.4)

$$k_{j}=\oint\limits_{\mathbf{B_{j}}}dp\,,\quad\omega_{j}=\oint\limits_{\mathbf{B_{j}}}dq\,.$$

An alternative useful representation for the wave numbers and frequencies in terms of holomorphic differentials is:

(2.5)

$$k_{j}=-4\pi i\varkappa_{N,j}\,,\qquad\omega_{j}=-4\pi i\left[\tfrac{1}{2}\sum\limits_{k=0}^{N}(\alpha_{k}+\bar{\alpha}_{k})\ \varkappa_{N,j}+\varkappa_{N-1,j}\right],\quad j=1,\dots,N,$$

where $\varkappa_{j,k}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$ are the coefficients of the normalized holomorphic differentials found from the system

(2.6)

$$\sum\limits_{i=1}^{N}\varkappa_{m,k}\oint\limits_{\hat{\gamma}_{k}}\frac{z^{m}}{R(z)}dz=\delta_{mk}\,,\quad m,k=1,\dots,N.$$

Here $\delta_{mk}$ is the Kronecker symbol.

The generating equation (2.1) has several fundamental consequences:

(i) By multiplying (2.1) by $(z-\alpha)^{3/2}$ and letting $z\to\alpha_{j}$ one obtains the diagonal system (1.5) with the characteristic speeds $V_{j}$ given by

(2.7)

$$V_{j}=\frac{dq}{dp}{\bigg{|}}_{z=\alpha_{j}},\quad\overline{V}_{j}=\frac{dq}{dp}{\bigg{|}}_{z=\bar{\alpha}_{j}}.$$

(ii) By expanding (1.5) near $z=\infty$ we obtain an infinite series of averaged local conservation laws of the form $\partial_{t}P_{j}({\boldsymbol{\alpha}},{\boldsymbol{\bar{\alpha}}})=\partial_{x}Q_{j}({\boldsymbol{\alpha}},{\boldsymbol{\bar{\alpha}}})$, where $P_{j}$ are the averaged densities of the NLS conservation laws and $Q_{j}$ the corresponding averaged fluxes. Any $2N+2$ of these conservation laws are independent.

(iii) By integrating (1.5) over each of the $N$ $\mathbf{B}$-cycles we obtain, on using (2.4), $N$ equations for conservation of waves

(2.8)

$$\frac{\partial}{\partial t}k_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})=\frac{\partial}{\partial x}\omega_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}}),\qquad j=1,\dots,N\,.$$

Since (2.8) must be consistent with (1.5) we obtain a compact and physically insightful representations for $V_{j}$’s as nonlinear group velocities is (see [8], [9], [11]).

(2.9)

$$V_{j}^{(N)}=\frac{\partial\omega_{i}}{\partial\alpha_{j}}/\frac{\partial k_{i}}{\partial\alpha_{j}},\qquad\hbox{for any}\quad i=1,\dots,N.$$

We note that equations (2.8) represent the consistency conditions in the formal averaging procedure [36, 37] as well as in the the WKB-type multiple-scale expansions [24], [37], [6] leading to the same Whitham system (1.5). Within this (general) modulation theory framework the wavenumbera $k_{j}$ and the frequences $\omega_{j}$ in the modulated wave are defined as the derivatives of the phase $\eta_{j}$:

(2.10)

$$k_{j}=(\eta_{j})_{x},\quad\omega_{j}=(\eta_{j})_{t},\quad j=1,\dots,N.$$

Clearly, the definitions (2.10) by Clairaut’s theorem are consistent with the wave conservation equation (2.8).

The Whitham system (1.5) can be integrated using the Tsarev generalized hodograph transform [32]. This method was originally developed for hyperbolic systems of hydrodynamic type but is equally applicable to elliptic systems. Tsarev’s result in the application to our present problem can be formulated as follows. Any local non-constant solution of the modulation system (1.5) for a given genus $N$ is given in an implicit form by the system of $N$ algebraic equations with complex coefficients

(2.11)

$$x+V_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})t=w_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}}),\quad x+\overline{V}_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})t=\bar{w}_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}}),\qquad j=1,2,\dots,N,$$

where the characteristic speeds $V_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})\equiv V^{(N)}_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$ are given by (2.9), (2.5). The $2N$ unknown complex functions $w_{j}$, $\bar{w}_{j}$, $j=1,2,\dots,N$ satisfy the system of linear partial differential equations

(2.12)

$$\frac{\partial_{\alpha_{j}}w_{k}}{w_{k}-w_{j}}=\frac{\partial_{\alpha_{j}}V_{k}}{V_{k}-V_{j}}\quad\hbox{and \ c.c.};\quad j,k=1,2,\dots,N,\quad k\neq j\,,$$

where $\partial_{\alpha_{j}}\equiv\frac{\partial}{\partial\alpha_{j}}$. System (2.12) is overdetermined but compatible for the NLS case studied here owing to the integrability of fNLS being preserved under the Whitham averaging [6].

We now note that, for the solution $q(x,t;\varepsilon)$ of the semi-classical fNLS to have an asymptotic representation in the form of the modulated finite-band potential locally depending on $N$ “torus” phases $\varepsilon^{-1}\eta_{j}=\varepsilon^{-1}(k_{j}x+\omega_{j}t+\eta_{j}^{0})$ (see (1.4)), and at the same time to satisfy the general kinematic conditions (2.10) one must require that the “initial phases” $\eta_{j}^{0}$ are not constants but depend on $x,t$ via the branch points $\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}}$. To this end we introduce the modulation phase shift functions $\Upsilon_{j}(\boldsymbol{\alpha,\bar{\alpha}})$ by $\eta_{j}^{0}=-\Upsilon_{j}(\boldsymbol{\alpha,\bar{\alpha}})$ so that the normalized phases $\eta_{j}$ assume the form (see [10], [11])

(2.13)

$$\eta_{j}=k_{j}x+\omega_{j}t-\Upsilon_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}}).$$

Then the definition of the local wavenumber in (2.10) implies

(2.14)

$$k_{j}=\frac{\partial(k_{j}x+\omega_{j}t-\Upsilon_{j})}{\partial x},\quad j=1,2,$$

which yields

(2.15)

$$\frac{\partial k_{j}}{\partial\alpha_{m}}x+\frac{\partial\omega_{j}}{\partial\alpha_{m}}t=\frac{\partial\Upsilon_{j}}{\partial\alpha_{m}}\,,\quad\hbox{and c.c.},\quad j,m=1,2,\dots,N.$$

provided $\partial\alpha_{j}/\partial x\neq 0$, $\partial\bar{\alpha}_{j}/\partial x\neq 0$, $j=1,2$. Here $k_{j}({\boldsymbol{\alpha}},{\boldsymbol{\bar{\alpha}}})$ and $\omega_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$ are defined by (2.5), and $\Upsilon_{j}$’s are yet to be found. Note that the second condition (2.10) leads to the same set of equations (2.15). As we shall see, only half of the equations (2.15) are independent, so it is sufficient to consider either $j=1$ or $j=2$. We also note that equations (2.15) admit a compact and elegant representation in the form of the stationary phase conditions:

(2.16)

$$\frac{\partial\eta_{j}}{\partial\alpha_{m}}=0,\qquad\frac{\partial\eta_{j}}{\partial\bar{\alpha}_{m}}=0,\quad j,m=1,2,\dots,N.$$

For given functions $\Upsilon_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$ equations (2.15) fully define the modulations $\boldsymbol{\alpha}(x,t)$, $\boldsymbol{\bar{\alpha}}(x,t)$ (assuming invertibility of (2.15), which is not guaranteed a priori). Comparing equations (2.15) with the hodograph solution (2.11) and using the representation (2.9) for the characteristic speeds $V_{j}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$ in (2.11) one readily makes the identification

(2.17)

$$w_{m}=\frac{\partial_{\alpha_{m}}\Upsilon_{j}}{\partial_{\alpha_{m}}k_{j}}\quad\hbox{and c.c.},\quad j,m=1,2,\dots,N\,.$$

Now we observe that formula (2.17) must yield the same function $w_{m}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$ for all values of $j$. This is a consequence of the consistency of the genus $N$ Whitham modulation system with $N$ “extra” conservation laws (2.8) (the same argument was used to establish the ‘nonlinear group velocity’ representation (2.9) for the characteristic speeds of the Whitham modulation system). Thus, it is sufficient to consider any $N$ of the equations (2.15) for any given $j$.

Summarizing, the integration of the Whitham equations reduces to the determination of the “modulation phase shift” vector function $\boldsymbol{\Upsilon}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$. As we shall see, this function naturally arises in the RHP construction, thus enabling one to circumvent the complicated matching regularization procedure necessary for the determination of the dependence of $\boldsymbol{\Upsilon}(\boldsymbol{\alpha},\boldsymbol{\bar{\alpha}})$ on the fNLS initial conditions within the Whitham modulation theory framework. Another feature of the RHP analysis enhancing the modulation theory is that it reveals the precise mechanism of the genus change across a breaking curve. We note that within the Whitham modulation theory the genus change determination is essentially an “educated guess” process which must be confirmed by the construction of the full modulation solution.

It is well known (see, for example [39]) that the RHP (1.3) has a unique solution $m(z)$ that has asymptotics $m(z)={\bf 1}+{m_{1}\over z}+O(z^{-2})$ as $z\rightarrow\infty$, and that the solution to the NLS (1.1) is given by $\psi(x,t,\varepsilon)=-2(m_{1})_{12}$, where $(m_{1})_{12}$ denotes the $(1,2)$ entry of matrix $m_{1}$. In the case when $r_{0}(z,\varepsilon)$ has analytic continuation into the upper halfplane, the RHP for $m(z)$ is simplified by factorizing the jump matrix

(3.1)

$$V=\begin{pmatrix}1+r\bar{r}&\bar{r}\cr r&1\cr\end{pmatrix}=\begin{pmatrix}1&\bar{r}\\ 0&1\end{pmatrix}\begin{pmatrix}1&0\\ r&1\end{pmatrix}=V_{-}V_{+}~~,$$

and, thus, “splitting” the jump condition (1.3) into two jumps: one with triangular jump matrix $V_{+}$ along some contour $\Gamma^{+}$ in the upper halfplane $\bar{\mathbb{C}}^{+}$ and the other with triangular jump matrix $V_{-}$ along some contour $\Gamma_{-}$ in the lower halfplane $\bar{\mathbb{C}}^{-}$ (here we assume that ${\mathbb{R}}$ is included in $\bar{\mathbb{C}}^{\pm}$). Contours $\Gamma^{\pm}$ are deformations of ${\mathbb{R}}$. Due to the Schwarz symmetry of the ZS problem, contours $\Gamma^{\pm}$ can be chosen to be symmetric to each other with respect to the real axis, and we can restrict our attention only to the contour $\Gamma^{+}\subset\bar{\mathbb{C}}^{+}$. For simplicity, we assume $\Gamma$ to be a simple, smooth (except for a finitely many points) contour without self-intersections.

The central idea of the Deift-Zhou nonlinear steepest descent method for asymptotic analysis of RHPs is proper factorization of a jump matrix accompanied by the proper deformation of the contour $\Gamma$. That is why we assumed that $r_{0}(z;\varepsilon)$ has analytic continuation from ${\mathbb{R}}$ into $\bar{\mathbb{C}}^{+}$. Another essential element in the nonlinear steepest descent analysis is the concept of $g$-function. We use the $g$-function $g(z)=g(z;x,t)$ to define a transformation, that reduces the RHP (1.3) to another RHP that, in the small $\varepsilon$ limit, has piece-wise constant (in $z$) jump matrices. In some sense, the $g$ function can be compared with the oscillatory phase function in the WKB method, which reduces a singularly perturbed ODE (like, for example, (1.2)) to a system that can be solved by a power series in the small parameter $\varepsilon$ (or a fractional power of $\varepsilon$).

The $g$ is build for a particular solution of the fNLS (1.1), given by the corresponding scattering data. Thus, the $g$ function will be defined by a Schwarz symmetric function $f_{0}(z)$, that can be associated with a scaled logarithm of $r_{0}(z)$. We assume $f_{0}(z)$ to be analytic, or at least, piece-wise analytic, in some Schwarz symmetrical domains that contain $\Gamma^{+}$ and $\Gamma^{-}$ respectively. The meaning of piece-wise analyticity will be addressed below.

Because of Schwarz symmetry, $f_{0}(z)$ may have a purely imaginary jump $2i\Im f_{0}$ on ${\mathbb{R}}$. Depending on whether or not $\Im f_{0}(z)\equiv 0$ on some interval of ${\mathbb{R}}$, $f_{0}(z)$ may have either one or two analytic components in some region containing ${\mathbb{R}}$. In general, $f_{0}(z)$ may also depend on $\varepsilon$.

In order to control the growth of $r(z,\epsilon)$ on the contour $\Gamma$, we split it into a number of arcs of 2 different types: main arcs (bands) and complementary arcs (gaps). Main arcs are always bounded whereas complementary arcs are either unbounded (we do not consider them as they do not affect the $g$-function) or bounded arcs that neighbor main arcs at each endpoint. The number and the positions of the endpoints of these arcs depend on the initial potential as well as on a particular point $(x,t)$ of the physical (space-time) variables. The number of bounded complementary arcs (called simply complementary arcs) does not exceed the number $n\in{\mathbb{N}}$ of main arcs. Because of Schwarz symmetry, each main arc either has a complex conjugate (with anti complex conjugate orientation) or is Schwarz-symmetrical itself, and in this case it crosses the real axis. The same is true for complementary arcs. We use notations $\gamma_{m,j}$, $\gamma_{c,j}$ for the j-th pair of complex conjugate main and complementary arcs respectively (that also include single self-symmetrical arcs), whereas $\gamma^{\pm}_{m,j}$, $\gamma^{\pm}_{c,j}$ denote parts of $\gamma_{m,j}$, $\gamma_{c,j}$ in the upper or lower half-planes respectively, see Figure 1, Left. Here we discuss the general setting of the problem for the $g$-function.

The complex valued Schwarz symmetrical $g$-function is defined as satisfying the following jump and analyticity conditions:

(3.5)		$\displaystyle g_{+}+g_{-}$	$\displaystyle=f+W_{j}~~\mbox{ on the main arc $\gamma_{m,j}$, $j=0,\cdots,n$}$
(3.6)		$\displaystyle g_{+}-g_{-}$	$\displaystyle=\Omega_{j}~~\mbox{ on the complementary arc $\gamma_{c,j}$, $j=0,1,\cdots,n$ }$
(3.7)		$\displaystyle g(z)$	is analytic in ${\bar{\mathbb{C}}}\setminus\gamma$,

where the function

(3.9)

$$f(z)=f_{0}(z)+2zx+2tz^{2}$$

is a given input to the problem and the real constants $W_{j}$ and $\Omega_{j}$ are to be determined (we take $W_{0}=0$). Note that conditions (3.5) form a scalar RHP for the $g$-function $g(z)$.

By the Sokhotski-Plemelj formula,

(3.10)

$$g(z)={{R(z)}\over{2\pi i}}\sum_{j=0}^{n}\left[\int_{\gamma_{m,j}}{{f(\zeta)+W_{j}}\over{(\zeta-z)R_{+}(\zeta)}}d\zeta+\int_{\gamma_{c,j}}{{\Omega_{j}}\over{(\zeta-z)R(\zeta)}}d\zeta\right]~.$$

Here the radical $R(z)=\sqrt{\prod_{j}(z-\alpha_{j})}$ has branchcuts $\gamma_{m,j}$, where the product is taken over all the endhpoints $\alpha_{j}$ of $\gamma_{m,j}$, $j=0,1,\cdots,n$,

Let us assume for simplicity that each $\gamma_{m,j}$, $j=1,\cdots,n$, consists of two arcs and $\gamma_{m,0}$ is a single arc. Other possible configurations can be considered similarly. Then the total number of endpoints is $4n+2$ and the hyperelliptic Riemann surface $\mathcal{R}$ of $R(z)$ has the genus $N=2n$ and there are no more than $2n$ complementary arcs $\gamma_{c,j}$, $j=1,\dots,n$. We fix the branch of $R(z)$ by the requirement

(3.11)

$$\lim_{z\rightarrow\infty}\frac{R(z)}{z^{N+1}}=1$$

on the main sheet of $\mathcal{R}$. Due to the analyticity of $f$, the integrals over $\gamma_{m,j}$, $\gamma_{c,j}$ in (3.10), can be deform into the corresponding loop integrals $\hat{\gamma}_{m,j}$, $\hat{\gamma}_{c,j}$ on $\mathcal{R}$ around $\gamma_{m,j}$, $\gamma_{c,j}$ respectively. Thus,

(3.12)

$$g(z)=O(z-a_{k})^{\frac{1}{2}}$$

around any endpoint $\alpha_{k}$ of a main arc. Here and henceforth ee will refer to the endpoints $\alpha_{j}$ as branchpoints (of $\mathcal{R}$).

The scalar RHP (3.5) for $g(z;x,t)$ is completely defined if we know the genus $N$ of $\mathcal{R}$, the branchpoints $\alpha_{k}$ and the real constants $W_{j},\Omega_{j}$. From the point of view of solving RHP (3.5), the exact location of the main arcs $\gamma_{m,j}$ is not crucial, as, due to the analyticity of $f(z)$, they can be deformed without changing $g$.

The real constants in the RHP (3.5) are determined by the requirement that $g$ is analytic at $z=\infty$ and will be discussed in Section 5 below. What conditions determine the genus $N$ and the branchpoints? Without going into the details of $g$-function mechanism, we will state that the function

(3.13)

$$h(z;x,t)=2g(z;x,t)-f(z;x,t)$$

must satisfy the following “sign distribution” conditions in ${\mathbb{C}}^{+}$:

(3.14)		$\displaystyle\Im h<0$	on both sides of each main arc $\gamma_{m,j}$, $j=0,1,\cdots,n$,
(3.15)		$\displaystyle\Im h>0$	along each complementary arc $\gamma_{c,j}$, $j=1,\cdots,n$.

The corresponding inequalities in the lower half-plane, due to Schwarz symmetry, have the opposite signs.) Note that according to (3.5), (3.13), $h_{+}(z)+h_{-}(z)=2W_{j}$ on $\gamma_{m,j}$, so that all the main arcs lie on zero level curves of $\Im h(z)$. Conditions (3.14) can be violated in no more than a finite number of points. If this is the case, the corresponding value of $(x,t)=(x_{b},t_{b})$ is called a breaking point.

At any branchpoint, where a main and a complementary arcs meet (movable branchpoint), conditions (3.14), combined with (3.12), imply that

(3.17)

$$\Im h(z;x,t)=O(z-a_{k})^{\frac{3}{2}}$$

Equations (3.17), are known in the RHP literature as the modulation equations (although they essentially are solutions of the (differential) Whitham equations), determine the location of all movable (with $x,t$) branchpoints. Fixed branchpoints (“hard edges”) are known to appear in some non-analytic cases (see, for example, [17], for the box initial data). If two or more branchpoints collide at some $\alpha$, we will have

(3.18)

$$\Im h(z;x,t)=o(z-a_{k})^{\frac{3}{2}}$$

instead of (3.17). Now, the only remaining question is: how is the genus $N$ defined?

Usually, the genus $N$ is defined for some special values of $x,t$, say, for the initial data $t=0$. As we then continuously deform external parameters $x,t$, the branchpoints $\alpha_{j}$ move according to (3.17), pulling (deforming) main and complementary arcs of the contour $\gamma=\gamma(x,t)$ with them. The genus will be preserved under such deformation until a breaking point $(x_{b},t_{b})$ is reached. A regular breaking point is attained by changing the topology of the zero level of $\Im f(z;x,t)$, called pinching: the required inequalities (3.14) failed in one or several points on the main and/or complementary arcs. Such points are called double points. Another cause of breaking points is interaction of the Riemann-Hilbert contour of $g$ (the collection of the main and complementary arcs) with singularities of $f_{0}(z)$, or, in the case of piece-wise analytic $f_{0}(z)$, the interaction of the Riemann-Hilbert contour with various elements of $f_{0}(z)$.

The continuation principle ([29]) states that, in the case of a regular breaking point $(x_{b},t_{b})$, we can continue deformation of $x,t$ (with sign conditions satisfied) past $(x_{b},t_{b})$ with the appropriate change of the genus $N$.

The scalar RHP (3.5) and the modulation equations (3.17) is the starting point of our analysis. Various forms of Whitham equations from Section 2 and conservation equations can be derived from (3.5) and (3.17).

In the case of $n=0$ there are no real constants in the RHP (3.5). Then the modulation equations (3.17), defining the branchpoints of $\mathcal{R}(x,t)$, can be written as

(4.1)

$$\oint_{\hat{\gamma}_{m}}{{f(\zeta)}\over{(\zeta-\alpha)R(\zeta)_{+}}}d\zeta=0~,$$

and its complex conjugate, where $\alpha=\alpha_{0}$, $f$ is given by (3.9) and $\hat{\gamma}_{m}$ is a negatively (clockwise) oriented loop around the main arc (band) that intersects ${\mathbb{R}}$ only at some point $\mu$. Here $R(z)=\sqrt{(z-\alpha)(z-\bar{\alpha})}$.

Clearly, (4.1) and (3.9) yield

(4.2)

$$2x+2(\Re\alpha+\alpha)t=-\frac{1}{2\pi i}\oint_{\hat{\gamma}_{m}}{{f_{0}(\zeta)}\over{(\zeta-\alpha)R(\zeta)_{+}}}d\zeta~$$

together with its complex conjugate. Comparison of (4.2) and the hodograph equation (2.11) yields

(4.3)

$$V(\alpha,\bar{\alpha})=\frac{3}{2}(\alpha+\bar{\alpha}),~~~w(\alpha,\bar{\alpha})=-\frac{1}{4\pi i}\oint_{\hat{\gamma}_{m}}{{f_{0}(\zeta)}\over{(\zeta-\alpha)R(\zeta)_{+}}}d\zeta$$

and their complex conjugates, where the velocity coincides with (1.6). In the genus zero case there is only one Tsarev equation (2.12)

(4.4)

$$\frac{\partial_{\bar{\alpha}}w}{w-\bar{w}}=\frac{\partial_{\bar{\alpha}}V}{V-\bar{V}}$$

and its complex conjugate, that can be shown true by immediate calculation. Thus, we proved that $\alpha,\bar{\alpha}$ satisfy the corresponding genus zero Whitham equations (1.6). In fact, (4.1) and its complex conjugate represent the hodograph solution to the Whitham equations (1.6).

According to (3.10) and (3.13), in the genus zero case

(4.5)

$$h(z)=\frac{R(z)}{2\pi i}\oint_{\hat{\gamma}_{m}}{{f(\zeta)}\over{(\zeta-z)R(\zeta)_{+}}}d\zeta,$$

as there are no complementary arcs and the constants $W_{0}$ can be choosen as zero. Using

(4.6)

$$\frac{\partial}{\partial\alpha}\frac{R(z)}{R(\zeta)}=\frac{1}{2}\frac{\zeta-z}{(z-\alpha)(\zeta-\alpha)}\frac{R(z)}{R(\zeta)}$$

and (4.5), we obtain

(4.7)

$$\frac{\partial}{\partial\alpha}h(z)=\frac{R(z)}{2(z-\alpha)}\cdot\left.\frac{h(\zeta)}{R(\zeta)}\right|_{\zeta=a}$$

and its complex conjugate. Since, according to (4.1), the last fraction is equal to zero, we obtain

(4.8)

$$\frac{\partial}{\partial\alpha}h(z)\equiv 0\Leftrightarrow{\rm modulation~equations}.$$

According to (4.8) and (3.9), we have

(4.9)

$$h(z)=h(z;x,t,\alpha)=xh_{x}(z)+th_{t}(z)+h_{0},~~~{\rm where}~~~h_{0}=\frac{R(z)}{2\pi i}\oint_{\hat{\gamma}_{m}}{{f_{0}(\zeta)}\over{(\zeta-z)R(\zeta)_{+}}}d\zeta$$

and

(4.10)

$$h_{x}(z)=2R(z),~~~~~~~~~~~~~~~h_{t}(z)=2(z+a)R(z)$$

were calculated in [30]. Here $a=\Re\alpha$ and $h_{x},h_{t}$ are partial derivatives in the corresponding variables, which, according to (4.8), coincide with “full” derivatives on the solutions of modulation equations $\alpha=\alpha(x,t)$. Obviously,

(4.11)

$$\partial_{x}h_{t}\equiv\partial_{t}h_{x},$$

and one can check directly that this condition is equivalent to to the modulation equation (4.1). Moreover, $dh_{x}$, $dh_{t}$ are meromorphic differential of the second kind on $\mathcal{R}$ with the only poles at $\infty^{\pm}$ and, according to (2.2),

(4.12)

$$dp=dh_{x},~~~dq=dh_{t}.$$

Thus, $p=h_{x}(z),q=h_{t}(z)$ are primitive functions for the meromorphic differentials $dp,dq$, and $h(z)$ is a potential for $p,q$. So, we come to the conclusion that the topology of zero level curves of the imaginary part of the potential for $p,q$ that is defined by a particular solution, determines the change of genus of Whitham equations.

Consider the case when the (generally piece-wise) analytic function $f_{0}(z)$ does not have a jump on some interval $I\subset{\mathbb{R}}$. Let us assume for the moment that all the branchpoints $\alpha_{j}$ in the upper half-plane have distinct real parts Then, deforming the main and complementary arcs, as shown in Figure 1, Right, we can reduce the RHP (3.5) for $g(z)$ to the equivalent RHP with jumps along the vertical contours $\tilde{\gamma}_{j}$ with the branchpoints $\bar{\alpha}_{j},\alpha_{j}$, $j=0,\dots,N$:

(5.1)		$\displaystyle g_{+}+g_{-}$	$\displaystyle=f+\tilde{\eta}_{j}~~\mbox{ on the vertical contour ${\tilde{\gamma}}_{j}$, $j=0,\cdots,N$}$
(5.2)		$\displaystyle g(z)$	is analytic in ${\bar{\mathbb{C}}}\setminus\cup{\tilde{\gamma}}_{j}$,

where the real constants $\tilde{\eta}_{2j}$ are expressed as differences of constants $W,\Omega$ on the neighboring main and complementary arcs. For example, we have $\tilde{\eta}_{2j}=W_{j}-\Omega_{j+1}$, $\tilde{\eta}_{2j+1}=W_{j+1}-\Omega_{j+1}$ for the configuration, shown on Figure 1, Right. The $g$ functions, defined by the RHPs (3.5) and (5.1) coincide outside the region, encircled by the “old” main and complementary arcs and “new” vertical cuts (called bands), see Figure 1, Right. Inside this region, the “new” $g(z)$ coincides with: the “old” $g(z)$ up to appropriate real constant $\Omega_{j}$ on the positive (left) side of any two neighboring bands ${\tilde{\gamma}}$; the values of the “old” $g(z)$ from the second sheet of $\mathcal{R}$ up to appropriate real constant $W_{j}$ on the negative (right) side of any two neighboring bands ${\tilde{\gamma}}$ (note the opposite orientation of the neighboring bands). Thus, the modulation equations (3.17) remains valid for the “new” $g$, so that either RHP (3.5) or RHP (5.1) can be used to calculate the branchpoints $\alpha_{j}$, the Whitham equations, conservation equations, etc.

Introducing now $\tilde{g}=g-\frac{\tilde{\eta}_{0}}{2}$, we obtain the RHP

(5.4)		$\displaystyle\tilde{g}_{+}+\tilde{g}_{-}$	$\displaystyle=f+\eta_{j}~~\mbox{ on the band ${\tilde{\gamma}}_{j}$, $j=0,\cdots,2n$}$
(5.5)		$\displaystyle\tilde{g}(z)$	is analytic in ${\bar{\mathbb{C}}}\setminus\gamma$,

where $\eta_{2j}=\tilde{\eta}_{2j}-\tilde{\eta}_{0}$. In the configuration of Figure 1, Right, $\eta_{2j}=W_{j}-\Omega_{j+1}+\Omega_{1}$, $\eta_{2j+1}=W_{j+1}-\Omega_{j+1}+\Omega_{1}$. The constants $\boldsymbol{\eta}=(\eta_{1},\dots,\eta_{2n})$ are exactly the constants from the argument of the multi-phase nonlinear wave solution (given through the Riemann Theta function) of the model problem, see [30], Section 8. They are also introduced in (1.4) as fundamental phases. With a mild abuse of notation, we will use $g(z)$ instead of $\tilde{g}(z)$ below and use $\mathcal{R}$ to denote the hyperelliptic Riemann surface with branchcuts along ${\tilde{\gamma}}_{j}$. Here and henceforth we also assume that all the bands are oriented upwards. This orientation does not change any jump conditions (5.4).