Discontinuous shock solutions of the Whitham modulation equations as dispersionless limits of traveling waves

TL;DR

This paper links discontinuous solutions of the Whitham modulation equations, called Whitham shocks, to zero dispersion limits of traveling waves in higher order dispersive PDEs, providing new insights into nonlinear wave behavior.

Contribution

It introduces a novel method to identify and construct admissible shock solutions of the Whitham equations as limits of traveling waves in higher order dispersive PDEs.

Findings

Whitham shocks correspond to zero dispersion limits of traveling waves.

Numerical simulations show heteroclinic traveling waves limit to Whitham shocks.

Admissible shocks are used to solve the generalized Riemann problem.

Abstract

Whitham modulation theory describes the zero dispersion limit of nonlinear waves by a system of conservation laws for the parameters of modulated periodic traveling waves. Here, admissible, discontinuous, weak solutions of the Whitham modulation equations--termed Whitham shocks--are identified with zero dispersion limits of traveling wave solutions to higher order dispersive partial differential equations (PDEs). The far-field behavior of the traveling wave solutions satisfies the Rankine-Hugoniot jump conditions. Generally, the numerically computed traveling waves represent heteroclinic connections between two periodic orbits of an ordinary differential equation. The focus here is on the fifth order Korteweg-de Vries equation. The three admissible one-parameter families of Whitham shocks are used as solution components for the generalized Riemann problem of the Whitham modulation equations. Admissible KdV5-Whitham shocks are generally undercompressive, i.e., all characteristic families pass through the shock. The heteroclinic traveling waves that limit to admissible Whitham shocks are found to be ubiquitous in numerical simulations of smoothed step initial conditions for other higher order dispersive equations including the Kawahara equation (with third and fifth order dispersion) and a nonlocal model of weakly nonlinear gravity-capillary waves with full dispersion. Whitham shocks are linked to recent studies of nonlinear higher order dispersive waves in optics and ultracold atomic gases. The approach presented here provides a novel method for constructing new traveling wave solutions to dispersive nonlinear wave equations and a framework to identify physically relevant, admissible shock solutions of the Whitham modulation equations.

Tables

Table 1. Table 1: Three distinct Whitham shock loci and right ( + + ) characteristic velocities for ( 55 ) when the right periodic wave φ + subscript 𝜑 \varphi_{+} is approximated by a weakly nonlinear Stokes wave ( 19 ), ( 20 ).

${\overline{u}}_{-}$	$V$	$k_{+}$	${\lambda }_1^{(+)}$	${\lambda }_2^{(+)}$	${\lambda }_3^{(+)}$
0.2963	$-0.0628$	0.3936	-0.7271	0.1858	1.0416
0.3657	0.0119	0.4775	-0.3768	0.4809	0.5360
0.9724	0.4219	0.8082	0.0586	1.6775	2.5904

Table 2. Table 2: Three distinct Whitham shock loci and right characteristic velocities using numerically computed periodic traveling waves. Compare with the approximate loci in Table 1 .

	${\overline{u}}_{-}$	$V$	$k_{+}$	${\lambda }_1^{(+)}$	${\lambda }_2^{(+)}$	${\lambda }_3^{(+)}$
$▲$	0.2522	$-0.1133$	0.3173	$-0.1462$	$-0.0780$	0.4056
$🌑$	0.3479	$-0.0071$	0.4496	0.6071	$-0.0237-0.1194i$	$-0.0237+0.1194i$
$■$	0.9726	0.4213	0.8080	0.0623	1.6472	2.5723

Equations

$$u_t+u{u}_x+u_{xxx}=0.$$

$$u_t+u{u}_x+u_{xxxxx}=0.$$

$$U_T+U{U}_X+{ϵ}^4{U}_{XXXXX}=0,$$

$$u_t+u{u}_x+\alpha u_{xxx}+u_{xxxxx}=0,$$

$$u_t+u{u}_x+\mathcal{K}\ast u_x=0,\mathcal{K}(x)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\frac{\omega (k)}{k}{e}^{ikx}dk,$$

$$(K\ast u_x)(x,t)={\int }_{\mathbb{R}}K(x-y)u_y(y,t)dy.$$

$${\omega }_{\mathrm{ww}}(k)$$

$$\sim k+\frac{3B-1}{6}{k}^3+\frac{(19-15B(3B+2))}{360}{k}^5+\dots ,0<k\ll 1$$

$$u(x,0)=\{\begin{array}{ll}\Delta & x<0\\ 0& x\ge 0\end{array},\Delta >0$$

$$u(x,t)=f(\xi ),\xi =x-ct$$

$$-c{f}^{\prime }+f{f}^{\prime }+f^{(5)}=0.$$

$$f^{\prime \prime \prime \prime }+\frac{1}{2}{f}^2-cf-A/2=0,$$

$$f^{\prime \prime \prime }{f}^{\prime }-\frac{1}{2}{\left(f^{\prime \prime }\right)}^2+\frac{1}{6}{f}^3-\frac{c}{2}{f}^2-\frac{A}{2}f+B=0,$$

$${\mathbf{\Phi }}^{\prime }=L\mathbf{\Phi }+R(\mathbf{\Phi }),L=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ c& 0& 0& 0\end{array}\right),R(\mathbf{\Phi })=\left(\begin{array}{c}0\\ 0\\ 0\\ \frac{1}{2}(A-{\Phi }_1^2)\end{array}\right),$$

$${\mathbf{U}}_0^{(\pm )}={\left(\begin{array}{cccc}{U}_0^{(\pm )}& 0& 0& 0\end{array}\right)}^T,U_0^{(\pm )}=c\pm \sqrt{c^2+A},$$

$$L+\frac{\partial \mathbf{R}}{\partial \mathbf{U}}({\mathbf{U}}_0^{(\pm )})=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ c-U_0^{(\pm )}& 0& 0& 0\end{array}\right),$$

$$\begin{array}{rl}{\lambda }_1^{(\pm )}& =-{\mu }^{(\pm )},{\mathbf{r}}_1^{(\pm )}={\left(1,-{\mu }^{(\pm )},{\mu }^{(\pm )},-{\left({\mu }^{(\pm )}\right)}^3\right)}^{\mathrm{T}},\\ {\lambda }_2^{(\pm )}& =-i{\mu }^{(\pm )},{\mathbf{r}}_2^{(\pm )}={\left(1,-i{\mu }^{(\pm )},-{\mu }^{(\pm )},i{\left({\mu }^{(\pm )}\right)}^3\right)}^{\mathrm{T}},\\ {\lambda }_3^{(\pm )}& =i{\mu }^{(\pm )},{\mathbf{r}}_3^{(\pm )}={\left(1,{\mu }^{(\pm )},-{\mu }^{(\pm )},-i{\left({\mu }^{(\pm )}\right)}^3\right)}^{\mathrm{T}},\\ {\lambda }_4^{(\pm )}& ={\mu }^{(\pm )},{\mathbf{r}}_4^{(\pm )}={\left(1,{\mu }^{(\pm )},{\mu }^{(\pm )},{\left({\mu }^{(\pm )}\right)}^3\right)}^{\mathrm{T}},\end{array}$$

$$f(\xi )=\phi (\theta ;\overline{u},a,k),\phi (\theta +2\pi ;\overline{u},a,k)=\phi (\theta ;\overline{u},a,k),$$

$$k,$$

$$\overline{u}=\frac{1}{2\pi }{\int }_0^{2\pi }\phi (\theta )d\theta ,$$

$$a=\underset{\theta \in [0,\pi ]}{\max }\phi (\theta )-\underset{\theta \in [0,\pi ]}{\min }\phi (\theta ),$$

$$u(x,t)\to u(x-u_{0}t,t)+u_0,$$

$$u\to bu,x\to b^{-1/4}x,t\to b^{-5/4}t$$

$$\begin{array}{rl}\phi (\theta ;\overline{u},a,k)& =\overline{u}+a\tilde{\phi }(a^{1/4}\theta -a^{1/4}k\overline{u}t;a^{-1/4}k),\\ c(\overline{u},a,k)& =a\tilde{c}(a^{-1/4}k)+\overline{u},\end{array}$$

$$\phi =\overline{u}+a{\phi }_1+a^2{\phi }_2+\dots ,c=c_0+a^2{c}_2+\dots .$$

$$\phi (\theta ;\overline{u},a,k)$$

$$c(\overline{u},a,k)$$

$$(u{)}_t+{\left(\frac{1}{2}{u}^2+u_{xxxx}\right)}_x$$

$${\left(\frac{1}{2}{u}^2\right)}_t+{\left(\frac{1}{3}{u}^3+u{u}_{xxxx}-u_x{u}_{xxx}+\frac{1}{2}{u}_{xx}^2\right)}_x$$

$$k={\theta }_x=S_X,\omega =-{\theta }_t=-S_T$$

$$u(x,t)=\phi (\theta ;\overline{u},a,k)+ϵu_1(\theta ,X,T)+\mathcal{O}({ϵ}^2),$$

$$\frac{\partial }{\partial x}=k\frac{\partial }{\partial \theta }+ϵ\frac{\partial }{\partial X},\frac{\partial }{\partial t}=-\omega \frac{\partial }{\partial \theta }+ϵ\frac{\partial }{\partial T}.$$

$$\overline{F}=\frac{1}{2\pi }{\int }_0^{2\pi }Fd\theta .$$

Full text

• July 2019

1 Introduction

Nonlinear dispersive wave equations are host to numerous solutions ranging from steady periodic traveling waves and solitary waves to unsteady dispersive shock waves (DSWs). The canonical model equation that includes third order (long wave) dispersion is the Korteweg-de Vries (KdV) equation

[TABLE]

In their seminal paper, Gurevich and Pitaevskii (GP) [1] construct the KdV DSW by solving a Riemann problem consisting of step initial conditions for the KdV-Whitham modulation equations. The KdV-Whitham modulation equations are a first order, quasi-linear system of three conservation laws for the slow evolution of a KdV nonlinear periodic traveling wave’s parameters [2]. Despite the apparent contradiction of applying the Whitham modulation equations—obtained by averaging conservation laws of the KdV equation under the assumption of long wavelength modulations—to discontinuous initial data, the GP DSW modulation solution is not a discontinuous shock but rather consists of a self-similar, centered rarefaction wave solution (a 2-wave [3]) between two constant levels. In a series of papers, Lax, Levermore [4, 5, 6, 7] and Venakedis [8, 9] utilized the inverse scattering transform to prove that the KdV-Whitham system (and their multiphase generalizations [10]) describe the zero dispersion limit of the KdV equation. The limit was shown to be weak in the $L_{2}$ sense and thus provides a rigorous justification for Whitham’s approach.

Because the KdV-Whitham equations are strictly hyperbolic [11], a generic class of initial data will necessarily lead to singularity formation and multivalued solution regions. Singularity formation in the KdV-Whitham equations is regularized by adding an additional phase to the multiphase modulation solution. The simplest case is breaking of the dispersionless Hopf equation (the zero-phase KdV-Whitham equation) for the mean flow that is regularized by the addition of an oscillatory region that emanates from a point of gradient catastrophe in the $x$-$t$ plane and is described by the one-phase KdV-Whitham equations for modulated periodic waves [12, 13, 14]. The boundaries of this modulated, one-phase region matches continuously to the zero-phase solution, which generalizes the GP step problem to smooth initial conditions. Singularity formation in the one-phase Whitham equations similarly leads to $x$-$t$ regions that are described by the two-phase Whitham equations [10] and match continuously to the solution of the one-phase Whitham equations [15]. Generally, compressive waves in the KdV-Whitham equations that lead to gradient catastrophe are regularized in terms of global, continuous expansion waves. This is also the case for other integrable systems such as the defocusing nonlinear Schrödinger equation [16] so that one can regularize singularity formation in the Whitham equations by an appropriate choice of degenerate initial conditions for sufficiently large phase modulations that yield global, continuous solutions [17, 18, 19].

Despite our understanding of singularity formation in the Whitham equations for integrable systems, there remains an outstanding question regarding discontinuous shock solutions. Do discontinuous, weak solutions of the Whitham equations have any meaning vis-a-vis the governing partial differential equation (PDE) and the physical system that it models? One standard approach in the conservation law community is to identify discontinuous shock solutions of hyperbolic systems as admissible if they are the pointwise limit of traveling wave (TW) solutions to the dissipatively or dissipative-dispersively perturbed conservation law as dictated by the microscopic physics of the problem [20, 21]. The TW profiles are heteroclinic, equilibrium-to-equilibrium solutions of a stationary ordinary differential equation (ODE). This approach leads to both admissible compressive Lax shocks [22, 23] and admissible nonclassical undercompressive shocks [24]. But the Whitham modulation equations are the zero dispersion limit of a conservative PDE so that there is no justification for introducing dissipative perturbations to the Whitham equations.

Instead, we expand the collection of TWs and prove that if TW solutions of a conservative PDE consist of heteroclinic connections between an equilibrium or a periodic orbit and another periodic orbit exist, their far-field behavior satisfies the Rankine-Hugoniot relations for discontinuous, shock solutions of the Whitham modulation equations. This enables us to define admissible shock solutions of the Whitham modulation equations—Whitham shocks—as (weak) zero dispersion limits of oscillatory, heteroclinic TW solutions. The solution of the governing PDE corresponding to a Whitham shock consists of two disparate oscillatory waves that are connected by a transition region occurring over the length scale of a single oscillation, the dispersive coherence length [3].

A phase plane analysis shows that the KdV equation (1) does not admit heteroclinic TW solutions [25]. Consequently, the Whitham shocks studied here are not admissible solutions of the KdV-Whitham equations. Instead, we consider a different class of nonlinear dispersive equations.

While KdV is a universal model of DSWs in convex media [3], modifications to it are required when higher order, e.g., short wavelength, or large amplitude, effects occur. The KdV equation is then modified to include higher order dispersive and/or nonlinear terms that more accurately model the underlying physics. This manuscript focuses on the implications of higher order and nonlocal dispersive effects on both steady, traveling wave solutions as well as unsteady DSWs. A canonical model incorporating higher order dispersion and weak nonlinearity is the fifth order Korteweg-de Vries equation (KdV5)

[TABLE]

The dispersionless limit of the KdV5 equation (2) which will be of use throughout this manuscript is made explicit via the hydrodynamic variable transformation $X=\epsilon x$, $T=\epsilon t$, and $U(X,T;\epsilon)=u(X/\epsilon,T/\epsilon)$. Substitution of these variables into (2) yields

[TABLE]

so that the singular dispersionless limit is achieved by passing $\epsilon\to 0$.

In a typical multi-scale asymptotic expansion for weakly nonlinear, long waves, the dispersion relation expansion for small wavenumber involves successively smaller terms. In this scenario, cubic dispersion typically dominates, which is the reason the KdV equation is said to be universal. However, when successive terms are comparable, e.g., third and fifth order dispersion near an inflection point [26], then higher order dispersion is operable. The KdV5 equation (2) is an example in which fifth order dispersion dominates over all others, but other well known examples also exist. A long wave, nonlinear model that incorporates the competition between third and fifth order dispersion is the Kawahara equation [27]

[TABLE]

where $\alpha\in\mathbb{R}$ is a parameter. Higher order dispersive effects can also occur when full, or nonlocal dispersion is included. A class of models in this vein is

[TABLE]

where $\omega(k)$ is the linear dispersion relation and

[TABLE]

This model was originally proposed by Whitham as a weakly nonlinear, fully dispersive model of water waves [28, 23] and has since been called the Whitham equation, not to be mistaken for the aforementioned Whitham modulation equations. Since then, this model has been used in applications to understand the propagation of water waves with surface tension [29]. The dispersion relation for gravity-capillary water waves is

[TABLE]

so that the Fourier transform (5) and convolution (6) should be considered in the distributional sense [30, 31]. For Bond number $B$ (the ratio of surface tension to gravity forces) near $1/3$, higher order dispersion is important. Generally, whenever linear dispersion exhibits an inflection point for nonzero wavenumber $k$, higher order dispersive effects take precedence [26].

Dispersive systems with higher order dispersion occur in a variety of applications, for example, gravity-capillary water waves [32, 33, 34] and flexural ice sheets [35, 36, 37], nonlinear optical systems [38, 39, 40, 41], and Bose-Einstein condensates [42] all exhibit a change in dispersion curvature for sufficiently short waves. Nonlinear PDEs with higher order or nonlocal dispersion admit solutions wholly different from their lower order dispersive counterparts. Examples include multi-mode periodic waves (traveling waves in which two wavenumbers are resonant) [43, 44, 45], solitary waves consisting of multiple pulses [46], solitary waves with decaying oscillatory tails [27, 47], and solitary waves accompanied by co-propagating small, but finite, amplitude oscillations spatially extending to infinity [48, 49, 50].This rich zoology of traveling waves (TWs) can be attributed to additional degrees of freedom inherent in the higher order ODE that results from a TW ansatz. In this paper, we identify new TW solutions to further expand this zoology.

The increased variety of solutions to higher order dispersive PDEs is not limited to steady solutions. The unsteady dynamics that result from dispersively regularized gradient catastrophe, most simply embodied by the step initial data

[TABLE]

for, e.g., equations (2) and (4), are notably different from lower order dispersive systems such as KdV (1). Previous numerical and asymptotic studies have jointly identified the following distinct, unsteady regimes for Riemann problems of higher order or nonlocal PDEs:

DSW implosion: above a critical jump height, the DSW’s harmonic wave edge becomes modulationally unstable and a multiphase structure emerges [51, 52]. 2. 2.

Radiating DSW (RDSW): a perturbed convex DSW with exponentially small radiation emitted from the solitary wave edge. This results from a linear resonance between the solitary wave and shorter linear waves [41, 26]. 3. 3.

Traveling DSW (TDSW): a rapid transition from an equilibrium level to a comoving nearly periodic wavetrain slowly transitions to a constant level via a partial DSW from the intermediate periodic wavetrain to a background constant [53, 26, 54]. See Figure 1 for an example TDSW solution of the KdV5 equation. 4. 4.

Crossover DSW: can be thought of as a combination of the RDSW and TDSW, where the linear resonance begins to grow to a finite amplitude wavetrain and the wave dynamics are more complex [41, 26].

The original motivation for this manuscript was to fully describe the TDSW solution in Fig. 1 in terms of modulation theory. As we show in Sec. 6, the TDSW can be described in terms of a shock-rarefaction solution to the KdV5-Whitham and the Kawahara-Whitham modulation equations. However, during our research, a broader theme emerged and revealed a host of novel heteroclinic and homoclinic TW solutions as well as a Whitham modulation theory framework to interpret them.

A recent preprint considers the so-called generalized Riemann problem for the Whitham modulation equations of the Serre-Greene-Naghdi (SGN) equations that model fully nonlinear shallow water waves [55]. In particular, the numerical evolution of initial conditions consisting of one periodic wave that rapidly transitions to another periodic wave was shown as evidence that discontinuous, weak solutions of the SGN-Whitham equations can be realized. Emerging from this class of initial data were structures that resemble genuine traveling wave solutions of the governing PDE. The dispersion relation for the SGN equations is $\omega(k)=k/\sqrt{k^{2}/3+1}$, which does not exhibit an inflection point. Consequently, the numerical results obtained in [55] and their interpretation as shock solutions of the Whitham equations suggest that Whitham shocks persist in nonlocal, strongly nonlinear regimes in addition to the weakly nonlinear cases studied here.

This paper is organized as follows. In Section 2, we discuss TW solutions of the KdV5 equation. Homoclinic TW solutions are used to derive the KdV5-Whitham modulation equations in Section 3. We also study the properties of the Whitham modulation equations and identify a class of 2-wave curves that we use to solve the KdV5 equation (2) with Riemann data (7). In Section 4, we study shock solutions of the modulation equations. Here, we also prove that heteroclinic TWs consisting of disparate, co-propagating periodic waves necessarily satisfy the Rankine-Hugoniot jump conditions for the Whitham equations. We use this result to define admissible Whitham shocks as the zero dispersion limit of KdV5 heteroclinic TW solutions. In Section 5, we use the Whitham shock locus to compute a rich variety of heteroclinic and homoclinic TW solutions. We discuss extensions and applications of the newly developed theory to the Riemann problem (7) for the KdV5 equation (2), the Kawahara equation (4), and the Whitham equation (5) in Sec. 6. In Section 7, we conclude the manuscript and postulate new and related problems to the present work.

2 KdV5 periodic traveling wave solutions

Traveling waves are sought in the form

[TABLE]

where $c$ is the phase velocity. The profile function $f$ satisfies the ODE

[TABLE]

Integration yields the fourth order profile ODE

[TABLE]

where $A\in\mathbb{R}$ is a constant of integration. This ODE can be integrated again to yield the energy-type integral

[TABLE]

where $B\in\mathbb{R}$ is again a constant of integration. We can rewrite the ODE (9) as the first order system

[TABLE]

where $\Phi_{n}=\mathrm{d}^{n-1}f/\mathrm{d}\xi^{n-1}$, $n=1,2,3,4$, are the components of the vector $\mathbf{\Phi}$.

2.1 Linearization

The fixed points $\mathbf{U}(\xi)=\mathbf{U}_{0}$ of the ODE (11) are

[TABLE]

so that $c^{2}>-A$ for real equilibria. Linearization about these fixed points results in the matrices

[TABLE]

whose eigenvalues $\lambda_{j}^{(\pm)}$ come in quartets with corresponding right eigenvectors $\mathbf{r}_{j}^{(\pm)}$

[TABLE]

where $\mu^{(\pm)}=(c-U_{0}^{(\pm)})^{1/4}$. We observe that the fixed point $\mathbf{U}_{0}^{(-)}$ has $U_{0}^{(-)}=c-\sqrt{c^{2}+A}<c$ so that the complex conjugate eigenvalues $\lambda_{2,3}^{(-)}=\mp i(c^{2}+A)^{1/8}$ lie on the imaginary axis whereas $\lambda_{1,4}^{(-)}=\mp(c^{2}+A)^{1/8}$ are opposite and lie on the real axis. The eigenvalues $\lambda_{j}^{(+)}$ are $\pi/4$ rotations of $\lambda_{j}^{(-)}$ in the complex plane for $j=1,2,3,4$.

This linear analysis suggests the existence of small amplitude periodic orbits near $\mathbf{U}_{0}^{(-)}$, with one-dimensional stable and unstable manifolds in the directions $\mathbf{r}_{1}^{(-)}$ and $\mathbf{r}_{4}^{(-)}$, respectively. These manifolds are what enable the periodic-heteroclinic-to-periodic TWs numerically computed in this paper. The existence of periodic orbits for (11) was demonstrated in [56] and related results establishing the existence and stability of wavetrains in the related Kawahara equation (4) are established in [57]. This motivates our investigation of heteroclinic connections between differing periodic orbits.

First, we consider $2\pi$-periodic traveling wave solutions to (2) in the form

[TABLE]

where $\theta=k\xi$ is the phase variable and $\varphi$ possesses three distinct parameters identified in Figure 2 and defined as

[TABLE]

and wave frequency $\omega=ck$, $c=c(\overline{u},a,k)$.

2.2 Scaling and Galilean symmetries

KdV5 (2) admits the Galilean

[TABLE]

and scaling

[TABLE]

symmetries. Therefore, the three-parameter family of periodic traveling waves can be generated from the one-parameter family $\left\{(\tilde{\varphi}(\theta;\tilde{k}),\tilde{c}(\tilde{k}))\right\}_{\tilde{k}\geq 0}$ by

[TABLE]

for any $a>0$, $k>0$, $\overline{u}\in\mathbb{R}$. Note that $\tilde{\varphi}(\theta;\tilde{k})=\varphi(\theta;0,1,\tilde{k})$ and $\tilde{c}(\tilde{k})=c(0,1,\tilde{k})$.

2.3 Approximate and numerical computations of periodic waves

We now obtain the periodic solutions to (8) via a weakly nonlinear Stokes frequency shift calculation [23] and through numerical computations. The periodic solutions will be used to obtain the Whitham modulation system.

To this end, we seek an approximate form for the periodic wave $\varphi(\theta)$ and its phase speed $c$ as series expansions in the small, finite amplitude parameter $a$

[TABLE]

Inserting the asymptotic approximation and equating like powers of $a$, we obtain the periodic wave and phase velocity up to $\mathcal{O}(a^{2})$

[TABLE]

Here, the restriction on the wavenumber $k\gg a^{1/4}$ is required so that the asymptotic series remains well ordered. We remark that the asymptotic series can be rescaled to arbitrary $a$ via the symmetries (16), (17) so that the requirement $a\ll 1$ can be formally relaxed but we must maintain the short wave requirement $k\gg a^{1/4}$ in order to respect asymptotic ordering.

We now present computations of periodic solutions to the profile ODE (8) where we take $A=0$ via a Galilean shift. The solutions are computed via a pesudospectral projection onto a Fourier basis. A solution to the nonlinear system for the periodic wave’s Fourier coefficients is carried out with a Newton iteration. Details of this computation are contained in Appendix A.

Examples of periodic traveling wave solutions and their corresponding nonlinear dispersion curve are given in Fig. 3. Small $\tilde{k}$ corresponds to a well separated, solitary wave train. One metric for the validity of the Stokes approximation is the accuracy of the phase velocity (20), which compares favorably with the numerically computed phase velocity for $\tilde{k}\gtrsim 0.4$.

3 KdV5-Whitham modulation equations

The Whitham modulation equations describe the dispersionless limit of nonlinear wavetrains. Equivalently, these equations describe slow modulations of periodic traveling waves. The modulation equations are a system of first order, quasi-linear conservation laws for the parameters $(\overline{u},a,k)$ that we now derive. Equation (2) admits the two conservation laws

[TABLE]

The Whitham equations can be obtained by a period-average of the conservation laws, evaluated on the periodic TW manifold [2]. To see this, we introduce the rapid phase variable $\theta=S(X,T)/\epsilon$ satisfying

[TABLE]

where we recall $X=\epsilon x$ and $T=\epsilon t$ are the long space and slow time scales with $0<\epsilon\ll 1$ as in (3). We then seek an asymptotic solution in the form

[TABLE]

with a leading order periodic traveling wave whose parameters $\overline{u}$, $a$, $k$ depend on the slow variables $X$ and $T$ and each of the $u_{n}$ are $2\pi$-periodic in $\theta$. Consequently, the derivatives in (21), (22) transform according to

[TABLE]

We now introduce averaging of the quantity $F[u(x,t);X,T]$ according to

[TABLE]

Then, substituting the asymptotic series (23) for $u$ into $F$ results in

[TABLE]

owing to the fact that the period average of an exact $\theta$ derivative is zero for periodic functions. Applying the averaging operator (24) to the conservation laws (21), (22) with the asymptotic series (23) for $u$, we obtain

[TABLE]

subject to $\mathcal{O}(\epsilon)$ corrections that vanish in the zero dispersion limit. The modulation system is closed by the phase compatibility condition $S_{XT}=S_{TX}$, yielding the conservation of waves

[TABLE]

Thus, we conclude that the KdV5-Whitham modulation equations (25), (26), and (27) describe the dispersionless limit of the KdV5 equation (3). If initial data for Eqs. (25), (26), and (27) are invariant to the hydrodynamic scaling transformation $X\to bX$, $T\to bT$—e.g., the step initial data considered later—then an equivalent alternative to the zero dispersion limit is the long time limit of Eqs. (25), (26), and (27) with the unscaled independent variables $x=X/\epsilon$, $t=T/\epsilon$.

The conservation laws (25), (26) can be simplified via integration by parts

[TABLE]

These identities, along with a return to the unscaled variables $x,t$, imply that the Whitham equations in conservative form can be written compactly

[TABLE]

This will be convenient of our purposes because it eliminates the small dispersion parameter $\epsilon$ from the problem.

3.1 Properties of the Whitham equations

3.1.1 Weakly nonlinear regime

We approximate the averaging integrals in the Whitham equations (28)–(30) with the Stokes wave (19) and its dispersion relation (20) yielding

[TABLE]

accurate to $\mathcal{O}(a^{2})$. Inserting these approximate integrals into the averaged conservation laws (28)–(30), the weakly nonlinear KdV5-Whitham equations in conservative form are approximately

[TABLE]

The non-conservative, quasi-linear form is

[TABLE]

A standard perturbation calculation gives the approximate eigenvalues

[TABLE]

and eigenvectors

[TABLE]

of the flux Jacobian matrix in (33). The quasi-linear system (33) is genuinely nonlinear if $\mu_{j}=\nabla\lambda_{j}\cdot\mathbf{v}_{j}\neq 0$ for all $j$. By inspection and direct calculation, we find that the weakly nonlinear Whitham equations in the asymptotic regime $0<a\ll 1$, $k\gg a^{1/4}$ are strictly hyperbolic and genuinely nonlinear.

3.1.2 Strongly nonlinear regime

In the strongly nonlinear regime, we rewrite the modulation equations (28)–(30) in non-conservative form in terms of the modulation variables $\mathbf{q}=[\overline{u},a,k]^{\rm T}$

[TABLE]

where

[TABLE]

which is valid provided the matrix inverse exists. For all of the computations performed here, $\mathcal{A}$ is well-defined. The symmetries (16), (17) of solutions to the KdV5 equation can be used to directly compute the integral dependence on the mean $\overline{u}$ and amplitude $a$. We define the averaging integrals

[TABLE]

These integrals can be written explicitly in terms of period averages over the one-parameter family of periodic solution $\tilde{\varphi}(\theta;\tilde{k})$ using (18)

[TABLE]

where we have introduced the $\tilde{\cdot}$ variables to denote averaging integrals of unit amplitude, zero mean periodic solutions.

We now investigate the hyperbolicity/ellipticity and genuine nonlinearity of the modulation equations. These criteria depend solely upon the eigenvalues $\lambda_{j}$ and eigenvectors $\mathbf{r}_{j}$ of $\mathcal{A}$ that satisfy

[TABLE]

where $I$ is the $3\times 3$ identity matrix. In general, we expect three eigenpairs $\left\{(\lambda_{j},\mathbf{r}_{j})\right\}_{j=1}^{3}$ that depend on the values $(\overline{u},a,k)$. By the symmetries (16), (17), it is enough to determine their dependence on $\tilde{k}$ only.

Weak hyperbolicity of the modulation equations is a necessary condition for the modulational stability of periodic waves [23, 58]. Numerical computation of the characteristic velocities depicted in Figure 4 illustrate hyperbolicity of the modulation equations for a range of $\tilde{k}$. When $0.41\lesssim\tilde{k}\lesssim 0.65$ and $0.25\lesssim\tilde{k}\lesssim 0.3$, two complex conjugate characteristic velocities with nonzero imaginary part emerge, so that corresponding periodic waves are modulationally unstable. It is further shown in Fig. 4 that the weakly nonlinear calculations for the characteristic velocities Eq. (34)–(36) are in excellent agreement with fully nonlinear calculations for $\tilde{k}\gtrsim 0.7$.

In a similar manner, we utilize the eigenvalues and eigenvectors to determine regions of genuine nonlinearity for the modulation equations (28)–(30) by computing the quantity

[TABLE]

For waves of unit amplitude and zero mean, we find points of non-genuine nonlinearity where $\mu_{1}\approx\mu_{2}\approx 0$, which are identified by the red points in the lower panel of Fig. 4 (a). Our careful numerical computations (see Appendix) suggest that both the first and second characteristic fields are not genuinely nonlinear at these select points $\tilde{k}\in\{0.167,0.185,0.222,0.370\}$.

In strictly hyperbolic regions, we compute wave curves that parameterize self-similar, simple wave solutions to the modulation equations (40). These wave curves will be applied to the solution of Riemann problems in Sec. 6. The $j$-th wave curve is the integral curve in the direction of $\mathbf{r}_{j}$ [21]. We introduce the self similar parameterization $s=s(x,t)$ so that $\mathbf{q}=\mathbf{q}(s)$ and, along $s=\lambda_{j}$,

[TABLE]

Illustrative 2-wave curves are displayed in Figs. 4(b) and (c) and correspond to the shaded regions in Fig. 4(a) (lower panel) on the left and right respectively. In these figures, we project the three dimensional curve $(\overline{u}(s),a(s),k(s))$ onto both the $\overline{u}$-$a$ plane and the $\overline{u}$-$k$ plane. We normalize the wave curves to emanate from the zero mean, unit amplitude state $\overline{u}=0$, $a=1$. For this visualization, we should note that $\overline{u}(s)$ is a monotonically decreasing function in $s$.

4 Whitham shocks: Abstract Setup

Now that we have obtained the Whitham modulation equations and described their structure, we consider shock solutions for the conservation laws (28)–(30) and their connections to heteroclinic traveling wave solutions of the KdV5 equation.

Shock solutions of the modulation equations (28)–(30) take the form of a moving discontinuity

[TABLE]

In order for the shock (44) to be a weak solution of (28)–(30), the left $(\overline{u}_{-},a_{-},k_{-})$, right $(\overline{u}_{+},a_{+},k_{+})$ states and the shock velocity $V$ satisfy the Rankine-Hugoniot jump conditions

[TABLE]

where $\mathbf{P},\mathbf{Q}\in\mathbb{R}^{3}$ are the averaged densities and fluxes in Eqs. (28)–(30) and $\llbracket\cdot\rrbracket$ represents the difference in the left, right values. Consequently, the jump conditions take the form

[TABLE]

where $\varphi_{\pm}=\varphi(\theta;\overline{u}_{\pm},a_{\pm},k_{\pm})$ are the right ($+$)/left ($-$) periodic traveling waves that compose the Whitham shock. We summarize this in the following.

Definition 1** (Whitham shock).**

A discontinuity (44) moving with velocity $V$ is called a Whitham shock when it is a weak solution of the Whitham modulation equations in conservative form (28)–(30). That is, the left $(\overline{u}_{-},a_{-},k_{-})$ and right $(\overline{u}_{+},a_{+},k_{+})$ states satisfy the Rankine-Hugoniot jump conditions (45)–(47).

In order to determine admissibility of Whitham shocks, we consider smooth heteroclinic TW solutions to the KdV5 equation (2) that asymptote to distinct periodic waves as $x\to\pm\infty$

[TABLE]

with each far-field periodic wave ($\varphi_{\pm}$) characterized by the parameters $(\overline{u}_{\pm},a_{\pm},k_{\pm})$ and phases $\theta_{\pm}$.

Our main result states that KdV5 heteroclinic traveling wave solutions exhibiting the far field behavior (48) necessarily lie on the Whitham shock locus of states satisfying (45)–(47).

Theorem 1**.**

Suppose $f(\xi)$ with speed $c$ ($\xi=x-ct$) satisfies (8), hence is a traveling wave solution of the differential equation (2) such that

[TABLE]

uniformly as $\xi\to\pm\infty$, then the Rankine-Hugoniot relations (45)–(47) for the KdV5-Whitham equations (28)–(30) are satisfied by the far-field periodic waves $\varphi_{\pm}$ with parameters $\overline{u}_{\pm}$, $a_{\pm}$, $k_{\pm}$. The traveling wave speed is the shock speed $c=V$ and coincides with the phase velocities $c=c_{\pm}$.

Proof.

By assumption, there exists some TW profile $f(\xi)$ that solves equation (8) with the far-field behavior described by the boundary conditions (48). The factors $\theta_{\pm}$ are appropriate phase shifts so that the traveling wave matches the far-field periodic orbits $\varphi_{\pm}$. We define the averaging operators $\bar{\mathcal{L}}^{\pm}$ acting on a periodic, integrable function $F(\xi)$ as

[TABLE]

Since $f\to\varphi_{\pm}$ uniformly, we can compute for $n=1,2,3$

[TABLE]

and

[TABLE]

Since $f$ is a traveling wave solution to (8) with two periodic wave limits, the phase speed of each must be the same, i.e.,

[TABLE]

This condition immediately implies the jump condition (47) with $V=c$ from the conservation of waves. All TW solutions with speed $c$ admit the first and second integrals (9) and (10). We now apply the operator $\bar{\mathcal{L}}^{\pm}$ to the first integral

[TABLE]

where we used $\bar{\mathcal{L}}^{\pm}[f^{\prime\prime\prime\prime}]=0$ by virtue of periodicity. Equating the two relations in (50) to eliminate $A$ gives

[TABLE]

which is the first jump condition (45) when we identify $V=c$. Applying the operator $\bar{\mathcal{L}}^{\pm}$ to the second integral (10) results in

[TABLE]

where we used $\bar{\mathcal{L}}^{\pm}[f^{\prime\prime\prime}f^{\prime}]=k_{\pm}^{4}\overline{\varphi_{\pm}^{\prime\prime\prime}\varphi_{\pm}^{\prime}}$. Integrating by parts and applying (51) simplifies Eq. (52) to

[TABLE]

The third and final Rankine-Hugoniot condition (46) with $V=c$ is found by subtracting the $+$ and $-$ instances of Eq. (53) to eliminate $B$

[TABLE]

thereby completing the proof. ∎

Theorem 1 motivates the following.

Definition 2** (admissibility).**

A KdV5 Whitham shock (recall Definition 1) is admissible if there exists a heteroclinic traveling wave solution with speed $c=V$ to the KdV5 equation (2) satisfying the profile equation (8) and the boundary conditions (48) corresponding to the left ($-$)/right ($+$) states of the Whitham shock (44) and $V=c_{+}=c_{-}$.

In the following section, we provide extensive numerical computations of heteroclinic TW solutions that support the existence of admissible Whitham shocks for the KdV5 equation.

5 Whitham shocks

In this section, we study admissible KdV5 Whitham shocks via increasing levels of complexity. First, we consider the case of a shock solution (44) to the modulation equations where one far-field state degenerates to zero wavenumber, i.e., a solitary wave. These results are then generalized to Whitham shocks in which two periodic waves satisfying the jump conditions propagate with identical phase velocities. The section culminates with computations of two co-propagating Whitham shocks where the corresponding TW solutions are homoclinic, localized oscillatory patterns on an oscillatory or uniform background.

5.1 Solitary wave to periodic

We consider the case where the left periodic wave degenerates to a solitary wave ($k_{-}\to 0$) and the right periodic wave is of unit amplitude and zero mean ($a_{+}=1$, $\overline{u}_{+}=0$). The associated Whitham shock (44) is

[TABLE]

In the solitary wave limit ($k\to 0$) of the modulation equations (28)–(30), the averaged quantities are

[TABLE]

Therefore, the jump conditions (45)–(47) with a solitary wave on the left are

[TABLE]

An illuminating calculation using the Stokes wave approximation (19) for $\varphi_{+}$ leads to explicit formulae. In this case, the jump conditions (56)–(58) are solved by

[TABLE]

where $\overline{u}_{-}^{2}$ is a root of the polynomial

[TABLE]

which has three positive solutions $\overline{u}_{-}^{2}\in\{0.08777,0.1337,0.9455\}$. The positive square roots are inserted into (59) to obtain the wavenumber $k_{+}$, shock velocity $V$, and the right characteristic velocities from Eqs. (34)–(36) for three distinct shock loci, all summarized in Table 1. The left solitary wave amplitude $a_{-}$ can be recovered from the solitary wave speed-amplitude relation by equating it to the shock velocity $c(\overline{u}_{-},a_{-},0)=V$. The left characteristic speeds for the left solitary wave state are $\lambda_{1}^{(-)}=\overline{u}_{-}$, $\lambda_{2}^{(-)}=\lambda_{3}^{(-)}=V$. The reason that two of the characteristic velocities for the left state are the same is that two of the three modulation equations coincide in the solitary wave limit $k\to 0$, which is a well-known property of the Whitham equations [3]. The cases where $\overline{u}_{-}^{2}<0$ or $\overline{u}_{-}<0$ in Eq. (60) can be dismissed because these choices result in a complex value for $\overline{u}_{-}$ or $k_{+}$ in Eq. (59).

We test the accuracy of the approximate Whitham shock loci reported in Table 1 by solving the jump conditions (56)–(58) with the family of periodic traveling wave solutions obtained numerically in Sec. 2.3. Results from the numerical computations are given in Table 2. The Whitham shock locus with the largest root of Eq. (60) is well-approximated to three digits (denoted by $\blacksquare$). The corresponding right characteristic velocities on this locus are accurate to one or two digits. The two Whitham shock loci corresponding to the two smaller roots of Eq. (60) are less accurately approximated. The reason for this is the Stokes approximation restriction $k_{+}\gg a_{+}^{1/4}=1$ that requires a sufficiently large wavenumber $k_{+}$. The Whitham shock locus with the largest root $\overline{u}_{-}$ also exhibits the largest wavenumber $k_{+}$, hence is expected to be a better approximation to the true Whitham shock locus, although good agreement is notable given that $k_{+}<1$. We remark that a numerical search did not yield any other Whitham shock loci.

Two of the characteristic velocities for the locus denoted $\CIRCLE$ are complex, therefore we expect this locus to correspond to unstable Whitham shocks. We investigate stability questions in Sec. 5.3.

We now compute heteroclinic traveling wave solutions whose zero dispersion limit correspond to a Whitham shock from each of the loci reported in Table 2, hence demonstrating Whitham shock admissibility. See the Appendix for computational details. The obtained solutions are depicted in Fig. 5. All three TW solutions visually look similar for $x<0$, which corresponds to the left solitary wave $(\overline{u}_{-},a_{-},0)$ of the Whitham shock (55). To investigate this further, we compare this portion of the heteroclinic TW solutions with solitary wave solutions that move with the shock velocity $V$ on the background $\overline{u}_{-}$ (the solitary wave amplitude $a_{-}$ is obtained from $c(\overline{u}_{-},a_{-},0)=V$). Figure 6 consists of the heteroclinic TW solutions depicted in Fig. 5 overlaid with the left solitary wave (dashed red) and the right periodic wave $\varphi(\theta;0,1,k_{+})$ (dashed blue) that form the corresponding admissible Whitham shock. Both the left solitary waves and the right periodic waves are visually indistinguishable from the heteroclinic TW in their respective regions of validity. This corroborates our formulation and interpretation of the zero dispersion limit of heteroclinic TW solutions as admissible Whitham shocks where the left wave is a solitary wave that rapidly transitions to a co-propagating finite amplitude periodic wave.

Admissible Whitham shock solutions and characteristics

[TABLE]

corresponding to the loci $\blacktriangle$ and $\blacksquare$ in Table 2 are shown in Figures 7 and 8, respectively. Both Whitham shocks with real characteristic velocities are weakly compressive in the first characteristic family $\Gamma_{1}$ because $\lambda_{1}^{(+)}<V=\lambda_{1}^{(-)}$, while the second characteristic family is weakly expansive $\lambda_{2}^{(-)}=V<\lambda_{2}^{(+)}$, and the third characteristic family passes through the Whitham shock, decelerating $V<\lambda_{3}^{(+)}<\lambda_{3}^{(-)}$ for $\blacktriangle$ and accelerating $V<\lambda_{3}^{(-)}<\lambda_{3}^{(+)}$ for $\blacksquare$ [21]. Consequently, we refer to this class of Whitham shock solutions as weak 1-shocks.

The degeneration of the periodic wave to a solitary wave on the left state in Fig. 5 allows us to obtain three additional admissible Whitham shock loci by reflecting the initial data (55) and the heteroclinic traveling wave across $x=0$. These solutions are weakly compressive in the second characteristic family, $\Gamma_{2}$, hence are called weak 2-shocks. Further implications of this observation are discussed in Sec. 5.3.

5.2 Periodic to periodic

Let us continue our computation of the jump conditions Eqs. (45)–(47) where we now consider shock solutions of the modulation equations corresponding to one periodic wave connected to another. We slightly modify the Riemann problem (55) from the previous section where the left state is assumed to have zero wavenumber. In this section, we scale the Whitham shock (44) so that the left state consists of an arbitrary periodic wave, and the right state consists of a periodic wave with zero mean, unit amplitude and arbitrary wavenumber. This parameter set results in the shock solution

[TABLE]

where $V=c_{\pm}$. The Rankine Hugoniot jump relations (45)–(47) are a nonlinear system of three equations that relate the four remaining parameters: $k_{+}$, $\overline{u}_{-}$, $a_{-}$, and $k_{-}$. We use $\overline{u}_{-}$ as the continuation parameter to obtain the one-parameter family of Whitham shock loci

[TABLE]

Attempts to solve the jump conditions (45)–(47) approximately by using the Stokes wave approximation (19) yield no nontrivial asymptotically valid solutions, aside from the solitary wave to periodic shocks discussed previously. As a result, we rely on numerical continuation along the parameterized curves $(k_{+},a_{-},k_{-})(\overline{u}_{-})$ where we start from the three Whitham shock loci calculated in Section 5.1 so that $k_{-}=0$ and $\overline{u}_{-}$, $k_{+}$, $V$ are initialized from Table 2 ($a_{-}$ satisfies $c(\overline{u}_{-},a_{-},0)=V$). We numerically continue solutions to the jump conditions (45)–(47) with Matlab’s fsolve function for decreasing values of $\overline{u}_{-}$. Continuation is terminated when the value $k_{+}(\overline{u}_{-})<10^{-3}$ is reached, indicating that the oscillatory wavetrain in the left far-field is nearly a solitary wave. Figures 9(a), (c) and (e) present the three Whitham shock loci. We also compute heteroclinic TW solutions for select $\overline{u}_{-}$, seeded with left/right states from Whitham shock loci, and plot them in Figures 9(b), (d), and (f). Consequently, we find that all three periodic to periodic Whitham shock loci are admissible.

Dynamic stability of admissible Whitham shocks is determined by the hyperbolicity of the Whitham modulation equations (recall Fig. 4). Whitham shocks with real characteristic velocities are identified by the white background in Fig. 9 (a), (c), and (e) while the gray background denotes complex characteristic velocities $\mathrm{Im}\lambda_{1}\neq 0$ and $\mathrm{Im}\lambda_{2}\neq 0$ and therefore a modulationally unstable heteroclinic TW. This will be investigated further in Sec. 5.3.

In Figure 10, we plot the numerically computed characteristic curves $\Gamma_{j}$, $j=1,2,3$ and the shock trajectory for the periodic to periodic Whitham shock in Fig. 9(a) where $\overline{u}_{-}=0.6$ (the corresponding heteroclinic TW is shown in Fig. 9(b)). We observe that all characteristic families pass through the shock and satisfy

[TABLE]

Shocks that are not compressive in any characteristic family violate the Lax entropy conditions. In the conservation law community, non-Lax shocks can be identified as admissible when they are the limit of vanishing dissipative-dispersive heteroclinic TWs [24, 20]. In such cases, they are referred to as undercompressive. Here, we find that the vanishing dispersion limit of heteroclinic periodic to periodic TWs generally converge to undercompressive Whitham shocks. To prove this, consider Fig. 11 where the scaled characteristic velocities $(\lambda_{j}-\overline{u})/a$, $j=1,2$ and shock velocity $(V-\overline{u})/a$ are plotted in the scaled coordinate $\tilde{k}=a^{-1/4}k$. Since $\lambda_{3}>\lambda_{2}$, we see that all admissible Whitham shocks with real characteristic velocities exhibit the undercompressive relations (63) except for the very narrow band of waves $0.662\lesssim ka^{-1/4}\lesssim 0.672$. Because this band is so narrow in parameter space, we do not investigate these solutions any further. Undercompressive shocks were first described for general $2\times 2$ systems in [59] and have since been observed, for example, in fluid dynamics [60].

5.3 Two shock solutions

In Sections 5.1 and 5.2, we established the admissibility of KdV5-Whitham shocks by computing heteroclinic traveling wave solutions to the KdV5 equation (2). Each locus of admissible Whitham shocks depicted in Fig. 9(a), (c), (e) include, up to scaling and Galilean symmetries, pairs of shocks with reflected $+$ and $-$ states. The implication of this observation is that it is possible to arrange for co-propagating Whitham shocks with the same equilibrium on the left and right, respectively, as well as the periodic wave on the right and left, respectively. Can we compute homoclinic TWs consisting of two co-propagating Whitham shocks? Furthermore, can this construction be extended to solutions tending to a periodic wave in the far-field? We affirmatively answer these questions now.

We formulate the double Whitham shock problem as a locally periodic wave that terminates after a finite number of oscillation periods and transitions to a different periodic wave, possibly a solitary wave, in the far field. To this end, let us consider solutions of the form

[TABLE]

where the subscript i refers to the inner periodic wave with $n$ oscillations and subscript o denotes the outer periodic wave in the far field. Following the normalization in Sec. 5.2, we have normalized the inner periodic solution to have zero mean and unit amplitude.

Figure 12 depicts the structure of an example homoclinic TW solution in the physical plane overlaid with co-propagating weak 1-shock and weak 2-shock solutions that satisfy the jump conditions for zero wavenumber in the far field, i.e., $k_{\rm o}=0$. In general, the far-field periodic wave can be chosen from a continuous set of nonzero wavenumbers (cf. Fig. 9 for the range of possible values, equating $k_{\rm o}=k_{-}$).

Utilizing the three Whitham shock loci in Table 2, the scaling symmetries (16) and (17), and double Whitham shock data (64), we compute homoclinic TW profiles and initialize the KdV5 equation with these solutions plus small, band-limited noise (see Appendix for computational details). The time evolution of the perturbed homoclinic TW solutions then allows us to corroborate our hypothesis that traveling waves comprised of modulationally unstable portions will be unstable. Indeed, in Fig. 13, we depict the evolution of three distinct homoclinic solutions with $n=25$ oscillation periods. Note that our computational approach enables for the construction of homoclinic solutions with any number of oscillation periods. The intermediate periodic wave in Fig. 13(b) lies on the Whitham shock locus $\CIRCLE$ in Table 2, which exhibits complex characteristic velocities. As expected, the homoclinic TW solution is unstable. In Figs. 13 (a) and (c), the homoclinic oscillatory defect is comprised of periodic waves on the Whitham shock loci $\blacktriangle$ and $\blacksquare$ in Table 2, respectively, where the modulation equations are strictly hyperbolic. Perturbed homoclinic traveling waves are numerically stable over a long integration time.

These homoclinic TWs are similar to solutions obtained in the investigation of reversible Hamiltonian systems [61, 46], though we have presented an alternative approach by which multi-pulse solutions with an arbitrary number of peaks can be computed from the Whitham shock loci in Table 2.

Similarly, homoclinic TWs consisting of distinct inner and outer periodic waves can be constructed from Whitham shock loci shown in Fig. 9. Numerical computations of initially perturbed homoclinic traveling waves of this type along with their long time evolution are displayed in Fig. 14, again revealing numerically stable evolution for those solutions with all real characteristic velocities—Fig. 14(a,d)—and numerically unstable evolution for those double Whitham shocks that exhibit complex characteristic velocities—Fig. 14(b,c).

6 Applications and extensions

In this section, we discuss applications of the theory presented in this manuscript and various physical systems to which it applies. We construct the modulation solution for the TDSW depicted in Fig. 1 for the KdV5 equation (2) and discuss how the same construction can be carried out for the Kawahara equation (4), where these structures were numerically coputed in [26]. This section concludes with numerical experiments simulating the Riemann problem (7) for the Whitham equation [23, 28], which describes weakly nonlinear, fully dispersive gravity-capillary water waves.

6.1 Dam break problem for the modulation equations: Traveling dispersive shock waves

Several recent works have considered step initial data for systems with higher order dispersion [26, 54, 53, 40, 41]. Among each of these studies, one key observation from numerical simulations was made: the nonlinear structure arising from smoothed step initial data developed into a wave with two distinct regions, depicted in Fig. 1. We now identify the two distinct regions as follows. At the left, trailing edge of the TDSW, there is a heteroclinic traveling wave corresponding to an admissible Whitham shock solution of the modulation equations (cf. Sec. 5.1). The heteroclinic traveling wave terminates at a partial DSW described by a continuous, simple wave solution of the modulation equations (cf. Sec. 3.1.2) at the leading edge. Because the modulated leading edge initiates from an oscillatory wavetrain, rather than a solitary wave, it is termed a partial DSW. Similar partial DSWs have been observed in the initial boundary value problem of nonlinear systems with lower order dispersion [62, 63]. We now explain this Whitham shock-rarefaction modulation solution to the intial value problem

[TABLE]

Numerical simulations of (65) in Fig. 15 show two distinct regions in the solution. The boundaries of these regions are defined by identifying the following velocities: the trailing edge Whitham shock velocity, $V_{\rm s}$, the intermediate velocity $V_{i}$ joining the TW to the partial DSW and $V_{+}$ the leading, harmonic edge velocity where the partial DSW terminates when the modulation amplitude vanishes $a\to 0$. If we denote the parameters of the periodic portion of the heteroclinic wave by $(\overline{u}_{r},a_{r},k_{r})$ and the leading harmonic wave edge wavenumber by $k_{+}$, then $V_{i}=\lambda_{2}(\overline{u}_{r},a_{r},k_{r})$ and $V_{+}=\lambda_{2}(0,0,k_{+})$. The TDSW traveling wave portion bears a striking resemblance to the traveling wave that is identified by the symbol $\blacksquare$ in Fig. 5. The Whitham shock locus for this solution in Table 2 is well approximated using a Stokes wave in Table 1, scaled by (16), (17) so that $\overline{u}_{-}\to 1$ and $\overline{u}_{+}$ is a to-be-determined parameter. The intermediate periodic wave properties are then given as a one parameter family

[TABLE]

where $\zeta\approx 0.9724$ is the positive square root of the largest, positive real root of the quartic polynomial (60), and $\overline{u}_{\rm r}$ is the mean to be determined by matching to a simple wave describing the transition from the oscillatory wavetrain to the leading, constant level, $\overline{u}=0$ propagating with velocity $V_{+}$. The computation of the self-similar simple wave curve is described in Sec. 3.1.2, denoted by $\mathbf{q}(x/t)$. The Whitham shock locus, 2-wave curve and the far-field mean values completely determine the Whitham shock-rarefaction modulation solution

[TABLE]

This shock-rarefaction solution can be viewed as the higher dispersion analog of the shock-rarefaction solution for the shallow water dam break problem or the gas dynamics piston problem.

In Figure 16, we compare the Whitham shock-rarefaction solution of the modulation equations with numerical simulations of smoothed step initial data. Numerically extracted values of the modulation parameters exhibit excellent agreement with the Whitham shock-rarefaction modulation solution. The small oscillations in the modulation variables that occur along the intermediate equilibrium state are due to higher order dispersive effects that are not captured by leading order modulation theory.

6.2 Traveling waves in gravity-capillary water waves

In this section, we consider two models of gravity-capillary water waves that incorporate higher order dispersive effects. We recall the linear dispersion relation for right-moving water waves on a constant depth, normalized to $1$, with surface tension is given in nodimensional form by

[TABLE]

In this section, we are motivated by the physical problem in which $B$ is sufficiently close to $1/3$, so that $B-1/3\sim k^{2}$ when $k\ll 1$ so that third and fifth order dispersion are in competition.

6.2.1 Kawahara equation

A weakly nonlinear, long wave approximation for water waves with Bond number $B$ sufficiently close to but less than $1/3$ is the Kawahara equation (4), where $\alpha=1$[32]. Previous numerical simulations of the Riemann problem (7) for eq. (4) show that for sufficiently large values of the initial jump $\Delta$, the result is a TDSW [26]. Motivated by these numerical experiments, we now construct traveling wave solutions of the underlying PDE by computing the Whitham shock loci of the modulation system using the weakly nonlinear Stokes approximation.

The Kawahara Whitham modulation equations are obtained by averaging the first two conservation laws

[TABLE]

where $\varphi=\varphi(\theta)$ is a periodic, traveling wave solution to (4). Based on numerical computations in [26], we expect that a heteroclinic TW solution exists and connects two disparate far-field waves in which the leftmost state is a solitary wave with mean $\overline{u}_{-}$. The jump conditions are

[TABLE]

We use the leading, periodic solution in the form of a Stokes wave [26]

[TABLE]

with approximate phase velocity

[TABLE]

Inserting the Stokes expansion (69) into the jump conditions (66)–(67), we then solve the algebraic system numerically for a fixed value of $\overline{u}_{-}=\Delta$ and the ensuing wave parameters are then compared to numerical solutions of the Riemann problem (7). In Figure 17, we compare the periodic wave parameter values from the jump conditions using the Stokes approximation (69) and those obtained from numerical computations of heteroclinic traveling waves using Matlab’s bvp5c. Details of the numerical method

Since we are using the Stokes approximation for the leading periodic wave, we only expect agreement with numerical simulations for sufficiently small wave amplitude. However, the results depicted in Fig. 17 demonstrate that even for quite large values of $a_{+}$, the Stokes approximation remains very accurate. Notice further that the numerical computations of traveling waves, whose parameters are depicted by the red dots, cease for $\Delta\lesssim 0.58$. Here, numerical fixed point computations of traveling waves do not converge, yet the Whitham shock locus continues well past this point. We interpret these computations as evidence of inadmissible Whitham shocks for the Kawahara equation. The mere existence of a Whitham shock is not a sufficient condition to guarantee the existence of a traveling wave solution.

6.2.2 Gravity-capillary Whitham equation

In 1967, Whitham proposed a method to capture the full dispersive properties of a physical system by using a convolution kernel so that the linear dispersion relation of a weakly nonlinear model is prescribed [28]. The so-called Whitham equation takes the form

[TABLE]

where $\mathcal{K}=\mathcal{F}^{-1}c(k)$, where $c(k)$ is the phase velocity. For applications to right-moving gravity-capillary water waves, we take $c(k)=\sqrt{\frac{(1+Bk^{2})\tanh k}{k}}$.

Solutions similar to those appearing in the KdV5 equation (2) and the Kawahara equation (4) are expected to also persist for the Whitham equation when the Bond number $B$, is sufficiently close to $1/3$ and third and fifth order dispersion are in balance. We revisit our original motivation for investigating traveling wave solutions to the KdV5 equation (2) by conducting direct numerical simulations of the Riemann problem (7) for (71) with $\Delta=0.25$. The results of the simulations at $t=2000$ are shown in Figure 18 for different values of the Bond number, $B$. The simulations indicate that for $B$ near $1/3$ (e.g., $0.3\leq B\leq 0.375$), where third order dispersion is weak, familiar TDSW-looking dynamics emerge (cf. Fig. 1). The trailing edge resembles a nearly uniform traveling wave and the leading portion resembles a partial DSW terminating into the intermediate nonlinear wavetrain. These numerical simulations suggest that the Whitham shocks described in this manuscript are, in fact, quite general. With that being said, the numerical simulations depicted in Figure 18 for $B=0.4$ and $B=0.425$ no longer resemble TDSWs, instead they look like classical, convex DSWs. In this regime, lower order dispersion dominates, preventing the higher order dispersive balance that appears to be needed for the existence of admissible Whitham shocks. A future avenue of research will be to understand the transition from solutions discussed throughout this manuscript to the classical convex, KdV-type DSW solutions arising from step initial data.

7 Conclusions

In this manuscript, we considered the fifth order KdV equation (2) and its novel traveling wave solutions. The heteroclinic TW solutions studied here consist of two periodic waves that co-propagate and transition between one another on a length scale commensurate with the length scale of a single wave period. Utilizing the framework of Whitham modulation theory that describes the zero dispersion limit of nonlinear oscillatory wavetrains via the Whitham system of conservation laws, we establish that these heteroclinic traveling wave solutions limit to discontinuous shock solutions of the Whitham modulation equations, i.e., Whitham shocks. Within this framework, we prove that a heteroclinic traveling wave solution of the KdV5 equation (2) connecting disparate periodic waves necessarily satisfies the Rankine-Hugoniot jump conditions of the Whitham equations in conservative form. The jump conditions reveal admissible Whitham shock loci—far-field periodic wave parameters that are then used as candidate traveling wave solutions of the governing PDE.

Numerical computations of the characteristic velocities in tandem with the novel traveling wave profiles allow us to specify the characteristic structure of Whitham shocks. These computations reveal that, but for a negligibly small portion of parameter space, admissible Whitham shocks are undercompressive, i.e., each characteristic family passes through the shock and therefore violates the Lax entropy condition. The characteristic velocities of the modulation equations also reveal the stability of heteroclinic traveling wave solutions. Since the hyperbolicity of the Whitham modulation equations is a necessary condition for the stability of periodic solutions [58], the stability of traveling wave solutions considered in this manuscript are determined by the stability of the periodic waves that constitute them, which can be readily checked by examining the Whitham shock velocities.

A prominent contributing feature to the emergence of novel heteroclinic traveling waves and admissible Whitham shocks found in this manuscript is the presence of higher order dispersive effects. For example, in scalar problems exhibiting nonlocal dispersion, such as the Whitham equation in Sec. 6, we demonstrate numerical evidence of the ubiquity of Whitham shocks in other systems. To further analytically study these heteroclinic traveling waves connecting two far-field periodic waves, one promising route is to investigate integrable higher order nonlinear, dispersive systems [64] so that heteroclinic TWs may be understood in the context of the inverse scattering transform. For instance, the mathematical structure underlying these traveling waves could be eludicated though a detailed study of the Lax equation [65, 66]–an integrable equation with higher order dispersive terms.

Of related importance is the extensive literature describing spontaneous, localized pattern formation in dissipative systems, where the prototypical example is the Swift-Hohenberg equation that leads to a fourth order ODE for stationary solutions that resemble the TW ODE considered here [67, 68, 69, 70, 71, 72]. Here, stationary solutions of the governing equation were computed in which an equilibrium state spontaneously transitions to a localized periodic state and back to equilibrium, the stationary, dissipative analog of the homoclinic TWs limiting to double Whitham shocks computed here. The results have since been extended to consider steady solutions in which a large, localized periodic pattern persists on a small amplitude oscillatory background [73]. A promising avenue for further study is the consideration of modulation dynamics in the Swift-Hohenburg equation. Of course, the ensuing dynamics and stability of these solutions will differ from the dynamics studied here because the regularizing mechanisms in each case—dissipation or dispersion respectively—are wholly different. From this point of view, TW solutions occurring in higher order dispersive systems could serve as a bridge between conservative Hamiltonian systems and dissipative pattern forming systems.

This work was supported by NSF grant DMS-1816934. The authors gratefully acknowledge valuable discussions with Michael Shearer and Noel Smyth.

Appendix A

In this appendix, we outline numerical methods used throughout this manuscript.

First, we describe the method used to compute periodic traveling wave solutions to (2). We seek a one parameter family of solutions, $\tilde{\varphi}$ parameterized by an arbitrary wavenumber $\tilde{k}$ and fixed amplitude $\tilde{a}=1$ and mean $\tilde{\overline{u}}=0$. For clarity in notation we will continue to denote functions, variables, and parameters associated with the one parameter family with tildes “$\sim$”. The one parameter family of solutions is computed through a pseudospectral Fourier projection of the differential equation (8) with $A=0$ without loss of generality and where $c=\tilde{c}$ is an eigenvalue of the differential operator and used as a continuation variable.

On the interval $\xi\in[-\pi,\pi)$ with periodic boundary conditions the periodic wave is approximated by

[TABLE]

which yields a nonlinear system of algebraic equations for the Fourier coefficients $a_{n}$. The nonlinear term is numerically computed in physical space while the linear terms are computed directly via the fast Fourier transform. The resulting nonlinear system is then solved via Matlab’s fsolve function. The number of Fourier modes, $N$, is chosen sufficiently large so that Fourier modes decay to $O(10^{-16})$. The nonlinear solution can then be rescaled using the symmetries (16), (17) so that the solution is of unit amplitude and zero mean. The solution library following scaling to unit amplitude and zero mean consists of approximately $6000$ periodic waves non-uniformly spaced on the grid $\tilde{k}\in[0.002,1.25]$.

We compute the KdV5-Whtiham modulation equations (40) by numerically averaging the periodic solutions of the KdV5 equation (2) over their periods. The averaged integrals, $\tilde{I}_{j}$ and $\tilde{J}_{2}$ in (41) are computed using the spectrally accurate trapezoidal rule. We use a finite difference scheme on the unequally spaced grid in $\tilde{k}$ to compute the appropriate gradients of $\tilde{I}_{j}$ and $\tilde{J}_{2}$ which are then used to build the matrix $\mathcal{A}$ appearing in Eq. (40). The eigenvalues and eigenvectors of $\mathcal{A}$ are determined directly using Matlab’s eig function. Next, we compute the quantities $\mu_{j}$ in Eq. (42) on the grid $\tilde{k}$ using the aforementioned finite differencing scheme.

The numerical computation of heteroclinic traveling waves is implemented using the iterative Newton-conjugate gradient method [74, 75]. The numerical implementation relies upon periodic boundary conditions and projection onto a Fourier basis. By considering the even reflection of the Whitham shock at the peak of the right oscillatory wave $\varphi_{+}$, the resulting extended structure then satisfies the periodic boundary condition requirement for the numerical method. An example of the initial guess for the iterative solver is shown in Fig. 19, where intermediate periodic wave wavenumber is determined by the jump conditions. Outside the oscillatory region in the initial guess, the solution abruptly transitions to the known far-field constant value.

Homoclinic traveling waves corresponding to double Whitham shock solutions of the modulation equations are computed using the fifth order collocation method bvp5c in Matlab with periodic boundary conditions imposed on the ODE (8). The solutions are used as initial contions to KdV5 (2), which is numerically evolved with a pseudospectral spatial discretization via the FFT and an integrating factor fourth order Runge-Kutta method [76] with time step in the range $\Delta t\in[10^{-4},10^{-3}]$.