# Order-disorder transition in repulsive self-propelled particle systems

###### Abstract

We study the collective dynamics of repulsive self-propelled particles. The particles are governed by coupled equations of motion that include polar self-propulsion, damping of velocity and of polarity, repulsive particle-particle interaction, and deterministic dynamics. Particle dynamics simulations show that the collective coherent motion with large density fluctuations spontaneously emerges from a disordered, isotropic state. In the parameter region where the rotational damping of polarity is strong, the systems undergoes an abrupt shift to the absorbing ordered state after a waiting period in the metastable disordered state. In order to obtain a simple understanding of the mechanism underlying the collective behavior, we analyze binary particle scattering process. We show that this approach correctly predicts the order-disorder transition at dilute limit. The same approach is expanded for finite densities, although it disagrees with the result from many-particle simulations due to many-body correlations and density fluctuations.

## I Introduction

Following the seminal works by Vicsek et al. Vicsek et al. (1995) and Toner and Tu Toner and Tu (1995, 1998), the interest of physicists in collective motion of self-propelled particle (SPP) systems has been growing in last two decades. The main aim of the field is to obtain an unified understanding over non-trivial ordering behaviors seen in a group of living organisms in various length scales, from flocks of birds and schools of fish Ballerini et al. (2008); Cavagna et al. (2010); Katz et al. (2011); Bialek et al. (2012); Gautrais et al. (2012); Lopez et al. (2012) to tissue cell migration in monolayers Szabó et al. (2006); Poujade et al. (2007); Szabó et al. (2010); Anon et al. (2012); Serra-Picamal et al. (2012) and bacterial colonies Zhang et al. (2010); Peruani et al. (2012); Wensink et al. (2012) and even to the intracellular structure such as actin filaments and microtubules Schaller et al. (2010); Sumino et al. (2012). A well established approach in the field is the Vicsek-type SPP model Vicsek et al. (1995); Grégoire and Chaté (2004); Chaté et al. (2008). It assumes the particles drive themselves with a constant speed while adjusting their direction of motion parallel to their neighbors’ velocities. It has been shown that the non-equilibrium character of the systems leads to development of long-range order and giant density fluctuations, which are unusual in equilibrium systems Toner and Tu (1995, 1998); Aditi Simha and Ramaswamy (2002); Ramaswamy et al. (2003); Chaté et al. (2006).

Experiments using vibrated grains and driven colloids Narayan et al. (2007); Aranson et al. (2008); Deseigne et al. (2010, 2012); Buttinoni et al. (2013) suggest that the similar properties can also be found in non-living systems. However, the Vicsek model is not likely to illustrate the microscopic nature of these systems, where particles do not “compute” the average of their neighbors’ velocities. It is rather natural to assume that the elements interact with repulsion through excluded volume. Pedestrian movement, another major target of SPP studies, is also known to be well described by models that assume short-ranged repulsive force from neighboring people Helbing and Molnár (1995); Helbing et al. (2000). Therefore, it is of great importance to capture dynamics of SPP systems where repulsive interactions are dominant.

Previously, Hanke et al. Hanke et al. (2013) proposed a model of soft, repulsive active particles, explored its collective behavior over parameters, and found that a polarized state emerge at certain parameter region. However, their model has some inconsistencies at microscopic level when compared to actual granular or colloidal matter. First, the particles are not always as soft and strongly overdamped as the model requires them to be in order to achieve the collective motion. Second, the origin of the noise is unclear, since thermal fluctuation is not relevant in the length scale we deal with.

Here we consider a simple model of repulsive SPP that can be applied to various systems from motile cells, active colloids and grains to pedestrian crowds. This generalized model obeys deterministic equations of motion and interacts with other particles by short-range repulsion. We explore the phase diagram and discuss the nature of the phase transition. In section 4 and 5, we analyze the emergence of collective motion from binary scattering process.

## Ii Model

We consider disk particles of equal radius in a two-dimensional continuous surface. The dynamics is governed by the following deterministic equations Hiraoka et al. (2015):

(1) | |||

(2) |

Eq. (1) describes the Newtonian equation of motion with velocity . The first term of the rhs is the self-propelling force of fixed magnitude along the direction of “polarity”, an internal degree of freedom defined by an unit vector . The second term is dissipation proportional to , whose strength is determined by . The interaction force is given by binary, short-ranged repulsion: here we assume linear Hookean contact, if (in contact) and otherwise. Eq. (2) describes the time evolution of the polarity . When deviates from the direction of the velocity, it is rotated by a torque proportional to the deviation with a coefficient . Thus the equations include damping term for both translational and rotational degrees of freedom. The former is underdamped, while the latter is overdamped.

Note that each parameter gives a different characteristic time: is the time scale that a particle at the stationary speed takes to run its own diameter; is the relaxation time of speed; is the characteristic time during which two colliding particles are in contact; is the relaxation time of polarity. Without loss of generality, we set unit of length and time as and , and obtain rescaled equations.

## Iii Numerical results

### iii.1 Ordering behavior

We perform particle dynamics simulations of a system in a square box of size with periodic boundaries. Initial positions and polarities are randomly assigned; overlaps between particles are reduced prior to each run. Unless otherwise specified, we discuss results for and . Under such choice of parameter values, the elasticity is large enough to avoid unphysical situations where particles in contact pass through each other. A fourth-order Runge-Kutta method is employed for the numerical integration. The time step size is chosen to be sufficiently small compared to both and , either of which defines the shortest time scale of the dynamics in the system.

As shown in Fig. 1, the system exhibits polar ordering and large density fluctuation, the two characteristics also seen in the Vicsek model. The polar order is characterized by the average normalized velocity,

(3) |

If the system is in random state , while for a perfectly ordered state. In order to quantify the density fluctuation, we divide the system into small cells of size and take the standard deviation of the local packing fraction,

(4) |

where denotes the local packing fraction in the cell .

Fig. 1 shows a typical ordering behavior. Initially, the system is randomized in position and in polarity, so both and have small values. When the ordering starts, relatively dense and locally ordered regions appear and grow to form several clusters. The clusters eventually merge with each other until all the particles move into an identical direction. The time series of and that of display a simultaneous increase.

### iii.2 Phase diagram

We explore the behavior of the model by implementing a set of simulations with different packing fractions and rotational damping parameter . We find that the globally aligned state emerges in the regime where system is dense and rotational damping is weak, while the disordered, isotropic state persists if we set the packing fraction small and the damping parameter large. For , where polarities are never rotated from initial randomized condition, the system maintains a trivial disordered state; however, a small but finite value of leads to a slow ordering.

We construct the phase diagram by performing a set of runs (typically 8 to 16) with different initial configurations for a certain simulation time . If polar order, namely , is established for one or more runs, the parameter set is classified as a part of ordered region; otherwise, it is in the disordered phase. We choose and and depict Fig. 2.

In the ordered region near the phase boundary, the system maintains the disordered state (small ) until it suddenly transits to the ordered state (). The lifetime of the disordered state, which we refer to as the waiting time, varies depending on the initial configuration. As we increase , tends to be longer and, eventually, ordering behavior does not take place within the simulation time for any realizations.

We look into the distributions that follows for each parameter set . In Fig. 3, the second central moment and third central moment are plotted against the average . In the large waiting time regime, namely where , two moments satisfy and , which is expected for exponential distribution. The fact that the waiting time follows an exponential distribution implies that ordering events occur as Poisson processes involving nucleation phenomena.

Next, we discuss how the phase diagram would be changed if we choose other system sizes and time scales. Finite-size effects are investigated by changing the number of particles while keeping the volume fraction fixed. In small systems, strongly depends on ; in larger systems, however, the increase becomes insignificant (Fig. 4, inset). We confirm that is sufficiently large to avoid the finite-size effects for the time scale that we deal with.

displays a rapid increase as a function of (Fig. 4). It is difficult to identify the function that fits the growth, but the phenomenological double-well potential picture, which will be described in the next subsection, implies that ordering behavior can occur in any parameter region if we wait long enough. If this is the case, the phase classification inevitably depends on the observation time. Fortunately, the rapid increase, which seems to be exponential or faster, also ensures that the results with an moderately long simulation time provide a good approximation of the results with longer time scales. In other words, the phase diagram would be changed only slightly even if an order of magnitude longer observation time is employed.

### iii.3 Dynamics of order parameter

In this subsection, we show that the dynamics of the order parameter is described as a motion in a double-well potential. This phenomenological picture is consistent with the Poissonian character of the ordering process. It also provides a possible explanation for the divergence of the waiting time: In the ordered regime far from the boundary, the system quickly evolves from an unstable disordered state to an ordered state; near the boundary, however, the disordered state becomes metastable and a nucleation process is necessary to escape from it.

When the system is in the disordered state, the global polarization vector fluctuates around . The microscopic origin of the fluctuation is change in polarity of the particles caused by collisions to their neighbors. We expect that the collision events are uncorrelated to each other, and the time development of can be treated as Brownian motion in the potential. Assuming that the potential is harmonic around , we expect that constitutes a two-dimensional Ornstein-Uhlenbeck processUhlenbeck and Ornstein (1930). The stationary distribution of the radial component should be given by a Rayleigh distribution,

(5) |

where the scale parameter is the characteristic amplitude of the fluctuation. In fact, the observed distribution can be fitted reasonably by a Rayleigh distribution (Fig. 5(a)), providing an estimate of as

(6) |

Aside from the fluctuation amplitude, the escape rate is determined by the “position” of the barrier in order parameter space, which can be estimated by preparing a system in which a fraction of particles are given the same initial polarity so that initial polarization has a finite value. The initial positions are uniformly distributed both for aligned and unaligned particles. The double-well potential picture suggests that the time development of the system depends on which side of the barrier the initial state is situated: Systems with smaller than a certain value relax to the disordered state and those with larger than the same threshold goes to the ordered state. The threshold indicates the position of the barrier (Fig. 5(b)).

Decreasing in the disorder region towards the transition point, the threshold decreases (Fig. 5(c)) and the fluctuation amplitude increases (Fig. 5(d)). These results imply that the escape rate increases and the average time until spontaneous polarization occurs will be shorter.

Due to the absence of noise, a fully-ordered system does not evolve back to a disordered state; the ordered state is an absorbing state. We also study the stability of the ordered state numerically. The initial polarities are set to be aligned except for a small fraction ( to of all particles) with random directions. The initial positions are, again, uniformly distributed both for aligned and unaligned particles. In all the realizations, the system quickly relaxes to a fully ordered state, even if we choose a parameter set deep in the disordered phase, such as . This result suggests that the ordered state is stable against perturbation.

## Iv Binary scattering at dilute limit

In this section, we give a simple explanation to understand the mechanism that underlies the ordering behavior shown in the previous section by focusing on the binary particle collision process, a method also employed in Hanke et al. (2013). Here we assume that the system is dilute () so that the collisions are uncorrelated with each other, and that both the velocity and the polarity are fully relaxed before every collision.

Let us consider a binary scattering process between particle and (Fig. 6). Since we assume the rotational invariance, the geometry of the moment of contact is fully specified by two scalar parameters: the impact parameter , where , and the relative angle . The impact parameter shows the perpendicular offset of the two bodies’ center of mass from head on collision. If the collision is head on whereas it is a miss if .

Instantaneous alignment of the two particles is characterized by two-particle polarization,

(7) |

which corresponds to Eq. (3) with . We measure the post-collisional two-particle polarization at a point where the polarities and the velocities are fully relaxed, and compare it to the pre-collisional polarization . The increment indicates the magnitude of parallel alignment caused by the scattering process.

Assuming the system is homogenous and isotropic, two particles should collide in the relative angle of with a probability proportional to the relative velocity ; the impact parameter should be uniformly distributed. The average tendency of binary alignment, as a function of , is then obtained by taking an integral weighted by the “scattering cross section,”

(8) |

where is a normalization constant.

The result shown in Fig. 7 indicates that the alignment tendency has a maximum at . For , which corresponds to the regime where angular relaxation is slow, goes to zero. For large , namely , has a negative value.

Qualitative explanation is as follows. If the rotational damping is weak, the polarity of two particles remain unchanged, so the directions of motion will eventually restored to the original direction. In contrast, if the damping is strong, the polarity rotates itself quickly to follow the change in the direction of motion. It is with an intermediate value of the damping strength that the motion of two bodies align parallel, due to the competing effect of repulsive collision and subsequent rotational damping.

This argument is consistent with the results obtained from the many-particle simulations: First, the ordering in many-body system is the fastest in the parameter region that maximizes the value of ; second, the transition point in a dilute system is in agreement with the zero-crossing point (Fig. 2). These agreements imply that the onset of the collective motion arises from iteration of binary scattering, although one will have to take into account many-body correlation for the late stage of the ordering process, where the isotropic and uncorrelated conditions no longer hold.

## V Binary scattering at finite densities

The analysis in the previous section is only valid for dilute limit, where the velocity and the polarity of the particles are fully relaxed before every collisional event. At finite densities, however, the relaxation may not complete between collisions (Fig. 6).

For the parameter values we choose ( and ), in the disordered phase, meaning that relaxation of polarity is much faster than that of velocity. We check this by measuring the speed and the deviation between direction of velocity and polarity of all the particles directly from the many-particle simulation. Fig. 8 shows that the distribution of is considerably broad compared to the polarity distribution, which has a narrow peak at the point where . We will therefore take into account the speed distribution while safely neglecting the angular deviation in the following.

We apply a similar binary scattering study as in the previous section, except that is not defined at the point where the relaxation is complete, but at the point where the particle reaches the distance of the mean free path (mfp) away from the binary contact.

To calculate the mfp, we assume that every particle moves at speed . The magnitude of relative velocity between two particles is

(9) |

If collisions are uncorrelated to each other, the mean relative velocity is

(10) |

Since the scattering cross section and the number density , the mfp is given as

(11) |

Suppose particles and of speed and come into contact and are scattered. We assume that the speed of the particles independently obeys an identical distribution . The speeds evolve to and when the particles travels the distance of mfp after the collision. The post-collision speed distribution can be written as

(12) |

where denotes the scattering operator.

If the system is in the disordered phase, the distribution is stationary, so the following self-consistent condition should be satisfied:

(13) |

where is the stationary distribution. According to the Perron-Frobenius theoremBerman and Plemmons (1979), is the eigenfunction of operator which corresponds to the eigenvalue one.

is numerically derived by mapping pre-collision speeds onto post-collision speeds . We divide the -space into bins and simulate a binary scattering process for the representative value for each bin to obtain the matrix . Applying the power iteration method, we yield as the eigenvector of .

Once is known, the averaged increment in the two particle polarization is calculated as

(14) |

where and is a normalization constant.

Again, by numerically calculating the point where crosses zero, we obtain the estimated phase boundary. However, the estimation deviates from the results of many-body simulations, as shown in Fig. 2. The binary scattering approach is based on three assumptions: (i) the collisions are uncorrelated to each other; (ii) the motion of the particles is isotropic; (iii) the system is homogeneous, i.e., the density does not fluctuate. We suppose that these assumptions basically contribute to suppress the emergence of ordering behavior and that the deviation with the results from many-body simulations stems from the failure of one or more of the assumptions. A more detailed analysis of this point, particularly a closer look into the local quantities, is left for future work.

## Vi Summary

In this work we consider a simple model of repulsive self-propelled particle systems. Numerical simulations indicate that the system exhibits polar-ordering and a transition between ordered phase and disordered phase; the phase diagram is depicted as a function of the packing fraction and the rotational damping strength. We show that the slowdown of the dynamics at the proximity of the phase boundary is caused by metastability of the disordered state, and that transition to a polarized state requires a nucleation process. The physics that induces the emergence of the collective motion is the iteration of binary collision and the rotational relaxation afterwards. This argument is supported by the fact that the binary scattering analysis predicts the transition point from the many-particle simulation in the dilute limit. For finite densities, the approach deviates from the actual boundary because the disordered state is unstable against the many-particle correlations, anisotropy of the collisions, and the density fluctuation in the system. Recently, much effort has been devoted to derive a hydrodynamic description directly from the microscopic dynamics of SPPBertin et al. (2006); Baskaran and Marchetti (2008); Peshkov et al. (2014); Ihle (2014). Applying the kinetic approaches to the model we present will be an interesting direction of future work.

This model is consistent with some other models presented in previous studies. (i) Soft particle models for epithelial cell migration developed in Szabó et al. (2006); Henkes et al. (2011) has a first order equation of motion instead of Eq. (1), namely , where is a constant speed. This model is equivalent to ours when we take the overdamped limit, with . (ii) The model mentioned in Vicsek and Zafeiris (2012) and explored in Hanke et al. (2013) does not have polarity degree of freedom and employs the equation of motion as follows: , which means self-propulsion is always directed towards the direction of motion. This model corresponds to another limit in our model, . (iii) The social force model for pedestrian movement Helbing et al. (2000) is given by following dynamics: , where denotes the individual’s desired direction. This equation of motion is analogous to case in our model, where polarity of particles are never changed, although the particular form of the short-range repulsive interaction is different. Hence, our model can be regarded as a generalized model that bridges the above three models. While these previous models incorporate noise as a control parameter, we show that the phase transition is realized even in the absence of noise. In conclusion, our findings elucidate the universality of collective ordering behavior in repulsive SPP systems.

## Acknowledgments

T. H. thanks Koji Oishi for helpful discussions. The authors also appreciate comments and suggestions from anonymous reviewers, which contributed greatly to improve the paper. This work was supported by CREST program of Japan Science and Technology Agency.

## References

- Vicsek et al. (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
- Toner and Tu (1995) J. Toner and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995).
- Toner and Tu (1998) J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998).
- Ballerini et al. (2008) M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Anim. Behav. 76, 201 (2008).
- Cavagna et al. (2010) A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale, Proc. Natl. Acad. Sci. U.S.A. 107, 11865 (2010).
- Katz et al. (2011) Y. Katz, K. Tunstrom, C. C. Ioannou, C. Huepe, and I. D. Couzin, Proc. Natl. Acad. Sci. U.S.A. 108, 18720 (2011).
- Bialek et al. (2012) W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Proc. Natl. Acad. Sci. U.S.A. 109, 4786 (2012).
- Gautrais et al. (2012) J. Gautrais, F. Ginelli, R. Fournier, S. Blanco, M. Soria, H. Chaté, and G. Theraulaz, PLoS Comput. Biol. 8, e1002678 (2012).
- Lopez et al. (2012) U. Lopez, J. Gautrais, I. D. Couzin, and G. Theraulaz, Interface Focus 2, 693 (2012).
- Szabó et al. (2006) B. Szabó, G. J. Szöllösi, B. Gönci, Z. Jurányi, D. Selmeczi, and T. Vicsek, Phys. Rev. E 74, 061908 (2006).
- Poujade et al. (2007) M. Poujade, E. Grasland-Mongrain, A. Hertzog, J. Jouanneau, P. Chavrier, B. Ladoux, A. Buguin, and P. Silberzan, Proc. Natl. Acad. Sci. U.S.A. 104, 15988 (2007).
- Szabó et al. (2010) A. Szabó, R. Unnep, E. Méhes, W. O. Twal, W. S. Argraves, Y. Cao, and A. Czirók, Phys. Biol. 7, 046007 (2010).
- Anon et al. (2012) E. Anon, X. Serra-Picamal, P. Hersen, N. C. Gauthier, M. P. Sheetz, X. Trepat, and B. Ladoux, Proc. Natl. Acad. Sci. U.S.A. 109, 10891 (2012).
- Serra-Picamal et al. (2012) X. Serra-Picamal, V. Conte, R. Vincent, E. Anon, D. T. Tambe, E. Bazellieres, J. P. Butler, J. J. Fredberg, and X. Trepat, Nat. Phys. 8, 628 (2012).
- Zhang et al. (2010) H. P. Zhang, A. Be’er, E.-L. Florin, and H. L. Swinney, Proc. Natl. Acad. Sci. U.S.A. 107, 13626 (2010).
- Peruani et al. (2012) F. Peruani, J. Starruß, V. Jakovljevic, L. Søgaard-Andersen, A. Deutsch, and M. Bär, Phys. Rev. Lett. 108, 098102 (2012).
- Wensink et al. (2012) H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Proc. Natl. Acad. Sci. U.S.A. 109, 14308 (2012).
- Schaller et al. (2010) V. Schaller, C. Weber, C. Semmrich, E. Frey, and A. R. Bausch, Nature 467, 73 (2010).
- Sumino et al. (2012) Y. Sumino, K. H. Nagai, Y. Shitaka, D. Tanaka, K. Yoshikawa, H. Chaté, and K. Oiwa, Nature 483, 448 (2012).
- Grégoire and Chaté (2004) G. Grégoire and H. Chaté, Phys. Rev. Lett. 92, 025702 (2004).
- Chaté et al. (2008) H. Chaté, F. Ginelli, G. Grégoire, F. Peruani, and F. Raynaud, Eur. Phys. J. B 64, 451 (2008).
- Aditi Simha and Ramaswamy (2002) R. Aditi Simha and S. Ramaswamy, Physica A 306, 262 (2002).
- Ramaswamy et al. (2003) S. Ramaswamy, R. Aditi Simha, and J. Toner, Europhys. Lett. 62, 196 (2003).
- Chaté et al. (2006) H. Chaté, F. Ginelli, and R. Montagne, Phys. Rev. Lett. 96, 180602 (2006).
- Narayan et al. (2007) V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 (2007).
- Aranson et al. (2008) I. S. Aranson, A. Snezhko, J. S. Olafsen, and J. S. Urbach, Science 320, 612 (2008).
- Deseigne et al. (2010) J. Deseigne, O. Dauchot, and H. Chaté, Phys. Rev. Lett. 105, 098001 (2010).
- Deseigne et al. (2012) J. Deseigne, S. Léonard, O. Dauchot, and H. Chaté, Soft Matter 8, 5629 (2012).
- Buttinoni et al. (2013) I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Phys. Rev. Lett. 110, 238301 (2013).
- Helbing and Molnár (1995) D. Helbing and P. Molnár, Phys. Rev. E 51, 4282 (1995).
- Helbing et al. (2000) D. Helbing, I. Farkas, and T. Vicsek, Nature 407, 487 (2000).
- Hanke et al. (2013) T. Hanke, C. A. Weber, and E. Frey, Phys. Rev. E 88, 052309 (2013).
- Hiraoka et al. (2015) T. Hiraoka, T. Shimada, and N. Ito, in Proceedings of the International Conference on Social Modeling and Simulation, plus Econophysics Colloquium 2014, edited by H. Takayasu, N. Ito, I. Noda, and M. Takayasu (Springer International Publishing, Cham, 2015) pp. 243–253.
- Uhlenbeck and Ornstein (1930) G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. 36, 823 (1930).
- Berman and Plemmons (1979) A. Berman and R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences (Academic Press, 1979).
- Bertin et al. (2006) E. Bertin, M. Droz, and G. Grégoire, Phys. Rev. E 74, 022101 (2006).
- Baskaran and Marchetti (2008) A. Baskaran and M. C. Marchetti, Phys. Rev. E 77, 011920 (2008).
- Peshkov et al. (2014) A. Peshkov, E. Bertin, F. Ginelli, and H. Chaté, Eur. Phys. J. Spec. Top. 223, 1315 (2014).
- Ihle (2014) T. Ihle, Eur. Phys. J. Special Topics 223, 1293 (2014).
- Henkes et al. (2011) S. Henkes, Y. Fily, and M. C. Marchetti, Phys. Rev. E 84, 040301 (2011).
- Vicsek and Zafeiris (2012) T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012).