Dispersive Riemann problem for the Benjamin-Bona-Mahony equation

Long time dynamics of the smoothed step initial value problem or dispersive Riemann problem for the Benjamin-Bona-Mahony (BBM) equation $u_t + uu_x = u_{xxt}$ are studied using asymptotic methods and numerical simulations. The catalog of solutions of the dispersive Riemann problem for the BBM equation is much richer than for the related, integrable, Korteweg-de Vries equation $u_t + uu_x + u_{xxx} =0.$ The transition width of the initial smoothed step is found to significantly impact the dynamics. Narrow width gives rise to rarefaction and dispersive shock wave (DSW) solutions that are accompanied by the generation of two-phase linear wavetrains, solitary wave shedding, and expansion shocks. Both narrow and broad initial widths give rise to two-phase nonlinear wavetrains or DSW implosion and a new kind of dispersive Lax shock for symmetric data. The dispersive Lax shock is described by an approximate self-similar solution of the BBM equation whose limit as $t \to \infty$ is a stationary, discontinuous weak solution. By introducing a slight asymmetry in the data for the dispersive Lax shock, the generation of an incoherent solitary wavetrain is observed. Further asymmetry leads to the DSW implosion regime that is effectively described by a pair of coupled nonlinear Schr\"{o}dinger equations. The complex interplay between nonlocality, nonlinearity and dispersion in the BBM equation underlies the rich variety of nonclassical dispersive hydrodynamic solutions to the dispersive Riemann problem.


Introduction
The unidirectional regularized shallow water equation, often called the Benjamin-Bona-Mahony (BBM) equation [4], was originally introduced by Peregrine [41] to numerically model the development of shallowwater undular bores. The canonical representation of the BBM equation as an asymptotic model for shallow water waves has the form u t + u x + ε(uu x − αu xxt ) = 0, (1.1) where ε 1 and α = O(1) are constants. Equation (1.1) shows that, to leading order in ε, u t = −u x so that (1.1) is asymptotically equivalent to the Korteweg-de Vries (KdV) equation Despite the asymptotic equivalence of the KdV and BBM equations, their mathematical properties are drastically different. In particular, the BBM equation (1.1) yields more satisfactory short-wave behavior, owing to its regularization of the unbounded growth in the frequency, phase and group velocity values present in the linear dispersion relation for the KdV equation. This leads to a number of advantages of the BBM equation in the context of well-posedness and computational convenience. On the other hand, the BBM equation lacks integrability, which is a prominent feature of the KdV equation with many remarkable consequences including an infinite number of conserved quantities and the existence of solitons, localized solutions exhibiting elastic interactions. In addition to the particular application of the KdV and the BBM equations as asymptotic models for weakly nonlinear shallow water waves and other dispersive media, these equations exemplify two qualitatively different ways to regularize a scalar nonlinear conservation law, acutely displayed by their normalized versions: and dropping the tildes.
In this paper, we study the dispersive regularization of an initial step-like transition for the BBM equation (1.4). This kind of dispersive Riemann problem was introduced and studied for the KdV equation (1.3) by Gurevich and Pitaevskii [25] using Whitham modulation theory [51]. The most prominent feature of the KdV dispersive Riemann problem is the occurrence of a dispersive shock wave (DSW), a smooth, expanding nonlinear wavetrain that replaces the discontinuous, traveling shock solution of the hyperbolic conservation law u t + ( 1 2 u 2 ) x = 0-the Hopf equation-subject to the initial data when u − > u + (the Lax entropy condition) and the shock speed c satisfies the Rankine-Hugoniot jump condition c = 1 2 (u − +u + ). If u − < u + , the initial discontinuity (1.6) for the Hopf equation evolves into an expanding, continuous rarefaction wave (RW), which retains its general structure under KdV dispersive regularization subject to smoothing of its weak discontinuities and small-amplitude oscillations at the left corner of the classical RW. DSWs and RWs are the only possible KdV dispersive regularizations of the Riemann problem. Because qualitatively similar DSWs and RWs have been found to arise in a variety of "KdV-like" dispersive regularizations of hyperbolic conservation laws, such DSWs and RWs are referred to as classical or convex.
It turns out that the BBM dispersive Riemann problem exhibits a number of features that are markedly different from the KdV dispersive Riemann problem, i.e., are nonclassical. To elucidate them, we first note that the consideration of discontinuous Riemann data (1.6) for a dispersive equation leads to anomalous features, such as the generation of waves with unbounded phase and group velocities in the KdV equation [6]. One way to avoid this unphysical behavior, is to introduce smoothed Riemann data, e.g., where the parameter ξ represents the characteristic width of the initial transition. If ξ → 0, (1.7) converges to the discontinuous data (1.6). For applications, the consideration of a nonlinear dispersive equation with smoothed step data of the type (1.7) constitutes a physically meaningful dispersive Riemann problem. Henceforth, we will refer to the data (1.7) as a smoothed step. Drawing upon the theory of Riemann problems for hyperbolic conservation laws, one may hypothesize that the leading order, long-time asymptotic solution of the dispersive Riemann problem is independent of the value of ξ. Indeed, the evolution of a smooth, compressive step with u − > u + in the KdV equation is asymptotically described by a self-similar solution of the Whitham modulation equations, which does not involve the initial step's width or shape and only depends on the boundary parameters u ± where x → ±∞ [1]. However, solutions of BBM Riemann problems turn out to be quite sensitive to the value of ξ, leading to qualitative differences in the wave patterns generated by the evolution of the smoothed step (1.7) with the same values of u ± but different values of ξ. This can be understood by noting that the BBM equation (1.4) can be put in the integro-differential form [4] u t (x, t) = 1 4 +∞ −∞ sgn(x − y)e −|x−y| u(y, t) 2 dy, (1.8) explicitly reflecting its nonevolutionary and nonlocal character. The factor e −|x−y| in (1.8) introduces the intrinsic BBM nonlocality length scale = 1. Consequently, the regimes of the initial value problem (1.4), (1.7) characterized by ξ 1 and ξ 1 are expected to be qualitatively different. The nonevolutionary character of the BBM equation also results in nonconvexity of the linear dispersion relation, with zero dispersion points for certain values of the wavenumber and background. Generally, we find that the complex interplay between nonlocality, nonlinearity and dispersion in the BBM equation gives rise to remarkably rich nonclassical dynamics that are revealed in the study of the Riemann problem.
It is worth commenting briefly on the relevance of the BBM equation (1.1) and its reduced version (1.4) for physical applications. While (1.1) is a generic asymptotic model of weakly nonlinear, dispersive waves, it is not an asymptotically resolved equation. The small parameter cannot be scaled out of the equation while maintaining its asymptotic validity. For example, the transformation (1.5) breaks the asymptotic long wave assumption. However, there are geophysical scenarios in which the reduced BBM equation (1.4) is a reasonable physical model. An example scenario is the circular, free interface between two viscous, Stokes fluids whose long wave evolution is well-approximated by the conduit equation u t + uu x − uu txx + u t u xx = 0 [39,34]. Other examples include a class of models for magma transport in the upper mantle u t + (u n ) x − (u n (u −m u t ) x ) x = 0 where n ∈ [2,5], m ∈ [0, 1] are parameters [43,48] and channelized water flow beneath glaciers u t + (u α (1 − (sgn(u t )|u t u −1 | 1/n ) x ) β ) x = 0 where n ≥ 1, α > 1, β > 0 are parameters [47]. The BBM equation (1.4) captures the quadratic hydrodynamic flux 1 2 u 2 x and nonlocal, nonevolutionary dispersion −u xxt inherent to these physical models in certain parameter regimes. Therefore, the BBM equation in its reduced, normalized form (1.4) represents a significant, physically inspired dispersive regularization of the Hopf equation that is distinct from the KdV dispersive regularization (1.3).
In this paper, we report on wave structures and phenomena associated with the Riemann problem for the BBM equation. Solutions of (1.4), (1.7) include familiar wave patterns such as classical RWs and DSWs that also appear in the dispersive Riemann problem for the KdV equation. However, the BBM dispersive Riemann problem additionally exhibits a variety of nonclassical waves that do not appear as solutions to KdV-type equations. These include dynamics that depend strongly upon the smooth step transition width ξ and those that do not. When ξ 1, two-phase linear wavepackets accompany the RWs and DSWs. When u + = −u − , the BBM equation (1.4) admits a stationary, discontinuous weak solution that is stable according to the Lax entropy condition [31,32] if 0 < u − . For smooth initial data (1.7) in which u + = −u − < 0, the numerical solution evolves toward the Lax shock solution but is accompanied by increasingly short waves. These waves are described by an approximate self-similar solution to the BBM equation that we term a dispersive Lax shock. However, the solution is observed to be structurally unstable to slightly asymmetric initial conditions (u + = −u − ), numerically evolving into an incoherent collection of small amplitude, short waves and solitary waves. Remarkably, for u + = −u − > 0, for which the discontinuous solution violates the Lax entropy condition and is an expansive shock wave solution of the Hopf equation, the solution of the corresponding initial value problem (1.4), (1.7) with 0 < ξ 1 exhibits a smooth solution that approximates the discontinuity. In this case, the smooth expansion shock decays algebraically in time, giving way to a RW as shown by an asymptotic analysis in [15]. In that paper, it was also observed that when the initial data fail to be symmetric, specifically (u + = −u − ,), then the expansion shock may be accompanied by a train of one or more solitary waves. Finally, for arbitrary ξ > 0, we identify a regime of the BBM Riemann problem that gives rise to DSW implosion in which a nonlinear two-phase interaction develops from the small amplitude edge of the DSW. DSW implosion was previously predicted and observed in the conduit and magma equations [35].
As we have noted, the BBM equation can be viewed as a prototypical model for wave phenomena that arises in other physical systems described by nonevolutionary, nonlinear dispersive wave equations. While DSW implosion has been observed in the conduit and magma equations, other intriguing BBM wave patterns such as expansion shocks, solitary wave shedding, and dispersive Lax shocks await their realization in physical systems. These nonclassical wave patterns require new mathematical approaches for their interpretation. We utilize detailed numerical simulation, asymptotic methods, and properties of the BBM equation in order to provide a relatively complete classification of solutions of the BBM dispersive Riemann problem.
The paper is organized as follows. Owing to the classification's richness and complexity, we begin with a high-level description of our findings in Section 2. We review the long-time asymptotic description of linear waves in Section 3 as they play a prominent role in BBM dispersive Riemann problem dynamics. This is followed by Section 4 in which rarefaction waves, dispersive Lax shocks, and expansion shocks are studied. This leads naturally to the description and analysis of the shedding of solitary waves in Section 5. Section 6 is devoted to DSW solutions, DSW implosion that includes a new two-phase description of the dynamics, and the generation of incoherent solitary wavetrains. We conclude in Section 7 with a discussion and future outlook.

The BBM Riemann problem: summary of solutions and analysis
With the scaling of variables, u → αv, x → βX, t → γT, Thus (1.4) is invariant under the change of variables (2.1) if αγ = β = ±1. In particular the BBM equation remains unchanged after (u, x) → (−u, −x). We therefore restrict the study of the dispersive Riemann problem (1.4), (1.7) to initial data satisfying u − + u + ≥ 0. In this section and throughout the paper, we refer to Figures 1 and 2, which concisely capture the BBM dispersive Riemann problem classification. The figures include a partitioning of the u + ≥ −u − half plane for the smoothed step initial data (1.7) and corresponding representative numerical simulations of each qualitatively distinct wave pattern that emerges during the course of BBM (1.4) evolution. Figure  1 corresponds to the broad, slowly varying smoothed initial step in which the transition width parameter ξ 1 (ξ = 10 in all numerical simulations). Figure 2 represents solutions for which the smoothing is narrow and sharp, corresponding to ξ 1 (ξ = 0.1 in all numerical simulations). Some wave patterns persist for both large and small transition width, but many do not. We do not investigate the subtleties of the crossover ξ ∼ 1 in which the transition width coincides with the nonlocality length, choosing instead to utilize scale separation for ξ 1 and ξ 1 where asymptotic methods are available.
Large transition width ξ 1, Figure 1 When u + > max(µu − , −u − ) where µ = e −3/2 /4, convex RWs (region a) and DSWs (region b) are generated, qualitatively similar to those generated by KdV (1.3) evolution. In Sec. 6.2, we use the DSW fitting approach to make quantitative predictions for the DSW's properties (edge velocities, harmonic edge wavenumber, soliton edge amplitude, and DSW structure near the harmonic edge). We also identify the loss of convexity at the harmonic edge (where a zero linear dispersion point is attained) as the progenitor for nonclassical wave patterns when u + is below the line µu − .  For initial data sufficiently close to the line u + = µu − (region c), DSW implosion occurs. A weakly nonlinear, two-phase modulation theory using coupled nonlinear Schrödinger equations is developed in Sec. 6.3, the two wave envelopes associated with modulations of short and long waves, respectively. We find that there is a spatial redistribution of energy from the long waves in a partial DSW on the right to a short wave wavepacket on the left. The partial DSW is typically accompanied by some number of depression envelope solitary waves. Very similar DSW implosion dynamics are also observed in the small transition width regime ξ 1 (see region f in Fig. 2). There exists a stationary, discontinuous compressive shock solution of the BBM equation when u + = −u − < 0 (region e). In Sec. 4.1, we show that smoothed step initial data for this case exhibits the usual nonlinear self-steepening but is also accompanied by the continued production of shorter and shorter waves during evolution. In a neighborhood of x = 0, an oscillatory overshoot appears that quickly saturates to a constant magnitude. An approximate self-similar solution in the form u(x, t) = g(xt) is obtained that describes these oscillations as t → ∞. This feature is reminiscent of Gibbs' phenomenon in the theory of Fourier series and linear dispersive partial differential equations with discontinuous initial data [7]. We refer to this wave pattern as a dispersive Lax shock. However, if the initial data is perturbed asymmetrically so that −u − < u + < σ 3 (ξ)u − (region d), then the evolution exhibits a large number of solitary waves and highly oscillatory wavepackets that do not exhibit a discernible, coherent pattern. We refer to this as the incoherent solitary wave regime and describe it in Subsection 6.4. The dispersive Lax shock and the incoherent solitary wave patterns also occur for small transition width ξ 1 (see regions h and g in Fig. 2).
Small transition width ξ 1, Figure 2 Generally, the features observed for large transition width in Fig. 1 persist into the small transition width regime shown in Fig. 2 but with modification. Rarefaction waves are accompanied by a decaying expansion shock in regions a and b, depression solitary waves propagating to the left in regions b and c, and small amplitude linear waves propagating to the left in regions c and d.
From the ideal expansion shock solution for symmetric data (u − = −u + < 0) (reviewed in Sec. 4.3) into region a for asymmetric data, a short-time analysis of the solution reveals an asymmetry that favors the generation of a depression wave for x < 0 (Sec. 5.2). When the threshold u + = −σ 1 (ξ)u − is exceeded, the asymmetry produces a solitary wave that is shed from the expansion shock. Greater asymmetry generates the shedding of more depression solitary waves until u + ≥ 0 and an expansion shock is no longer generated. Nevertheless, solitary wave shedding persists until u + = σ 2 (ξ)u − . Solitary Wave shedding thresholds σ 1 (ξ) and σ 2 (ξ) are analyzed in Sec. 5.3.
When u − ≥ 0 (regions c-e), the small transition width gives rise to the generation of small amplitude wavepackets that accompany the RW or DSW. In the absence of solitary wave shedding, we can conveniently partition the initial smoothed step in the Fourier domain to prescribe initial data for the linearized BBM equation. In Sec. 3, long-time asymptotic analysis of the solution integral yields a precise prediction for three distinct wavepacket regimes corresponding to one-phase, two-phase, and Airy modulations. The analysis of DSW implosion (regions c and f in Figs. 1 and 2, respectively) can be viewed as a weakly nonlinear generalization of the linear wavepacket analysis.
The dynamics of DSW implosion (region f ), incoherent waves (region g), and the dispersive Lax shock (region h) are qualitatively similar to the respective regions in the large transition width case identified in Fig. 1.
A general conclusion from our analysis, in both the large and small transition width regimes, is that the loss of linear dispersion convexity underlies the generation of nonclassical wave patterns.

Linear wavetrains and Airy modulation near the zero dispersion point
In this section, we investigate the asymptotic structure of the approximately linear wavetrains generated in BBM dispersive Riemann problems with smoothed step initial data (1.7) in the parameter domain (u − , u + , ξ) defined by: 0 < u − , µ < u + /u − < σ 2 (ξ), where µ = e −3/2 /4 and σ 2 (ξ) is defined later in equation (5.31); see regions d and e in Fig. 2. The significance of σ 1,2 will be clarified later. Depending on the relative values of u − and u + , the initial step generates either a RW (u − < u + ) or a DSW (u − > u + ). In both cases, an approximately linear wavetrain develops on the background u = u − behind the RW or DSW and occupies an expanding region x ∈ [x − 0 (t), x + 0 (t)], where the right boundary x + 0 coincides with the trailing edge of the RW or DSW and x − 0 will be determined. In this section, we use linear theory to analyze the corresponding wavetrains, hence we refer to them as "linear wavetrains", with an understanding that this is an approximation.  The triangle is a partitioning of the boundary data in the u + + u − ≥ 0 half plane. Corresponding representative numerical simulations for ξ = 0.1 in each case are shown. The table is a description of the partitioning and wave patterns. µ = e −3/2 /4 ≈ 0.056, σ 2 (ξ) = 1 + √ 6βξ 1/2 + 27β/8ξ 5/2 + · · · , β ≈ 9.5, σ 3 (ξ) is generally unknown but empirically σ 3 (0.1) ≈ −0.70. σ 1 (ξ) is generally unknown but empirically σ 1 (0.1) ≈ 0.68.

Linear wavetrains in the BBM Riemann problem
We start by linearizing the BBM equation (1.4) about a constant background u > 0, u(x, t) = u + ϕ(x, t), |ϕ| u, and obtain , we obtain the BBM linear dispersion relation: Then the group velocity is given by the first derivative of ω 0 and the dispersion sign is defined by the sign of the second derivative The dispersion relation is nonconvex and the dispersion sign is zero for k = 0, k = √ 3 or u = 0. We stress that it is the non-convexity of the linear dispersion relation (3.2) that gives rise to crucial differences between the BBM and KdV dispersive Riemann problems. The KdV linear dispersion relation is strictly convex for k > 0.
The group velocity (3.3) has a minimum when k = √ 3, and a maximum at k = 0, Consequently, a BBM linear wavetrain that develops from a localized initial disturbance considered in isolation on the constant background u is confined to the region x ∈ [s min t, s max t] in the asymptotic regime t 1, although some rapidly decaying oscillations are possible outside this region. We now place this in the context of the dispersive Riemann problem where the linear wavetrain is generated on the background u = u − as part of the nonlinear-dispersive regularization of smoothed step initial data. We consider two cases u − < u + and u − > u + separately, distinguishing between a linear wave trailing a RW or a DSW (see Fig. 2(d,e)).
(i) First, consider the case u − < u + that leads to a RW. As we shall see later in Sec. 4.2, the trailing edge of the RW is located at x = u − t = s max (u − )t. Provided u + /u − < σ 2 (ξ), no solitary wave propagates in the region x < s max t (cf. Sec. 4.2). Thus the two kinds of structures generated in this class of BBM dispersive Riemann problems-the linear wavetrain for x ∈ [s min (u − )t, s max (u − )t], and the RW for x ∈ [u − t, u + t]-do not overlap and can be described separately, except in the vicinity of the point x = u − t where the two solutions could be matched. A linear wavetrain can also be generated for u + /u − ≥ σ 2 but we do not consider this case because one or more solitary waves may accompany the wavetrain, complicating its analysis.
(ii) For the case u − > u + , which leads to a DSW, we show in Sec. 6 that the trailing edge of the DSW is located at x = s − t < s max (u − )t so that the trailing linear wavetrain is necessarily attached to the DSW. For u + /u − > µ, we have s min (u − ) < s − and our main concern here will be the description of the linear wavetrain located outside the vicinity of the DSW trailing edge point x = s − t, where the two modulated wavetrains could be matched.
Summarizing, we assume that for dispersive Riemann problems with µ < u + /u − < σ 2 (ξ), the generation and propagation of a small-amplitude wavetrain is governed by the linearized BBM equation (3.1) with u = u − , and its dynamics are fully decoupled from the nonlinear DSW or RW dynamics. In other words, we suppose that the initial profile's Fourier spectrum can be partitioned into two essentially distinct components. The short-wave (large wavenumber) component of the spectrum is assumed to evolve according to the linear wave equation (3.1). The long-wave (small wavenumber) component of the spectrum is assumed to evolve according to the nonlinear BBM equation (1.4). Consequently, the linear wave evolution is superposed with the nonlinear evolution. We justify our partition of the initial Fourier spectrum by comparison of our linear analysis with numerical simulations. Since the Fourier transform of the initial condition (1.7) is given by the formula, the partitioning assumption implies that the Fourier transform of the linear wave ϕ(x, t) at t = 0 is approximated byφ , ∀k ∈ R. (3.8) The initial condition (3.8) simplifies for a step (ξ → 0) toφ 0 (k) = (u + − u − )/ik, k > k 0 . Note that |φ 0 (k)| is a rapidly decaying function of ξ, and so no linear wavetrain is expected to be generated for ξ 1. This agrees with Fig. 1, which shows the evolution for the dispersive Riemann problem with ξ = 10.
The solution of (3.1), is given by For notational convenience, we suppress the dependence of χ on x t .

Stationary phase analysis
Consider the asymptotic regime t → ∞ with x/t fixed. We can evaluate the integral (3.9) using the stationary phase method, cf. for instance [51]. Fast oscillations average out in the integral, and the main contribution comes from the neighborhood of the stationary points k n defined as solutions of χ (k) = 0, namely (3.10) A first approximation of ϕ(x, t) is obtained by substituting the Taylor expansions: Equation (3.10) implies that k n depends on s ≡ x/t. The computation of (3.9) thus reduces to the integration of a Gaussian, and the corresponding solution ϕ only depends on the number of stationary points and their respective positions in k-space: Real solutions of (3.10) come in pairs (−k n , k n ) and, sinceφ 0 (−k) =φ * 0 (k) and ω 0 (−k, u − ) = −ω 0 (k, u − ), the sum in (3.12) simplifies to [51]: The stationary phase equation (3.10) divides space-time into two regions. (a) In the region s max > s = x/t > 0, denoted hereafter as Region I, (3.10) has only one positive solution so the series (3.13) has just one term. This is comparable to what we would obtain with a convex linear dispersion relation, and we refer to this region as the convex linear regime. In fact, for k 1, we obtain the linear dispersion relation ω 0 (k, u) ∼ u(k − k 3 ), which is the linear dispersion relation for the KdV equation.
(b) In the region s min < s = x/t < 0, denoted hereafter as Region II, (3.10) has k 1 (s) and an additional positive solution so that the series (3.13) has two terms. The coexistence of two waves with wavenumbers k 1 and k 2 at the same position x within Region II is a direct consequence of the nonconvexity of the dispersion relation (3.2). The modulation of the linear wave in this region dramatically differs from a convex-dispersion, KdVtype linear modulation. In Sec. 6.3 below, we show that the coexistence of two dominant wavenumbers also persists in the nonlinear regime, leading to the phenomenon of DSW implosion. The upper panel of Fig. 3 displays the comparison between the wavenumber of the waves obtained numerically from the full partial differential equation (PDE) simulation, shown in the lower panel, and the graph of x given by (3.10) as a function of k = k n for a fixed value of t, with k n , n = 1, 2 given by the formulas (3.14), (3.15). The two regions I and II are treated separately using the stationary phase approximation (3.13). A third region, denoted III in Fig. 3 and defined hereafter as x s min t, is resolved with an expansion near x = s min t. The solution is well-approximated by a one-phase linear wave in region I and a two-phase linear wave in region II. The variation of the waves' wavenumber extracted from the numerical solution is represented by black disks in the upper panel; the procedure to extract wavenumbers in region II is detailed in Appendix A.2. The blue (red) dashed line corresponds to the variation of the analytical solutions k 1 (x, t) (k 2 (x, t)).
(Region I) In the convex linear regime of Region I, the stationary phase approximation yields This solution is not valid close to the so-called caustic trajectory x = u − t = s max t, where the two dominant wavenumbers k 1 and −k 1 coalesce at k 1 = 0, (cf. for instance [9]) and ∂ kk ω 0 (k 1 , u − ) = 0. At this edge of Region I, one needs to consider the next order in the expansion of χ(k) in the Taylor expansion (3.11).
Besides, as pointed out earlier, a RW or a DSW develops close to the point x/t = u − or x/t = s − < u − (respectively), and the linear wave has to be asymptotically matched with the hydrodynamic state. Such a higher order derivation, which can be achieved by following the matched asymptotics procedure developed for the KdV equation [22,33], is not important for us here because we have a well-defined initial value problem (3.1), (3.8). The comparison between the approximate solution (3.17) and the numerical solution of the dispersive Riemann problem in region I is displayed in Fig. 4.
(Region II) In Region II, for which x ∈ [s min t, 0], the modulated wave corresponding to the large wavenumber branch k 2 coexists with the wave from the lower wavenumber branch k 1 , leading to a modulated beating pattern. The stationary phase approximation then yields the superposition In this region, the expansion (3.11) is insufficient. We shall denote the union of this special region s ∼ s min and s < s min as s s min and call it Region III. In order to describe the oscillations in this region, we expand χ(k) close to the inflection point k = √ 3 (3.19) Substituting this expansion into eq. (3.9), we obtain where Ai(y) is the Airy function defined by [2] Ai(y) The solution (3.21) represents a harmonic cosine wave modulated by the Airy function. Note that the Airy approximation (3.21) is also valid for x/t < s min , where the Airy modulation describes an exponential decay of the wave's amplitude: Ai(y) ∼ y −1/4 exp(− 2 3 y 3/2 ) as y → +∞, which in our case translates to (x − s min t) → −∞.

Uniform Airy approximation
Using the technique originally developed in [8], we derive a uniform approximation that is valid in the union of regions II and III whose limiting behavior is given by (3.21) when x s min t and (3.18) when x s min t. We summarize here the derivation detailed in [9]. First, we suppose that the solution in region II accords with the ansatz where φ 0 ∈ C, φ 1 ∈ C, A ∈ R and ζ ∈ R are functions of (x, t) that are to be determined. By selecting the exponential's phase ψ = A − κζ + κ 3 /3, we will ensure that it has the same stationary points as χ(k) for x close to s min t. Note that, although only the leading order termφ 0 (k n ) in the Taylor expansion of ϕ 0 (k) near k = k n in the approximations leading to (3.13) and (3.20) was used, incorporating the first order termφ 0 (k n )(k − k n ) in the analysis is essential for constructing a uniform approximation. We first determine the functions φ 0 , φ 1 , A and ζ by evaluating the integral (3.23) using the stationary phase method. For simplicity, we consider the computation for ζ > 0; the computation is similar in the case ζ < 0, cf. [9]. The phase ψ has two stationary points: Close to the stationary points ± √ ζ, the phase is given by Substituting the Taylor expansion of ψ into (3.23), we obtain (cf. previous stationary phase computations) Since the solution (3.23) should be asymptotically valid in all of region II, we impose that (3.26) identify with (3.18) for x s min t. The identification imposes Then, direct integration of (3.23) yields: where Ai(y) is given by (3.22) and Ai (y) ≡ dAi/dy. By construction, (3.28) is asymptotic to (3.18) when , such that the uniform approximation (3.28) asymptotically matches with the Airy approximation (3.21) when x ∼ s min t (recall that α 1 in the long time regime). Note that the uniform approximation (3.28) is also valid in the region x < s min t, or equivalently ζ < 0, where the analytical expressions (3.14) and (3.15) are complex.

III
II I Note that, according to (3.8), the linear wave depends on the sign of u + −u − and the green line represents the field −ϕ(x, t) with u − > u + in order to be compared to the solution ϕ(x, t) with u − < u + .
To summarize, the linear wave generated by the dispersive Riemann problem is approximated by the piecewise-defined function where ϕ I and ϕ uni are defined in (3.17) and (3.28), respectively. This stationary phase approximation compares well with the numerical solution, as displayed in Fig. 4. An important consequence of the approximate initial data (3.8) for ϕ is that the linear wave solution ϕ(x, t) is valid for both polarities of the initial jump u − < u + and u − > u + , i.e., for either the generation of waves accompanying a RW or the generation of waves accompanying a DSW, respectively. Moreover, the linear waves accompanying a RW, ϕ RW (x, t), with jump u + − u − > 0 are related to the linear waves accompanying a DSW with the same jump magnitude u − − u + > 0 by a minus sign: ϕ DSW (x, t) = −ϕ RW (x, t). Figure 5 displays the numerical solutions for two different initial conditions u − < u + and u − > u + with a common difference magnitude |u + − u − |. Although the two numerical solutions are comparable, the discrepancy between them increases with increasing x; the discrepancy is larger in Region I. In fact, the small amplitude assumption at the core of linear theory is not well-fulfilled since, initially, |ϕ(x, 0)| = O(|u − −u + |). It is, therefore, remarkable that linear theory still predicts the evolution of the fast oscillations with excellent accuracy, even if the small amplitude assumption is initially violated. The discrepancy in Fig. 5 is due to different matching conditions for the linear wave with the leftmost, trailing edge of the RW or DSW. This asymptotic matching has been investigated for the KdV equation in [33,22].
Finally, we note that the asymptotic behavior (3.21) near the zero dispersion point is general and can be obtained for other wave equations exhibiting a non-convex linear dispersion relation. For example, it was obtained (with appropriate modifications due to a different form of the dispersion relation) in [49] for linear wavepackets in the Gardner-Ostrovsky equation.

Shocks and rarefactions
In this section, we consider shock wave and rarefaction wave solutions of the underlying conservation law

Stationary Lax shocks
We first consider the solution when u − = −u + = 1. Both the Hopf u t + uu x = 0 and the BBM equation (1.4) admit the stationary, discontinuous, weak solution or stationary Lax shock. The shock wave is compressive in the sense that the Hopf characteristics (with speed u ± = ∓1) on each side of the shock located at x = 0 propagate into the shock [32]. Although this BBM solution has been acknowledged elsewhere, e.g. [15], it has remained mostly a curiosity. We investigate the dynamics when (4.1) is smoothed at t = 0 according to (1.7) with u ± = ∓1. The numerical solution u(x, t) is displayed in Fig. 6 at different times. Evolution initially leads to the usual nonlinear self steepening (t = 10) due to broad initial data (ξ = 10). This long wave, Hopf equationlike evolution, is soon accompanied by oscillations (t = 20) when the solution is steep enough and the BBM dispersive term u xxt is important. Dispersive regularization of shock formation typically results in a DSW (cf. Sec. 6), with an expanding collection of oscillations. The long-time dynamics of a generic DSW consists of finite amplitude waves that are bounded in wavenumber [16]. But here, oscillations of increasing, apparently unbounded wavenumber develop as time increases, concentrating near the origin. The last panel of Fig. 6 (t = 300) shows that the large t behavior for both broad and narrow initial smoothed step profiles approach approximately the same structure, suggesting the existence of some self-similar asymptotic configuration. We remark that a previous numerical study of BBM has depicted similar behavior for broad, periodic initial data [21], associating it with discontinuity formation.
To further investigate these dynamics, we introduce the short space ρ = x/ and long time T = t scaling where 0 < 1 is a small parameter representing the characteristic length scale of the oscillations in a neighborhood of the origin. With this scaling, the BBM equation (1.4) becomes Expanding u in as u = U 0 (ρ, T ) + U 1 (ρ, T ) + · · · , we obtain the leading order equation This equation can be integrated once in ρ, yielding a nonlinear Klein-Gordon equation In general, α is an arbitrary function of T . Motivated by the numerical observations in Fig. 6, we seek a self-similar solution that must be independent of the small but otherwise arbitrary parameter , thus can only depend on ρT = xt: In order for g to satisfy a well-defined ordinary differential equation (ODE) If the integration constant is set to α = 1/2, then h(η) = 0 is a stable fixed point of the homogeneous equation (4.7) with the general, Bessel-type solution h(η) Consequently, the large η asymptotics of the Bessel solution demonstrate that the sought solution to eq. (4.6) with α = 1/2 exhibits algebraic, oscillatory decay to the requisite boundary condition for some constants C and φ.
Since the symmetric, smoothed step initial data considered here is an odd function of x, we seek a solution to (4.6) that is an odd function of its argument so that g(0) = 0. Evaluating eq. (4.6) at η = 0 also determines g (0) = −1/2. The initial conditions uniquely determine the solution of eq. (4.6) with α = 1/2 that decays to −1 according to (4.8). We compute the numerical solution to this initial value problem and display the result in Fig. 7. Motivated by the linearized equation (4.7), we also plot in 7 the function −1 , which provides a good empirical fit to g(η) when 2 √ η 8. By the asymptotics of the Bessel function J 0 , we can infer from Fig. 7 that C ≈ π −1/2 and φ ≈ π/4 + 0.7 in eq. (4.8).
In Fig. 8, the numerical solution of the dispersive Riemann problem for both narrow and broad initial data approaches the self-similar solution g(xt) (with an odd extension) for large t. The narrow case in 8a leads to the development of oscillations with wavelengths that are shorter than the self-similar profile whereas the broad case in 8b exhibits longer wavelength oscillations than g(η). In both cases, there is no indication of discontinuity formation in finite time.
Since the odd-extended self-similar solution g(η) → ±1 as η → ∓∞, we observe that it converges to the stationary Lax shock (4.1) as t → ∞ for each fixed x (4.10) Figure 7: Self-similar dispersive Lax shock profile (solid, black) and a Bessel function approximation (dashed, red).
although the convergence is not uniform in x. The self-similar, oscillatory profile g(xt) describes how the stationary Lax shock develops as t → ∞. Based on these observations, we call the profile g(xt) a dispersive Lax shock. This asymptotic solution represents a completely new type of dispersive shock structure.
The simulations in Figs. 6 and 8 exhibit behavior reminiscent of Gibbs phenomenon in the theory of Fourier series. Perhaps a more apt comparison is to initial value problems for linear dispersive PDE in which the data is discontinuous [7]. For a broad class of constant coefficient, linear dispersive PDE, the convergence of the solution as t → 0 + for initial data with a discontinuity is not uniform. Moreover, it exhibits Gibbs phenomenon whereby the solution has an overshoot that, in proportionality to the jump, converges to the Wilbraham-Gibbs constant In contrast, the dispersive Lax shock g(xt) converges to the discontinuity (4.10) as t → ∞. While the oscillations for η sufficiently far from 0 are linear and Bessel-like, the compression of these oscillations is an inherently nonlinear process. Narrower oscillations appear as the initial smoothed step profile steepens and approaches a jump discontinuity. In contrast to the aforementioned linear initial value problems, we observe an oscillatory overshoot that can be numerically estimated from the self-similar profile g(xt) to be, in proportion to the jump height, a bit more than twice the Wilbraham-Gibbs constant g. Numerical simulations are limited by the discretization and method utilized. The structure and selfsimilar scaling of g(xt) will, for large enough t, exceed the resolution of any fixed discretization scheme. Our analysis of this dispersive Riemann problem provides a means to interpret solutions in the vicinity of dispersive Lax shocks.
We show in Sec. 6.3 that the solution of the dispersive Riemann problem for slightly asymmetric boundary conditions 0 < u − + u + u − is drastically different from the symmetric case investigated here. In particular, we will show that solitary waves are shedded in an incoherent fashion on top of a compressive shock structure.

Rarefaction waves
From now on, we consider the case u + = 1. In Fig. 9, we show numerical solutions of (1.4), (1.7) with u − < 1. Dispersive effects are negligible for smoothed step initial data (ξ 1), and the asymptotic solution of the dispersive Riemann problem is a smoothed modification of the RW solution   u − = −1, the RW develops with an expansion shock that is described in Sec. 4.3 below. The evolution of asymmetric smoothed step initial data (having −1 < u − < 1) displays richer structure (cf. Fig. 9a). The interplay between long and short waves in the linear regime (clearly visible in Fig. 9a with u − = 1/2) is analysed in Sec. 3, and the emergence of a train of solitary waves is investigated in Sec. 5. We show in Sec. 5.2 that the development of expansion shock and solitary wavetrain structures for ξ → 0 is linked to the atypical early time evolution of the dispersive Riemann problem.

Expansion shocks
with ξ 1. In [15], an asymptotic solution is derived in detail, using ξ as a small parameter, with long-time scaling T = ξt. Here, we outline the description in [15].
The inner structure of the solution, near x = 0, is described using the same short space scaling χ = x/ξ as that utilized for the description of the dispersive Lax shock in Sec. 4.1. However, the small parameter ξ in this case is fixed by the choice of initial data (1.7). These scalings introduced into the BBM equation (1.4) lead to , we obtain the leading order equation In addition to the self-similar solution to this equation that we obtained for the dispersive Lax shock, here we obtain a separated solution describing the expansion shock. Equation (4.16) with initial condition U 0 (χ, 0) = tanh (χ) corresponding to (4.14), admits the separated solution Thus the inner solution of the dispersive Riemann problem (1.4), (4.14) for ξ 1 is given by The outer structure of the solution is described using the long space scaling X = ξx, leading to Looking for u in the formŨ 0 (X, T ) + ξŨ 1 (X, T ) + O(ξ 2 ), we obtain the Hopf equation at leading order The general solution of this equation This formula is continuously matched with the far-field conditions: lim x→±∞ u = ±1, leading to the outer solution, Combining the descriptions (4.18) and (4.22), we obtain the uniformly valid asymptotic solution of (1.4), (4.14) for ξ 1 As remarked in [15], the accuracy of the asymptotic solution (4.23) is excellent, cf. Fig. 10. Note that the expansion shock structure converges to the rarefaction solution (4.13) as t → ∞. We also observe in Fig. 10 that the expansion shock persists for slightly asymmetric boundary conditions with −1 < u − < 0, and u + = 1, see for instance the solution for u − = −0.7. The outer solution for such conditions is a slight modification of (4.22) which is continuous at x = u − t. Thus the truncated solution reads:

Solitary wave shedding
In this section, we study the generation of solitary waves observed in numerical solutions. Throughout this section, we keep u + = 1 and consider −1 < u − < 1, with ξ 1 (regions a, b, c, and d in Fig. 2). Numerical results, described in Sec. 5.1, exhibit a train of solitary waves for a range of values of u − . However, for x > 0, the solutions have either a rarefaction wave accompanied by a region I linear wave (described in Sec. 3), or a portion of the expansion shock solution (described in Sec. 4.3). There are no solitary waves generated in the region x > 0. We can therefore restrict attention to x < 0.
In Sec. 5.2, we explain the appearance of solitary waves by analyzing the early time behavior of the solutions. In Sec. 5.3, we use the energy equality satisfied by smooth solutions, to establish the values of u − for which a solitary wave is generated from those for which the solution has no solitary wave.

Numerical Solutions
We start with a qualitative description of the numerical solutions shown in Fig. 12 for various values of u − < 1 and fixed ξ = 0.  (1.4) is given by [40] u( where c is the wave speed and u is the background value. The wave amplitude is a = 3(c − u) yielding the speed-amplitude relation c = c(a, u) = u + a/3.  If u = u − = 0, the linear dispersion relation (3.2) is ω 0 (k, u) = 0, and only a train of solitary waves is generated for x < 0 as in Fig. 12b. Similarly, when u − < 0, although a linear wave can propagate on a strictly negative background, its group velocity satisfies ∂ k ω 0 > u − so that no linear waves remain in the region x < u − t for t 1. The solitary wavetrain then propagates from the tail of the truncated expansion shock solution (4.25), as depicted in Fig. 12b. No solitary wave emission is observed when u − < u − ] is depicted in Fig. 13. We empirically observe that when u − ∈ (u (a) , 0] and multiple solitary waves are emitted, the amplitude a(x, t) of the solitary wavetrain asymptotically satisfies c(a, u − ) = x/t where c is the speed-amplitude relation (5.2). A detailed study of the solitary wavetrain using, for example, Whitham modulation theory [51], is left as a subject for future work.

Early time evolution
The unusual generation of linear waves and solitary waves can be qualitatively explained by investigating the early time dynamics of the dispersive Riemann problem. We consider in this section the small-time expansion of the solution where u 0 (x) is given by (1.7) with u + = 1. Substituting (5.3) into (1.4), we obtain at the zeroth order in t where ≡ d/dx. The small parameter condition tu 1 (x) u 0 (x) yields the boundary conditions: lim x→±∞ u 1 (x) = 0. Let G(x) be the Green's function satisfying, The solution of (5.4) is then Note that u 1 (x) (and higher order terms) can be directly obtained from the integro-differential form of the BBM equation (1.8) by iteration, which is part of the contraction mapping argument used to prove local existence of solutions to the BBM initial value problem [4]. The correction u 1 (x) is bounded by an estimate that decreases with increasing ξ, so that u 1 is negligible for large ξ.
On the other hand, in the ideal case ξ → 0, the initial condition is the step given by (1.6), and (5.7) simplifies to leading to a decrease in u(x, t) for sufficiently small t > 0. This introduces a dip in the graph of u(x, t) for x < 0. The small amplitude assumption tu 1 (x) 1 − u − (recall that 1 − u − represents the initial smoothed step height) then reduces to t 1/(1 + u − ). Figure 14a shows good agreement between the approximate solution (5.3) with (5.9) and the corresponding numerical solution of the dispersive Riemann problem when ξ = 0.1. While the initial smoothed step is rarefying in the region x > 0, a dip develops  immediately in the region x < 0. As a result, the initial smoothed step transition persists, contrasting with numerical solutions for ξ = 10 ( Fig. 14b) that more closely resemble the classical RW, with no dip developing. Figure 15 displays different evolution scenarios of smoothed step initial data for ξ = 0.1 depending on the value of u − . If u

The threshold for generation of a negative solitary wave
In this subsection, we derive the conditions for and properties of the emitted solitary wave in the simplest case where only one solitary wave is shedded, cf. Fig. 13. We also determine the threshold u − = u Integration of (5.10) and (5.11) between the limits x = −L and x = +L yields 14) Because of the term u 2 x , the initial energy E(t = 0), obtained by substituting the initial smoothed step (1.7) in (5.15), strongly depends on the transition width ξ, which plays a crucial role in the appearance of negative solitary waves, as was demonstrated numerically in Fig. 9.
Motivated by observations in Sec. 5.1, we suppose that when u − is very close to u − , the solution u(x, t) can be approximated for large t by where b and c remain to be determined. Expressions for u RW and u sol are given in (4.13) and (5.1), respectively, and ϕ(x, t) is a correction term. We suppose here that the RW is not centered and denote b the value of the shift. The incorporation of the RW shift is crucial for the precise determination of the threshold for solitary wave emission and the accompanying solitary wave's speed c. Figure 17a shows the necessity of incorporating a vertical shift −b/t of the centered RW solution u RW (x, t) in order to achieve good agreement with the numerical simulation. We also hypothesize that the dominant contribution to the correction ϕ(x, t) is a linear wavetrain, analogous to that studied in Sec. 3. Note that the solitary wave and the linear wavetrain correction propagate on the background u = u − , fixed by the boundary condition. Thus, the condition c < s min ensures that the solitary wave is negative and that it is wellseparated from the RW and the linear wavetrain for t 1. This assumption is consistent with numerical computations. We now show that the formula (5.16) represents a solution of the dispersive Riemann problem (1.4),(1.7) for a certain value of the solitary wave speed c depending on ξ 1 and u − . In a first approximation, we suppose that the correction ϕ(x, t) has a negligible mass and energy contribution: We evaluate the mass N (t) and the energy E(t) in two different ways: i) by direct substitution of (5.17) in (5.13) and (5.15) and ii) by solving (5.12) and (5.14) with the initial condition (1.7) and the boundary conditions lim x→±∞ u(x, t) = u ± , lim x→±∞ u xt (x, t) = 0. Equating the two expressions for N (t) and E(t) then yields an approximation to the solitary wave speed c and the RW shift b.   18) and the direct substitution of the ansatz (5.17) into (5.15) yields We choose L and t such that: L t 1. The contribution of the RW to the energy then simplifies to Since L 1 and u sol (x, t; c) decays exponentially to 0 when |x| → ∞, we can allow L → ∞ in the integral for N sol and E sol (ii) Now we consider the conservation laws (5.12) and (5.14). Since L t, the RW and the solitary wave have not reached the positions x = ±L at the time t and the solution u(x, t) decays rapidly to u ± for |x| ∼ L. In particular, we have u ∼ u ± and |u xt | u 2 ± at x = ±L. Thus the conservation laws read where u 0 is the smoothed step initial data (1.7). We calculate Since L ξ, the initial energy E(0) simplifies to which combine into a single equation for the solitary wave speed c Since ξ 1, we have the approximation The critical value u  Discrepancy between the analytical results and the numerics in Fig. 18 is due to neglecting the linear corrective term ϕ in the ansatz (5.16). Because it separates from the RW and solitary wave, the contribution of the linear wave to the mass and the energy are dx. comparable. Thus, neglecting ϕ in (5.16) overestimates the solitary wave mass and energy. An improved estimate for the solitary wave amplitude is possible if the variation of ϕ(x, t) for x ∈ (s min t, u − t) can be determined. While we have obtained an asymptotic description of linear wavetrains that are generated by smoothed step dispersive Riemann problems (cf. eq. (3.29)), our analysis did not include the generation of a solitary wave, which is a significantly more complex problem. In particular, the initial data for the linear wavetrain (3.8) would require modification due to initial solitary wave-linear wavetrain interaction. Therefore, we have not found a simple analytical estimation of N ϕ and E ϕ . We emphasize that both computations presented in this section predict solitary wave emission for u − close to the threshold u − σ 2 (ξ) = 1. Indeed, we observe in the previous section that the numerical solution involves the generation of multiple solitary waves for sufficiently small values of |u − |, which is not taken into account by the ansatz (5.16) considered here. For example, multiple solitary waves are generated for u − < 0.13 when ξ = 0.1.
A similar argument could be used to explain the emission of solitary waves for u (1) − ≤ u − < 0, where the solution could be approximated by with u aexp given by (4.25) and ϕ a corrective term. Contributions to the mass and the energy from the correction term are of the same order as the contributions from the solitary wave. The determination of ϕ is, as demonstrated in the previous computation, necessary to compute the solitary wave's amplitude. Typical variations of ϕ(x, t) within the interior region of the asymmetric expansion shock x ∈ (u − t, t) are displayed in Fig. 17b. Contrary to the previous computation, we have not found a simple ansatz that can describe the variation of ϕ(x, t). The determination of the solitary wave's amplitude for u − < 0 is left as a subject for future work.

Dispersive Shock Waves
DSWs are expanding modulated nonlinear wavetrains regularizing wavebreaking singularities in dispersive hydrodynamics, a conservative counterpart of regularizing classical shock waves in viscous fluid dynamics. The shock structure of a DSW is more complex than the shock structure of a viscous shock wave. In particular, a DSW cannot be described by a traveling wave solution of the nonlinear wave equation [11]. We refer the reader to the recent review [16] for the principal ideas and applications of DSW theory. In this section, we analyze the DSWs generated in the BBM dispersive Riemann problem for |u + | < u − ; see regions b and e in Figs. 1 and 2, respectively. We first consider the region 0 < µu − < u + < u − , where µ = e −3/2 /4. This restriction on the initial data (the formula for µ is derived later) is related to DSW admissibility conditions specific to the BBM equation. As an aside, we note that for the KdV equation (1.3), the only requirement on the Riemann data for a DSW solution is u − > u + . Although the structure of the BBM DSW within the admissibility region (6.1) is qualitatively similar to that of the KdV DSW, the quantitative description of BBM DSWs requires a modified approach due to non-integrability of the BBM equation. The DSW fitting method introduced in [12] (see also [16]) enables the analytical determination of the DSW edge speeds in non-integrable dispersive-hydrodynamic systems. A recent extension of DSW fitting [10] also enables partial determination of the interior DSW structure. In Section 6.2, we apply the extended DSW fitting method to the BBM equation and compare the obtained analytical results with numerical simulations of the smoothed dispersive Riemann problem (1.7). In the complementary region −u − < u + < µu − of Riemann data, DSW fitting fails due to the development of nonlinear two-phase oscillations at one of the DSW's edges. This phenomenon, termed DSW implosion, is a nonlinear counterpart of the two-phase interference pattern occurring in linear wavetrains near the zero dispersion point, described earlier in Section 3. Numerical simulations and the onset of DSW implosion were first studied in [35] for the viscous fluid conduit equation. In Section 6.3, we develop an analytical framework for the structure of DSW implosion by deriving a coupled NLS equation for a nonlinear superposition of two BBM Stokes waves that asymptotically describe the nonlinear twophase dynamics of the imploded region.

Classical, convex DSWs: general properties, fitting relations and interior modulation
It is convenient to describe properties of DSWs in the framework of a fairly general scalar dispersive equation where f (u) is the hyperbolic flux and D[u] is a dispersion operator, giving rise to a real valued linear dispersion relation ω = ω 0 (k, u). For the BBM equation, f (u) = 1 2 u 2 , D[u] = u xt and ω 0 (k, u) is given by (3.2).
The rapidly oscillating structure of DSWs-see, e.g., Fig. 19-motivates the use of asymptotic, WKBtype, methods for their analytical description. One such method, known as Whitham modulation theory [51] (see also [28]), is based on the averaging of dispersive hydrodynamic conservation laws over nonlinear periodic wavetrains, leading to a system of first order quasilinear PDEs, known as the Whitham modulation system. Classical DSW theory has been developed for KdV-type equations. More generally, it is useful to define the dispersive hydrodynamic equations (6.2) to be classical or convex if the associated Whitham modulation system has the properties of strict hyperbolicity and genuine nonlinearity; see [16]. For the scalar equation (6.2), convexity of the flux, i.e., f (u) = 0, and of the linear dispersion relation, i.e., ∂ kk ω 0 (k, u) = 0 for k = 0 and u ∈ R, are defining conditions for convexity of the dispersive hydrodynamic equation [16]. This convexity property provides for the existence of a class of solutions to the Whitham modulation system that describe the DSW structure.
For the BBM equation, the hyperbolic flux is indeed convex, but the linear dispersion relation is not, because ∂ kk ω 0 (k, u) = 0 at the zero dispersion points k = √ 3 and u = 0 (additionally k = 0), see (3.4). As a result, the classical, convex DSW regime is realized only for a restricted domain of Riemann data u ± given by (6.1), and established below.
Whitham theory has proven particularly effective for the description of DSWs in both integrable and non-integrable systems. If the dispersive hydrodynamics are described by an integrable equation such as the KdV equation, the associated Whitham system can be represented in a diagonal, Riemann invariant form [50,20,28]. This fact enabled Gurevich and Pitaevskii (GP) [25] to construct an explicit modulation solution for a DSW generated by a dispersive Riemann problem for the KdV equation. The GP construction is based on a self-similar, rarefaction wave solution of the KdV-Whitham system of three modulation equations for three slowly varying parameters that locally describe the oscillatory wave. For our purposes, it is convenient to choose the mean u, the wavenumber k and the amplitude a as the three parameters, all expressible in terms of the Riemann invariants. The interior shock structure of a DSW is then described by a self-similar, centered solution of the modulation equations whose values, (u, k, a)(x, t) for t 1, lie on a 2-wave rarefaction curve connecting the edge points (u, k, a) = (u − , k − , 0) and (u, k, a) = (u + , 0, a + ). These edge points are associated with a harmonic, small-amplitude wave (a → 0) and a solitary wave (k → 0), respectively. The two integrals associated with the 2-wave curve generate a family of Riemann problems that are parametrized by (u − , u + ), i.e. k − = k − (u − , u + ) and a + = a + (u − , u + ). We emphasize that the GP, self-similar solution results from a genuine Riemann problem for the Whitham modulation equations, i.e. discontinuous initial data, in contrast to an initial smoothed step profile like (1.7) posited for the dispersive Riemann problem. For convex dispersive hydrodynamic problems, this distinction is not so important because the long-time evolution of a smoothed step profile is described by the GP solution. But for nonconvex problems such as BBM, there are subtleties associated with the smoothing length scale ξ-possibly beyond the application of leading order Whitham modulation theory-that need to be considered.
A DSW in a general convex, scalar (i.e., unidirectional) dispersive hydrodynamic system such as (6.2) has a multi-scale structure that is similar to the KdV DSW. It consists of an oscillatory transition between two non-oscillatory-e.g., slowly varying or constant-states: one edge is associated with a solitary wave that is connected, via a slowly modulated periodic wavetrain, to a harmonic, small-amplitude wave at the opposite edge. The relative position (left/trailing or right/leading) of the soliton and harmonic edges determines the DSW orientation d ∈ {±1}, where d = +1 if the solitary wave is on the right/leading edge. DSW orientation is determined by the dispersion sign as d = −sgn[∂ kk ω 0 (k, u)] [16,11]. The DSW polarity p ∈ {±1} is defined by the polarity of the solitary wave edge, where p = +1 if the solitary wave is a wave of elevation relative to its background. DSW polarity is given by p = −sgn [∂ kk ω 0 f (u)]. The DSW shown in Fig. 19 has d = p = +1.
By direct evaluation of the BBM linear dispersion relation (3.2), we observe that the BBM equation is expected to support DSWs with positive orientation/polarity d = p = +1 and negative orientation/polarity d = p = −1. These two kinds of DSWs map to each other through the BBM invariant transformation x → −x, u → −u. Hence, it is sufficient to only consider DSWs with d = p = +1.

DSW fitting relations
For non-integrable dispersive equations such as the BBM equation, diagonalization of the associated Whitham modulation system in terms of Riemann invariants is generally not possible. Thus, the explicit determination of the Whitham system's simple wave solution corresponding to a DSW is problematic, even though its existence requires only strict hyperbolicity and genuine nonlinearity in the 2-wave characteristic family. One can, however, explicitly determine key observables associated to each DSW edge for integrable Figure 19: KdV-type DSW in convex dispersive hydrodynamics with negative dispersion. The DSW fitting relations determine the soliton edge velocity s + and amplitude a + as well as the harmonic edge velocity s − and wavenumber k − . and non-integrable equations provided the dispersive hydrodynamics are convex. These observables include the DSW edge speeds s ± and their associated wave parameters-the harmonic edge wavenumber k − and the soliton edge amplitude a + , as functions of the initial data u ± . The determination of these observables represents the fitting of a DSW to the long-time dynamics of dispersive Riemann initial data. The DSW fitting method proposed in [12] (see also [16]) is based on a fundamental, generic property: the Whitham modulation equations admit exact reductions to a set of common, much simpler, analytically tractable equations in the limits of vanishing amplitude and vanishing wavenumber, which correspond to the harmonic and soliton DSW edges, respectively. DSW fitting has also been successfully applied to the propagation of a broad localized pulse into a constant state [14,29,36].
We now formulate the set of DSW fitting relations for an initial step in the mean u, a dispersive hydrodynamic analogue of the classical Rankine-Hugoniot jump conditions for viscously regularized shocks. As outlined above, the DSW fitting relations specify the speeds s − and s + of the DSW edges and the associated wave parameters k − and a + in terms of the Riemann data (u − , u + ). Since we focus on positive DSW polarity and orientation, the leftmost trailing edge, propagating with constant speed s − , is associated with the vanishing amplitude harmonic wave, while the rightmost leading edge, propagating with speed s + , is associated with the solitary wave. Thus, s + > s − . In what follows, the Riemann step data (u − , u + ) of (1.7) are fixed in the regime of interest (6.1), where we assume the existence of the integral 2-wave curve of the modulation system connecting the harmonic and soliton edges. The DSW fitting relations are summarized in three steps (see [12,16] for additional details): (i) The DSW harmonic edge is the modulation characteristic x = s − t with speed s − = s − (u − , u + ) determined by the linear group velocity for the edge mean u = u − and wavenumber k = k − : Given u − and u + , the harmonic edge wavenumber k − is determined by solving the initial value problem , k(u + ) = 0, (6.4) for K(u, u + ) ≡ k(u), and setting k − = K(u − , u + ).
(ii) The DSW soliton edge is the characteristic x = s + t whose speed s + = s + (u − , u + ) coincides with the solitary wave velocity for the edge mean u = u + and amplitude a = a + s + = c(a + , u + ), (6.5) where c(a, u) is the velocity of a solitary wave with amplitude a on the background u. For the description of the soliton edge, it is instructive to make the change of variables (u, a) → (u,k), where the conjugate wavenumberk is defined implicitly by c(a, u) =ω 0 (k, u) k , (6.6) andω 0 (k, u) is the conjugate dispersion relation defined in terms of the linear dispersion relation according toω 0 (k, u) = −iω 0 (ik, u). Given u − , u + , the valuek + at the soliton edge of the DSW is determined by solving the initial value problem ,k(u − ) = 0 (6.7) forK(u − , u) ≡k(u), and settingk + =K(u − , u + ). The soliton edge speed s + and the solitary wave amplitude a + are then found from (6.5), (6.6), withk =k + , u = u + . Note that we have used a slightly modified presentation of DSW fitting as compared to Refs. [12,16] in order to explicitly highlight the dependence of the edge velocities s ± and edge wave parameters k − , a + on the initial Riemann data u − , u + .
(iii) Admissibility conditions: The DSW fitting relations are subject to the existence of a 2-wave curve of the Whitham modulation equations continuously connecting (u, k, a) = (u − , k − , 0) to (u, k, a) = (u + , 0, a + ). The following admissibility conditions are necessary for the existence of a 2-wave solution [12,26] Causality:

DSW interior modulation
It was shown in [10] that the DSW fitting procedure can be extended into the DSW interior x > s − t adjacent to the harmonic edge, where the wave amplitude is small. This is done by exploiting the overlap domain for the asymptotic applicability of weakly nonlinear modulated (Stokes) waves in the Whitham modulation equations (describing slow modulations of arbitrary amplitude waves) and the nonlinear Schrödinger equation (describing slow modulations of finite but small amplitude waves). A brief outline of the results of [10] relevant to this paper is presented below.
The weakly nonlinear regime of the DSW |u(x, t) − u − | u − is well described by a slowly modulated Stokes wave where k − = K(u − , u + ) is determined by DSW fitting, A(x, t) ∈ C and B(x, t) ∈ R are slowly varying fields in comparison to the fast oscillations of the carrier wave. Note that B(x, t) describes the induced mean flow. The crest to trough amplitude a(x, t) and the wavenumber k(x, t) = k − + v(x, t) of the Stokes wave are related to the complex envelope A(x, t) by The dynamics of A(x, t) are governed by the nonlinear Schrödinger (NLS) equation where β − = β(k − , u − ) = 1 2 ∂ kk ω 0 (k − , u − ) and γ − = γ(k − , u − ) are determined as part of the standard multiple scales derivation of the NLS equation (see for instance [5]). Note that the slow variation of the background u(x, t) is solely due to the induced mean flow and is constrained to the variation of the amplitude according to is also determined as part of the NLS derivation. The NLS equation describes a stable modulation iff β − γ − < 0 (defocusing regime). Also, the derivation requires that β − does not vanish, which is the case in the convex modulation regime considered here. For the description of a DSW's structure, the NLS equation (6.12) is considered with boundary conditions specified at the harmonic edge a(x, t) = 0, k(x, t) = k − . (6.14) As a result, the universal DSW modulation near the trailing edge is given by the following special selfsimilar rarefaction wave solution of (6.12) The formulae (6.15) represent a first-order approximation of the DSW's structure near its harmonic edge. A more accurate description of the DSW modulation, which asymptotically extends deeper into the DSW's interior, is achieved by including higher-order terms in the derivation of the NLS equation. The resulting HNLS equation (see Appendix) subject to the boundary conditions (6.14) yields a rarefaction wave solution for a and k that includes terms O( x t − s − ) 2 in addition to the linear terms of eq. (6.15).

BBM DSW fitting
We now apply the DSW fitting relations (6.3)-(6.7) and the universal modulation formulae (6.15) (and its HNLS improvement) to the description of BBM DSWs generated within the region of Riemann data u − , u + defined by the admissibility conditions (6.8), (6.9). The BBM equation satisfies the prerequisites for DSW fitting [12,16]: periodic travelling waves parameterised by three independent integrals of motion exist (with harmonic wave and solitary wave limits) and BBM has at least two conservation laws (1.4) and (5.11). The speed-amplitude relation for solitary waves is given by (5.2). We initially consider the Riemann data (1.6) with u − > u + > 0 (other configurations will be considered later).
We now verify the DSW fitting admissibility conditions (6.8), (6.9). The first two causality conditions (6.8) for BBM assume the form s − < u − , s + > u + and are readily verified providedk + < 1, which is true so long as u − > u + > 0. Next, a calculation yields While the first convexity condition in (6.9) is always satisfied, the second condition requires k − = √ 3, which, after substitution in (6.18), gives (u − /u + ) = 4e 3/2 . The critical line u + = µ u − , µ = e −3/2 /4 ≈ 0.056 (6.27) corresponds to the zero dispersion point k − = √ 3 and delimits the applicable region for the DSW fitting method 0 < µu − < u + . (6.28) Another calculation yields (6.29) Equation (6.23) implies that 0 <k + < 1 within the region (6.28). Then the third and fourth convexity conditions in (6.9) are readily verified to hold true for all u − , u + in (6.28). It follows that within the admissibility region (6.28), k − < √ 3 so the dispersion sign is Consequently, the DSW orientation and polarity are d = p = 1, consistent with our original assumption s + > s − , which ensures the third causality condition (6.8) within the region (6.28) for µ defined in (6.27). Figure 20a shows excellent agreement between the analytical expressions (6.19) and (6.24) and the velocities of the DSW edges extracted from numerical simulations. Additionally, Fig. 20b   The interior BBM DSW modulation in the vicinity of the harmonic edge is given by (6.15). It is determined by the coefficients β(k, u) and γ(k, u) of the NLS equation (6.12) for slowly varying modulations of the weakly nonlinear periodic traveling wave of the BBM equation with k = k − , u = u − . The derivation of the NLS equation is standard and we only present here the expressions for β and γ .
Within the admissible region (6.28), k − < √ 3 and so β(k − , u − )γ(k − , u − ) < 0, thus ensuring modulational stability of the DSW's harmonic edge. Additionally, the function b(k, u) determining the variation of the induced mean flow (6.13) near x = s − t is given by . (6.32) A comparison between the modulation solution (6.15), (6.13), and the DSW structure extracted from numerical simulations is displayed in Fig. 21. The figure also displays the modulation solution obtained using the HNLS approximation, see Appendix B and [10] for details. The HNLS modulation solution provides the best agreement with the numerical DSW structure near the harmonic edge and deeper into the DSW interior. Near the soliton edge, there is a departure from the HNLS description, which is to be expected because the (H)NLS descriptions are based upon modulations of a weakly nonlinear Stokes wave, which does not describe the DSW in the strongly nonlinear, or solitary wavetrain regime.  The wavenumber and mean value predictions in Figures 21b, 21c visually exhibit better agreement with the numerically extracted DSW parameters than the DSW amplitude in 21a. There is some deviation in the numerically extracted oscillation envelope amplitude depicted in Fig. 21a near the trailing edge from predicted linear behavior a ∝ x t − s − in (6.15). This is not surprising given that our analysis in Section 3 demonstrates that the transition width for smoothed step initial data significantly alters the DSW's oscillatory structure in the neighborhood of the DSW harmonic edge. In fact, a fundamental discrepancy between the DSW modulation solution and the numerical oscillation behavior is known and even occurs for the KdV equation [22]. In contrast to the well-defined DSW soliton edge, the DSW harmonic edge discrepancy from leading order modulation theory leads to ad hoc approaches for comparing the harmonic edge velocity and wavenumber with numerical simulations. Nevertheless, modulation theory provides a reasonable prediction of the DSW's amplitude, wavenumber, and mean structure near the DSW harmonic edge. Figure 22 shows that the solution of the dispersive Riemann problem is no longer a single-phase DSW when u + < µu − . A beating pattern develops close to the DSW small-amplitude edge that is clearly depicted in the examples for which u + ≤ −0.3. This phenomenon has been previously identified [35] as DSW implosion in which the group velocity of linear waves at the harmonic edge experiences a minimum, i.e. at a zero dispersion point. In this section, we investigate the solution for u + < µu − in greater detail. In the numerical examples presented in this section, ξ = 10 unless otherwise stated. The variation of a typical solution when u + = −0.5 in Fig. 22 shows that 3 distinct modulation regions exist: a finite amplitude, stable, single-phase wave in region I, a beating pattern in region II and a small amplitude wave in region III; the boundaries between regions I and II denoted x = x b , and regions II and III denoted x = x a , separate the spatial regions where the wave is single-phase and two-phase.

Region I: stable single-phase wave
We first consider the modulation of the wave in region I. The variations of the crest to trough amplitude a(x, t), the wavenumber k(x, t), and the mean u(x, t) obtained numerically are displayed in Fig. 23 (where u 1 = u in region I, as detailed in Sec. 6.3.2). Region I can be subdivided into two subregions: a region denoted Ia where a(x, t), k(x, t), andū(x, t) are approximately constant, apart from localized, smallwidth modulations, and a region denoted Ib where a(x, t), k(x, t), and u(x, t) vary monotonically. The modulation in region Ib is similar to the DSW modulation investigated in Secs. 6.1 and 6.2 where the wavenumber k(x, t) monotonically decreases until it reaches the solitary wave limit k = 0, while the   amplitude monotonically increases toward the soliton edge. However, the DSW fitting of Sec. 6.1 does not apply to this region since the admissibility condition (6.28) is not fulfilled, and the parameters of the soliton edge cannot be obtained using this method. In fact, the wave structure in region Ib can be identified as a partial DSW, where the nonlinear modulated wave does not reach the small-amplitude, harmonic limit. Instead, it terminates at some nonzero amplitude a 0 , meanū 0 , and wavenumber κ 0 that are matched to the right edge of region Ia. These types of structures occur in initial-boundary value problems for nonlinear dispersive PDEs (see, e.g. [38]) and, in particular, are realized in problems involving resonant or transcritical shallow-water flows past localized bathymetry [24,44,13]. We do not present a quantitative description of partial BBM DSWs here as this would require the development of a full nonlinear modulation theory for the BBM equation, a significant undertaking that deserves a dedicated study.
The specific modulation of the wave in region Ia is a new feature that only develops when σ 3 u − < u + < µu − . The localized depletions of the amplitude a(x, t) and the corresponding augmentations of the wavenumber k(x, t) propagate on the "background" (a, k) = (a 0 , κ 0 ) at constant velocity; we have, for instance in the example of Fig. 23 where u + = −0.5, (a 0 , κ 0 ) ≈ (2.19, 1.19). In the following, we investigate the modulation of the field u(x, t) using a weakly nonlinear NLS description that was introduced in the previous sections. The nonlinear wave can be approximated in the weakly nonlinear regime by the Stokes wave (cf. for instance (6.10)) u(x, t) ∼ u(x, t) + A(x, t)e i(κ0x−ω0(κ0,u−)t) + c.c. , (6.33) where u(x, t) ∼ u − is a slowly varying field. The crest to trough amplitude a and the wavenumber k are given by a(x, t) = 4|A(x, t)|, k(x, t) = κ 0 + ∂ x arg A. (β(κ 0 , u − )γ(κ 0 , u − ) < 0): We first consider the nonlinear plane wave solution a = a 0 = const, k = κ 0 = const. As a result, the weakly nonlinear wave in region Ia is approximately showing that weakly nonlinear interaction shifts the linear frequency ω 0 (κ 0 , u − ) by −∆ω 0 . Figure 24a displays the frequency shift computed numerically for the example u + = −0.5, and shows reasonable agreement with the NLS equation description (6.35). Equation (6.38) also yields u = 0.77 whereas the average extracted numerically is u ≈ 0.67 as shown in Fig. 23. A better analytical approximation of ∆ω 0 and u necessitates higher order terms in the NLS description.  This phase-shift can be observed in the numerical solution but differs from the analytical formula (6.40). This indicates that the wave in region I is governed by the NLS equation only to a first approximation, and the description of traveling waves of the envelope A(x, t) necessitates a higher order description such as the HNLS equation (B.2). We leave this comparison as a subject for future work. The emergence of envelope solitary waves is reminiscent of the shedding of solitary waves investigated in Sec. 5 because solitary waves are only generated for a sufficiently large jump |u − − u + |. Figure 22 shows that the implosion of the DSW occurs in regions II and III, which broadly overlap with the two-phase region x ∈ (s min t, 0) identified in the linear analysis of Sec. 3. In this section, we focus on region II. We show that u(x, t) displays a beating pattern that can be approximated by a superposition of two single-phase, weakly nonlinear waves

Region II: stable two-phase wave
were u 1 and u 2 correspond, respectively, to slow and fast oscillating waves. In order to determine u 1 and u 2 from our numerical simulations, we compute the spatial Fourier transform of u II (x, t), which is the restriction of u(x, t) to region II. The example u + = −0.5 is shown in Fig. 25 depicting multiple maxima of the Fourier transform. The two largest spectral contributions correspond to the dominant wavenumber κ 1 of the slowly oscillating wave u 1 and the dominant wavenumber κ 2 > κ 1 of the fast oscillating wave u 2 . Other maxima correspond to harmonics that result from weakly nonlinear interactions such as second harmonic generation at 2κ 1 and 2κ 2 , and two wave interaction κ 2 ± κ 1 . We then identify the small-wavenumber content of the Fourier spectrum (k < k c ) as the Fourier transform of u 1 (x, t) and the large-wavenumber content as the Fourier transform of u 2 (x, t), where k c is the minimum identified for κ 1 < k < κ 2 (k c ≈ 2.8 in the example depicted in Fig. 25). The variation of the fields u 1 and u 2 is obtained by computing the inverse Fourier transforms of the band limited functions The division into small and large wavenumber bands of the Fourier spectrum was examined in Sec. 3 for linear waves. There, k c = √ 3 separated the small and large wavenumber bands. In the weakly nonlinear regime, k c = √ 3 because the second harmonic 2κ 1 can be larger than √ 3 (2κ 1 2.6 in the example presented in Fig. 25). The variation of u 1 (x, t) and u 2 (x, t) extracted from the numerical solution for u − = −0.5 and their modulation parameters in region II are depicted in Fig. 23. Note that the numerical method presented here to extract the two single-phase waves only yields an approximation to u 1 and u 2 . For instance, the higher harmonic κ 2 − κ 1 contributes to the Fourier spectrum of u 2 in the example considered in Fig. 25, even though it corresponds to the nonlinear interaction between u 1 and u 2 . In practice, the quasi-monochromatic aspect of the extracted waves u 1 and u 2 (cf. Fig. 23) indicates that the numerical method presented above yields a good approximation. We now make the ansatz that the two modulated wavetrains u 1 and u 2 can be modeled by two Stokes waves u 1 (x, t) = A 1 (x, t)e i(κ1x−ω0(κ1,u−)t) + c.c. + u(x, t), u 2 (x, t) = A 2 (x, t)e i(κ2x−ω0(κ2,u−)t) + c.c., (6.43) where A 1 (x, t), A 2 (x, t) are complex valued and u(x, t) is real valued, all of which are slowly varying. The dominant wavenumbers of the Stoke waves satisfying κ 1 < √ 3 and κ 2 > √ 3 are extracted numerically using the procedure already described. Note that by the definition of u 1 and u 2 , the slow modulation of the background u(x, t) is necessarily incorporated in the variation of u 1 (x, t) so that u 2 (x, t) has zero mean. The crest to trough amplitudes a j and wavenumbers k j are a j (x, t) = 4|A j (x, t)|, k j (x, t) = κ j + ∂ x arg A j , j = 1, 2. (6.44) The modulation equations for the complex envelopes A 1 and A 2 are two coupled NLS equations obtained with a standard multiple-scale computation. The coefficients β(κ, u − ) and γ(κ, u − ) are given in (6.31) and b(κ, u − ) in (6.32). The nonlinear interaction between the two components A 1 and A 2 is characterized by the coefficient ν(κ, λ, u − ) = κ(1 + λ 2 )(9 + 6κ 2 + κ 4 + 21λ 2 + 13κ 2 λ 2 + κ 4 λ 2 + 9λ 4 + 3κ 2 λ 4 + λ 6 ) (1 + κ 2 )λ 2 (3 + λ 2 )(3 + κ 2 − κλ + λ 2 )(3 + κ 2 + κλ + λ 2 )u − > 0. (6.46) The modulation description (6.43), (6.45) is the weakly nonlinear extension of the modulated two-phase linear wave solution (3.18), which applies when the two wavenumbers k 1 and k 2 are distinct. As indicated in Fig. 23, Region II can be divided into two subregions: region IIa where the wave parameters are approximately constant, and region IIb where the variations of a j (x, t), j = 1, 2 are approximately monotone.
We first focus on the modulation in region IIb where the amplitudes a 1 (x, t) and a 2 (x, t) are, respectively increasing and decreasing functions of x. Ultimately a 2 (x, t) = 0 at the boundary with region I (x = x b ) and the amplitude of the slow varying wave a 1 (x, t) matches with a(x, t) = a 0 in region Ia. We thus extend the definition of u 1 and u 2 to region I by setting: The modulation close to x = x b can thus be uniformly described, across the boundary between the region IIb and Ia, by the simplified coupled equations The approximation of (6.45) invoked here is |A 2 | |A 1 |, as suggested by Fig. 23. The modulation of A 1 is independent of A 2 and is governed by the, modulationally stable, defocusing NLS equation as described in Sec. 6.3.1. The modulation of A 2 is a linear Schrödinger equation with potential −ν(κ 2 , κ 1 , u − )|A 1 (x, t)| 2 . Consequently, although the wavenumber κ 2 is larger than √ 3, the modulation dynamics of u 2 remain stable in the vicinity of x = x b . This highlights the stabilizing effect of the slowly-oscillating component u 1 (x, t) on the fast-oscillating component u 2 (x, t) in region IIb. We do not carry out a detailed analysis of the modulations A j in region IIb, leaving that for a separate, dedicated study.
We now consider region IIa for which the modulation variables are approximately constant where we identify the modulation wavenumber k j with the constant, dominant wavenumber κ j extracted from the numerical simulation so that ∂ x arg A j = 0. The values of the parameters a j are also extracted from the numerical solution. In Fig. 23, we have (a 1 , k 1 ) ≈ (1.85, 1.31) and (a 2 , k 2 ) ≈ (0.77, 4.75). Substituting (6.49) into the coupled NLS equations (6.45), we obtain the modulation of the complex envelopes in the form of two plane waves As a result, the two-phase weakly nonlinear wave in region IIa for the BBM equation is approximately showing that weakly nonlinear interactions between the two wave modes shift the linear frequencies ω 0 (κ 1 , u − ) and ω 0 (κ 2 , u − ) by −∆ω 1 and −∆ω 2 , respectively. Figure 24b displays good agreement between the weakly nonlinear frequency predictions (6.38) and the temporal Fourier spectrum of the BBM numerical simulation at x = −300 for the case u + = −0.5. Equation (6.51) also yields u = 0.86 whereas the average extracted numerically in Fig. 23 is u ≈ 0.72. A better approximation of ∆ω 1 , ∆ω 2 and u could be achieved with higher order terms in the NLS description. It is important to point out that, in the absence of two-phase coupling, the modulation A 2 satisfies the scalar NLS equation (6.35) with κ 0 → κ 2 . Since κ 2 > √ 3, β(κ 2 , u − )γ(κ 2 , u − ) > 0 so that the equation is of focusing type and the plane wave solution (6.36) (with a 0 → a 2 ) for A 2 is linearly unstable to infinitesimal, long wavelength perturbations, also known as modulationally unstable. However, the nonlinear coupling between the modulations in the coupled NLS equations (6.45) can stabilize the plane wave solution (6.50). For the specific example shown in region IIa of Fig. 23, we have verified that this is indeed the case with a calculation described in Appendix C.

Region III: linear wave
Since the carrier wavenumber of the field u(x, t) is larger than √ 3 in region III-cf. Fig. 23 where (u 1 , u 2 ) = (u − , u − u − ) in region III-u represents a fast oscillating wave according to the criterion introduced in Sec. 6.3.2, and we set (6.52) Contrary to the neighborhood of the boundary x = x b , the fields u 1 and u 2 are rapidly changing at the boundary x = x a . The amplitude of u(x, t) in region III is small and its modulation can be described by the linear theory introduced in Sec. 3. The bottom row of Fig. 23 shows that the modulation of the wavenumber k(x, t) obtained numerically lies on the stationary point k 2 (x/t) from linear theory given by (3.15).

Summary of DSW implosion
Here we summarize our description of DSW implosion for dispersive Riemann problems with σ 3 u − < u + < µu − in regions c and f of Figs. 1 and 2, respectively (σ 3 is defined in the next subsection). The entire coherent wave structure arises due to compressive initial data that cannot be resolved by an admissible DSW alone. Instead, both short and long waves are generated that represent a nonlinear generalization of the two-phase, linear wave modulation investigated in Sec. 3. From right to left, the solution is described by modulated long waves in a partial DSW (region I) that transitions to region II where long and short waves coexist, significantly influencing one another. This two-phase modulation abruptly terminates at a short wavelength, linear wavepacket (region III) propagating on the background u − , hence the entire structure transitions from modulated, long nonlinear waves to modulated, short linear waves.
Similar to region I identified in linear theory, the nonlinear wave in region I is a modulated one-phase wave with a small wavenumber (k < √ 3). The one-phase wave in region Ib is a partial DSW, a carryover from dispersive Riemann problems with µu − < u + < u − in which DSWs are admissible. Its wave parameters a, k, and u are strictly monotone functions of x in this region. Adjacent to the partial DSW is a stable nonlinear plane wave with multiple localized, depression modulations of the wave envelope (dark envelope solitary waves) that stably propagate in region Ia.
The wavetrain in region Ia matches smoothly to a nonlinear beating pattern in region II at x = x b . In a first approximation, u(x, t) is the superposition of two modulated, weakly nonlinear, monochromatic (Stokes) waves u 1 (x, t) and u 2 (x, t), each with distinct dominant wavenumbers κ 1 < √ 3 and κ 2 > √ 3. This is a generalization of the superposition of two modulated, linear wavetrains with distinct dominant wavenumbers (see region II described in Sec. 3). The weakly nonlinear components u 1 and u 2 are approximately governed by coupled NLS equations, as demonstrated by their nonlinear frequency shifts. Although a weakly nonlinear wave of dominant wavenumber larger than √ 3 is modulationally unstable, the component u 2 (x, t) is shown to be stable when nonlinearly coupled to u 1 (x, t).
Region III displays a principal difference from linear theory. In the linear theory, the two wavenumbers k 1 and k 2 coalesce at √ 3 when x = −u − /8t, and the corresponding linear wave u(x, t) continuously changes from a beating two-wave structure to an Airy wave. In the dispersive Riemann problem considered here, the modulation changes very rapidly, in fact, discontinuously on the modulation length scale at x = x a . The adjacent wave in region III is a linear, one-phase wave with wavenumber k 2 > √ 3. The discontinuity at x = x a displayed in Fig. 23 resembles a so-called Whitham shock wave structure that was recently identified and investigated for several nonlinear dispersive PDEs in [46]. Higher order dispersion was shown to be essential for the generation of Whitham shocks. The coupled NLS equations (6.45) we have used to effectively describe region II do not incorporate higher order dispersion.
We have focused in this subsection on DSW implosion resulting from broad initial data ξ = 10. Therefore, this nonclassical wave pattern emerges during the course of nonlinear evolution and does not require the imposition of a small scale feature in the initial data. In fact, DSW implosion for narrow data ξ 1 results in essentially the same nonlinear wave coherent structure. See Fig. 2f.

Incoherent solitary wavetrain
As the initial smoothed step height is increased through the DSW implosion regime that is described in the previous subsection, a new feature emerges when −u − < u + < σ 3 (ξ)u − (cf. Fig. 1d): the emission of solitary waves. This effect is thoroughly described in section 5 for expansive initial data. We now examine solitary wave emission accompanying compressive data. This subsection is somewhat speculative and brief because we base our description of the dynamics solely upon numerical simulation. Figure 26 depicts the numerical solution for narrow 26a and broad 26b, smoothed step initial data. In the vicinity of u + = −0.7, we observe a depression solitary wave emitted to the left edge of the two-phase oscillatory region II (recall Fig. 25) by t = 1000. Other simulations, not shown, exhibit solitary wave emission at later times for u + ≈ −0.65. As the right edge u + is reduced below −0.7, more solitary waves, including elevation solitary waves propagating to the right, are emitted. These simulations depict a gradual transition as u + is decreased from DSW implosion to the emission of solitary waves of an increasingly incoherent character. As u + approaches −u − , the solitary waves exhibit different amplitudes and are emitted intermittently as shown in the simulation animations that are included as supplementary material to this manuscript. While the simulations for different smoothed step width ξ shown in Figs. 26a and 26b exhibit differences, their qualitative character is the same. This stands in contrast to the case of solitary wave shedding for expansive data where we showed in section 5 that the initial smoothed step width must be sufficiently narrow, i.e. the initial data must include a sufficiently small scale feature. We remark that qualitatively similar dynamics were observed numerically for broad, localized initial and periodic data in [21].
These pseudospectral simulations were highly resolved. The numerical method is different from that utilized in other sections of the paper (see Appendix A.1) and is described in [11]. The grid spacing ∆x = 0.00343 and timestep ∆t = 0.002 ensured that the conserved energy E(t) in eq. (5.15) remains numerically within 10 −10 of its predicted value E(t) = 1 3 (u 3 + − u 3 − )t + E(0) (eq. (5.14)) for the simulated times t ∈ [0, 1000] and domain x ∈ [−900, 900]. In all cases, the relative magnitude of the Fourier coefficients |û(k, t)/û(0, t)| decays to 10 −15 for sufficiently large k. Furthermore, the large spatial domain ensures that the boundary deviations from u ± remain below 10 −15 in magnitude for the simulated times.
One important finding from the simulations in Figure 26 is that even a slight asymmetry in the initial smoothed step (u + = −0.99u − ) leads to a significant departure from the symmetric dispersive Lax shock that is generated when u + = −u − (recall section 4.1). The animations for this case that are included as supplementary material are particularly revealing in this regard. Thus, the dispersive Lax shock lacks robustness to small changes in the boundary conditions.

Conclusions and Outlook
In this paper, we have developed a classification of solutions to the dispersive Riemann problem for the BBM equation. The BBM equation (1.1) as a mathematical model for unidirectional, weakly nonlinear dispersive shallow water wave propagation is asymptotically equivalent to the celebrated KdV equation while providing more satisfactory short-wave/high-frequency behavior in the sense that the linear dispersion relation is bounded for the BBM equation, but unbounded for the KdV equation. However, the bounded dispersion relation is nonconvex, a property that gives rise to a number of intriguing features that are markedly different from those found in the KdV equation. Some of these features exemplify phenomena previously observed in other nonlinear dispersive equations, but some are new, providing the motivation for the study of the reduced BBM equation (1.4) as a distinct dispersive regularization of the Hopf equation.
The main feature of the BBM equation (1.4) that distinguishes it from the KdV equation (1.3) is the non-evolutionary dispersive term −u xxt . This term introduces an intrinsic nonlocality length scale and results in a qualitative and quantitative dependence of solutions on the comparison between and the typical spatial scale ξ associated with the initial condition. The nonlocality length coincides with the BBM coherence length (the unit amplitude solitary wave width) defined by the balance of nonlinearity and dispersion. In the normalization of (1.4), = 1. The nonlocal nature of the BBM equation (1.4) imparts a bounded linear dispersion relation ω = ω 0 (k, u) that is nonconvex, exhibiting zero-dispersion points ∂ kk ω 0 (k, u) = 0 for particular values of the wavenumber and background. Generally, we find that the complex interplay between nonlocality, nonlinearity and nonconvex dispersion in the BBM equation gives rise to remarkably rich dynamics for the smoothed dispersive Riemann problem with initial data (1.7) in the form of a monotone transition between two constant states u − at x → −∞ and u + at x → +∞. We identify regions in the plane of boundary states (u − , u + ) where BBM solutions exhibit qualitatively different behaviors. Additionally, the evolution of data within these regions often depends significantly upon the smoothing of the initial data characterized by its transition width ξ, thus introducing a further distinction between various regions. To capture the smoothing-dependent features, we construct two separate dispersive Riemann problem classifications: for ξ 1 and for ξ 1 shown in Figures 1 and 2, respectively. A detailed guide to the wave patterns presented in these figures is contained in Section 2, which summarizes our findings.
Emergent wave phenomena for the BBM dispersive Riemann problem can be roughly split into two categories: classical and nonclassical. Classical phenomena include dispersive shock waves (DSWs) and rarefaction waves, which are also observed in convex, KdV-type dispersive hydrodynamics, hence termed classical. The nonclassical features are due to nonconvex dispersion and include the generation of twophase linear wavetrains, expansion shocks, solitary wave shedding, dispersive Lax shocks (in a conservative dispersive equation!), DSW implosion and the generation of incoherent solitary wavetrains. Some of these striking nonclassical features are present only in the evolution of narrow initial steps, ξ 1. We stress that all the observed nonclassical phenomena involve short wave features that are formally not applicable in the traditional asymptotic regime of the original BBM equation (1.1) in which u xxt ∼ −u xxx .
The description of a diverse range of linear and nonlinear wave patterns emerging from the BBM dispersive Riemann problem has required an equally diverse set of approaches-analytical and numerical-to their analysis and description. We analyze linear wavetrains generated near the zero dispersion point using the stationary phase method and matched asymptotics, while solitary wave shedding is tackled by a combination of short time analysis and energy arguments. DSWs in the convex propagation regime, defined by strict hyperbolicity and genuine nonlinearity of the associated Whitham modulation system, are analyzed using the extended DSW fitting method while DSW implosion dynamics are shown to be approximately described by the coupled NLS equation derived from BBM via multiple-scale expansions for weakly nonlinear, two-phase waves. All analytical results are accompanied by detailed numerical simulations. The numerical methods and approaches are presented in the Appendix.
Another interpretation of these results incorporates the distinction between long and short waves, evident in the aforementioned difference between the linear dispersion relations for KdV (1.3) and BBM (1.4) in which the phase and group velocities, although asymptotically equivalent for long waves, are bounded for BBM and unbounded for KdV in the short wave regime. While linear theory predicts that short waves in BBM are feeble, in the words of the equation's namesake [4], nonlinear wave propagation does not suffer from this restriction [37]. The prevalence of nonlinear, short-long wave interaction in the dispersive Riemann problem for BBM is striking. When a short scale is encoded in the initial condition, new features arise in the solution, including expansion shocks, the emission of solitary waves, and 2-phase linear wavepackets. Short wave effects also arise during the course of nonlinear evolution, independent of the spatial scale of the initial data. Dispersive shock wave implosion, the dispersive Lax shock, and incoherent solitary wavetrains are notable examples. In fact, the self-similar dispersive Lax shock u(x, t) = g(xt) exhibits short waves that are unbounded in wavenumber as t → ∞. While a slight asymmetry in the boundary conditions arrests the unbounded wavenumber production in the dispersive Lax shock, the resulting intermittent solitary wave emission suggests a conservative counterpart to the chaos observed in numerical simulations of a damped-driven BBM equation [42]. None of these short wave effects occur in the KdV equation.
Some of the results obtained in this paper carry over to other dispersive hydrodynamic systems. The observed prominent dependence of the long-time dynamics on the initial transition width points to a leading order effect that has mostly been ignored in previous DSW studies. The analysis of twophase linear wavetrains near a zero dispersion point is general, and the coupled NLS description of DSW implosion can be applied to the conduit and magma equations where this effect is known to occur. Expansion shocks and solitary wave shedding have been observed in the regularized Boussinesq shallow water system [17] albeit in a nonphysical, short-wave regime. Still, it would be fair to say that many of the striking nonconvex wave regimes generated in the BBM dispersive Riemann problem await their realization in other systems.
Our work poses a number of interesting questions that will hopefully be addressed in future research. These include a better understanding of solitary wave shedding and the intermittency associated with the incoherent solitary wave regime, stability analysis of dispersive Lax shocks, and a full analytical description of the nonlinear two-phase structures and envelope solitary waves associated with DSW implosion. and ϕ 2 should only hold in the (x, t)-region where k 1 (x, t) √ 3 k 2 (x, t). In practice we see in Fig. 3 that the two linear waves are still distinct even if k 1 and k 2 are close to √ 3. Note that we also obtain numerically k 1 ∼ k 2 ∼ √ 3 for x s min t even if no linear waves can propagate in the region. This is a numerical artifact and one can check that ϕ 1 and ϕ 2 are out of phase in this region such that ϕ = ϕ 1 + ϕ 2 = 0 as expected.

B HNLS description B.1 HNLS equation
The NLS equation displays both a dispersive effect term A xx ∝ ε 3 and a weakly nonlinear term |A| 2 A ∝ ε 3 , where ε is the small parameter describing the slow spatio-temporal dependence and small amplitude of the envelope A(x, t), cf. [5]. The next order of the envelope dynamics of the weakly nonlinear wave: u(x, t) ∼ u 0 + A(x, t)e i(κ0x−ω0(κ0,u0)t) + c.c. + B(x, t), (B.1) is commonly called higher order NLS (HNLS) equation, includes the terms: A xxx ∝ ε 4 , |A| 2 A x ∝ ε 4 and A 2 A * x ∝ ε 4 (cf. our paper and references therein). The HNLS equation for the BBM equation reads: iA t + i∂ k ω 0 A x + βA xx + γ|A| 2 A + iδA xxx + iλ|A| 2 A x + iνA 2 A * x = 0, where the coefficients β, γ and b are given in (6.31), and a and v are defined in (6.11). δ, λ, ν and b 2 are function of the parameters κ 0 , u 0 given by:

B.2 DSW modulation
The universal DSW modulation near the trailing edge is given by the following special self-similar rarefaction wave solution of (B.2), cf. [10]: (B.4) Figure 28: Linear dispersion relation v(q) for perturbations to the plane wave solution (6.50) corresponding to Fig. 23. The exact roots of (C.1) (solid) and their approximation (C.2) (dashed) are shown. All are purely real, implying modulational stability of region IIa.

C Stability of nonlinear two-phase solution
In [23], it was shown that the stability of the plane wave solution (6.50) to the coupled NLS equations (6.45) boils down to determining the roots v of the quartic polynomial where we have introduced the notation v gj = ∂ k ω 0 (κ j , u − ), β j ≡ β(κ j , u − ), γ j ≡ γ(κ j , u − ), and ν ij ≡ ν(κ i , κ j , u − ) for brevity. A root v of the quartic (C.1) is the phase velocity and q is the wavenumber, both of the infinitesimal plane wave perturbation ∝ exp(iq(x − vt)) of the nonlinear solution (6.50). Consequently, if all four roots of (C.1) are real for every q ∈ R, then the weakly nonlinear two-phase wavetrain (6.50) is modulationally stable. In general, there are no simple, explicit expressions for the roots of (C.1). However, a simplification of the general calculations in [23] can be made when the two group velocities v g1 , v g2 are close so that (C.1) can be approximated by taking (v − v gj ) 2 → (v −v g ) 2 wherev g = 1 2 (v g1 + v g2 ) is the average of the group velocities. Then (v −v g ) 2 in the approximation of (C.1) satisfies a quadratic equation. Introducing the difference of the group velocities ∆v g = 1 2 (v g2 − v g1 ) and assuming |∆v g /v g | 1, the four roots of (C.1) are approximately 2(v −v g ± ∆v g ) 2 = −Q 1 (q) − Q 2 (q) ± (Q 2 (q) − Q 1 (q)) 2 + 4R. (C.2) Then the necessary and sufficient condition for stability is the real, non-negativity of −Q 1 (q) − Q 2 (q) ± (Q 2 (q) − Q 1 (q)) 2 + 4R. Since Q 1 (q) < 0 and Q 2 (q) < 0, it suffices to consider only the long wavelength case q = 0. For the simulation shown in Fig. 23, we calculate − Q 1 (0) − Q 2 (0) ± (Q 2 (0) − Q 1 (0)) 2 + 4R ≈ 0.0111 ± 0.0109 > 0, (C. 3) hence the two-phase plane wave solution (6.50) is stable. A plot of the linear dispersion relation v = v(q) for perturbations to (6.50) as roots of (C.1) and their approximation (C.2) for the extracted simulation parameters in Fig. 23 is shown in Fig. 28.