Deterministic particle approximation of the Hughes model in one space dimension

In this paper we present a new approach to the solution to a generalized version of Hughes' models for pedestrian movements based on a follow-the-leader many particle approximation. In particular, we provide a rigorous global existence result under a smallness assumption on the initial data ensuring that the trace of the solution along the turning curve is zero for all positive times. We also focus shortly on the approximation procedure for symmetric data and Riemann type data. Two different numerical approaches are adopted for the simulation of the model, namely the proposed particle method and a Godunov type scheme. Several numerical tests are presented, which are in agreement with the theoretical prediction.


Introduction
The mathematical modeling of human crowds has gained significant attention in the last decades, see [7,12,22,23,24,30,32,34,35]. In 2002, Roger L. Hughes [24] proposed a 'thinking fluid' approach to this subject, which consists in modeling a human crowd as a continuum medium, with Eulerian velocity computed via a nonlocal constitutive law of the overall distribution of pedestrians. Such a nonlocal dependence is encoded in a weighted distance function computed at a quasi-equilibrium regime, obeying an eikonal equation with right-hand side depending on the density of pedestrians. In its multidimensional formulation, the model reads Here x denotes the position variable, t is the time, and ρ = ρ(t, x) ≥ 0 represents the averaged crowd density. The model is typically posed on a bounded domain x ∈ Ω ⊂ R 2 , with boundary conditions depending on the position of the exits in ∂Ω, see [24,10] for further details. The map ρ → v(ρ) represents the absolute value of the velocity, and is assumed to be decreasing, as lower velocities correspond to higher densities. The function φ in (1) models the pedestrians perception of the exits, and is computed as a weighted distance encoding the overall distribution of the crowd. It may be interpreted as an estimated exit time for a given distribution of pedestrians. More precisely, pedestrians choose their direction upon evaluation of a running cost function ρ → c(ρ). Such a function is assumed to be increasing, since pedestrians choose their path toward the closest exit avoiding densely crowded regions.
Hughes' approach addresses specifically the movement of pedestrians in highly populated domains. In [25], the model has been suggested in particular to provide assistance in the annual Muslim Hajj and, quoting from [25], 'to locate barriers that actually increase the flow of pedestrians above that when there are no barriers present'. Such a phenomenon, analogous to a similar one known in vehicular traffic modeling, is often referred to in the literature as Braess' paradox for pedestrians, see [8,34]. Hughes' model is way far from being the only significant PDE-based approach to pedestrian movements, see e.g. [12,29]. Compared to other existing models in the literature, in addition to its potential applications listed in [25], this model has the undeniable advantage of combining classical overcrowding avoidance effects with singular direction-switching phenomena, which occur typically in a densely crowded region with limited information on the location of the exits. Clearly, the global perception of the domain implied by the eikonal equation in (1) is not totally realistic in most situations. However, the model is very flexible and can be easily improved to localize the range of perception, see [11]. Another argument in favour of (1) has been recently provided in [10], in which a full-transient version of this model has been justified as a mean field game.
In this paper we focus on the Hughes model [24] in its one-dimensional formulation coupled with homogeneous Dirichlet boundary conditions ρ(t, −1) = ρ(t, 1) = 0, φ(t, −1) = φ(t, 1) = 0, t > 0, and initial datum In this version of the model, two exits are located at the edges of the interval (−1, 1). It will be assumed that a maximal value of ρ is prescribed, which we set equal to 1 for simplicity. A preliminary formal maximum principle proved in [20] shows that ρ never exceeds the range [0, 1] ifρ takes values in [0, 1]. We shall assume what follows.
An important example is see [24,10] and Example 1.2 below. The boundary conditions (2b) are justified as follows. If we denote the efficiency β ∈ [0, 1] of the exits by the ratio between the escape velocity allowed by the exit and the maximum pedestrian velocity, we have that the flux of pedestrians through the exit is proportional to β, i.e. we obtain the Robin type condition condition where the term φ x /|φ x | gets cancelled when considering the outgoing normal direction of the flux. This is equivalent to imposing the exit velocity Since we assume that the velocity is related to the density via a monotonic function, then imposing the velocity is equivalent to imposing the density. Hence, at variance with what happens in diffusion problem in which Robin boundary conditions become in the limit either Dirichlet or Neumann condition, in our case we simply obtain a Dirichlet conditions on the density. For simplicity, we always considered the case of maximum efficiency β = 1, i.e. v = v max , which in turn gives zero boundary conditions on the density. Note however that in the cases in which the characteristics are exiting the domain, no boundary conditions have to be assigned on the boundary. A discussion about this issue is provided later on in this section. By a simple duality argument, a similar zero Dirichlet condition is obtained for φ, see the discussion in [10]. A fully satisfactory existence theory for (2) is still missing. A mathematical theory in this setting was first addressed in [18], in which the eikonal equation ∇φ = c(ρ) was replaced by two regularised versions involving a Laplacian term. In both cases, existence of weak solutions was proved, and the uniqueness of entropy solutions was provided in one of the two cases. [18] was the first rigorous attempt at the mathematical theory for the Hughes model. The use of a regularised eikonal equation implies a C 1 regularity of the velocity field in the continuity equation for the density ρ. Such a property does not hold in the non-regularized case (2a), and the numerical simulations performed in [18] showed evidence of discontinuities in the velocity field, expressing sudden changes of direction in the movement of the pedestrians.
A first attempt at a rigorous mathematical treatment of (2a) (without any regularization) was performed independently in [1] and [20]. In both papers, the authors provided a complete description of the behaviour of solutions for short times in the case of Riemann-type initial data on an interval. In [1] a solution semigroup was defined for short times for piecewise constant initial conditions. As a special case, solutions with direction-switching phenomena were produced analytically. In [20] the Riemann problem was analyzed more thoroughly, and the optimization of the evacuation time in that case was studied. Numerical simulations were performed in [18] and in [20], essentially based on Godunov schemes in both papers.
The results in [1] laid the basis for tackling the convergence of a wave-fronttracking scheme [9,15] for (2), see the numerical simulations in [21]. Such an approach was only partly successful, and led to the results in [2], which show existence of entropy solutions when the initial condition is well-separated, i.e. is yielding the formation of two distinct groups of pedestrians moving towards the two exits, with the emergence of a vacuum region in between, persisting until the total evacuation of the domain. A more general set of initial conditions was never covered so far, the main difficulty lying in the very large number of possible case studies in the interactions between classical and nonclassical waves in the wave-front-tracking algorithm.
As already observed in [1,2,20], the one-dimensional Hughes model (2a) can be reformulated in terms of the density ρ as follows: Here, ξ(t) plays the role of a turning point at which pedestrians switch their target exit in view of the overall crowd distribution. The reduction of (2a) to (4) is soon explained. In order to have a unique physically relevant (semi-concave) viscosity solution to the Dirichlet problem for φ, the derivative φ x can change its sign just once from positive to negative. The turning point ξ(t) is then defined as the point x at which x → φ(t, x) reaches its maximum. ξ depends on time because of the ρdepending term on the right hand side of the eikonal equation in (2a). Indeed, ξ(t) depends non-locally on ρ(t, ·). The second equation in (4) is recovered after a simple integration on the eikonal equation in (2a), upon the assumption that x → φ(t, x) is continuous at ξ(t). Notice that the flux F is possibly discontinuous along the turning curve x = ξ(t). For a detailed study of a Rankine-Hugoniot condition along the turning point we refer to [1,20].
Example 1.1. The simplest choice for the cost function is c ≡ 1. In this case the obvious solution to the eikonal equation in (2a) is the distance function φ(t, x) = dist(x, {−1, 1}) for all t ≥ 0. As a consequence, the ratio φ x /|φ x | is equal to −1 for x < 0 and to 1 for x > 0, and the model is equivalent to the scalar conservation law with discontinuous flux ρ t + sign(x) ρ v(ρ) As a result, Hughes' model reduces to the Lighthill-Whitham-Richards (LWR) model for traffic flow [28,33] with negative velocity on (−1, 0) and positive velocity on (0, 1). Roughly speaking, pedestrians are moving toward the closest exit regardless of the overall distribution. This is the typical behaviour of a crowd in a panic state, see [20].
see [1,2,10,18,24,25,36]. In this case, in order to bypass the technical issue of c blowing up at ρ = 1 in the eikonal equation of (2a), the initial conditionρ is assumed to take values in [0, 1 − δ], for some constant δ > 0 sufficiently small. This assumption, together with the maximum principle, ensures that the cost (5) computed along any solution of (2) is well defined. Rewriting Hughes' model in the form (4) highlights the non-local mechanism by which the overall distribution of pedestrians in the domain biases the choice of the exit path for each single pedestrian. Since c is increasing w.r.t. ρ, such a mechanism penalises regions with high density in the choice of the nearest exit. In the relevant case in (5), c(ρ) gets very large when ρ is close to 1, and the weighted distance function φ gets very steep in this case.
The reformulation of (2a) in the form (4) clearly suggests that Hughes' model can be seen as a two-sided first order LWR model for traffic flow [28,33], with the turning point ξ(t) splitting the whole interval (−1, 1) into two subintervals. Hinted by the new particle approach developed in [17], in which nonnegative entropy solutions to nonlinear scalar conservation laws ρ t + [ρ v(ρ)] x = 0 are rigorously approximated by the empirical measure of the follow-the-leader ODE particle system, in which particle i satisfies an equation of the forṁ where m denotes the mass assigned to each particle, in the present paper we propose a new deterministic particle approach to (4), which can be roughly formulated as follows. Here we just outline our particle method and refer to section 2.1 for its detailed application to approximating the solution of (4) via atomization. For a fixed N ∈ N and a fixed mass m > 0, let −1 <x 0 < . . . <x N < 1 be the initial positions of N + 1 particles in (−1, 1). The initial positionξ of the discrete turning point is determined uniquely by the identity which replaces the integral identity in (4), with the density ρ replaced by a discretized particle density, and with the assumption that a vacuum region is placed aroundξ and near the boundary. We observe that here we assume zero density around the turning point and outside the interval [x 0 ,x N ], which yields (6) in view of the assumption c(0) = 1. The pair (I 0 ,ξ) is uniquely determined by condition (6) andx I0 ≤ξ <x I0+1 . Assuming for simplicity thatx I0 =ξ, i.e. thatx I0 <ξ <x I0+1 , we let the particlesx 0 , . . . ,x I0 move towards the left exit via a backward follow-the-leader system, and the particlesx I0+1 , . . . ,x N move towards the right exit via a forward one. More precisely, with x i (0) =x i for all i = 0, . . . , N . The discrete turning point ξ(t) is required to obey ξ(0) =ξ and the identity The above formula is used as long as all points {x i } N i=0 belong to the interval [−1, 1], and is suitably modified by removing the particles leaving the domain and adjusting the length of the first and last interval (see the subsection 2.1 for the details). Clearly, system (7a) must be restarted once ξ(t) crosses one of the two neighbor particles x I0 (t), x I0+1 (t). In this case, a detailed analysis comparing the speed of the neighbor particle withξ at the crossing time is needed in order to define the post-crossing evolution consistently. We shall refer to system (7) as the FTL-Hughes model.
The goal of this paper is twofold.
• First of all we use the FTL-Hughes model (7) as a way to prove existence of weak solutions to the continuum model (2). We perform this task under analogous assumptions to those required in [2], namely in situations in which a vacuum region is generated around x = ξ(t), but with a lighter proof and under more general assumptions on the cost function c (only the cost function c(ρ) = 1/v(ρ) is covered in [2]). As a byproduct of our result, we obtain a rigorous deterministic many-particle approximation result for a class of solutions to (2). • As a second task, we implement the FTL-Hughes scheme (7) numerically in more complicated situations yielding particles crossing the turning point. We shall compare these simulations to existing ones in the literature [18,20] based on Godunov's scheme, and show that our deterministic particle approach is able to capture mass transfer phenomena around the turning point.
Next we provide our notion of solution to (2) and state our main results. For notational simplicity, we introduce the maximal initial density ρ max ∈ (0, 1] such that c is well defined on [0, ρ max ]. Then the initial datumρ is assumed to be in Since we are dealing with zero Dirichlet boundary data, assuming strict concavity of ρ → f (ρ) = ρ v(ρ) on [0, 1], it is easily seen that the proper boundary condition set in [6] is Tr(ρ) ∈ [0,ρ], whereρ ∈ (0, 1) is uniquely determined by f (ρ) = 0, see condition (V) above. Such a condition expresses the fact that characteristics are pointing outward the domain (both in x = −1 and in x = 1) if and only if Tr(ρ) ∈ [0,ρ] on those two points. It is well known (see [19]) that the boundary problem with Dirichlet conditions for a nonlinear conservation law is solved by considering a Riemann problem involving the boundary datum and the trace Tr(ρ) at the boundary. Since the boundary datum is zero and since the flux is concave, the Riemann solver at the boundary only consists of rarefaction waves, which leave the domain if Tr(ρ) <ρ and are centered on the value ρ =ρ at the boundary points if Tr(ρ) ≥ρ. Now, a solution of the Cauchy problem with a general L ∞ initial datum with support on [−1, 1] and without boundary conditions within the wavefront-tracking algorithm [9,15] easily shows that ρ achieves only values in [0,ρ) outside [−1, 1]. Therefore, the solution via wave-front-tracking without boundary data coincides with the one in which the boundary datum ρ = 0 is imposed every time a wave interacts with the boundary. Hence, by using the convergence results in [26], we recover that actually no boundary conditions have to be prescribed in our case, and the solution ρ to (2) can be constructed by coupling the integral relation in (4) with the continuity equation on the whole R restricted to the interval [−1, 1].
For this reason, also due to the fact that we shall only consider solutions with a vacuum region around the turning point, our notion of entropy solution will simply be obtained by 'merging' the two traffic equations on [−1, ξ(t)] and [ξ(t), 1] respectively. Observe that the strong traces of the solution at the boundary points exist due to the genuine non-linearity of the flux ( [31,37]) and must satisfy In this paper we only deal with boundary conditions with maximal efficiency β = 1 in (3). We shall address the case β < 1 in a future work. We also plan to impose (local and non-local) point constraints on the flow at the two boundary points to model the presence, for instance, of exit doors, see [3,4,5,12,13,14] or [34,Chapter 6] for more details.
We now state our notion of solution, which is motivated by the presence of a persistent vacuum region around the turning point for the class of solutions considered in this paper.
For any initial datumρ in S, there exists a unique entropy solution ρ to (2) in the sense of Definition 1.1, such that ρ(t, ·) ∈ S for all t > 0. Such a solution is obtained as a strong L 1 limit of the discrete density R(t, x) constructed in subsection 2.1 via the FTL-Hughes particle system (7).
In the next theorem we state our main result, which deals with a class of small initial data in BV. For further use, we define the function which is strictly decreasing in view of assumption (C) above. We then set Theorem 1.3. If ρ max < 1 and the initial datumρ ∈ BV((−1, 1) then there exists a unique entropy solution ρ to (2) in the sense of Definition 1.1 defined globally in time. Such a solution is obtained as a strong L 1 limit of the discrete density R(t, x) constructed in subsection 2.1 via the FTL-Hughes particle system (7).
We remark that the ρ max < 1 assumption above is essential in order to have the right-hand-side in the inequality (12) strictly positive.
The proofs of Theorems 1.2 and 1.3 are carried out in subsections 2.2 and 2.3 respectively.
The paper is structured as follows. In section 2 we cover the analytical theory. In particular, in subsection 2.1 we set up the deterministic particle scheme for (2), namely we construct the approximation to (4) via the FTL-Hughes model (7). In subsections 2.2 and 2.3 we prove our analytical results, with Theorem 1.3 being our main result. In section 3 we show the numerical simulations of our particle methods in simple Riemann-type initial conditions. In particular, we compare our tests with previous ones performed in [18,21]. The numerical section shows evidence that the two methods, though conceptually different, are in very good agreement when using a large number of particles. We stress here that, although the analytical results concerning our deterministic particle method are restricted to cases in which particle separate into two distinct sets keeping the same direction for all times, the numerical simulations also cover cases with direction switching, showing that the particle approach works also in those cases. Our study is restricted to the one-dimensional case. Possible extensions to multidimensional cases rely on the extension of the results in [17] to multi dimensional scalar conservation laws, which is not available in the literature. We set where spt denotes the support. We recursively definē The above defines a set of N + 1 particles −1 ≤x 0 <x 1 < . . . <x N −1 <x N ≤ 1 with the property In particular, this implies that For future reference, we denote the local discrete initial densities and introduce the discretized initial densityR : We now define the initial approximated turning pointξ n via the formula Clearly, the above formula definesξ n ∈ (−1, 1) uniquely. The next step is the definition of the evolving particle scheme. Roughly speaking, ξ n splits the set of particles into left and right particles, the former moving according to a backward follow-the-leader scheme, the latter according to a forward one. By a slight modification of the initial condition, we may always assume thatξ n does not coincide with any of the particles. Then, there exists I 0 ∈ {0, . . . , N } such that ξ n ∈ (x I0 ,x I0+1 ). We then set We consider the corresponding discrete density Notice that the above density has been set to equal zero outside the particle region [x 0 (t), x N (t)) and around the turning point, namely in [x I0 (t), x I0+1 (t)). The latter in particular is simply due to a consistency with the numerical simulations, in which the computation of the turning point is made simpler in this way. This simplifying assumption will introduce a small error m = 2 −n M in the total mass. In view of the above notation, (15) can be re-written as Notice that R I0+1/2 does not bias the movement of any of the particles, and this is an argument in favour of the ansatz R I0+1/2 (t) = 0. For future reference, we compute for any i ∈ {0, . . . , .
In particular, the time derivative of the discrete density R i+1/2 satisfies The (unique) solution to the system (15) is well defined until the turning point does not collide with a particle. Similarly, the expressions (18) for the time derivative of the discrete densities only hold until the first collision of a particle with the turning point. We note that the density R I0+1/2 (t) is equal to zero until the turning point collides with a particle.
In view of the discussion before Definition 1.1 on the solution of the continuum problem at the boundary, we shall not impose any boundary condition to the particle system (15), and we shall follow the movement of each particle whether or not they are in [−1, 1]. Hence, the discrete densities R i+1/2 (t) are (in principle) defined for all t > 0 and for all i ∈ {0, . . . , N − 1}. We remark that a careful analysis of the many particle approximation for the Dirichlet boundary value problem for a scalar conservation law is performed in [16].
The approximated turning point ξ n (t) is implicitly uniquely defined by c(R(t, y)) dy, where R : R + × R → [0, ρ max ] is the discretized density defined by Clearly ξ n (t) belongs to (−1, 1) for any t ≥ 0. We underline that ξ n (0) does not necessarily coincide withξ n . We compute now the time derivative of ξ n (t).
Proof. By (19) we have that Take the time derivative of the above equatioṅ (17) and (16) we obtain Finally, by applying (15) we conclude the proof.

2.2.
Proof of Theorem 1.2. Fix an initial datumρ in S. In [2] it is proven that the function ρ ∈ S defined for x ∈ [0, 1) as the entropy solution to the initial-boundary value problem if x ∈ (0, 1), is an entropy solution of (2) in the sense of Definition 1.1, with ξ ≡ 0. Moreover, as already pointed out in the introduction, the entropy solution to the above initialboundary problem coincides with the entropy solution to the Cauchy problem whereρ e is obtained fromρ by extending it to zero outside (−1, 1). Now, adopting the atomization procedure described in subsection 2.1, it is clear that the discrete turning point will be located at zero for all times due to the symmetry of the particle system. Therefore, the particles will split into two sets which will not change in time, with the two particles nearest the turning point getting further and further away from each other. Hence, by the results in [17], we obtain that the FTL-Hughes model (7) converges to the entropy solution of the above Cauchy problem as m goes to zero and N to infinity, and this concludes the proof. Proposition 2.1. Let L be as defined in (10) and C be as defined in (11).
then no particles change direction and the problem (15) admits a global-in-time solution.
Proof. Directly from Lemma 2.1 we deduce that Indeed, it is easy to obtain from the estimate proved in Lemma 2.1 that Moreover, by the discrete maximum principle proven in [17,Lemma 1], as long as no collisions occur between a particle and the turning point, we have thaṫ In conclusion we proved that the particles do not belong to the cone where ξ n evolves. As a consequence, no particle crosses the turning curve, namely no particle changes direction. Hence problem (15) can be split in two follow-the-leader problems, for which global existence follows from [17, Lemma 1]. Moreover, the BV contraction estimate proven in [17,Proposition 5] implies that the TV c (R) ≤ L TV(ρ), and this concludes the proof.
In order to prove that ρ is the unique entropy solution to (2) according to Definition 1.1, we only need to prove that the initial conditionρ is achieved as t 0 at least in a weak measure sense. Now, it is immediately seen that R(0, x) only differs from the initial atomizationR(x) on the interval [x I0 ,x I0+1 ). Moreover, recalling the definition ofR in (13) and of R in (20), we have x I 0R (y) dy and the last term above clearly tends to zero as n → +∞. SinceR(x) converges tō ρ in the sense of measures as n → +∞, we have that also R(0, x) approachesρ as n → +∞. This concludes the proof of Theorem 1.3.

2.4.
Particle approximations of Riemann data. In this subsection we consider the particle approximation of an initial condition with ρ − < ρ + . In order to simplify the computations, we first fix the number of particles [−1, 0) to be equal to N − , and the number of particles on (0, 1] equal to N + , with N + > N − . We call N . = N + + N − . We then fix the particle mass m > 0 such that m < 1/N + , and obtain the two Riemann values We set the initial N + 1 particle positions as follows The initial position of the discrete turning point is computed as In order to avoid the situation in whichξ coincides with the initial position of a particle, we remove that particle from the system in such a case. Lemma 2.1 implieṡ The assumptions on c and v clearly imply thatξ(0) < 0. Therefore, the avoidance of a collision between the turning point ξ(t) and a neighbor particle is guaranteed under the assumption It is immediately seen that condition (23) is satisfied when ρ + = ρ − and ρ + < 1.

Numerical simulations
In this section we present some numerical simulations of the particle method described above, and compare the tests with classical methods available in the literature. More precisely, we solve the particle system (15) using the Runge-Kutta MATLAB solver ODE23 with initial mesh sizes automatically determined by the total number of particles N and the initial density values. We then compare the particle scheme simulations with solutions of (2a) obtained via Godunov scheme.
An important remark has to be stated about the boundary conditions. As described in section 2.1, we do not impose any boundary condition on the particle method. The two leading particles x 0 and x N move with maximal velocity towards the opposite directions. For the Godunov method the boundary conditions are assigned as follows. We create two extra ghost cells, one just at the left of −1 and one juts at the right of 1. In all our numerical computation we set the value of ρ in those cells to zero, to mimic 'perfect exits'. In practice, any value that results in characteristic speeds pointing out of the domain will provide the solution, since the upwind scheme will not make use of such values (except to check that the characteristic actually point out of the domain). This procedure is used for convenience, since it simplifies the implementation of boundary conditions a lot.
Particular attention is devoted to the turning point evolution in the particle simulations, obtained by discretizing (19). Since no boundary conditions are imposed for the particle method, particles are free to exit the domain following the evolution of the two leaders, whereas only the particles still inside the domain bias the evolution of the turning point.
We test the atomization algorithm described in section 2.1 considering some examples introduced in the literature for various initial dataρ, see [18] and [20]. In each of the example reported, the choice for the cost function is , v(ρ) = 1 − ρ, and we show time evolution of the discrete density R(t, x) (20) in the domain Ω = (−1, 1). In Figure 1 and Figure 2 we consider two constant initial conditions ρ(x) = 0.25 andρ(x) = 0.6 respectively. In the latter case, as we expect, the simulation shows the formation of two rarefaction waves centered at x = −1 and x = 1, with trace R(t, ±1) ∼ 1/2. In the former case the rarefaction waves exit the boundary, hence we see three constant states inside Ω. Figure 1. Evolution of R(t, x) with initial dataρ(x) = 0.25 at times t = 0, t = 0.5 and t = 1. In the particle simulations the blu dots represent particles positions, whereas the red line is the discretized density. The magenta vertical line describes the turning point evolution.
In Figure 3 we consider the Riemann problem with initial condition for which the condition (23) in section 2.4 is satisfied. In this case, a vacuum region is formed around the turning point. In Figure 4 we consider the initial datum for which condition (23) is not satisfied. In this case, the turning point collides with its left neighbour particle. In order to compare our method with the tests performed in [18,21], in Figure 5 we consider the three-step initial condition As shown in Figure 4 and Figure 6, examples (25) and (26), exhibit the typical mass transfer phenomenon occurring when the turning point is not surrounded by a vacuum region. In such a case, particles are crossing the turning point ξ(t), and a non-classical shock is formed around ξ(t), see [1,Remark 5]. These examples show that the particle method introduced in section 2.1 has good chances to rigorously approximate the continuum solution to Hughes' model under less restrictive assumptions than the ones considered in theorems 1.2 and 1.3. Indeed, a key requirement in order to have the particle scheme in section 2.1 well posed and consistent in the limit is the avoidance of wild oscillations of the particle trajectories around the turning point. Such a behaviour is ensured in the examples (25) and (26), and we believe that it can be obtained in a less restrictive setting for the initial conditions than the one we imposed in Theorems 1.2 and 1.3. In support of this conjecture we performed several simulations with various initial conditions, and there is no evidence of such strange behavior. Of course a more detailed analysis is needed to rigorously prove convergence of the scheme.
In all the examples we set N = 200 particles and plot the particle positions and the discrete densities. We highlight the interaction with the boundary and the evolution of the turning point. (1) Given ρ, solve the eikonal equation.
(2) Given the solution to the eikonal equation, solve the non-linear scalar conservation law using the Godunov, Lax-Friedrichs or ENO schemes (see [18] and [20]).
A comparison betweeen the first-order particle method and a classical Godunov scheme is showed in Figure 7, where we plot the solution of the previous simulations (at time t = 1) both with the FTL-Hughes particle method and with the Godunov scheme. We set the spatial discretization according to the number of particles N = 200 and the time discretization with ∆t = 5 × 10 −2 . Note that the time step is the same for both methods and is selected so as to respect the CFL condition for the Godunov method. Empirically, the observed time step restriction for the FTL-Hughes method (7) is much less severe than for the Godunov method applied to the Eulerian descriprion of the flow. However a detailed stability analysis of the method is beyond the scope of the present paper, and will be subject of future investigation. It is evident from Figure 8 that the two methods, though conceptually different, produce solutions which differ by a very small error. A more refined analysis of the error will be produced in a future work.

Conclusions
The FTL-Hughes particle scheme introduced in section 2.1 has been proven to rigorously converge to weak entropy solutions to Hughes' model in 1d with zero Dirichlet conditions (2) in cases for which a vacuum region is generated around the turning curve x = ξ(t) for all times, i.e. no mass transfer occurs through the turning point ξ(t). At the particle level, this means that particles initially split into two sets moving towards distinct directions, with no change of direction for all times.
The numerical simulations in section 3 supported this result, and show evidence that the proposed particle method is consistent also in cases in which particles do cross the turning point thus switching direction. The observation of such a phenomenon motivates further studies on the rigorous analytical study of the particle method, in particular to show (for the first time) an existence result for (2) with large initial data. Even the global-in-time well-posedness of the particle scheme is difficult in such cases. Indeed, even in some simple symmetric cases the continuation of the particle scheme after collision with the turning point cannot be prescribed by actually crossing the turning point itself. The authors will continue the study of this problem in a future work. equazioni alle derivate parziali nella matematica applicata. GR was partially supported by ITN-ETN Marie Curie Actions ModCompShock -'Modelling and Computation of Shocks and Interfaces'. MDR acknowledges hospitality by the Gran Sasso Science Institute in L'Aquila in a period in which this manuscript was being written.