Shock bifurcation and emergence of diffusive solitons in a nonlinear wave equation with relaxation

A hyperbolic generalization of Burgers' equation, which includes relaxation, is examined using analytical and numerical tools. By means of singular surface theory, the evolution of initial discontinuities (i.e. shocks) is fully classified. In addition, the parameter space is explored and the bifurcation experienced by the shock amplitude is identified. Then, by means of numerical simulations based on a Godunov-type scheme, we confirm the theoretical findings and explore the solution structure of a signaling-type initial-boundary-value problem with discontinuous boundary data. In particular, we show that diffusive solitons (or Taylor shocks) can emerge in the solution, behind the wavefront. We also show that, for certain parameter values, a shock wave becomes an acceleration wave in infinite time, an unexpected result that is the exact opposite of the well-known phenomenon of finite-time acceleration wave blow-up. Finally, the ‘red light turning green’ problem is re-examined.


Introduction
At one end of the spectrum of nonlinear wave phenomena we find the classical solitons [1]. These extremely stable solitary waveforms, whose permanence is due to a balance between nonlinear and dispersive effects, are best known for their remarkable ability to interact in a particle-like manner. At the other end is wave chaos [2], i.e. deterministic wave motion that exhibits an extremely sensitive dependence on the initial data. Somewhere in between lie singular surfaces [3], i.e. propagating wavefronts, across which some physical quantity of interest suffers a jump discontinuity (or jump for short).
Followed closely by acceleration waves [3]- [10], the best known subclass of singular surfaces are, of course, shock waves [3,11]. The classical definition of a shock wave (or shock for short) comes to us from acoustics. There, a shock is defined as a propagating surface across which at least one of the flow variables suffers a jump (see [12], p 218). Of course, shocks arise not only in acoustics, but in many other areas of science: e.g. the evolution of black hole cosmologies [13];Čerenkov radiation [14]; certain flows of Maxwell fluids [15,16], where the wavefronts are known as vortex sheets ([17], p 155); second-sound (i.e. thermal wave) phenomena [5,7,8], [18]- [21]; and traffic flow modeled under the continuum assumption [11,22,23], to mention just a few.
In this paper, we examine the evolution of kinematic 3 shock waves under a nonlinear partial differential equation (PDE) first put forth (but not analyzed) by Lighthill and Whitham in their seminal work on the kinematic-wave theory of traffic flow (see [24], equation (21)). Our aim here is not to present a study of traffic flow per se but, rather, to gain a better understanding of the range of wave phenomena captured by this model. To this end, it is useful to first (briefly) review the basic theory of traffic flow involving diffusion.
Assuming one-dimensional (1D) flow along the x-axis, the simplest continuum-based model of traffic flow that takes into account vehicular diffusion is given by (see, e.g., [22,23]) 3 Here, ρ = ρ(x, t) is the traffic density, which carries units of vehicles per unit distance; q = q(x, t) is the flux, which carries units of vehicles per unit time; the positive constant ν, often referred to as the diffusion coefficient, denotes the diffusivity; and the positive constants v m and ρ s denote the maximum speed in the limit of ρ → 0 and the saturation value of the density, which is such that 0 ρ ρ s , respectively. Furthermore, implicit in (1) is the assumption that the flow of traffic occurs on a long, straight stretch of highway with no exits or on-ramps. Eliminating q between these two equations yields the well-known PDE which, of course, is Burgers' equation [22].
Since it is based on Fick's law (embodied by the second equation in (1)), Burgers' equation, like the well-known diffusion (or heat) equation, is a parabolic PDE. Therefore, it too suffers from the 'paradox of instantaneous propagation' [18], also known as the 'paradox of heat conduction' [25], and thus violates the principle of causality. With regard to traffic flow, this physically unrealistic behavior implies that drivers and their vehicles react instantaneously to changes in the density and its gradient.
To account for both vehicular diffusion and the reaction time of drivers and their vehicles, Jordan [26] turned to the well-known Maxwell-Cattaneo flux law [18]- [21] from hyperbolic heat conduction theory. Specifically, he replaced Fick's law with the constitutive relation where the positive constant τ 0 represents the reaction time of drivers and their vehicles to changes in the traffic flow. Consequently, in place of Burgers' equation (i.e. equation (2)), he obtained (what he termed) the hyperbolic Burgers' equation (HBE): which, for the purposes of the present investigation, we immediately recast as the system of balance laws Here, c 0 = √ ν/τ 0 is the characteristic (or 'sound') speed, and we observe that (5) is a special case of a semilinear hyperbolic system. Except for the recent studies carried out in [26,27], equation (4) appears to have received little attention since it first appeared in [24], over fifty years ago. Thus, the aim of the present paper is to expand upon these earlier works and present an in-depth analysis of shock phenomena described by this nonlinear PDE. In particular, we fully classify the evolution of shocks using singular surface theory, study the stability and bifurcation of the shock amplitude, and provide numerical simulations in support of our analytical findings. We also note how, depending on the parameter values involved, the HBE assumes the character of the wellknown PDEs that comprise its special/limiting cases. Finally, conclusions are stated, the 'red light turning green' problem is briefly re-examined, and we end by pointing out a number of unresolved issues that we find worthy of further investigation.

Shock analysis
Suppose that at time t = 0, the shock amplitude, [[ρ]], is nonzero. Here, [[F]] ≡ F − − F + denotes the jump in a function F(x, t) across the right-propagating (planar) wavefront x = (t), where F ± ≡ lim x→ (t) ± F(x, t) are assumed to exist and the ± superscripts correspond to the regions ahead of and behind , respectively. Then, referring the reader to [26] for the details, it can be shown that ] is the dimensionless shock amplitude, satisfies the following (quadratic) Bernoulli equation: Above, the 1D displacement derivative, δ/δt, denotes the time-rate-of-change measured by an observer traveling with ; µ 0 = (2τ 0 ) −1 − β 0 + 2β 0 ρ + c /ρ s , where the constant ρ + c ∈ [0, ρ s ] denotes the value of ρ immediately ahead of ; and β 0 = v m c 0 /(2ν).
Solving (6) using the standard (i.e. linearizing) substitution S = 1/s, its exact solution is readily found to be [26] where S(0)( = 0) denotes the value of S at t = 0. Additionally, the constant α * , known as the critical initial amplitude [3], is given by where Ma = v m /c 0 is a Mach number, ∈ [0, 1] is a constant, and here and henceforth we let ρ + c = ρ s to simplify our notation. From a strictly mathematical standpoint, (7) indicates that the evolution of S can occur in any one of the following ten ways.  For α * = 0: In cases (iv), (viii) and (x), the breakdown time t ∞ (> 0) is given by where we observe that sgn(α * ) = sgn(µ 0 ) since β 0 > 0. Here, it should be noted that (only) cases (i)-(iv) were discussed in [26], wherein the focus was limited to α * > 0. And we note for future reference that a compression (resp. expansion) shock is such that S(t) > 0 (resp. S(t) < 0) [3].

Shock bifurcation analysis
Although we have an exact expression for S, it is nevertheless instructive to investigate the steady-state behavior of the shock amplitude using qualitative methods. To this end, let us return to (6) and examine the stability characteristics of its two equilibrium solutionsS = 0 andS = −α * , which correspond to the roots of the quadratic equation Now, a standard stability analysis of (10) and its roots reveals that S undergoes a transcritical bifurcation (see, e.g., Strogatz [28]). Here, however, the bifurcation occurs not at a point, as is often the case in textbook examples, but rather along the curve which is obtained by setting α * = 0 in (8) and solving for . What is interesting about this type of bifurcation is the 'swapping' of stability between the two equilibria, which in this case occurs across the (bifurcation) curve given in (11). Observe thatS = 0 andS = −α * are stable and unstable, respectively, when α * > 0, under which fall (only) cases (i)-(iv). On the other hand, when α * < 0 just the opposite is true, i.e.S = 0 and S = −α * are, respectively, unstable and stable, and (only) cases (v)-(viii) can occur. If, however, α * = 0, then the two roots of our quadratic have coalesced at the origin (i.e. zero has become a root of multiplicity two) and only cases (ix) and (x) are possible sinceS is restricted to the bifurcation curve. Equation (12), below, summarizes the stability of the equilibrium solutions in terms of and the Mach number: where we recall that ∈ [0, 1] and Ma > 0.

The MUSCL-Hancock scheme
To numerically solve the system given by (5), we employ the MUSCL-Hancock scheme [29,30], which is a high-resolution Godunov-type scheme. The latter combines a slope-limited extrapolation of the cell-interface values in the spirit of van Leer's MUSCL [29]- [31] with a predictor-corrector time-stepping to achieve a non-oscillatory second-order-accurate (in space and time) scheme. Hence, let us consider a generic system of balance laws: where Q = Q(x, t), F and R are column vectors containing the conserved quantities, their corresponding fluxes and sources, respectively. Note that the system of conservation laws resulting from taking R ≡ 0 is assumed to be hyperbolic. Given the domain [0, L], we begin the numerical solution by constructing a cell-centered grid with N cells; i.e. let x = L/N be the mesh size, thereupon the cells' centers are located at Then, the average of Q over the ith cell at time t = t n (i.e. at the nth time step), denoted by Q n i and assigned to be the value of Q at that cell's center, is the quantity we seek to compute. Finally, by t n we denote the variable time step size.
We approximate the conservation law in (13) via the unsplit (shock-capturing) conservative discretization [29,31] where F n+1/2 i∓1/2 are the fluxes through the left and right interface of the ith cell, respectively, and R n+1/2 i is an approximation of the cell-average of the source term. Also, an n + 1/2 superscript indicates that a quantity is based on a mid-time-step prediction, as discussed below. Furthermore, we henceforth leave the temporal superscript n on Q understood.
The discretization given in (14) in only first-order accurate in space and in time. To obtain second-order accuracy in space we perform slope-limited linear extrapolation of the cellinterface values of Q based upon the cell-center values. In other words, we let be the values of Q at the ith cell's left and right interfaces, respectively. The slope Q i in each cell is given by where θ ∈ [1, 2] is a parameter, and mm(·) is the (generalized) 'minmod' limiter: By limiting the slopes, we obtain a nonlinear scheme, so that it can be both second-order accurate and non-oscillatory (in fact, total variation diminishing [30,31]), i.e. the scheme does not introduce spurious oscillations in the solution.
Upon extrapolating the cell-interface values of Q, we perform a predictor half-time step to obtain the evolved cell-interface values: 7 Then, we set up a local Riemann problem (RP), i.e. a Cauchy problem with discontinuous initial data, on each cell-interface of the grid. And we continue by computing the fluxes through the cell-interfaces from the solutions of the local RPs. However, in order to decrease the computational time and the complexity of the algorithm, we obtain an approximate solution to the RP. Due to the simplicity of the advective part of the HBE (recall (5)), the well-known Lax-Friedrichs flux [30,31], i.e.
is as good as any. Finally, for the cell-average of the source term, we take Note that using the evolved cell-interface values in (19) and (20), and then computing Q n+1 i via (14), constitutes a corrector time step, resulting in a scheme that is (formally) second-order accurate in time.
The stability of any explicit finite-difference method, such as the one described above, is contingent upon the satisfaction of the Courant-Friedrichs-Lewy (CFL) condition, i.e. the time step must be such that for this problem, where C CFL is the CFL number (0 < C CFL 1). In practice, C CFL is determined empirically.
In 1D, there are two types of physical boundary conditions: transmissive and reflective. We implement these boundary conditions via two ghost cells at each end of the grid (see, e.g. [30,31] for more details). In our simulations below, we employed transmissive boundary conditions at both ends of the computational domain.

Numerical results and discussion
Since the study carried out in [26] relied exclusively upon analytical methods, one of the primary aims of the present investigation is to employ advanced numerical techniques to generate accurate simulations that illustrate not only the evolution of the shock amplitude S(t), but, more importantly, that of the HBE's solution profile in the region behind the shock front , i.e. within the interval 0 < x < c 0 t. (See also [31], example 17.7, wherein a physically questionable special case of (5) is numerically solved.) To this end, values for the physical constants must be specified. Here, we set ν = 0.1655 m 2 s −1 and τ 0 = 0.0629 s. These values correspond to the kinematic viscosity and relaxation time, respectively, of the fluid listed as 'Vistanex' 4 in Joseph et al ( [15], table 4). According to these authors, this liquid behaves as a non-Newtonian fluid of the Maxwell type, i.e. a particular type of viscoelastic fluid. We have chosen these values for ν and τ 0 , which yield c 0 ≈ 1.622 m s −1 , along with the initial and boundary conditions given below in (22) so that the reader might compare our simulations with those carried out in [16]. Therein, the equation of motion is the damped wave equation (DWE), which is the linearized version of the HBE that results from letting v m → 0, and the focus is on Stokes' first problem for a Maxwell fluid. Here, we point out that the quantities ρ, x, ν, τ 0 , c 0 , ρ − and ρ + in the present work correspond to u, y, ν, λ, c, V and zero, respectively, in [16]. Additionally, it should be noted that in this particular (Maxwell) fluid flow analogy, (4) corresponds to the equation of motion obtained by adding the (hypothetical) body force-related term v m λ −1 (1 − 2u/V m )u y , where the constant V m denotes the maximum value of |V |( = 0), to the right-hand side of ( [16], equation (2.1)), i.e. the DWE, and then multiplying through by λ.
Having obtained realistic numerical values for the primary physical parameters, we now proceed to numerically solve (5) subject to the following initial and boundary conditions: Here, H (·) denotes the Heaviside unit step function and ρ ± ∈ [0, ρ s ] are constants. Note that our initial condition implies that ρ + c ≡ ρ + (recall (6)), and also that S(0) = (ρ − − ρ + )/ρ s . In addition, for simplicity, we take v m ∝ c 0 in the simulations below. The results are shown in figures 1-4, which are color-coded using the following convection: red corresponds to the solution at t = 0.125 s, blue to t = 0.25 s, green to t = 0.5 s, orange to t = 1.0 s, violet to t = 2.0 s, gray to t = 4.0 s, black to t = 6.0 s and pink to t = 8.0 s. (For those readers viewing this paper in black and white, see the appendix.) Furthermore, though L and N vary depending on the simulation, the mesh size is fixed at x = 5 × 10 −4 m. Additionally, we picked θ = 2.0 as the generalized minmod limiter's parameter, and we used a CFL number of C CFL = 0.9. Finally, all simulations were carried out for ρ s = 1.0 and all limits stated in the figure captions are taken as t → ∞. Figure 1, which corresponds to case (ix), shows the solution of a variation on the 'red light turning green' problem ( [11], example 11.2). Here, we see how the inclusion of vehicular diffusion and the reaction time of the drivers affects the solution of this well-known problem. Specifically, it is clear that the HBE propagates the initial jump from the boundary into the half-space x > 0, with S(t) decaying to zero as t → ∞. In contrast, classical traffic flow theory, which is based on the Lighthill-Whitham-Richards (LWR) model [27], i.e. (2) in the limit of ν → 0, dictates that a rarefaction wave, which is continuous and connects the states ρ = ρ + and ρ = ρ − , emerges from the boundary at x = 0 ( [11], figure 11.21). Though seemingly different, the solutions of the HBE and LWR model for this problem are, in fact, very similar, in the long-time limit. To understand this connection, we observe that, based on the nearly piecewiselinear shape of the profiles shown in figure 1, the following ad hoc approximation can be easily constructed: Thus, figure 1 suggests that the discontinuous solution of the HBE is evolving into the continuous rarefaction wave of the classical theory; in other words, under case (ix), the HBE assumes the character of the LWR model as t → ∞. It should also be noted that while ρ evolves toward a continuous function of x under case (ix), ρ x does not; in fact, it is clear from (23) that [[ρ x ]] ≈ (ρ + − ρ − )/(c 0 t) for t sufficiently large. This means that the decaying shock wave shown in figure 1 also gives rise to an acceleration wave as t → ∞. Such behavior is, of course, the exact opposite of the well-known phenomenon of acceleration wave blow-up in finite time (see, e.g., [3,4]). Figure 2 shows both the decay of the shock amplitude to zero (left panel) and its growth to the value |α * | (right panel), which corresponds to cases (ii) and (vi), respectively. In both simulations, we observe what can only be described as the emergence of a diffusive soliton [1], or Taylor shock [11], behind . Note that much like the classical solitons, which result from the balance of nonlinearity by dispersion, diffusive solitons emerge from the balance of nonlinearity by diffusion [1]. Indeed, it is shown in [26] that (4) possesses a Taylor shock traveling-wave solution (TWS) of the form where 0 ≤ ρ 1 < ρ 2 ≤ ρ s , ξ = x − vt and v = v m [1 − (ρ 1 + ρ 2 )/ρ s ] is the wave speed of the TWS, l = l b [1 − (v/c 0 ) 2 ] is its 'shock thickness' [11,22] with l b = 4νρ s [v m (ρ 2 − ρ 1 )] −1 being the shock thickness of the TWS of Burgers' equation (i.e. equation (2)). Note that for case (ii) (left panel) we have ρ 1 = ρ − and ρ 2 = ρ + , while for case (vi) (right panel) we have ρ 1 = ρ − and ρ 2 = α * .
To support the conjecture that parts of these two solution profiles are evolving into Taylor shocks, we have also plotted ρ −1 s f (ξ − vφ), where φ is a constant temporal shift, in   figure 2 using the '×' symbols. As expected, the agreement between the (shifted) Taylor shock ρ −1 s f (ξ − vφ), where we note that φ is positive for both cases shown in figure 2, and the similar structure that emerges in the numerical solution was found to improve as the latter ages and approaches a steady state. In fact, the Taylor shock that eventually emerges under case (ii) is, qualitatively, very similar to the TWS of the classic Burgers' equation; see [23]. More importantly, however, the emergence of a diffusive soliton under case (vi) (right panel of figure 2) clearly illustrates the two 'faces' of this damped, nonlinear, hyperbolic PDE; specifically, the HBE's remarkable ability to support both viscous (i.e. continuous) and inviscid (i.e. discontinuous) shocks at the same time. Figure 3 depicts the decay of a compression shock (left panel) and the propagation of a constant amplitude shock (right panel) which, respectively, correspond to cases (i) and (vii). In the former, we see that both S(t) and ρ − x go to zero as t → ∞, and so an acceleration wave cannot form. Since this behavior is very similar to that of the DWE (see [16], figure 1), it appears that the effects of nonlinearity are negligible under case (i). In contrast, the right panel of figure 3 shows the HBE behaving exactly like the linear wave equation ρ tt = c 2 0 ρ x x , the solution of which, subject to both (22) and case (vii), is ρ(x, t) = ρ + + ρ s |α * |H (t − x/c 0 ). Evidently, the advection terms effectively cancel out the damping term, ρ t , under case (vii). Figure 4, our last simulation, presents the decay of an initial jump to the stable equilibrium value |α * | under case (v). Clearly, a Taylor shock does not form behind in this case, contrary to case (vi), where the initial jump grows to the value |α * |. Figure 4 also highlights the fact that, for these parameter values, |α * | is a nonzero stable equilibrium.

Conclusions
Using both analytical and numerical methods, we have carried out a study of the HBE, equation (4), which is a hyperbolic generalization of Burgers' equation, in the context of shock phenomena. The singular surface results presented in [26] on the growth/decay of shocks have been extended to include the cases for which α * 0. A bifurcation analysis of the shock amplitude solution was also performed and the question of shock stability was addressed. Additionally, a number of high-resolution simulations of shock propagation, using the MUSCL-Hancock scheme to numerically solve the HBE subject to (22), were presented and discussed. Based on an analysis of these findings, we report the following: • For α * = 0 (resp. α * = 0), S(t) is an exponential (resp. algebraic) function of time (see also section 6.2).
• The shock amplitude S undergoes a transcritical bifurcation along the curve given in (11), which lies in the plane α * = 0 of (Ma, , α * )-space (see section 3). • Under case (ix), both rarefaction and acceleration waves form as the shock amplitude decays to zero; i.e. the HBE assumes the character of the LWR model, a first order PDE, as t → ∞ (see figure 1).
• Under cases (ii) and (vi), diffusive solitons, which evolve into the TWS given in (24), emerge behind (see figure 2). Also, the TWS that eventually emerges under case (ii) is, qualitatively, very similar to the TWS of (2), the classic Burgers' equation (see the left panel of figure 2).
• Under case (i), the HBE exhibits a behavior very similar to that of the DWE, a linear PDE (compare the left panel of figure 3 with [16], figure 1).
• Under cases (iii) and (vii), the HBE behaves exactly like the linear wave equation; however, only case (vii) corresponds to a stable, nonzero equilibrium solution of (6), namely, S = |α * | (as depicted by the right panel of figure 3 and in figure 4).

Final remarks on the 'red light turning green' problem
It is interesting to note that, while the behaviors seen in figure 1, the right panel of figure 2, and the left panel of figure 3 are obviously very different, all three simulations correspond to some version of the 'red light turning green' problem. Mathematically, these differences follow from the fact that cases (ix), (vi) and (i) 'lie' on, below and above, respectively, the bifurcation curve (i.e. (11)). Physically, however, these differences are due to the Mach number, i.e. the ratio v m /c 0 . Observing that the values of Ma taken in these three simulations correspond to v m = c 0 , v m > c 0 and v m < c 0 , respectively, and recalling that c 0 = √ ν/τ 0 , these cases admit the following physical interpretations: (I) v m = c 0 , meaning that v m is equal to the speed at which propagates, corresponds to the situation in which drivers have exactly the critical reaction time τ 0 = νv −2 m [26], i.e. drivers and the highway are perfectly 'balanced' on the bifurcation curve, but S → 0 only algebraically; (II) v m > c 0 implies that drivers are not capable of taking full advantage of the highway, i.e. they react too slowly, and the shock maintains its initial amplitude forever; and (III) v m < c 0 represents a situation in which drivers react more quickly than the highway demands, and thus S → 0 exponentially in time.
More importantly, these three cases can also be considered in terms of highway safety. In particular, the presence of a shock on the highway indicates that the driving conditions are less than ideal (likely dangerous) because traffic cannot 'spread out' evenly (e.g. there is a traffic jam or a slow-reacting driver ahead). Thus, in this respect, (III) is the 'safest' and (II) is the least safe of these three situations. In other words, when drivers have a slow reaction time (i.e. (II), where τ 0 > νv −2 m ) there is a greater danger on the road since the resulting shock never decays! (As a point of reference, Green [32] notes that reaction times typically fall in the range 0.7-1.5 s.)

Future work
While the present study has addressed a number of issues not covered in [26], it has also given rise to a number of new questions that we feel deserve further study. For example: (A) Can an analytical expression for the temporal shift φ be determined? (B) What is the behavior of the shock amplitude in the case where ρ + is a function of time rather than a constant? (C) How does a density-dependent diffusivity, i.e. where ν = ν(ρ), affect the growth/decay of shock waves? (D) What other physical phenomena might be described by the HBE?
In figures 1-4, the solution profiles are shown at different times using different colors. Here, for the benefit of those readers who have access to only a black and white version of this paper, we list the times at which the solution profiles are shown in each frame of each figure.