Nonlinear Schr\"odinger equations and the universal description of dispersive shock wave structure

The nonlinear Schr\"odinger (NLS) equation and the Whitham modulation equations both describe slowly varying, locally periodic nonlinear wavetrains, albeit in differing amplitude-frequency domains. In this paper, we take advantage of the overlapping asymptotic regime that applies to both the NLS and Whitham modulation descriptions in order to develop a universal analytical description of dispersive shock waves (DSWs) generated in Riemann problems for a broad class of integrable and non-integrable nonlinear dispersive equations. The proposed method extends DSW fitting theory that prescribes the motion of a DSW's edges into the DSW's interior, i.e., this work reveals the DSW structure. Our approach also provides a natural framework in which to analyze DSW stability. We consider several representative, physically relevant examples that illustrate the efficacy of the developed general theory. Comparisons with direct numerical simulations show that inclusion of higher order terms in the NLS equation enables a remarkably accurate description of the DSW structure in a broad region that extends from the harmonic, small amplitude edge.


INTRODUCTION
There has been a surge of interest recently in the subject of dispersive hydrodynamics and, in particular, dispersive shock waves (DSWs) (see Refs. 1 and 2 and references therein). This has largely occurred thanks to the growing recognition of the fundamental nature and ubiquity of DSWs in physical applications: from shoaling tsunami waves [3][4][5] and internal undular bores in the ocean [6][7][8] and atmosphere 9,10 to nonlinear diffraction patterns and optical shocks in laser beam propagation, [11][12][13][14][15] quantum shocks in superfluids, [16][17][18][19] and nonlinear spin wave propagation in magnetic thin films. 20 On the other hand, the study of DSWs has revealed a number of challenging mathematical problems in the context of both integrable and nonintegrable nonlinear wave equations.
A DSW is an expanding, modulated nonlinear wavetrain that connects two disparate hydrodynamic states ( Figure 1). It can be viewed as a dispersive counterpart to the dissipative, classical shock. Hydrodynamic wave breaking singularities in dispersive media are generically resolved by DSWs. A DSW has a distinct multiscale structure consisting of an oscillatory transition between two nonoscillatoryfor example, slowly varying or constant-states: one edge is associated with a solitary wave or soliton (for convenience, we use the term soliton regardless of the integrability of the governing equation) that is connected, via a slowly modulated periodic wavetrain, to a harmonic, small-amplitude wave at the opposite edge. The relative position (left/trailing or right/leading) of the soliton and harmonic edges determines the DSW orientation , found in terms of the curvature of the linear dispersion relation as = −sgn[ 0 ( , 0 )]. 2,21 Here, = 0 ( , 0 ) is the frequency of a small amplitude wave with wavenumber that propagates on the mean flow background 0 . The DSW shown in Figure 1 has = 1 because the solitary wave is on the rightmost, leading edge. The shock structure of a DSW-an unsteady oscillatory wavetrain-is more complex than the stationary shock structure of a viscous shock wave. In particular, a DSW cannot be described by a traveling wave (ODE) solution of the nonlinear wave equation. 21 The rapidly oscillating structure of DSWs motivates the use of asymptotic, WKB-type, methods for its analytic description. One such method, known as Whitham modulation theory 22,23 (see also, Ref. 24), is based on the averaging of dispersive hydrodynamic conservation laws over nonlinear periodic wavetrains leading to a system of first-order quasilinear partial differential equations. Whitham theory has proved particularly effective for the description of DSWs in both integrable and nonintegrable systems. If the dispersive hydrodynamics are described by an integrable equation such as the Korteweg-de Vries (KdV) or nonlinear Schrödinger (NLS) equation, the associated Whitham system can be represented in a diagonal, Riemann invariant form. 22,24,25 This fact enabled Gurevich and Pitaevskii (GP) 26 to construct an explicit modulation solution for a DSW generated by x F I G U R E 1 DSW expanding oscillatory structure in convex dispersive hydrodynamics with negative dispersion a Riemann problem for the KdV equation. The GP construction is based on a self-similar, rarefaction wave solution of the KdV-Whitham equations. This modulation solution describes the interior shock structure of a DSW and reveals a monotonic change in the nonlinear wave's wavenumber, mean, and amplitude as the DSW is spatially traversed.
For nonintegrable dispersive equations, diagonalization of the associated modulation systems in terms of Riemann invariants is generally not possible, often presenting an insurmountable obstacle to the explicit determination of the Whitham system's simple wave solution, although its existence requires only strict hyperbolicity and genuine nonlinearity. 68 Consequently, the analytical description of a DSW's interior structure has so far been limited to integrable systems or a detailed analysis of the Whitham modulation system in certain limiting regimes on a case-by-case basis. 27,28 One can, however, explicitly determine key observables associated with each DSW edge, even for nonintegrable equations. These observables include the DSW edge speeds and their associated wave parameters-the harmonic edge wavenumber and the soliton edge amplitude. The determination of these observables represents the fitting of a DSW to the long-time dynamics of piecewise constant, initial Riemann data. The DSW fitting method proposed in Ref. 29 (see also, Ref. 2) is based on a fundamental, generic property: the Whitham modulation equations admit exact reductions to a set of common, much simpler, analytically tractable equations in the limits of vanishing amplitude and vanishing wavenumber, which correspond to the harmonic and soliton DSW extremes ( Figure 1). Therefore, the DSW fitting method can be viewed as a universal dispersive hydrodynamic analog to the Rankine-Hugoniot conditions for dissipative, viscous shocks. The key advantage of the method is that it involves neither the determination nor the analysis of the full Whitham system because the required zero-amplitude and zero-wavenumber reductions are available directly and are ultimately determined by the nonlinear, hyperbolic flux and the linear dispersion relation of the dispersive hydrodynamics. The method greatly expands the scope of DSW analysis as it is not reliant on integrability of the governing nonlinear dispersive equation via inverse scattering theory. It has been successfully applied to many nonintegrable dispersive hydrodynamic systems. See, for example, Refs. 28,30-37. Restrictions to the method's applicability are related to possible violations of genuine nonlinearity (convexity) and strict hyperbolicity of the modulation system. 33,34 Once the parameters of the leading and trailing edges have been determined by DSW fitting, wave modulation in the vicinity of these edges can be, in principle, determined by expanding the full Whitham system for small amplitudes near the harmonic edge and small wavenumbers near the soliton edge. 27,28 Such an asymptotic consideration, however, has a number of serious drawbacks due to the need to derive and analyze the full modulation system. Apart from being a potentially daunting technical task, this presents a major disadvantage to the whole procedure as it is system specific. It would therefore be highly desirable to have a more direct, widely applicable method for the determination of the DSW structure including modulation near the harmonic and soliton edges, which would complement and extend the existing general DSW fitting procedure.
In this paper, we develop a general method for the determination of the universal nonlinear DSW modulation-the DSW structure-near the harmonic edge. This asymptotic modulation provides crucial information about the variation of the amplitude in the DSW (ie, the envelope of the oscillatory wavetrain) as well as other physical DSW parameters such as the wavenumber and mean flow. The modulation is universal because it is derived from the NLS equation, a universal model of weakly nonlinear, modulated dispersive wavetrains 38 (see also, Ref. 39). The present work takes advantage of the asymptotic overlap region in the vicinity of the DSW harmonic edge between the semiclassical, long-wave limit of the NLS equation and the small amplitude limit of the full Whitham modulation equations. The commonalities and differences between Whitham modulation theory and the NLS equation have been widely discussed in the literature (see, eg, Refs. 2 and 40) but to the best of our knowledge, the overlap regime for the applicability of the two descriptions has never been used in practice, except very recently in Ref. 41. While the Whitham equations describe modulations of fully nonlinear wavetrains, the NLS description is advantageous in the weakly nonlinear regime because it incorporates higher order dispersive effects of the wave envelope that are not accounted for in leading order Whitham theory.
We use the parameters obtained by DSW fitting applied to the harmonic edge as input information for a standard, small amplitude, multiple scales expansion that leads to the NLS equation for the wave's envelope and phase. The specific information relevant to dispersive hydrodynamics consists of the NLS' nonlinear and dispersion coefficients. The universal asymptotic modulation in a DSW is then found as a special vacuum rarefaction simple wave solution of the NLS equation in the long-wavelength, dispersionless limit, which is similar to the solution to the shallow water equations for the classical dam break problem with a dry downstream bed.
We consider several representative integrable and nonintegrable examples to illustrate the efficacy of the developed general theory. Comparisons with direct numerical simulations show that the accuracy of the asymptotic description improves dramatically when higher order terms of the NLS equation are taken into account in the so-called HNLS equation. The HNLS equation was first derived in the nonlinear optics context, 42,43 but it is a universal equation that also arises in other applications including fluid dynamics 44,45 and plasma physics. 46 We observe that in all considered examples, the vacuum rarefaction simple wave solution of the semiclassical, dispersionless HNLS equation provides a remarkably accurate description of the DSW modulation, and therefore the DSW structure, in a broad vicinity of the harmonic edge. Finally, we show that convexity of the linear dispersion relation for the original dispersive hydrodynamics plays a key role in the determination of DSW stability via the focusing/defocusing character of the associated NLS equation.
The paper is organized as follows. We begin in Section 2 with a brief outline of the necessary elements of DSW modulation theory and, in particular, the DSW fitting method for the determination of the DSW harmonic edge in scalar dispersive hydrodynamic systems. Section 3 develops an asymptotic, multiple scales expansion in the vicinity of the DSW harmonic edge that leads to the NLS equation. This is used to derive the universal, first-order approximation of the DSW modulation as a vacuum rarefaction simple wave solution of the long-wave, dispersionless limit of the NLS equation. We then extend the first-order analysis by including higher order terms in the multiple scales expansions in Section 4. This results in the HNLS equation for which we find the appropriate simple wave solution in the longwavelength limit. Section 5 is devoted to applications to several representative examples. The examples include the KdV equation, the conduit equation that models the interfacial dynamics of a rising, buoyant, viscous fluid within another miscible, high viscosity contrast fluid, 33,47 and the Serre equations for fully nonlinear shallow water waves. 28,48,49 The latter two equations are nonintegrable. We complete the paper with a summary and discussion of further directions in Section 6. Appendices A, B and C detail the multiple scales derivations of the NLS and HNLS equations for the conduit equation and the Serre system. Appendix D is devoted to a brief description of numerical methods used for simulations.

DISPERSIVE SHOCK WAVES: MODULATION THEORY
In this section, we outline the elements of DSW modulation theory that are necessary for developments in subsequent sections. A detailed exposition can be found in Refs. 2 and 21.

Modulation equations and the matching regularization of the Riemann problem
We consider scalar dispersive hydrodynamics generically described by the equation with ′′ ( ) ≠ 0. The dispersive operator (generally integro-differential) is assumed to have the property that Equation (1) admits the real-valued linear dispersion relation = 0 ( , 0 ) with long-wave expansion for small-amplitude waves proportional to ( − ) and propagating on a constant (or slowly varying) background = 0 . We will initially assume that the dispersion relation is purely convex or concave, so that 0 ≠ 0. Suppose that Equation (1) has a three-parameter periodic traveling wave solution and at least two local conservation laws. These basic requirements are quite generic and are satisfied by many dispersive hydrodynamic equations arising in applications. When the dispersive hydrodynamics admit the above properties, we say that they are of KdV type.
We will consider the evolution of Riemann step initial data for Equation (1). Our consideration will be based upon the fundamental assumption that the initial step (3) is regularized in the long-time limit by the emergence of three distinct regions in the -upper half space-time plane so that the solution is given by two constant states = − and = + that are separated by an expanding DSW region ( Figure 2). Within this region, the solution has an oscillatory structure described by a modulated, locally periodic wavetrain that exhibits a solitary wave at one edge and a vanishing amplitude linear wave at the opposite edge (recall Figure 1). This asymptotic structure of the Riemann problem solution has been rigorously recovered for a number of integrable equations (see, eg, Refs. 50 and 51). For nonintegrable equations such as the Serre system 28 or conduit equation, 33 the existence of a modulated, periodic, single phase wave structure of a DSW is a plausible assumption, which can be inferred from numerical simulations. We assume that Equation (1) admits a three parameter family of periodic, traveling wave solutions ( , ) = ( ), where = − , being the wavenumber and the wave frequency, so that ( + 2 ) = ( ). It is convenient to use the period mean = (2 ) −1 ∮ d , the amplitude = max − min , and the wavenumber as a basic parameter set, that is, ( ) ≡ ( ; , , ); all other physical parameters, such as the frequency or the mean square 2 are functions of the basic triple ( , , ). We also assume that the solution ( ; , , ) has two asymptotic limits: (i) when → 0, it transforms into a linear wave on the background = with the dispersion relation = 0 ( , ); (ii) when → 0, it transforms into an exponentially decaying solitary wave. Examples of dispersive equations whose periodic solutions satisfy the above properties include KdV, modified KdV, the conduit equation, and others. We now consider a solution of Equation (1) represented by the 2 -periodic traveling wave with slow ( , )-dependence of ( , , ). This slowly varying traveling wave is characterized by the generalized phase ( , ) so that the local wavenumber and frequency are given by = and = − , respectively. The variations of ( , , ) satisfy the Whitham modulation equations, 23 which can be obtained by applying multiple scales expansions or, equivalently, by averaging two independent conservation laws of (1) over the periodic family and completing the system with the consistency equation = that yields wave number conservation + = 0. The same set of modulation equations can be derived via an averaged variational principle. 52 Assuming nonvanishing Jacobians, the Whitham system can be represented as a system of quasilinear first-order equations where A( , , ) ∈ ℝ 3×3 is a matrix. We initially assume hyperbolicity so that the eigenvalues 1 ≤ 2 ≤ 3 of A are real and the eigenvectors { | A = , = 1, 2, 3} form a basis in ℝ 3 . In the context of a DSW that is described by a modulated periodic wave solution, the Whitham equations (4) are subject to free boundary (matching) conditions 26,29 thus ensuring continuity of the mean flow at the unknown DSW edges = ± ( ). Outside the DSW region − ( ) ≤ ≤ + ( ), the solution is given by = − for < − ( ) and = + for > + ( ). Here, for specificity, we have assumed a positive DSW orientation ( Figure 1) so that the harmonic edge is trailing, = − ( ), and the soliton edge is leading, = + ( ). We also assume concave flux, ′′ ( ) > 0, which implies the admissibility or causality condition − > + for a compressive DSW. 2,21 The cases of positive dispersion that yield either a negative DSW orientation, = −1, or convex flux ′′ ( ) < 0 (which implies − < + for compressive DSW formation) can be considered in a similar fashion. The hydrodynamic scaling invariance → , → of both the modulation equations (4) and initial conditions (3), together with hyperbolicity, necessitates a self-similar, simple wave modulation solution. To satisfy the matching conditions (5), the modulation solution must be a 2-wave rarefaction curve 2,29 given by where 1,2 are integrals of the Whitham system (4) on the solution curve. The integrals are parametrized by the Riemann data ± , that is, 1 ( − , 0, − ) = 0 and 2 ( + , + , 0) = 0 determine the trailing edge harmonic wavenumber − and the leading edge soliton amplitude + . For the KdV equation solution, Equation (7) was found explicitly by GP 26 in terms of Riemann invariants that are available for the KdV-Whitham system owing to its complete integrability; 25 in fact, Whitham was able to determine the Riemann invariants explicitly via an involved, direct calculation 22 (see also, Refs. 24 and 53). The DSW edge speeds ± are constant and follow from the modulation solution (7) in the → 0 (harmonic, trailing edge) and → 0 (soliton, leading edge) limits by evaluating the linear group velocity and soliton phase velocity, respectively. The ( , )-contour plot of the GP solution to the Riemann problem for the KdV equation + + = 0 illustrates the described modulation theory setting in Figure 2.
The above outlined construction of the DSW modulation solution for scalar equations (1) can be extended to systems describing bidirectional dispersive hydrodynamics, for example, the Serre shallow water equations, the generalized NLS equation, and other systems. See, for example, Refs. 28, 29, 34, and 54.

The determination of the harmonic edge: DSW fitting
Modulation systems (4) obtained by averaging dispersive hydrodynamic equations (1) exhibit an important general property: they admit exact reductions to lower order quasilinear systems in the harmonic ( → 0) and soliton ( → 0) limits (recall that these two limits correspond to special wave regimes realized at the DSW edges 2,27,29 ). Importantly, these exact reductions are often available explicitly, without the necessity to derive the full modulation system. Another fundamental fact is that in the Riemann problem, the DSW edges are characteristics when the modulation system (4) is hyperbolic. As a result, the speeds ± of the harmonic and soliton edges can be obtained from the analysis of the reduced modulation systems. The corresponding technique proposed in Ref. 29 is sometimes referred to as the DSW fitting method.
Determining the harmonic edge via the DSW fitting method is particularly simple. Indeed, the harmonic reduction ( = 0) of the modulation system (4) can be shown to be universally represented in the physically transparent form Then, the DSW harmonic edge speed coincides with the linear group velocity for the edge parameters = − , = − : ( ) being the harmonic edge wavenumber locus function, which is determined as follows. Let the value = + at the DSW soliton edge be fixed. Then, the function ( ) is found from the ODE The ODE (10) follows by integrating the differential associated with the group velocity characteristic of Equation (8)  The determination of the soliton edge is analogous, although it involves some nontrivial change of variables which we do not describe here (see Refs. 2 and 29). The extension of scalar DSW fitting to bidirectional, Eulerian dispersive hydrodynamic systems has been developed in Refs. 29 and 34. The DSW fitting procedure is subject to certain admissibility conditions derived from the requirement of monotonicity for the relevant real characteristic velocity along the integral curve (7), that is, the modulation system must be genuinely nonlinear and strictly hyperbolic along the entire integral curve. 33,34 Therefore, the DSW fitting construction is not reliant upon the integrability of the dispersive hydrodynamic evolution equation (1), it only requires strict hyperbolicity and genuine nonlinearity of the modulation system (4).

SMALL AMPLITUDE DSW REGIME AND THE NLS EQUATION
We will be interested in the region of a DSW adjacent to the harmonic edge = − , where the oscillation amplitude is relatively small. One can, in principle, expand the Whitham equations (4) in powers of for ≪ 1 and solve them by seeking a solution in powers of the amplitude, subject to the matching condition (5) at = − . This programme has been to some extent realized in Refs. 27 and 28 for several nonintegrable dispersive hydrodynamic systems, including the equations for ionacoustic and magneto-acoustic waves in collisionless plasma and the Serre equations for fully nonlinear dispersive shallow water waves. In all cases, the full modulation system (4) (or its bi-directional generalization) was derived, and then reduced to an abstract, model Whitham system for and (see Ref. 23, ch. 16.15) involving an effective nonlinear frequency correction 2 ( , − ) 2 to the linear dispersion relation 0 ( , − ). As a result, the first-order DSW modulation near the harmonic edge was determined, including the linear growth of the amplitude with distance. This approach, however, has a major drawback, due to its reliance on the full modulation system in each case, while only the small amplitude expansion is actually used.
Here, instead of deriving the full modulation system and making subsequent small amplitude expansions, we perform an appropriate small amplitude expansion directly on the original system and then derive modulation equations for weakly nonlinear, periodic (Stokes) waves. 23 Slow modulations of almost monochromatic Stokes waves for a broad class of nonlinear dispersive systems are known to be described by the NLS or its higher order versions. 38,55,56 Consequently, both Whitham modulation theory and an NLS description can be used in the inner vicinity of the DSW harmonic edge. Moreover, because a DSW is described by a rarefaction wave solution of the Whitham equations, the counterpart NLS description will reduce to finding a simple wave solution to the dispersionless limit of the corresponding NLS equation or one of its higher order versions.
It is instructive to start with an outline of the standard derivation of the NLS equation. See, for example, Ref. 55 for examples and further details. Let be a small parameter characterizing the wave amplitude. We seek the solution of the dispersive hydrodynamic equation (1) in the form of an asymptotic expansion about the constant 0 for a nearly monochromatic wave with the dominant carrier wavenumber where Substituting expansion (11) into Equation (1) and collecting terms in powers of , we obtain the linear dispersion relation = 0 ( , 0 ) at the first order. To eliminate secular terms at ( 2 ), we require that the complex wave envelope moves with the group velocity The NLS equation arises as the condition for removal of secular terms at ( 3 ) and has the form where ( , 0 ) = 1 2 0 . We also obtain the variation of the mean − 0 ∼ 2 1 ( , 0 )| | 2 as a byproduct of the ( 3 ) calculation. Here, the factors 1 and have no general expressions.
Although the above outlined derivation is standard, it can be quite laborious, especially for systems. The difficult part of the derivation is the determination of the nonlinear coefficient ( , 0 ), but this computation can be readily performed with the aid of a symbolic algebra package such as Mathematica. See Appendix B for an outline of the calculations for the Serre system. In fact, is precisely the sought for nonlinear frequency correction 2 ( , ) mentioned earlier that is obtained in a weakly nonlinear analysis of the Whitham equations. Equations for 1 and 2 can be combined into a single equation for the unscaled envelopẽ( , ) defined bỹ Hence, one has the following substitution rules: The envelope of the wave packet = 0 +̃( , ) ( − 0 ) + c.c. is then governed by the equation: The sign of the product determines the NLS type: if > 0, the equation is attractive or focusing and describes the envelope of a modulationally unstable wave while for < 0, it is repulsive or defocusing and describes the envelope of a modulationally stable wave.
To apply the NLS equation to the description of the DSW harmonic edge vicinity, we assume in (11): where the dependence of − on the Riemann data ± is obtained by DSW fitting. We introduce the Madelung transformatioñ= √ , = to cast the NLS equation (17) in dispersive-hydrodynamic form where ∼ − − , and √ = |̃| ∼ ∕4 in the DSW context (we recall that = max − min ). The DSW modulation solution is a rarefaction curve of the Whitham equation, thus the relevant NLS solution must also be a rarefaction wave described by the long-wave, dispersionless limit. Neglecting the dispersive term in (18), we obtain The characteristic velocities are ± √ −2 , thus the system is hyperbolic if < 0 and elliptic if > 0, consistent with the defocusing and focusing character of the NLS equation (13), respectively. We assume for now that < 0, so that the system (19) is equivalent to the shallow water equations.
We now need to specify boundary conditions for the dispersionless NLS equation (19) at the DSW harmonic edge. This is done using the GP matching conditions (5) and the DSW fitting data. In the original modulation variables, we have from (5) Note that, unlike the free boundary in (5), the boundary in (20) is known from DSW fitting. Translating (20) into the variables of (19), we arrive at a boundary value problem for the vacuum rarefaction wave We are now looking for a self-similar rarefaction wave solution of the shallow water equations (19) subject to the boundary conditions (21). There are two such nontrivial solutions-the fast and slow waves. The correct one is chosen by the admissibility condition that the wavenumber decreases monotonically as the DSW is traversed from the harmonic to the soliton edge, that is, ∕ < 0 or, equivalently, ∕ < 0. Then, the required solution of (19) is given by Using the dispersionless NLS solution (22) and the expansions (11), we recover the leading order behaviors of the amplitude and the wavenumber near the DSW harmonic edge in terms of the dispersion and nonlinearity coefficients and of the associated NLS equation (17) We also recover the variation of the mean: Equation (23) is the universal description of the DSW envelope (with "martini-glass" shape 2,57 ) near the harmonic edge for systems with convex dispersion ( ≠ 0). Solution (23) is valid when < 0, which is the hyperbolicity condition for the dispersionless NLS system (19) and can be interpreted as a necessary condition for DSW modulational stability. For nonconvex dispersion relations, the sign of can change and the system may exhibit an unstable behavior described by the focusing NLS equation where > 0. An example of such behavior has been recently reported in Refs. 33 and 41, where it was shown that nonconvexity of the linear dispersion relation for the conduit equation implies an elliptic regime for the associated Whitham equations in a certain region of parameter phase space. This gives rise to modulationally unstable dynamics near the DSW harmonic edge that can be described by the focusing NLS equation in the small amplitude regime.
For systems with nonconvex dispersion, the study of DSW behavior near the zero dispersion point = 0 necessitates inclusion of higher order terms in the associated NLS equation. It turns out that the inclusion of such terms is beneficial even outside the nonconvex, zero-dispersion regime, as we now demonstrate.

HIGHER ORDER NLS APPROXIMATION
We can include higher order nonlinear/dispersive effects in the asymptotic expansion (11) by introducing a third, slower time scale 3 = 3 and assuming 1 = ( , 1 , 2 , 3 ) ( − ) + c.c. The cancellation of secular terms at ( 4 ) then gives where ( , 0 ) = − 0 ∕6 and * is the complex conjugate of . Once again, the laborious part of the computation consists of finding the coefficients ( , 0 ) and ( , 0 ). Written in the moving reference frame ( , ) (cf. Appendix A for the derivation of the unscaled equation), where = − 0 ( , 0 ) (recall Equation (17)), the resulting equation is the higher order NLS, or HNLS, equatioñ initially derived as an improvement to the standard NLS equation for signal propagation in optical fibers; 43 it also arises in geophysical fluid dynamics 44  Similar to the previous section, we cast the HNLS equation (25) in dispersive hydrodynamic form using the Madelung transform̃= √ ∫ . Upon neglecting dispersive terms, we obtain the following quasilinear system for long waves As expected, the system (26) is equivalent to the shallow water equations (19) when = = = 0. In the context of DSWs, the dispersionless limit (26) of the HNLS equation should be considered with the same vacuum rarefaction conditions (21) augmented by the DSW admissibility inequality ∕ < 0. Before we proceed with the integration of system (26), let us briefly discuss its structure. The characteristic velocities are (27) and the associated right eigenvectors are implying that the system (26) is hyperbolic in the region ( , ) where > 0, the discriminant is positive > 0, and ≠ ( − ) so that ± are independent. In the small amplitude regime ( , | | ≪ 1), we recover the standard hyperbolicity condition < 0 because , ≠ 0. Consequently, the DSW modulation near the harmonic edge is determined by the similarity solution of (26), where = −sgn is the DSW orientation and the dependence ( ) is determined by the rarefaction curve As a by-product of the multiscale expansion to order ( 4 ), we obtain a higher order correction to the mean value (recall (24)) described by the new expression where 1 has been obtained at the previous order ( ( 3 )), and 2 is determined as a by-product of the ( 4 ) computation.
An exact analytical solution of the simple wave ODE (30) is generally not available, thus one can resort to numerical solution. However, because Equations (29) and (30) can formally be applied only to the small-amplitude vicinity of the DSW harmonic edge-that is, when and are small-we seek a series solution near = 0. The asymptotic solution of (29) and (30) in the form of a series in powers of ∕ yields general expressions for the approximate DSW amplitude and wavenumber modulations As expected, the first-order terms in (32) agree with the NLS result (23). We note that it is implicit in the expansions (32) that | ( − , − )| = (1), which is the case for dispersive hydrodynamic equations with a convex dispersion relation. For systems with nonconvex dispersion such as the Benjamin-Bona-Mahony equation 58 or the conduit equation, 33,47 the behavior near the zero-dispersion point, ( − , − ) = 0 captured by the HNLS equation (25) requires a separate consideration, which is beyond the scope of the present paper.
For convex dispersive hydrodynamics, the second-order approximation (32) formally delivers the same accuracy of the DSW harmonic edge vicinity as the HNLS equation (25) itself. However, comparisons with results of direct numerical simulations of the Riemann problem for the example dispersive hydrodynamic equations in the next section show that the simple wave solution (29) and (30) of the full dispersionless HNLS (26) is more accurate than the expansion (32). Remarkably, the accuracy holds over a significant portion of the DSW, where the amplitude is not small and the HNLS description, let alone the expansions (32), is not expected to be applicable.

NLS DESCRIPTION OF DISPERSIVE SHOCK WAVES: EXAMPLES
We now demonstrate the effectiveness of the developed general approach by applying it to several specific dispersive hydrodynamic equations and comparing the results with direct numerical simulations of the corresponding Riemann problems.

Korteweg-de Vries equation
As a first example, we consider the KdV equation with Riemann initial data (3). The aim here is to compare the results of the developed asymptotic approach with the known GP modulation solution. 26 The KdV linear dispersion relation is The multiple scales asymptotic expansion of (33) leading to the NLS equation for KdV weakly nonlinear wavepackets is standard and can be found in the literature, see, for example, Refs. 55 and 59. The coefficients in (13) and (24) are We also derive the coefficients of the HNLS equation (25) and the associated higher order correction of the mean value (31) 59 The trailing edge wavenumber is readily obtained from DSW fitting (see Ref. 29 and section 2.2). Solving the ODE (10), we obtain − = ( − ) in the form which yields the harmonic edge velocity Substituting (36) and (37) into (31) and (32), we obtain the second-order expansions which agree with the corresponding expansions of the exact GP solution. 26 The comparison between the exact simple wave solution (29) and (30) of the dispersionless HNLS equation (26) with coefficients (35) and (39), the asymptotic expansion (39)- (41), and the direct numerical solution of the KdV Riemann problem is displayed in Figure 3. The numerical method used in the simulations is detailed in Appendix D. Although the NLS description is formally limited to the small-amplitude regime in the vicinity of the DSW harmonic edge, the KdV DSW amplitude is almost linear for ∕ ∈ [ − , + ], thus the term linear in ∕ in (39) fits almost the entire DSW with good accuracy. It is not the case for the wavenumber and the mean flow: good agreement between numerics and the first-order (NLS) approximation is observed only in the vicinity of the trailing edge, but the second-order expansions (40) and (41) exhibit better agreement with the direct KdV numerics over a broader ∕ interval. The agreement further improves when the full solution (29) and (30) of the dispersionless HNLS equation is used where all modulation parameters fit almost the entire DSW with good accuracy, with the exception of some vicinity of the leading edge, where → 0 and → + logarithmically in ∕ , 26 behavior that cannot be captured by the expansions (40) and (41).
We should make an important comment regarding the comparison of the numerically computed DSW oscillations near the harmonic edge with the results of modulation theory. We can see from Figure 3, left panel, that there is some deviation of the envelope profile near the trailing edge from linear behavior ∝ ( ∕ − − ) that is predicted by modulation theory. In particular, the amplitude of the oscillations is not exactly zero for ∕ ≤ − . This discrepancy between the DSW modulation solution and the exact oscillation behavior is known, having been studied for the KdV equation in detail in Ref. 60, where it was shown that the envelope amplitude difference between the numerical KdV solution and the modulation theory solution in the region ∕ < − decreases roughly as −1∕3 , likely due to higher order dispersive effects that are not captured by the long-wave approximation utilized here. However, modulation theory typically provides a very satisfactory prediction of the amplitude growth near the DSW harmonic edge even for relatively moderate times.

Conduit equation
We now consider the conduit equation at time and vertical spatial coordinate , separating a light, viscous fluid rising buoyantly through a heavier, more viscous, miscible fluid at small Reynolds numbers. 33,47 Equation (42) has convex hyperbolic flux ( ) = 2 and linear dispersion relation which is nonconvex as 0 can change sign. The coefficients of the NLS equation (13)  ) , . (44) We see that, while the nonlinearity coefficient is always positive, the dispersion coefficient (and therefore the parameter ) can change sign, thus the parameter space ( , 0 ) of conduit Stokes waves is split into two domains-which correspond to the defocusing and focusing NLS regimes-that are separated by the line = √ 3∕ 0 . In the context of DSWs, the line − ( − , + ) = √ 3∕ − in the − -+ phase plane of Riemann data (3) separates the regimes of DSW harmonic edge stability and instability. Here, is the conduit DSW harmonic edge wavenumber obtained from DSW fitting. 33 The condition − < √ 3∕ − or, equivalently, + ∕ − > 5∕32 ≈ 0.156 is the DSW fitting admissibility condition whose violation was associated in Ref. 33 with a gradient catastrophe for the wavenumber and a subsequent DSW implosion-the formation of a two-phase region near the trailing edge. Within the NLS description of DSW modulations developed here, the above admissibility condition is naturally interpreted as a conduit DSW modulational stability condition. The plots of DSWs for stable and unstable regimes are presented in Figures 4(A) and (B), respectively.
We now compare the predictions of the (H)NLS analysis for the conduit DSW modulation with direct numerical simulations within the admissible range of Riemann data with + ∕ − > 5∕32 that produces a stable DSW. Within this region, < 0, thus the DSW orientation = +1, and the harmonic edge is the trailing one, see Figure 4. The coefficients , , in the conduit-HNLS equation (25) and the coefficient The comparisons between the dispersionless HNLS vacuum rarefaction simple wave solution (29) and (30) with coefficients (44) and (A.6), its asymptotic expansions (23) and (32), and the numerical solution of the Riemann problem for the conduit equation are shown in Figure 5. As in the previous cases, the full simple wave solution exhibits very good agreement with the direct numerical solution over a broad DSW region, while the first-and second-order approximations work satisfactorily only in a relatively narrow vicinity of the harmonic edge.

Serre equations
The presentation until now has emphasized scalar dispersive hydrodynamic equations in the form (1). Our methodology, however, can be applied to systems of dispersive hydrodynamic equations. As an example, we now consider the Serre system modeling fully nonlinear shallow water waves 48,49 + ( ) = 0, We defer technical calculations to Appendices B and C. Here, is the total depth of the fluid and is the depth averaged horizontal velocity. The linear dispersion relation of (46) for small amplitude waves propagating on the background ( 0 , 0 ) has the form We assume Riemann initial data for (46) subject to an additional constraint (a 2-wave rarefaction curve) that ensures a simple wave, 2-DSW resolution of (48) corresponding to the fast "+" mode in the dispersion relation (47). 2,28 Due to scaling and Galilean symmetries of (46), we can assume + = 1, + = 0. While the Serre system (46) is not integrable, it satisfies the prerequisites of the DSW fitting method except for the loss of genuine nonlinearity in a certain parameter regime, discussed further below. The corresponding analysis has been carried out in Ref. 28, demonstrating excellent agreement with direct numerical simulations. An implicit expression for the harmonic (trailing) edge wavenumber − of a 2-DSW as a function of the Riemann data (48), obtained in Ref. 28 by integrating a bi-directional generalization of the ODE (10), has the form Here, Δ = − ∕ + and − = + 0 ( − , − , − ). The result (50) is valid as long as − ∕ − < 0 (DSW fitting admissibility 34 ) leading to the condition Δ < Δ ≈ 1.43. 28 When Δ = Δ , the Whitham modulation system loses genuine nonlinearity at harmonic trailing edge.
The multiple scales asymptotic expansions for the Serre equations (46) that lead to the NLS equation (13) are carried out in Appendix B. In these expansions, the envelopes of the small amplitude oscillations of ( , ) and ( , ) are proportional to each other: , The coefficients of the NLS equation (13) where = 0 . Here, 1, , 1, are the coefficients in the mean flow expansions = 0 + 1, and = 0 + 1, . The NLS equation for the component̃is obtained by combining the NLS equation for̃and the proportionality relation (51).
Going to ( 4 ), we obtain the coefficients , , in the HNLS equation (25), as well as the coefficients 2, , 2, in the second-order expansions of the mean flow , in terms of , 0 , 0 . These are presented in Appendix C (see formulae (C.5) and (C.6)). In this higher order description, the envelopes and̃are no longer proportional to each other, and the coefficient in the HNLS equation has to be derived separately for each component. To apply the HNLS equation to the description of the vicinity of the Serre DSW harmonic edge, we set 0 = − , 0 = − , = − , where the dependence of − on the Riemann data (48) is obtained numerically from the implicit equation (50).
The comparison between the simple wave solution (29) and (30) with parameters given by (52), (C.5), and (C.6), and the numerical solution of the Riemann problem for the Serre equations is displayed in Figure 6. Also shown are the curves corresponding to the first-order (23) and second-order (32) approximations of the full solution (29) and (30). The comparisons are made for the depth and velocity amplitudes and in the DSW, the wavenumber , the mean depth , and the mean velocity . We can see that the full simple wave solution (29) and (30) of the dispersionless HNLS equation provides a more accurate description of the DSW modulation than the first-and second-order approximations. Similar to the KdV case, the full simple wave solution of the dispersionless HNLS equation provides a good approximation of the nonlinear wave modulation over a significant portion of the DSW, well beyond the formal applicability of the small amplitude (H)NLS approximation. On the other hand, we can see that, in contrast to the KdV case, the second-order approximation (32) develops quite strong deviation from the actual modulation for moderate values of ( ∕ − − ).
A final comment on the comparison concerns the already discussed generic discrepancy between the simple wave modulation DSW solution and direct numerical solution of the Riemann problem in the vicinity of the harmonic edge (see the discussion at the end of Section 5.1). This discrepancy is more pronounced for the Serre equations than for the KdV equation although the overall agreement with the dispersionless HNLS solution is still quite good.
We note in conclusion that the (H)NLS description of 1-DSWs is completely analogous to the 2-DSW description, with condition (49) replaced by the rarefaction curve for a 1-wave and using the "−" mode in the dispersion relation (47). Combined with known rarefaction wave solutions for the classical shallow water equations and the bi-directional DSW fitting construction, 29,34,54 the theory developed here gives the DSW structure for the full Serre equation Riemann problem.

CONCLUSION AND DISCUSSION
In this work, we have developed an efficient, universal approach for the analytical description of the interior structure of a DSW that extends the previously developed DSW fitting method 29 for the DSW edge speeds. The key element of the extension is the realization that the DSW modulation described by an expansion fan solution of the Whitham modulation equations can be universally approximated, in the vicinity of the weakly nonlinear harmonic edge, by a special vacuum rarefaction solution of the shallow water equations. The connection between the original dispersive hydrodynamics and the approximating shallow water system occurs via a long-wave, dispersionless limit of the NLS equation for weakly nonlinear, narrow-band Stokes waves, whose parameters are determined by DSW fitting when the NLS equation is of defocusing type. The NLS type (defocusing or focusing) determines DSW stability properties. The developed approach is particularly attractive for applications as it allows one to avoid a potentially complex, full Whitham modulation analysis of DSWs in favor of the more straighforward and standard NLS theory.
The efficacy of the developed approach is demonstrated by several representative examples including the KdV equation, the Serre shallow water equations, and the viscous fluid conduit equation, the latter two systems being nonintegrable. In all considered cases, it is shown that the inclusion of higher order terms in the NLS equation dramatically improves agreement between the approximate modulation solution and the numerical solution of the original dispersive Riemann problem. The proposed method has broad implications for DSW analysis in nonintegrable systems, where exact methods based on inverse scattering theory are not available. One interesting perspective is to use the NLS approximation for the analytical description of multiphase modulations that are symptomatic of DSW implosions (see Ref. 33 and Section 5.2 of this paper). In this context, the description will necessarily depend upon dispersive terms in the NLS equation, which do not play a role in the classical expansion fan DSW solutions considered in the present paper. Further, the improved DSW description in Section 4, based on the higher order NLS equation, provides a general mathematical framework for DSW analysis in systems with nonconvex dispersion, which are currently under active investigation. 61,62 We also envisage intriguing connections with the multisymplectic theory of universal dispersive deformations of the Whitham equations near coalescing characteristics, precisely the configuration that occurs at the DSW harmonic edge; see Refs. 63 and 64 and references therein.
Probably, the most appealing extension of the developed harmonic edge DSW structure theory would be to find its counterpart in the vicinity of the DSW soliton edge (utilizing small asymptotics), and to construct a universal, matched, uniformly valid asymptotic solution for the entire DSW modulation. The cancellation of the secular terms at ( 2 ) and ( 3 ) in the expansion in gives Equations (12) and (13), respectively, and the cancellation of the secular term at ( 4 ) gives As a by-product of the ( 4 ) expansion, we also obtain a higher order correction of the mean which reads: We define the unscaled envelopẽ( , ) by: implying the new substitution rule (cf. Equation (15)): Combining Equations (12), (13), and (A.2) with the substitution rule (A.5), we obtain the HNLS equation (25) for the unscaled envelopẽ( , ).
Applying the above algorithm to the conduit equation (42), we find the coefficients of the HNLS equation (25) to be:

APPENDIX B: DERIVATION OF THE NLS EQUATION FOR THE SERRE SYSTEM
We detail in this Appendix the multiple scales expansions for the Serre equations (46) leading to the NLS equation (13). Although the derivation of the NLS equation for systems is standard (see, eg, Ref. 65), the computation can be rather cumbersome because of the vectorial nature of the system; this difficulty can be overcome using symbolic computations, as we have done in this case (we used Mathematica). Similar to the multiple scales asymptotic expansions for scalar equations, we look for the solution in the form where Ψ ∈ ℂ 2 is a complex two-component vector. Substituting (B.1) into the Serre system (46) and collecting the ( ) terms, we get where = 0 is a convenient parameter and ( , 1 , 2 ) is a scalar. The kernel of is spanned by Assuming that Ψ is given by (B.5) and = + 0 ( , 0 , 0 ), the second order of the asymptotic expansion reads where we drop the dependences of the dispersion relation + 0 by convenience. The vectors 1 and 2 are given by: Because det ( , − + 0 ) = 0, a compatibility condition is necessary to solve Equation (B.7) (cf. for instance, Ref. 65), thus we impose the orthogonality requirement + 0 ⋅ 1 = 0. This condition is satisfied if the wave packet propagates with the group velocity + 0 ( , 0 , 0 ): Providing that (B.8) is respected, one solution of (B.7) is: where ( , 1 , 2 ) and ( , 1 , 2 ) are two unknown fields that remain to be determined. and are necessary for the consistency of the asymptotic expansion, as we will see at the next order (cf. Equations (B.11)).
The solution (B.9) is not unique and the general solution of (B.7) reads as: Ξ 2 + [ ( , 1 , 2 ) ( − + 0 ) + c.c.], where ( , 1 , 2 ) belongs to the kernel of ( , − + 0 ). This additional term does not modify the NLS equation (13) in the end, but it plays an important role in higher order descriptions (cf. Appendix C). In practice, we choose ( , 1 , 2 ) such that one of the two components of Ξ 2 proportional to ( − + 0 ) is equal to 0 (which is already the case here). Finally, if we substitute Ξ 2 by the solution (B.9), the ( 3 ) of the expansion reads We do not present the coefficients 2 and 3 for the second and third harmonic terms at this stage because they are not needed for the derivation of the NLS equation. However, these terms are needed to solve (B.10), and ultimately derive the HNLS equation at the next order. Because (0, 0) = 0, the constant term 0 should be equal to 0. This condition is respected for ( , 1 , 2 ) and ( , 1 , 2 ) given by:  We note that the coefficient in front of̃2̃ * in (C.4) is not proportional to as one might expect if one considered the inadequate definitioñ= . Nonetheless, to the first order, both̃= and definition (C.2) yield the same NLS equation for̃which can be simply obtained by substituting in (16)̃bỹ∕ .
We now present, without derivation, the coefficients for the HNLS equation (25)  and apply the same algorithm.