1. Introduction
The modeling of vehicular traffic flow via mathematical equations is a key building block in traffic simulation, state estimation, and control. Important ways to describe traffic flow dynamics are microscopic/vehiclebased
[53, 46, 5], cellular [45, 11], and continuum models. This last class is the focus of this paper, particularly: inviscid macroscopic models [43, 54, 61, 50, 51, 37, 4] that describe the spatiotemporal evolution of the vehicle density (and other field quantities) via hyperbolic conservation laws. Other types of continuum models exist as well, including gaskinetic [26, 52, 28], dispersive [35, 34], and viscous [32, 33] models. Hyperbolic models do not resolve zones of strong braking, but rather approximate them by traveling discontinuities (shocks) whose dynamics are described by appropriate jump conditions [15]. Macroscopic models play a central role in traffic flow theory and practice because:
Computationally, a macroscopic description is a natural framework to upscale millions of vehicles to a celltransmission model [11]
with much fewer degrees of freedom.

Societally, traffic descriptions that do not resolve individual vehicles are desirable for privacy and data security.
In this work, we focus on the laneaggregated description of traffic flow dynamics on uniform highways without any road variations, let alone intersections or bottlenecks. The reason is that even in this simple scenario, real traffic flow tends to develop complex nonlinear dynamics, particularly the phantom traffic jam phenomenon [31, 24]: initially uniform flow develops (under small perturbations) into nonlinear traveling waves, called jamitons [19]. This occurrence of instabilities and waves without discernible reason has been demonstrated and reproduced experimentally [58, 56]. While these features can be reproduced in microscopic carfollowing models, a key goal is to capture these nonequilibrium phenomena via macroscopic models (to facilitate the model advantages described above).
The archetype macroscopic model is the LighthillWhithamRichards (LWR) model [43, 54]
(1) 
that describes the evolution of the vehicle density where is the road position and is time. The fundamental diagram (FD) function , where the equilibrium velocity function is the bulk flow velocity as a function of density, is motivated by the 1935 measurements by Greenshields [22], and many types of FD have been proposed [43, 20, 61, 47, 11]. As a matter of fact, real FD data exhibits a substantial spread in the congested regime [31]. More complex traffic models capture this spread [9, 55, 17, 16], but the LWR model does not. Yet, due to its simplicity it nevertheless is widely used. Moreover, as we highlight below, it also is motivated as a reduced equation for more complex models.
Another critical shortcoming of the LWR model is that it cannot reproduce the phantom traffic jam phenomenon: being a firstorder model, it exhibits a maximum principle, and thus small perturbations to a uniform solution cannot amplify (instead, they turn into Nwaves and decay). In this work, we focus on secondorder models that augment the vehicle density by an independent field variable for the bulk velocity , and describe their evolution via a balance law system, specifically: a hyperbolic conservation law system with a relaxation term in the velocity equation. Due to conservation of vehicles, the density always evolves by the continuity equation, . In turn, the velocity equation encodes the actual modeling of the vehicle dynamics and interactions. The PayneWhitham (PW) model [50, 64]
(2) 
was the first secondorder model proposed. Here is the desired velocity function, and is the relaxation time that determines how fast drivers adjust to their desired velocity . The traffic pressure models preventive driving. Even though the PW model does capture traffic waves accurately [19, 55], it is generally rejected [12] due to spurious shocks that overtake vehicles from behind; and other hyperbolic models are preferred (see below). However, the fundamental structure of a hyperbolic system with a relaxation in the second equation, is common to all models of interest in this study.
Models with the structure described above possess a critical phase transition. If the
subcharacteristic condition (SCC) is satisfied, then uniform flow is stable [63, 64, 44, 7]. Conversely, when it is violated, uniform flow is unstable and nonlinear traveling wave solutions exist [41, 29, 48, 19, 55]. The SCC is defined as follows. Let be the two characteristic speeds of the hyperbolic part of the model, and let be the characteristic speed of the reduced equation (1) (with ), which arises in the formal limit ; in which relaxes infinitely fast to . Then the SCC is: .The case of the SCC satisfied is well studied [63, 64, 44, 7, 42]. In particular, it is related to positive diffusion when conducting a ChapmanEnskog expansion of the model [35, 25]. In contrast, this paper focuses on understanding the behavior and stability of solutions when the SCC is violated.
This paper is organized as follows. In §2, we introduce the equations. Then we characterize the nature of the instabilities to uniform flow, and the traveling wave solutions that then arise: the jamitons. In §3, a systematic computational study of the stability of jamitons is conducted. Those results then motivate a stability analysis of those nonlinear traveling waves, presented in §4. We close with a discussion and a broader outlook in §5.
2. Macroscopic Traffic Models with Instabilities and Traveling Waves
While the general results and methodologies apply to a wide class of secondorder models with relaxation (including the PW model (2) and generic secondorder models [38, 16]), we focus this study on the inhomogeneous AwRascleZhang (ARZ) model [4, 66]. In nonconservative form it reads as
(3) 
where is called the hesitation function. We assume that: is strictly decreasing, is strictly concave, is strictly increasing, and is strictly convex. In particular these assumptions yield a hyperbolic system, which has no waves that overtake vehicles (the 2waves are contacts) [4]. While originally proposed in homogeneous form, the addition of the relaxation term [21] allows for the violation of the SCC.
In the homogeneous ARZ model, the field can be interpreted as a convected quantity moving with the flow (the hesitation function reduces the empty road velocity by ). Hence, the conserved variables are and , and the conservative form of the equations is
(4) 
with associated RankineHugoniot jump conditions
(5) 
Here denotes the jump of the variable across the discontinuity, and is the speed. In addition, the Lax entropy conditions [15] must be satisfied. Specifically: one family of characteristics goes through the discontinuity, while the other converges into it (for a shock), or is parallel to it (for a contact). In particular, the assumptions on made below (3) guarantee that the entropy conditions are equivalent to: the shocks are compressive (i.e., as vehicles go through a shock, the density increases) and move slower than the vehicles [55].
The characteristic speeds of (4) are:
(6) 
where the is genuinely nonlinear (associated with shocks and rarefactions), while the is linearly degenerate (associated with contacts).
2.1. Specific model functions
While the analysis and general results derived below hold for generic models (4), the computational study and the illustrative graphs are presented for a specific choice of model functions. As in [55], we choose , , and construct the fundamental diagram function
that is a smoothed version of the NewellDaganzo triangular flux [47, 11]. The parameters are chosen , , and to have the function fit real sensor data [55]. Hence . Moreover, we choose , and the relaxation time s. Note that these values are for a single lane. When considering multilane traffic, realistic values result by scaling and by the number of lanes.
2.2. Linear stability of uniform flow
Before analyzing the stability of nonlinear waves, we discuss important aspects regarding the stability of uniform flow, i.e., base state solutions of (3) in which and are constant in space and time. The linear stability analysis itself is a wellestablished normal models analysis [32, 19], and we briefly outline the key steps. Consider infinitesimal wave perturbations (where is the wave number and the complex growth rate) of the base state ,
substitute the perturbed solution and into (3), and consider only constant and linear terms. This leads to the system
(7) 
for the perturbation amplitudes, where , , and . Nontrivial solutions can only exist if the matrix in (7) has vanishing determinant, which requires
where satisfies . Writing in terms of its real and imaginary part yields the two equations and , which then leads to the following quadratic equations for :
(8) 
Here and . The positive solution of (8), as a function of , is
(9)  
(10) 
This function has the following properties:

.

, which follows from (9) and the asymptotic () formula:

It is strictly monotonic if , i.e., it is strictly increasing if and strictly decreasing if . This fact follows from (9), because the sign of the term determines the slope of : if , it is constant; and if the term is positive (negative), the function goes up (down) with .
The growth rate of normal modes is
Linear stability, i.e., , is equivalent to (only the negative root of could cause positive growth). Hence, stability holds exactly if , or equivalently , or equivalently
(11) 
This last condition is exactly what the subcharacteristic condition (SCC) [63, 64] yields as well [55]: the LWR characteristic speed, lies in between the two ARZ characteristic speeds, and , exactly if (11) holds.
To recap, for the inhomogeneous ARZ model (3), there are exactly two possibilities: Either the stability condition (the SCC) (11) holds; then all basic wave perturbations have nonpositive growth rates, and solutions are linearly stable. Or (11) is violated; then all waves grow. Moreover, the rate of growth is an increasing function of the wave number , that has , and approaches (as ) the asymptotic growth rate
Figure 1 shows the growth rate functions for the specific model given in §2.1, with stable base states in panel 0(a) and unstable base states in panel 0(b). In the latter, one can clearly see the strict increase of with , and the asymptotic limit . Panel 0(c) shows a plot of the asymptotic growth rate as a function of .
Clearly, base states that satisfy (11) are wellbehaved. However, with regards to modeling phantom traffic jams and jamitons, we are particularly interested in base states that violate (11). These require some more careful discussion. While instabilities to uniform states are ubiquitous in science and engineering, having a growth rate that is increasing for all wave numbers is unusual. The much more common scenario (for example, fluid instabilities moderated by viscosity or surface tension [13]) is that medium wave length are unstable and short waves (i.e., large) are stable again, yielding a critical wave number of maximal growth. In that case, one can argue that out of infinitesimal perturbations, in which all wave lengths are present, the linearized dynamics will single out the ones with dominant growth. Hence, the wave number will be selected to first enter the nonlinear regime.
However, arguments of that type do not work for (3) because, as we have shown, its growth function has no maximum. Rather, the shorter the waves in the perturbation, the faster their growth. It should be stressed that despite this behavior, the linearized model for (3) is mathematically wellposed: for any final time , the amplification of normal modes is bounded by . Still, from an application perspective, properly answering the question of which wave lengths dominate once an amplified perturbation leaves the linear regime, is important; but it is more challenging than in the usual situation.
While the PDE model (3) has no maximum wave number, reality does, namely the vehicle scale. Specifically, wave numbers beyond a , given by the minimum spacing between vehicles, have no practical meaning. One possible way to exclude features on such unphysically short length scales is to add a small amount of viscosity to the ARZ model (3), as in KernerKonhäuser [32, 33] for the PW model (2). In Fig. 0(b), this would change the functions to drop off once gets close to the vehicle scale. Similarly, the numerical discretization of the PDE (3) on grids that are never finer than the vehicle scale will produce a wave number cutoff via numerical viscosity of the method [39].
Another possibility (employed here in §3) is to consider small perturbations, rather than infinitesimal perturbations, and provide a model for the noise. Specifically, we argue that on real roads, perturbations of all wave lengths will act: due to small variations in road features, wind, etc.; and due to variabilities across vehicles. The simplest such noise model is one where all wave numbers appear with equal amplitudes, and perturbations with do not occur.
Because the growth function tends to have a plateau near (see Fig. 1), this linear growth/noise model will yield that all wave numbers near but below will be amplified to reach the nonlinear regime at the same time. This is not unrealistic, as it means that noise close to the vehicle scale will dominate before systematic nonlinear wave effects kick in.
As a final remark we wish to point out that once solutions of the ARZ model (3) leave the linear regime (around a uniform base state), the nonlinear dynamics tend to turn those vehiclescale waves into oscillations with shocks that then collide and merge to form nonlinear wave structures of much smaller amplitude to wavelength ratios. However, those nonlinear transient dynamics are extremely complicated, and this insight is merely based on our observations from numerous highly resolved computations (like those done in §3). What we will study, though, is the stability of true traveling wave solutions of (3) (jamitons) in the situation when the SCC (11) is violated (see §4).
2.3. Traveling wave analysis and jamitons
Before studying waves, it is important to stress that macroscopic models (without explicit lane changing) can equivalently be written in Lagrangian variables. In (4) the equations are cast in Eulerian variables and . The Lagrangian formulation, as used in [21, 55], employs the variables and , where is the (continuous) vehicle number, defined so that , and is the specific traffic volume, i.e., the road length per vehicle. In these variables the ARZ model reads as
(12) 
where and . The assumptions on the model functions in Eulerian variables (, , , ) translate to the following assumptions in Lagrangian variables: , , , and . For simplicity, we now omit the hats, unless explicitly required for clarity. The characteristic speeds of (12) are and , and the associated RankineHugoniot shock jump conditions are
(13) 
where is the propagation speed of the shock in the Lagrangian variables (in the Eulerian frame is the flux of vehicles through the shock). Note that, for contact discontinuities, the conditions are: and .
Below, we are going to employ both types of (equivalent) descriptions of the ARZ model. Eulerian (4) for the computational study of the nonlinear model in §3, and Lagrangian (12) for the jamiton stability analysis in §4.
Jamiton solutions can now be constructed via the Zel’dovichvon NeumannDöring (ZND) theory [18]. One starts out with a traveling wave ansatz. In Eulerian variables, one seeks for solutions , of (4) that depend on the single variable . In Lagrangian variables, one considers solutions , of (12), where . Here is the traveling wave speed in the road frame, while the Lagrangian wave speed relates to the mass flux of vehicles through the wave.
We start with the Lagrangian formulation [55]. The traveling wave ansatz leads to
(14)  
(15) 
Equation (14) yields that
(16) 
where is a constant of integration. Using (16) to substitute by in (15), we obtain the scalar firstorder jamiton ODE
(17) 
where the two functions and are defined as
Because and , the denominator in (17) has exactly one root, the sonic value (occurring at the sonic point), such that . The ODE (17) can be integrated through if the numerator in (17) has a simple root at as well. This leads to the ChapmanJouguet condition [18]
which yields a relationship between the constants and as follows:
One therefore has a oneparameter family of smooth traveling wave solutions, parameterized by , each being solutions of (17).
Into these smooth profiles shocks can be inserted that move with the same speed . The first condition in (13) implies that the quantity is conserved across the shock (in addition to being conserved along the smooth parts by (16)). And both conditions in (13) together imply that is conserved across shocks. Hence, when integrating (17), one can at any value insert a shock that jumps to a value with and continue integrating (17) from there. Moreover, for those shocks to satisfy the Lax entropy conditions [64], one can only jump downwards, i.e., . This, in turn requires that the smooth jamiton profile must be an increasing function. Using L’Hôpital’s rule in (17) at the sonic point yields that
which means exactly that the SCC is violated. In other words, as shown in [55], jamiton profiles with shocks can exist if and only if the SCC is violated.
The construction in Eulerian variables is analogous, albeit a bit more technical (cf. [19]). The traveling wave ansatz leads to
Integrating the first equation yields , which allows one to substitute via and vice versa. The second equation becomes the jamiton ODE
where . The ChapmanJouguet condition (matching roots of numerator and denominator) leads to the relations: and . Shock and entropy conditions are then implemented analogous to the Langrangian situation.
With these rules, jamiton solutions can be constructed (in either choice of variables). For a given choice of (and thus uniform propagation speed), any pattern of solutions to (17) connected by shocks (satisfying the above conditions) results in a traveling wave solution. The jamitons between any two shocks can be arbitrarily short (with a small variation around ), or may be arbitrarily long. In fact, it is not even required for the jamitons between shocks to have the same length (see [19, 55] for visualizations of jamiton profiles).
While all of these constitute feasible traveling wave solutions of the ARZ model (4), it does not mean that all such profiles would be dynamically stable under perturbations. In fact, both numerical evidence (see §3) as well as intuition dictate that neither very short, nor very long jamitons should be stable. The former because they can be thought of as a small (sawtooth) perturbation of the constant state (which is unstable because the SCC is violated, see §2.2); and the latter because their long tail will itself be close to a constant which, if that state violates the SCC, will be dynamically unstable. In other words, too short jamitons merge and have longer waves form between them; and long jamitons have new instabilities grow in their tails. It is only the middle range of jamitons (not too short and not too long) that is expected to be dynamically stable; and only those should arise in actual practice.
3. Computational Study of Jamiton Stability
To understand the dynamic stability of jamitons, we conduct a systematic study of the ARZ model (4) via direct numerical computation. After constructing a periodic jamiton as outlined in §2.3, we insert that profile as an initial condition into a numerical scheme (§3.1) and investigate whether the profile is maintained under small perturbations (§3.2).
3.1. Numerical scheme for the ARZ model with relaxation term
The ARZ model (4) is a system of hyperbolic conservation laws with a relaxation term. The hyperbolic part of the system can be solved using a finite volume scheme based on an approximate Riemann solver [40]. To find the numerical flux at the cell boundaries, we use the HLL approximate Riemann solver [23], which guarantees that the numerical fluxes satisfy the entropy condition [36]. Given the grid cell , where is the cell size, let
denote the approximate solution (cell average) in cell and the numerical flux at the boundary between cells and , respectively, at time (th time step).
A numerically robust treatment of the relaxation term is achieved by treating it implicitly, resulting in the semiimplicit update rule
This ensures stability even when is small. Note that, because the implicit term appears only in the equation and because it is linear in , the formally semiimplicit numerical scheme is actually fully explicit and the update step can be conducted in two substeps:

Update the component explicitly:

Now, with known from the first step, update
3.2. Results on the stability of jamitons
Using the numerical scheme described above, we conduct a computational investigation of the stability of jamitons (of the ARZ model (4) with the specific model functions and parameters described in §2.1
). Specifically, we classify the jamitons as follows: Evolve the solution up to some large final time, while regularly adding small perturbations. Then a jamiton is classified as stable if the jamiton profile is (within a tolerance) maintained at the final time, and unstable otherwise.
To classify a given jamiton (of length , with sonic density , upstream density , and speed ), we set up a periodic domain of length with initial conditions , i.e., the initial profile is four consecutive jamitons with shocks in between. We discretize using 10,000 grid cells, and run the numerical scheme (from §3.1) up to 3,000 (seconds; we omit units below).
During the numerical solution process, a small smooth perturbation is added to the vehicle velocity field in each step. The perturbation in the th step is
where the
are normally distributed random numbers with mean zero and standard deviation 1. As in the EulerMaruyama method, the additive noise is scaled with
. The value is chosen so that the highest frequency mode has a period that is not below the vehicle length , i.e.,. In other words, we have white noise exactly until the vehicle scale, which is wellresolved by the numerical scheme. Finally, the noise scale is
for , and for . The rationale for this larger initial “thermal noise” is, like in probabilistic optimization techniques, to make it easier for the solutions to escape their initial configuration in case it is only mildly unstable.Once the solution at is found, we first determine the number of shocks. If that number is not equal to 4, we immediately classify the jamiton as unstable. Otherwise, we check the jamiton speed by plotting the points for in the fundamental diagram (FD), and calculate as the least squares best fit slope of these data points (see [55] for the reason why is the slope in the FD). If , we classify as unstable. Otherwise, we classify as stable.
This process is now conducted (and run in parallel on a HPC cluster) for 980 different jamitons that are sampled as follows. First we sample 35 values of equidistant in the interval where the SCC is violated. Then, for each , we pick 28 values of in , where is the upstream density corresponding to the infinite jamiton [55].
The results of this classification are displayed in Fig. 2. Each of the four panels shows the same results, but in four different “phase planes”. Each jamiton is uniquely determined by two parameters: (i) the sonic density or equivalently the wave speed ; and (ii) the downstream shock density , or equivalently, the average density across the jamiton, or equivalently, the jamiton length . Panels 1(a), 1(c), and 1(d) have the on the horizontal axis, and , , and , respectively, on the vertical axis. Panel 1(b) displays vs. . In each quantity except , the jamiton region (where the SCC (11) is violated) spans an interval. The dashed brown curve corresponds the zerolength jamiton limit (in which ), while the solid dark blue curve represents the limit of infinitely long jamitons. Inside that jamiton domain, the 980 investigated jamitons are displayed as colored dots: stable jamitons are light blue; unstable jamitons are red. Note that the void regions visible in Panel 1(a) (top left), Panel 1(b) (bottom left), and Panel 1(d) (bottom right), also possess jamitons that were not simulated due to the sampling strategy of the 980 examples.
The results display intriguingly clear patterns: there appear to be two smooth curves inside the jamiton region that separate the stable from the unstable jamitons. Specifically, there are two unstable regions separated by a stable region: short jamitons which perturbations cause to coalesce into bigger ones (a “merging” instability); and long jamitons in which the long tail is linearly unstable and sheds growing waves (a “splitting” instability). This last characterization of these two mechanisms is based on observing the timeevolution of the computations, as well as the stability analysis below.
4. Stability Analysis of Jamiton Solutions
We now move towards a mathematical analysis of the dynamic stability of jamitons. For this we switch to the Langrangian variables introduced in §2.3. Consider a given jamiton , with sonic specific volume , and Lagrangian length (which is actually the number of vehicles in the jamiton) . We start by writing the (Lagrangrian) ARZ model (12) in the frame of reference of this jamiton, which has a propagation speed . Thus we introduce the variables (the same variable used in §2.3 to construct the jamitons) and the nondimensional time (for consistency with the scaling used for ). Because of that last choice, any instability growth rate computed with these variables needs to be scaled by to recover physical units.
In the coordinates defined above, equations (12) become
(18) 
This system is in conservative form, with conserved quantities and . The characteristic speeds of (18) are
(19) 
The RankineHugoniot shock jump conditions associated with (18) are
(20) 
where is the shock speed in the – frame. Contacts require and .
4.1. Perturbation system for singlejamiton waves
We now formulate a linear perturbation system of (18). There are two fundamental differences to the linear perturbation analysis for uniform flow presented in §2.2. First, because the jamiton profile is nonconstant, we obtain a variable coefficient linear system. Second, because the jamiton contains a shock, we must introduce a perturbation to the shock’s position as an additional variable (a variable not needed for perturbations of smooth solutions). As we will see below in more detail, both aspects render this analysis significantly more complicated than the one in §2.2.
Here we consider the stability of periodic jamiton profiles with one shock per period, under periodic perturbations. Note that this setup excludes the possibility of jamitons merging by means of adjacent shocks approaching each other. Hence, we only study the “splitting instability” for long jamitons, not the “merging instability” for short jamitons (see §3.2).
Consider a periodic jamiton profile , of length between shocks, and write it as , — a solution of (18) on with the shock placed at 0. Now write and , where and are infinitesimal perturbations. Substituting into (18) yields the linear system for and :
(21) 
We also need to track the infinitesimal perturbation of the shock position . We do so by implementing the RankineHugoniot conditions (20) in a way consistent with solving (18) on with periodic boundary conditions. This then generates boundary conditions for (21). The first equation in (20) yields
Expanding this equation, ignoring terms beyond , and using that and , we obtain
The second equation in (20) becomes
Again, ignoring terms beyond and using that , we get
In this setup the bracket notation denotes . Therefore, we have derived the following variablecoefficient linear model for and on , with boundary conditions that involve the shock position perturbation :
(22)  
(23) 
We conduct two further simplifications to the model. First, we transform it to characteristic form by writing it in terms of the Riemann variables and . Second, we replace the shock perturbation variable by a Robin b.c. for the PDE, as follows. Differentiating the boundary conditions with respect to time yields
Using the fact that , we obtain Robin boundary conditions for the PDE. Altogether, we obtain the following system
(24) 
with boundary condition
(25) 
The coefficients are computable from the jamiton functions as
where
4.2. Qualitative characterization of the jamiton perturbation system
We now adopt a short notation for the jamiton perturbation system (24), with b.c. (25), by writing and in place of of and , and introducing coefficient functions to obtain:
(26) 
with b.c. . The characteristic speed is constant and positive. In turn, vanishes at the sonic point , and is negative (positive) for (). Hence, the only ingoing characteristic is at , for (consistent with a single b.c.). The function crosses from negative to positive at as well, and it is always negative for ; it may or may not cross back to negative for some . Finally, everywhere. Figures 3 and 4 display the functions and characteristic curves, respectively, for an example jamiton.
Qualitatively, the solutions of (26) behave as follows. Being an advectionreaction system, its solutions are generally wavelike in nature. Waves enter the field at and are transported with the field to the right with constant speed , while being dampened by the term and modified (via the field) through the term. Likewise, the field constantly feeds into the field via the
Comments
There are no comments yet.