Shock waves in dispersive Eulerian fluids

# Shock waves in dispersive Eulerian fluids

## Abstract

The long time behavior of an initial step resulting in a dispersive shock wave (DSW) for the one-dimensional isentropic Euler equations regularized by generic, third order dispersion is considered by use of Whitham averaging. Under modest assumptions, the jump conditions (DSW locus and speeds) for admissible, weak DSWs are characterized and found to depend only upon the sign of dispersion (convex or concave) and a general pressure law. Two mechanisms leading to the breakdown of this simple wave DSW theory for sufficiently large jumps are identified: a change in the sign of dispersion, leading to gradient catastrophe in the modulation equations, and the loss of genuine nonlinearity in the modulation equations. Large amplitude DSWs are constructed for several particular dispersive fluids with differing pressure laws modeled by the generalized nonlinear Schrödinger equation. These include superfluids (Bose-Einstein condensates and ultracold Fermions) and “optical fluids”. Estimates of breaking times for smooth initial data and the long time behavior of the shock tube problem are presented. Numerical simulations compare favorably with the asymptotic results in the weak to moderate amplitude regimes. Deviations in the large amplitude regime are identified with breakdown of the simple wave DSW theory.

###### pacs:
03.75.Kk, 03.75.Lm, 05.45.Yv, 42.65.Sf, 47.40.Nm

## 1 Introduction

Nonlinear wave propagation in dispersive media with negligible dissipation can lead to the formation of dispersive shock waves (DSWs). In contrast to classical, viscous shock waves which are localized, rapid jumps in the fluid’s thermodynamic variables, DSWs exhibit an expanding oscillatory region connecting two disparate fluid states. A schematic depicting typical left-going DSWs for positive and negative dispersion fluids is shown in figure 1. These structures are of particular, current interest due to their recent observation in superfluidic Bose-Einstein condensates (BECs) of cold atomic gases [1, 2, 3, 4, 5] and nonlinear photonics [6, 7, 8, 9, 10, 11, 12, 13]. Dispersive shock waves also occur in a number of other dispersive hydrodynamic type systems including water waves [14] (known as undular hydraulic jumps or bores), two-temperature collisionless plasma [15] (called collisionless shock waves), and fluid interfaces in the atmosphere [16, 17, 18] and ocean [19].

The Whitham averaging technique [20, 21] is a principle analytical tool for the dispersive regularization of singularity formation in hyperbolic systems; see e.g., the review [22]. The method is used to describe slow modulations of a nonlinear, periodic traveling wave. Given an order nonlinear evolution equation, implementation of the method requires the existence of a -parameter family of periodic traveling wave solutions , with period and phase . Additionally, the evolution equation must admit conserved densities and fluxes , corresponding to the conservation laws

 Extra open brace or missing close brace (1.1)

Assuming slow spatio-temporal evolution of the wave’s parameters , the conservation laws are then averaged over a period resulting in the modulation equations

 (1L∫L0Pi[ϕ(θ;p)]\rmdθ)t+(1L∫L0Qi[ϕ(θ;p)]\rmdθ)x=0,i=1,…,n−1. (1.2)

The Whitham modulation equations are completed by the addition of the conservation of waves to (1.2)

 kt+ωx=0,k=2π/L=θx,c=−θt, (1.3)

a consistency condition () for the application of modulation theory. The Whitham equations are a set of first order, quasi-linear partial differential equations (PDEs) describing the slow evolution of the traveling wave’s parameters .

As laid out originally by Gurevich and Pitaevskii [23], a DSW can be described by the evolution of a free boundary value problem. The boundary separates the oscillatory, one-phase region, described by the Whitham equations, from non-oscillatory, zero-phase regions, described by the dispersionless evolution equation. The regions are matched at phase boundaries by requiring that the average of the one-phase solution equals the zero-phase solution. Thus, the free boundary is determined along with the solution. There are two ways for a one-phase wave to limit to a zero-phase solution. In the vicinity of the free boundary, either the oscillation amplitude goes to zero (harmonic limit) or the oscillation period goes to infinity, corresponding to a localization of the traveling wave (soliton limit). The determination of which limiting case to choose at a particular phase boundary requires appropriate admissibility criteria, analogous to entropy conditions for classical shock waves.

Riemann problems consisting of step initial data are an analytically tractable and physically important class to study. For a system of two genuinely nonlinear, strictly hyperbolic conservation laws, the general solution of the Riemann problem consists of three constant states connected by two self-similar waves, either a rarefaction or a shock [24, 25]. This behavior generalizes to dispersive hydrodynamics so, borrowing terminology from classical shock theory, it is natural to label a left(right)-going wave as a 1(2)-DSW or 1(2)-rarefaction. See figure 1 for examples of -DSWs where the corresponds to positive or negative dispersion. For a DSW resulting from the long time evolution of step initial conditions, the oscillatory boundaries are straight lines. These leading and trailing edge speeds can be determined in terms of the left and right constant states, analogous to the Rankine-Hugoniot jump conditions of classical gas dynamics. Whitham modulation theory for DSWs was initially developed for integrable wave equations. Integrability in the context of the modulation equations [26] implies the existence of a diagonalizing transformation to Riemann invariants where the Riemann problem for the hyperbolic modulation equations could be solved explicitly for a self-similar, simple wave [23, 27]. The two DSW speeds at the phase boundaries coinciding with the soliton and harmonic limits are the characteristic speeds of the edges of the simple wave. Thus, the dispersive regularization of breaking in a hydrodynamic system is implemented by the introduction of additional conservation laws (the Whitham equations) that admit a global solution. An important innovation was developed by El [28] whereby the DSW’s trailing and leading edge speeds could be determined without solving the full set of modulation equations. The Whitham-El DSW construction relies on the existence of a simple wave solution to the full, strictly hyperbolic and genuinely nonlinear modulation equations, but does not require its complete determination, hence analytical results are available even for non-integrable equations.

In this work, the one-dimensional (1D) isentropic Euler equations are regularized by a class of third order dispersive terms, modeling several of the aforementioned physical systems. The time to breaking (gradient catastrophe) for smooth initial data is numerically found to fall within bounds predicted by the dispersionless Euler equations. In order to investigate dynamics post-breaking, the long time resolution of the Riemann problem is considered. The DSW locus relating upstream, downstream flow configurations and the DSW speeds for admissible weak shocks are determined explicitly for generic, third order dispersive perturbations. The results depend only upon the sign of the dispersion () and the general pressure law assumed. A fundamental assumption in the DSW construction is the existence of an integral curve (simple wave) of the Whitham modulation equations connecting the upstream and downstream states in an averaged sense. Explicit, verifiable criteria for the breakdown of the simple wave assumption are given. The regularization for large amplitude DSWs depends upon the particular form of the dispersion. Thus, DSWs are explicitly constructed for particular pressure laws and dispersive terms of physical origin including a generalized Nonlinear Schrödinger (gNLS) equation modeling superfluids and nonlinear optics. Comparisons with a dissipative regularization are presented in order to highlight the differences between viscous and dispersive shock waves.

The outline of this work is as follows. Section 2 presents the general dispersive Euler model and assumptions to be considered, followed by section 3 outlining the gNLS and other dispersive Eulerian fluid models. Background sections 4 and 5 review the theory of the hyperbolic, dispersionless system and the Whitham-El method of DSW construction, respectively. A detailed analysis of DSW admissibility criteria and the breakdown of the simple wave assumption are undertaken in section 6 followed by the complete characterization of admissible weak DSWs in section 7. The theory of the breaking time for the gNLS model is shown to agree with numerical computation in section 8. Large amplitude DSWs for the gNLS equation are studied in section 9. The manuscript is completed by conclusions and an appendix on the numerical methods utilized.

## 2 Dispersive Euler Equations and Assumptions

The 1D dispersive Euler equations considered in this work are, in non-dimensional form

 ρt+(ρu)x=0,(ρu)t+[ρu2+P(ρ)]x=[D(ρ,u)]x,−∞

where is a fluid density, is the velocity, and is the (conservative) dispersive term. Formally setting gives the hydrodynamic approximation, valid until gradient catastrophe when the dispersion acts to regularize the singular behavior. The dispersionless, hyperbolic, isentropic Euler equations are known as the -system whose weak solutions to the Riemann problem are well-known [29, 25]. Here, the long time behavior of the dispersively regularized Riemann problem is analyzed. By the formal rescaling , ,

 ρT+(ρu)X=0,(ρu)T+[ρu2+P(ρ)]X=ε2[D(ρ,u)]X,−∞

the long time () behavior of the dispersive Euler equations in the independent variables is recast as a small dispersion () problem. Due to the oscillatory nature of the small dispersion limit, it is necessarily a weak limit as shown rigorously by Lax, Levermore, and Venakides for the Korteweg-deVries equation (KdV) [30, 31, 32, 33]. In this work, the multiscale Whitham averaging technique will be used to study the behavior of the dispersive Euler equations (2.1) for sufficiently large time and long waves.

The Whitham-El DSW simple wave closure method [28] is used to construct DSWs under the following assumptions.

• (sound speed) The pressure law is a smooth, monotonically increasing function of , for so that the speed of sound

 c0=c(ρ0)≡√P′(ρ0),

is real and the local Mach number

 M0≡|u0|c(ρ0),

is well-defined. It will also be assumed that the pressure is convex for so that .

• (symmetries) Equations (2.1) admit the Galilean invariance:

 D(ρ,u−u0)(x−u0t,t)=D(ρ,u)(x,t),

 D(ρ,−u)(−x,t)=−D(ρ,u)(x,t), (2.2)

so that (2.1) are invariant with respect to , .

• (dispersive operator) The dispersive term is a differential operator with of second order in spatial and/or mixed partial derivatives such that the system (2.2) has the real-valued dispersion relation

 ω=u0k±ω0(k,ρ0), (2.3)

with two branches found by linearizing about the uniform background state , with small amplitude waves proportional to . The appropriate branch of the dispersion relation is fixed by the in (2.3) with for , . The dispersion relation has the long wave expansion

 ω0(k,ρ0)=c0k+μk3+o(k3),k→0,μ≠0. (2.4)

The sign of the dispersion is for . Using (2.4) and the convexity or concavity of as a function of , one finds

 sgn(ω0(k,ρ0)k)k=sgn∂2ω0∂k2(k,ρ0).

Therefore, positive dispersion corresponds to increasing phase and group velocities with increasing while negative dispersion leads to decreasing phase and group velocities.

• (Whitham averaging) Equations (2.1) are amenable to Whitham averaging whereby a DSW can be described by a slowly varying, single-phase traveling wave. This requires

1. The system possesses at least three conservation laws. The mass and momentum equations in (2.1) account for two. An additional conserved quantity is required.

2. There exists a four parameter family of periodic traveling waves parametrized by, for example, the wave amplitude , the wavenumber , the average density , and the average velocity limiting to a trigonometric wave for small amplitude and a solitary wave for small wavenumber. In the cases considered here, the periodic traveling wave manifests as a solution of the ordinary differential equation (ODE) where is smooth as it varies over three simple, real roots. Two roots coincide in the small amplitude and solitary wave limits.

• (Simple wave) The Whitham-El method requires the existence of a self-similar simple wave solution to the four Whitham modulation equations (the averaged conservation laws and the conservation of waves). For this, the modulation equations must be strictly hyperbolic and genuinely nonlinear.

Assumption A1 provides for a modulationally stable, hydrodynamic long wave limit. The symmetry assumptions in A2 are for convenience and could be neglected. As will be demonstrated in section 3, A3 is a reasonable restriction still allowing for a number of physically relevant dispersive fluid models. The assumptions in A4 and A5 allow for the application of the Whitham-El method. While the assumptions in A4 are usually verifiable, A5 is often assumed. Causes of the breakdown of assumptions A3 (unique dispersion sign) and A5 (genuine nonlinearity) are identified and associated with extrema in the DSW speeds as either the left or right density is varied.

The nonstationary DSW considered here is the long time resolution of an initial jump in the fluid density and velocity, the Riemann problem

 u(x,0)={u1x<0u2x>0,ρ(x,0)={ρ1x<0ρ2x>0, (2.5)

where , .

## 3 Example Dispersive Fluids

The dispersive Euler equations (2.2) model a number of dispersive fluids including, among others, superfluids and optical fluids. The particular model equations described below were chosen because they incorporate different pressure laws and allow for different signs of the dispersion, key distinguishing features of Eulerian dispersive fluids and their weak dispersion regularization.

### 3.1 gNLS Equation

The generalized, defocusing nonlinear Schrödinger equation

 \rmiψt=−12ψxx+f(|ψ|2)ψ,f(0)=0,f(ρ)>0,ρ>0, (3.1)

or gNLS, describes a number of physical systems. For example, the “polytropic superfluid”

 f(ρ)=ρp,p>0, (3.2)

corresponds to the cubic NLS when that describes a repulsive BEC and intense laser propagation through optically defocusing (normal dispersion) media. The model (3.2) with describes a zero temperature Fermi gas near unitarity [34, 35] which is of special significance as recent experiments have been successfully interpreted with both dissipative [36] and dispersive [37] regularizations. Moreover, the regime describes the so-called BEC-BCS transition in ultracold Fermi gases [38]. The quintic NLS case, , models three-body interactions in a BEC [39, 40]. A BEC confined to a cigar shaped trap exhibits effective 1D behavior that is well-described by the non-polynomial nonlinearity [41, 42]

 f(ρ)=2√1+γρ−2γ,γ>0, (3.3)

here scaled so that , . In spatial nonlinear optics, photorefractive media corresponding to [43, 44]

 f(ρ)=ρ1+γρ,γ>0, (3.4)

is of particular interest due to recent experiments exhibiting DSWs [6, 8, 7, 45, 46]. For , the leading order behavior of (3.3) and (3.4) correspond to the cubic NLS.

The complex wavefunction can be interpreted in the dispersive fluid context by use of the Madelung transformation [47]

 ψ=√ρ\rme\rmiϕ,u=ϕx. (3.5)

Using (3.5) in (3.1) and equating real and imaginary parts results in the dispersive Euler equations (2.1) with

 P(ρ)=∫ρ0~ρf′(~ρ)\rmd~ρ,c(ρ)=√ρf′(ρ),x=14[ρ(lnρ)xx]x=ρ2[(√ρ)xx√ρ]x. (3.6)

The dispersive regularization of (2.1) corresponds to the semi-classical limit of (3.1), which, in dimensional units corresponds to for quantum many body systems. In applications, the dispersive regularization coincides with a strongly interacting BEC or a large input optical intensity.

Assumption A1 restricts the admissible nonlinearity to those satisfying

 f′(ρ)>0,(ρf′(ρ))′>0,ρ>0, (3.7)

which is realized by (3.2), (3.3) generally and for (3.4) when . Assumptions in A2 are well-known properties of the gNLS equation [48]. Assumption A3 is clear from (3.6) and the dispersion relation is

 ω0(k,ρ)=k√c2+k2/4∼ck+18ck3,|k|≪c. (3.8)

The dispersion is positive because for , .

Inserting the traveling wave ansatz

 ρ=ρ(x−Vt),u=u(x−Vt),

into (2.1) with (3.6) and integrating twice leads to

 u=V+Aρ, (3.9a) (ρ′)2=8[ρ∫ρρ1f(~ρ)\rmd~ρ+Bρ2+Cρ−A22]≡G(ρ). (3.9b)

It is assumed that has three real roots related to the integration constants , , and so that, according to a phase plane analysis, a periodic wave exists with maximum and minimum densities and , respectively. The fourth arbitrary constant, due to Galilean invariance, is the wave speed . In addition to mass and momentum conservation, an additional energy conservation law exists [49] which reads

 E≡ρu22+ρ2x8ρ+∫ρ0f(~ρ)\rmd~ρ,Et+{u[E+P(ρ)]}x=14[uρxx−(ρu)xρxρ]x,

hence the assumptions in A4 are satisfied. The hyperbolicity of the Whitham equations can only be determined by their direct study. The genuine nonlinearity of the system will be discussed in section 6. It will be helpful to note the solitary wave amplitude/speed relation which results from the boundary conditions for a depression (dark) solitary wave

 u0≡lim|ξ|→∞u(ξ),ρ0≡lim|x|→∞ρ(ξ),ρmin≡minξ∈Rρ(ξ).

A phase plane analysis of (3.9b) implies that the roots of satisfy , resulting in the solitary wave speed satisfying

 (s−u0)2=2ρmin(ρ0−ρmin)2[(ρ0−ρmin)f(ρ0)−∫ρ0ρminf(~ρ)\rmd~ρ]. (3.9j)

The soliton profile can be determined by integration of (3.9b).

Dispersive shock waves for the gNLS equation have been studied for the pure NLS case [50, 51] as well as in 1D photorefractive media [52] and the cubic-quintic case [53, 54]. A general DSW analysis will be presented in section 9.

### 3.2 Other Systems

The gNLS equation exhibits positive dispersion. Two additional examples are briefly given here with negative dispersion.

Two-temperature collisionless plasma: The dynamics of the ionic component of a two-temperature unmagnetized plasma [55] satisfy the dispersive Euler equations with

 P(ρ)=ρ,c(ρ)=1,D(ρ,u)=12ϕ2x−ϕxx,−ϕxx=ρ−\rmeϕ.

The electronic potential introduces nonlocal dispersion with dispersion relation

 ω0(k,ρ)=k√1+k2/ρ∼k−k32ρ0,|k|≪1.

It can be shown that , thus the system exhibits negative dispersion.

This system has been analyzed in detail [50] and satisfies assumptions A1-A4. Large amplitude dispersive shock waves were constructed in [28] under the assumptions of A5.

Fully nonlinear shallow water: Shallow waves in an ideal fluid with no restriction on amplitude satisfy the generalized Serre equations (also referred to as the Su-Gardner or Green-Naghdi equations) [56, 57, 58, 59] with

 Missing or unrecognized delimiter for \right (3.9k)

The density corresponds to the free surface fluid height and is the vertically averaged horizontal fluid velocity. The Bond number is proportional to the coefficient of surface tension. The dispersion relation is

 ω0(k,ρ)=k(ρ1+σk21+ρ2k2/3)1/2∼√ρ(k+3σ−ρ26k3),k→0.

The sign of the dispersion changes when corresponding to the critical values

 σcr=ρ23orkcr=1ρ(3+3√1+ρ2/σ)1/2.

The critical value expresses the fact that shallow water waves with weak surface tension effects, , exhibit negative dispersion for sufficiently long waves () and support elevation solitary wave solutions. Strong surface tension, , corresponds to positive dispersion and can yield depression solitary waves. Assumptions A1-A4 hold [59]. DSWs in the case of zero surface tension were studied in [60].

The Serre equations (3.9k) with and a model of liquid containing small gas bubbles [61] can be cast in Lagrangian form to fit into the framework of “continua with memory” [62]. The Whitham modulation equations for these dispersive Eulerian fluids were studied in [63]. Explicit, sufficient conditions for hyperbolicity of the modulation equations were derived.

The properties of DSWs for these systems will be discussed in section 10.

## 4 Background: Dispersionless Limit

The analysis of DSWs for (2.1) requires an understanding of the dispersionless limit

 ρt+(ρu)x=0,(ρu)t+[ρu2+P(ρ)]x=0, (3.9a)

corresponding to . Equations (3.9a) are the equations of compressible, isentropic gas dynamics with pressure law [64]. They are hyperbolic and diagonalized by the Riemann invariants (see e.g. [65])

 r1=u−∫ρc(ρ′)ρ′\rmdρ′,r2=u+∫ρc(ρ′)ρ′\rmdρ′, (3.9b)

with the characteristic velocities

 λ1=u−c(ρ),λ2=u+c(ρ), (3.9c)

so that

 ∂rj∂t+λj∂rj∂x=0,j=1,2. (3.9d)

By monotonicity of

 g(ρ)=∫ρc(ρ′)ρ′\rmdρ′,

the inversion of (3.9b) is achieved via

 u=12(r1+r2),ρ=g−1(12(r2−r1)). (3.9e)

In what follows, an overview of the properties of equations (3.9a) is provided for both the required analysis of DSWs and for the comparison of classical and dispersive shock waves.

### 4.1 Breaking Time

Smooth initial data may develop a singularity in finite time. The existence of Riemann invariants (3.9b) allows for estimates of the breaking time at which this occurs. In what follows, Lax’s breaking time estimates [66] are applied to the system (3.9d) with smooth initial data.

Lax’s general approach for hyperbolic systems is to reduce the Riemann invariant system (3.9d) to the equation

 z′=−a(t)z2,z(0)=m, (3.9f)

along a characteristic family, , and then bound the breaking time by comparison with an autonomous equation via estimates for and in terms of initial data for , .

Following Lax [66], integration along the 1-characteristic family in (3.9f) leads to , and satisfies

 ∂h∂r2=∂λ1∂r2λ1−λ2.

By direct computation with eqs. (3.9b), (3.9c), one can verify the following

 h=12ln[c(ρ)ρ], (3.9ga) a=\rme−h∂λ1∂r1=c(ρ)+ρc′(ρ)2c(ρ)[ρc(ρ)]1/2, (3.9gb) z=\rmeh∂r1∂x=∂r1∂x[c(ρ)ρ]1/2. (3.9gc)

The initial data for and are assumed to be smooth and bounded so that they satisfy

 r1––≤r1(x,t)≤¯¯¯¯¯r1,r2––≤r2(x,t)≤¯¯¯¯¯r2,0≤t

Assuming (non-vacuum conditions), then (3.9gb) implies and is decreasing along the 1-characteristic. Bounds for are defined as follows

 A=minρ∈RAa,B=maxρ∈RBa, (3.9gi)

where and are intervals related to the bounds on the initial data, chosen shortly. The initial condition is chosen as negative as possible

 x0=argminx∈Rz(0)=argminx∈R[∂u∂x−c(ρ)ρ∂ρ∂x][c(ρ)ρ]1/2∣∣ ∣∣t=0, (3.9gja) m=minx∈Rz(0)=[∂u∂x−c(ρ)ρ∂ρ∂x][c(ρ)ρ]1/2∣∣ ∣∣(x,t)=(x0,0). (3.9gjb)

These estimates lead to the following bounds on the breaking time

 −1mB≤tbr≤−1mA. (3.9gjk)

It is still necessary to provide the intervals and . The possible values of and in (3.9gh) and the monotonicity of the transformation for in (3.9e) suggests taking the full range of possible values . However, this choice does not provide the sharpest estimates in (3.9gjk). The idea is to use the fact that is constant along 1-characteristics. The choice for in (3.9gjb) suggests taking and allowing to vary across its range of values. While the optimal is associated with this characteristic, it does not necessarily provide the optimal estimates for or . A calculation shows

 ∂a∂r1=−18c3(ρc)1/2(c2−4ρcc′−2ρ2cc′′+3ρ2c′2).

It can be verified that for the example dispersive fluids considered here. In this case, any characteristic with can cause to increase, leading to a larger and a tighter bound in (3.9gjk). If , then may decrease, leading to a smaller and a tighter bound on . Combining these deductions leads to the choices

 RA=[g−1(r2––−r1(x0,0)2),g−1(¯¯¯¯¯r2−r1––2)], (3.9gjla) RB=[g−1(r2––−¯¯¯¯¯r12),g−1(¯¯¯¯¯r2−r1(x0,0)2)], (3.9gjlb)

when .

In summary, given initial data satisfying (3.9gh), the point and are determined from (3.9gja) and (3.9gjb). If , then there is no breaking. Otherwise, after verifying , the sets and are defined via (3.9gjla) and (3.9gjlb) leading to and in (3.9gi). The breaking time bounds are given by (3.9gjk). A similar argument integrating along the 2-characteristic field yields another estimate for the breaking time . The only changes are in (3.9gja) and (3.9gjb) where the minus sign goes to a plus sign and the choices for and reflect . These results will be used to estimate breaking times for dispersive fluids in section 8.

### 4.2 Viscous Shock Waves

It will be interesting to contrast the behavior of dispersive shock waves for (2.1) with that of classical, viscous shock waves resulting from a dissipative regularization of the dispersionless equations. For this, the jump and entropy conditions for shocks are summarized below [65, 25].

The Riemann problem (2.5) for (3.9a) results in the Hugoniot locus of classical shock solutions

 u2=u1±{[P(ρ2)−P(ρ1)](ρ2−ρ1)ρ1ρ2}1/2. (3.9gjlm)

The () corresponds to an admissible 1-shock (2-shock) satisfying the Lax entropy conditions when the characteristic velocity () decreases across the shock so that (). Weak 1-shocks connecting the densities and , exhibit the shock speed

 v(1)∼u1−c1−12(c1ρ1+c′1)Δ,0<Δ≪ρ1. (3.9gjln)

While the Riemann invariant exhibits a jump across the 1-shock, it is third order in so is approximately conserved for weak shocks. Weak, steady (non-propagating) 1-shocks satisfy the jump conditions

 Δ∼2ρ1c1c1+ρ1c′1(M1−1),M2∼1−(M1−1),0

where are the Mach numbers of the up/downstream flows and . The upstream flow indexed by is supersonic and the downstream flow is subsonic, this behavior also holding for arbitrary amplitude shocks.

### 4.3 Rarefaction Waves

Centered rarefaction wave solutions of (3.9a) exhibit the following wave curves connecting the left and right states

 1−rarefaction:u1 =u2+∫ρ2ρ1c(ρ)ρ\rmdρ,ρ1>ρ2, (3.9gjlpa) 2−rarefaction:u1 =u2−∫ρ2ρ1c(ρ)ρ\rmdρ,ρ2>ρ1, (3.9gjlpb)

where admissibility is opposite to the shock wave case. The characteristic velocities increase across a rarefaction wave. Since rarefaction waves are continuous and do not involve breaking, the leading order behavior of dispersive and dissipative regularizations for (3.9a) are the same. A dispersive regularization of KdV [67, 68] shows the development of small amplitude oscillations for the first order singularities at either the left or right edge of the rarefaction wave with one decaying as and the other . The width of these oscillations expands as [23] so that their extent vanishes relative to the rarefaction wave expansion with .

### 4.4 Shock Tube Problem

Recall that the general solution of the Riemann problem consists of three constant states connected by two waves, each either a rarefaction or shock [24, 25]. The shock tube problem [64] involves a jump in density for a quiescent fluid . The solution consists of a shock and rarefaction connected by a constant, intermediate state . For the case , a 1-shock connects to a 2-rarefaction via the Hugoniot locus (3.9gjlm) (with ) and the wave curve (3.9gjlpb), respectively. For example, a polytropic gas with gives the two equations

 1−shock:um=−[(κργm−κργ1)(ρm−ρ1)ρmρ1]1/2,2−rarefaction:um=−2(κγ)1/2γ−3[ρ(γ−1)/22−ρ(γ−1)/2m]. (3.9gjlpq)

Equating these two expressions provides an equation for the intermediate density and then the intermediate velocity follows.

## 5 Background: Simple DSWs

The long time behavior of a DSW for the dispersive Euler model (2.1) was first considered in [28]. In this section, the general Whitham-El construction of a simple wave led DSW for step initial data is reviewed. This introduces necessary notation and background that will be used in the latter sections of this work.

Analogous to the terminology for classical shocks, a 1-DSW is associated with the characteristic family of the dispersionless system (3.9a) involving left-going waves. In this case, the DSW leading edge is defined to be the leftmost (most negative) edge whereas the DSW trailing edge is the rightmost edge, these roles being reversed for the 2-DSW associated with the characteristic family. There is a notion of polarity associated with a DSW corresponding to its limiting behavior at the leading and trailing edges. The edge where the amplitude of the DSW oscillations vanish (the harmonic limit) is called the linear wave edge. The soliton edge is associated with the phase boundary where the DSW wavenumber (the soliton limit). Thus, the soliton edge could be the leading or trailing edge of the DSW, each case corresponding to a different DSW polarity. The polarity is generally determined by admissibility criteria and typically follows directly from the sign of the dispersion [50], as will be shown in section 6. The DSW construction for 1-DSWs is outlined below. A similar procedure holds for 2-DSWs.

Assuming the existence of a DSW oscillatory region described by slow modulations of the periodic traveling wave from assumption A4, three independent conservation laws are averaged with the periodic wave. The wave’s parameters , the average density, , the average velocity, , the generalized (nonlinear) wavenumber, and , the wave amplitude are assumed to vary slowly in space and time. The averaging procedure produces three first order, quasilinear PDEs. This set combined with the conservation of waves, ( here is the generalized, nonlinear frequency), following from consistency of wave modulations, results in a closed system for the modulation parameters, the Whitham modulation equations. As originally formulated by Gurevich and Pitaevskii [23], the DSW free boundary value problem is to solve the dispersionless equations (3.9a) outside the oscillatory region and match this behavior to the averaged variables and from the Whitham equations at the interfaces with the oscillatory region where (soliton edge) or (linear wave edge). These GP matching conditions correspond to the coalescence of two characteristics of the Whitham system at each edge of the DSW. Assumption A5 can be used to construct a self-similar, simple wave solution of the modulation equations connecting the soliton edge with the linear wave edge via an integral curve so that the two DSW boundaries asymptotically move with constant speed, the speeds of the double characteristics at each edge. In the Whitham-El method, the speeds are determined by the following key mathematical observations [28]

• The four Whitham equations admit exact reductions to quasi-linear systems of three equations in the (soliton edge) and (linear wave edge) regimes.

• Assuming a simple wave solution of the full Whitham equations, one can integrate across the DSW with explicit knowledge only of the reduced systems in the or planes of parameters, thereby obtaining the DSW leading and trailing edge speeds.

This DSW closure method is appealing because it bypasses the difficult determination and solution of the full Whitham equations. Furthermore, it applies to a large class of nonintegrable nonlinear wave equations. Some nonintegrable equations studied with this method include dispersive Euler equations like ion-acoustic plasmas [28], the Serre equations with zero surface tension [60], the gNLS equation with photorefractive [52] and cubic/quintic nonlinearity [69], and other equations including the Miyata-Camassa-Choi equations of two-layer fluids [70].

Simple DSWs are described by a simple wave solution of the Whitham modulation system which necessitates self-similar variation in only one characteristic field. Using a nontrivial backward characteristic argument, it has been shown that a simple wave solution requires the constancy of one of the Riemann invariants (3.9b) evaluated at the left and right states [28]. Then a necessary condition for a simple DSW is one of

 1−DSW:u2=u1−∫ρ2ρ1c(ρ)ρ\rmdρ,ρ2>ρ1, (3.9gjlpa) 2−DSW:u2=u1+∫ρ2ρ1c(ρ)ρ\rmdρ,ρ1>ρ2. (3.9gjlpb)

1-DSWs (2-DSWs) are associated with constant () hence vary in the () characteristic field. Equations (3.9gjlpa), (3.9gjlpb) can be termed DSW loci as they are the dispersive shock analogues of the Hugoniot loci (3.9gjlm) for classical shock waves. It is worth pointing out that the DSW loci correspond precisely to the rarefaction wave curves in (3.9gjlpa) and (3.9gjlpb). However, the admissibility criteria for DSWs correspond to inadmissible, compressive rarefaction waves where the dispersionless characteristic speed decreases across the DSW. The coincidence of rarefaction and shock curves does occur in classical hyperbolic systems but is restricted to a specific class, the so-called Temple systems [71] to which the dispersionless Euler equations do not belong.

Recall from section 4.2 that across a viscous shock, a Riemann invariant is conserved to third order in the jump height. Since the DSW loci (3.9gjlpa), (3.9gjlpb) result from a constant Riemann invariant across the DSW, the DSW loci are equal to the Hugoniot loci (3.9gjlm) up to third order in the jump height.

### 5.1 Linear Wave Edge

The integral curve of the Whitham equations in the (linear wave edge) plane of parameters reduces to the relationships , and the ODE

 \rmdk\rmd¯¯¯ρ=ω¯¯ρ¯¯¯u(¯¯¯ρ)−c(¯¯¯ρ)−ωk, (3.9gjlpq)

where the average velocity is constrained by the density through a generalization of (3.9gjlpa)

 ¯¯¯u(¯¯¯ρ)=u1−∫¯¯ρρ1c(ρ′)ρ′\rmdρ′. (3.9gjlpr)

The negative branch of the linear dispersion relation in (2.3), , is associated with left-going waves, hence a 1-DSW. Using (3.9gjlpr) and (2.3), (3.9gjlpq) simplifies to

 \rmdk\rmd¯¯¯ρ=ck/¯¯¯ρ+ω0¯ρc−ω0k. (3.9gjlps)

Equation (3.9gjlps) assumes , an exact reduction of the Whitham equations only at the linear wave edge. Global information associated with the simple wave solution of the full Whitham equations is obtained from the GP matching condition at the soliton edge by noting that the modulation variables satisfy for some depending on whether the soliton edge is leading or trailing, independent of the soliton amplitude . Thus (3.9gjlps) can be integrated in the plane with the initial condition , to associated with the linear wave edge, giving the wavenumber of the linear wave edge oscillations. This wavepacket’s speed is then determined from the group velocity .

Based on the Riemann data (2.5), the integration domain for (3.9gjlps) is . The initial condition occurs at either the leading edge where or the trailing edge where . The determination of the location of the linear wave edge, leading or trailing, is set by appropriate admissibility conditions discussed in section 6. The solution of (3.9gjlps) with initial condition at will be denoted so that one of

 k(ρj;ρj)=0,j=1,2, (3.9gjlpt)

holds for (3.9gjlps). Evaluating the solution of (3.9gjlps) at the linear wave edge , gives the wavenumber of the linear wavepacket at the leading (trailing) edge when (). The associated 1-DSW linear wave edge speed is denoted and is found from the group velocity evaluated at

 vj(ρ1,ρ2)=ωk[k(ρj;ρ3−j),ρj]=¯¯¯u(ρj)−ω0k[k(ρj;ρ3−j),ρj],j=1,2. (3.9gjlpu)

### 5.2 Soliton Edge

An exact description of the soliton edge where involves the dispersionless equations for the average density and velocity as well as an equation for the soliton amplitude . While, in principle, one can carry out the simple DSW analysis on these equations, it is very convenient to introduce a new variable called the conjugate wavenumber. It plays the role of an amplitude and depends on , , and thus is not a new variable. The resulting modulation equations at the soliton edge then take a correspondingly symmetric form relative to the linear edge. The speed of the soliton edge is determined in an analogous way to the linear edge by introducing a conjugate frequency

 ~ω(~k,¯¯¯ρ)=−\rmiω(\rmi~k,¯¯¯ρ)=¯¯¯u(¯¯¯ρ)~k+\rmiω0(\rmi~k,¯¯¯ρ)=¯¯¯u(¯¯¯ρ)~k+~ω0(~k,¯¯¯ρ). (3.9gjlpv)

The conjugate wavenumber plays the role analogous to an amplitude so that corresponds to the linear wave edge where . Integrating the ODEs for a simple wave in the plane results in

 \rmd~k\rmd¯¯¯ρ=c(¯¯¯ρ)~k/¯¯¯ρ+~ω0¯ρc(¯¯¯ρ)−~ω0~k, (3.9gjlpw)

the same equation as (3.9gjlpq) but with conjugate variables. It is remarkable that the description of the soliton edge so closely parallels that of the linear wave edge. The initial condition is given at the linear wave edge where . As in (3.9gjlpt), the solution with zero initial condition at is denoted according to

 ~k(ρj;ρj)=0,j=1,2. (3.9gjlpx)

Then the soliton speed , is the conjugate phase velocity evaluated at the conjugate wavenumber associated with the soliton edge

 sj(ρ1,ρ2)=~ω[~k(ρj;ρ3−j),ρj]~k(ρj;ρ3−j)=¯¯¯u(ρj)−~ω0[~k(ρj;ρ3−j),ρj]~k(ρj;ρ3−j). (3.9gjlpy)
###### Remark 1

In a number of example dispersive fluids studied in this work and elsewhere, the transformation to the scaled phase speed

 α(¯¯¯ρ)=ω0[k(¯¯¯ρ),¯¯¯ρ]c(¯¯¯ρ)k(¯¯¯ρ),

of the dependent variable in (3.9gjlps) and the analogous transformation

 ~α(¯¯¯ρ)=~ω0[~k(¯¯¯ρ),¯¯¯ρ]c(¯¯¯ρ)~k(¯¯¯ρ),

for (3.9gjlpw) are helpful, reducing the ODEs (3.9gjlps) and (3.9gjlpw) to simpler and, often, separable equations for and .

### 5.3 Dispersive Riemann Problem

The integral wave curves (3.9gjlpa), (3.9gjlpb) and the DSW loci (3.9gjlpa), (3.9gjlpb) can be used to solve the dispersive Riemann problem (2.5) just as the wave curves and the Hugoniot loci are used to solve the classical Riemann problem [25]. In both cases, the solution consists of two waves, one for each characteristic family, connected by a constant intermediate state. Each wave is either a rarefaction or a shock.

In contrast to the classical case, the integral wave curves and the DSW loci are the same for the dispersive Riemann problem. It is the direction in which they are traversed that determines admissibility of a rarefaction or a DSW. This enables a convenient, graphical description of solutions to the dispersive Riemann problem as shown in figure 2. Solid curves (——) correspond to example 1-wave curves (3.9gjlpa), (3.9gjlpa) and the dashed curves (-  -  -  -) correspond to example 2-wave curves (3.9gjlpb), (3.9gjlpb). The arrows provide the direction of increasing dispersionless characteristic speed for each wave family. Tracing an integral curve in the direction of increasing characteristic speed corresponds to an admissible rarefaction wave. The decreasing characteristic speed direction corresponds to an admissible DSW. The solution to a dispersive Riemann problem is depicted graphically by tracing appropriate integral curves to connect the left state with the right state . There are multiple paths connecting these two states but only one involves two admissible waves. This is shown by the thick curve in figure 2. Since the 1-wave curve is traced in the negative direction to the intermediate, constant state , this describes a 1-DSW. The 2-wave curve is then traced in the positive direction to the right state, describing a 2-rarefaction. The 1-DSW is admissible because . Since the characteristic speed is monotonically increasing, the 2-rarefaction is admissible (). The direction of the curve connecting the left and right states was taken from left to right. The opposite direction describes an inadmissible 1-rarefaction connected to an inadmissible 2-DSW.

The example shown in figure 2 corresponds to the dispersive shock tube problem consisting of an arbitrary jump in density for a quiescent fluid . Such problems have been studied in a number of dispersive fluids, e.g. [72, 60, 73, 70]. The determination of proceeds by requiring that the left state lie on the 1-DSW locus (3.9gjlpa)

 um=−∫ρmρ1c(ρ)ρ\rmdρ. (3.9gjlpz)

The second wave connects to the right state via the 2-rarefaction wave curve (3.9gjlpb)

 um=−∫ρ2ρmc(ρ)ρ\rmdρ. (3.9gjlpaa)

Equating (3.9gjlpz) and (3.9gjlpaa) leads to

 ∫ρmρ1c(ρ)ρ\rmdρ−∫ρ2ρmc(ρ)ρ\rmdρ=0,

which determines the intermediate density . The intermediate velocity follows.

This construction of the wave types and the intermediate state is independent of the sign of dispersion and the details of the dispersive term in (2.1), depending only upon the pressure law . For example, polytropic dispersive fluids with (e.g., gNLS with power law nonlinearity (3.2) and gSerre (3.9k)) yield the intermediate state (previously presented in [74])

 ρm=[12(ρ(γ−1)/21+ρ(γ−1)/22)]2/(γ−1),um=2(κγ)1/23−γ[ρ(γ−1)/21−ρ(γ−1)/22]. (3.9gjlpab)

This prediction will be compared with numerical computations of gNLS in section 9.2.