Velocity alignment promotes motility-induced phase separation

We study the phase behavior of polar Active Brownian Particles moving in two-spatial dimensions and interacting through volume exclusion and velocity alignment. We combine particle-based simulations of the microscopic model with a simple mean-field kinetic model to understand the impact of velocity alignment on the motility-induced phase separation of self-propelled disks. We show that, as the alignment strength is increased, approaching the onset of collective motion from below, orientational correlations grow, rendering the diffusive reorientation dynamics slower. As a consequence, the tendency of particles to aggregate into isotropic clusters is enhanced, favoring the complete de-mixing of the system into a low and high-density phase.

Active systems made of self-propelled particles have been (and still are) the object of intense research [1][2][3]. Such an interest is mainly due to the fact that, unlike passive systems, active systems are driven out-of-equilibrium at the level of a single constituent, giving rise to a rich and novel phenomenology. In particular, activity may trigger non-equilibrium large-scale behavior like clustering and phase separation in the absence of cohesive forces, or the spontaneous emergence of collective directed motion. For instance, bacterial colonies self-organize into coherently moving swarms [4][5][6], self-propelled Janus colloids form dynamical clusters [7][8][9][10][11][12] and actin filaments selforganize into dense traveling structures at high density [13,14]. Similarly, large-scale directed motion has been observed in monolayers of polar grains and colloidal rollers suspensions [15,16].
In order to gain understanding on the general mechanisms underlying such collective phenomena, several simplified model systems have been introduced and analyzed in depth. In particular, the role played by the nature (or symmetry) of the microscopic interactions between agents on phase transitions and structure formation has been extensively discussed. Among these minimal models, the Active Brownian Particle (ABP) and Vicsek model, together with their variants, occupy a privileged place. The ABP model [17][18][19][20][21][22][23][24] introduces the simplest interaction mechanism between isotropic particles, i.e. volume exclusion, while the Vicsek model [2,[25][26][27][28][29][30][31] considers a local 'ferromagnetic' velocity alignment mechanisms between otherwise non-interacting polar agents. As such, they constitute the prototypical, and therefore most studied, models of isotropic and polar active matter. In the ABP model, the tendency to form dense clusters arises from the 'mere' combination of self-propulsion and excluded volume interactions. These clusters grow as self-propulsion or density is increased and eventually lead to a complete phase separation, the so-called Motility-Induced phase separation (MIPS) [32] originally reported in the context of run-andtumble particles [33]. In the Vicsek model, the direction of self-propulsion of a particle tends to align towards the mean direction of its neighbors. Such an alignment interaction triggers an order-disorder phase transition [2,25,26] characterized by the emergence of different kinds of traveling structures [27][28][29][30].
In experimental set-ups, interactions between components obviously involve more ingredients than the simplified ones captured by the ABP or Vicsek models. Therefore, the impact of the interplay between these two paradigmatic interactions, excluded volume and velocity alignment, on the large-scale behavior of active systems p-1 arXiv:1807.07497v1 [cond-mat.soft] 19 Jul 2018 appears as a natural question which has motivated the introduction of several models [34][35][36][37][38][39][40]. However, such a question has been mainly addressed from continuum models where excluded volume effects are captured through an effective coarse-grained density-dependent self-propulsion velocity [37,39]. While this mesoscopic description is wellgrounded for isotropic active matter, for which it can be systematically derived from the microscopic dynamics [20,21,41,42], such a connection between the microscopic and mesoscopic description is still lacking in the presence of alignment.
Here, in order to clarify the role played by alignment on the phase behavior of ABP, we keep a microscopic description of the interactions. We show that, below the onset of polar order (or flocking), the self-trapping mechanism responsible of MIPS is maintained in the presence of velocity-alignment. We show that alignment favors the emergence of MIPS, in agreement with the continuum formulation in [39]. Our simulations, complemented by a simple kinetic argument based on [19], allow us to understand the microscopic origin of such behavior: alignment induces orientational correlations which slow down the random reorientation of the particles, thus increasing the persistence of their motion.
Model and methods. -We study a two dimensional system of N self-propelled particles in a L × L box with periodic boundary conditions [43]. Each particle at position r i (t) = (x i , y i ) is self-propelled with constant speed v 0 along the direction given by n i = (cos θ i , sin θ i ). The over-damped Langevin equations governing the evolution of the system arė where η i and ν i are two independent white Gaussian noises with zero mean and unit variance: η α i (t)η β j (t ) = δ ij δ αβ δ(t − t ) and ν i (t)ν j (t ) = δ ij δ(t − t ) (α, β denoting cartesian coordinates). The noise η i represents a thermal bath at temperature T . The diffusivity D 0 and mobility µ verify D 0 = µk B T . Particles are subjected to both excluded volume and velocity alignment interactions. The inter-particle force, F i = − j =i ∇ i U (r ij ), accounts for short-range repulsions which we model by a repulsive potential of the form U (r) = u 0 (σ/r) 12 , with an upper cut-off at 3σ. Velocity alignment is introduced as a torque in terms of the sum of the phase difference and its strength is controlled by the 'ferromagnetic' coupling parameter K ≥ 0 [27,37,39,40,[43][44][45]. The sum in eq. 2 runs over the particles in the vicinity of particle i, denoted ω i , defined by the interaction range, R, which is chosen to be slightly larger than σ. This mimics an alignment mechanism which is of non-steric origin, meaning that particles do not need to be in contact in order to align, as in [34,35,38]. As such, it is closer to the alignment mechanism of flocking birds or hydrodynamically coupled micro-swimmers [2,16,46].
It is convenient at this stage to identify the relevant set of dimensionless parameters. The parameters σ, τ = D −1 θ and u 0 provide the natural units of length, time and energy of the model, respectively. Here, we fix R = 2σ and D θ = 3D 0 /σ 2 . Then, we are left with the reduced coupling parameter g = K 4πσ 2 D θ quantifying the strength of the alignment interaction with respect to angular diffusion; the Péclet number Pe = v0 σD θ and the mean density φ = ρπ/4 (where ρ = N/L 2 ). The ratio Γ = u0µ σv0 accounting for the competition between the softness of the potential and the self-propulsion strength should, in general, be considered as well. However, in the present study we focus on a stiff potential which diverges at contact and in a regime of (small enough) v 0 which ensures that repulsive forces always dominate over self-propulsion at short distances.
For g = 0 our model reduces to the paradigmatic ABP model, for which it is known that at high enough density and self-propulsion velocity it undergoes MIPS [18,19,[21][22][23][24]. For u 0 = 0, in turn, the model is equivalent to a variant of the Vicsek model in continuous time [37,39,43,45], which exhibits a transition towards collective motion above a critical value of the coupling g. Thus, our model includes the main ingredients of two paradigmatic models of active matter, introducing a direct coupling between the velocity and the density of aligning active particles.
Here, we numerically solve the equations of motion 1 and 2 using a Euler integrator with time step ∆t = 10 −3 . We simulate systems of N = 4000 particles at fixed φ = 0.4, D θ = 0.005 and µ = 1. We explore the phase behavior of the system in the g−Pe plane by varying v 0 from 0 to 0.3 (thus implying Pe ∈ [0 : 60]) and K from 0 to 0.6 (thus g ∈ [0, 5]). To analyze the long-time behavior of the system, we let it evolve from a homogeneous initial state until it reaches stationarity. The measurements reported have been obtained by averaging over 10 4 independent steadystate configurations.
Flocking. -We start our study by analyzing the onset of collective motion, or flocking transition, induced by the presence of velocity alignment interactions. As it is manifest from the model equation (2), the tendency of the particles to align their direction of motion competes with rotational noise. Thus, at small g, one expects noise to dominate, destroying the tendency to locally order. At higher g, the strength of the interaction might be large enough to overcome the random reorientation of the particles and eventually lead to large-scale ordering in the space of velocities. Global order in the system can thus be characterized by the average polarization, P = || 1 N i n i || (where * denotes an average over different steady-states). As shown in Fig. 1, the system orders as we increase g, from where we can infer the emergence of directed motion (flocking): a macroscopic p-2 Velocity alignment promotes motility-induced phase separation fraction of the system moves, in average, along a preferred direction. We identify the onset of flocking g c by P (g = g c ) 0.35. As shown in Fig. 1, the onset of flocking seems, within our numerical accuracy, independent of the choice of Pe and equal to g c ≈ 0.34. This result is in agreement with previous particle-based simulation studies [40] and with theoretical treatments that consider the coupling between velocity and density at a coarse-grained level [37,39]. In [37], the authors use a mesoscopic description of excluded volume effects, encoded in a densitydependent velocity, and find that the onset of flocking in their effective theory is located at g c = 2 4πσ 2 ρ , here, g c = 0.312, independently of Pe. This predicted value is reasonably close to our numerical estimate, although we find it slightly above. As noted in [37], this is consistent with the fact that the theoretical prediction is based on a linear stability analysis which can only access the limit of stability of the disordered state or spinodal line.
Phase separation. -In the absence of alignment (g = 0), particles slow down due to collisions, resulting in the decrease of particles' velocity with increasing local density. This creates a positive feed-back by which more particles accumulate in denser regions and thus, slow down, eventually triggering a complete phase separation between a dense, slow-swimming phase and a dilute, fastswimming one. This feed-back mechanisms leading to MIPS is usually described at a coarse-grained level by an effective density-dependent velocity v(ρ), that decreases sufficiently fast with the density [20,21,32,33,47].
We identify the onset of MIPS by the emergence of a large dense cluster (large in the sense that it represents a finite fraction of the total system). We define a cluster as a connected set of particles distant of less than 1.4σ (corresponding with the location of the first peak of the pair distribution function). Once clusters have been identified, we can compute the fraction of particles belonging to the largest one, Ψ. We thus use Ψ as a phenomenological order parameter to identify the onset of phase separation [21,23,40]. As shown in Fig. 2, Ψ becomes finite as Pe is increased. To systematically locate the onset of phase separation we consider the value of the Pe for which the fraction of particles belonging to the macroscopic cluster is Ψ = 0.35. According to this criterion, in the absence of alignment (g = 0), MIPS occurs for Pe > Pe * = 31 ± 2.
We focus now on the impact of alignment on the aggregation of ABP reported above. We proceed by slowly varying the coupling parameter g and computing the fraction of particles belonging to the largest cluster as a function of Pe, see Fig. 2. Since the flocking phase transition is approached from below (g < g c ), the system does not develop any global net polarization in this regime. However, as shown in Fig. 2, as the particle's tendency to align increases, the emergence of a large structure takes place at lower values of Pe. This means that it is easier for particles to aggregate in the presence of alignment. Is this aggregation phenomena in the presence of velocityalignement dominated by the feedback mechanism behind MIPS? Or is it rather controlled by the mechanisms responsible of structure formation in the form of traveling fronts in polar active matter?
In order to answer these questions and better understand the aggregation mechanism for 0 < g < g c , we analyze the dependence of particles' velocity with local density. We compute the average velocity projected in the direction of self-propulsion as a function of the density v(φ). We proceeded by subdividing the simulation box in cells and measure v j = 1 Nj i∈cjṙ i · n i , where N j is the number of particles in cell c j and φ j = πN j /(4c 2 j ).
p-3  The data obtained after averaging over 10 4 configurations is shown in Fig. 3, where the local mean velocity in the self-propulsion direction, v, is plotted as a function of the local density.
Below the onset of flocking, v(φ) decreases with increasing local density, indicating that particles are slowed down due to crowding effects, as expected from the MIPS scenario. However, above the flocking phase transition, directed motion emerges, such that a macroscopic (finite) fraction of the system travels coherently at the same mean velocity. In this regime the system is globally oriented and, thus, a macroscopic fraction of particles in the system do not block each other when they collide, but rather mutually align, giving rise to clusters that move coherently [40]. We conclude that, below the flocking phase transition, the particles' aggregation arises from the competition between self-propulsion and steric effects, through an effective density-dependent velocity. Therefore, the kinetic self-trapping mechanism responsible of MIPS survives in systems of polar ABP below the onset of flocking and it is responsible for the emergence of large clusters (as shown in Fig. 2). However, in contrast with mesoscopic approaches that consider a given functional decay of v(ρ) (usually exponential) [37,39], it is evident from the data Fig. 3 that alignment alters the decay of v(φ). Indeed, the linear decay of velocity with density of isotropic ABP [47] is maintained at finite g for low enough densities, but shifts to a slower, yet linear, decay at higher φ (see discussion below).
From the data shown in Fig. 2, we extract the dependence of the MIPS threshold as a function of the alignment strength, denoted Pe * (g) and shown by green symbols in the phase diagram Fig. 4. The vertical dashed line in the phase diagram shows the onset of flocking g c . Above it, in the ordered region, the formation of clusters cannot be attributed to MIPS but to the directed motion of polar selfpropelled particles (flocking phase denoted F). Below the onset of flocking, the system is controlled by the physics of ABP: at low self-propulsion the system is in a homogeneous gas phase (denoted H), while for Pe > Pe * (g) the system exhibits MIPS (phase separated region, denoted PS). The MIPS threshold Pe * (g) decreases with g: alignment favors phase separation. As illustrated by the snapshots Fig. 4, at g = 0, the orientation of neighboring particles, represented by a color code, is completely decorrelated. As g increases, approaching the flocking phase transition from below, orientational correlations grow, as evidenced by the growing size of the (polarized) singlecolor domains of particles. The growth of such orientational correlations with g is the crucial feature that explains why velocity-alignment facilitates MIPS.
To quantify local orientational correlations we compute the correlation function C s (r) = n i (t s ) · n j (t s ) |ri−rj |=r (3) in the steady state (i.e. t s τ ). In Fig. 5(a) we show C s for several couplings across the onset of flocking. As expected, spatial correlations grow as g increases, indicating the appearance of polarized domains. Above the onset of flocking, C s does not decay to zero, indicating that the system is globally ordered. To quantify the impact of the growth of spatial correlations on the orientational dynamics of particles, we also compute orientational correlations in time, captured by the self-correlation function As expected, time correlations decay exponentially and the relaxation becomes slower for g > 0 [see Fig. 5(b)], meaning that the orientation of particles has a longer memory in the presence of alignment or, put in different words, their dynamics is more persistent and their random reorientation slower. From a kinetic point of view, the slowdown of the orientational dynamics due to the appearance of polarized domains should have an impact on the rate of absorption and evaporation of particles in clusters. Indeed, when two particles collide they get stuck until they reorient in a time scale that, as shown in Fig. 5(b), is larger in the presence of alignment. Following [19], the rate of absorption of particles in a cluster can be expressed as where ρ l is the number density of the low-density gas background, and the rate of evaporation where D θ (g) is the effective rotational diffusion coefficient which, of course, generically depends on g. The parameter κ quantifies the total number of particles that leave the cluster in each escape event. The rate of evaporation of particles, k out , depends on the time D θ (g) −1 needed for a particle in the surface to randomly reorient, and pick an orientation that allows it to move away from the cluster. Since the orientations are correlated for longer times as alignment strength increases, k out will decrease with g, while the rate of absorption k in remains unchanged. As a result, the aggregation of particles into clusters is enhanced, favoring the system's phase separation.
In order to account for the dependence of the onset of MIPS with g, we extend the kinetic approach given by Redner et al. [19] to our case with velocity alignement. In the steady-state, the rate of absorption and evaporation should be equal, leading to the gas density The average density of the system can be expressed as where f is the fraction of particles belonging to the dense phase (i.e. the largest cluster) and ρ c its number density. We identify the onset of MIPS with the limit of coexistence f = 0, such that eq. 8 together with eq. 7 leads to an expression of the onset of MIPS as a function of the coupling strength: All the g-dependence is encoded in the effective rotational diffusivity. Since, as already discussed, D θ decreases with g, this simple argument allows us to rationalize why the onset of MIPS is shifted towards lower Pe in the presence of alignment.
In order to make a quantitative comparison between our simulation results and the prediction from the kinetic mean-field model, we use the onset of MIPS extracted from the simulations at g = 0 as an input in eq. 9, with D θ = D θ , to estimate κ ≈ 5. We further assume κ to be independent of g and measure D θ from the long time behavior of the angular mean-square displacement This quantity is plotted in Fig. 5 (c), where it is shown that D θ decreases with g, as expected from the previous discussion of C(t). From the measurement of D θ , we estimate the onset of MIPS, Pe * (g), using eq. 9. Such estimation is shown in the phase diagram Fig. 4   criterion on Ψ (see Fig. 2) are in considerable agreement and show the same overall behavior: phase separation can be triggered by a mechanism analog to MIPS of apolar ABP in the presence of alignment. As illustrated in the sequence of snapshots Fig. 4, increasing g can drive MIPS at constant Pe. As shown in Fig. 3, the decrease of v(φ) below the onset of flocking shows two clearly distinct regimes. In the dilute regime, the decay is independent of the alignment interaction between particles (as expected since in dilute regions particles are on average distant of more than the alignment interaction range). Conversely, at high packing fraction, the decay is slower in the presence of alignment, showing that the latter has an important impact on the dynamics of the system, even below the onset of flocking. Strikingly, the linear decay slope of v(φ) at high densities is, up to numerical accuracy, independent of the precise value of g. There are therefore two well-defined distinct regimes, indicating a change on the large-scale properties of the system as the density is varied. This suggests that such crossover is associated to the emergence of a MIPS, which, in the presence of alignment, is characterized by a dense phase with enhanced orientational correlations.
Conclusion and Discussion. -We find that polar ABP subjected to local alignment interactions undergo a MIPS-like phase separation, in that the aggregation mechanism below the onset of polar order is controlled by the competition between self-propulsion and steric effects. Moreover, the formation of non-polar clusters is favored when increasing the alignment interaction, thus fostering the complete phase separation of the system. Such behavior is understood as arising from the enhancement of orientational correlations which largely suppresses the random reorientation of the particles. The microscopic character of our approach allows to explore the kinetic mechanisms responsible of clustering in systems of ABP with velocity alignment and critically asses the appropriateness of the MIPS picture in this context. The understanding of the two different regimes in the decay of the velocity we identified from a coarse-grained theory calls for further developments, since the mesoscopic descriptions in terms of a single-mode density-dependent velocity do not correctly capture the competition between self-propulsion and excluded volume in the presence of alignment. * * * D.L. acknowledges received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie (IF) grant agreement No 657517. I.P. acknowledges MINECO and DURSI for financial support under projects FIS2015-67837-P and 2017SGR-884, respectively. All the authors acknowledge support from the COST Action MP1305-Flowing Matter.