Kinetic description of collision avoidance in pedestrian crowds by sidestepping

In this paper we study a kinetic model for pedestrians, who are assumed to adapt their motion towards a desired direction while avoiding collisions with others by stepping aside. These minimal microscopic interaction rules lead to complex emergent macroscopic phenomena, such as velocity alignment in unidirectional flows and lane or stripe formation in bidirectional flows. We start by discussing collision avoidance mechanisms at the microscopic scale, then we study the corresponding Boltzmann-type kinetic description and its hydrodynamic mean-field approximation in the grazing collision limit. In the spatially homogeneous case we prove directional alignment under specific conditions on the sidestepping rules for both the collisional and the mean-field model. In the spatially inhomogeneous case we illustrate, by means of various numerical experiments, the rich dynamics that the proposed model is able to reproduce.


Introduction
The complex dynamical behaviour of large pedestrian crowds has always fascinated researchers from various scientific fields. Academic studies began in earnest in the last century, starting with empirical observations in the early 1950s and continuing with the development of models in the field of applied physics. In recent years, applied mathematicians have become increasingly interested in analytic aspects as well as computational challenges related to simulation and calibration. Modelling the intricate individual dynamics and understanding the emergence of complex macroscopic phenomena, such as the formation of directional lanes or collective patterns, is an active area of research nowadays.
Different classes of models have been proposed in the literature -either at the micro-, mesoor macroscopic level. More recently also multiscale descriptions to model the impact of single individuals on the dynamics of a larger crowd have been discussed. Microscopic dynamics are based on the change of the individual position and possibly velocity due to interactions with others walkers and with the surrounding. The social force model [22], which is based on Newton's laws of motion, is among the most prominent and successful models in this class. But also cellular automata, see for example [9], or stochastic optimal control approaches, such as the one in [23], have been studied to model the individual behaviour. In kinetic models, the evolution of the crowd distribution with respect to the microscopic position and velocity of the pedestrians is described by a Boltzmann-type equation, in which interactions are included in the so-called "collision kernel"; see for example [19]. Macroscopic PDEs describe instead the crowd as a continuum with density, of which they study the evolution in space and time, see for example [7,14,36]. More recently multiscale approaches, which allow one to model the interactions of single individuals, for example leaders, with a large crowd have been proposed in [2,13,16]. For a detailed review of mathematical models of pedestrian dynamics we refer to the work by Cristiani and co-workers [17, Chapter 4].
Experiments confirm that collision avoidance is one of the main driving forces in pedestrian dynamics, see [25]. Individuals actively anticipate the behaviour of others and try to avoid collisions while maintaining their desired direction. The deviations from the desired direction result from stepping aside -either to the right or the left. In a force based model, such as the social force model proposed by Helbing and co-workers, sidestepping can be included by a force which depends on the estimated collision time and acts perpendicular to the vector connecting the two involved individuals. In cellular automata approaches sidestepping in bidirectional flows can be accounted for by enhancing the transversal transition rates (with respect to the desired direction). Based on these different microscopic collision avoidance mechanisms, various meso-and macroscopic models have been proposed and studied in the literature. At the mesoscopic level these interactions often result in complex collisional operators, see for example the work of Degond and co-workers [19]. Burger et al. [8] studied a minimal macroscopic model for bidirectional flows with sidestepping, which resulted in the segregation of the two flows into separate lanes.
In the present work we propose instead a minimal kinetic description, which is based on the assumptions that individuals try to avoid collisions while moving in their desired direction, see also [20]. We would like to mention that the collision avoidance plays an important role also in kinetic models of vehicular traffic. For detailed information we refer to [26,31]. This paper is organised as follows: in Section 2 we discuss microscopic collision mechanisms. In Section 3 we formulate a Boltzmann-type kinetic model for collision avoidance based on the minimal concepts of desired direction and sidestepping. In Section 4 we identify conditions for directional alignment in the case of spatially homogeneous pedestrians for both the Boltzmann-type model and its hydrodynamic counterpart under the grazing collision limit. We also show various numerical experiments which confirm our theoretical findings. Finally, in Section 5 we illustrate the dynamics of the proposed model for more general problems under spatially inhomogeneous interaction rules.

Microscopic collision mechanisms
We start by discussing various mechanisms of pedestrian collision at the microscopic scale. We consider a crowd of N > 1 pedestrians, whose position and velocity are denoted by x i ∈ R 2 and v i ∈ R 2 , i = 1, . . . , N , respectively. The dynamics of the pedestrians are driven by two goals: moving with a desired velocity or in a desired direction and, at the same time, trying to avoid collisions with others.
We assume that the reaction of a pedestrian to nearby walkers depends on his/her prediction of the time before a possible collision [25]. This time, the so-called time to collision, can be estimated by extrapolating the pedestrian positions on the basis of the current velocities. Specifically, we imagine that the individuals continue to walk with their current velocity, hence along straight paths, until they collide. Furthermore, since we regard pedestrians as point particles, we assume that a collision between two of them, say i and j, occurs once they get closer than a certain threshold γ > 0 and we fix the time to collision t = t ij at the instant in which their distance is exactly γ, cf. Figure 1. Therefore we obtain t ij by solving the equation with respect to t with the constraints t ∈ R, t ≥ 0, which yields where · denotes the inner product in R 2 while the discriminant ∆ is Notice that in the case |x i − x j | > γ a future collision can occur, i.e. (1) admits non-negative real solutions, only if the states ( While the first condition guarantees ∆ ≥ 0, hence that the solutions to (1) are real, the second condition ensures that the solutions to (1) are non-negative. In this case, the smaller one is retained as time to collision. In order to reduce the analytical and computational complexity of the model (considering that for two-dimensional position and velocity the dimension of the state space is 4), it may be convenient to assume that all the pedestrians walk at a constant speed: In particular, working with dimensionless variables, we set V 0 = 1 with reference to the standard walking speed of a pedestrian in normal conditions, i.e. O(1 m/s). This assumption implies: where i, j are the unit vectors of the coordinate axes in R 2 and θ i ∈ R is the angle giving the orientation of the velocity v i . Hence the microscopic state of the ith pedestrian is fully characterised by the pair (x i , θ i ) ∈ R 3 , see Figure 2, which reduces the dimension of the state space by one, cf. [1].
Remark 2.1. Because of the 2π-periodicity of the orientation, the angle θ i can actually be taken in any bounded interval I ⊂ R of length 2π. Throughout this paper we set The interval I is centred around a preferential angle α d ∈ [−π, π), which corresponds to the desired direction of the pedestrians, cf. again Figure 2.
where the discriminant is now given by In view of a statistical description of the ensemble of pedestrians, the time to collision can be used to define a probability of collision P between two individuals with states (x i , θ i ) and (x j , θ j ). For instance: where τ > 0 is a constant. Notice that 0 < P ≤ 1 for t ij ≥ 0 as expected and furthermore P → 1 when t ij → 0 + while P → 0 when t ij → +∞. Often, however, the qualitative asymptotic trends of the collision dynamics are studied in spatially homogeneous conditions, i.e. when pedestrians are so well mixed that their statistical distribution can be considered independent of the space variable. In this case, the function (5) has to be approximated by a function of the angles θ i , θ j only. In this paper we consider: i.e. we basically assume that P is proportional to the distance between the angles, which is linked to the incidence of the walking directions, taking into account the 2π-periodicity of the orientation of the velocities in R 2 .
In Figure 3 we compare the functions (5), (6) for two pedestrians located in x i = (0, 0) and in x j = (0, 1) with θ i = θ variable in the range [−π, π) and either θ j = π or θ j = − π 4 . In both cases, (6) can be seen as a piecewise linear approximation of (5) in the range of θ where conditions (2) are met. If instead θ is such that either t ij = 0 or t ij = +∞ then the spatially inhomogeneous and homogeneous probabilities of collision may differ consistently, because the spatially homogeneous model does not consider the relative position of the interacting pedestrians.

Kinetic modelling
Based on the microscopic modelling discussed before we now frame some principles of collision avoidance in a kinetic context. This approach allows us to make the modelling more amenable to mathematical analysis and to the understanding of group-wise dynamics.
As already discussed, we assume that the pedestrians walk at a constant speed, hence we describe their microscopic state by means of their position in space x ∈ R 2 and the orientation θ ∈ I of their velocity. We denote by f = f (t, x, θ) : R + × R 2 × I → R + the kinetic distribution function, which is such that f (t, x, θ) dx dθ is the fraction of pedestrians who at time t are in the infinitesimal volume dx centred in x with an orientation in [θ, θ + dθ]. We assume that f obeys the following Boltzmann-type equation: where the collisional term at the right-hand side describes binary interactions between pairs of pedestrians. Explicitly: with (x, θ * ), (y, ϕ * ) and (x, θ), (y, ϕ) representing the pre-and post-interaction states, respectively, and |J| the determinant of the Jacobian J of the transformation from the former to the latter. Notice that binary interactions are assumed to modify only the orientation of the velocity but not the position of the interacting individuals, which instead changes because of the transport term at the left-hand side of (7). In order to avoid the computation of the Jacobian J in Q it is customary to use the weak formulation of (7), which is formally obtained by multiplying the equation by a sufficiently smooth and x-compactly supported test function ψ : R 2 × I → R and integrating with over x, and θ: This form is easier to handle and can provide information about the evolution of some macroscopic quantities of the system. To characterise the collisional term we propose an interaction rule of the form where h, k ∈ Z are chosen to ensure that θ, ϕ ∈ I. Therefore they depend in general on θ * , ϕ * ∈ I. In (10), α d is the desired angle introduced in Remark 2.1 whereas the deviation angle α c ∈ [−π, π) corresponds to the lateral displacement of the pedestrians, when they step aside to avoid collisions. A pedestrian steps to the left when α c > 0 and to the right when α c < 0 (with respect to their current direction of motion). Finally, the term P (x, y, θ * , ϕ * ) is the probability of collision introduced in Section 2.

The spatially homogeneous problem
In the space homogeneous case pedestrians are assumed to be well mixed so that their distribution function f is constant in x, i.e. f = f (t, θ). In other words, the statistical distribution of the angle θ is the same at every point x.
The weak formulation (9) of the kinetic equation then becomes: We point out that, in contrast to (9), here the test function ψ : I → R does not have to be either smooth or compactly supported. Furthermore, in the following we assume that the test functions are extended to R by 2π-periodicity, so as to avoid unneccessary technicalities caused the periodic boundary conditions. In particular, taking the constant test function ψ ≡ 1 we see that the macroscopic density is constant in time because the right-hand side of (11) vanishes. This means that the crowd density can be regarded as a parameter, say ρ, of the model fixed by the initial condition. By dropping the transport term v · ∇ x f in the equation, the assumption of space homogeneity allows one to study the genuine effect of the microscopic interactions, which are now expressed as and P given by (6). In more detail, we assume corresponds to a function which belongs to (6), while a(ρ) ∈ [0, 1] models the influence of the collective state of the crowd on the individual interaction rules. Considering a dimensionless ρ ∈ [0, 1] referred to a typical congestion density O(4 ped/m 2 ), cf. [32,38], we suggest: (i) a(ρ) ∝ ρ to model a probability of collision proportional to the number of pedestrians in the crowd; (ii) a(ρ) ∝ ρ(1 − ρ) to account for a small probability of collision in the case of spread out groups (ρ ≈ 0) and in very crowded situations (ρ ≈ 1), when pedestrians tend to be passively dragged by the flow. The highest probability of collision arises at intermediate congestion levels (ρ ≈ 1 2 ). 4.1. Asymptotic alignment under binary interactions. Now we investigate the conditions under which interactions lead pedestrians to move in the same direction. We call this alignment, which is a form of consensus [4,10,27]. We say that there is alignment when the kinetic distribution f converges asymptotically in time to a distribution of the form ρδ α for some α ∈ I, where δ α is the Dirac delta centred at θ = α. In fact this means that, in the long run, the mass concentrates at θ = α.
We shall use the following notations and definitions: • Let M ρ + (I) denote the set of the positive measures on I with total mass ρ ≥ 0; • Let W 1 (µ, ν) denote the 1-Wasserstein distance [5] for µ, ν ∈ M ρ + (I), given by where Π(µ, ν) is the set of the transference plans from µ to ν, i.e. the measures on I 2 with marginals µ and ν. First of all, we show that α d is the only possible direction of alignment. Proposition 4.1. A distribution of the form ρδ α ∈ M ρ + (I), with ρ > 0 and α ∈ I, is a steady solution to (11) with interaction rule (12) and interaction probability (13) if and only if α = α d .
Next we study under which conditions the system aligns in the desired direction α d starting with a possibly generic initial distribution f 0 ∈ M ρ + (I). Theorem 4.2. Let ρ > 0 and P be given by (13). Define moreover Proof. We define for ρ > 0 and for every t > 0 the average distance of the crowd from the desired direction Proving that lim t→+∞θ (t) = 0 implies the thesis, because f (t, θ) ⊗ ρδ α d (ϕ) is a possible transference plan in Π(f (t), ρδ α d ) and therefore In particular, from the interaction rules (10) we have In fact the value of h ∈ Z is chosen such that C (θ * , ϕ * ) + 2hπ ∈ I, whereas in general one may have C (θ * , ϕ * ) ∈ I. In view of this and using the interaction probability (13) together with the triangular inequality we obtain Plugging this into (16) and using the definition ofθ we get In order to integrate the previous differential inequality we define . Going back toθ we deduce that From (14) we have which, in view of (18), implies thatθ(t) → 0 for t → +∞. In order to complete the proof we examine the cases a(ρ) = 0,θ 0 = 0. In the former we have µ = 0 and λ = −1, thus from (18) we deduceθ(t) ≤θ 0 e −ρt → 0 for t → +∞. In the latter we deduce insteadθ(t) = 0, i.e. f (t, θ) = ρδ α d (θ), for all t ≥ 0, which concludes the proof.  (14) are met we infer the following estimate forθ: where λ, µ are given in (17). In Section 4.2 we compare this estimate with the numerical values ofθ(t) computed by simulating the kinetic model with a Monte Carlo algorithm.   Figure 4 illustrates the regions of alignment in this particular case.
In Figure 5, we fix ρ = 1 2 and α c = π 5 , which comply with (21), and we represent the mean velocity of the pedestrians at two successive times, namely t = 10 −2 after one iteration of the  algorithm and t = 10 after 10 3 iterations. Note that Figure 5 shows only the subdomain [0, 1] × [0, 1], to illustrate the alignment at α d = 0 better.
Next we analyse the rate of convergence in time to the desired direction by plotting an MC approximation, sayθ MC , ofθ cf. (15). Following [29, Chapter 3] we takē where N is the number of particles used in the MC simulations and the Θ i 's are N values of the angle θ sampled from the (MC-approximated) probability distribution 1 ρ f (t, θ). In Figure 6 we compare the mapping t →θ MC (t) with the theoretical estimate (19) for the choices (ρ, α c ) = 1 2 , π 5 on the left and (ρ, α c ) = 1 3 , 3 5 π on the right (notice that also the latter pair satisfies (21)). In both examples the simulations converge faster than the theoretical rate.
In case (i), (21) predicts alignment for |α c | ≤ π 4 but the numerical simulations show alignment also for α c = 3 10 π, 2 5 π. Note that Theorem 4.2 provides only sufficient conditions to characterise the asymptotic trend of the system. On the other hand, the numerical simulations confirm that alignment is not possible for every value of α c . Specifically, we observe thatθ MC always converges to some possibly non-zero asymptotic value, nonetheless, due to Proposition 4.1, this does not indicate a possible alignment in a direction different from α d . It means instead that interactions may, in some cases, preserve or even increase the initial "disorder" of the pedestrian orientations. In fact we see that for α c ≥ 7 10 π the value ofθ MC either remains almost constant or increases in time.
In case (ii), (21) predicts alignment for ρ ≤ 10 19 ≈ 0.52 but the numerical simulations show alignment also for ρ < 9 10 . In case (iii) the conditions (14) of Theorem 4.2 are always violated, because for ρ = 1 it results a(ρ) = 1 and subsequentlyθ 0 < 0. In spite of this, the numerical simulations show that in some cases the alignment is possible, particularly whenθ is sufficiently close to α d = 0. This can be understood by noticing that forθ = α d = 0 the mean of (22) is preciselyθ, for f 0 is symmetric in I. Therefore on average the pedestrians are initially not too far from their desired direction. However, such an asymptotic alignment is soon lost by varyingθ, henceθ 0 .
The tests (i)-(iii) confirm that, although the hypotheses of Theorem 4.2 are possibly overlyrestrictive, the theoretical predictions capture qualitatively the role played by the parameters of the model in the alignment process.
4.3. The quasi-invariant direction limit. In this section we study the quasi-invariant direction limit, namely a form of grazing collision limit [37], see also [29,34], in order to extract the principal part of the interaction rules encoded in the Boltzmann collisional operator. This procedure allows us to derive a hydrodynamic-type model closer to the macroscopic scale, in which binary interactions are replaced by mean-field interactions defining a transport velocity in the space of the microscopic states.
The basic idea is to assume that pedestrian interactions get simultaneously weaker and more frequent. To this purpose we begin by rewriting (12) as where > 0 is a dimensionless parameter measuring the strength of the interaction. For = 1 we recover precisely (12). At the same time, we scale the time in the Boltzmann equation (11) as t/ , meaning that the interaction rate becomes O(1/ ), so that we finally have where now we take ψ ∈ C 2 (I) with ψ(±π + α d ) = 0, still 2π-periodic on the whole real line.
As we are interested in small values of , in view of (23) and of the said periodicity of ψ we can expand for someθ ∈ I depending on θ * and θ. Let us define for brevity; plugging into (24) we obtain where, considering that |P(θ * , ϕ * )| ≤ π because both angles α d − θ * and α c belong to [−π, π) while 0 ≤ P (θ * , ϕ * ) ≤ 1, the remainder R( ) satisfies with · ∞ denoting the usual ∞-norm. Therefore in the limit → 0 + we find where, recalling also (13), we define Equation (25) corresponds to the weak form of the PDE which is a non-linear conservation law with non-local flux.

4.4.
Asymptotic alignment under mean-field interactions. Like for the Boltzmann-type model, we are interested in the steady state solutions of (27) to find out under which conditions this equation predicts the emergence of group-wise alignment. Since (27) has been derived as an approximation of (11) on a time scale much larger than that of the binary interactions (recall the scaling t → t/ ), it is reasonable to expect that the large-time behaviour of its solutions, including the conditions for alignment, approximates the asymptotic trends of (11). To begin with, by a straightforward calculation we show, like in Proposition 4.1, that also for (27) the only possible alignment corresponds to α d . Since G(0) = 0 we further have H[ρδ α ](α) = ρ(α d − α), cf. (26). This equation has to hold for every ψ, therefore α = α d .
Next we claim that, at least in a suitable density-dependent range of values of the deviation angle α c , the group-wise alignment at θ = α d is the configuration reached asymptotically in time by the mean-field interaction model regardless of its initial configuration. In other words, we can give (sufficient) conditions on ρ, α c ensuring that ρδ α d is a globally attractive equilibrium to (27).
To tackle this we introduce the characteristic line of (27) issuing from θ ∈ I as the mapping t → γ t (θ) satisfying where the dot over a variable means henceforth time derivative. In practice, γ t (θ) is the direction at time t > 0 of a pedestrian who initially (t = 0) moved in the direction θ. We point out that γ t (θ) is understood "modulus 2π", meaning that, as a function of θ, it maps I into itself at every time. Since (27) describes a transport at speed H[f ] of the initial kinetic distribution f 0 ∈ M ρ + (I) in the state space I, the (weak) solution f can be written formally by pushing f 0 forwards in time along the characteristics: (29) f (t) = γ t #f 0 ∈ M ρ + (I), t > 0. We study the characteristics of (27) by the following lemma, which establishes the continuous dependence of the mapping t → γ t (θ) on the angle θ ∈ I. Lemma 4.5. Set C := 1 + |αc| π with α c ∈ [−π, π). For all θ 1 , θ 2 ∈ I it holds: Proof. First we notice that, using (29), we can rewrite the transport speed H[f ] as We fix now θ 1 , θ 2 ∈ I. Writing (28) first for θ = θ 1 then for θ = θ 2 and taking the difference of the two equations we get d dt Let us temporarily call I the second term at the right-hand side. Since γ t (θ 1 ) ∈ I for all t > 0 by 2π-periodicity, it results γ t (θ 1 ) − α d + α c ∈ [−π + α c , π + α c ), hence where C is the constant defined in the statement of the lemma. Moreover, from (13) we see that the function G is Lipschitz continuous with Lipschitz constant equal to 1 π , therefore on the whole we have

This implies on one hand
and on the other hand d dt thus, setting for ease of notation, Let Λ(t) := t 0 λ(s) ds. Multiplying the chain of inequalities above by e −Λ(t) and integrating in time gives Invoking now Gronwall's inequality we deduce that which, after substituting the definitions of u(t), u 0 in terms of γ t , γ 0 , concludes the proof.
Next we establish a sufficient condition on the deviation angle α c such that all the characteristics of (27) converge in time to the desired angle α d . Proposition 4.6. Let α c ∈ [−π, π) be such that Then lim Proof. We preliminarily observe that, under the stated constraint on α c , the constant C introduced in Lemma 4.5 is such that a(ρ)(1 + C) < 1. Moreover, for θ, ϕ ∈ I it results |θ − ϕ| < 2π, thus from Lemma 4.5 we get on the whole uniformly with respect to θ, ϕ. We can rephrase this by saying that for all > 0 there exists a time t 0 > 0, independent of θ, ϕ, such that |γ t (θ) − γ t (ϕ)| < for all t > t 0 . Since G is continuous with G(0) = 0, we conclude that there also exists t 0 > 0 such that G(|γ t (θ) − γ t (ϕ)|) < for all t > t 0 and all θ, ϕ ∈ I.
Let us now consider the equation of the characteristics: For t >t := max{t 0 , t 0 } we obtain the following estimate: We set u(t) := γ t (θ) − α d and observe thatu(t) =γ t (θ), multiply the chain of inequalities above by e ρt and integrate over the time interval [t, t] to obtain In the limit t → +∞ we find which, since > 0 was arbitrary, implies the thesis.
Then lim for all f 0 ∈ M ρ + (I). Proof. Proceeding like in Theorem 4.2 we have hence, by dominated convergence and Proposition 4.6, which concludes the argument.
Remark 4.8. Conditions (14) and (31) resemble each other, thereby confirming the expectation that the asymptotic behaviour of the mean-field interaction equation should approximate that of the binary interaction model. However, a remarkable difference, which makes it not immediate to compare the two conditions in general, is that the former depends on the initial condition f 0 through the quantityθ 0 while the latter is independent of f 0 .
We can get a clue on the relationship between the two conditions by fixing the special case f 0 (θ) = ρ 2π , θ ∈ I, which corresponds to the initially most "disordered" situation with pedestrians following homogeneously all possible directions. Conditions (14) specialise then in (21) and can be compared with (31) for different choices of the function a(ρ).
In Figure 8 we consider the cases a(ρ) = ρ and a(ρ) = 3 2 ρ(1 − ρ) already discussed at the very beginning of Section 4. In the former we see that either condition may be more or less restrictive depending on the considered density range. In the latter we observe instead that the mean-field interaction model may guarantee alignment in a range of values of α c invariably larger than the corresponding one of the binary interaction model. Finally, in both cases we notice that the density range for unconditional alignment, i.e. alignment for every α c ∈ [−π, π), is typically larger in the case of the mean-field interaction model.
We should stress that the considerations proposed in Remark 4.8 apply a priori to the two types of model, but alignment in situations not encompassed by either (14) or (31) can still be observed a posteriori. In fact these conditions are on one hand quite general, in that they are not too bound to the specific form of the initial distribution f 0 , but on the other hand, as already observed from the numerical simulations of the binary interaction model, cf. Section 4.2, only sufficient. In the case of the mean-field interaction model, we prove that particular initial distributions might produce alignment also in cases not explicitly covered by Theorem 4.7. Proof. The equation of the characteristics for the given f 0 iṡ In particular, for θ = α we obtainγ which, together with γ 0 (α) = α, gives γ t (α) = α d + (α − α d )e −ρt . Using this in the calculations of the proof of Theorem 4.7 yields finally whence the thesis follows.
We conclude that the extra microscopic information brought by the kinetic distribution function f 0 may be essential to get a full understanding of the alignment dynamics in particular cases. 4.5. Numerical simulations. Next we illustrate the time-asymptotic evolution of the solution to (27) by various numerical experiments. To discretise the equation we use the semi-Lagrangian scheme described in Appendix A.2, which has good stability and accuracy properties even in the case of non-local fluxes.
We consider the uniform initial distribution (20) with density ρ = 1 3 and we take α d = 0, a(ρ) = ρ, α c = π 3 . In this case, consensus towards α d is asymptotically expected because, as shown in the left panel of Figure 8, these parameters fall in the range of unconditional alignment. The corresponding time evolution of the kinetic distribution function f is displayed in the top row of Figure 9.
Changing the initial condition, for instance using the function (22) with various choices of the parametersθ, σ 2 , does not affect the asymptotic behaviour of the system, cf. the bottom row of Figure 9. This is indeed consistent with the theoretical prediction of Theorem 4.7.
Similarly to the numerical tests presented in Section 4.2, we also perform a parametric study of the trends of the mean-field interaction model for large times by computing a Finite-Element-type approximation, sayθ FEM , of the quantityθ(t). Specifically, we compute the integral in (15) by the rectangle method over a grid of points {θ i } M i=1 ⊂ I with step ∆θ := 2π M +1 , see Appendix A.2 for further details:θ In Figure 10 we illustrate the evolution of the mapping t →θ FEM (t) for: (i) ρ = 1 2 and setting α c = k 10 π for k = 1, . . . , 10, cf. Figure 10 (left); (ii) α c = π 5 and setting ρ = k 10 for k = 1, . . . , 10, cf. Figure 10 (right). Case (i) is completely not covered by condition (31); case (ii), instead, is covered for ρ < 5 11 ≈ 0.45. In both cases, however, we do not find any numerical evidence that the alignment fails for some choices of the parameters. As already discussed, this is not in contrast with condition (31), which indeed is only sufficient.

The spatially inhomogeneous problem
In the spatially inhomogeneous case one considers a generic distribution in space of the pedestrians, which can vary from point to point x ∈ R 2 and over time. Therefore the kinetic distribution function fully depends on both state variables (x, θ), which implies that the density of the crowd: is no longer a constant parameter of the model. However, from (7)- (8) we see that the total mass of the system, namely the quantity I R 2 f (t, x, θ) dx dθ = R 2 ρ(t, x) dx, is conserved in time because I R 2 Q(f, f )(t, x, θ) dx dθ = 0 for all t > 0. In the following we consider the case I R 2 f (t, x, θ) dx dθ = 1, hence f can be understood as a probability density.  In this section we explore computationally the space inhomogeneous problem by means of the binary collision model, cf. (7)- (8) or (9), which we simulate using a Nanbu-like Monte Carlo particle method, see Appendix A.1 for details. We show that the microscopic model of walking behaviour discussed in Sections 2, 3, though much simplified with respect to other sophisticated models such as e.g. [19], is nonetheless able to reproduce various emergent macroscopic patterns while being more realistic from the phenomenological point of view than the interaction models based on position-dependent repulsion potentials.
In the case studies addressed here we consider a two-dimensional bounded spatial domain, specifically the square Ω := [− L 2 , L 2 ] ⊂ R 2 with edge length L > 0 and periodic boundary conditions. Table 1 states all the relevant parameters used in the simulations.

Test 1 -Alignment.
In the first test we consider a similar situation as for the space homogeneous model in Section 4. We choose an initial condition of the form such that I Ω f 0 (x, θ) dx dθ = 1, which corresponds to a crowd concentrated mainly in a horizontal middle stripe of Ω with walking direction uniformly distributed in I, see Figure 11a. Like in the space homogeneous case, the pedestrians tend to align group-wise to the desired direction in finite time (cf. Figure 11c). However, the simulation shows that such a consensus is faster where the density is higher (cf. Figure 11b). This can indeed be inferred also from the results of the space homogeneous case, cf. the central panel of Figure 7. Moreover, we notice that during the transition from the initial condition to the consensus the individuals tend to move from high to low density areas due to sidestepping for collision avoidance (cf. the upper and lower white stripes in Figure 11b). Such an effect appears to be top/bottom-symmetric despite the fact that the deviation angle α c is constant across the domain (in this case α c > 0, cf. Table 1, corresponding to leftwards stepping).
This may explain the early non-monotone trend of the mapping t →θ MC (t), cf. Figure 11d, which in the present case is computed as a Monte Carlo approximation of In particular, the steep initial rise might be due to the fact that most individuals in the top area of the domain have to span anticlockwise all the walking directions in the interval I before finding one free from collisions with neighbouring pedestrians. The subsequent exponential-like decay of θ MC is instead consistent with the observations made in the spatially homogeneous case.

Test 2 -Counterflow.
A very popular problem in crowd dynamics is the counterflow, i.e. the case of two groups walking towards each other. In order to simulate this situation we consider an extension of model (7) to two groups of pedestrians, each described by its own distribution function f p = f p (t, x, θ) which satisfies the equation The additional collisional term Q(f p , f q ) at the right-hand side takes into account the interactions of either group with the opposite one. We assume that the two groups differ only in the desired directions, which we take to be oriented rightwards and leftwards, respectively. Hence each collisional term in (32) implements an interaction rule of the form (10) but with α d replaced by α p d , cf. Table 1. We consider for both groups an identical initial distribution which is uniform in θ ∈ I and in x ∈ − L 2 , L 2 × −η L 2 , η L 2 ⊂ Ω, 0 < η ≤ 1, and is such that I Ω f p 0 (x, θ) dx dθ = 1. Hence the two crowds share initially the same area of the domain and are well mixed therein (Figure 12a).
The simulation shows that the two groups fully segregate, see Figure 12. A similar behaviour has been ovserved in the literature for different types of models, cf. e.g. [8,15,17,28]. Furthermore, within either lane each group aligns to its desired direction. As a matter of fact, this is the optimal way for pedestrians to avoid collisions with others.
In the early stages of the segregation (Figure 12b) the lanes are thinner and well separated by an area in which pedestrians with opposite desired directions are still mixed. In the long run, when either group aligns to its desired direction (Figure 12c), the lanes become wider and the groups clearly segregate. Figure 12. Test 2 -Lane formation in pedestrian counterflow. For pictorial purposes, in order to represent the two crowds in the same picture, the blue density is plotted on a negative scale. In white areas where the velocity is nonzero the two densities take nearly the same values.

Test 3 -Crossing flows.
We now consider two groups of pedestrians having perpendicular desired directions. As initial conditions we prescribe the distribution functions which are analogous to (33) but for the fact that here the smaller edge of the initial stripe is in the x q -direction. Hence the group p = 1 is distributed horizontally and walks rightwards while the group p = 2 is distributed vertically and walks upwards, see Figure 13a and cf. Table 1.
The snapshots in the top row of Figure 13 show that the pedestrians anticipate the interactions by starting to step aside slightly before they reach the area in which the two streams cross. As a result, the macroscopic flows deviate from their initial orthogonal directions thereby giving rise to quite realistic patterns.
In order to better investigate the dynamics in the crossing area, we consider at last two groups of pedestrians, with the same perpendicular desired directions as before, that are now uniformly distributed in the whole spatial domain Ω. As Figure 13d shows, when they begin to interact we observe the emergence of segregation with stripe formation, a phenomenon well-documented in the experimental literature [21,33,39]. The slope of the stripes corresponds to an angle between α 1 d and α 2 d , which in the present case seems to be close to However, further analytical study of the model is needed to possibly confirm this guess.

Appendix A. Numerical tools
In this appendix we provide more details about the numerical schemes that we used for simulating the Boltzmann-type collisional model (both homogeneous and inhomogeneous in space) and the hydrodynamic mean-field model deduced from the former in the quasi-invariant direction limit.
A.1. Nanbu-like Monte Carlo algorithm. The Nanbu algorithm is a particle method belonging to the family of the Monte Carlo numerical methods for the approximate solution of collisional kinetic equations. We used it for producing the simulations presented in Sections 4.2 and 5.
Here we follow the approach presented in [3,29], see also [6], to approximate the spatially inhomogeneous Boltzmann-type equation (7) with microscopic states (x, θ) ∈ R 3 and v = V 0 (cos θi+ sin θj). The reader can easily adapt this description to the spatially homogeneous case, when the transport term v · ∇ x f drops. The algorithm is based on a splitting procedure: a first collisional step is followed by a transport step according to ∂ t f (t, x, θ) = Q(f, f )(t, x, θ) (34a) ∂ t f (t, x, θ) + v · ∇ x f (t, x, θ) = 0. (34b) If we assume that, up to normalisation, the initial datum f 0 = f 0 (x, θ) is a probability density, i.e. I R 2 f 0 (x, θ) dx dθ = 1, then, since the total mass is conserved, f (t, x, θ) is in turn a probability density in (x, θ) for all t > 0. Hence we can rewrite the operator Q involved in the collisional step as, cf. (8), where Q + (f, f )(t, x, θ) := x, θ * )f (t, y, ϕ * ) dy dϕ is the so-called gain operator, which implements the interaction rule. Still because of mass conservation, or equivalently I R 2 Q(f, f )(t, x, θ) dx dθ = 0, also the gain operator Q + (f, f )(t, x, θ) is a probability density in (x, θ) for all t > 0. Equation (34a) is discretised in time by the forwards Euler scheme to get (35) f n+1 (x, θ) = (1 − ∆t)f n (x, θ) + ∆tQ + (f n , f n )(x, θ), where ∆t > 0 is the time step and f n is an approximation of f at time t n := n∆t, n = 0, 1, 2, . . . .

Algorithm 1
Nanbu-like algorithm for (7) 1: fix N > 1, ∆t ∈ (0, 1] 2: define the random variable ξ ∼ Bernoulli(∆t) 3: for n = 0, 1, 2, . . . do 4: sample N particles with states {(X n i , Θ n i )} N i=1 from the distribution f n 5: for i = 1 to N do Under the constraint ∆t ≤ 1, from (35) we see that f n+1 is a convex linear combination of f n and Q + (f n , f n ), hence it is in turn a probability density. This is the key to give (35) the following probabilistic interpretation in terms of the underlying microscopic particle system: to obtain samples distributed according to f n+1 we have to sample either from f n , with probability 1 − ∆t, or from Q + (f n , f n ), with probability ∆t. As a matter of fact, in (35) f n is related to the event that no binary interactions take place during the time ∆t, so that the kinetic distribution does not change in the time step n → n + 1, while Q + (f n , f n ) is related to the event that an interaction occurs. Therefore, after sampling from f n a certain number N of particles with states {(X n i , Θ n i )} N i=1 , with probability 1−∆t we set Θ n+1 i = Θ n i , i.e. we leave the direction of the ith particle unchanged. Conversely, with probability ∆t we: (i) select uniformly a second particle (X n j , Θ n j ), j ∈ {1, . . . , N }, with j = i; (ii) update Θ n i to Θ n+1 i according to the binary interaction rule (10) applied to the ith and jth particles. In the subsequent transport step we update the positions of the sampled particles. Owing to (34b), we have: X n+1 i = X n i + V 0 (cos Θ n i i + sin Θ n i j) ∆t, i = 1, . . . , N.
The whole method is summarised in Algorithm 1 in the form of a pseudo-code.
A.2. Semi-Lagrangian scheme for conservation laws with non-local flux. In recent works, see e.g. [18,24], semi-Lagrangian (SL) schemes have been proved efficient in approximating conservation laws. Since these numerical methods are not classical and need some adaptations to our case, we report here a short description of their derivation for the reader's convenience. In particular, following [12,11], see also [30,35], we describe a SL scheme for the numerical approximation of the following problem, cf. Section 4.5: To derive the scheme, we begin by fixing a time step ∆t > 0 through which we define the discrete time instants t n := n∆t, n = 0, 1, 2, . . . . Then we multiply the first equation in (36) by a sufficiently smooth test function ψ : I → R with ψ(±π + α d ) = 0 and integrate over I × [t n , t n+1 ] using integration by parts with respect to θ: