The BGK approximation of kinetic models for traffic

We study spatially non-homogeneous kinetic models for vehicular traffic flow. Classical formulations, as for instance the BGK equation, lead to unconditionally unstable solutions in the congested regime of traffic. We address this issue by deriving a modified formulation of the BGK-type equation. The new kinetic model allows to reproduce conditionally stable non-equilibrium phenomena in traffic flow. In particular, stop and go waves appear as bounded backward propagating signals occurring in bounded regimes of the density where the model is unstable. The BGK-type model introduced here also offers the mesoscopic description between the microscopic follow-the-leader model and the macroscopic Aw-Rascle and Zhang model.


Introduction
There are mainly three modeling scales in the mathematical description of vehicular traffic flow. The microscopic one is based on the prediction of the trajectory of each single vehicle via systems of ordinary differential equations. The macroscopic one is based on the assumption that traffic flow behaves like a fluid where individual vehicles cannot be identified. Therefore, the flow is represented by a density function and evolved in space and time via partial differential equations. The last one is the mesoscopic scale which offers an intermediate description of traffic between the microscopic and the macroscopic scale. In fact, mesoscopic or kinetic models are characterized by a statistical description of the microscopic states of vehicles but, at the same time, still provide the macroscopic aggregate representation of traffic flow e.g. linking collective dynamics to pairwise interactions among vehicles at a smaller microscopic scale.
Prototype examples of microscopic models are given by the optimal velocity model by Bando [3] and the follow-the-leader model [15]. The first and classical example of a macroscopic model for traffic was introduced by Lighthill and Whitham [30] and independently by Richards [44]. This model is based on the single conservation law for the density assuming that the dynamics is governed by a velocity function which is a function of the density only. An extension of this is provided by the second-order macroscopic model introduced by Aw and Rascle [2] and by Zhang [48]. It is based on a system of two partial differential equations describing also the evolution of the speed of the flow. There is a proved link between common microscopic and macroscopic models. More precisely, the follow-the-leader model approaches to the Lighthill-Wihtham-Richards model and to the Aw-Rascle-Zhang model in the limit of infinitely many vehicles [1,10,11,20,21].
This work is mainly concerned with the analysis of traffic flow models at kinetic level. However, we point-out that throughout the paper the study is performed by linking models at different scales. Classical kinetic models for vehicular traffic [26,27,37,40,41] are based on a Boltzmann-type equation that, in classical kinetic theory, describes the statistical behavior of a system of particles in a gas. Other type of kinetic formulations are available for modeling traffic flow and they also allow to simplify the complex structure of the integro-differential Boltzmann-type equation. Among them, we mention Vlasov-Fokker-Planck type models in which the integral structure of the collision kernel is replaced by differential operators, obtained also by means of suitable time scaling, e.g. see [18,22], and discrete-velocity models, e.g. see [9,12,19].
We study spatially non-homogeneous Boltzmann-type kinetic models for traffic and in particular the BGK approximation, originally introduced by Bhatnagar, Gross and Krook [4] for mesoscopic models of gas particles. The BGK equation replaces the Boltzmann-type kernel with a relaxation towards the equilibrium distribution of the full kinetic equation. For this reason, the BGK model provides a good description of the Boltzmann-type equation in the limit of many interactions, since in this regime the kinetic distribution is close to equilibrium. However, when the interaction frequency is extremely high to immediately drive the kinetic distribution towards the equilibrium one, the kinetic equation converges to the solution of the Lighthill-Whitham-Richards model. Since the flow is at equilibrium, non-equilibrium phenomena are not described in this regime. If the interaction frequency is still large but there is some residual between the kinetic and the equilibrium distribution, then the kinetic model is able to reproduce non-equilibrium waves. We are interested in the study of the stability of the BGK model in the latter case. Using a Chapman-Enskog expansion we show that the BGK equation for traffic flow problems yields an advection-diffusion equation having a negative diffusion coefficient in dense traffic. This analysis is closely related to the investigation of the sub-characteristic condition for relaxation systems [6,23]. This implies produces possible unbounded growth of non-equilibrium waves. This behavior is investigated at the particle description and at the hydrodynamic limit of the BGK model.
Due to the possible instability of the kinetic model, here we aim to derive a novel BGK equation which is stable or weakly-stable in dense traffic. As in Section 3.2, we define the BGK model being stable if the sub-characteristic condition is always verified and weakly-stable if the sub-characteristic condition is violated in a proper subdomain of the admissible values of the density. As discussed in [46], a weakly-stable model for traffic is not necessarily unrealistic: the instability occurring in the regime can be regarded as model for the appearance of moving unstable but bounded waves, as e.g. stop-and-go waves. The boundness of these instabilities is numerically observed in Section 3.3. In [46] the analysis of the sub-characteristic condition has been performed for the Aw-Rascle and Zhang model, showing that it could be either satisfied or weakly-satisfied, depending on the choices of the equilibrium speed or flux and of the pressure function. Based on this study, our goal is to develop a BGK-type equation which has as hydrodynamic limit the Aw-Rascle and Zhang model so that it may be at least weakly-stable. The new kinetic model for traffic is derived by starting from the follow-the-leader model that can be seen as a discretization of the Aw-Rascle and Zhang model. Therefore, we do not only introduce a well-posed or weakly wellposed kinetic equation but we also contribute in deriving the mesoscopic description of common models for traffic flow. The diagram in Figure 1.1 schematically synthesizes the contribution of this work.
The manuscript is organized as follows. In Section 2 we briefly review the homogeneous kinetic model for traffic flow introduced in [43]. In Section 3 we introduce the non-homogeneous version of the kinetic model and we study non-equilibrium phenomena by using the classical BGK approximation of the full kinetic equation. We show that model is unstable, namely it is able to produce unbounded growth of small perturbations, via Chapman-Enskog expansions and hydrodynamic limits. In Section 4 we propose a way to derive an alternative formulation

Microscopic scale
Mesoscopic scale Macroscopic scale Spatially homogeneous model (Section 2): The solution of Q[f, f ] = 0 provides the distribution at equilibrium Mf (v; ρ) Classical BGK model (Section 3.2): Sub-characteristic condition Payne-Whitham type model (Section 3.2): Sub-characteristic condition U eq (ρ) 2 ρ 2 < 0 never satisfied Bando's model (Section 3.2): It does not consider the interaction term and it is unstable if > 1 2|U eq (ρ)| Follow-the-leader model (Section 4.1): Modified BGK model (Section 4.2): The sub-characteristic condition is conditionally verified in dense traffic Aw-Rascle and Zhang model (Section 3.3): of a BGK-type model for traffic flow. This approach has also the advantage of prescribing the mesoscopic step between the microscopic follow-the-leader model and the macroscopic Aw-Rascle and Zhang model. Finally, we discuss the results and possible perspectives in Section 5.

Review of a homogeneous kinetic model
Recent homogeneous kinetic model for traffic flow have been discussed e.g. in [43,45]. The model proposed therein is used as starting point for the discussion in the following. However, the particular choice below is not required for the shown results later. General kinetic traffic models of the form read is the mass distribution function of the flow, i.e. f (t, x, v)dxdv gives the number of vehicles in [x, x + dx] with velocity in [v, v + dv] at time t > 0. Moments of the kinetic distribution function f allow to define the usual macroscopic quantities for traffic flow yielding density, flux and mean speed of vehicles, respectively, at time t and position x. The space V := [0, V M ] is the space of microscopic speeds. We suppose that V is limited by a maximum speed V M > 0, which may depend on several factors, and that the density is also limited by a maximum density ρ M . Throughout the work we will take dimensionless quantities, thus V M = 1 and ρ M = 1.
In [43] the spatially homogeneous kinetic model corresponding to (2.1) has been discussed: where we consider binary interactions, in which a test vehicle with velocity v * interacts with a field vehicle with speed v * , emerging with speed v as a result of the interaction and with probability given by P. The kinetic model for traffic flow studied in [43] can be generalized to the following probability distribution accounting for stochastic speed transitions, see [45]: is a decreasing function of the density giving the probability of accelerating, ∆ a and ∆ b are the acceleration and the braking parameters, respectively. We can think of ∆ a as the instantaneous physical acceleration of a vehicle. The parameter ∆ b instead corresponds to an uncertainty in the estimate of the field vehicle speed, either one's own, or the one of the vehicle ahead, as in the first and the second line respectively of (2.6). Note that the model is continuous across the line v * = v * and that mass conservation holds true since In the following we will consider ∆ b = 0 and ∆ a = ∆v constant, so that one recovers exactly the model in [43].
In [43] a result on the existence and uniqueness of a particular class of stationary solutions has been established: Modeling with kinetic equations such as (2.1) is driven by the interaction kernel Q[f, f ] and for this reason homogeneous kinetic models are widely studied for traffic flow problems. Using the steady-state of the kinetic model it is possible to define the flux and the mean speed of vehicles at equilibrium as Notice that Theorem 2.1 states that M f depends only on a finite number of speeds, that it is parameterized by the local density ρ and that, for each fixed value of ρ, it is unique. Therefore, the maps ρ → Q eq (ρ) and ρ → U eq (ρ) are actually functions prescribing a uniquely defined correspondence between the density and the quantities (2.7). These relations provide the so-called simulated fundamental diagrams of traffic which are used in order to validate the model. In Figure 2.1 we show the fundamental diagrams of homogeneous kinetic model (2.4)-(2.6) that are computed by using the equilibrium distribution prescribed by Theorem 2.1. The top row shows the correspondence ρ → Q eq (ρ) and the bottom row shows ρ → U eq (ρ). We take P (ρ) = 1 − ρ and we consider the case of N = 2 speeds (left panels) and N = 3 speeds (right panels) with ∆v = 1 N −1 . In Figure 2.2 we show that the fundamental diagrams can fit experimental flux diagrams. This is for instance the case of the experimental measurements provided by the Minnesota Department of Transportation in 2013 and tuning the parameters of the model as N = 3 speeds, ∆v = 1 2 and P (ρ) = 1 − ρ 1 4 . Observe that the flux at maximum density is not zero since we compute the diagram up to a value of the density ρ such that ρ < ρ M . This is due to the observation that experimental data contain a residual movement even in the congested phase. The model is able to describe the phase transition and to catch the two regimes of traffic compared with experimental data. For low values of ρ, more precisely for ρ < ρ crit , where ρ crit is the solution of P (ρ) = 1 2 , the flux is linear in ρ. This is the phase of free flow, when U eq (ρ) = V M . For larger values of ρ, i.e. ρ larger than the critical density ρ crit , the role of the interactions increases, and the flux decreases. This corresponds to the congested phase of traffic flow, see also [24]. Remark 2.2. Thanks to Theorem 2.1 we have that the equilibrium distribution M f has a fixed dependence on v, and depends on x and t only through the local density ρ(x, t). However, for ∆ b = 0, it was shown in [43] that unstable equilibria can be found, for which M f depends not only on ρ, but also on the initial distribution f (t = 0, x, v). Those equilibria are unstable, in the sense that they disappear and drive to the solutions given by Theorem 2.1 when the initial data are perturbed. If the braking uncertainty ∆ b > 0, it was shown in [45] that only stable equilibria are found, they still depend on a finite number of speeds. In this work we still focus on the case ∆ b = 0 and consider initial distributions f (t = 0, x, v) such that only stable equilibria occur.
3 Uncontrolled instabilities in standard non-homogeneous kinetic models

Motivation: non-equilibrium regimes
The properties of the fundamental diagrams derived in Section 2 are the result of the modeling of the microscopic interactions. Although the model is able to reproduce typical features of the experimental diagrams, we point-out that the simulated diagrams are only single valued. This is a direct consequence of the uniqueness of the equilibrium distribution provided by Theorem 2.1. This means that the model cannot reproduce and explain the scattering behavior of the experimental measurements, at least not at equilibrium. Mathematical models being able to describe multivalued fundamental diagrams were also studied and they link the observed scattered data to several reasons: for instance the presence of heterogeneous components, such as different drivers [31,33] or populations of vehicles [42], or the deviation of microscopic speeds at equilibrium [13], or the presence of non-equilibrium phenomena such as stop and go waves [25,46,49]. In this work we will focus on the last motivation and consider non-equilibrium phenomena that might lead to complex traffic phenomena as e.g. stop and go waves. Such phenomena cannot be described by spatially homogeneous models but there is the need to reconsider the space dependence and thus to take into account the full kinetic equation (2.1).

Non-homogeneous kinetic models
In the non-homogeneous case, the relaxation parameter plays a crucial role since it balances the weight between the convection and the source term. In analogy to gas dynamics, where the interaction frequency is given by the Knudsen number, different values of lead to different  regimes. The situation is synthesized in Table 3.1. If we allow for = 0, i.e. we suppose that the interactions are so frequent to instantaneously relax f to the local equilibrium distribution M f , we are in the so-called equilibrium flow regime where (2.1) is solving the conservation law for the density with Q eq (ρ) defined in (2.7). Equation (3.1) is a first-order macroscopic model for traffic flow [30,44]. Instead, we expect that if is small, but not vanishing, then we are in a regime where the kinetic equation (2.1) is solving the continuity equation (3.1) with a diffusive perturbation of order O( ). For 1 we are in the transitional regime where the counterpart of the kinetic equation is given by extended continuum hydrodynamic equations. In traffic flow, this is the case of the Aw-Rascle [2] and Zhang [48] model. For > 1, but not too large, we are in the kinetic regime or rarefied gas regime and finally for → ∞ we obtain the regime of the collision-less kinetic equation.
Here, we will consider a type of approximation for the full non-linear Boltzmann-type equation which formally holds for small values of , namely the BGK approximation. This technique is based on replacing the Boltzmann-type kernel in (2.1) by a relaxation term: where M f (v; ρ) is the distribution of the microscopic speeds v at equilibrium, corresponding to the density ρ at time t and position x derived from the spatially homogeneous equation (2.4). Equation (3.2) assumes that the equilibrium distribution is known. For small values of , the BGK model provides the same behavior of the full kinetic equation (2.1). In contrast, (3.2) is only an approximation of (2.1) in the regime with large.
Remark 3.1 (Backward propagating information). Let us consider fixed and 1, i.e. the regime in which the transport term in (3.2) is stronger then the relaxation term. It is well known that in this regime, since the velocities in traffic are non-negative, the information propagate in the wrong direction for congested traffic situations. In this phase, in fact, we expect to see backward propagation of signals but the transport part move them forward. This phenomenon is shared both by a full kinetic model (2.1) and the BGK equation (3.2). For this reason, several models with the ability of overcoming this drawback were introduced in the mathematical literature [12,27]. This effect could be linked to a fixed choice of . Instead, in analogy with the Knudsen number in gas models, should be a decreasing function of the density. In this way, is large in the free flow phase since the interactions are less frequent and the convective term rules the dynamics. On the contrary, since decreases when density increases, i.e. interactions among vehicles are dominant, we expect that relaxation is faster when ρ is high. It is clear from the fundamental diagrams that, in the regimes in which f is close to the local equilibrium M f , signals move backward in congested phases. See also the left panel of Figure 3.1 where we show the characteristic speeds of the kinetic model when f is close to equilibrium for several values of the speed number N . In other words we claim that the regimes described in Table 3.1 depend on the phase of traffic.
Numerical scheme. We briefly describe the numerical discretization of BGK using a first-order implicit-explicit (IMEX) scheme in time [34,39] by treating the convective term explicitly and the collision term implicitly. This allows to treat the stiffness for small values of . Moreover, observe that, due to the linearity of the source term with respect to f , the scheme remains still explicit.
Let us consider the computational domain (t, Here, the numerical parameter δv is a sub-multiple of the physical parameter ∆v defined in Section 2. In this way, N , i.e. the number of discrete speeds defining the equilibrium distribution M f , see Theorem 2.1, is a sub-multiple of N v . Although in the transient speeds that are not spaced by ∆v give a contribution to the dynamics, while at equilibrium only the N speeds spaced by ∆v give a non-zero contribution, here we always consider the simple case with δv = ∆v. For piece-wise constant approximation of cell-averages in space we for i = 1, . . . , N x , j = 1, . . . , N v and n ≥ 0. Observe that the explicit knowledge of the equilibrium distribution M f , given by Theorem 2.1, is crucial for the performance of the numerical scheme. The quantity F can be one of the known Lipschitz continuous, conservative and consistent numerical flux function. The asymptotic preserving [14,35] result holds for the scheme (3.3).
Proposition 3.2 (Asymptotic preserving). Scheme (3.3) is asymptotic preserving in the sense that it solves a global first-order discretization of the macroscopic equation (3.1) with the flux function Q eq (ρ) given by (2.7) when → 0.
This scheme is an explicit conservative, consistent and stable discretization of the macroscopic equation (3.1) with the flux function Q eq (ρ) given by (2.7), provided F(·, ·) is a conservative, consistent and stable numerical flux function for the non-linear conservation law (3.1). Hence, we choose the time-step δt using the CFL condition applied to the non-linear conservation law (3.1) and in view of the nonlinear transport on macroscopic level we employ the Lax-Friedrichs flux: Non-equilibrium effects. We show numerically that the multivaluedness in fundamental diagrams can be reproduced as a result of non-equilibrium effects. The top rows of Figure 3 The dots in the right bottom panel contain the evolution of the flux ρu defined in (2.3) which is drawn as a function of ρ. Note that, at equilibrium, the fundamental diagram Q eq (ρ) is a function of ρ. When the system is not in equilibrium, the flux ρu depends not only on ρ, but also on the distribution f . Again, the green dots are obtained in the initial stages of the simulation, and they fade into blue as the solution evolves. If the flow were in equilibrium, as in the case = 0, then f coincides with the equilibrium distribution M f (uniquely determined by ρ). Thus the dispersion in these plots is a measure of deviations from equilibrium. Hence, the dispersion of the points in blue may reflect the multivaluedness of the fundamental diagram In all simulations we consider N = N v = 3 speeds and so the physical acceleration is ∆v = 1 2 . The spatial domain is [0, 1] and discretized by N x = 400 cells. The final time is T M = 0.2. We consider an initial smooth perturbation in the density with periodic boundary conditions. The initial kinetic distribution f (t = 0, x, v) is the uniform distribution over the velocity space with density ρ 0 (x). In Figure 3.2, ρ min = 0.1 and ρ max = 0.3. The perturbation in the density is below the critical density ρ crit = 1 2 . Thus here the whole solution moves towards the right for = 0 because in this regime the flux is linear with respect to the density. In the case = 10 −1 the solution also moves to the right but with very little deformation. As consequence, the flux diagram shows a very small dispersion, which means that the flow is almost at equilibrium.
If we choose ρ min = 0.6 and ρ max = 0.9, the initial perturbation occurs on a dense traffic flow, see is no longer satisfied. These are purely non-equilibrium effects, as shown also in the flux diagram. The equilibrium fundamental diagram is shown with the red line, and it can be seen bisecting the ρu values computed during the evolution of the case = 10 −1 .
The last smooth test we show has a large initial pulse, with ρ min = 0.1 and ρ max = 0.9. The results are in Figure 3.4. The low density part of the wave on the right accelerates and propagates towards the right. The high density part of the wave propagates backward, increasing the slope of the profile. Thus a shock forms, propagating towards the left. This structure appears quickly in the case = 0 while it is slower for = 10 −1 where however we see the appearance of a new maximum. Again the solution with = 10 −1 exhibits non equilibrium waves which are represented by the scattered values in the flux diagram. Remark 3.3 (Solutions of Riemann problems). As already observed, in the relaxation limit the full kinetic equation (2.1) and the BGK approximation (3.2) solve the conservation law (3.1) with Q eq computed using the equilibrium distribution of the homogeneous kinetic model. We notice that the fundamental relation Q eq defined in (2.7) has two main properties.
1. It is continuous: in fact, from Theorem 2.1 it is clear that Thus, there is no need to treat the macroscopic model (3.1) with special techniques for conservation laws with discontinuous fluxes.
2. It is not convex: in fact, it is linear in the free flow regime (i.e. when ρ ≤ ρ crit ) and it is convex in the congested regime (i.e. when ρ > ρ crit ). See the right panel in Figure 3.1 which shows the behavior of the second derivative of Q eq for different values of the speed number.
The non-convexity of the flux-density diagram represents the most interesting property since it implies that the wave structure at equilibrium is non-trivial, because even a simple Riemann problem initial data will result in more than a single wave. For a more detailed discussion we refer to [5,29].
Unbounded instabilities. The results of the previous simulations show that non-homogeneous models are important to describe non-equilibrium phenomena. Note that, in the case of the perturbation in dense traffic, the non-equilibrium effects provides new maxima and minima that apparently grow in time. These results suggest that the appearance of stop and go waves on highways could be interpreted and reproduced as non-equilibrium effects.
Here, we prove that the BGK model (3.2) has the drawback that it could lead to unbounded unstable waves in the congested regime of traffic. In fact, let us consider the same setting of the simulations provided in Figure 3.3, namely a perturbation of the density in the congested regime. In Figure 3.5 we show the solutions for a larger final time with = 10 −1 in the left panel and = 10 −2 in the right panel. We stop the simulations as soon as the density profile exceeds the maximum density. We observe that the BGK model propagates the initial profile backward producing unstable waves that do not remain bounded. This phenomenon will be investigated below using the Chapman-Enskog expansion.
Consider the Chapman-Enskog expansion of the BGK equation (3.2) as first order perturbation in of f around the equilibrium distribution M f (v; ρ). This leads to the advection-diffusion equation where the diffusion coefficient µ(ρ) is given by  In fact, considering the expansion plugging the approximation of f into the BGK equation (3.2) and integrating with respect to the velocity we obtain If µ(ρ) < 0 then the advection-diffusion equation is ill-posed and therefore has solutions with unbounded growth even starting from small perturbations. In the case of the kinetic model (3.2), the sign of the diffusion coefficient depends on the equilibrium distribution M f . The request and, since Q eq (ρ) is the fundamental diagram at equilibrium, this condition requires that the square of the characteristic velocities is bounded by the variation of the kinetic energy in each regime. The condition (3.7) depends only on the equilibrium distribution M f which is given by the homogeneous kinetic model and is the result of the modeling of the microscopic interactions. In Figure 3.6 we show the sign of the diffusion coefficient (3.6) in the case of the equilibrium distribution computed in [43], and recalled in Section 2, with two and three speeds corresponding to ∆v = 1 and ∆v = 1 2 , respectively. We observe that µ(ρ) ≥ 0 if ρ ≤ ρ crit = 1 2 , while µ(ρ) < 0 if ρ > ρ crit = 1 2 . This result could explain the unbounded growth of perturbations in the density, numerically observed in Figure 3.5.
In the following result, we prove that the instability of the model does not depend on the choice of the equilibrium distribution but it occurs for any suitable equilibrium of kinetic traffic models.
More specifically, for the case of two microscopic speeds and in view of the previous proposition, the model can be seen as the kinetic relaxation system of Jin and Xin [23] that cannot satisfy the sub-characteristic condition since v > 0 and Q eq (ρ) < 0 in the congested regime.
Remark 3.5. We observe that conditions (3.8) do occur in traffic flow. In particular, the first one, Q eq (ρ) < 0, is always verified in the congested phase of traffic. The second condition instead requires that the variance of the microscopic speeds decreases as the density increases. This is a reasonable assumption in traffic since we expect that the freedom in choosing a microscopic velocity reduces when traffic becomes dense.
In view of the results provided by the Chapman-Enskog analysis, we state the following definition which establishes that the BGK model is unstable. In order to further investigate the instability of the BGK model when applied to traffic flow problems, we study its hydrodynamic limit. The second-order macroscopic model derived from the BGK equation (3.2) is obtained without closing the continuity equation with the distribution at equilibrium and computing the evolution equation for the first moment. Then we have System (3.9) is closed and thus solvable when a closure law for the second moment of the kinetic distribution f is prescribed.
1. If we consider the following approximation with h(ρ) being an increasing function of the density, then from (3.9) it is possible to recover the Payne and Whitham model [38] second-order model for traffic flow: It is well-know that model (3.10) suffers of the following drawback (see for instance the works by Daganzo [8]): the characteristic velocities of (3.10) are λ 1 = u − h (ρ) and λ 2 = u + h (ρ). Thus one characteristic speed is always larger than the speed u of vehicles.

If we instead approximate
that is again the Payne and Whitham model without the anticipation term. However, the model still suffers of the following drawbacks. First, the system is not strictly hyperbolic since the eigenvalues are λ 1 ≡ λ 2 = u leading to a lack of strict hyperbolicity. The possible unbounded growth of small perturbations is highlighted by the fact that the sub-characteristic condition for the model is −ρ 2 U eq (ρ) > 0 which is never satisfied.
3. It is easy to show that if we approximate around the equilibrium distribution, i.e.
Remark 3.7. The particle dynamics resulting from the BGK model (3.2) are described by the ODE systemẋ where i is the label of a vehicle. System (3.12) represents the microscopic model by Bando [3] which is know to be stable provided < 1 2|U eq (ρ)| , i.e. if the relaxation is large.
The last results of this paragraph prove that properties of the BGK model studied via Chapman-Enskog expansion occurs also in its second-order hydrodynamic limit and in its particle formulation. Due to this consideration we aim to derive kinetic models representing the mesoscopic scale of suitable second-order macroscopic models. Therefore, in the following section we investigate properties of the second-order macroscopic model [2,48,46].

The case of the Aw-Rascle and Zhang model
The prototype example of second-order macroscopic traffic models is given by the Aw and Rascle [2] and Zhang [48] (ARZ) model where the conservation law for the density ρ(t, x) of vehicles at time t ≥ 0 and position x is coupled with an additional equation for the evolution of the macroscopic speed of the flow u(t, x). We focus on the non-homogeneous version of the ARZ model which writes in conservation form where q is the second conservative variable and it is defined as q := ρ(u + h(ρ)). The function h = h(ρ) is a strictly increasing function of the density and it is called hesitation function or traffic pressure although it is homogeneous to a velocity. The quantity is a time which rules the relaxation speed of the flux ρu to the equilibrium flux Q eq (ρ) which is a given function of the density, i.e. the fundamental diagram. Here Q eq is not necessarily the one given in (2.7). The ARZ model was introduced in order to overcome the drawbacks of the second-order macroscopic model introduced by Payne and Whitham and recalled in equation (3.10). The characteristic velocity of the ARZ model are λ 1 = u − ρh (ρ) ≤ u = λ 2 and thus information propagates at most with the vehicles' velocity.
If is small, but not vanishing, (3.14) approaches the advection-diffusion equation (3.5) where the diffusion coefficient µ(ρ) is given by In fact, using a Chapman-Enskog expansion, the coefficient (3.15) is obtained by considering a first-order expansion of the speed u = U eq (ρ)+ u 1 around the equilibrium velocity function U eq (ρ) Finally, plugging the expansion of u just obtained into the first equation of (3.14) we get the advection-diffusion equation (3.5) with coefficient µ(ρ) given by (3.15). The condition µ(ρ) > 0 provides the so-called sub-characteristic condition [6,23]. For the ARZ model µ(ρ) > 0 is satisfied if 0 > U eq (ρ) > −h (ρ). (3.16) The sub-characteristic condition guarantees the stability of solutions of the ARZ model based on initial conditions where ρ(t, x) = ρ and u = U eq (ρ). In fact, the LWR model (3.1) and the ARZ model (3.14) share solutions of the type (ρ, u) and in this case the sub-characteristic condition (3.16) can be rewritten in terms of the characteristic velocities of the two models computed in (ρ, u) as The above condition implies that the characteristic velocity of the LWR model must lie between the two characteristic velocities of the ARZ model. Whitham [47] proved that this inequality is verified if and only if solutions of the second-order macroscopic model with initial condition (ρ, u) are linearly stable, i.e. small perturbations to (ρ, u) decay in time [6]. We stress the fact that condition (3.16) strongly restricts the possible choice of U eq and h. For instance, the Greenshields' closure U eq (ρ) = 1 − ρ and the hesitation function h(ρ) = Cρ γ , C, γ > 0 given in [2] violate (3.16) if ρ γ−1 ≤ 1 Cγ . In Figure 3.7 we numerically analyze the sign of the diffusion coefficient µ(ρ) of the ARZ model given in (3.15) for several choices of the equilibrium velocity and the hesitation function. We observe that standard choices lead to regimes where the diffusion coefficient is negative in the congested regime and consequently the sub-characteristic condition is not verified. This produces  given in [1] with pressure h(ρ) = 2ρ.      instabilities of uniform flow as discussed above. However, for suitable choices of U eq and h, regimes where µ(ρ) < 0 can be properly contained in [0, ρ M ], guaranteeing the existence of a bounded density regime. This result was already mentioned and analyzed in [46] where the authors link the appearance of such regimes to the presence of nonlinear traveling wave solutions, called jamitons, and identified as the stop and go waves observed in real traffic situations.
Bounded waves. Numerical evidence suggests that if the instability of the model is restricted on a set (ρ 1 , ρ 2 ) properly contained in [0, ρ M ] guarantees that unstable waves occurring in this regime remain bounded. In fact, the growth of a wave is smoothed when the solution enters in the regime [0, ρ M ] \ (ρ 1 , ρ 2 ) where the diffusion is positive.
We point-out that this consideration is purely heuristic and based on numerical investigations. To the best of our knowledge there are few contributions on the analysis of advection-diffusion equations with a chancing sign diffusion coefficient, see e.g. [7,16,28,32].
In Figure 3.8 we numerically study the evolution of a density perturbation with the ARZ model as did for the BGK model in the previous section. Here we consider the case in which the bump in the density intersect a regime where the diffusion is negative. In particular, we consider the equilibrium flux Q eq and the equilibrium velocity U eq prescribed by the three-velocity kinetic distribution at equilibrium, see (2.7). The pressure or hesitation function is chosen as h(ρ) = 2ρ. This is the situation described in Figure 3.7f where µ(ρ) < 0 for ρ ∈ (0.5, 0.66). The initial condition is given by (3.4) with ρ min = 0.55 and ρ max = 0.9. The space domain is [−1, 1] and discretized with N x = 400 cells. The top-left panel shows the density profile at different instants in time. We observe that unstable waves appear in the regime where the diffusion coefficient is negative and they move backward. The top-right panel shows a zoom in the region where the unstable waves are rising and the bottom-left panels shows the evolution of the macroscopic speed. The growth is bounded as confirmed in the bottom-right panel where we show the density profile for a longer final time and using periodic boundary conditions.
In the following, we consider models that allow for a Chapman-Enskog expansion as advectiondiffusion equation with a negative diffusion coefficient in a subset (ρ 1 , ρ 2 ) properly contained in [0, ρ M ]. Those are not necessarily unrealistic, since the violation of the sub-characteristic condition lead to instabilities that in turn can be regarded as models for stop and go waves. These waves are however unbounded if the regime of negative diffusion is not properly contained in [0, ρ M ] as in the case of the classical BGK model.

Derivation of a modified BGK-type model for traffic flow
In this section we address the issue of deriving a modified BGK-type equation for traffic flow satisfying the following properties. The starting point for the derivation of the BGK-type model are microscopic follow-the-leader (FTL) models. The FTL model is proved to converge to the ARZ model in the macroscopic limit [1].
The kinetic model derived with this approach will also automatically guarantee Property 2. Moreover, since we have shown that the ARZ model may have a negative diffusion coefficient in a small density regime, see for instance Figure 3.7d-3.7f, we expect that the BGK-type equation will also satisfy Property 1.

Microscopic follow-the-leader model
In [1] a first-order semi-discretization in space of the Lagrangian formulation of the ARZ model is equivalent to the microscopic follow-the-leader model: and where x i and v i are the microscopic position and velocity of the vehicle i at time t ∈ R + , U eq (ρ i ) is a target speed depending on the value of the local density ρ i = ∆X xi+1−xi . The quantity ∆X is the typical length of a vehicle. In the following computations we assume ∆X normalized to 1. The constants C γ > 0, γ > 0 and the relaxation time are given parameters. The link between the microscopic and the macroscopic model is also studied in [10,11,17].
Following [1], we consider an alternative but equivalent formulation of the model (4.1) based on the quantity w i := v i + p(ρ i ). The function p = p(ρ i ) is the so-called traffic pressure and satisfies p(ρ) ≥ 0, p (ρ) > 0. In terms of the new variable w i , model (4.1) readṡ with K given by (4.2). Model (4.3) is identical to (4.1) but the introduction of w i allows us to rewrite the acceleration equation as a relaxation step.Note that the choice of the function p(ρ) is determined by the modeling of K. Since a different interaction term leads to different pressure functions, we keep consider the case of a general function K. Next, we introduce stochastic particle interactions in terms of the desired speed w i that is an approximation of (4.3). Then, we derive the corresponding kinetic equation for the evolution of the distribution of the desired velocity via the relaxation algorithm, see [36,Section 4.2.2]. Finally, we rewrite the kinetic equation in terms of the actual microscopic speeds of vehicles and we verify that Property 1 and 2 hold true for the new kinetic model for traffic.

A BGK-type model for traffic
Let g = g(t, x, w) : R + ×R×W → R + be the kinetic distribution function such that g(t, x, w)dxdw gives the number of vehicles at time t ∈ R + with position in [x, x + dx] ⊂ R and desired speed in [w, w + dw] ⊂ W . The desired speed w is assumed to be the speed that drivers want to keep in "optimal" situations. We define W := [w min , +∞) the space of the microscopic desired speeds where w min > 0 may be interpreted as the minimum speed limit in free-flow conditions.
The macroscopic density, i.e. the number of vehicles per unit length, at time t and position x is defined by and we define the macroscopic quantity q(t, x) := W w g(t, x, w)dw (4.5) In order to derive the evolution equation for the kinetic distribution g = g(t, x, w), we first prescribe the update of the microscopic states in a probabilistic interpretation that has a close link with the microscopic particle model (4.3). Let us assume a set of N particles (x 0 i , w 0 i ) is given at time t 0 characterized by the position x 0 i and the desired speed w 0 i . We can think of (x 0 i , w 0 i ) as samples from the distribution g(t = 0, x, w). The update of the microscopic states (x i , w i ) at time t 0 + ∆t is generated by the following algorithm.
Step 1. Reconstruct the density ρ i = 1 xi+1−xi at time t Step 3. For each i, with probability ∆t replace w 0 i with a velocity sample from a distribution M g with density ρ i ; otherwise set w i = w 0 i .
So far, the distribution M g = M g (w; ρ) only has to fulfill the requirement W M g (w; ρ)dw = ρ(t, x).
In fact, we observe that under this assumption the probabilistic microscopic algorithm prescribes an average speed which is equal to the output velocity prescribed by a first-order time discretization of (4.3).
According to [36,Section 4.2.2], we interpret the interactions introduced above as a stochastic particle method solving This equation is still a BGK-type equation since the collision kernel is linear and describes the relaxation of g towards a given distribution M g parameterized by the density ρ.
The deviation is controlled by p(ρ) in such a way that p(0) = 0 so that v = w, i.e. drivers can behave as they wish in the vacuum regime, and p(ρ M ) ≤ w min in order to ensure the bound v ≥ 0.
In terms of f = f (t, x, v) we define the density as 8) and the flux of vehicles as (ρu)(t, Moreover, for each fixed density ρ, we have that Hence, V vM f (v; ρ)dv = ρU eq (ρ).

Properties of the BGK-type model
We verify that Property 1 and Property 2 hold true for the kinetic equation (4.6).
Property 2 -Macroscopic limits. Due to the derivation from a FTL model this property should hold. We can also explicitly show that Property 2 holds by computing the hydrodynamic limit of (4.6). The continuity equation is obtained by integration in w space as where q is defined by (4.5). If one assumes that g is at equilibrium , g = M g , we obtain the classical first-order model for traffic flow Instead, the second-order moment equation gives the system x, w)dw − p(ρ)q(t, x) = 1 (ρU eq (ρ) + ρp(ρ) − q(t, x)) . System (4.12) is not closed since the knowledge of the second moment of g is required. In particular, the ARZ model (3.13) is recovered by closing as W w 2 g(t, x, w)dw ≈ q(t, x) 2 ρ(t, x) (4. 13) and taking the function p(ρ) = h(ρ).

Conclusions
In this paper we studied the stability of the classical BGK formulation of a kinetic model for traffic flow. Via Chapman-Enskog expansion we have shown that the BGK model leads to an advection-diffusion equation with a negative diffusion coefficient for a density greater than the critical one. This leads to unbounded growth of waves.
The issue has been addressed by deriving a modified BGK-type equation as mesoscopic description of the classical microscopic follow-the-leader model. The new model is of ARZ-type in the hydrodynamic limit and, as direct consequence, it result to be conditionally well-posed. The instabilities occurring in the regime of high densities have now a bounded growth and can be regarded as model for stop and go waves. Further, with the new kinetic model, we have also introduced the mesoscopic description between the microscopic follow-the-leader model and the macroscopic ARZ model.