Self-propelled hard disks: implicit alignment and transition to collective motion

We show that low density homogeneous phases of self propelled hard disks exhibit a transition from isotropic to polar collective motion, albeit of a qualitatively distinct class from the Vicsek one. In the absence of noise, an abrupt discontinuous transition takes place between the isotropic phase and a fully polar absorbing state. Increasing the noise, the transition becomes continuous at a tri-critical point. We explain all our numerical findings in the framework of Boltzmann theory, on the basis of the binary scattering properties. We show that the qualitative differences observed between the present and the Vicsek model at the level of their phase behavior, take their origin in the complete opposite physics taking place during scattering events. We argue that such differences will generically hold for systems of self-propelled particles with repulsive short range interactions.

Self-propelled particles borrow energy from their environment and convert it to translational motion. Depending on their density and interaction details, energy taken up on the microscopic scale can be converted into macroscopic, collective motion. A transition from an isotropic phase to a polar, collective phase is observed. Most theoretical knowledge about such active matter and its transitions has been developed for Vicsek-like models: point particles travel at constant speed, to represent self-propulsion, and their interaction rules comprise both explicit alignment and noise, to account for external or internal perturbations. In addition, self-diffusion noise may be incorporated in the individual motion. Typically, the alignment strength and noise intensity control the transition between the isotropic and the polar phases [1][2][3][4][5][6].
In order to describe more closely experiments with bacteria [7][8][9] or with self-propelled actin filaments [10], other models incorporates the alignment mechanism through the steric interaction of elongated objects [11][12][13]. The interaction accounts for the body-fixed axis of these particles, which is in some sense more microscopic than the purely effective (mesoscopic) alignment rules in Vicsek-like models. Steric interaction is not the only microscopic interaction that can lead to collective behaviour. An alternative way to use dynamics for aligning (hard or soft) circular disks is to incorporate inelasticity or softness in the interaction rules [14][15][16][17]. Not surprisingly, also these models exhibit a transition to a globally aligned phase. A model experimental system of self-propelled disks has been realized by vertically vibrating millimeter-sized circular metal disks with two different "feet" [18,19]. When isolated, such particles advance in the direction of their body axis, with a reasonably high persistence length. When sufficiently crowded, they start to align to each other and move along in dense clusters. Successful modelling of this experiment, in the form of Newton's equations for rigid disks with a force acting on the center and a torque turning it around, plus slightly inelastic collisions and some active noise included, have demonstrated that purely repulsive mechanical interaction with no explicit alignment, were sufficient ingredients to observe a transition to collective motion [20]. The nature of this effective alignment remains however mysterious. Is inelasticity a required ingredient? Or, as suggested in [18][19][20], does it come from repeated recollisions of particle pairs, during which the velocity vector repeatedly converges against the body axis, but also the body axis has time to relax towards the timeaverage of the velocity axis. Given that the experiment is very much the one closest to a simple hard-disk liquid and lends itself to a description in terms of statistical physics, understanding precisely its phase transitions, the effective alignment mechanism, and the role of the noise would be a major step towards a theoretical description of active liquids in general.
In this letter we provide such an understanding, restricting ourselves to the analysis of spatially homogeneous  Figure 1: (a) An isolated self-propelled particle converges to its stationary state where velocity v and polarityn are parallel. (b) A single binary "scattering event" can consist of many harddisk elastic collisions. (c) Stable phases in the absence of noise. Between isotropic ( ψ ≈ 0) and polar ( ψ = 1) phase is a discontinuous transition. phases in the low density limit. To do so we proceed in three steps. (i) We perform molecular dynamics of the model equations for purely elastic interactions, with and without noise: In the absence of noise, the system exhibits a strongly first order transition from the isotropic to the collective motion phase (see Fig. 1c). Above a finite level of noise, the transition becomes second order -a tricritical point exists. This establishes the phase behaviour which we will explain from theoretical considerations. (ii) We analyze the model equations on the grounds of the Boltzmann equation, by making use of a recently proposed observable p · δp which quantifies the non-conservation of momentum [21]. This observable allows to span the bridge from the microscopic dynamics, in particular binary collisions such as depicted in Fig. 1b, to macroscopic order parameters. From a direct sampling of all possible binary scattering events, we obtain an excellent quantitative prediction of our numerical findings. (iii) We scrutinize the very peculiar dynamics of a collision between two self propelled disks and explain the specific shape of the scattering function that was obtained numerically in (ii). We further find that recollisions are not necessary for the observed alignment, contrary to our previous belief.
Model of self-propelled hard disks. The model consists of N hard disks in a square box of size L×L, with periodic boundary conditions. The density is ρ = N/L 2 . Particles, being self-propelled, relax to a stationary speed v 0 . As units of length and time we choose the diameter d 0 of the particles and d 0 /v 0 , respectively. A particle i has coordinates r i , velocity v i , and a body axis given by the unit vectorn i (see Fig. 1a). Between collisions, it evolves according to the equations The competition between the self-propulsionn and the viscous damping −v in Eq. (1b) lets the velocity relax ton on a timescale τ v . Similarly, in Eq. (1c), the polarityn undergoes an overdamped torque that orients it toward v on a timescale τ n . Interactions between particles are elastic hard-disk collisions which change v but notn. After such a collision, v andn are not collinear, and the particles undergo curved trajectories which are either interrupted by another collision (Fig. 1b), or the particles reach their stationary state, where v =n and the trajectory is straight at a speed v 0 = 1 (Fig. 1a). The final direction of v (equal to that ofn) depends on the parameter which can be understood as the persistence of the polarityn. Linearizing the evolution equations around the stationary state, one can show that the final polar angle is given by the weighted average of the initial angles, (θ n + αθ v )/(1 + α). When α 1,n is practically always directed along v.
On top of the deterministic trajectories given by the Eqs. (1), we add some angular noise by the following procedure. Given a time step δt → 0, we rotate v i andn i by the same angle η i (t), distributed normally with zero mean and variance 2Dδt, where the constant D ≥ 0 fixes the level of the angular noise. Noises of different particles are statistically independent. We choose δt much smaller than all other timescales in the dynamics. The relevant parameter to characterize the angular noise is then D/λ, where λ = 4ρ/π is the characteristic scattering rate of the system, which is proportional to the density [21].

Molecular dynamics (MD) simulations.
We now establish the phase behaviour of the model for N particles. MD simulations were performed at τ v = 4 with N = 1000 or N = 4000, focusing on the dilute regime ρ 1 (see below for a discussion of the effect of τ v ). We are thus left with two microscopic parameters, namely α and D/λ. Also, the system size is chosen not too large, in order to keep the system spatially homogeneous, which we have checked by visual inspection. We measured the order parameter ψ(t) = i v i (t) /N , which is of order 1/ √ N for the isotropic state and of order one for the polar state.
Let us first look at the case without angular noise, D/λ = 0. We initialized simulations from random isotropic conditions and waited for the isotropic state to eventually destabilize. When a stationary state was reached, we started to average the order parameter over time, ψ . As shown in Fig. 1c, we found the isotropic state to be stable at low values of α, whereas it becomes unstable at larger values, in favour of a polar state. Between the two phases, an abrupt discontinuous transition takes place at α * . Quite remarkably, in the whole polar phase the dynamics converges to ψ = 1, where particles are all strictly parallel. Further, choosing some random state with ψ ≈ 1 as initial condition, we found that the polar state ψ = 1 is stable for all α > 0, in particular also when α < α * .
In Fig. 2a, we show again the (now rescaled) order parameter in the isotropic state, this time for different densities. For a given density, the data for different values of N collapse, showing that finite-size effects are under control. In all cases we observed convergence to the fully polar state beyond the points shown. The theoretical framework used below adds the line for ρ → 0, drawn in black. Increasing density here clearly favours the polar state. The departure of the transition due to density effects, 1 − α * (ρ)/α * (0), is plotted in Fig. 2b. Adding angular noise to the trajectories quite changes the picture. During simulations we first increased D/λ quasi-statically and then decreased it again. The transition was measured by looking at the maximum of the fluctuations of the order parameter among many realizations of the dynamics. The resulting phase diagram at density ρ = 10 −2 is shown in Fig. 2c. In agreement with intuition,the isotropic state is always stable at strong enough angular noise. When decreasing the noise we pass into the polar phase, but the nature of this transition can be either discontinuous or continuous. For values α < α c ≈ 0.157 the transition has some hysteresis, as indicated in the upper panel of Fig. 2d. The discontinuous nature of the transition is thus robust when adding angular noise, and the phase areas in Fig. 2c overlap, presenting an area which we could call coexistence region if the system were not homogeneous. For α > α c , the hysteresis is no longer observed at our level of numerical precision 1 , as can be seen in the lower panel of Fig. 2d. At α c , the coexistence zone vanishes into a single line of transition (tricritical point [22][23][24]). Interestingly, once in the polar phase, further increasing α leads to a re-entrant transition towards the isotropic phase.
Kinetic theory framework. We now rationalize these numerical observations in the context of kinetic theory, using the properties of the binary scattering, following Refs. [17,21]. We summarize here the procedure described in Ref. [21]. In the dilute regime, the mean free-flight time λ −1 is long enough so that particles have mostly reached their stationary velocity v 0 before interacting with another particle (τ v , τ n λ −1 ). The binary scattering of selfpropelled particles does not conserve momentum, and that is why a polar state can emerge from an isotropic initial condition. Assuming molecular chaos and choosing the von Mises distribution as an ansatz for the angular distribution of the velocities, one can write down an evolution equation for the order parameter ψ. This equation can then be expanded, up to order ψ 3 to study the stability of the isotropic phase [21]: where the coefficients are given by In the stationary state, the left-hand size of Eq. (3) vanishes. The transition line is obtained by solving the equation µ(α * ) = D/λ for α * , while the sign of ξ(α * ) at the transition tells whether it is continuous or discontinuous. The coefficient µ is exact within the assumptions of kinetic theory, while ξ should depend on the ansatz used for the angular distribution. Both coefficients are an average over all pre-scattering parameters, as given in Eq. (6), where b is the impact parameter and ∆ is the angle between the incoming particles' velocities. The averaging needs not be done over the norms of the velocities, since those are fixed to v 0 = 1. We have checked explicitly that this assumption holds very well in the numerical simulations in the dilute regime. In Eqs. (4) and (5), p is the pre-scattering momentum of the two colliding particles, and δp is the change of their momentum by the scattering event. In contrast with Ref. [17], the above equations explicitly identify the forward component of the momentum change p · δp, as the proper quantity to describe the alignment in a binary scattering event. The predictions thus depend only on the microscopic details of the deterministic dynamics through the scattering function p · δp(b, ∆), which we now explore for the model (1).
Binary scattering. Remind that a scattering event can consist of several, if not many, hard-disk collisions, such as depicted in Fig. 1b. Before the first collision, both particles are taken to be in the stationary state, i.e. v =n.
After the collision, we integrate Eqs. (1) numerically until another collision possibly occurs. The binary scattering is considered over when both particles have again reached their stationary state and are heading away from each other. Repeating the procedure in the whole range of initial parameters (b, ∆) yields the scattering function p · δp(b, ∆). Figure 3 shows this function for several α, as well as its integral over b and the full integrals yielding the coefficients µ and ξ.
Let us stress that p · δp is not changed by the collision itself, which conserves momentum. All (dis)alignment must here come from the relaxation of the post-collisional value of |v| to unity. Fig. 3a-d shows that scattering at low angular separation, small ∆, always creates forward momentum. In other words, two nearly parallel particles that interact become even more parallel, which gives rise to an effective alignment, p · δp > 0. On the other hand, for small enough α, particles that enter in interaction frontally (∆ ≈ π) tend to disalign, except for special symmetry such as b ≈ 0. Increasing α favours aligning scattering events until eventually only aligning events remain. This is best summarized by integrating out all parameter dependence except the incoming angle, as is plotted in Fig. 3e,f. The coefficients µ and ξ are then obtained by integrating over ∆, with weights prescribed by Eqs. (4) and (5). Their dependences on the microscopic parameters of the dynamics, α and τ v , are shown in Fig. 3g,h. In the absence of noise, the transition occurs for α = α * such that µ(α * ) = 0; ξ(α * ) is negative, hence the discontinuous transition. When angular noise is added, the transition is obtained by solving the equation µ(α * ) = D/λ. From the shape of the curve µ(α), one obtains two values α ± , with α − → α * and α + → ∞ when D → 0. As D/λ is increased, ξ(α − ) eventually becomes positive and the transition turns continuous at a tricritical point (α c , D c /λ). Note that ξ(α + ) > 0: the re-entrant transition is always continuous. Regarding the role of τ v , α − is practically independant on τ v , while α + increases when τ v decreases. The theoretical predictions are shown as solid lines in Fig. 2c. The agreement with the MD simulations data for density ρ = 10 −2 is excellent. The small shift of the measured transition lines to the left with respect to the theoretical one comes from finite-density effects.
Finally, we also learn from the examination of the scattering maps that, in the absence of noise, the polar phase ψ = 1 is actually an absorbing phase [25]: this is because all binary scattering events at small ∆ have p · δp > 0. When all particles in a system are sufficiently parallel, binary scattering events can only align the system more. This is true for all α, and most remarkably for α → 0.
Altogether, our kinetic theory description, using the von Mises ansatz for the angular distribution, captures quantitatively all the phenomenology reported in the numerical simulations at low enough density. It however relies on the numerical evaluation of the scattering maps. In the last part of the paper, we would like to provide some intuition on the origin of the peculiar form of these maps. Also, we will elucidate the role of the multiple collisions which can take place during a scattering event.
When we take the limit τ n → 0, the vectorn has no persistence at all. After two particles collide elastically, then rotate to their respective v instantaneously, the two particles cannot collide a second time and the post-collisional relaxation of |v| to unity occurs on a straight trajectory. We take advantage of this simplification to compute the scattering map in the case α = 0. Because velocities have equal modulus, we can use the nice visualisation of an elastic collision as a rotation in the center-of-mass frame (Fig. 4a,b). One can then find an analytic expression 2 for p · δp(b, ∆), which is plotted in Fig. 4c. This plot nicely completes the series of varying α from Fig. 3a-d. The essential structure of the scattering maps can be accounted for by the sole ingredients present in the case α = 0. Inelastic collisions and persistence ofn are not essential. Conversely, non-conservation of momentum due to the relaxation of the particle speeds to unity is crucial.
At values α > 0, the vectorn has non-vanishing persistence, which results in curved trajectories and possible recollisions. This could lead to the effective alignment mechanism which was proposed to be at the root of the collective motion in the experiment of vibrated polar disks [18][19][20]. Concerning the influence of recollisions in the scattering, we ask whether they contribute significantly to µ and ξ. In particular, how numerous are they depending on the scattering geometry (b, ∆), what is their statistical weight and how much do they affect the scattering map p · δp(b, ∆)? One can get an intuitive picture of the recollision mechanism looking at a simple case, the symmetric (b = 0) binary scatterings in Fig. 5a. A first collision makes v rotate from being equal ton to pointing away from the other particle. On the following trajectory, v andn relax towards each other, on their respective timescales τ v and τ n . For small enough τ v , the persistence ofn allows for a recollision to take place, and both vectors get closer to each other from one collision to the next. The highly symmetric case depicted in Fig. 5a is however so degenerate that it has either one or infinite collisions 3 . Unfortunately, the analytic treatment of asymmetric collisions is much too technical to gain physical insight from it and we prefer coming back to the numerical sampling of the scattering events.
For α = α − , as shown in Fig. 5b (top row), recollisions can be numerous, but take place only in a small region of the parameter space (b, ∆), the smaller the larger τ v . We measure their relative contribution to p · δp by Υ = 1 2 sin(∆/2)(p · δp − p · δp| 1 )/|p · δp|, where p · δp| 1 is the contribution of the first collision only. As shown in Fig. 5b (bottom row), even for the smallest value τ v = 0.04, this contribution is locally less than 5 %. For larger τ v , the contribution is vanishingly smaller. Thus, recollisions has no practical influence on the transition at α = α − . Conversely, as shown in Fig. 5c (top row), for α = α + , the number of recollisions is not so large, but they occur for wider ranges of scattering parameters and their relative contribution to p · δp is significant, see Fig. 5c (bottom row). This is all the more true for smaller τ v and is responsible for the dependance of µ and ξ on τ v in this regime.
Discussion. Let us recast our main findings and the understanding we obtained from the above analysis. A homogeneous system of self propelled hard disks, obeying the minimal deterministic Eqs. (1), exhibits a sharp discontinuous transition from an isotropic to an absorbing polar collective motion state, when increasing the persistence of the body axis vectorn. Adding angular noise, the transition becomes continuous via a tricritical point. At one collision to the next. For α < 1, there is no recollision. For α > 1, an infinite number of recollisions occur, with a collision rate that converges to a positive constant (α < 2) or diverges (α > 2). fixed noise level, increasing further the persistence of the body axis vector, a re-entrant continuous transition from the polar state to the isotropic state occurs.
All these macroscopic behaviours can be understood by the study of the scattering maps p · δp, which describe the level of alignment of all possible scattering events. The most striking feature of these maps is that scattering at low angular separation always creates forward momentum. This is even true for α = 0, for which, we could prove it explicitly from a mechanical analysis of the collisions. This feature explains the stability of the polar phase in the absence of noise. For small α, the dis-alignment, present in the frontal collisions is strong enough to stabilize the isotropic state. For α > α * , this does not hold anymore, and the isotropic state becomes unstable. The coexistence of stability of the two states for α < α * guarantees the transition to be discontinuous.
The effect of noise is most easily understood from the shape of the curves µ(α) and ξ(α). The key point here is that in the present scheme of approximation -Boltzmann equation plus von Mises Ansatz -the computation of ξ is identical to that of µ, apart from the additional asymmetric factor (1/2 − cos ∆). This does not rely on the specificity of the microscopic model. As a result, the curve ξ(α) is similar to the curve µ(α), with a shift towards larger α. Adding noise simply shifts the transition towards larger α: at some point ξ changes signs at the transition, which becomes continuous, hence the tricritical point.
Finally, the re-entrant transition from polar back to isotropic at large α (while increasing α) is a direct consequence of the non-monotonous shape of µ(α). At moderate α, the effective alignment µ increases with the persistence ofn. However, a too large persistence reduces the alignment, because the scattering event does not reorientn significantly, and no memory of the collision is kept: v relaxes to practically the samen as before the scattering. Some alignment still takes place but is most easily destroyed by small amounts of noise. Interestingly, while the recollisions play no role at the first transition from the isotropic to the polar state, they here increase the effective alignement and thereby postpone the re-entrant transition to larger α.
As can be seen from the above discussion, the results do not depend on the details of the scattering maps, as long as binary scattering at low angles produces effective alignment. The partially integrated scattering function (Fig. 3e,f) is then positive for small angles. Models which fall into this class are (i) the here described hard disks with v andn, (ii) the inelastically colliding hard disks with only v [21], and (iii) the soft disks described in Ref. [17]. All three are self-propelled in the sense that their freeflight speed relaxes to v 0 .
It is important to note that the distinctive behaviour of the above class of system is markedly different from the original Vicsek model, in which binary scattering at low angle effectively disaligns. In the Vicsek model, the partially integrated scattering function [21, Fig. 2a] behaves precisely in the opposite way of the one here. It is therefore unlikely that a coarse-graining of self-propelled hard disks (or any of the above three models) would yield the Vicsek model. Altogether, the present new class of models, including the self-propelled hard particles, could well play for the theory of "simple active liquids" the role that hard spheres play for the statistical mechanics of gases.