Evolution of initial discontinuities in the DNLS equation theory

The problem of classification of wave structures evolving from initial discontinuities has played important role since the classical paper of Riemann [1]. Complemented by the jump conditions of Rankine [2] and Hugoniot [3, 4], it provided a prototypical example of formation of shocks in dispersionless media with small viscosity, and the full classification of possible wave patterns evolving from initial discontinuities with general initial data in adiabatic flows of ideal gas was obtained by Kotchine [5]. These results were generalized to the class of so-called genuinely nonlinear hyperbolic systems (see, e.g., [6, 7]), however, situation beyond this class is much more complicated and suffers from ambiguity of possible solutions. One of the methods to remove this ambiguity is introduction of small viscosity into equations followed by taking the limit of zero viscosity. This approach seems very natural from physical point of view since it provides some information on the inner structure of viscous shocks. At the same time, there exists another method of regularization of hydrodynamics-like equations, namely, introduction of small dispersion. Although in this case the limit of zero dispersion does not lead to the same shock structure, this approach is of considerable interest since, on one side, it is related with the theory of dispersive shock waves (DSWs) that finds a number of physical applications (see, e.g., review article [8] and references therein) and, on another hand, there are situations when the regularized equation belongs to the class of completely integrable equations and therefore it admits quite thorough investigation including even cases of non-genuinely nonlinear hyperbolic systems.

The simplest example of dispersive nonlinear evolution equation is apparently the famous KdV equation and in this case the solution of the Riemann problem is extremely simple: Gurevich and Pitaevskii showed [9] with the use of Whitham modulation theory [10] that there are only two possible ways of evolution of initial discontinuity—it can evolve into either rarefaction wave or DSW whose parameters can be expressed in explicit analytical form by solving the Whitham equations. This result was obtained without explicit use of the complete integrability of the KdV equation [11], but its extension to the NLS equation became possible [12] only after derivation of the Whitham modulation equations [13, 14] by the methods based on the inverse scattering transform for the NLS equation [15] which means its complete integrability. It was shown in [12] that the NLS equation evolution of any initial discontinuity leads to a wave pattern consisting of a sequence of building blocks two of which are represented by either the rarefaction wave or the DSW, and they are separated by plateau, or vacuum, or two-phase self-similar solution close to unmodulated nonlinear periodic wave. The rarefaction waves are here self-similar simple wave solutions of the dispersionless limit of the NLS equation (i.e., of the shallow water equations) and DSW is described by a self-similar solution of the Whitham modulation equations. In total, there are six different possible wave patterns that can evolve from a given initial discontinuity. Similar classification of wave patterns was also established for the dispersive shallow water Kaup-Boussinesq equation [16, 17].

For classification of wave patterns arising in solutions of the Riemann problem of the KdV or NLS type, it is important that the corresponding dispersionless limits (Hopf equation or shallow water equations) are represented by the genuinely nonlinear hyperbolic equations. If it is not the case, then the classification of the KdV-NLS type becomes insufficient and it was found that it should include new elements—kinks or trigonometric dispersive shocks—for mKdV [18] and Miyata-Camassa-Choi [19] equations. The mKdV equation is a modification of KdV equation and it also describes a unidirectional propagation of wave with a single field variable, so it can be considered as a simplest example of non-convex dispersive hydrodynamics. In spite of its relative simplicity, the full classification of the wave patterns in the solution of the Riemann problem is much more complicated than that in the KdV equation case and it was achieved in [20] for the Gardner equation (related with the mKdV equation) with the use of Riemann invariant form of the Whitham modulation equations obtained in [21]. These results were adapted to mKdV equation in [22] and for this equation the Whitham modulation equations were obtained by the direct Whitham method in [23]. Instead of two possible patterns in KdV case, in the mKdV-Gardner case we have eight possible wave structures which depend now not only on the sign of the jump at the discontinuity, but also on the values of wave amplitudes at both its sides. No similar classification has been obtained yet for two-directional waves although important partial results were obtained in [19] for the Miyata-Camassa-Choi equation. However, this equation is not completely integrable and although the principles of such a classification are the same for completely integrable and non-integrable equations and in non-integrable case some essential information can be derived with the use of El's method [24], we prefer here to turn first to the case of completely integrable derivative nonlinear Schrödinger (DNLS) equation when more complete study is possible.

Thus, in this paper, we shall give full solution of the Riemann problem for evolution of initial discontinuities in the theory of the DNLS equation

$$\begin{eqnarray}&&{\rm{i}}{{\rm{\Psi }}}_{t}+\displaystyle \frac{1}{2}{{\rm{\Psi }}}_{{xx}}-{\rm{i}}{(| {\rm{\Psi }}{| }^{2}{\rm{\Psi }})}_{x}=0.\end{eqnarray} \tag{ 1 }$$

This equation appears in the theory of nonlinear Alfvén waves in plasma physics (see, e.g., [25] and references therein) and in nonlinear optics (see, e.g., [26] and references therein). Its complete integrability was established in [27, 28], periodic solution and Whitham modulation equations were derived in [29, 30]. Partial solution of the Riemann problem was obtained in [33], however, only in the sector of the NLS equation type structures. Here we develop the method which permits one to predict the wave pattern arising from any given data for an initial discontinuity. The method is quite general and it was applied to the generalized NLS equations [34] with Kerr-type cubic nonlinearity added to (1), what is important for nonlinear optics applications, and to the Landau–Lifshitz equation for magnetics with easy-plane anisotropy [35]. Here we develop a similar theory for the equation (1).

In many situations, it is convenient to transform the DNLS equation (1) to the so-called hydrodynamic form what is achieved by means of the substitution

$$\begin{eqnarray}&&{\rm{\Psi }}(x,t)=\sqrt{\rho (x,t)}\exp \left({\rm{i}}{\int }^{x}u(x^{\prime} ,t){\rm{d}}x^{\prime} \right).\end{eqnarray} \tag{ 2 }$$

After separation of real and imaginary parts, this equation is easily reduced to the system

$$\begin{eqnarray}&&{\rho }_{t}+{\left[\rho \left(u-\displaystyle \frac{3}{2}\rho \right)\right]}_{x}=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{u}_{t}+{{uu}}_{x}-{(\rho u)}_{x}+{\left(\displaystyle \frac{{\rho }_{x}^{2}}{8{\rho }^{2}}-\displaystyle \frac{{\rho }_{{xx}}}{4\rho }\right)}_{x}=0.\end{eqnarray} \tag{ 4 }$$

These equations can be interpreted as hydrodynamic form of the DNLS equation with equation (3) playing the role of the continuity equation and equation (4) of the Euler equation for a fluid with depending on the flow velocity 'pressure' ρu and 'quantum pressure' represented by the last term. However, one should keep in mind that we are dealing with an anisotropic medium where the flux of mass in (3) does not coincide with the momentum density. As a result, the conservation of momentum equation takes the form

$$\begin{eqnarray}&&{[\rho (u-\rho )]}_{t}+{\left[\rho {u}^{2}-3u{\rho }^{2}+2{\rho }^{3}-\displaystyle \frac{1}{4}\rho {(\mathrm{ln}\rho )}_{{xx}}\right]}_{x}=0.\end{eqnarray} \tag{ 5 }$$

This feature of the DNLS equation, which in our case means that the 'right' and 'left' directions of wave propagation cannot be exchanged by an inversion operation x → −x, can be illustrated by the linear approximation.

Let us consider linear waves propagating along the background flow (ρ₀, u₀), that is ρ = ρ₀ + ρ', u = u₀ + u', where $| \rho ^{\prime} | \ll {\rho }_{0}$, $| u^{\prime} | \ll {u}_{0}$. Linearization with respect to small variables ρ', u' yields the system

$$\begin{eqnarray}&&{\rho }_{t}^{{\prime} }+({u}_{0}-3{\rho }_{0}){\rho }_{x}^{{\prime} }+{\rho }_{0}{u}_{x}^{{\prime} }=0\\ &&{u}_{t}^{{\prime} }+({u}_{0}-{\rho }_{0}){u}_{x}^{{\prime} }-\displaystyle \frac{1}{4{\rho }_{0}}{\rho }_{{xxx}}^{{\prime} }=0.\end{eqnarray} \tag{ 6 }$$

Looking for the plane wave solution $\rho ^{\prime} ,u^{\prime} \propto \exp [{\rm{i}}({kx}-\omega t)]$, we find that it exists if only the dispersion law

$$\begin{eqnarray}&&\omega (k)=k[{u}_{0}-2{\rho }_{0}\pm \sqrt{{\rho }_{0}({\rho }_{0}-{u}_{0})+{k}^{2}/4}]\end{eqnarray} \tag{ 7 }$$

is fulfilled. In the limit of small wave vectors k we find

$$\begin{eqnarray}&&\omega (k)\approx ({u}_{0}-2{\rho }_{0}\pm \sqrt{{\rho }_{0}({\rho }_{0}-{u}_{0})})k\pm \displaystyle \frac{{k}^{3}}{8\sqrt{{\rho }_{0}({\rho }_{0}-{u}_{0})}}.\end{eqnarray} \tag{ 8 }$$

As we see, there are two modes of propagation of linear waves with different absolute values of propagation velocities even for medium at rest with u₀ = 0: the initial disturbance decays to two wave packets propagating with different absolute values of group velocities.

Another important feature of the dispersion law (7) is that it leads to modulationally unstable modes with complex ω for u₀ > ρ₀. In this paper, we shall confine ourselves to the stable situations only.

The above properties of the wave propagation in the DNLS equation theory are preserved in the weakly nonlinear cases, that is if we take into account weak nonlinear effects in the above modes with ρ' small but finite. Before proceeding to this task, we shall consider in the next section the dispersionless dynamics when the dispersion effects are completely neglected.

The nonlinear and dispersive effects have the same order of magnitude, when in equations (3), (4) we have ${u}^{2}\sim \rho u\sim {\rho }_{{xx}}/\rho$, hence the last term in equation (4) can be neglected if the variables ρ and u change little on distances Δx ∼ 1/ρ. In this dispersionless approximation, the flow is governed by the equations

$$\begin{eqnarray}&&{\rho }_{t}+{\left[\rho \left(u-\displaystyle \frac{3}{2}\rho \right)\right]}_{x}=0,\\ &&\quad {u}_{t}+{{uu}}_{x}-{(\rho u)}_{x}=0\end{eqnarray} \tag{ 9 }$$

or

$$\begin{eqnarray}&&{\left(\begin{array}{c}\rho \\ u\end{array}\right)}_{t}+{\mathbb{A}}{\left(\begin{array}{c}\rho \\ u\end{array}\right)}_{x}=0,\quad {\mathbb{A}}=\left(\begin{array}{cc}u-3\rho & \rho \\ -u & u-\rho \end{array}\right).\end{eqnarray} \tag{ 10 }$$

The characteristic velocities of this system

$$\begin{eqnarray}&&{v}_{\pm }=u-2\rho \pm \sqrt{\rho (\rho -u)}\end{eqnarray} \tag{ 11 }$$

coincide, naturally, with the phase velocities ${\omega /k| }_{k\to 0}$ for the dispersion laws (7) in the long wave limit. The system (10) of first-order equations can be easily transformed to a diagonal form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {r}_{+}}{\partial t}+{v}_{+}\displaystyle \frac{\partial {r}_{+}}{\partial x}=0,\quad \displaystyle \frac{\partial {r}_{-}}{\partial t}+{v}_{-}\displaystyle \frac{\partial {r}_{-}}{\partial x}=0\end{eqnarray} \tag{ 12 }$$

for the Riemann invariants

$$\begin{eqnarray}&&{r}_{+}=u/2-\rho +\sqrt{\rho (\rho -u)},\quad {r}_{-}=u/2-\rho -\sqrt{\rho (\rho -u)},\end{eqnarray} \tag{ 13 }$$

with the velocities (11) expressed in terms of the Riemann invariants as

$$\begin{eqnarray}&&{v}_{+}=\displaystyle \frac{3}{2}{r}_{+}+\displaystyle \frac{1}{2}{r}_{-},\quad {v}_{-}=\displaystyle \frac{1}{2}{r}_{+}+\displaystyle \frac{3}{2}{r}_{-}.\end{eqnarray} \tag{ 14 }$$

If the solution of equation (12) is known, then the physical variables ρ, u are given by the expressions

$$\begin{eqnarray}&&\rho =\displaystyle \frac{1}{2}{(\sqrt{-{r}_{+}}\pm \sqrt{-{r}_{-}})}^{2},\qquad u=\pm 2\sqrt{{r}_{+}{r}_{-}},\end{eqnarray} \tag{ 15 }$$

where both Riemann invariants are negative: ${r}_{-}\leqslant {r}_{+}\leqslant 0$.

The Riemann invariants (13) and the characteristic velocities (11) are real for ρ ≥ u (ρ ≥ 0 by definition), that is the inequalities ρ ≥ 0, ρ ≥ u define the hyperbolicity domain in the plane (u, ρ) of physical variables. Besides that, it is extremely important that the Riemann invariant r₊ reaches its maximal value r₊ = 0 along the ρ-axis where u = 0. It means that its dependence on the physical variables is not monotonous. We say that the ρ-axis u = 0 cuts the hyperbolicity domain into two monotonicity regions u < 0 and u > 0. Correspondingly, the dependence of the physical variables on the Riemann invariants is not single-valued—it is two-valued in our case of a single maximum of r₊, if the solution of our hydrodynamics equations crosses the axis u = 0. As we shall see, this leads to important consequences in classification of wave structures evolving from initial discontinuities.

Now we turn to derivation of the evolution equations for weakly nonlinear waves with small dispersion.

The linear modes correspond to flows with fixed relationship between ρ' and u' and generalizations of these waves to the nonlinear regime are simple waves with one of the Riemann invariants r_± constant. In the leading order, when the nonlinear and dispersive corrections are accounted in their main approximations, we can add their effects in the resulting evolution equations. The small dispersive effects are described by the last terms in the dispersion laws (8) that can be transformed to the differential equations for ρ' by the replacements $\omega \to {\rm{i}}{\partial }_{t}$, $k\to -{\rm{i}}{\partial }_{x}$:

$$\begin{eqnarray}&&{\rho }_{t}^{{\prime} }+({u}_{0}-2{\rho }_{0}\pm \sqrt{{\rho }_{0}({\rho }_{0}-{u}_{0})}){\rho }_{x}^{{\prime} }\mp \displaystyle \frac{1}{8\sqrt{{\rho }_{0}({\rho }_{0}-{u}_{0})}}{\rho }_{{xxx}}^{{\prime} }=0.\end{eqnarray} \tag{ 16 }$$

Therefore it is enough to consider now the weak nonlinear effects neglecting the dispersion. To simplify the notation, we shall consider waves propagating along a uniform quiescent background with ρ = ρ₀, u = u₀ = 0.

4.1. Kortweg-de Vries mode

At first we shall consider waves with r₊  = const, and it is easy to find that far enough from a localized wave pulse this Riemann invariant vanishes and the identity r₋ = 0 is fulfilled with the accuracy up to the first order of small quantities ρ' and u. Consequently, the equation for r₊ is already satisfied with this accuracy and for the waves of density ρ' we can substitute u = 0 into dispersionless expressions (11) and (13) for r₋ and v₋, correspondingly, to find

$$\begin{eqnarray*}&&{r}_{-}\approx -2({\rho }_{0}+\rho ^{\prime} ),\qquad {v}_{-}\approx -3({\rho }_{0}+\rho ^{\prime} ).\end{eqnarray*}$$

Thus, dispersionless Hopf equation for this mode obtained from (12) reads

$$\begin{eqnarray*}&&{\rho }_{t}^{{\prime} }-3({\rho }_{0}+\rho ^{\prime} ){\rho }_{x}^{{\prime} }=0,\end{eqnarray*}$$

and addition of dispersion term from (16) for lower sign yields the KdV equation

$$\begin{eqnarray}&&{\rho }_{t}^{{\prime} }-3({\rho }_{0}+\rho ^{\prime} ){\rho }_{x}^{{\prime} }-\displaystyle \frac{1}{8{\rho }_{0}}{\rho }_{{xxx}}^{{\prime} }=0.\end{eqnarray} \tag{ 17 }$$

Solution of the Riemann problem for this equation has very simple Gurevich-Pitaevskii type [9].

4.2. Modified Korteweg–de Vries mode

In the mode with r₋ = −2ρ₀ = const we have to make calculations with accuracy up to the second order with respect to ρ'. The condition r₋ = −2ρ₀ gives us the relationship

$$\begin{eqnarray*}&&u\approx 2\rho ^{\prime} -\displaystyle \frac{\rho {{\prime} }^{2}}{{\rho }_{0}},\end{eqnarray*}$$

and its substitution into expressions (11) and (13) for r₊ and v₊ yields with the same accuracy

$$\begin{eqnarray*}&&{r}_{+}\approx -\displaystyle \frac{\rho {{\prime} }^{2}}{2{\rho }_{0}},\qquad {v}_{+}\approx -{\rho }_{0}-\displaystyle \frac{\rho {{\prime} }^{2}}{{\rho }_{0}}.\end{eqnarray*}$$

Hence equation (12) for r₊ reduses to the dispersionless equation for the density

$$\begin{eqnarray*}&&{\rho }_{t}^{{\prime} }-\left({\rho }_{0}+\displaystyle \frac{\rho {{\prime} }^{2}}{{\rho }_{0}}\right){\rho }_{x}^{{\prime} }=0,\end{eqnarray*}$$

and addition of dispersion term from (16) for upper sign yields the mKdV equation

$$\begin{eqnarray}&&{\rho }_{t}^{{\prime} }-\left({\rho }_{0}+\displaystyle \frac{\rho {{\prime} }^{2}}{{\rho }_{0}}\right){\rho }_{x}^{{\prime} }+\displaystyle \frac{1}{8{\rho }_{0}}{\rho }_{{xxx}}^{{\prime} }=0.\end{eqnarray} \tag{ 18 }$$

For this mode the solution of the Riemann problem [20, 22] is much more complicated and this fact suggests that the Riemann problem for the DNLS equation must differ considerably from that for the NLS equation [12]. To find this solution, we have to obtain the periodic solutions in convenient for us form parameterized by the Riemann invariants of the Whitham modulation equations and to derive these modulation equations. Actually, that was done in [29, 30], however, for completeness we shall reproduce here briefly these results with some improvements.

The finite-gap integration method (see, e.g., [31]) of finding periodic solutions is based on possibility of representing the DNLS equation (1) as a compatibility condition of two systems of linear equations with a spectral parameter λ

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}\left(\begin{array}{c}{\psi }_{1}\\ {\psi }_{2}\end{array}\right)={\mathbb{U}}\left(\begin{array}{c}{\psi }_{1}\\ {\psi }_{2}\end{array}\right),\quad \displaystyle \frac{\partial }{\partial t}\left(\begin{array}{c}{\psi }_{1}\\ {\psi }_{2}\end{array}\right)={\mathbb{V}}\left(\begin{array}{c}{\psi }_{1}\\ {\psi }_{2}\end{array}\right),\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}{\mathbb{U}}=\left(\begin{array}{cc}F & G\\ H & -F\end{array}\right),\quad {\mathbb{V}}=\left(\begin{array}{cc}A & B\\ C & -A\end{array}\right),\end{eqnarray} \tag{ 20 }$$

with

$$\begin{eqnarray}F & = & -2{\rm{i}}{\lambda }^{2},\quad G=2\lambda {\rm{\Psi }},\quad H=2\lambda {{\rm{\Psi }}}^{* },\\ A & = & -{\rm{i}}(4{\lambda }^{4}+2{\lambda }^{2}| {\rm{\Psi }}{| }^{2}),\quad B=4{\lambda }^{3}{\rm{\Psi }}+\lambda ({\rm{i}}{{\rm{\Psi }}}_{x}+2| {\rm{\Psi }}{| }^{2}{\rm{\Psi }}),\\ C & = & 4{\lambda }^{3}{{\rm{\Psi }}}^{* }-\lambda ({\rm{i}}{{\rm{\Psi }}}_{x}^{* }-2| {\rm{\Psi }}{| }^{2}{{\rm{\Psi }}}^{* }).\end{eqnarray} \tag{ 21 }$$

The compatibility condition of linear systems (19),

$$\begin{eqnarray}&&{{\mathbb{U}}}_{t}-{{\mathbb{V}}}_{x}+[{\mathbb{U}},{\mathbb{V}}]=0,\end{eqnarray} \tag{ 22 }$$

where [·, ·] is a commutator of matrices, is equivalent to the DNLS equation.

If we denote as (ψ₁, ψ₂) and $(\bar{{\psi }_{1}},{\bar{\psi }}_{2})$ two basis solutions of linear systems (19) and introduce a matrix of 'squared basis functions'

$$\begin{eqnarray}{\mathbb{W}}=\left(\begin{array}{cc}-{\rm{i}}f & g\\ h & {\rm{i}}f\end{array}\right),\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&f=-\displaystyle \frac{{\rm{i}}}{2}({\psi }_{1}{\bar{\psi }}_{2}+{\psi }_{2}{\bar{\psi }}_{1}),\quad g={\psi }_{1}{\bar{\psi }}_{1},\quad h=-{\psi }_{2}{\bar{\psi }}_{2},\end{eqnarray} \tag{ 24 }$$

then equations for these functions can be written in matrix form

$$\begin{eqnarray}&&{{\mathbb{W}}}_{x}=[{\mathbb{U}},{\mathbb{W}}],\qquad {{\mathbb{W}}}_{t}=[{\mathbb{V}},{\mathbb{W}}].\end{eqnarray} \tag{ 25 }$$

It is known that the characteristic polynomial

$$\begin{eqnarray}&&\det ({iw}\cdot {\mathbb{I}}-{\mathbb{W}})=-{w}^{2}+{f}^{2}-{gh}\end{eqnarray} \tag{ 26 }$$

does not depend on t and x (in our simple case it can be checked by a simple calculation and the general proof of this theorem can be found, e.g., in appendix B of [32]). Hence, it defines the curve

$$\begin{eqnarray}&&{w}^{2}=P(\lambda ),\qquad P(\lambda )={f}^{2}-{gh},\end{eqnarray} \tag{ 27 }$$

where $P(\lambda )$ depends on $\lambda$ only.

Periodic solutions are distinguished by the condition that $P(\lambda )$ be a polynomial in $\lambda$, and then the structure of the matrix elements (21) suggests that f, g, h must also be polynomials in $\lambda$. The simplest one-phase solution corresponds to the polynomials f, g, h in the form

$$\begin{eqnarray}&&f={\lambda }^{4}-{f}_{1}{\lambda }^{2}+{f}_{2},\quad g=\lambda ({\lambda }^{2}-\mu /2){\rm{\Psi }},\quad h=\lambda ({\lambda }^{2}-{\mu }^{* }/2){{\rm{\Psi }}}^{* }.\end{eqnarray} \tag{ 28 }$$

The functions f₁(x, t), f₂(x, t), μ(x, t) and ${\mu }^{* }(x,t)$ are unknown yet, but we shall see soon that μ(x, t) and ${\mu }^{* }(x,t)$ are complex conjugate, whence the notation. Then the polynomial $P(\lambda )$ can be written as

$$\begin{eqnarray}&&P(\lambda )=\displaystyle \prod _{i=1}^{4}({\lambda }^{2}-{\lambda }_{i}^{2})={\lambda }^{8}-{s}_{1}{\lambda }^{6}+{s}_{2}{\lambda }^{4}-{s}_{3}{\lambda }^{2}+{s}_{4},\end{eqnarray} \tag{ 29 }$$

where s_i are symmetric functions of the four zeroes ${\lambda }_{i}^{2}$ of the polynomial,

$$\begin{eqnarray}&&{s}_{1}=\displaystyle \sum _{i}{\lambda }_{i}^{2},\quad {s}_{2}=\displaystyle \sum _{i\lt j}{\lambda }_{i}^{2}{\lambda }_{j}^{2},\quad {s}_{3}=\displaystyle \sum _{i\lt j\lt k}{\lambda }_{i}^{2}{\lambda }_{j}^{2}{\lambda }_{k}^{2},\quad {s}_{4}={\lambda }_{1}^{2}{\lambda }_{2}^{2}{\lambda }_{3}^{2}{\lambda }_{4}^{2},\end{eqnarray} \tag{ 30 }$$

and the identity (27) yields the conservation laws

$$\begin{eqnarray}{s}_{1} & = & 2{f}_{1}+\nu ,\quad {s}_{2}={f}_{1}^{2}+{{sf}}_{2}+\nu (\mu +{\mu }^{* }),\\ {s}_{3} & = & 2{f}_{1}{f}_{2}-\nu \mu {\mu }^{* },\quad {s}_{4}={f}_{2}^{2},\end{eqnarray} \tag{ 31 }$$

where $\nu =| {\rm{\Psi }}{| }^{2}$. This system permits one to express f₁(x, t), f₂(x, t), μ(x, t) and ${\mu }^{* }(x,t)$ as functions of ν:

$$\begin{eqnarray}&&{f}_{1}=({s}_{1}-\nu )/2,\qquad {f}_{2}=\pm \sqrt{{s}_{4}},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\mu ,{\mu }^{* }=\displaystyle \frac{1}{8\nu }(4{s}_{2}\pm 8\sqrt{{s}_{4}}-{(\nu -{s}_{1})}^{2}+i\sqrt{-{ \mathcal R }(\nu )}),\end{eqnarray} \tag{ 33 }$$

where the polynomial

$$\begin{eqnarray}{ \mathcal R }(\nu ) & = & {\nu }^{4}-4{s}_{1}{\nu }^{3}+(6{s}_{1}^{2}-8{s}_{2}\pm 48\sqrt{{s}_{4}}){\nu }^{2}\\ & & -(4{s}_{1}^{3}-16{s}_{1}{s}_{2}+64{s}_{3}\pm 32{s}_{1}\sqrt{{s}_{4}})\nu +{(-{s}_{1}^{2}+4{s}_{2}\pm 8\sqrt{{s}_{4}})}^{2}\end{eqnarray} \tag{ 34 }$$

is called a resolvent of the polynomial $P(\lambda )$ since its zeroes ν_i are related with the zeroes ${\lambda }_{j}$ of $P(\lambda )$ by symmetric formulae: the upper signs (+) in (34) corresponds to the zeroes

$$\begin{eqnarray}{\nu }_{1} & = & {(-{\lambda }_{1}+{\lambda }_{2}+{\lambda }_{3}+{\lambda }_{4})}^{2},\quad {\nu }_{2}={({\lambda }_{1}-{\lambda }_{2}+{\lambda }_{3}+{\lambda }_{4})}^{2},\\ {\nu }_{3} & = & {({\lambda }_{1}+{\lambda }_{2}-{\lambda }_{3}+{\lambda }_{4})}^{2},\quad {\nu }_{4}={({\lambda }_{1}+{\lambda }_{2}+{\lambda }_{3}-{\lambda }_{4})}^{2},\end{eqnarray} \tag{ 35 }$$

and the lower signs (−) in equation (34) correspond to the zeroes

$$\begin{eqnarray}{\nu }_{1} & = & {(-{\lambda }_{1}+{\lambda }_{2}+{\lambda }_{3}-{\lambda }_{4})}^{2},\quad {\nu }_{2}={({\lambda }_{1}-{\lambda }_{2}+{\lambda }_{3}-{\lambda }_{4})}^{2},\\ {\nu }_{3} & = & {({\lambda }_{1}+{\lambda }_{2}-{\lambda }_{3}-{\lambda }_{4})}^{2},\quad {\nu }_{4}={({\lambda }_{1}+{\lambda }_{2}+{\lambda }_{3}+{\lambda }_{4})}^{2}.\end{eqnarray} \tag{ 36 }$$

This can be proved by a simple check of the Viète formulae. In both cases the zeroes are ordered according to ν₁ ≤ ν₂ ≤ ν₃ ≤ ν₄ for λ₁ ≤ λ₂ ≤ λ₃ ≤ λ₄ ≤ 0.

From the components

$$\begin{eqnarray}&&{g}_{x}=2{\rm{i}}{Gf}+2{Fg},\qquad {g}_{t}=2{\rm{i}}{Bf}+2{Ag}\end{eqnarray} \tag{ 37 }$$

of the matrix equation (25) at $\lambda ={\mu }^{1/2}$ we find that μ satisfies the equations

$$\begin{eqnarray}&&{\mu }_{x}=4{\rm{i}}\sqrt{P({\mu }^{1/2})},\quad {\mu }_{t}=8{\rm{i}}(2{f}_{1}+\nu )\sqrt{P({\mu }^{1/2})}=2{s}_{1}{\mu }_{x},\end{eqnarray} \tag{ 38 }$$

where we have used the first equation (31). Consequently, μ depends on the phase ξ = x − Vt only, where $V=-2{s}_{1}=-2\sum {\lambda }_{i}^{2}$. Then the variable ν also depends on ξ only. Substitution of $g=\lambda ({\lambda }^{2}-\mu ){\rm{\Psi }}$ into the first equation (37) gives ${{\rm{\Psi }}}_{x}=-4{\rm{i}}{\rm{\Psi }}({f}_{1}-\mu )$, so that ${\nu }_{x}=4{\rm{i}}\nu (\mu -{\mu }^{* })$, and, with the use of (33), we obtain equation for ν,

$$\begin{eqnarray}&&\displaystyle \frac{d\nu }{d\xi }=\sqrt{-{ \mathcal R }(\nu )},\quad \xi =x-{Vt},\quad V=-\sum {\lambda }_{i}^{2}=\displaystyle \frac{1}{4}\sum {\nu }_{i}.\end{eqnarray} \tag{ 39 }$$

The real solutions of this equation correspond to oscillations of ν within the intervals where $-{ \mathcal R }(\nu )\geqslant 0$. We shall discuss two possibilities separately.

(A) At first we shall consider the periodic solution corresponding to oscillations of ν in the interval

$$\begin{eqnarray}&&{\nu }_{1}\leqslant \nu \leqslant {\nu }_{2}.\end{eqnarray} \tag{ 40 }$$

Standard calculation yields after some algebra the solution in terms of Jacobi elliptic functions:

$$\begin{eqnarray}&&\nu ={\nu }_{2}-\displaystyle \frac{({\nu }_{2}-{\nu }_{1}){\mathrm{cn}}^{2}(\theta ,m)}{1+\tfrac{{\nu }_{2}-{\nu }_{1}}{{\nu }_{4}-{\nu }_{2}}{\mathrm{sn}}^{2}(\theta ,m)},\end{eqnarray} \tag{ 41 }$$

where it is assumed that ν(0) = ν₁, and

$$\begin{eqnarray}&&\theta =\sqrt{({\nu }_{3}-{\nu }_{1})({\nu }_{4}-{\nu }_{2})}\,\xi /2,\qquad m=\displaystyle \frac{({\nu }_{4}-{\nu }_{3})({\nu }_{2}-{\nu }_{1})}{({\nu }_{4}-{\nu }_{2})({\nu }_{3}-{\nu }_{1})},\end{eqnarray} \tag{ 42 }$$

$\mathrm{cn}$ and $\mathrm{sn}$ being Jacobi elliptic functions. The wavelength of the oscillating function (41) is

$$\begin{eqnarray}&&L=\displaystyle \frac{4K(m)}{\sqrt{({\nu }_{3}-{\nu }_{1})({\nu }_{4}-{\nu }_{2})}}=\displaystyle \frac{K(m)}{\sqrt{({\lambda }_{3}^{2}-{\lambda }_{1}^{2})({\lambda }_{4}^{2}-{\lambda }_{2}^{2})}},\end{eqnarray} \tag{ 43 }$$

where K(m) is the complete elliptic integral of the first kind.

In the limit ${\nu }_{3}\to {\nu }_{2}$ ($m\to 1$) the wavelength tends to infinity and the solution (41) acquires the soliton form

$$\begin{eqnarray}&&\nu ={\nu }_{2}-\displaystyle \frac{{\nu }_{2}-{\nu }_{1}}{{\cosh }^{2}\theta +\tfrac{{\nu }_{2}-{\nu }_{1}}{{\nu }_{4}-{\nu }_{2}}{\sinh }^{2}\theta }.\end{eqnarray} \tag{ 44 }$$

This is a 'dark soliton' for the variable ν.

The limit $m\to 0$ can be reached in two ways.

• (i)  
  If ${\nu }_{2}\to {\nu }_{1}$, then the solution transforms into a linear harmonic wave

  $$\begin{eqnarray}&&\nu \cong {\nu }_{2}-\displaystyle \frac{1}{2}({\nu }_{2}-{\nu }_{1})\cos (k\xi ),\quad k=\sqrt{({\nu }_{3}-{\nu }_{1})({\nu }_{4}-{\nu }_{1})}.\end{eqnarray} \tag{ 45 }$$
• (ii)  
  If ${\nu }_{4}={\nu }_{3}$ but ${\nu }_{1}\ne {\nu }_{2}$, then then we arrive at the nonlinear trigonometric solution:

  $$\begin{eqnarray}&&\nu ={\nu }_{2}-\displaystyle \frac{({\nu }_{2}-{\nu }_{1}){\cos }^{2}\theta }{1+\tfrac{{\nu }_{2}-{\nu }_{1}}{{\nu }_{3}-{\nu }_{2}}{\sin }^{2}\theta },\quad \theta =\sqrt{({\nu }_{3}-{\nu }_{1})({\nu }_{3}-{\nu }_{2})}\,\xi /2.\end{eqnarray} \tag{ 46 }$$

If we take the limit ${\nu }_{2}-{\nu }_{1}\ll {\nu }_{3}-{\nu }_{1}$ in this solution, then we return to the small-amplitude limit (45) with ${\nu }_{4}={\nu }_{3}$. On the other hand, if we take here the limit ${\nu }_{2}\to {\nu }_{3}={\nu }_{4}$, then the argument of the trigonometric functions becomes small and we can approximate them by the first terms of their series expansions. This corresponds to an algebraic soliton of the form

$$\begin{eqnarray}&&\nu ={\nu }_{2}-\displaystyle \frac{{\nu }_{2}-{\nu }_{1}}{1+{({\nu }_{2}-{\nu }_{1})}^{2}{\xi }^{2}/4}.\end{eqnarray} \tag{ 47 }$$

(B) In the second case, the variable ν oscillates in the interval

$$\begin{eqnarray}&&{\nu }_{3}\leqslant \nu \leqslant {\nu }_{4}.\end{eqnarray} \tag{ 48 }$$

Here again, a standard calculation yields

$$\begin{eqnarray}&&\nu ={\nu }_{3}+\displaystyle \frac{({\nu }_{4}-{\nu }_{3}){\mathrm{cn}}^{2}(\theta ,m)}{1+\tfrac{{\nu }_{4}-{\nu }_{3}}{{\nu }_{3}-{\nu }_{1}}{\mathrm{sn}}^{2}(\theta ,m)}\end{eqnarray} \tag{ 49 }$$

with the same definitions (42) and (43) of θ, m, and L. In this case we have $\nu (0)={\nu }_{4}$. In the soliton limit ${\nu }_{3}\to {\nu }_{2}$ ($m\to 1$) we get

$$\begin{eqnarray}&&\nu ={\nu }_{2}+\displaystyle \frac{{\nu }_{4}-{\nu }_{2}}{{\cosh }^{2}\theta +\tfrac{{\nu }_{4}-{\nu }_{2}}{{\nu }_{2}-{\nu }_{1}}{\sinh }^{2}\theta }.\end{eqnarray} \tag{ 50 }$$

This is a 'bright soliton' for the variable ν.

Again, the limit $m\to 0$ can be reached in two ways.

• (i)  
  If ${\nu }_{4}\to {\nu }_{3}$, then we obtain a small-amplitude harmonic wave

  $$\begin{eqnarray}&&\nu \cong {\nu }_{3}+\displaystyle \frac{1}{2}({\nu }_{4}-{\nu }_{3})\cos (k\xi ),\quad k=\sqrt{({\nu }_{3}-{\nu }_{1})({\nu }_{3}-{\nu }_{2})}.\end{eqnarray} \tag{ 51 }$$
• (ii)  
  If ν₂ = ν₁, then we obtain another nonlinear trigonometric solution,

  $$\begin{eqnarray}&&\nu ={\nu }_{3}+\displaystyle \frac{({\nu }_{4}-{\nu }_{3}){\cos }^{2}\theta }{1+\tfrac{{\nu }_{4}-{\nu }_{3}}{{\nu }_{3}-{\nu }_{1}}{\sin }^{2}\theta },\quad \theta =\sqrt{({\nu }_{3}-{\nu }_{1})({\nu }_{4}-{\nu }_{1})}\,\xi /2.\end{eqnarray} \tag{ 52 }$$

If we assume that ${\nu }_{4}-{\nu }_{3}\ll {\nu }_{4}-{\nu }_{1}$, then we reproduce the small-amplitude limit (51) with ν₂ = ν₁. On the other hand, in the limit ${\nu }_{3}\to {\nu }_{2}={\nu }_{1}$ we obtain the algebraic soliton solution:

$$\begin{eqnarray}&&\nu ={\nu }_{1}+\displaystyle \frac{{\nu }_{4}-{\nu }_{1}}{1+{({\nu }_{4}-{\nu }_{1})}^{2}{\xi }^{2}/4}.\end{eqnarray} \tag{ 53 }$$

For both cases (35), (36) we have the identities

$$\begin{eqnarray}&&m=\displaystyle \frac{({\nu }_{4}-{\nu }_{3})({\nu }_{2}-{\nu }_{1})}{({\nu }_{4}-{\nu }_{2})({\nu }_{3}-{\nu }_{1})}=\displaystyle \frac{({\lambda }_{4}^{2}-{\lambda }_{3}^{2})({\lambda }_{2}^{2}-{\lambda }_{1}^{2})}{({\lambda }_{4}^{2}-{\lambda }_{2}^{2})({\lambda }_{3}^{2}-{\lambda }_{1}^{2})}.\end{eqnarray} \tag{ 54 }$$

The importance of this form of periodic solutions of our equation is related with the fact that the parameters λ_j, connected with ν_i by the formulae (35), (36), can play the role of Riemann invariants in the Whitham theory of modulations.

In modulated waves the parameters ${\lambda }_{i}$ become slowly varying functions of the space and time variables and their evolution is governed by the Whitham modulation equations. Whitham showed in [10] that these equations can be obtained by averaging the conservation laws of the full nonlinear system over fast oscillations (whose wavelength L changes slowly along the total wave pattern). Generally speaking, in cases where the periodic solution is characterized by four parameters, this averaging procedure leads to a system of four equations of the type ${\nu }_{i,t}+{\sum }_{j}{v}_{{ij}}({\nu }_{1},{\nu }_{2},{\nu }_{3},{\nu }_{4}){\nu }_{j,x}=0$ with 16 entries in the 'velocity matrix' v_ij. However, for the case of completely integrable DNLS equation, this system of four equations reduces to a diagonal Riemann form for the Riemann invariants ${\lambda }_{i}$'s, similar to what occurs for the usual Riemann invariants of non-dispersive waves (see equation (12)). We shall derive the modulation Whitham equations by the method developed in [30, 31].

First of all, we notice that with the use of (22) and (37) it is easy to prove the identity

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left(\sqrt{P(\lambda )}\cdot \displaystyle \frac{G(\lambda )}{g(\lambda )}\right)-\displaystyle \frac{\partial }{\partial x}\left(\sqrt{P(\lambda )}\cdot \displaystyle \frac{B(\lambda )}{g(\lambda )}\right)=0,\end{eqnarray} \tag{ 55 }$$

where we have introduced under the derivative signs the constant on periodic solutions factor $\sqrt{P(\lambda )}$ to transform the identity (27) to the form

$$\begin{eqnarray*}&&{\left(\displaystyle \frac{f}{\sqrt{P(\lambda )}}\right)}^{2}-\displaystyle \frac{g}{\sqrt{P(\lambda )}}\cdot \displaystyle \frac{h}{\sqrt{P(\lambda )}}=1,\end{eqnarray*}$$

so that the right-hand side is independent of the variations of ${\lambda }_{i}$ in a modulated wave, hence the densities and fluxes in the conservation laws can change due to modulations only, as it should be, and any changes caused by λ-dependent normalization of the f, g, h-functions are excluded. We shall use the equation (55) as the generating function of the conservation laws of the DNLS equation: a series expansion in inverse powers of $\lambda$ gives an infinite number of conservation laws of this completely integrable system.

Substitution of equations (22) and (37) into (55) and its simple transformation gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{\sqrt{P(\lambda )}}{{\lambda }^{2}-\mu }\right)-\displaystyle \frac{\partial }{\partial x}\left[\sqrt{P(\lambda )}\left(2+\displaystyle \frac{{s}_{1}}{{\lambda }^{2}-\mu }\right)\right]=0,\end{eqnarray} \tag{ 56 }$$

Averaging of the density and of the flux in this expression over one wavelength L

$$\begin{eqnarray}&&L=\oint \displaystyle \frac{{\rm{d}}\mu }{4\sqrt{-P({\mu }^{1/2})}}\end{eqnarray} \tag{ 57 }$$

according to the rule

$$\begin{eqnarray*}&&\langle \{\ldots \}\rangle ={\int }_{0}^{L}\{\ldots \}\displaystyle \frac{{\rm{d}}x}{L}=\displaystyle \frac{1}{L}\oint \{\ldots \}\displaystyle \frac{{\rm{d}}x}{{\rm{d}}\mu }{\rm{d}}\mu =\displaystyle \frac{1}{L}\oint \{\ldots \}\displaystyle \frac{{\rm{d}}\mu }{4\sqrt{-P({\mu }^{1/2})}}\end{eqnarray*}$$

yields the generating function of the averaged conservation laws:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left[\displaystyle \frac{\sqrt{P(\lambda )}}{L}\oint \displaystyle \frac{{\rm{d}}\mu }{4({\lambda }^{2}-\mu )\sqrt{-P({\mu }^{1/2})}}\right]\\ &&\quad -\,\displaystyle \frac{\partial }{\partial x}\left[\sqrt{P(\lambda )}\left(2+\displaystyle \frac{{s}_{1}}{L}\oint \displaystyle \frac{{\rm{d}}\mu }{4({\lambda }^{2}-\mu )\sqrt{-P({\mu }^{1/2})}}\right)\right]=0.\end{eqnarray} \tag{ 58 }$$

The condition that in the limit $\lambda \to {\lambda }_{i}$ the singular terms cancel yields

$$\begin{eqnarray}&&\oint \displaystyle \frac{d\mu }{4({\lambda }_{i}^{2}-\mu )\sqrt{-P({\mu }^{1/2})}}\cdot \displaystyle \frac{\partial {\lambda }_{i}^{2}}{\partial t}\\ &&\quad -\,\left(2L+{s}_{1}\oint \displaystyle \frac{d\mu }{4({\lambda }_{i}-\mu )\sqrt{-P({\mu }^{1/2})}}\right)\cdot \displaystyle \frac{\partial {\lambda }_{i}^{2}}{\partial x}=0.\end{eqnarray} \tag{ 59 }$$

From the definition (57) of L one obtains

$$\begin{eqnarray*}&&\oint \displaystyle \frac{d\mu }{4({\lambda }_{i}^{2}-\mu )\sqrt{-P({\mu }^{1/2})}}=-2\,\displaystyle \frac{\partial L}{\partial {\lambda }_{i}^{2}},\end{eqnarray*}$$

which makes it possible to cast equation (59) in the form of a Whitham equation for the variables ${\lambda }_{i}$:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\lambda }_{i}}{\partial t}+{v}_{i}\displaystyle \frac{\partial {\lambda }_{i}}{\partial x}=0,\end{eqnarray} \tag{ 60 }$$

where the Whitham velocities v_i are given by

$$\begin{eqnarray}&&{v}_{i}=-{s}_{1}+\displaystyle \frac{L}{\partial L/\partial {\lambda }_{i}^{2}},\quad {\rm{for}}\quad i=1,2,3,4.\end{eqnarray} \tag{ 61 }$$

The values ${\lambda }_{i}$ of the spectral parameters are well-defined Riemann invariants of the Whitham system of modulation equations, however, they do not suit well enough to the problems with matching of modulated cnoidal waves and smooth dispersionless solutions. Therefore it is more convenient to define new set of Whitham invariants by using simple fact that any function of a single argument ${\lambda }_{i}$ is also a Riemann invariant. We define the new Riemann invariants by the formulae

$$\begin{eqnarray}&&{r}_{i}=-2{\lambda }_{i}^{2},\quad i=1,2,3,4.\end{eqnarray} \tag{ 62 }$$

They are negative and ordered according to ${r}_{1}\leqslant {r}_{2}\leqslant {r}_{3}\leqslant {r}_{4}\leqslant 0$ for ${\lambda }_{1}\leqslant {\lambda }_{2}\leqslant {\lambda }_{3}\leqslant {\lambda }_{4}\leqslant 0$. The parameters ν_i of the periodic solutions of the DNLS equation are expressed in terms of r_i as

$$\begin{eqnarray}{\nu }_{1} & = & {(-\sqrt{-{r}_{1}}+\sqrt{-{r}_{2}}+\sqrt{-{r}_{3}}+\sqrt{-{r}_{4}})}^{2}/2,\\ {\nu }_{2} & = & {(\sqrt{-{r}_{1}}-\sqrt{-{r}_{2}}+\sqrt{-{r}_{3}}+\sqrt{-{r}_{4}})}^{2}/2,\\ {\nu }_{3} & = & {(\sqrt{-{r}_{1}}+\sqrt{-{r}_{2}}-\sqrt{-{r}_{3}}+\sqrt{-{r}_{4}})}^{2}/2,\\ {\nu }_{4} & = & {(\sqrt{-{r}_{1}}+\sqrt{-{r}_{2}}+\sqrt{-{r}_{3}}-\sqrt{-{r}_{4}})}^{2}/2,\end{eqnarray} \tag{ 63 }$$

or

$$\begin{eqnarray}{\nu }_{1} & = & {(-\sqrt{-{r}_{1}}+\sqrt{-{r}_{2}}+\sqrt{-{r}_{3}}-\sqrt{-{r}_{4}})}^{2}/2,\\ {\nu }_{2} & = & {(\sqrt{-{r}_{1}}-\sqrt{-{r}_{2}}+\sqrt{-{r}_{3}}-\sqrt{-{r}_{4}})}^{2}/2,\\ {\nu }_{3} & = & {(\sqrt{-{r}_{1}}+\sqrt{-{r}_{2}}-\sqrt{-{r}_{3}}-\sqrt{-{r}_{4}})}^{2}/2,\\ {\nu }_{4} & = & {(\sqrt{-{r}_{1}}+\sqrt{-{r}_{2}}+\sqrt{-{r}_{3}}+\sqrt{-{r}_{4}})}^{2}/2.\end{eqnarray} \tag{ 64 }$$

The phase velocity and the wavelength are given by

$$\begin{eqnarray}&&V=\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{4}{r}_{i},\quad L=\displaystyle \frac{2K(m)}{\sqrt{({r}_{4}-{r}_{2})({r}_{3}-{r}_{1})}},\quad m=\displaystyle \frac{({r}_{4}-{r}_{3})({r}_{2}-{r}_{1})}{({r}_{4}-{r}_{2})({r}_{3}-{r}_{1})}.\end{eqnarray} \tag{ 65 }$$

The Whitham modulation equations read

$$\begin{eqnarray}&&\displaystyle \frac{\partial {r}_{i}}{\partial t}+{v}_{i}\displaystyle \frac{\partial {r}_{i}}{\partial x}=0,\quad i=1,2,3,4,\end{eqnarray} \tag{ 66 }$$

where the Whitham velocities v_i are given by

$$\begin{eqnarray}&&{v}_{i}=-{s}_{1}-\displaystyle \frac{L}{2\partial L/\partial {r}_{i}},\quad {\rm{for}}\quad i=1,2,3,4,\end{eqnarray} \tag{ 67 }$$

and substitution of L from (65) gives after simple calculation the following explicit expressions

$$\begin{eqnarray}{v}_{1} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{4}{r}_{i}-\displaystyle \frac{({r}_{4}-{r}_{1})({r}_{2}-{r}_{1})K(m)}{({r}_{4}-{r}_{1})K(m)-({r}_{4}-{r}_{2})E(m)},\\ {v}_{2} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{4}{r}_{i}+\displaystyle \frac{({r}_{3}-{r}_{2})({r}_{2}-{r}_{1})K(m)}{({r}_{3}-{r}_{2})K(m)-({r}_{3}-{r}_{1})E(m)},\\ {v}_{3} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{4}{r}_{i}-\displaystyle \frac{({r}_{4}-{r}_{3})({r}_{3}-{r}_{2})K(m)}{({r}_{3}-{r}_{2})K(m)-({r}_{4}-{r}_{2})E(m)},\\ {v}_{4} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{4}{r}_{i}+\displaystyle \frac{({r}_{4}-{r}_{3})({r}_{4}-{r}_{1})K(m)}{({r}_{4}-{r}_{1})K(m)-({r}_{3}-{r}_{1})E(m)},\end{eqnarray} \tag{ 68 }$$

where K(m) and E(m) are complete elliptic integrals of the first and second type, respectively.

In a modulated wave representing a dispersive shock wave, the Riemann invariants change slowly with x and t. The dispersive shock wave occupies a space interval at whose edges two of the Riemann invariants are equal to each other. The soliton edge corresponds to r₃ = r₂ (m = 1) and at this edge the Whitham velocities are given by

$$\begin{eqnarray}{v}_{1} & = & \displaystyle \frac{1}{2}(3{r}_{1}+{r}_{4}),\quad {v}_{4}=\displaystyle \frac{1}{2}({r}_{1}+3{r}_{4}),\\ {v}_{2} & = & {v}_{3}=\displaystyle \frac{1}{2}({r}_{1}+2{r}_{2}+{r}_{4}).\end{eqnarray} \tag{ 69 }$$

The opposite limit m = 0 can be obtained in two ways. If r₃ = r₄, then we get

$$\begin{eqnarray}{v}_{1} & = & \displaystyle \frac{1}{2}(3{r}_{1}+{r}_{2}),\quad {v}_{2}=\displaystyle \frac{1}{2}({r}_{1}+3{r}_{2}),\\ {v}_{3} & = & {v}_{4}=2{r}_{4}+\displaystyle \frac{{({r}_{2}-{r}_{1})}^{2}}{2({r}_{1}+{r}_{2}-2{r}_{4})},\end{eqnarray} \tag{ 70 }$$

and if r₂ = r₁, then

$$\begin{eqnarray}{v}_{1} & = & {v}_{2}=2{r}_{1}+\displaystyle \frac{{({r}_{4}-r)}^{2}}{2({r}_{3}+{r}_{4}-2{r}_{1})},\\ {v}_{3} & = & \displaystyle \frac{1}{2}(3{r}_{3}+{r}_{4}),\quad {v}_{4}=\displaystyle \frac{1}{2}({r}_{3}+3{r}_{4}).\end{eqnarray} \tag{ 71 }$$

From these equations it is clear that at the edges of the oscillatory zone the Whitham equation for two Riemann invariants coincide with those for the dispersionless equations, that is the oscillatory zone can match at its edges with smooth solutions of the dispersionless equations.

Now we are ready to discuss the key elements from which consists any wave structure evolving from an initial discontinuity.

Our aim in this paper is to develop the method of derivation of the asymptotic solution of the DNLS evolution problem for a discontinuous step-like initial conditions

$$\begin{eqnarray}\rho (x,0)=\left\{\begin{array}{ll}{\rho }^{L} & \mathrm{for}\ x\,\lt \,0\\ {\rho }^{R} & \mathrm{for}\ x\,\gt \,0\end{array}\right.\qquad u(x,0)=\left\{\begin{array}{ll}{u}^{L} & \mathrm{for}\ x\,\lt \,0\\ {u}^{R} & \mathrm{for}\ x\,\gt \,0\end{array}\right..\end{eqnarray} \tag{ 72 }$$

As we shall see, evolution of this step-like pulse leads to formation of quite complex wave structures consisting of several simpler elements of simple wave type with only one Riemann invariant changing. Therefore we shall first describe these elements in the present section.

7.1. Rarefaction waves

For smooth enough dependence of wave parameters on x and t, we can neglect the dispersion effects and use the dispersionless equations derived in section 3. First of all, we notice that the system (12) has a trivial solution for which r₊  = const and r₋ = const. We shall call such a solution a 'plateau' because it corresponds to a uniform flow with constant density and flow velocity given by (15).

The initial conditions (72) do not contain any parameters with dimension of time or length. Therefore solutions of equation (12) can depend on the self-similar variable ζ = x/t only, that is ${r}_{\pm }={r}_{\pm }(\zeta )$, and then this system reduces to

$$\begin{eqnarray}&&({v}_{+}-\zeta )\displaystyle \frac{{{dr}}_{+}}{d\zeta }=0,\quad ({v}_{-}-\zeta )\displaystyle \frac{{{dr}}_{-}}{d\zeta }=0.\end{eqnarray} \tag{ 73 }$$

Evidently, these equations have solutions with one of the Riemann invariants constant and the other one changes in such a way, that the corresponding velocity equals to ζ = x/t. To be definite, let us consider the solution

$$\begin{eqnarray}&&{r}_{+}={r}_{+}^{0}=\mathrm{const},\qquad {v}_{-}=\displaystyle \frac{1}{2}{r}_{+}^{0}+\displaystyle \frac{3}{2}{r}_{-}=\zeta =\displaystyle \frac{x}{t}.\end{eqnarray} \tag{ 74 }$$

Consequently, ${r}_{-}$ depends on x/t as

$$\begin{eqnarray}&&{r}_{-}(x,t)=-\displaystyle \frac{1}{3}{r}_{+}^{0}+\displaystyle \frac{2}{3}\cdot \displaystyle \frac{x}{t},\end{eqnarray} \tag{ 75 }$$

and according to equation (15) the physical variables are given by

$$\begin{eqnarray}{\rho }_{\pm }(x,t) & = & \displaystyle \frac{1}{2}{\left(\sqrt{-{r}_{+}^{0}}\pm \sqrt{\displaystyle \frac{1}{3}{r}_{+}^{0}-\displaystyle \frac{2}{3}\cdot \displaystyle \frac{x}{t}}\right)}^{2},\\ {u}_{\pm }(x,t) & = & \pm 2\sqrt{-{r}_{+}^{0}\left(\displaystyle \frac{1}{3}{r}_{+}^{0}-\displaystyle \frac{2}{3}\cdot \displaystyle \frac{x}{t}\right)}.\end{eqnarray} \tag{ 76 }$$

Here the single solution (74) of equations written in Riemann form yields two solutions (76) in physical variables which we distinguish by the indices ±. These rarefaction waves match to the plateau solutions at their left and right edges. At both edges the invariant r₊ has the same value ${r}_{+}={r}_{+}^{0}$ whereas we have ${r}_{-}={r}_{-}^{L}$ at the left boundary and ${r}_{-}={r}_{-}^{R}$ at the right boundary. Correspondingly, the above two solutions match to the values of the density

$$\begin{eqnarray}(i)\quad {\rho }_{+}^{L} & = & \displaystyle \frac{1}{2}{(\sqrt{-{r}_{+}^{0}}+\sqrt{-{r}_{-}^{L}})}^{2},\quad {\rho }_{+}^{R}=\displaystyle \frac{1}{2}{(\sqrt{-{r}_{+}^{0}}+\sqrt{-{r}_{-}^{R}})}^{2},\\ ({ii})\quad {\rho }_{-}^{L} & = & \displaystyle \frac{1}{2}{(\sqrt{-{r}_{+}^{0}}-\sqrt{-{r}_{-}^{L}})}^{2},\quad {\rho }_{-}^{R}=\displaystyle \frac{1}{2}{(\sqrt{-{r}_{+}^{0}}-\sqrt{-{r}_{-}^{R}})}^{2},\end{eqnarray} \tag{ 77 }$$

and similar formulae can be written for the flow velocities ${u}_{\pm }^{L,R}$. The edge points propagate with velocities

$$\begin{eqnarray}&&{s}^{L}={v}_{-}({r}_{+}^{0},{r}_{-}^{L})=\displaystyle \frac{1}{2}{r}_{+}^{0}+\displaystyle \frac{3}{2}{r}_{-}^{L},\quad {s}^{R}={v}_{-}({r}_{+}^{0},{r}_{-}^{R})=\displaystyle \frac{1}{2}{r}_{+}^{0}+\displaystyle \frac{3}{2}{r}_{-}^{R}.\end{eqnarray} \tag{ 78 }$$

Since ${r}_{-}^{R}\lt {r}_{+}^{R}={r}_{+}^{0}$, we always have ${s}^{R}\lt {r}_{+}^{0}$.

In a similar way we obtain the second solution

$$\begin{eqnarray}&&{r}_{+}=-\displaystyle \frac{1}{3}{r}_{-}^{0}+\displaystyle \frac{2}{3}\cdot \displaystyle \frac{x}{t},\quad {r}_{-}={r}_{-}^{0}=\mathrm{const},\end{eqnarray} \tag{ 79 }$$

hence

$$\begin{eqnarray}{\rho }_{\pm }(x,t) & = & \displaystyle \frac{1}{2}{\left(\sqrt{\displaystyle \frac{1}{3}{r}_{-}^{0}-\displaystyle \frac{2}{3}\cdot \displaystyle \frac{x}{t}}\pm \sqrt{-{r}_{-}^{0}}\right)}^{2},\\ {u}_{\pm }(x,t) & = & \pm 2\sqrt{-{r}_{-}^{0}\left(\displaystyle \frac{1}{3}{r}_{-}^{0}-\displaystyle \frac{2}{3}\cdot \displaystyle \frac{x}{t}\right)}.\end{eqnarray} \tag{ 80 }$$

In this case we have

$$\begin{eqnarray}&&{r}_{-}^{0}\leqslant {s}^{L}=\displaystyle \frac{3}{2}{r}_{+}^{L}+\displaystyle \frac{1}{2}{r}_{-}^{0}\lt {s}^{R}=\displaystyle \frac{3}{2}{r}_{+}^{R}+\displaystyle \frac{1}{2}{r}_{-}^{0}\lt \displaystyle \frac{1}{2}{r}_{-}^{0}.\end{eqnarray} \tag{ 81 }$$

Diagrams of the Riemann invariants for these rarefaction wave solutions are shown in figure 1: the case (a) corresponds to equations (74), (75) and the case (b) to equation (79). Corresponding plots of densities are demonstrated in figure 2 by thick lines together with plateau distributions at the edges of the rarefaction waves. Dashed thick lines show both branches of the solutions (76) and (80). It is worth noticing that the edge velocities are determined by the Riemann invariants only and do not depend on the choice of the branch into which the Riemann invariants are mapped.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is useful to give another graphic representation of the rarefaction waves. From definition (13) of we Riemann invariants we find that they are constant along parabolas

$$\begin{eqnarray}&&\rho =-\displaystyle \frac{1}{2r}{\left(\displaystyle \frac{u}{2}-r\right)}^{2}\end{eqnarray} \tag{ 82 }$$

in the (u, ρ)-plane, where r is the value of the corresponding Riemann invariant. If a rarefaction wave corresponds to ${r}_{+}=\mathrm{const}$, then both its left and right points must lie on the same parabola shown in figure 3(a) by a blue line. These points can be represented as crossing points of this blue parabola with other two parabolas that represent curves with constant ${r}_{-}^{L}$ and ${r}_{-}^{R}$ and are shown by red lines. We have two pairs of 'left' and 'right' points and obtain, consequently, two types of rarefaction waves described by the diagram figure 1(a). These transitions ${L}_{a}\to {R}_{a}$ and ${L}_{b}\to {R}_{b}$ correspond to different signs in the formulas (76). As we see, both transitions give the growth of ρ with increase of ζ in agreement with the plots in figure 2(a). In a similar way, the situations corresponding to the diagram figure 1(b) with constant Riemann invariant r₋ are represented by the parabolas shown in figure 3(b). Now transitions from the 'left' points to the 'right' ones give the growth of ρ in one case and its decrease in another case, as it is shown in figure 2(b). It is important to notice that according to equation (76) these transitions connect the points with the same signs of u, that is they do not intersect the ordinate axis separating the monotonicity regions. Thus, these rarefaction waves connect the states belonging to the same regions of monotonicity of the Riemann invariants. In the next sections we shall generalize this graphical representation to other wave structures what will be quite helpful in classification of possible wave structures evolving from initial discontinuities.

Zoom In Zoom Out Reset image size

Both solutions for ρ₋ can describe flow of liquid into vacuum—in case (76) from left to right and in case (80) from right to left. It is worth noticing a curious particular solution for u^L = 0, when ${r}_{+}^{0}=0$, r₋ = 2x/(3t) and we get ρ = −x/(3t), u = 0. It is easy to see that dispersionless equation (9) admit such a solution.

Considered here wave structures satisfy the conditions (a) ${r}_{+}^{L}={r}_{+}^{R}$, ${r}_{-}^{L}\lt {r}_{-}^{R}$ or (b) ${r}_{+}^{L}\lt {r}_{+}^{R}$, ${r}_{-}^{L}={r}_{-}^{R}$. It is natural to ask, what happens if we have the initial conditions satisfying opposite inequalities, and to answer this question we have to consider the DSW structures.

7.2. Cnoidal dispersive shock waves

The other two possible solutions of equation (73) are sketched in figure 4, and they satisfy the boundary conditions (a) ${r}_{+}^{L}={r}_{+}^{R}$, ${r}_{-}^{L}\gt {r}_{-}^{R}$ or (b) ${r}_{+}^{L}\gt {r}_{+}^{R}$, ${r}_{-}^{L}={r}_{-}^{R}$. In the dispersionless approximation these multi-valued solutions are nonphysical. However, following to Gurevich and Pitaevskii [9], we can give them clear physical sense by understanding r_i as four Riemann invariants of the Whitham system that describe evolution of a modulated nonlinear periodic wave. Naturally, now r_i are the self-similar solutions of the Whitham equation (66), that is of the equations

$$\begin{eqnarray}&&({v}_{i}-\zeta )\displaystyle \frac{{{dr}}_{i}}{d\zeta }=0,\qquad i=1,2,3,4,\end{eqnarray} \tag{ 83 }$$

which are obvious generalizations of (74):

$$\begin{eqnarray}(a)\quad {r}_{1} & = & {r}_{-}^{R},\quad {r}_{3}={r}_{+}^{R},\quad {r}_{4}={r}_{-}^{L},\quad {v}_{2}({r}_{-}^{R},{r}_{2},{r}_{-}^{L},{r}_{+}^{L})=\zeta ,\\ (b)\quad {r}_{1} & = & {r}_{-}^{R},\quad {r}_{2}={r}_{+}^{R},\quad {r}_{4}={r}_{+}^{L},\quad {v}_{3}({r}_{-}^{R},{r}_{+}^{R},{r}_{3},{r}_{+}^{L})=\zeta ,\end{eqnarray} \tag{ 84 }$$

where the last relations determine implicitly dependence of r₂ and r₃, correspondingly, on ζ. Sketches of these solutions are shown in figure 4. Velocities of the edges of the oscillatory zone whose envelopes are described by the solutions of the Whitham equations are given by

$$\begin{eqnarray}(a)\quad {s}^{L} & = & {v}_{2}({r}_{-}^{R},{r}_{-}^{R},{r}_{-}^{L},{r}_{+}^{L})=2{r}_{-}^{R}+\displaystyle \frac{{({r}_{+}^{L}-{r}_{-}^{L})}^{2}}{2({r}_{+}^{L}+{r}_{-}^{L}-2{r}_{-}^{R})},\\ {s}^{R} & = & {v}_{2}({r}_{-}^{R},{r}_{-}^{L},{r}_{-}^{L},{r}_{+}^{L})=\displaystyle \frac{1}{2}({r}_{-}^{R}+2{r}_{-}^{L}+{r}_{+}^{L}),\\ (b)\quad {s}^{L} & = & {v}_{3}({r}_{-}^{R},{r}_{+}^{R},{r}_{+}^{R},{r}_{+}^{L})=\displaystyle \frac{1}{2}({r}_{-}^{R}+2{r}_{+}^{R}+{r}_{+}^{L}),\\ {s}^{R} & = & {v}_{2}({r}_{-}^{R},{r}_{-}^{R},{r}_{-}^{L},{r}_{+}^{L})=2{r}_{-}^{R}+\displaystyle \frac{{({r}_{+}^{L}-{r}_{-}^{L})}^{2}}{2({r}_{+}^{L}+{r}_{-}^{L}-2{r}_{-}^{R})},\end{eqnarray} \tag{ 85 }$$

correspondingly.

Zoom In Zoom Out Reset image size

If we substitute the solutions (84) into formulae (63) and (63), then we determine the dependence of the parameters ν_i on ζ. There are two possibilities shown in figure 5: the diagram figure 4(a) is mapped by both sets of formulae (63) and (64) into the type figure 5(i), whereas the diagram figure 4(b) is mapped by the formulae (63) into the type figure 5(ii) and by the formulae (64) into the type figure 5(i).

Zoom In Zoom Out Reset image size

The solutions obtained here are interpreted as formation of cnoidal dispersive shock waves evolving from initial discontinuities with such a type of the boundary conditions. Indeed, equation (41) upon substitution of obtained ν_i yields the plots shown in figure 6(a)I,II and figure 6(b)II, whereas equation (49) yields the plot figure 6(b)I. We summarize the appearing possibilities in the following list:

• Figure 4(a) $\mathop{\longrightarrow }\limits^{(63)}$ Figure 5(i) ⟶ figure 6(a) plot I
• Figure 4(a) $\mathop{\longrightarrow }\limits^{(64)}$ Figure 5(i) ⟶ figure 6(a) plot II
• Figure 4(b) $\mathop{\longrightarrow }\limits^{(63)}$ Figure 5(i) ⟶ figure 6(b) plot II
• Figure 4(b) $\mathop{\longrightarrow }\limits^{(64)}$ Figure 5(ii) ⟶ figure 6(b) plot I

As we see, each solution of Whitham equations expressed in terms of Riemann invariants is mapped into two different DSW structures which satisfy the boundary conditions compatible with given values of the Riemann invariants. Such a behavior is typical for non-convex dispersive hydrodynamics and has already been discussed in simpler situation of mKdV (Gardner) equation in [20].

Zoom In Zoom Out Reset image size

This two-valued connection of Riemann invariants with solutions in terms of physical variables is similar to the situation described above for the rarefaction waves: the diagram figure 1(a) yields two decreasing with ζ density distributions shown in figure 2(a) whereas the diagram figure 1(b) yields decreasing and increasing distributions shown in figure 2(b). These two types of wave structures will serve us as building blocks appearing in evolution of arbitrary initial discontinuity. It is clear that these cnoidal DSWs are described by the same diagrams of figure 3 as the rarefaction waves, but with inverted 'left' and 'right' points. Hence, the cnoidal DSWs still connect the states belonging to the same regions of monotonicity of the dispersionless Riemann invariants. But there must be waves that connects the states at opposite sides of the ρ-axis u = 0 in the (u, ρ)-plane and they also appear as elementary building blocks which are described by the self-similar solutions of the Whitham equations. We shall turn to this type of waves in the next subsection.

7.3. Trigonometric (contact) dispersive shock waves

At first we shall consider the situation in which the Riemann invariants have equal values at both edges of the shock, i.e., when ${r}_{-}^{L}={r}_{-}^{R}$, ${r}_{+}^{L}={r}_{+}^{R}$. In this case we obtain a new type of wave structure which we shall call a contact dispersive shock wave since it has some similarity with contact discontinuities in the theory of viscous shock waves (see, e.g., [36]). For this situation, the parabolas corresponding to ${r}_{-}^{L}=\mathrm{const}$ and ${r}_{-}^{R}=\mathrm{const}$ in figure 3(a) coincide with each other and cnoidal DSWs disappear. Instead, there appears the path connecting the identical left and right states labeled by the crossing points of two parabolas as is shown in figure 7. Such waves can arise only if the boundary points are located on the opposite sides of the line u = 0, i.e. in the different regions of monotonicity.

Zoom In Zoom Out Reset image size

Along the arc of the parabola connecting the points P₁ and P₂ the two biggest Riemann invariants must be equal to each other, r₃ = r₄, and at the left soliton edge they must equal to their boundary value ${r}_{3}={r}_{4}={r}_{2}={r}_{+}^{L}={r}_{+}^{R}$. Hence, we arrive at the diagram of the Riemann invariants shown in figure 8. Along this solution we have m = 0 and the solutions of the Whitham equations is determined by the formula

$$\begin{eqnarray}&&{v}_{3}={v}_{4}=2{r}_{4}+\displaystyle \frac{{({r}_{+}^{L}-{r}_{-}^{L})}^{2}}{2({r}_{+}^{L}+{r}_{-}^{L}-2{r}_{4})}=\zeta \end{eqnarray} \tag{ 86 }$$

from which we obtain

$$\begin{eqnarray}&&{r}_{3}={r}_{4}=\displaystyle \frac{1}{4}[{r}_{+}^{L}+{r}_{-}^{L}+\zeta +\sqrt{{({r}_{+}^{L}+{r}_{-}^{L}-\zeta )}^{2}+2{({r}_{+}^{L}-{r}_{-}^{L})}^{2}}].\end{eqnarray} \tag{ 87 }$$

At the left soliton edge we have ${r}_{4}={r}_{+}^{L}$ and at the right small-amplitude edge r₄ = 0. Therefore equation (70) yields velocities of the edges:

$$\begin{eqnarray}&&{s}^{L}=\displaystyle \frac{3}{2}{r}_{+}^{L}+\displaystyle \frac{1}{2}{r}_{-}^{L},\qquad {s}^{R}=\displaystyle \frac{{({r}_{+}^{L}-{r}_{-}^{L})}^{2}}{2({r}_{+}^{L}+{r}_{-}^{L})}.\end{eqnarray} \tag{ 88 }$$

The sign of the square root in equation (87) is chosen in such a way that this formula gives ${r}_{4}={r}_{+}^{L}$ at the left edge with $\zeta =(3{r}_{+}^{L}+{r}_{-}^{L})/2$.

Zoom In Zoom Out Reset image size

As one can see from equation (63), in this case ν₃ = ν₄ and ν oscillates in the interval ${\nu }_{1}\leqslant \nu \leqslant {\nu }_{2}$. Then equation (46) yields the plot shown in figure 9(a) with dark algebraic solitons at the left soliton edge. In case of equation (64) we have ν₁ = ν₂, hence ν oscillates in the interval ${\nu }_{3}\leqslant \nu \leqslant {\nu }_{4}$, and equation (52) yields the plot figure 9(b) with bright algebraic solitons at the soliton edge. Again the same solution of the Whitham equations represented by a single diagram figure 8 is mapped into two different wave structures.

Zoom In Zoom Out Reset image size

7.4. Combined shocks

Now we turn to the last elementary wave structures connecting two plateau states. They can also be symbolized by single parabolic arcs between two points in the (u, ρ)-plane. This type of paths is illustrated in figure 10 and obviously it is a generalization of the preceding structure. In this case, the boundary points are also located in different monotonicity regions. One of the Riemann invariants still remains constant (${r}_{-}^{L}={r}_{-}^{R}$), however, the boundary values of the other Riemann invariant are different: we have ${r}_{+}^{L}\lt {r}_{+}^{R}$ in case (a) and ${r}_{+}^{L}\gt {r}_{+}^{R}$ in case (b). The corresponding diagrams of Riemann invariants are shown in figure 11. As we see, in case (a) the oscillating region located between two plateaus consists of two subregions—one with four different Riemann invariants, what corresponds to a cnoidal DSW, and another one with r₃ = r₄, what corresponds to a trigonometric DSW, and there is no any plateau between them. Thus, this diagrams leads to a combined wave structure of 'glued' cnoidal and trigonometric DSWs. This structure is illustrated in figure 12(a). At the soliton edge the cnoidal DSW matches with the left plateau and the edge with m = 0 it degenerates into the trigonometric shock. Velocities of the edge points are equal to

$$\begin{eqnarray}{s}_{1}^{L} & = & {v}_{3}({r}_{-}^{L},{r}_{+}^{R},{r}_{+}^{R},{r}_{+}^{L})=\displaystyle \frac{1}{2}({r}_{-}^{L}+2{r}_{+}^{R}+{r}_{+}^{L}),\\ {s}_{2}^{L} & = & {v}_{3}({r}_{-}^{L},{r}_{+}^{R},{r}_{+}^{L},{r}_{+}^{L})=2{r}_{+}^{L}+\displaystyle \frac{{({r}_{+}^{R}-{r}_{-}^{R})}^{2}}{2({r}_{+}^{R}+{r}_{-}^{R}-2{r}_{+}^{L})},\\ {s}^{R} & = & {v}_{3}({r}_{-}^{L},{r}_{+}^{R},0,0)=\displaystyle \frac{{({r}_{+}^{R}-{r}_{-}^{R})}^{2}}{2({r}_{+}^{R}+{r}_{-}^{R})}.\end{eqnarray} \tag{ 89 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In a similar way, in case (b) we have a single trigonometric DSW region glued with a rarefaction wave, as is shown in figure 12(b). In this case the edge velocities are given by

$$\begin{eqnarray}{s}_{1}^{L} & = & \displaystyle \frac{1}{2}(3{r}_{+}^{L}+{r}_{-}^{L}),\\ {s}_{2}^{L} & = & \displaystyle \frac{1}{2}(3{r}_{+}^{R}+{r}_{-}^{R}),\\ {s}^{R} & = & {v}_{3}({r}_{-}^{L},{r}_{+}^{R},0,0)=\displaystyle \frac{{({r}_{+}^{R}-{r}_{-}^{R})}^{2}}{2({r}_{+}^{R}+{r}_{-}^{R})}.\end{eqnarray} \tag{ 90 }$$

In both cases, the oscillatory wave is described by the formula (49) or its limit (52) with oscillations of ν in the interval ${\nu }_{3}\leqslant \nu \leqslant {\nu }_{4}$.

Now, after description of all elementary wave structures arising in evolution of discontinuities in the DNLS equation theory, we are in position to formulate the main principles of classification of all possible wave structures.

Classification of possible structures is very simple in the KdV equation case when any discontinuity evolves into either rarefaction wave, or cnoidal DSW [9]. It becomes more complicated in the NLS equation case [12] and similar situations as, e.g., for the Kaup-Boussinesq equation [16, 17], where the list consists of eight or ten structures, correspondingly, which can be found after simple enough inspection of available possibilities which are studied one by one. However, the situation changes drastically when we turn to non-convex dispersive hydrodynamics: even in the case of unidirectional Gardner (mKdV) equation we get eight different patterns (instead of two in KdV case) due to appearance of new elements (kinks or trigonometric and combined dispersive shocks), but these patterns can be labeled by two parameters only and therefore these possibilities can be charted on a two-dimensional diagram. In our present case the initial discontinuity (72) is parameterized by four parameters u^L, ρ^L, u^R, ρ^R, hence the number of possible wave patterns considerably increases and it is impossible to present them in a two-dimensional chart. Therefore it seems more effective to formulate the principles according to which one can predict the wave pattern evolving from a discontinuity with given parameters. Similar method was used [34, 35] in classification of wave patterns evolving from initial discontinuities according to the generalized NLS equation and the Landau–Lifshitz equation for easy-plane magnetics or polarization waves in two-component Bose–Einstein condensate.

We begin with the consideration of the classification problem from the case when both boundary points lie on one side of the axis u = 0 separating two monotonicity regions in the $(u,\rho )$-plane. At first we shall consider situation when the boundary points lie in the left monotonicity region u < 0. We show in figure 13(a) the two parabolas corresponding to the constant dispersionless Riemann invariants ${r}_{\pm }^{L}$ related with the left boundary state. Evidently, they cross at some point $L({u}^{L},{\rho }^{L})$ representing the left boundary. These two parabolas cut the left monotonicity region into six domains labeled by the symbols A, B, ..., F. Depending on the domain, in which the point R with coordinates $({u}^{R},{\rho }^{R})$, representing the right boundary condition, is located, one gets one of the six following possible orderings of the left and right Riemann invariants:

$$\begin{eqnarray}{\rm{A}}: & & \quad {\lambda }_{-}^{R}\lt {\lambda }_{+}^{R}\lt {\lambda }_{-}^{L}\lt {\lambda }_{+}^{L},\\ {\rm{B}}: & & \quad {\lambda }_{-}^{R}\lt {\lambda }_{-}^{L}\lt {\lambda }_{+}^{R}\lt {\lambda }_{+}^{L},\\ {\rm{C}}: & & \quad {\lambda }_{-}^{L}\lt {\lambda }_{-}^{R}\lt {\lambda }_{+}^{R}\lt {\lambda }_{+}^{L},\\ {\rm{D}}: & & \quad {\lambda }_{-}^{R}\lt {\lambda }_{-}^{L}\lt {\lambda }_{+}^{L}\lt {\lambda }_{+}^{R},\\ {\rm{E}}: & & \quad {\lambda }_{-}^{L}\lt {\lambda }_{-}^{R}\lt {\lambda }_{+}^{L}\lt {\lambda }_{+}^{R},\\ {\rm{F}}: & & \quad {\lambda }_{-}^{L}\lt {\lambda }_{+}^{L}\lt {\lambda }_{-}^{R}\lt {\lambda }_{+}^{R}.\end{eqnarray} \tag{ 91 }$$

All these six domains and corresponding orderings yield six possible wave structures evolving from initial discontinuities. Let us consider briefly each of them.

Zoom In Zoom Out Reset image size

In case (A) two rarefaction waves are separated by an empty region. Evolution of Riemann invariants and sketch of wave structure are shown in figure 14(A).

Zoom In Zoom Out Reset image size

In case (B) two rarefaction waves are connected by a plateau whose parameters are determined by the dispersionless Riemann invariants ${r}_{\pm }^{P}$ equal to ${r}_{-}^{P}={r}_{-}^{R}$ and ${r}_{+}^{P}={r}_{+}^{L}$. Here left and right 'fluids' flow away from each other with small enough relative velocity and rarefaction waves are able now to provide enough flux to create a plateau in the region between them (see figure 14(B)).

In case (C) we obtain a dispersive shock wave on the left side of the structure, a rarefaction wave on its right side and a plateau in between (see figure 14(C)).

In case (D) we get the same situation as in the case (C), but now the dispersive shock wave and rarefaction wave exchange their places (see figure 14(D)).

In case (E) two DSWs are produced with a plateau between them. Here we have a collision of left and right fluids (see figure 14(E)).

In case (F) the plateau observed in the case (E) disappears. It is replaced by a nonlinear wave which can be presented as a non-modulated cnoidal wave (see figure 14(F)).

The possible structures for this part of the (u, ρ)-plane coincide qualitatively with the patterns found in similar classification problem for the nonlinear Schrödinger equation [12] and this case was already studied in [33].

If we turn to consideration of the classification problem for the case when both boundary points lie to the right of the line u = 0, then we get the diagram in the $(u,\rho )$-plane shown in figure 13(b). We see that the parabolas cut again this right monotonicity region into six domains. For this case the Riemann invariants can have the same orderings (91) as in the previous case. Depending on the location of the right boundary point in a certain domain, the corresponding wave structure will be formed. Qualitatively these structures coincide with those for the previous case.

At last, we have to study the situation when the boundary points lie on different sides of the line u = 0, that is in different monotonicity regions. As we have seen in the previous section, in this case new complex structures, namely, combined shocks, appear. It is easy to see that if the left boundary corresponds to the point in the left monotonicity region, then we get again six wave patters, and if it correspond to the point in the right monotonicity region, we get six more patterns, twelve in total. In principle, they can be considered as generalizations of those shown in figure 14 with simple elements (rarefaction waves and cnoidal DSWs) replaced by combined shocks. Instead of listing all possible patterns, we shall formulate the general principles of their construction and illustrate them by a typical example. This will provide the method by which one can predict the wave pattern evolving from any given initial discontinuity.

For given boundary parameters, we can construct the parabolas corresponding to constant Riemann invariants ${r}_{\pm }^{L,R}$: each left or right pair of these parabolas crosses at the point L or R representing the left or right boundary state's plateau. Our task is to construct the path joining these two points, then this path will represent the arising wave structure. We already know the answer for the case when the left and right points lie on the same parabola, see, e.g., figure 7. If this is not the case and the right point R lies, say, below the parabola ${r}_{-}^{L}=\mathrm{const}$, see figure 15(a), then we can reach R by means of more complicated path consisting of two arcs of parabolas LP and PR joined at the point P. Evidently, this point P represents the plateau between two waves represented by the arcs. At the same time, each arc corresponds to a wave structure discussed in the preceding section. Having constructed a path from the left boundary point to the right one, it is easy to draw the corresponding diagram of Riemann invariants. To construct the wave structure, we use the formulae connecting the zeroes ν_i of the resolvent with the Riemann invariants r_i and expressions for the solutions parameterized by ν_i. This solves the problem of construction of the wave structure evolving from the initial discontinuity with given boundary conditions. In fact, there are two paths with a single intersection point that join the left and right boundary points and we choose the physically relevant path by imposing the condition that velocities of edges of all regions must increase from left to right.

Zoom In Zoom Out Reset image size

For example, let us consider the case ρ^L = 0.45, u^L = −4, ρ^R = 1.6, u^R = 1.5 which corresponds to figure 15(a) with the transition $L\to R\in B(u\gt 0)$. In this case ${r}_{-}^{L}=-3.87$, ${r}_{+}^{L}=-1.035$, ${r}_{-}^{R}=-1.25$, ${r}_{+}^{R}=-0.45$ and we see that the arc PB of the parabola with ${r}_{-}^{R}=\mathrm{const}$ in the above transition crosses the axis u = 0 as is illustrated in figure 15(a). Thus, we arrive at the diagram of Riemann invariants shown in figure 15(b). Consequently, at the left edge we have a standard rarefaction wave (the arc LP does not cross the axis u = 0) and at the right edge the combination of a trigonometric shock with a rarefaction wave. Between these waves we get a plateau characterized by the Riemann invariants ${r}_{-}^{P}={r}_{-}^{R}$ and ${r}_{+}^{P}={r}_{+}^{L}$. This plateau is represented by a single point P in figure 15(a). The rarefaction waves are described by the formulas (76) (left wave) and (80) (right wave) with 'minus' sight chosen in them. The profile of the oscillatory wave structure can be obtained by substitution of the solution of the Whitham equations

$$\begin{eqnarray*}&&{r}_{1}={r}_{-}^{R},\quad {r}_{2}={r}_{+}^{R},\quad {v}_{3}({r}_{-}^{R},{r}_{+}^{R},{r}_{3},{r}_{3})={v}_{4}({r}_{-}^{R},{r}_{+}^{R},{r}_{3},{r}_{3})=\zeta ,\end{eqnarray*}$$

into equation (52) with ν_i given by equation (64). The velocities of the edge points are equal to

$$\begin{eqnarray*}{s}_{1}^{L} & = & \displaystyle \frac{1}{2}({r}_{+}^{L}+3{r}_{-}^{L}),\quad {s}_{2}^{L}=\displaystyle \frac{1}{2}({r}_{+}^{L}+3{r}_{-}^{R}),\\ {s}_{1}^{R} & = & \displaystyle \frac{1}{2}(3{r}_{+}^{L}+{r}_{-}^{R}),\quad {s}_{2}^{R}=\displaystyle \frac{1}{2}(3{r}_{+}^{R}+{r}_{-}^{R}),\quad {s}_{3}^{R}=\displaystyle \frac{{({r}_{+}^{R}-{r}_{-}^{R})}^{2}}{2({r}_{+}^{R}+{r}_{-}^{R})}.\end{eqnarray*}$$

The resulting wave pattern is shown in figure 16. It is easy to see, that it represents a deformation of the plot figure 14(B): due to crossing the axis u = 0 the right rarefaction waves acquires a tail in the form of trigonometric DSW. It should be stressed that appearance of such a tail is impossible in the theory of dispersive shock waves in the NLS equation case.

Zoom In Zoom Out Reset image size

In a similar way we can construct all twelve possible wave patterns for this type on the boundary conditions.

In this paper, we have developed the Whitham method of modulations for evolution of waves governed by the DNLS equation. The Riemann problem of evolution of an initial discontinuity is solved for this specific case of non-convex dispersive hydrodynamics. It is found that the set of possible wave structures is much richer than in the convex case (as, e.g., in the NLS equation theory) and includes, as structural elements, trigonometric shock combined with rarefaction waves or cnoidal dispersive shocks. Evolution of these trigonometric shocks is described by the degenerate limits of the Whitham modulation equations. In the resulting scheme, one solution of the Whitham equations corresponds to two different wave patterns, and this correspondence is provided by a two-valued mapping of Riemann invariants to physical modulation parameters. Thus, the algebraic resolvents introduced in [29] for effectivization of periodic solutions of integrable equations turns out to be crucially important also for establishing the relations between Riemann invariants and modulation parameters of periodic solutions. To determine the pattern evolving from given discontinuity, we have developed a graphical method which is quite flexible and was also applied to other systems with non-convex hydrodynamics—generalized NLS equation for propagation of light pulses in optical fibers [34] and Landau–Lifshitz equation for dynamics of magnetics with uniaxial easy-plane anisotropy [35]. Comparison of analytic theory with numerical solutions of the corresponding evolution equations demonstrates very good accuracy of the Whitham modulation approach even for not very long time of evolution.

The developed theory opens the door to two possible future developments. First, the principles of classification of the wave patterns are based on the study of the dispersionless limit which is not related with the integrability property of the evolution equation under consideration. Therefore it can be applied to the non-integrable equation with less detailed study of dispersive shock waves evolved from the initial discontinuity. Second, the small dissipative perturbation can change classification of the wave patterns, as it was found for the mKdV-Burgers case in [22]. Therefore the presented here theory can be generalized on the case of perturbed integrable equation with the use of the corresponding generalization of the Whitham theory obtained in [37]. These two possible generalizations can find a number of applications to concrete physical problems.

This work was partially supported by RFBR grant 16-01-00398. I am grateful to S K Ivanov for useful discussions.