"Irregularization"of Systems of Conservation Laws

We explore new ways of regulating defect behavior in systems of conservation laws. Contrary to usual regularization schemes (such as a vanishing viscosity limit), which attempt to control defects by making them smoother, our schemes result in defects which are \textit{more singular}, and we thus refer to such schemes as"irregularizations". In particular, we seek to produce \textit{delta shock} defects which satisfy a condition of \textit{stationarity}. We are motivated to pursue such exotic defects by a physical example arising from dislocation dynamics in materials physics, which we describe.

A remarkable feature of systems of conservation laws is the tendency of smooth initial data to give way to non-smooth defects after only finite time. Examples include shocks in hydrodynamics, cracks and dislocations in crystalline solids, and traffic jams in continuum models of traffic flow. Such defects are important in determining the qualitative behavior of a system, but they pose various theoretical and practical challenges, stemming from the fact that the time evolution of defects cannot be determined solely from the evolution of the ambient continuum.
The physical origin of this difficulty can be understood as follows: Continuum equations of motion are derived from more fundamental, microscopic laws by discarding the "high-order" terms which are irrelevant at macroscopic scales. Thus, systems with different microphysics can yield the same continuum equations, and conversely the continuum equations alone do not specify the microphysics. Defects are inherently microscopic, and thus their behavior cannot be inferred from the continuum laws.
The ambiguity of defect behavior in systems of conservation laws is reflected in the mathematical theory by non-uniqueness of solutions containing defects. To obtain uniqueness, one must impose additional conditions upon the solution, presumably corresponding to the microphysics not accounted for by the continuum equations. Common conditions of this sort include entropy conditions and any of various types of regularization. In the former case, inequalities are imposed in analogy with the thermodynamical principle that entropy must always increase, while in the latter case the conservation laws themselves are modified by the addition of small terms that serve to smooth out defects. We will return in more detail to both of these cases in the next section.
For the most common kinds of conservation laws (notably, those of hydrodynamics), the "proper" choice for solutions with defects turns out to be relatively tame, with defects consisting only of jump discontinuities and being described by regularization with a vanishingly small viscosity term. Some more exotic systems feature different classes of defects which require other types of regularization; for example, superfluids are known 10 to exhibit dispersive shock waves which are described by a very small dispersive term. Still another type of defect (which is a primary subject of this paper) is a delta shock, which consists of a delta function-like spike, possibly riding atop a regular (viscous) shock. This latter type of defect is much more poorly understood than the former two and differs in some important regards: Firstly, there is no particular regularization associated with delta shocks (for some systems, delta shocks may even be contained in the vanishing viscosity solutions 5 ). Secondly, there are few physically-realizable systems which are currently known to exhibit delta shocks. Various types of shocks evolving from smooth, sinusoidal initial conditions, depending on different regularizations of the Hopf equation: viscous shock (blue), dispersive shock (green), and delta shock (red).
The present paper was motivated by the discovery of arXiv:1506.05743v1 [math.AP] 18 Jun 2015 delta shocks in a system arising from materials physicscontinuum dislocation dynamics (CDD)-which to our knowledge provides the best physical interpretation of delta shocks currently available. These delta shocks have presented substantial challenges to physicists interested in simulating CDD numerically, which has led to a desire for new numerical methods tailored to better handle such exotic defects. We explore in this paper (in a mostly informal way) perspectives on delta shocks which may be of use to practitioners interested in simulating them numerically. (We do not, however, explicitly construct any numerical schemes here.) In particular, we describe a condition of stationarity which forces the formation of delta shocks, and we construct a new kind of "irregularization" which also produces delta-like defects (though we do not rigorously prove these to be true delta shocks).
The outline of this paper is as follows: In the first section we provide some background on classical shocks. In the second section we describe delta shocks and discuss various ways in which they can arise in systems of conservation laws, and we introduce the condition of stationarity and the irregularization mentioned above. Then in the third section we describe in more detail the system which initiated the preceding work-continuum dislocation dynamics.

I. BACKGROUND
In this section, we review basic facts about systems of conservation laws and standardize our notation. A system of conservation laws (or simply "conservation law") is a partial differential equation (PDE) of the form where u : Ω × R + → R m (for Ω ⊂ R n ), u = u(x, t), is the unknown, and F : R m → R n is a given flux function.
The conserved quantity u we refer to as "mass" in the following discussion. In the case where u is scalar valued (m = 1), we refer to eq. (1) as a scalar conservation law. Our working example of a conservation law will be the Hopf equation (a.k.a. the "inviscid Burgers' equation"), which is a scalar conservation law in one spatial dimension (n = 1). The Hopf equation is one of the simplest examples of a PDE which develops discontinuities in finite time 4 . Such discontinuities and other non-smooth features appearing in solutions to PDE we refer to collectively as defects. In the presence of discontinuities in u, the derivatives appearing in eqs. (1)(2) are no longer well defined, and so we must generalize our notion of "solution" to a PDE in order to make sense of the conservation law at later times. Such a generalization is referred to as a weak solution. Several distinct notions of weak solution exist, such as integral solutions, viscosity solutions, variational solutions, and various types of regularized solutions, to name a few 4 . Not all forms of weak solution are applicable to all PDE's. The types of weak solution which are most relevant for systems of conservation laws are regularized solutions and integral solutions. In the former case, the PDE is modified by the addition of a small term that serves to smooth out the defects; for the Hopf equation, such a modification looks like where η ∈ R controls the strength of the regularization.
To recover a solution to the original system of conservation laws, we take a limit as η → 0. Note that, in principle, the form of the regularizing term should correspond to the microphysics of a given system-so for example, systems with viscosity have a second-order regularization, η∂ 2 x u, whereas inviscid systems (like plasmas or superfluids) typically have dispersive (third-order) regularization, η∂ 3 x u. Below, in eqs. (14)(15)(16), we will consider a different form of regularization altogether, which nonetheless bears a familial resemblance to eq. (3).
The main alternative kind of weak solution, the integral solution, proceeds instead by integrating the original PDE to remove the problematic derivatives. Integral solutions are not of direct interest to us here, but we point out one tangential relationship to our current discussion: Integral solutions are generally not unique, and to select a proper weak solution requires additional microphysical information. This information is usually specified as a relationship between the values of the solution u on either side of a defect; when the relationship takes the form of an inequality it is called an entropy condition 1,2 , and when it is an equality it is called a kinetic relation 3 . The stationarity condition we propose below is similar to a kinetic relation, but with the generalization that the speed of the defect can also figure into the equality.

II. DELTA SHOCKS
For our purposes, a delta shock (also written "δshock") may be defined informally as any defect which has a finite amount of mass concentrated at a point. We note in passing that there are also meaningful notions of δ -shocks (and δ (n) -shocks, more generally) corresponding to defects containing derivatives of delta functions 9 . Our numerical solutions (Section II C, Figs 2 and 3) appear to contain an upward and downward-pointing singularity -the sum being the strength of the delta shock and the difference a δ component. In this section we shall derive conditions for the evolution of the strength of the delta shock component.
Delta shocks are known to arise naturally in a number of systems of conservation laws-for example, the system considered by Tan, Zheng, & Zhang 5 , and the system In the former case, the delta shocks have been shown to be vanishing viscosity limits of the regularized system The mechanisms which cause delta shocks to form in systems of conservation laws are not fully understood 8 , though in some cases the necessity of delta shocks is obvious. We highlight two such cases here: Firstly, if a conserved quantity v is a derivative of another variable u, and the latter exhibits regular (viscous) shocks, then clearly v must contain a delta shock. For example, if we take a derivative of the Hopf equation (2) and set v : which is precisely the latter half of system (4). Thus (provided the initial data for system (4) satisfies v(x, t = 0) = u (x, t = 0)) we see that this system can be interpreted as the evolution of the conserved quantity of the Hopf equation, along with that of its derivative.
This interpretation is not the only available, nor is it necessarily the best. Indeed, the closely related system 5 where F is an arbitrary smooth & convex function and g is an arbitrary smooth & increasing function, also displays delta shocks, though they cannot in this case be interpreted as derivatives of regular shocks. A more general interpretation for the appearance of delta shocks in such conservation laws is as a consequence of incompatibility of defect motion and the flow of the conserved quantity, which we explore in the next two subsections.

A. Generalized Rankine-Hugoniot relation
The Rankine-Hugoniot relation is a statement of conservation of mass across a shock. For the case of a regular shock in one spatial dimension which moves with speed σ, mass conservation requires that the net flux into the shock balance the mass of the "hole" left as the shock moves away. In symbols, this reads where u l and u r are the values of u to the immediate left and right (resp.) of the shock (see Evans 4 for a more formal derivation of this relation). If we allow the possibility of delta shocks storing mass on the defect, it is now the combination of the hole behind the shock and the mass m of the delta which must balance the incoming flux. The modified Rankine-Hugoniot relation is thus Note that both eqs. (9) and (10) are vector identities in general, and may be interpreted as specifying the shock speed σ and delta mass m (or the derivative thereof) in terms of the local environment, given by u l and u r . In particular, this implies that eq. (9) is only satisfiable in general when u is scalar, since vector u leads to an overspecification in (9) of the single unknown σ. Eq. (10) has no such problem, since m has the same number of components as u, so that there is always one more unknown than there are equations in (10). This shows that, except in the presence of high degeneracy so that eq. (9) is soluble for σ, we expect multicomponent conservation laws to require delta shocks in order to conserve mass in the vicinity of a defect.
This provides our preferred explanation of the delta shocks in system (4): the classical Rankine-Hugoniot relation (9) for this system reads: which does not hold unless u l v r = u r v l (i.e. the necessary degeneracy for (9) to be satisfiable does not generally exist). A similar statement shows why the more general system (8) also displays delta shocks.

B. Stationarity
The Rankine-Hugoniot relations illustrate how multicomponent systems can develop delta shocks when a defect's motion cannot be simultaneously compatible with the flow of all components of the conserved quantity. We can also produce a similar effect even in scalar conservation laws by insisting that a defect move according to specified kinematics which are not compatible with the flow of the scalar conserved quantity. The simplest such example is to demand that a defect be stationary. This might be an appropriate kinematic constraint in, for example, an electrical network containing fuses, where a blown fuse is modeled as a stationary defect. A stationarity condition could also be physically appropriate for the continuum dislocation dynamics example we will describe below. This stationarity condition is analogous to a kinetic relation 3 , with the generalization that the speed of the defect can now figure in as well (rather than just the value of the local conserved quantity).
As an illustration, we consider the Riemann problem for the Hopf equation (2). Mass flows into the defect at x = 0 with rate F (0 − ) − F (0 + ) = 1 2 · 1 2 − 1 2 · 0 2 = 1 2 , and thus by the modified Rankine-Hugoniot relation (10) with σ = 0 we find that the mass of the delta at x = 0 grows at a rateṁ = 1 2 . Everywhere else the value of u is uniquely specified by projection of characteristics. Hence a reasonable candidate solution for all times t ≥ 0 is We note that this candidate solution is not an integral solution in the usual sense, as the square of a delta function (needed to make sense of an integral solution to eq. (2)) is undefined by classical distribution theory. (However, it is possible to overcome this technicality using generalized distributions 6 .)

C. Regularizations for stationarity
Our ad hoc construction of a solution (13) to the Hopf equation which satisfies the condition of stationarity does not lend itself well to numerical simulation. We thus find it advantageous to try to construct a "irregularization" of the Hopf equation which yields such a solution. Our attempted irregularization looks like: where the smooth & bounded function f ( , u x ) satisfies: For example, is such an irregularization. The intuition behind this is that when a shock develops in the Hopf equation, u x diverges to −∞, sending the modified flux 1 2 u 2 f ( , u x ) to zero. The quenching of the flux will hopefully cause a pileup of mass behind the defect. On the other hand, for small values of and away from defects, condition 3 above implies that f ≈ 1, and the original Hopf equation is recovered (approximately).
(This prescription for the irregularization is not as general as possible-for example, one also could add to the denominator of eq. (15) a non-vanishing polynomial in u, which would reflect some sort of self-interaction effect. However, the "minimal" irregularization of eq. (15) is sufficient to produce the defects we are interested in. Moreover, the effect is not peculiar to this particular form: We have also experimented with the irregularization and seen the same qualitative behavior. ) For numerical work, it is sometimes helpful to also add a small regularizing term to ensure solutions stay relatively smooth-e.g.
where k is typically 2 or 4.

FIG. 2: Irregularized Hopf equation. (Color online)
Simulation of the irregularized Hopf equation (16) starting with sinusoidal initial conditions, for two values of at time t = 1 (simulated with a conservative upwind differencing). Note that the singularity has both a positive and negative peak. Only the sum contributes to the delta shock (as a weak limit); the difference can be described as a derivative of a δ-function.

Does this succeed in making the defect stationary?
We begin to answer this by examining numerical results. Fig. (2) shows simulations of eq. (16) for = 0 & = 10 −5 on the unit interval with periodic boundary conditions, for parameters t = 1, k = 2, η = 10 −6 , and a grid spacing dx of .002.
Clearly the irregularization has substantially slowed the movement of the shock and thereby induced a deltalike spike to conserve mass. The extent to which this is a bona fide delta shock is a delicate question (we will thus cautiously refer to this feature as a delta spike). We mention three relevant considerations: 1. To begin with, the effect is inherently non-smooth, as for smooth solutions of eq. (15), we can prove the following maximum principle: Proposition 1. A solution u to eq. (15) on a compact domain satisfies max{u(·, t)} = max{u(·, 0)} (17) for all times t > 0 while the solution remains smooth.
We motivate this as follows: If we look at a curve y(t) such that u(·, t) attains a maximum at y(t), then we find Where the last expression vanishes because u x = 0 at the maximum y(t). This shows that the maximum value of u is unchanging in time. (This argument can be made more rigorous, avoiding the assumption that the maximum can be parameterized by a differentiable curve y(t), but we will not do so here.) Note that a similar result holds also in the case of eq. (16) with a 2nd order (viscous) regularization, with the modification that the maximum of u is then decreasing in time.
2. Examining the delta spike closely reveals an internal structure. To provide some perspective, we note first that eq. (15) admits two relevant exact solutions, one of the form a(t) + b(t)x (i.e. a straight line, moving in time), and the other a stationary solution of the form Experimentally, we find that the shape of the spike is well approximated by a profile of the form (19), as shown in fig. (3). More precisely, the form of the delta spike consists of two partial coshes, one positive and one negative, separated by an abrupt discontinuity. Away from the defect, the solution to eq. (15) approaches the aforementioned linear solution asymptotically in time. This eventual shape is analogous to the asymptotic "N-wave" profile of solutions to the regular Hopf equation. As discussed in the appendix, it is possible to derive analytic expressions for the parameters a(t), b(t), α, β describing this asymptotic form, and we find good agreement between simulation and these analytic values.

FIG. 3:
Cosh profile of delta spike. Comparison of a delta spike with a cosh profile of form (19). Note that the cosh profile shown here is not a best fit, but rather an analytically derived approximation, as discussed in the appendix. The delta spike comes from a simulation with parameters = 2.5e − 6, dx = 2.5e − 4, and initial data u(x, 0) = sin(2πx) + .1. Again, the positive and negative spikes must be added together to give the net weight of the delta shock.
An important feature of the cosh profile (19) is that its width scales like √ . To conserve mass, the height of the shock must therefore scale like 1/ √ . This suggests that, although the delta spike is not truly a delta shock for non-zero values of (indeed, the cosh profile (19) certainly does not qualify as a delta shock), as → 0 the spike does indeed converge to a delta function. Experimentally, this is what we seem to find (conditional upon some finessing of the numerical methods), as discussed in the next point.
3. Establishing numerical convergence of our simulations is tricky for two reasons: Firstly, since we are looking for convergence to a delta shock, we cannot hope to have convergence in e.g. the L 2 sense. Secondly, the equation (15) is highly unstable (not surprisingly, since the irregularization is designed specifically to produce a delta shock, which may be viewed as a sort of instability). After trying various finite-differencing schemes, we found that the best results came from a type of conservative upwind differencing, wherein u is differenced upwind but u x is not. Explicitly, the semi-discrete formulation of eq. (14) is where ∆ is the discrete forward difference operator (∆a i = a i+1 − a i ). All figures presented herein were produced using this discretization. (Remark: for the regularized eq. (16) we performed an operator splitting, evaluating the regularization term in Fourier space for stability.) Despite being the most successful numerical method we found, the stability properties of this discretization are not entirely satisfactory. In particular, we found that when the grid size became substantially smaller than the width of the delta spike (which is on the order of √ ), the delta spike would become unstable and break apart into multiple, smaller delta spikes. To overcome this difficulty, we let the mesh size dx and the width of the delta √ go to zero together, keeping dx/ √ fixed (recall that it was our intention to let → 0 anyhow, so as to recover a solution to the original Hopf equation). This technique appears to work, as seen in fig. (4): The delta spike gets narrower and higher, while remaining in essentially the same place, as dx, → 0 in this manner. Although our simulations seem to establish convergence to a delta shock, it is a very peculiar convergence wherein the strength of the irregularization must vanish together with the discretization length dx. The fact that the width of the delta spike must be kept close to dx suggests that there may be some sort of resonance between the irregularization of eq. (15) and the grid which drives the formation of the spike (especially in light of the maximum principle mentioned previously). Nonetheless, we can still hope to produce a meaningful weak solution for the Hopf equation in this way. We will not establish anything more rigorous at this time, though.

III. PHYSICAL MOTIVATION: DELTA SHOCKS AND DISLOCATION WALLS
Our interest in delta shocks is motivated by a tangible physical question. Why do dislocations in crystals form walls 20,21 ? Briefly, a crystal is a regular array of atoms. A dislocation line is a flaw in that array, such as the edge dislocation formed by the boundary of an extra plane of atoms, or the screw dislocation forming the central line where planes form a 'spiral staircase' 13 . These dislocations move to mediate plastic deformation when the crystal is bent, and form tangles that organize into walllike structures called cell walls 16 . In a continuum theory of dislocation dynamics, such walls must be described as delta shocks: their density scales as the inverse square of the lattice constant, which vanishes in the continuum limit.
One class of continuum dislocation dynamics theories 15,20 do form such delta shocks. In these theories, the dislocations move with a common velocity F. If one describes an incipient wall in the yz plane with a dislocation density that depends only on x, the dependence of this velocity on the dislocation density simplifies 14 to the Hopf equation: When the Hopf equation forms a step-like shock, the entrained dislocations pile up into a delta shock, forming a wall.
The central physical question is how that wall should evolve after it forms. If we regularize the Hopf equation (21) into Burger's equation (i.e. viscous regularization), the dislocation wall moves along with the Rankine-Hugoniot velocity. This produces a self-consistent, sensible model 18,21 , but one that is unsatisfying in two regards. First, the velocity of dislocation walls in crystals is determined by the microstructure of the walls and not by continuum properties. Secondly, generalizations of this continuum theory have differing velocities for the different components of the dislocation density -smearing the resulting walls.
On the other hand, physically at some junctions where two dislocations intersect, the point of intersection can be pinned in place (so-called sessile dislocation junctions 13 ). Dislocation walls can also be pinned by impurities that segregate to the boundary. This provides the motivation for our study in section (II B) of the Hopf equation (21) with stationary shocks.
We should note that the behavior of both our model in two and three dimensions and experimental dislocation systems under stress is more complex than simple wall singularities. There they produce complex cellular structures, which in our models and some experiments form self-similar, fractal morphologies. Our theories in higher dimensions show clear analogies 22 to the behavior of the Euler equation (the inviscid limit of Navier-Stokes). Nonetheless, these structures are complex, ramified wall-like entities, whose dynamics should be controlled by physics on the microscopic atomic scale, not by the continuum laws.

IV. CONCLUSIONS
We have seen that delta shocks arise naturally in a variety of systems of conservation laws, and their presence can often be traced to an incompatibility of the flow of the conserved quantity and the motion of defects. Although most of the existing examples of delta shocks do not have direct physical interpretation, the case of cell walls in crystals does appear to furnish such an example.
A recurring theme in our presentation has been the importance of developing numerical techniques capable of handling defects with more exotic behavior than traditional, viscous shocks. Ideally, one would like to be able to specify whatever defect behavior is believed appropriate for a given system and have numerical techniques which respect that behavior. In lieu of such a very general approach to simulating conservation laws, we must instead focus on specific classes of defect behavior which we believe are important. To this end, we propose further study of the stationarity condition introduced above.
We have concocted an "irregularization" (14-16) of the Hopf equation which appears numerically to exhibit delta shocks, suggesting a possible avenue to development of the aforementioned numerical methods for delta shocks. Though we are far from showing that these apparent deltas are bona fide, we nonetheless think that the irregularization and the defects it yields are interesting in their own right and merit continued study.

V. APPENDIX
We present here some further numerical and analytical observations (mostly without details or proofs) for the irregularization In particular, we sketch a systematic way to approximate the asymptotic form of solutions to this equation. Three useful characteristics of the regular Hopf equation can be shown to also hold for smooth solutions of this more general equation: 1. The maximum principle mentioned above.
2. Zero-crossings of a solution to eq. (22) remain fixed in time.
3. The total mass between any two zero-crossings is conserved.
The second and third items above allows us to predict (analytically, approximately) the asymptotic form of solutions to eq. (22). We sketch this analysis as follows: Shocks tend to form on downwards-going zero crossings, and away from shocks the solution becomes linear at large times. The overall profile of the solution thus is linear everywhere except at the shocks, where it is approximated by the cosh profile mentioned above. The stationarity of upwards-going zero crossings allows for analytical evaluation of the slope there, which asymptotically must equal the slope of the entire linear region of the solution. Explicitly, the slope u x at a zero crossing x zc satisfies the ODE 0 = ∂ t u x (x zc , t) + 1 2 u(x zc , t) 2 1 + u x (x zc , t) 2 xx = ∂ t u x (x zc , t) + u x (x zc , t) 2 1 + u x (x zc , t) 2 (23) which is obtained by differentiating eq. (22) w.r.t. x and noting that all terms with a factor of u vanish. This ODE (23) can be solved analytically, providing the asymptotic form of the linear part of the solution. The third item above allows us to evaluate the mass of the delta, which is just the difference between the total initial mass between two zero-crossings and the eventual mass of the linear solution on the same region. Examining the cosh profile of the delta mentioned above, we see that α must equal the height at which the delta begins, which can be determined from the linear solution. Thus α is determined, and then β can be determined from the known mass of the delta. Stringing together the above arguments allows a piecewise description of the solution across the entire domain. This analytic, approximate solution is found to agree reasonably well with simulation-for instance, the cosh profile shown in fig. (3) uses parameters calculated in this way.

VI. ACKNOWLEDGEMENT
We thank Randall LeVeque for helpful comments and guidance on hyperbolic conservation laws and numerical methods. We are especially indebted to Alexander Vladimirsky for his exceptionally insightful suggestions for the production of this paper. This work was supported by US DOE/BES grant DE-FG02-07ER46393.