Hydrodynamics of a Discrete Conservation Law
Patrick Sprenger Christopher Chong Emmanuel Okyere Michael Herrmann P. G. Kevrekidis Mark A. Hoefer
 



#1CC: #1
#1PS: #1
#1(MH: #1)
#1#1 







 


 


 





 
 








2cm






























e
i
ε

sech
Im
φ
∂


∇




ϵε
  d
𝐪
𝐅


    łℓdnℱ𝒱𝒮ℰ𝒲ℒcnsnd̂δ̂.06 lemmaLemma
    
    
    
    Department of Applied Mathematics, University of
      California Merced , Merced, California, 95343Department of Mathematics, Bowdoin College, Brunswick, Maine 04011Department of Mathematics, Bowdoin College, Brunswick, Maine 04011Institute of Partial Differential Equations, Technische Universität Braunschweig, 38106 Braunschweig, GermanyDepartment of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, Massachusetts 01003-4515Department of Applied Mathematics, University of
      Colorado, Boulder, Colorado 80309
      The Riemann problem for the discrete conservation law
      2 u̇_n + u^2_n+1 - u^2_n-1 = 0 is classified using
      Whitham modulation theory, a quasi-continuum approximation, and
      numerical simulations.  A surprisingly elaborate set of solutions to
      this simple discrete regularization of the inviscid Burgers'
      equation is obtained.  In addition to discrete analogues of
      well-known dispersive hydrodynamic solutions—rarefaction waves
      (RWs) and dispersive shock waves (DSWs)—additional unsteady
      solution families and finite time blow-up are observed.  Two
      solution types exhibit no known conservative continuum correlates:
      (i) a
      counterpropagating DSW and RW solution separated by a symmetric,
      stationary shock and (ii) an unsteady shock emitting two
      counter-propagating periodic wavetrains with the same frequency
      connected to a partial DSW or a RW.  Another class of solutions
      called traveling DSWs, (iii), consists of a partial DSW connected to
      a traveling wave comprised of a periodic wavetrain with a rapid
      transition to a constant. Portions of solutions (ii) and (iii) are
      interpreted as shock solutions of the Whitham modulation
      equations. 
    Hydrodynamics of a Discrete Conservation Law
        Mark A. Hoefer
        April 29, 2024
    ============================================
    
    
    
    § INTRODUCTION
    
    
    The hydrodynamics of the conservation law (inviscid Burgers' equation)
    
        u_t +  ( u^2  )_x = 0,
    
    with u(x,t) ∈, x,t ∈ is succinctly expressed by solutions of the
    Riemann problem that consists of Eq. (<ref>) for t > 0 subject
    to the initial condition
    
        u(x,0) =
          
            u_-     x ≤ 0 , 
        
            u_+     x > 0,
             x ∈ ,
    
    for u_±∈.  Solutions must be interpreted in a weak sense and
    depend intimately upon the regularization applied.  For the
    viscous regularization in which equation (<ref>) is
    modified to Burgers' equation u_t + (u^2)_x = ν u_xx where
    ν > 0, the weak solution of (<ref>)—either a moving
    discontinuity (shock) or rarefaction wave (RW)—is uniquely
    determined by considering the strong vanishing viscosity limit
    ν→ 0^+ of Burgers' equation <cit.>.  This results in
    the well-known Rankine-Hugoniot jump condition for the speed
    V = (u_- + u_+)/2 and Lax entropy condition u_+ < V < u_- of
    admissible shock solutions.  Regularization by more complex viscous
    terms (nonlinear, higher order, and viscous-dispersive) generally
    results in the same weak solution <cit.>.  An alternative
    dispersive regularization is the Korteweg-de Vries (KdV)
    equation u_t + (u^2)_x + ϵ^2 u_xxx = 0.  In this case, the
    zero dispersion limit ϵ→ 0 converges weakly in the sense
    that it satisfies the KdV-Whitham modulation equations corresponding
    to averaged conservation laws <cit.>.
    Gurevich and Pitaevskii recognized the physical importance of the
    asymptotic approximation obtained by considering small but nonzero
    ϵ and obtained the dispersive shock wave (DSW) solution of
    the KdV equation for (<ref>) when u_- > u_+ as a self-similar
    solution of the KdV-Whitham modulation equations <cit.>.  They
    also observed that when u_- < u_+, the KdV equation with small but
    nonzero ϵ is well-approximated by the same RW obtained by
    viscous regularization.  Dispersive shock waves are unsteady,
    modulated nonlinear wavetrains connecting two distinct levels
    <cit.>.  In contrast to viscous regularization,
    alternative dispersive regularizations can result in drastically
    different Riemann problem solution behavior in the small dispersion
    regime, particularly when higher order <cit.> or
    nonlocal <cit.> dispersive terms are considered.
    The multiscale dynamics of conservative nonlinear wave equations in
    the small dispersion regime are generally referred to as dispersive
    hydrodynamics <cit.>.
    
    In this paper, we study the dispersive hydrodynamics of the
    discrete regularization of Eq. (<ref>)
        / t u_n + 12  ( u_n+1^2 - u_n-1^2
             ) = 0,
    
    by solving the Riemann problem
    
        u_n(0) = 
                     u_-,       n ≤ 0 
        
                     u_+,       n > 0 ,
    
    for (<ref>) where n ∈, t∈, u=u_n(t) ∈.  Equation
    (<ref>) is the simplest centered differencing scheme for the
    hydrodynamic flux.  As recognized in the early days of computational
    fluid dynamics by Von Neumann, and later clarified by Lax,
    differencing schemes like (<ref>) introduce oscillations that
    require mitigation if one wishes to converge strongly to viscously
    regularized solutions that satisfy the conservation laws of fluid
    dynamics <cit.>.  In this paper, we consider
    Eq. (<ref>), subject to (<ref>), in its own right, divorced
    from the aim of approximating solutions of Eq. (<ref>).  The
    semi-discrete equation (<ref>) can be interpreted as the
    dispersive regularization
    
        u_T + i/ϵsin(-i ϵ∂_X)  (u^2
           ) = 0 ,
    
    by introducing T = ϵ t, X = ϵ n into
    Eq. (<ref>), where ϵ > 0 is the lattice spacing.  Here,
    isin(-i∂_x)=12(exp(∂_x)-exp(-∂_x))
    is the pseudo-differential operator for the centered discrete
    derivative and acts on a function f in the variable x via
    
        (isin(-i∂_x)f)(x)=12(f(x+1)-f(x-1)).
    
    Equation (<ref>) is similar to Whitham type evolutionary
    equations
    <cit.>
    except that it is further constrained to be band-limited.  For the
    lattice equation (<ref>) at time t, the support of the
    discrete-space Fourier transform of u_n(t) is [-π,π] due to
    the smallest length scale set by the lattice spacing, whereas the
    Fourier transform of quasi-continuum approximations is 
    not generally
    compactly supported.  As we will demonstrate, this fundamental
    property of lattice equations introduces new hydrodynamic solution
    features that do not appear within certain continuum limits of the
    model.
    
    Equation (<ref>) can be used to formally derive quasi-continuum
    approximations by using Padé approximants of
    i/ϵsin(-iϵ∂_X) for
    0 < ϵ≪ 1.  For example, the (1,3) Padé approximant
    
        i/ϵsin(-i ϵ∂_X) ≈ ( 1- ϵ^2/6∂^2_X )^-1∂_X ,
    
    leads to the Benjamin-Bona-Mahoney (BBM) equation absent the linear
    convective term <cit.>
        U_T + (U^2)_X - ϵ^2/6
          U_XXT = 0.
    
    This quasi-continuum approximation, inspired by the work of Rosenau on
    mass-spring chains <cit.>—see also the discussion
    in <cit.>—is expected to faithfully represent the
    long-wavelength behavior of the lattice model (<ref>). But it is
    no longer band-limited.  In fact, outside of RWs and DSWs, the Riemann
    problem solutions we obtain for the lattice equation (<ref>) bear
    no resemblance to the corresponding Riemann problem solutions of
    (<ref>) obtained in <cit.>.
    Nevertheless, the quasi-continuum model (<ref>) admits exact
    solitary and periodic traveling wave solutions that can be used to
    approximate corresponding solutions of the lattice model (<ref>).
    
    A particular feature of the band-limited lattice equation
    (<ref>), and others with centered differences, is the existence
    of stationary, period two (binary) oscillation solutions
    
        u_2n = α,    u_2n+1 = β,    n ∈,
    
    for any α,β∈.  Equation (<ref>) was studied in
    <cit.> using extensive numerical simulations for certain types
    of odd initial data and binary oscillations were found to play an
    important role.  Allowing for slow spatio-temporal modulations of this
    solution—α = α(X,T), β = β(X,T),
    T = ϵ t, X = ϵ n, 0 < ϵ≪ 1—Turner and
    Rosales <cit.> obtained the modulation equations
    
        α_T + (β^2)_X = 0,   β_T + (α^2)_X = 0 .
    
    Two important implications of the hydrodynamic-type equations
    (<ref>) are: (i) the equations (<ref>) are elliptic whenever
    αβ < 0 and (ii) the existence of discontinuous, shock solutions
    satisfying Rankine-Hugoniot jump conditions.  The ellipticity of
    (<ref>) was shown to provide a route through modulated binary
    oscillations to blow up of solutions of (<ref>).  
    
    
    
    While the topic of DSWs and dispersive hydrodynamics has been more
    extensively explored in the realm of continuum
    media <cit.>, as already implied by the above discussion,
    studies of the discrete realm have the potential to offer new and
    intriguing wave features.  In addition, the motivation for such
    explorations has significantly increased on account of a diverse range
    of corresponding applications. A central topic is the study of
    granular crystals and associated nonlinear
    metamaterials <cit.>, consisting typically of elastically
    interacting bead chains. There, a sequence of experimental efforts in
    simpler <cit.>, as well as in progressively more
    complex media, including dimers <cit.> and the more
    recent setup of hollow elliptic cylinders <cit.> have
    manifested the spontaneous emergence of DSWs under suitable loading
    conditions. However, this has not been the only setting where
    “effectively discrete” DSWs have experimentally emerged. Another
    example is in nonlinear optics where such structures have appeared in
    optical waveguide arrays <cit.>.  Finally, and quite
    recently, yet another setup has emerged, that of tunable magnetic
    lattices <cit.> in which ultraslow shock waves can arise and
    be experimentally imaged.
    
    Earlier interest in lattice shocks include the heyday of conservation
    laws and shock waves in the fifties and sixties when material
    scientists were interested in the compression of a solid by passage of
    a very strong shock wave through materials
    <cit.>.  Early numerical studies (molecular
    dynamics simulations) depicted what we now call a lattice DSW in the
    material's stress profile and recognized its unsteady character in a
    one-dimensional anharmonic chain <cit.>.  This contradicted
    the basic assumption of steadiness underlying the Rankine-Hugoniot
    jump conditions and led to some controversy in the field.  The DSW's
    leading edge was then identified with a homoclinic traveling wave
    solution (solitary wave) of a continuum approximation in
    <cit.>, later identified as a generic feature of
    DSWs in continuum media <cit.>.  The
    controversy continued for about 15 years until, in 1979,
    three-dimensional lattice simulations were shown to exhibit a
    transition from unsteady to steady shock fronts due to transverse
    strains <cit.>.  It is worth noting that the
    same transition from one-dimensional, unsteady (dispersive) to
    multi-dimensional, steady (effectively viscous) shock propagation, was
    recently observed in a completely different, ultracold atomic
    superfluid <cit.>.
    
    
    
    These works have motivated the present authors to revisit “lattice
    hydrodynamics” and the prototypical settings where DSW structures can
    arise in nonlinear dynamical lattices. Canonical examples of first
    order in time, quadratically nonlinear lattice nonlinear ordinary
    differential equations (ODEs) in a conservation law form were
    discussed extensively in the work of <cit.>. A subset of the
    present authors has recently revisited this class of models
    in <cit.> attempting to incorporate tools from Whitham
    modulation theory, bringing to bear both data-driven, as well as more
    theoretically-inspired quasi-continuum approaches to obtain an
    effective dimensional reduction, through an ODE description of the DSW states.
    Our aim here is to expand on this work, offering a more systematic
    classification of the possible solutions of such models by using
    Whitham modulation theory and, where appropriate,
    quasi-continuum approximation considerations supplemented by
    direct numerical simulations.
    
    We select the scale invariant, representative nonlinear example
    (<ref>) within the class of models of <cit.> and set up the
    corresponding discrete Riemann problem (<ref>).  Figure
    <ref> depicts our phase diagram as a partitioning of
    the parameter space (u_-,u_+) and identifies seven distinct
    solution behaviors. Some of these, such as the possibility of a RW or
    a DSW as well as that of blow up are to a certain degree expected or
    have been argued to be present previously <cit.>.  They are
    labeled RW, DSW, and ∞ in Fig. <ref>,
    respectively.  However, there are choices of initial conditions that
    yield less common dynamical outputs, some of which are genuinely
    discrete in nature with labels in Fig. <ref>
    identified parenthetically.  These include, for instance, a
    stationary, symmetric shock on its own (SS) or separating a RW and a
    DSW (DSW+SS+RW).  Another example is a traveling DSW (TDSW), which
    consists of a partial DSW connected to a heteroclinic
    periodic-to-equilibrium traveling wave.  Arguably the most complex
    structure encountered is an unsteady shock evolving between two
    distinct traveling waves featuring the same temporal frequency (US).
    In what follows, we explain our partitioning of the phase diagram
    <ref> into the regions pertaining to these different
    dynamical behaviors and we offer a set of tools that can be used to
    understand each one, as well as unveil some open directions for future
    exploration.
    
    There exists a family of Riemann problems that can be considered in the
    discrete setting by modifying the value u_0(0) in (<ref>).  For
    example, Turner and Rosales set u_0(0) = 0 and u_- = -u_+<cit.>. We primarily focus on the data (<ref>) in which
    u_0(0) = u_- resulting in the phase diagram of
    Fig. <ref>.  While changing u_0(0) does not affect
    the observed solution phases, it does change the phase boundaries.  We
    interpret this microscopic modification of the initial data impacting
    the macroscopic properties of solutions as an indication of
    non-uniqueness of the Riemann problem for the dispersive
    regularization (<ref>).
    
    It is important to distinguish our use of the term “shock” or
    “shock wave” from the classical notion of discontinuous weak
    solutions of inviscid Burgers' equation (<ref>).  We identify
    four classes of shock solutions to the discrete equation (<ref>)
    by prefacing each with a descriptor in Fig. <ref>.
    The simplest is the symmetric, stationary shock (SS) solution
    of the lattice
    
        u_n(t) =
          
            -u_0     n ≤ 0, 
        
            u_0     n > 0,
    
    where u_0 > 0.  The other shock solutions can be understood as
    special solutions of the first-order, quasi-linear Whitham modulation
    equations corresponding to Eq. (<ref>) that are described in
    Sec. <ref>.  The unsteady dispersive shock wave
    (DSW) is approximated by a nonlinear, periodic wavetrain modulated by
    a rarefaction wave solution of the Whitham modulation equations.  Note
    that the DSW does not satisfy the Rankine-Hugoniot jump conditions of
    the Whitham modulation equations.  On the other hand, the unsteady
    shock (US) is approximated by two periodic traveling waves that
    satisfy the jump conditions for the Whitham modulation equations.  The
    traveling dispersive shock wave (TDSW) consists of both an unsteady
    partial DSW—approximated by a rarefaction solution of the Whitham
    modulation equations—and a steady traveling wave—approximated by a
    periodic traveling wave and a solitary wave that satisfy the jump
    conditions for the Whitham modulation equations.  There is an
    important distinction between the traveling wave and US as
    discontinuous solutions of the Whitham modulation equations.  The
    phase speeds of the two periodic traveling waves in the US solution
    differ from one another and from the shock speed, which is zero.  On
    the other hand, the traveling wave solution consists of a single
    periodic traveling wave whose phase speed is the same as the shock
    speed.  For clarity, we summarize the four distinct uses of the term
    “shock” in this paper:
    
      *  the stationary lattice shock (SS) (<ref>);
    
      *  the unsteady dispersive shock wave (DSW) that is approximated by
      a rarefaction solution of the Whitham modulation equations;
    
      *  the unsteady shock (US) that is approximated by a discontinuous
      shock solution of the Whitham modulation equations;
    
      *  the traveling dispersive shock wave (TDSW) that is approximated
      by a shock-rarefaction solution of the Whitham modulation equations.
    
    
    Our presentation will be structured as follows.  In section
    <ref>, we present the model equations, as well as
    the principal setup and notation for our study.  In section
    <ref>, we focus on the Whitham modulation equation
    formulation for the discrete problem. We discuss the corresponding
    conservation laws and how their averaging can provide information for
    the DSW features of our model.  Section <ref> is dedicated
    to the systematic classification of our solutions in the different
    parametric regimes, accompanied by illustrative numerical computations
    of the different identified waveforms.  In section
    <ref>, we show how modification of the Riemann data
    (<ref>) at a single site can lead to drastically different
    solution behaviors.  Finally, in section <ref>,
    we summarize our findings and present our conclusions, as well as a
    number of open questions for further research into this budding theme.
    
    
    
    
    
    § MODEL EQUATIONS
    
    
    It will be beneficial to generalize Eq. (<ref>) and consider the
    discrete scalar conservation law <cit.>
    
        2 du_n/dt +  Φ'(u_n+1 )-  Φ'(u_n-1 ) =0,
    
    a discretization of the more general conservation law
    u_t + Φ'(u)_x = 0, where n∈, t∈, u=u_n(t) ∈ and
    the potential Φ(u) is assumed to be smooth with Φ'(u) a
    convex function of its argument Φ”'(u) 
    0. Equation (<ref>) possesses a Lagrangian and Hamiltonian
    structure <cit.>, yet it is first order only, making
    its analysis slightly more convenient when compared to classical
    nonlinear oscillators, such as those of the Fermi-Pasta-Ulam-Tsingou
    (FPUT) type <cit.>.  Besides serving as a prototype model for
    lattice DSWs, Eq. (<ref>) is also of interest for applications,
    such as in the description of traffic flow <cit.>; for a
    discussion of relevant models and their continuum limits see also
    Ref. <cit.>.
     
    In this paper, we primarily focus on the potential
    
        Φ(u) = u^3/3 .
    
    For this choice, the “mass”
        M(t) = ∑_n u_n(t)
    
    and “energy”
        E(t) = ∑_n Φ(u_n(t)),
    
    when well-defined, are conserved in the infinite lattice and in a
    finite lattice with periodic boundary conditions.  The linear
    dispersion relation for Eq. (<ref>) is
    
        ω_0(k,u̅) = Φ”(u̅) sin(k),    k,u̅∈,
    
    for linearized wave solutions of the form
    u_n(t) = u̅ + a e^i(kn-ω_0 t), |a| ≪ 1.  Throughout
    the manuscript we consider the Riemann, step initial data
    (<ref>).
    
    
    
    
    
    
    For numerical simulations, the infinite lattice is truncated by
    introducing N > 0 (even) to represent the number of lattice
    sites. The corresponding spatial domain is -N/2+1 < n ≤ N/2 and
    the simulation temporal domain is [0, T_f / ], where T_f is a
    fixed constant independent of ϵ= 1/N. We use free boundary
    conditions u_-N/2 = u_-N/2+1 and u_N/2 = u_N/2-1 in
    conjunction with the initial data Eq. (<ref>), and we choose
    domain sizes large enough that interactions with the boundary are
    negligible.  When investigating finite time blow up, we employ
    periodic conditions. This allows us to monitor if the rescaled
    quantities E(t) → E(t)/N and M(t) → M(t)/N are conserved
    (details in sec. <ref>).  A variational integrator is used
    for simulations, see <cit.>. Simulations were also
    carried out with a Runge-Kutta method to check for consistency which
    yielded negligible differences on the time scales considered in this
    paper (for cases that did not involve blow up features).   Due to the scaling symmetry
    t → a t, u_n → a u_n, for any nonzero a ∈ of equation
    (<ref>) subject to (<ref>), we can set either u_+ = 1 or
    u_- = 1 without loss of generality. Figure <ref>
    shows a classification of the zoology of solutions that arise from the
    Riemann problem. They include
    
    rarefaction waves (RWs) for u_+= 1 and u_-∈(0.18,1), 
    
    dispersive shock waves (DSWs) for u_+ ∈ (0,1) and u_- = 1, 
    
    solutions consisting of dispersive shock waves, stationary shocks and rarefaction waves (DSW + SS + RW) for u_+= 1 and u_-∈(-1,-0.26), 
    
    solutions consisting of traveling dispersive shock waves (TDSW) for u_+ ∈ (-0.724,0) and u_- = 1 ,
    
    unsteady shocks (US) for u_+= 1 and u_-∈(-0.26,0.18), 
    
    and blow up (∞) for u_+ ∈ (-1,-0.724) and u_- = 1.
    
    The region boundaries are approximate. In sec. <ref>, we
    provide a detailed analysis for each of the five solution types just
    described, starting first with the simplest, and moving through them
    gradually in terms of their complexity according to the table in
    Fig. <ref>.  We employ a number of tools for the
    study of these solutions, including direct numerical simulation,
    fixed-point iteration schemes, modulation theory, weak solutions, DSW
    fitting, and quasi-continuum modeling.  The details of these
    approaches will be given in the sections they are employed, with the
    exception of modulation theory. This analysis is slightly more
    involved, and thus has a dedicated section. Our intention in
    presenting these tools is to leverage this specific, but interesting
    in its own right, example in order to utilize a variety of techniques
    that may be of broader relevance to applications in other Hamiltonian
    nonlinear dynamical lattices. It would be of particular
    interest to identify similar phenomena or/and to leverage
    the techniques utilized herein in other dispersive,
    nonlinear lattice models.
    
    
    
     §.§ An alternative, integrable discretization
    
    
    Prior to describing the solutions of Eq. (<ref>) depicted in
    Fig. <ref>, we briefly comment on the alternative
    discretization
    
        u_n/ t + u_n ( u_n+1 - u_n-1 ) = 0 ,
    
    of Eq. (<ref>) subject to (<ref>).  This equation was studied
    in <cit.> where it was shown to exhibit
    DSWs and, for positive data, to be completely integrable by a
    transformation <cit.> to an equation related to the Toda
    lattice <cit.>.  We have performed numerical simulations of
    Eq. (<ref>) subject to the Riemann data (<ref>) and observe
    DSWs when u_- > u_+ ≥ 0, RWs when u_+ > u_- ≥ 0 and blow-up
    when u_+ and u_- exhibit opposite signs. Examples of numerical
    simulations of DSWs and RWs that emerge from strictly positive initial
    Riemann data are shown in Fig. <ref>
    
    Of course, complete integrability confers a great deal of mathematical
    structure.  Whitham modulation theory for the Toda lattice was
    developed in <cit.> while
    the inverse scattering transform for the Toda lattice with step-type
    initial data was developed in <cit.>; see also the
    recent discussion of Whitham theory applied to DSWs in the Toda
    lattice <cit.>.  Collectively, these works
    support our numerical observation that, for positive Riemann data,
    Eq. (<ref>) exhibits only RW and DSW solutions.  These Riemann
    problem solution behaviors are to be contrasted with those depicted in
    Fig. <ref>.  Although integrability of
    Eq. (<ref>) is lost for sign indefinite initial data, the only
    dynamics we numerically observe are indicative of blow up.  Thus, the
    discretization (<ref>) we focus on in this paper, admits a wider
    variety of dynamics than the integrable alternative of (<ref>).
    
    
    
    § WHITHAM THEORY
    
    
     §.§ Modulation equations for a continuum system
     
    
    In this section, we consider a continuum model system by introducing
    the interpolating function u(x,t) such that u(n,t) = u_n(t) for
    all n ∈ℤ. This allows us to represent the advance-delay
    operator in the discrete system (<ref>) as a pseudo-differential
    operator. The resulting continuum model is
    
        u_t + isin(-i∂_x) (Φ'(u)) = 0  ,
    
    Equation (<ref>) can be written in the Hamiltonian form 
    
        u_t = J δ H/δ u,
    
    where J = i sin(-i∂_x) is the antisymmetric operator and
    H = ∫Φ(u)  dx is the Hamiltonian. Periodic solutions of
    (<ref>) are of the form u(x,t) = φ(θ;𝐪)
    with phase θ = kx - ω t and parameters
    𝐪∈^3 (e.g., wavenumber, amplitude, mean).  They satisfy
    
        - ωφ_θ  + isin(-ik∂_θ)
          Φ'(φ) = 0 ,
    
    which is equivalent to the nonlinear advance-delay differential equation
    
        2 ω φ_θθ=Φ'(φ(θ+k))-Φ'(φ(θ-k))
    
    The solution theory of such nonlocal equations is rather intricate but
    the existence of a three-parameter family of traveling waves has been
    establshed in <cit.> by variational
    techniques. Integrating (<ref>) once with respect to θ
        -c_ pφ  +  sinc(-ik∂_θ) Φ'(φ)= A,
    
    where c_ p = ω/k is the phase speed and A is a
    real constant. The pseudodifferential operator is then interpreted as
    a multiplier on the Fourier coefficients of φ,
    
        sin(-ik∂_θ) Φ'(φ)    = ∑_nsin(n k) p̂_n e^in θ, 
        p̂_n     = 1/2π∫_0^2πΦ'(φ(θ))e^-inθ dθ.
    
    Equation (<ref>) possesses the two conserved quantities
    
        M(t)    = ∫ u dx 
        
        E(t)    = ∫Φ(u) dx,
    
    where the domain of integration is determined by the decay or
    periodicity of u.  We now seek the modulation equations for a
    periodic wave with the slowly varying ansatz
    
        u(x,t)= φ(θ;𝐪(X,T)) + ϵφ_1(θ,X,T) + ⋯ ,    X =ϵ x,    T
          =ϵ t,   
          0<ϵ≪ 1,
    
    in which the leading order term φ(θ;𝐪) is the
    periodic traveling wave solution satisfying (<ref>) with vector
    of parameters 𝐪 that varies on the slow scales X and T
    while φ is 2π-periodic in θ. We impose the
    generalized wavenumber and frequency relationships θ_x = k and
    θ_t = -ω along with their compatibility
    
        k_T + ω_X = 0 .
    
    The nonlocal operator acting on a modulated periodic function g(θ,X,T) ∈𝒞^1 has the multiple scale expansion 
    
        sin(-i∂_x) g    = sin(-i k∂_θ - iϵ∂_X) g
           ∼sin(-i k ∂_θ) g - i ϵ/2(cos(-i k ∂_θ) g_X + (cos(-i k ∂_θ) g)_X) + 𝒪(ϵ^2)
    
      The proof follows from the analyticity of sin(·). A detailed
      proof follows all of the ideas in <cit.>.
    
    
    We now average Eq. (<ref>) and its higher order conserved
    densities by introducing the averaging operator
    
        F[φ](X,T) = 1/2π∫_0^2π
          F[φ(θ;𝐪(X,T))] dθ ,
    
    where F[u(x,t)] = F(u,u_x,u_t,u_xt, …) is a local function of
    u and its derivatives.  If u is a multiscale
    function of the form u = g(θ,X,T), then
    
        ∂_t F   = -ω∂_θ F + ϵ∂_T F = ϵ∂_T F, 
        ∂_x F   = k ∂_θ F + ϵ∂_X F = ϵ∂_X F,
    
    by virtue of the fact that F = F[g] is periodic in θ so the
    period average of ∂_θ F[g] is zero. 
    We use Lemma
    <ref> to compute averages of, for example,
    sin(-i k ∂_θ)g for any g ∈ L_2([0,2π]) with the
    Fourier series g = ∑_n ĝ_n e^inθ:x
    
        sin(-i k ∂_θ)g    = ∑_n sin(n k) ĝ_n e^inθ
         = 0 .
    
    We now insert the multiple scales ansatz (<ref>) into the two
    conservation laws associated with (<ref>), and average. This
    procedure results in the system of conservation laws
    
        φ_T + Φ'(φ)_X     = 0 
        Φ(φ)_T + 1/2(∑_m cos(mk)
            p̂_m^2)_X     = 0 
        
            k_T + ω_X     = 0
    
    
    In the vanishing amplitude a → 0 limit, Eqs. (<ref>) and
    (<ref>) become the Hopf equation
    
        u̅_T + Φ”(u̅)u̅_X = 0,
    
    for the mean u̅ = φ, and the conservation of
    waves equation (<ref>) corresponds to linear wave modulation
    theory with frequency ω = ω_0 given by the linear
    dispersion relation (<ref>).
    
    The nonlinear modulation equations can alternatively be derived by
    employing Whitham's other method of an averaged Lagrangian functional,
    see for instance <cit.> and
    <cit.> for symplectic PDEs,
    <cit.> for an application to FPUT chains as the
    most prominent example of Hamiltonian lattices, and
    <cit.> for the discrete conservation law
    (<ref>). In this setting, the modulation equations take the form
    
        u_T+(E_ u)_X    =0
        
           k_T + (E_S)_X   =0
            
        
            
            S_T + (E_k)_X     = 0  .
    
    This is a system of Hamiltonian PDEs with density variables u, k,
    and S, which represent the wave mean, the nonlinear wave number, and a
    nonlocal auxiliary variable that might be regarded as a generalized
    wave momentum. Moreover, the equation of state E=E(u,k,S) describes
    the energy of a traveling wave and its partial derivatives provide the
    fluxes in (<ref>). The energy is also conserved
    according to the extra conservation law
    
        E_T+(12 E_u^2+E_kE_S )_X=0 ,
    
    which is implied by (<ref>) thanks to the chain rule.
    A closer look to the derivation of (<ref>) in <cit.> reveals that (<ref>) and (<ref>) correspond to (<ref>) and 
    (<ref>), respectively, while (<ref>) is the analogue to (<ref>).
    A complete 
    understanding of (<ref>) and (<ref>) is currently out of reach because we are not able to characterize the analytical properties of the constitutive relations since these depend in a very implicit and not tractable way on the three-dimensional solution sets of the nonlinear advance-delay-differential equation (<ref>). For instance, it is not even clear for which values of the parameters the Whitham system (<ref>) is hyperbolic or genuinely nonlinear. For this reason we do not work with the full lattice modulation equations directly but combine different approximation procedures with a careful evaluation of numerical data.
    
    
    
    
    
     §.§ Relation to the lattice dynamics
    
    
    Although neither analytical nor numerical solutions to the nonlinear modulation systems 
    (<ref>) or (<ref>) are available, we
    can extract important partial information from numerical simulations
    of initial value problems to (<ref>). The key observation is that
    the lattice ODE as well as an implied energy equation represent
    discrete counterparts of local conservation laws and transform under
    the hyperbolic scaling of space and time into first order PDEs. To see
    this, we fix a window functionχ that depends smoothly
    on the macroscopic variables (X,T), decays sufficiently fast, and
    has normalized integral. Using a shifted copy of χ, we are able
    to quantify the local moments of any microscopic observable near a fixed macroscopic point. For instance, the average
    
        ⟨ u_n⟩XT=^2∑_n̆∫ u_n̆(t̆ ) χ(t̆-T,n̆-X) dt̆ ,     T=  t ,     X=  n
    
    represents the mesoscopic space-time averages of u_n(t) near the macroscopic point (X,T). 
    
    
     Any bounded solution to (<ref>) satisfies in the hyperbolic scaling limit →0 the conservation laws
    
        ∂_T ⟨ u_n ⟩  + ∂_X ⟨Φ'(u_n) ⟩ =
            0, 
        ∂_T ⟨Φ(u_n)⟩ + ∂_X
            ⟨12 Φ'(u_n) Φ'(u_n+1)⟩ =0,
    
    provided that these are interpreted in a distributional sense.
     We only give an informal derivation but mention that an alternative and more elegant framework is provided by the theory of Young measures. The latter can also be applied to non-smooth window functions χ and reveals that the mesoscopic averages ⟨·⟩ can be expected to be independent of the particular choice for χ.
    Using the abbreviation p_n = Φ'(u_n), discrete integration by parts as well as the smoothness of χ we verify
    
        ⟨12  p_n+1-12  p_n-1⟩XT   =
            ^3∑_n̆∫ p_n̆(t̆ ) (12 χ(t̆-T,n̆--X)-
            12 χ(t̆-T,n̆+-X) ) dt̆
           =-
            ^3∑_n̆∫ p_n̆(t̆ ) ∂_X χ(t̆-T,n̆-X) dt̆+h.o.t.
           =
            - ∂_X ⟨ p_n⟩XT+h.o.t.
    
    and by similar computations we obtain
    
        ⟨u̇_n+1⟩XT =
            - ∂_T⟨ u_n⟩XT+h.o.t. .
    
    The asymptotic validity of (<ref>) is thus a direct consequence of the microscopic dynamics (<ref>), the definition of the bracket ⟨·⟩ in (<ref>), and the hyperbolic scaling. The lattice ODE (<ref>) implies
    with
    
        d/d tΦ(u_n)+ ( 12 Φ^'(u_n) Φ^'(u_n+1))-( 12 Φ^'(u_n-1) Φ^'(u_n) )=0 ,
    
    another discrete conservation law (in which the time derivative of a
    density is balanced by the discrete divergence of a flux quantity), so
    the second claim (<ref>) can be justified along the same
    lines.
    
    
    
    There is an important difference between the conservation laws in
    (<ref>) and (<ref>).  The PDEs in
    (<ref>) (and likewise those in
    (<ref>)) are derived under the hypothesis
    that the lattice solution can be approximated by a modulated traveling
    wave, see (<ref>), and the closure relations involve the
    (unknown) profile functions for traveling lattice waves as well as
    averages with respect to the scalar phase variable θ.
    Numerical simulations with well-prepared initial data (e.g., the
    Riemann initial data (<ref>)) indicate that the approximation
    assumption concerning the microscopic data is indeed satisfied but no
    rigorous proof is available, neither for the lattice (<ref>) nor
    for FPUT chains with convex interaction potential. The only exceptions are
    the few completely integrable cases but the details are still
    complicated and involve special coordinates related to the Lax
    structure. In particular, even for the lattice of Eq. (<ref>)
    and the Toda chain
    it is not easy to compute how the phase averages in the modulation
    equations depend on the traveling wave parameters.
    
    
    
    The status of (<ref>) is completely different. The two
    PDEs can be established under very mild assumptions (boundedness of
    lattice solutions) and by means of fundamental mathematical principles
    (such as integration by parts and compactness in the sense of Young
    measures). They reflect universal constraints for the
    macroscopic dynamics, do not require any a-priori knowledge on the
    fine structure of the microscopic oscillations, and hold for a large
    class of initial data (which might even be oscillatory or random).  
    
    Moreover, the
    mesoscopic space-time averages can easily by extracted from numerical
    data. In the simplest case, we use a straightforward box counting with
    space-time windows of microscopic length 1/√() (or 
    macroscopic length √()).  Of course, (<ref>)
    does not provide a complete set of macroscopic equations and without
    further information it is not clear whether or how the fluxes can be
    computed in a pointwise manner from the densities. The equations are
    nevertheless very useful since they allow us to derive and check
    partial information on the solution of the modulation equations from
    numerical data. In particular, in the context of modulated traveling
    waves, the PDEs (<ref>) and (<ref>)
    correspond to (<ref>) and (<ref>), respectively.
     
    
    
    
     §.§ Self-similar solutions
    
    
    The Whitham modulation equations (<ref>) are a system of
    conservation laws that can be compactly expressed in the form
    
        𝐏(𝐪)_T + 𝐐(𝐪)_X = 0,   𝐪 = [u̅,a,k]^ T,
    
    where the vectorial density 𝐏 and flux 𝐐 depend
    on the slowly varying parameters 𝐪 through integrals of the
    periodic orbit φ.  Equation (<ref>) can also be
    expressed in the form
    
        𝐪_T + 𝒜𝐪_X = 0,   𝒜 =
           ( ∂𝐏/∂𝐪 )^-1∂𝐐/∂𝐪,
    
    provided the inverse is nonsingular.  We will use solutions of the
    Whitham equations to approximate the long time dynamics of solutions
    to the Riemann problem (<ref>), (<ref>).  Consequently, it
    is natural to consider the Riemann problem
    
        𝐪(X,0) = 𝐪_-     X < 0
        𝐪_+     X > 0 
                   ,
    
    for the Whitham equations (<ref>) themselves.  Rarefaction
    (simple) wave solutions and discontinuous shock solutions of the
    binary oscillation modulation system (<ref>) were used in
    <cit.> to interpret various features of the numerical solutions.
    In this work, we will make use of rarefaction wave and discontinuous
    shock solutions of the more general Whitham modulation equations
    (<ref>).
    
    The invariance of the Riemann problem (<ref>), (<ref>)
    with respect to the hydrodynamic scaling X →σ X',
    T →σ T' for real σ 0 suggests seeking
    self-similar solutions in the form 𝐪 = 𝐪(ξ),
    ξ = X/T.  Equation (<ref>) possesses rarefaction waves
     satisfying <cit.>
        𝐪/ξ = 𝐫_i/∇λ_i
            ·𝐫_i,   𝐪(ξ_±) = 𝐪_±,
            ξ_- < ξ_+,
    
    where 𝒜𝐫_i = λ_i 𝐫_i and
    λ_i = ξ, provided the characteristic field is genuinely
    nonlinear ∇λ_i ·𝐫_i  0.  Since
    𝐪_± lie on the same, one-dimensional integral curve, they
    are constrained by two integral relations resulting from integration
    of the third order ODEs (<ref>).  Admissibility requires
    ξ_- < ξ_+.  The eigenvalues λ_i can be interpreted as
    speeds. For example, in the context of DSWs, the trailing edge speed
    is c_-=ξ_- and the leading edge speed c_+ = ξ_+.
    
    
    Another class of self-similar solutions are discontinuous shock
    solutions to the Whitham system (<ref>)
        𝐪(ξ) = 𝐪_-    ξ < V 
        𝐪_+    ξ > V 
                   ,
    
    where V is the velocity of the shock solution that satisfies the
    Rankine-Hugoniot jump conditions
    
        - V [[ 𝐏 ]] + [[ 𝐐 ]] = 0 .
    
    The brackets [[ · ]] denote the jump in its argument evaluated on
    the left and right triple 𝐪_± that parameterize distinct,
    steady periodic orbits φ_±. 
    
    For strictly hyperbolic, genuinely nonlinear Whitham modulation
    equations with negative linear dispersion
    (∂_k^2 ω_0 < 0), classical DSW solutions connecting
    the two constant states u_± are described by a rarefaction
    solution of (<ref>) in which λ = λ_2 is the middle
    characteristic speed and 𝐪_- = [u_-,0,k_-],
    𝐪_+ = [u_+,a_+,0].  The two constraints that result from
    integrating (<ref>) determine the trailing edge wavenumber k_-
    and speed ξ_- as well as the leading edge amplitude a_+ and
    speed ξ_+<cit.>.  Therefore, a classical DSW corresponds
    to a rarefaction wave solution of the modulation equations,
    not a shock solution.  For DSW construction, we will use the
    DSW fitting method, which leverages certain structural properties of
    the Whitham modulation equations under the assumptions of strict
    hyperbolicity and genuine nonlinearity in order to obtain k_-,
    a_+, and ξ_± by integrating a scalar ODE
    <cit.>.
    
    Whitham himself pondered the notion of discontinuous shock solutions
    to his eponymous equations <cit.>. But their utility was
    only recently discovered in <cit.> where shock solutions
    of the Whitham modulation equations for a fifth order KdV (KdV5)
    equation were deemed admissible if there exists a heteroclinic
    traveling wave solution connecting the corresponding left and right
    periodic orbits, each moving with the same speed as the shock.  Such
    traveling wave solutions are possible in higher order equations such
    as KdV5.  These Whitham shocks were used to solve the Riemann
    problem for KdV5 and, later, were investigated in the Kawahara
    equation <cit.>.  In this paper, we will show
    that similar Whitham shocks emerge as the traveling wave portion of
    the TDSW solution in Fig. <ref>.  We also provide
    analytical and numerical evidence of the existence of a new class of
    Whitham shocks, i.e., shock solutions of the Whitham modulation
    equations (<ref>) whose corresponding left and right
    periodic orbits possess the same frequency but different speeds than
    one another and the shock itself (see US in
    Fig. <ref>).
    
    
    
     §.§ Weakly nonlinear regime
    
    
    In the previous sections, we derived the modulation equations supposing
    the existence of a family of nonlinear periodic solutions. In the
    case where no known explicit periodic solution is available,
    it is useful to approximate the periodic solution with a truncated
    cosine series. The approximation via the Poincaré-Lindstedt method
    utilizes an asymptotic expansion of both the profile of the periodic
    solution and its frequency in the small amplitude parameter
    0 < a ≪ 1. The approximate periodic solution and its frequency are
    given, for a generic potential Φ by
    
        u    ∼u̅ + a/2cos(kn - ω t) +  a^2 sin (2 k) Φ^(3)(u̅)/16 Φ”(u̅)(2 sin (k) -sin (2 k) )cos(2(kn - ω t)) + o(a^2),
        ω   ∼Φ”(u̅) sin(k) + a^2 ω_2 + o(a^2),   ω_2 = 1/32sin (k) (Φ^(3)(u̅)^2/( (k)-1) Φ”(u̅)+Φ^(4)(u̅)) ,
    
    which maintain their asymptotic ordering so long as
    
        a^2/|k| ≪ 1   and  
          |a Φ^(3)(u̅)/Φ”(u̅)| ≪ 1,
    
    i.e., for Φ(u) = u^3/3, neither |k| nor |u̅| are too
    small.  Inserting (<ref>), (<ref>) into
    the modulation equations (<ref>), we obtain the weakly
    nonlinear Whitham modulation equations in conservative form by
    retaining terms up to O(a^2)
        u̅_T + (Φ'(u̅) + a^2/16Φ”'(u̅))_X     = 0, 
        (Φ(u̅) + a^2/16Φ”(u̅))_T 
        + (1/2Φ'(u̅)^2 + a^2/16(Φ”(u̅)^2 cos(k) + Φ'(u̅)Φ”'(u̅)))_X     = 0, 
        
        k_T + (Φ”(u̅) sin(k) + a^2/32sin (k)
          (Φ^(3)(u̅)^2/( (k)-1)
          Φ”(u̅)+Φ^(4)(u̅)))_X     = 0 .
    
    In the case of the cubic potential (<ref>), our focus here, the
    modulation equations are
    
        u̅_T + (u̅^2 + a^2/8)_X     = 0, 
        (u̅^3/3 + a^2/8u̅)_T 
            + (1/2u̅^4+ a^2/8u̅^2 (2
            cos(k) + 1))_X     = 0, 
        
            k_T + (2u̅sin(k) + a^2/16(sin(k)/( (k)-1) u̅))_X     = 0.
    
    Properties of the modulation equations can be elucidated by casting
    them in quasi-linear form
    𝐪̃_t + 𝒜̃𝐪̃_x = 0,
    where
    
        𝒜̃ =
          [                                                2 u̅                                                 1/8                                                   0;                                       4 a^2 cos (k)                                         2 u̅cos (k)                                    -2 a^2 u̅sin (k);                              2 sin (k) + a^2 ω_2,u̅                                                 ω_2 2u̅cos(k) + a^2
                                      ω_2,k ],
                                          𝐪̃ =
                                      [  u̅; a^2;   k ].
    
    A perturbation calculation gives the
    eigenvalues of the flux matrix 𝒜 to 𝒪(a)
        λ_3     = 2 u̅ +O(a^2), 
        λ_2     = 2u̅cos(k) +
                        a/2cos(k/2)√(2-cos(k)) +O(a^2),
            
        λ_1     = 2u̅cos(k) -
                        a/2cos(k/2)√(2-cos(k))  +O(a^2)
                        ,
    
    with the corresponding right eigenvectors
    
        𝐫_3   =  [
                 u̅tan(k/2), 0, 1  ]^T
                 +O(a^2),
        𝐫_2   =
                  [0, 0, √(2-cos (k)) ]^T
                 + a/4 [ (k/2), - 32 u̅sin(k/2), 0   ]^T
                 + O(a^2), 
        𝐫_1   = [0, 0, √(2-cos (k)) ]^T - a/4 [ (k/2), -32 u̅sin(k/2),
                 0   ]^T
                 + O(a^2).
    
        
    The quasilinear system is strictly hyperbolic if all of the
    eigenvalues are distinct, and real valued. To the order of the
    approximation given, the weakly nonlinear system is strictly
    hyperbolic provided a  0, k π, and
    
        |u| u_ cr,    u_ cr =
          acos ( k/2 ) 
            √(2-cos(k))/8 sin^2 ( k/2 ) .
    
    When u > u_ cr and a > 0, the eigenvalues are
    ordered λ_1 < λ_2 < λ_3.  When a = 0,
    Eqs. (<ref>) and (<ref>)
    coincide with the Hopf equation for the mean u̅ and the
    remaining equation corresponds to the conservation of waves from
    linear wave modulation theory.  When k = π,
    Eq. (<ref>) is identically satisfied.  While
    intuition might suggest that (<ref>) and
    (<ref>) are somehow related to the modulation
    equations for binary oscillations (<ref>), in fact, the
    asymptotic derivation breaks down.  For example, the period average of
    the weakly nonlinear solution (<ref>) is no longer
    u̅ but rather u̅ - a^2/8u̅ so that the
    density in (<ref>) does not correspond to the
    density in (<ref>).  When k = π, one should discard
    Eq. (<ref>) altogether in favor of the modulation
    equations for binary oscillations (<ref>), which apply beyond the
    weakly nonlinear regime considered here.
    
    
    
    § CLASSIFICATION OF SOLUTIONS
    
    
    From now onwards, we focus solely on the discrete equation
    (<ref>) (Eq. (<ref>) subject to (<ref>)).  Figure
    <ref> depicts seven qualitatively distinct solution
    families to the Riemann problem (<ref>), (<ref>) depending
    upon the parameters u_± in the initial data.  We now proceed to
    describe each of these solution families using a combination of
    numerical simulation, Whitham modulation theory, and quasi-continuum
    approximation.  The straight line boundaries between each solution
    family in Fig. <ref> are determined empirically (to
    two decimal digits accuracy) and some are explained by analytical
    considerations. By a possible reflection of the lattice n → -n and
    a rescaling of time, we can, without loss of generality, set either
    u_+ = 1 while varying u_- ∈ [-1,1] or set u_- = 1 while
    varying u_+ ∈ [-1,1].  Therefore, we can map out the phase diagram
    in the (u_-,u_+)  plane by traversing the top and right edges of the
    square [-1,1]^2.
    
    The special case in which u_+ = u_- is trivial but the case in which
    u_+ = -u_-  0 is the stationary shock (SS) solution
    (<ref>).  Otherwise, the solutions exhibit more complexity,
    which we now explore. We start with the simplest case first, and then
    work toward the richest, most complex scenario.
    
    
    
     §.§ Rarefaction waves (RWs)
    
    
    The simplest observed dynamical structure is the rarefaction wave
    shown as RW in Fig. <ref>.  Empirically, we find that
    they form when u_+= 1 and u_-∈(0.18,1).  The bifurcation at
    u_- = 0.18 will be described in section <ref>. The leading
    order RW behavior is given by the self similar solution
    (ξ = n/t = X/T)
    
        u_n(t) ∼u̅(ξ) =
          
            u_-    ξ≤ 2u_- 
        ξ/2     2u_-  < ξ≤ 2u_+  
         
            u_+     2u_+ < ξ,
    
    of the dispersionless equation (<ref>).  A favorable comparison
    of this profile with a numerical simulation is shown in
    Fig. <ref>. Because the data is expansive, the effect of
    dispersion manifests at higher order where a small amplitude,
    dispersive wavetrain is emitted from the lower, left edge of the RW.
    
    The slowest (most negative) group velocity is
    ∂_kω_0(π,u_-) = -2 u_-, which corresponds to an inflection
    point of the linear dispersion relation (<ref>). 
    Consequently,
    the leftmost edge of these small amplitude waves is expected to have
    an Airy profile whose decay estimate is proportional to t^-1/3,
    similar to the Fourier analysis carried out for linear FPUT chains
    <cit.>. The details of the linear wavetrain
    accompanying RWs and DSWs for the BBM equation were studied in
    <cit.>.  We follow a similar procedure by
    linearizing about the left initial state u_n = u_- + v_n to obtain 
    
        d/d tv_n + u_-(v_n+1 - v_n-1) = 0.
    
    The initial data (<ref>) then becomes
    
        f_n =
          
            0     n < 0, 
         u_+ - u_-     n ≥ 0 ,
    
    whose discrete-space Fourier transform is the distribution
    
        f̂(k) = ∑_n=-∞^∞ f_n
          e^-ink = (u_+ - u_-)  (1/1 - e^-ik +
          πδ(k)  ),    k ∈ (-π,π],
    
    where δ(k) is the Dirac delta.  To approximate the nonlinear
    equation (<ref>) by the linear equation (<ref>),
    one could seek solutions in which 0 < u_+ - u_- ≪ |u_+|+|u_-|.
    Alternatively, we follow <cit.> and consider
    scale separation in which the highest frequency components of
    (<ref>) are assumed to separate from the RW so that the initial
    data becomes
    
        v̂(k,0) =
          u_+ - u_-/1 - e^-ik    k_0 < |k| ≤π, 
        
            0    else ,
    
    for some 0< k_0 < π that is sufficiently far from the zero dispersion points, k = 0,π.  Then, the solution of the linear equation
    (<ref>) can be determined by taking the discrete-space
    Fourier transform v̂(k,t) = ∑_n v_n(t) e^-ink.  The
    solution of Eq. (<ref>) subject to (<ref>) is
    
        v_n(t) = 1/2π∫_-π^πv̂(k,0) e^iθ(k)t  dk     θ(k;n,t,u_-) = k n/t- 2u_-sin(k) .
    
    Quantitative information regarding the solution can be determined
    asymptotically for t →∞ with n/t fixed using the method of
    stationary phase <cit.>. The leading order behavior is
    determined by analyzing the integral (<ref>) near the
    stationary points, k_s where θ_k(k_s) = 0. Stationary points
    are therefore given by ± k_s where
    
        n/t = 2 u_- cos(k_s) ,
    
    for -2u_- < n/t < 2u_-. The leading order behavior in the vicinity
    of the stationary points is determined by expanding the integrand in
    (<ref>) about the stationary points k = ± k_s. When
    k_s π and |k_s| > k_0, the leading order behavior is
    
        v_n(t)    ∼1/√(2π t |∂_k^2 ω_0(k_s,u_-)|)(v̂(k_s,0) e^i k_s n - i ω_0(k_s,u_-) t + i
                   π/4 +  c.c.), 
            = 1/√(2π t |∂_k^2 ω_0(k_s,u_-)|) (u_+ - u_-)
                   (k_s2)sin(θ(k_s)t +
                   π/4 + k_s2) .
    
    The profile (<ref>) is compared with numerical
    simulations of the initial value problem in
    Fig. <ref> on the interval
    [-u_- t_f, u_-t_f] at a final simulation time of t = t_f =
    1000. The interval is chosen so that k_s ∈ (π/3,2π/3), i.e.,
    the truncation parameter k_0 = π/3 and we avoid the degenerate
    stationary points k_s = 0, π.  We observe that the linear profile
    (<ref>) is in good agreement with the numerical
    simulation. However, for larger initial jumps, the linear wave begins
    to deviate from the simulation.  This may be attributed to the
    emergence of stronger nonlinear effects not captured by the leading
    order asymptotics which require a larger truncation parameter k_0.
    
    
    To investigate the leftmost edge of the linear wave emitted from the
    RW, we modify our previous analysis and expand the phase in the
    integral (<ref>) about the inflection point
    k = π of the linear dispersion relation
    
        θ(k) ∼π n/t - (n/t + 2u_-)(π-k) + u_-(π-k)^3/3 +
          ⋯
    
    The expansion (<ref>) is inserted into the
    integral (<ref>). A calculation reveals that the leading
    order asymptotics in the vicinity of the ray n/t = -2u_- are given
    by
    
        v_n(t) ∼ - u_+-u_-/2(t u_-)^1/3cos(π n)
           Ai( - (n + 2 t u_-) (t u_-)^-1/3),
    
    where Ai(·) is the Airy function
    Ai(z) = 1/2π∫_ℝ e^i κ z + i
      κ^3/3  dκ . The Airy profile (<ref>)
    favorably compares with the two Riemann problem simulations depicted
    in Fig. <ref>, even for large u_+ - u_-.
    
    
    
    It is worth contrasting the observed RW dynamics with those of the
    quasi-continuum approximation in the BBM equation (<ref>) that
    was studied in <cit.>.  Qualitatively, the
    dynamics exhibited by the two models in overlapping regimes of the
    (u_-,u_+) plane of Riemann data are very similar.  Both equations
    exhibit large-scale dynamics that are well-approximated by the
    self-similar solution (<ref>) and its analogue for the BBM
    equation.  The details of the short-scale, emitted dispersive
    wavetrains are quantitatively different but, since both equations
    admit non-convex linear dispersion relations, they both exhibit Airy
    profiles with amplitude decay proportional to t^-1/3.
    
    The long-time dynamics produced by the lattice model (<ref>)
    significantly differs from those generated by its quasi-continuum BBM
    counterpart (<ref>) when either equation is strongly influenced by
    small scale effects.  The actual Riemann problems for the BBM equation
    studied in <cit.> were tanh-smoothed,
    monotone transitions between u_- and u_+, 
    a feature which introduces an
    external length scale characterizing the width of the initial
    transition.  When this width is larger than the 𝒪(1)
    oscillatory length scale (or 𝒪(ϵ) in
    Eq. (<ref>)), the BBM equation exhibits a RW for all
    |u_-| < u_+.  As shown in Fig. <ref>, RW generation
    on the lattice is limited to the region 0.18 u_- < u_+ < u_-,
    u_- > 0, with short-scale oscillatory dynamics occurring when
    -u_- < u_+ < 0.18 u_-, u_- > 0.  When the BBM initial transition
    width is sufficiently small, RW generation can be accompanied by the
    spontaneous generation of solitary waves and/or an expansion shock,
    features not observed in the lattice model.
    
    
    
     §.§ Dispersive shock waves (DSWs)
    
    
    For u_- = 1 and 
    u_+∈(0,1)
    the data is compressive and results
    in an expanding, modulated oscillatory wave train between the states
    u_- and u_+. This structure is called a dispersive shock wave; see
    the panel labeled DSW in Fig. <ref>. In
    <cit.>, DSWs were studied in Eq. (<ref>)
    with Φ(u) = u^2/2 + u^4/4 and u_- = 1, u_+ =0 using numerical
    simulations and a dimension reduction approach. In the following, we
    study DSWs in the system (<ref>) as the step value u_+ varies
    using two semi-analytical approaches, DSW fitting and a continuum
    model.  
    
    
    
      §.§.§ Approximation of the DSW harmonic and soliton edge
      speeds via DSW fitting
    
    
    In this section, we outline the method for fitting the macroscopic DSW
    properties (edge speeds, amplitudes, and wavenumbers) by applying the
    fitting method first introduced by El <cit.>; see also
    <cit.>.  This method was originally developed and justified for continuum PDEs where it has been extensively applied.  Since it only requires knowledge of the linear dispersion relation and the solitary wave amplitude-speed relation, it is straightforward to extend the method to the semi-discrete lattice equation (<ref>).
    
    An example DSW from a numerical simulation of the
    Riemann problem with u_+ = 0.5 and u_-=1 is shown in
    Fig. <ref>.  The DSW is comprised of a modulated,
    nonlinear wavetrain that terminates in two distinct limits: vanishing
    amplitude (called the harmonic edge) and vanishing wavenumber (called
    the soliton edge).  The modulation solution describing the DSW is the
    rarefaction solution (<ref>) of the Whitham modulation equations
    (<ref>) with 𝐪_- = [u_-,0,k_-]^T,
    𝐪_+ = [u_+,a_+,0]^T.  The harmonic edge wavenumber k_-
    and the soliton edge amplitude a_+, as well as their corresponding
    edge speeds ξ_- and ξ_+ are determined by integrating the ODE
    (<ref>), thus relating these macroscopic DSW properties to the
    initial data u_±. In numerical simulations, the amplitude of the
    DSW does not vanish exactly at the harmonic edge, so we define the
    location of the trailing edge by the intersection of the oscillatory
    envelope curves shown in Fig. <ref>.
    
    
    Although the DSW modulation is determined, in principle, by the
    aforementioned rarefaction solution of the Whitham modulation
    equations (assuming strict hyperbolicity and genuine nonlinearity of
    the second characteristic field), we do not have explicit 
    expressions relating
    the integrals in (<ref>) to the periodic orbit's
    parameters 𝐪 = [u̅,a,k]^T.  An alternative technique
    that allows one to obtain the DSW's edge properties is the DSW fitting
    method <cit.>.  This method assumes the existence of the
    rarefaction solution. For the sake of completeness, we will carry out
    the DSW fitting procedure for a generic potential Φ(u). The
    wavenumber of the DSW at the harmonic edge can be determined by solving
    the initial value problem
    
        k/u̅ = ∂_u̅ω_0(k,u̅)/Φ”(u̅) -
            ∂_kω_0(k,u̅) = Φ”'(u̅)sin
            k/Φ”(u̅)(1-cos k),    k(u_+) = 0 ,
    
    where ω_0 is the linear dispersion relation (<ref>).
    This ODE can be integrated by separation of variables to obtain
    
        k(u̅) = cos^-1 ( 2
          Φ”(u_+)-Φ”(u̅)/Φ”(u̅) ).
    
    The wavenumber at the DSW harmonic edge is k_- = k(u_-). The
    velocity of the harmonic edge is given by evaluating the linear group
    velocity at k_-
        c_- = ∂_kω_0(k_-,u_-) = 2Φ”(u_+) - Φ”(u_-).
    
    
    The velocity of the DSW at the soliton edge is calculated in a similar
    way. We begin by introducing the conjugate variables
    ω_0(k,u̅) = - i ω_0(i
    k,u̅) = Φ”(u̅) sinh(k), where
    k acts as a soliton amplitude parameter. The velocity of the
    DSW soliton edge is deduced by evaluating the solitary wave dispersion
    relation c_+ = ω(k_+,u_+)/k_+, where we
    find k̃_+ by solving the initial value problem
    
        d k/d u̅ =
          ∂_u̅ω̃_0/Φ”(u̅) -
          ∂_kω_0 = Φ”'(u̅)
          sinhk/Φ”(u̅)(1-coshk) ,      k(u_-) = 0.
    
    Integration results in
    
        k̃(u̅) = cosh^-1 ( 2 Φ”(u_-) -
              Φ”(u̅)/Φ”(u̅) ) ,
    
    and the soliton edge conjugate wavenumber
    k̃_+ = k̃(u_+).  The soliton edge velocity of the
    discrete DSW is given by
    
        c_+ =  ω̃(k̃_+,u_+)/k̃_+ =
          2/k̃_+√(Φ”(u_-)(Φ”(u_-)-Φ”(u_+)))
    
    Comparisons with numerical simulations of the Riemann problem are
    given in Fig. <ref> with the potential
    Φ(u) = u^3/3 and the initial data normalized so that u_- =
    1. 
    To estimate the velocity of the leading edge, we track the position of
    the right-most lattice site that is above the far-field initial data
    u_+ at the integer valued times in our numerical simulation. This
    time series data is fit with a line whose slope is the approximate
    velocity of the soliton edge of the DSW. To estimate the harmonic edge
    velocity, we produce a linear fit of the modulated wavetrain amplitude
    near the location of the harmonic edge, which is found by extracting
    peaks of the solution at output times. The intersection of this linear
    approximation with the constant level u_- is the approximate
    location of the DSW harmonic edge. The time interval of our numerical
    computations varies depending on u_+. For instance, when we take
    u_+ = 0.1, we approximate the solution at t ≈ 1000, while
    taking u_+ = 0.9 requires the longer time t ≈ 4000 for the
    solution to asymptotically develop. Upon varying u_+, we observe
    good agreement between the predictions of DSW fitting for the DSW
    harmonic and soliton edge velocities and numerical simulations for
    u_+ ≳ 0.5, while the predictions begin to deviate from what is
    observed in numerical simulations below this threshold.  The DSW
    fitting method is subject to the convexity conditions
    
        ∂ c_±/∂ u_- 0,   ∂
            c_±/∂ u_+ 0 .
    
    A direct calculation for the potential Φ(u̅) = u̅^3/3
    demonstrates that these convexity conditions are indeed satisfied.
    The numerical results suggest that the DSW fitting method provides an
    adequate approximate prediction of the discrete DSW edge properties.
    
    
    
    The harmonic edge of the DSW is accompanied by linear radiation much
    like the left edge of the RW in the previous section. To describe
    this, we apply the same approach as in Sec. <ref> to
    approximating the small amplitude linear waves that emanate from the
    DSW's harmonic edge. The only change is that, for the DSW,
    u_+ - u_- < 0 in Eq. (<ref>). A comparison is shown in
    Fig. <ref>.
    
    
    
      §.§.§ Approximations of the leading edge amplitude via a
      quasi-continuum model
    
    
    
    In this subsection, we go a bit deeper into the description of the
    solitary wave at the soliton edge of the DSW.  Like in the previous
    subsection, for u_+ ∈ (0,1) (where u_-=1 is fixed), we
    numerically solve the Riemann problem, generate a DSW and extract the
    amplitude of the soliton edge, and its speed. This is done by
    inspecting the time series of a node sufficiently far from the center
    of the lattice (we chose n=300, in which case we have observed the leading edge is developed) 
    and
    simply computing the amplitude as a_+ = maxu_n(t) - u_+. The
    speed is estimated by computing c_+ = 1/(t_n+1 - t_n) where t_n
    and t_n+1 are the time values where u_n and u_n+1 attain
    their maxima.  The blue open circles in
    Fig. <ref>(a) show the amplitude of the DSW soliton
    edge as a function of u_+. Since for each value of u_+, we compute
    both the speed c_+ and amplitude a_+, we also show a parametric
    plot of (c_+, a_+) parameterized by u_+ in
    Fig. <ref>(b).
    
    To confirm that the DSW soliton edge is indeed described by a solitary
    wave, we compute a numerical solitary wave solution of the lattice
    equation (<ref>) by using a fixed-point iteration scheme
    <cit.> to solve the advance delay equation
    
    
        2 c 𝒱^'(η)=Φ^'(u_+ + 𝒱(η+1))-Φ^'(u_+ + 𝒱(η-1)),
    
    which is a rescaled copy of (<ref>) and obtained by
    substituting
    
        u_n(t)=u_+ + 𝒱(η),     η = n-  c t
    
    into Eq. (<ref>).  While we are free to select values of c and
    u_+ when solving Eq. (<ref>), we select combinations of them
    according to the relationship extracted from the DSW soliton edge
    (i.e., the blue circles in Fig. <ref>(a)). Upon
    convergence of the scheme, we compute the amplitude of the resulting
    solitary wave, which is the maximum of the wave minus the background
    u_+.  The amplitudes of the solitary waves are shown as the solid
    red dots in Fig. <ref>(a,b). Note that these red
    dots fall almost exactly within the blue circles, indicating that the
    soliton edge of the DSW is indeed well-approximated by a solitary wave
    solution of the lattice equation.
    
    We can obtain an analytical approximation of discrete solitary waves
    by considering the BBM quasi-continuum approximation (<ref>)
    of the lattice dynamics (<ref>).
    Entering the moving frame U(X,T) = ϕ(X - c T) and integrating
    once, this PDE becomes the second order ODE
    
    
        c ϵ^2/6ϕ” = B +  c ϕ -  2ϕ^2
    
    where B is an arbitrary integration constant. This ODE can be solved
    using quadrature <cit.>.  In particular, the solitary wave
    with tails decaying to the background state u_+ has the form
    
    
        u_n(t) = ϕ(X - c T) = u_+ +
          (3c/2 - 3 u_+ )^2(√(3/2c- 3 u_+/c)  (n - c t))
    
    which assumes c > 2u_+. Note the maximum speed of linear waves on a
    background u_+ is 2u_+, implying the solitary waves travel faster
    than all linear waves, as expected.  In Eq. (<ref>), c and u_+
    can be chosen independently of each other but we once again select
    combinations of them according to the relationship extracting from the
    DSW soliton edge (the blue circles in
    Fig. <ref>(a)). The solid blue dots of
    Figure <ref>(a) show the quasi-continuum prediction
    of the amplitude a_+ = 3c_+/2 - 3 u_+, where it is seen that
    the approximation becomes better as the jump height decreases (i.e.,
    as u_+ → 1).  The quasi-continuum prediction of the
    amplitude is only semi-analytical as it relies on the numerically
    obtained relationship of u_+ and c_+ from the DSW soliton edge
    data. An analytical prediction can be derived by using the DSW soliton
    edge speed in Eq. (<ref>) of the previous subsection, which
    allows us to express the amplitude, a_+, of the DSW soliton edge in
    terms of just u_+ (or c_+):
    
        a_+ = 6√(1-u_+)/cosh^-1(2-u_+/u_+) - 3 u_+
    
    See the solid red line of Fig. <ref>(a,b) for a plot of
    this formula.
    
    Finally, the solitary wave profile given by Eq. (<ref>) matches the
    numerically computed solitary wave solution of Eq. (<ref>) quite
    well, especially for longer wavelength solutions.  See
    Fig. <ref>(c) for a comparison of the actual
    solitary wave (solid red dots) and quasi-continuum approximation
    (solid red line) for two example parameter sets.
    
    
      
    
    
    Because the quasi-continuum approximation of DSWs here and in
    <cit.> performs well, we briefly contrast the Riemann
    problems that result in DSWs for the lattice (<ref>) and BBM
    (<ref>) equations.  The DSW fitting technique was applied to
    the BBM Riemann problem in <cit.>.  In order to
    compare our results for the lattice with DSWs in the BBM
    equation (<ref>), we consider the initial transition 
    occurring
    between u_- = 1 and u_+ = 1-Δ for 0 < Δ≪ 1.  At the
    DSW harmonic edge, the characteristic wavenumber
    ϵ^2 K_-^2 = 4Δ + 26 Δ^2/9 + ⋯ and speed
    C_- = 2 - 4 Δ + 14 Δ^2/9 + ⋯ for BBM agree to
    𝒪(Δ) with the expansions
    k_-^2 = 4 Δ + 4/3Δ^2 + ⋯ and
    c_- = 2 - 4 Δ for the lattice.  Similarly, at the DSW soliton
    edge, the conjugate wavenumber
    ϵ^2 K̃_+^2 = 4 Δ + 10 Δ^2/9 + ⋯ and
    speed C_+ = 2 - 2 Δ/3 - 2 Δ^2/27 + ⋯ agree to
    𝒪(Δ) with the lattice:
    k̃_+^2 = 4 Δ + 8 Δ^2/3 + ⋯ and
    c_+ = 2 - 2 Δ/3 - 8 Δ^2/45 + ⋯.  The DSW's soliton
    edge amplitude in BBM is A_+ = 2 Δ - Δ^2/9 + ⋯
    whereas the prediction (<ref>) for the lattice expands
    as a_+ = 2 Δ - 4 Δ^2/15+⋯.  Note that to leading
    order, these predictions agree with the DSW edge characteristics of
    the KdV equation U_T + 2 UU_X + 13 ϵ^2 U_XXX = 0 for
    X = ϵ n, T = ϵ t, U(X,T) = u_n(t).  This is
    expected because BBM and KdV are asymptotically equivalent to leading
    order in the weakly nonlinear, long wavelength regime.
    
    
    In summary, the quasi-continuum approximation of lattice DSWs by DSWs
    in the BBM equation performs well for small initial jumps
    0 < Δ≪ 1 when the oscillation wavelengths are much larger
    than the lattice spacing.  The agreement to 𝒪(Δ) in
    the DSW properties is expected and, in fact, is a statement of
    universality of the KdV equation 
    
    
    as a weakly nonlinear, long wavelength model of dispersive
    hydrodynamics <cit.>.
    For sufficiently large Δ, the
    BBM DSW develops two-phase modulations near the trailing edge
    <cit.>.  In contrast, for Δ > 1, the
    lattice DSW bifurcates into a partial DSW connected to a traveling
    wave called a TDSW or exhibits blow up that we will describe in
    Secs. <ref> and <ref>, respectively.
    
    
    
    
    
    
     §.§ Dispersive shock wave  + stationary shock + rarefaction (DSW+SS+RW) 
    
    
    In this section we investigate the case where the initial step
    generates two unsteady waves: a leftward moving DSW and a right moving
    RW.  At the origin, there is a stationary shock (SS) joining symmetric
    states at the level u_n = ± u_0. A numerically computed example is
    depicted in Fig. <ref>. This class of solution is
    empirically found for u_+= 1 and u_-∈(-1,-0.26).  As will be
    shown, the bifurcation at u_- = -0.26 occurs when the DSW's harmonic
    edge exhibits zero velocity.
    
    
    
    The numerical simulation shown in Fig. <ref> suggests that
    the solution can be approximated for large t as:
    
        u_n(t) = 
                     u_-     n ≤ s_- t 
         
                     u_ DSW(n,t)     s_- t < n ≤ s_+ t 
         
                     -u_0     s_+ t < n ≤ 0
        
                     u_0     0  < n ≤ 2u_0 t
        
                     u_ RW(n,t)     2u_0 t < n  ≤ 2u_+ t
        
                     u_+     u_+ t < n
                    .
    
    The velocities s_- < s_+ give the motion of the DSW's soliton and
    harmonic edges, respectively.  Across all of the simulations
    performed, we found the following relation for the intermediate,
    symmetric states ± u_0 to hold to very high precision
    
        u_0 = u_+ - u_-/2 .
    
    This relation implies that the DSW and RW have the same jump height,
    albeit with opposite polarities.  We have been unable to
    mathematically justify this formula.  However, as we show in section
    <ref>, the value of the intermediate state u_0
    depends strongly on particular details of the initial data.  If the value
    of u_0(0) is changed, then the intermediate value u_0 differs from
    (<ref>).
    
    Utilizing the formula (<ref>), we can completely determine the
    velocities that divide the approximate solution (<ref>) into
    different wave patterns.  To determine the DSW edge velocities, we use
    Eqs. (<ref>) and (<ref>), which were derived under the
    assumption that u_- > u_+ > 0.  In the case of the solution
    (<ref>), the left (u_-) and right (-u_0) states are both
    negative.  Since the governing equation (<ref>) is invariant
    under the transformation u_n(t) → -u_-n(t), the DSW velocities
    are mapped as follows
    
        s_-(u_-,-u_0) = -c_+(u_0,-u_-),    s_+(u_-,-u_0) = -c_-(u_0,-u_-) .
    
    Then, using (<ref>), we find
    
        s_- = - 2/cosh^-1(u_+/|u_-|)√(u_+^2 - u_-^2)
          ,    s_+ = 3 u_- + u_+ .
    
    Figure <ref>(b) shows good agreement between the
      predicted velocities of the approximate solution (<ref>) and a
      numerical simulation when u_+ = 1, u_- = -0.5.
    
    The DSW remains detached from the stationary shock (SS) so long as the
    harmonic edge velocity s_+ remains negative.  From
    Eq. (<ref>), we predict that the DSW is no longer detached from
    the SS when 3 u_- + u_+ = 0, which, for u_+ = 1, occurs when
    u_- = -13.  As noted earlier, the bifurcation from the DSW +
    SS + RW to the US case is empirically identified as occurring when
    u_- = -0.26.  As shown in Fig. <ref>(b), this small
    discrepancy in the bifurcation value can be explained by the deviation
    of the computed DSW harmonic edge velocity from the DSW fitting
    prediction (<ref>).
    
    
    
    
     §.§ Traveling dispersive shock wave (TDSW)
    
    
    In this section, we consider the case where a partial DSW connects the
    level behind, u_-, to a periodic-to-equilibrium traveling wave
    solution to the level ahead u_+.  Although we do not directly
    compute it as a traveling wave solution of the discrete equation, it
    is interpreted as a heteroclinic connection between a periodic orbit
    with the constant level ahead u_+ based on an analysis of the
    numerical simulations of the Riemann problem.  Such heteroclinic
    solutions of continuum equations with higher (fifth) order dispersion
    were studied in <cit.> and were
    associated with so-called traveling dispersive shock waves (TDSWs)
    that emerge from an associated Riemann problem.
    
    For u_- = 1 and u_+ ∈ (-0.724,0), numerical simulations show a
    qualitatively similar solution pattern to that depicted in
    Fig. <ref> in which u_n(t) begins on the left with the
    value u_-.  It then progresses into an oscillatory wavetrain with
    increasing amplitude that resembles the leftmost portion of a DSW,
    called a partial DSW, that is then connected to a periodic traveling
    wave.  The periodic traveling wave is connected to the constant state
    ahead u_+ via an abrupt transition that moves with the same speed.
    Collectively, this partial DSW and traveling wave is referred
    to as the
    traveling DSW.
    
    
    
    The TDSW solution studied here is the discrete analogue of the TDSW
    studied for the KdV5 equation <cit.>. The terminology
    traveling dispersive shock wave unites the unsteady component of the partial DSW and the steady traveling wave component (a heteroclinic periodic-to-equilibrium solution) to which it is attached in this non-classical DSW.  Consequently, the entire TDSW structure is unsteady.
    
    
    
    
      §.§.§ Approximation of the traveling wave via the quasi-continuum model
    
    
    One can describe the periodic portion of the TW (see e.g. the interval
    n∈[-90,100]) of Fig. <ref>) using the continuum
    reduction presented in section <ref>. In particular, there is a
    three-parameter family of traveling wave solutions of the
    quasi-continuum BBM Eq. (<ref>), given by
    
    
        u_n(t) = r_1 + r_2 - r_3 + 2(r_3 - r_1) ^2( √(2(r_3 - r_1)/c)( n - c t),m) ,   m = r_2-r_1/r_3-r_1,    c = 2(r_1 + r_2 + r_3)/3.
    
    We treat r_1,r_2,r_3 as fitting parameters. After a sufficiently
    long time, the traveling wave in the numerical simulation forms, as in
    Fig. <ref>. We then isolate a small interval of that
    traveling wave, which is then fit to
    Eq. (<ref>). Figure <ref>(a) shows a
    comparison of the actual lattice dynamics at t=480 (blue markers)
    and quasi-continuum approximation (blue lines) with the step values
    u_-=1 and u_+=-0.16.
    
    
    Figure <ref>(b) shows the trajectory in the phase
    plane (u_400(t),u̇_400(t)) (two outermost lobes) and
    (u_140(t),u̇_140(t)) (two innermost lobes) for various
    values of u_- for the actual lattice dynamics (markers) and
    quasi-continuum approximation (lines). We show the phase plane for
    different values of n since the location of the traveling wave
    within the lattice is moving.
    
    The agreement is quite good throughout the interval of existence for
    these structures, but is best when the jump height is smallest
    (compare the blue and red trajectory of
    Figure <ref>(b)). The comparison of the
    frequency, amplitude, and mean parameters is shown in
    Figure <ref>(c).  These are computed via the
    following formulas with n fixed.  For u_-∈[-0.72,-0.28] the
    lattice index is fixed to n=140, for u_-∈(-0.28,-0.2] the index
    is n=400 and for u_-=0.08 the index is n=500.  The frequency is
    f = 1/T, where T is the period (computed as the peak-to-peak time
    of the trajectory); the mean is
    
        u = 1/T∫_I_T u_n(t)   dt
     where I_T is the time
    interval of one oscillation period. For the computations shown here,
    it is I_T = [480-T,480]. The amplitude is
    
        a = max_t∈ I_Tu_n(t) - min_t∈ I_T u_n(t).
    
    Once the best-fit values of r_1,r_2,r_3 are obtained, the wave
    parameters can be computed directly from Eq. (<ref>) as
    
        f = √(2 c (r_3 - r_1))/2K(m),   u = r_1 +
        r_2 - r_3 + 2(r_3-r_1) E(m)/K(m),    a=2(r_2 - r_1),
    
    where K(m) and E(m) are the complete elliptic integrals of the
    first and second kind, respectively. We note that while
    Eq. (<ref>) is able to describe the local periodic traveling
    wave dynamics of the TDSW structure, it does not admit solutions
    resembling the entire TDSW structure since heteroclinic
    periodic-equilibrium solutions do not exist for the planar ODE
    (<ref>). Such a description may be possible by using the (1,5)
    Padé approximant instead of the (1,3) approximant in (<ref>)
    to arrive at the 5th order model
    
        U_T + (U^2)_X - ϵ^2/6 U_XXT + 7
          ϵ^4/360 U_XXXXT = 0.
     Similar models have been shown to
    admit such solutions <cit.>. While this is an interesting
    topic for further study, we will not pursue the identification of such
    a heteroclinic orbit further herein.
    
    
    
    
      §.§.§ Modulation solution of the weakly nonlinear Whitham modulation equations
    
    
    
    In order to obtain the form of an approximate modulation solution
    𝐪 = [u̅,a,k]^T of the Whitham equations
    (<ref>) that describes the TDSW, we appeal to the
    structure of the TDSW evident in Fig. <ref>. An
    oscillatory wavetrain emerges from the left level u_- with
    increasing amplitude that saturates at a periodic traveling wave.  The
    traveling wave then abruptly transitions to the right level u_+.
    Guided by previous work on the traveling dispersive shock wave (TDSW)
    solutions of a fifth-order Korteweg-de Vries equation
    <cit.>, we make the self-similar modulation ansatz
    (ξ = X/T = n/t)
    
        𝐪(ξ) =
          𝐪_-    ξ < λ_2(𝐪_-) , 
        𝐪_ RW(ξ)
                            λ_2(𝐪_-)  ≤ξ <
                           λ_2(𝐪_ p) , 
        𝐪_ p   λ_2(𝐪_ p) ≤ξ <
                          λ_2(𝐪_+) , 
        𝐪_+
                            λ_2(𝐪_+)  ≤ξ ,
    
    for Eq. (<ref>) where λ_2 is the middle characteristic
    velocity.  Additionally, the constant states are
    
        𝐪_- = [u_-,0,k_-]^T,   𝐪_ p =
          [u_ p,a_ p,k_ p]^T,   𝐪_+ = [u_+,a_+,0]^T,
    
    and 𝐪_ RW(ξ) is the rarefaction solution (integral
    curve) of Eq. (<ref>) for the second characteristic field
    (λ_2,𝐫_2) that continuously connects the harmonic
    edge state 𝐪_- and the periodic traveling wave state
    𝐪_ p of the TDSW.  The discontinuity from
    𝐪_ p to 𝐪_+ satisfies the jump conditions
    (<ref>) for the Whitham modulation equations where
    V = λ_2(𝐪_+) is simultaneously the phase speed of the
    periodic traveling wave with parameters 𝐪_ p, the
    shock speed, and the phase speed of the solitary wave with parameters
    𝐪_+, i.e., it represents the TW component of the TDSW
    solution.  The modulation solution (<ref>) corresponds to a
    rarefaction-shock solution of the Whitham modulation equations.
    
    The five parameters (k_-,𝐪_ p,a_+) in the modulation
    solution (<ref>) could, in principle, be obtained by solving the
    full Whitham modulation equations (<ref>), but we lack
    explicit periodic traveling wave solutions.  Instead, we approximate
    the modulation solution (<ref>) in the weakly nonlinear regime
    by solving Eqs. (<ref>).
    
    
    
    
    A general feature of Whitham modulation systems in the weakly
    nonlinear regime is the 𝒪(a^2) mean induced by the finite
    amplitude modulated wavetrain <cit.>.  Absent mean changes
    due to initial/boundary data, the third order weakly nonlinear
    modulation system (<ref>) can be simplified, with the
    same order of accuracy, to a second order modulation system.  The
    procedure to do so is as follows.  The induced mean is represented by
    the ansatz
    
        u(X,T) = u_0 + a(X,T)^2 u_2(k(X,T)) + ⋯,
    
    where u_0 ∈ is a constant background.  This introduces the
    induced mean coefficient u_2(k) that gives rise to an effective
    nonlinear frequency shift ω̃_2 by inserting
    (<ref>) into Eq. (<ref>) and expanding as
    
        ω(k,u̅)    = ω_0(k,u_0) + a^2 ( ω_2(k,u_0) +
                                ∂_uω_0(k,u_0) u_2(k)  ) +
                                o(a^2), 
           ≡ω_0(k,u_0) + a^2 ω̃_2(k,u_0) +
                                o(a^2) .
    
    With this effective nonlinear frequency shift, weakly nonlinear wave
    modulations are generically described by the simplified system
    <cit.>
        (a^2  )_t +  ( ω_0,k a^2  )_x    = 0 ,
        
            k_t +   ( ω_0  )_x +  ω̃_2  (a^2
             )_x     = 0 ,
    
    provided
    
        u_2(k) = -1/16 u_0(1-cos k),   ω̃_2(k) =
          (cosk-2) (k/2)/16 u_0 .
    
    Note that the additional coupling term involving ω̃_2
    and its derivatives contribute at higher order in a.  The induced
    mean coefficient u_2(k) in Eq. (<ref>) is determined by
    compatibility of averaged mean, energy conservation laws
    (<ref>), (<ref>) with the
    induced-mean modulation system (<ref>), which asymptotically
    satisfies (F(k) a^2)_T + (F(k)ω_0,k a^2)_X = 0 for any
    differentiable F, in particular F(k) = u_2(k).
    
    Under the assumption of induced mean variation, we can analytically
    obtain the rarefaction solution 𝐪_ RW(ξ) in
    (<ref>) by solving Eq. (<ref>) for a RW and then
    inserting it into Eq. (<ref>).  For this, we express the
    induced-mean modulation system in Riemann invariant form
    
        ∂ r_±/∂ t + λ_±∂
          r_±/∂ x = 0,
    
    where 
    
        r_±   = a ∓1/2∫(
                    ω_0,kk/ω̃_2)^1/2 dk  
           = a ∓ 2√(2)  u_0 cos^-1(13(-1+2cos
                    k)),    u_0 >
                    0,     0 < k <
                    π,    a ≥ 0.
    
    The restriction to positive mean and wavenumber is due to the fact
    that we have selected the positive square root in (<ref>).  The
    characteristic velocities are
    
        λ_± = 2 u_0 cos k ±a/2cos ( k/2 )
          √(2 - cos k),
    
    so that λ_+ = λ_2 and λ_- = λ_1 in
    Eq. (<ref>).
    
    We can now solve for 𝐪_ RW(ξ) in (<ref>) by
    setting u_0 = u_- and taking the fast RW solution of
    Eq. (<ref>) that satisfies r_- = const and λ_+ = ξ.
    The constant slow Riemann invariant r_- implies
    
        cos^-1(13(-1 + 2cos k_-)) 
            = a_ p/ 4√(2) u_- +
            cos^-1(13(-1 + 2cos k_ p))  ,
    
      and the assumption of induced mean implies
      
        u_ p = u_- + u_2(k_ p) a_ p^2 .
    
      The RW profile is obtained by inverting λ_+ = ξ for
      a(ξ) and k(ξ) subject to the constraint
      
        cos^-1(13(-1 + 2cos k_-)) 
            = a(ξ)/ 4√(2) u_- +cos^-1(13(-1 +
              2cos k(ξ))) .
    
      Equations (<ref>) and (<ref>) are two conditions on the
      four unknown solution parameters
      (k_-,u_ p, a_ p, k_ p).  The other two
      conditions are obtained from the jump conditions (<ref>).
    
    
    
      The sharp transition from the periodic traveling wave
      𝐪_ p to the solitary wave ahead 𝐪_+ is
      achieved by a shock solution of the Whitham modulation equations.
      We obtain the jump conditions from the conservation laws
      (<ref>) and (<ref>) by assuming that the periodic
      traveling wave is in the weakly nonlinear regime (<ref>)
      and the level ahead is a solitary wave where k → 0, both with
      the same phase speed V:
      
        -V(u_ p - u_+) +  u_ p^2 +
            1/8a^2_ p - u_+^2
                = 0 
         
            -V(1/3u_ p^3 +
            1/8u_ p 
            a^2_ p - 1/3u_+^3) + 1/2u_ p^4
            + 1/4a_ p^2(u_ p^2cos(k_ p) + 1/2 u_
            p^2) - 1/2u_+^4 
                = 0 .
    
      The jump condition from the conservation of waves equation
      (<ref>) is satisfied because ω = k = 0 for the solitary
      wave and V = ω_ p/k_ p for the periodic traveling
      wave
      
        V = ω_0(k_ p,u_-) + a_ p^2 ω̃_2(k_
                p)/k_ p .
    
    
    Solving for a_ p and u_ p from (<ref>)
    and (<ref>), then inserting them into (<ref>),
    (<ref>) and using the phase velocity (<ref>)
    determines two nonlinear equations for k_- and k_ p.  We
    solve these equations numerically using standard root finding methods
    to obtain all the parameters of the shock-rarefaction modulation
    solution (<ref>) and compare it with numerical simulation in
    Figs. <ref> and <ref>.
    
    
    
    We observe in the figures that near the onset of the TDSW, the
    solution's mean, amplitude and frequency are accurately captured by
    the above theory. On the other hand, as the amplitude of the solution
    increases, the approximation loses quantitative efficacy.
    Nevertheless, the qualitative trend of the solution's properties are
    captured by the above analysis.
    
    
    
    
    
    
    
     §.§ Unsteady Shock (US)
    
    
    In section <ref>, we found that for u_+ = 1 and
    u_- ∈ (-1,-0.26), a stationary shock separated two counterpropagating
    waves, one a DSW, the other a RW.  When u_- = -0.26, the DSW no
    longer separates from the stationary shock.  Instead, the DSW+SS+RW is
    numerically observed to bifurcate into the unsteady generation of
    counterpropagating periodic waves that we term an unsteady shock (US)
    for u_- ∈ (-0.26,0.18).  When u_- exceeds 0.18, the US
    bifurcates into a RW, described in section <ref>.
    We now investigate the US.
    
    A plot of the solution for u_- = 0
    and u_+ = 1 is given in
    Figure <ref>. Two distinct counterpropagating
    periodic waves traveling with speed c_± emerge from the origin that then transition to the
    constant level u_- behind through a partial DSW and to u_+ ahead
    via a partial DSW and a RW. 
    
    
    
    While the velocities c_± are distinct, Figure
    <ref>(b) depicting the time series
    u_± 20(t) indicates that the two periodic waves have
    approximately the same temporal frequency, which we confirm below
    numerically to high precision.
    
    
    
    
      §.§.§ Poincaré description
    
    
    A spatio-temporal intensity plot of the US with u_- = -0.16 near n = 0 is shown
    in Fig. <ref>(a).  The
    counterpropagating traveling waves can clearly be identified. While
    the underlying wave parameters of the left- and right-moving waves
    will generally be different, they do share the same frequency.  Thus,
    the dynamics correspond to a time-periodic solution. Indeed,
    inspection of Fig. <ref>(c) confirms
    this, which shows the same intensity plot as panel (a), but with the
    solution sampled every T time units, where T is the period of
    oscillation. This is the Poincaré map of the dynamics.  With this
    sampling size, the solution appears to be constant, suggesting that
    the waveform is genuinely time-periodic. To demonstrate this further,
    we simulated the equations of motion on a small lattice n∈[-20,20]
    with boundary conditions given by the periodic solution, i.e., the
    left boundary is given by u_-20(t) and the right boundary is given
    by u_20(t). The dynamics upon initialization with the periodic
    solution u_n(450) are shown in
    Fig. <ref>(b), which can be hardly
    distinguished from the dynamics in (a).  The evolution remains
    periodic, as can be inferred from
    Fig. <ref>(d), which are the dynamics
    sampled every T seconds.  An avenue for potential further study,
    prompted by these findings, is the seeking of exact time-periodic
    solutions of the model and their corresponding Floquet analysis.
    
    
    
      §.§.§ Approximation of traveling waves via quasi-continuum model
    
    If considering the left-moving and right-moving waves as separate entities, we can once again apply the quasi-continuum
    reduction to describe the traveling wave using formula Eq. (<ref>). A comparison of the spatial profile
    of the left wave moving wave with u_+ = -0.16 at time t=480 is shown in Fig. <ref>(a).
    The phase plane for the left-moving waves (left lobes) and right-moving waves (right lobes) is shown in
    Fig. <ref>(b) in markers, with corresponding quasi-continuum approximations shown
    as solid lines. Like before, the reduction is best for smaller step heights (compare the blue and red 
    orbits in panel (b)). A comparison of the mean, amplitude and frequency of the simulation and quasi-continuum
    approximation as a function of u_- is shown in panel (c) for both the left- and right-moving waves. We can observe that the approximation provides a very
    adequate description of the relevant traveling patterns.
    
    
    
    
    
    
    
    
      §.§.§ Jump Conditions
    
    
    Figure <ref> shows two
    counterpropagating periodic traveling waves with a rapid transition
    between them in the vicinity of n = 0.  This motivates the
    hypothesis that these two waves satisfy the jump conditions obtained
    from the Whitham modulation system's conservation laws.  We denote the
    periodic traveling waves by φ_± for the left (-) and right
    (+) periodic waves, respectively. For a discontinuous, shock
    solution of the Whitham system at the origin, the corresponding jump
    conditions are Eq. (<ref>)
        ⟨Φ'(φ_-) ⟩ -  ⟨Φ'(φ_+) ⟩    = 0  
        ⟨Φ'(φ_-)𝒮Φ'(φ_-)
            ⟩ -  ⟨Φ'(φ_+)𝒮Φ'(φ_+)⟩    = 0  
        ω_- - ω_+     = 0,
    
    where S is the unit shift operator SR(η) = R(η + 1). To
    check if these jump conditions are indeed satisfied, we approximate
    the above averages using the numerical simulations. In particular, we
    let the structure come close to a periodic state (as in
    Fig. <ref>) and extract one period of
    motion at a particular node n. Let T_n be the period of node n
    and let I_T_n = [τ,τ+T_n] be the corresponding time interval from
    t = τ.  We then make the following approximations
    
    
        ⟨Φ'(φ) ⟩   ≈1/T_n∫_I_T_nΦ'(u_n(t)) dt =: f(n), 
        ⟨Φ'(φ_-)𝒮Φ'(φ_-)
            ⟩   ≈1/T_n∫_I_T_nΦ'(u_n(t))
                            Φ'(u_n+1(t)) dt =: g(n) ,
    
    for τ≫ 1 (we set τ = 240 in what follows).
    
    The first jump condition, Eq. (<ref>), is checked by
    comparing f(-5) to f(5), while the second jump condition
    Eq. (<ref>) is checked by comparing g(-5) with g(5).  The
    third jump condition is simply the difference in the frequency, which
    we estimate via 1/T_-5 and 1/T_5. Figure <ref> shows a
    plot of f(-5) (open blue circles), g(-5) (open blue squares) and
    1/T_-5 (open blue triangles), while the red dots show the
    corresponding quantities for n=5. Notice that each red point falls
    nicely into an open blue marker, indicating that the jump conditions
    are, up to some small numerical error, satisfied. The maximum residual
    over the interval of u_- values tested for the first jump condition
    was max_u_+ | f(-5) - f(5) | ≈ 0.001,
    whereas the maximum residual for the second condition was 0.0017 and
    the maximum residual for the third was 0.0008.
    
    The numerical evidence is a compelling indication that the US can be
    interpreted as a shock solution of the Whitham modulation equations.
    While shock solutions of the Whitham equations have been constructed
    previously <cit.>, their
    admissibility requires the existence of traveling wave solutions of
    the corresponding continuum PDE in which the phase velocities and
    shock velocity all coincide.  In the present case of the US for the
    lattice equation (<ref>), all three of these velocities differ
    but the frequencies are the same.  This suggests a new class of
    admissible shock solutions to the Whitham equations corresponding to
    time-periodic solutions of the lattice equation, an intriguing
    possibility for future work.
    
    
    
    
    
    
    We also note the similarity between the US and defect solutions of
    reaction-diffusion equations
    
    
    <cit.>. Because the waves in the US are in-phase (see
    Fig. <ref>), it most closely resembles a
    target pattern with a source from which waves are emanating, in the
    language of <cit.>.  The target pattern exhibits a Hopf
    bifurcation of the background state and a specific transition of the
    eigenvalues associated with the linearized operator about the
    background state.  It would be interesting to explore potential
    connections between the underlying diffusive regularization of the
    target pattern and the dispersive regularization of the US studied
    here.
    
    
    
     §.§ Blow up
    
    
    
    In section <ref>, we observed that a TDSW (a partial DSW
    connected to a traveling wave) is generated that transitions from the
    left level u_- = 1 to the right level u_+ ∈ (-0.724,0). By
    decreasing the value of u_+ below -0.724, a new wavetrain develops
    on the right level. Figure <ref> shows three example
    profiles. When u_+ is close to the transition value u_+=-0.724,
    the excitation on top of the right level u_+ is small, see Figure
    <ref>(a).  Decreasing u_+ further results in a
    larger amplitude wavetrain that resembles another TDSW; see
    Fig. <ref>(b).  Close to u_+ =- 1, the wavetrains
    on the left and right levels approach binary oscillations, as shown in
    Fig. <ref>(c).  Solution dynamics similar to those
    shown in Fig. <ref>(c), but with odd initial data,
    were studied extensively in <cit.>.
    In the work of <cit.>, it is claimed that the emergence of blow
    up for the odd initial data they considered is “almost always”
    associated with the emergence of regions of binary oscillations for
    which the upper and lower oscillatory envelopes have opposite sign.
    Such a scenario was observed to result in the loss of hyperbolicity in
    the modulation equations (<ref>) for binary oscillations and,
    consequently, a dynamical instability and thus exponential growth in a
    localized region of space that ultimately led to the finite time
    blow up of the wave pattern.  For the Riemann data considered here
    (<ref>), we numerically observe blow up in regions of the
    solution where binary oscillations are not apparent, which we now explore.
    
    
    We will start with a more detailed discussion of the structure in
    Fig. <ref>(b) with u_+ = -0.8, which is
    representative of many of the patterns found for u_+ ∈ (-1,-0.724)
    and u_- = 1. A zoom-in of the solution in the co-moving frame n-ct
    near the shock interface is shown in Fig. <ref> at three
    separate times (t ∈{340,345,350}) represented by different
    colors.  Both the solution on the lattice (dots) and its zero-padded
    Fourier interpolant (curves) are shown with two different speeds.  In
    Fig. <ref>(a) with speed c = c_1 = 0.229, the leftmost
    wave at the three distinct times overlap, suggesting that it moves in
    the steady frame n-c_1t.  Contrastingly, in
    Fig. <ref>(b) with speed c = c_2 = 0.15, the rightmost
    wave at the three distinct times overlap, suggesting it moves in the
    slower steady frame n-c_2 t.  
    
    
    
    
    The spatial profile at t = 350 and the Fourier transform
    û(k) of a 40-site window of the leftmost (rightmost) wave are
    shown in Fig. <ref>(a) in red (blue). Each wave has
    wavenumber concentration, and they are distinct (k≈ 2.76 for
    the left and k=2.51 for the right).  This solution is a
    generalization of the unsteady shock (US) studied in
    section <ref> in that two traveling waves are connected
    through a sudden jump. However, there are key differences.  First, the
    shock interface itself is moving as shown in Fig. <ref>.
    
    
    
    
    Second, the frequency for the leftmost and rightmost waves are not
    identical, which can be seen in Fig. <ref>(c) that shows
    the time series at a node located at the left wave (n=20) and the
    right wave (n=250).  Here, it is clear that the oscillation
    frequencies are distinct.  We conjecture that, much like the US, this
    wavetrain can be modeled by a discontinuous, weak solution of the
    Whitham equations satisfying the jump conditions (<ref>).
    Rather than pursue this further, we instead turn to an examination of
    the large t dynamics of the solution and, specifically, its eventual
    blow up.
    
    This structure is relatively coherent until about t=600, after which
    small disturbances in the leftmost wave develop. Disturbances are
    noticeable if one compares the time series for t∈(350,400) (shown
    in Fig. <ref>(c)) and for t∈(750,800) (shown in
    Fig. <ref>(d)). At about t=816 the solution appears to
    experience finite-time blow up. Figure <ref>(a) shows
    the blow up. The spatial profile close to (but before) the time of
    blow up is shown in Fig. <ref>(b). Notice the location of
    the blow up is spatially concentrated within the leftmost wave and that the
    wavenumber of the traveling wave where the instability seems to
    manifest is about k≈ 2.76, (i.e., not a binary oscillation).
    
    
    We conjecture that the observed blow up is due to an instability of
    the leftmost wave with wavenumber k≈ 2.76, and not due to
    numerical instability. A piece of evidence in this direction is that
    the mass and energy are conserved until times very close to the blow
    up time.  The solid gray and dashed lines of
    Figs. <ref>(c,d) show the mass M(t) and energy E(t),
    respectively. Note, in order for these quantities to be conserved, we
    employ periodic boundary conditions.  
    For the simulations in this
    subsection, we concatenate the initial condition Eq. (<ref>) with
    its reflection about the first site, leaving us with 2N total
    nodes. Thus, the relevant window of space for the plots is only the
    second half of the lattice. We define the lattice indices so that the
    initial (t=0) jump from u_- down to u_+ occurs at n=0. This
    makes the plots consistent with those in the previous sections.  We
    include the entire spatial window for the computation of E(t) and
    M(t). Note that in Figs. <ref>(c,d) the quantities M(t)
    and E(t) are indeed conserved, even as the waveform begins to break
    down, see t≈ 790 of panel Fig. <ref>(d). The gray
    solid and dashed lines of Fig. <ref>(a) show plots of
    M(t) and E(t), respectively, for times leading to the blow up
    itself. While the energy E(t) remains constant after the strong
    onset of instability at about t=785, the energy conservation breaks
    down for t>800, while M(t) remains conserved. The conservation of
    mass in the numerical scheme is not surprising, since by direct
    computation one sees that the variational integrator applied to
    Eq. (<ref>) with Φ'(u) = u^2 conserves the mass
    exactly <cit.>. The near
    conservation of energy for the variational integrator relies on the boundedness of
    the underlying numerical solution <cit.>. This will be clearly
    violated for solutions exhibiting 
    collapse-type phenomena
    and thus it is not surprising
    that the energy is not conserved in the numerical scheme close to the time of blow up.
    Thus, it seems the initial collapse
    is due to instability of the wave (energy remains conserved for t<790), but after experiencing sustained large amplitude oscillations, the numerical scheme may begin to exhibit additional numerical instabilities (since energy is not conserved for t>790). The blow up time
    found here is an approximation that depends on the particulars of the numerical scheme.
    
    
    For u_+ ∈ (-1,-0.724) and u_- = 1 we observe a similar blow up of the solutions,
    with the blow up time varying roughly between t=700 and t=2500, see Fig. <ref>(b).
    We define the solution as blown-up once max_n |u_n| exceeds a large threshold. We practically used 1000 as the threshold
    (the actual threshold makes little difference in Fig. <ref>(b)
    since the blow up occurs so quickly). For this figure panel, we used a lattice size of 2N=100,000
    and simulated until t=10,000. We observe blow up with u_+ = -0.725 (at about t=1500) but no blow up
    with u_+ = -0.724, even when simulating until t=10,000. This sharp transition between finite-time blow up
    and no blow up is further evidence that the blow up is due to an underlying instability of the waveform
    and not due to  numerical instability.
    
    The question of stability of traveling waves of this lattice is thus
    an important open question, meriting further investigation.  Based on
    these findings, it appears that traveling waves with wavenumbers
    other than k=π
    can lead to instabilities and finite-time blow up, a
    generalization of the findings in <cit.> where binary
    oscillations with k = π were identified as a primary instability
    mechanism.
    
    
    
    
    § NON-UNIQUENESS OF RIEMANN PROBLEMS
    
    
    
    
    As already mentioned in the introduction, the solutions to the Whitham system (<ref>) are non-unique for certain classes of initial data. To discuss this in greater detail, we start with some preliminary comments concerning the inviscid Burgers' PDE (<ref>), which describes the dispersionless, hyperbolic scaling limit of the lattice (<ref>) in the case of no oscillations and is a subsystem of (<ref>) that governs the dynamics in the case of zero wave amplitude or zero wavenumber. To obtain an elementary example for non-uniqueness, we impose the odd initial condition
    
    
        u(x,0)=sign(x)
    
    and notice that the rarefaction wave
    
    
        u(x,t)=U(x/t)    with    
        U(ξ)={[      -1     for -∞<ξ<-2;     ξ/2     for -2<ξ<+2;      +1     for +2<ξ<+∞;         ].
    
    is the unique self-similar solution according to the classical theory of hyperbolic conservation laws, see for instance <cit.>. The latter complements the PDE with an additional selection rules to exclude solutions that are considered to be unphysical. The most prominent example is the Lax condition for shocks but without such an admissibility criterion there exists a plethora of possibilities to fulfill the initial value problem (<ref>)+(<ref>). For instance, the formulas
    
    
        U(ξ)={[        -1       for   -∞<ξ<-2;       ξ/2       for -2<ξ<-2 μ;        -μ       for  -2 μ<ξ<0;        +μ       for  0<ξ<+2 μ;       ξ/2       for +2 μ<ξ<+2;        +1       for   +2<ξ<+∞;           ].
        
            and    
        U(ξ)={[        -1       for -∞<ξ<-1-μ;        -μ       for  -1-μ<ξ<0;        +μ       for  0<ξ<+μ+1;        +1       for   +μ+1<+∞;           ].
    
    provide two families of further self-similar solutions in dependence of the real parameter 0<μ<1 and μ>1, respectively, and contain (<ref>) as limiting case for μ=0. The corresponding  profile functions U are illustrated in Figure <ref> and combine a steady discontinuity at x=0 (which violates the Lax condition on both sides) with either two rarefaction waves or two Lax shocks. 
    
    
    
    
    In numerical simulations of the lattice (<ref>) we find a related family of solutions that differ only in microscopic details of the imposed initial conditions. As a typical but still elementary example we study the lattice initial data 
    
    
        u_n(0)={[  -1 for n<0;   η for n=0;  +1 for n>0;     ].
    
    where the free parameter η reflects that the macroscopic Riemann data (<ref>) can be realized in many different ways on the microscopic scale. Figure <ref> reveals that the numerical solutions for different choices of η correspond to distinct waves on the macroscopic scale. We always find a steady discontinuity at x=0 but its jump height as well as the other waves depend crucially on the value of microscopic parameter. In analogy to Figure <ref> we observe two rarefaction waves in the case of 0<η<1 but for η>1  we have to replace the Lax shocks by dispersive shock waves.
    
    
    It remains a challenging task to understand the macroscopic impact of small-scale fluctuations in the initial data. For instance, it seems that the parameter η in (<ref>) determines the jump height of the emerging steady discontinuity in Figure <ref> but we are not aware of any heuristic or even rigorous explanation thereof. Numerical simulations also indicate that microscopic details might not be relevant for Riemann data that are either positive or negative on both sides but also this must be investigated more thoroughly. Both issues are also intimately related to the properties of the Whitham system (<ref>) which can be regarded as an extension of the Burgers' equation (<ref>). At least for sign changing Riemann problems,  the lattice (<ref>) is able to produce an entire family  of Whitham solutions and it is very natural to investigate and classify them in terms of selection criteria and entropy inequalities. Moreover, other nonlinearities might produce further effects due to linearly degenerated states in the scalar first order PDE corresponding to (<ref>); see <cit.> for a related problem in FPUT chains whose force functions is increasing with inflection point.
    
    
    
    
    § CONCLUSIONS AND FUTURE CHALLENGES
    
    
    In the present work, we revisit the first order nonlinear dynamical
    lattice of Eq. (<ref>) through the lens of “lattice hydrodynamics” by
    providing a systematic analysis of the solutions to the canonical
    lattice Riemann problem (<ref>) for the case of quadratic flux
    (<ref>). Building on earlier works
    of <cit.>, we have characterized
    dispersive shock waves and rarefaction waves, familiar from continuum
    dispersive hydrodynamics.  But we also discover a variety of
    non-classical hydrodynamic-like solutions.  In addition to finite time
    blow up, recognized earlier in the work of <cit.>,
    
    we identify three additional dynamical regimes that are interpreted
    using a combination of numerical simulation, quasi-continuum
    approximation, and Whitham modulation theory. These regimes include
    the generation of a counterpropagating DSW and RW pair separated by a
    stationary, abrupt shock. Additionally, an abrupt, unsteady transition
    between two counterpropagating periodic waves with the same frequency
    is interpreted as a shock solution of the Whitham modulation
    equations.
    
    
    
    Finally, the phase diagram <ref> of solutions to the
    lattice Riemann problem is rounded out by a traveling DSW (TDSW),
    consisting of an unsteady partial DSW, connected to a
    heteroclinic periodic-to-equilibrium traveling wave.
    
    These elements of the lattice model's phenomenology provide an
    opportunity to develop different aspects of the mathematical analysis
    of lattice hydrodynamics. For instance, we adapt the DSW fitting
    method of <cit.> to lattice equations in order to characterize
    the expansion and edge properties of lattice DSWs.  We also follow up
    the work of <cit.>, providing abstract modulation
    equations for a quasi-continuum analogue of the lattice dynamical
    system. We connect these to the genuinely discrete modulation
    equations, made explicit in the weakly nonlinear regime by a
    Poincaré-Lindstedt expansion of periodic traveling wave solution
    profiles and their frequencies.  Hyperbolicity and self-similar
    solutions of the modulation equations are used to characterize lattice
    hydrodynamics. The quasi-continuum theory was also leveraged
    elsewhere, such as in characterizing the leading edge amplitude of the
    DSW and the periodic/traveling wave solutions of the discrete model.
    A truncation of the dynamics in wavenumber space allows for a
    linearized, large t analysis of small amplitude oscillations that
    accompany the RW and DSW solutions.
    
    Beyond the binary oscillations that have been studied previously, this
    work has identified new lattice hydrodynamic features that, so far,
    appear not to have a continuum, dispersive hydrodynamic parallel.  Of
    particular interest is the rapid, unsteady transition between two
    in-phase counterpropagating periodic traveling waves with the same
    frequency (US in Fig. <ref>), identified with a shock
    solution of the Whitham modulation equations.  This solution is born
    out of the bifurcation of another lattice solution particular
    to the lattice, the
    counterpropagating DSW, RW pair separated by a stationary shock (DSW+SS+RW
    in Fig. <ref>) when the DSW merges with the stationary
    shock.  The US represents an unsteady generalization of the steady
    periodic-to-equilibrium heteroclinic solutions of higher order
    continuum dispersive equations, themselves interpreted as admissible
    shock solutions of the corresponding Whitham modulation equations
    <cit.>.  This work points to a
    new class of admissible “Whitham shocks” in the lattice context.  It
    will be interesting to see if continuum solutions of this type also
    exist.
    
    The asymptotic and semi-analytical tools used in this work constitute
    an effective framework in which to investigate the hydrodynamics of
    other lattice dynamical systems for which dispersive shock wave
    phenomena may arise.  For example, the Riemann problem for FPUT chains
    generalizes to two wave families (second order in time) the problem of
    a single wave family studied here and has been shown to exhibit a
    variety of lattice hydrodynamic solutions including steady transition fronts <cit.>
    and DSWs <cit.>.
    
    A number of interesting open questions remain. Here, we list some of
    these. While the existence of periodic traveling waves has been proven
    <cit.>, their stability and the existence, as well as stability
    of a more general class of traveling waves consisting of
    periodic-to-periodic heteroclinic solutions are important problems to
    better understand lattice hydrodynamics, with implications for finite
    time blow up and the construction of the TDSW solution. In the
    DSW+SS+RW solution, we have not been able to provide a mechanism that
    leads to the particular selection of the intermediate constant
    u_0=(u_+-u_-)/2.  We suspect that this is due to microscopic details
    of the lattice equation that have been neglected in the analysis.  Of
    similar, unknown origin is the selection mechanism
    
    for shock solutions of the binary modulation
    equations (<ref>) investigated in detail in <cit.>.  On the
    other hand, the unsteady shock constitutes a genuine time-periodic
    orbit. Such periodic orbits are worthy of exploration in their own
    right, including from the perspective of dynamical systems by, for
    example, generalizing the spatial dynamics in dispersive
    <cit.> and dissipative <cit.> systems.
    The existence, stability and bifurcation analysis of such states are
    intriguing problems for future exploration.  
    
    
    
    
    
    
    As discussed in <cit.>, the model with quadratic flux
    (<ref>) analyzed herein is one among several
    
    
    possible discretizations of the inviscid Burgers' equation. A
    comparative study of different discretizations of the more general
    conservation law u_t + Φ'(u)_x = 0
    
    could lead to other lattice hydrodynamics and, in the case of
    integrable discretizations, be amenable to deeper mathematical
    analysis. These topics constitute some of the open problems
    
    within the theme of DSWs in one spatial dimension.  The study of
    higher-dimensional models appears to be wide open.
    
    
    Some of these topics are currently under consideration and will be
    reported in future publications.
    
    
    
    
    § ACKNOWLEDGMENTS
    
    The authors would like to thank the Isaac Newton Institute for
    Mathematical Sciences for support and hospitality during the programme
    Dispersive Hydrodynamics when work on this paper was undertaken (EPSRC
    Grant Number EP/R014604/1). This material is also based upon work
    supported by the US National Science Foundation under Grants
    DMS-2107945 (C.C. and E.O.), DMS-2204702 (P.G.K.) and PHY-2110030
    (P.G.K), DMS-2306319 (M.A.H.). This work began during P.S.'s time as INI-Simons Post Doctoral Research Fellow, and this author would like to thank the INI for support and hospitality during this fellowship, which was supported by Simons Foundation (Award ID 316017).
    
    
    
    unsrt
