Pattern phase transitions of self-propelled particles: gases, crystals, liquids, and mills

To understand the collective behaviors of biological swarms, flocks, and colonies, we investigated the non-equilibrium dynamic patterns of self-propelled particle systems using statistical mechanics methods and H-stability analysis of Hamiltonian systems. By varying the individual vision range, we observed phase transitions between four phases, i.e., gas, crystal, liquid, and mill-liquid coexistence patterns. In addition, by varying the inter-particle force, we detected three distinct milling sub-phases, i.e., ring, annulus, and disk. Based on the coherent analysis for collective motions, one may predict the stability and adjust the morphology of the phases of self-propelled particles, which has promising potential applications in natural self-propelled particles and artificial multi-agent systems.

The various collective patterns produced by self-propelled particles highly depend on inter-particle interactions. Reynolds [12] first proposed a widely-accepted interaction theory for natural biological flocks/ swarms, which comprises three heuristic rules: (1) Separation: avoid collisions with nearby group mates; (2) Alignment: match the velocity of group mates with moderate distances; (3) Cohesion: stay close to remote group mates. Later, Vicsek et al [1] introduced a flock model, the Vicsek model, based on inter-particle velocity alignment. Couzin et al [16] established a more comprehensive three-sphere model, and thereby reproduced three typical patterns of collective motions, i.e., congregation, flocking, and toruses. Inspired by the model of Couzin et al [16], Tanner et al [17] developed a model which leads to irregular collapse and fragmentation patterns. D'Orsogna et al [2] presented a method for predicting the stability and morphology of patterns based on a specific model developed for self-propelled particles. Zhang et al [18] investigated the dynamic mechanism of circular motion patterns of self-propelled particles, and developed a series of predictive control protocols [19] for such systems. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
© 2016 IOP Publishing Ltd and Deutsche Physikalische Gesellschaft A number of collective patterns have been observed in self-propelled particles, but the mechanism that governs inter-particle interactions to yield rich particle formations and behavioral patterns is poorly understood. One of the major unsolved issues is lacking of the capability to predict the eventual pattern and its stability merely based on the inter-particle interactions. Many physical, biological, and engineering multi-agent systems (MASs) have been investigated, and phase transitions have been revealed between different patterns [2,5,16,[20][21][22][23][24][25]. Among these investigations, Birnir [20] found the migratory (flocking) and stationary (circling) solutions of an ordinary differential equations derived from the Vicsek model [1]. Carrillo et al [21] discussed the existence condition of both single and double rotating mills. Single rotational structures are observed given general random initial conditions. However, the double rotating mills exist only with elaborately selected initial conditions, i.e., half particles rotate clockwise and half counterclockwise. Cheng et al [22] observed pattern transitions among gaseous, crystalline, and liquid flocks by slightly varying the zero crossing slope of the interparticle interaction function. Following the research line of the work of Gautrais et al [26], Calovi et al [27] exhibited a phase diagram of migratory and circular (milling) patterns of fish schools by varying the weights of two control terms. However, the collective patterns emerging from these models are still limited. In addition, a systematic investigation is still lacking in the phase transition mechanism underlying the abundant collective motional patterns.
In this paper, we found the phase transition mechanism between various geometrical patterns of self-driven particles by extensive numerical simulations on a 'minimal' model proposed by Vicsek et al [3,28]. The 'minimal' model achieves an ordered state merely by measuring the inter-particle relative distance, which minimizes the intra-system communication cost. Such a model is a double-integrator MAS that involves position, velocity, and acceleration. We obtained interesting phase-transition phenomena between gaseous, crystalline, liquid, and mill-liquid coexistence patterns by simply varying the vision range of self-propelled particles, rather than the weights of control terms [27]. Furthermore, we investigated the milling formation mechanism and identified sub-phase transitions, i.e., ring, annulus, and disk, according to the exponent of the repulsive function. In addition, the influence of noise is investigated by conducting simulation experiments with varying noise intensity levels. Finally, we analyzed the emergent pattern mechanism of self-propelled particles by energy dissipating theory. be the position and velocity of the particlei, respectively. We extend the 'minimal' model [3,28] by adding two parameters, α and β. The N self-propelled particles obey the stochastic ordinary equation below

Model
1, 2, , , x d is a standard Wiener process with amplitude σ,  i is the individual neighborhood, the scalars ( ) a = r d , d, α, a and b are the vision range, the expected inter-particle distance, the vision range parameter, the weight of the individual dynamics, and the weight of the inter-particle forces, respectively. Assume that the mass of each particle is one, the position-related pairwise force can be represented by Both the attractive and the repulsive inter-particle forces exist in ( )  F x ij for a > 1. That is because the repulsive force exists between agentsi and j provided that   < x d ij , while the attractive force exists between them otherwise (i.e., Specifically, when s = 0, model (1) becomes an ordinary differential equation and is noise-free. Model(1) is a typical Hamiltonian system with conservative energy, which belongs to a class of the canonical energy dissipative systems. To reproduce collective motions, model(1) can be generalized to a family of dynamics model [29] and g > 0. The 'minimal' model(1) is a special case of model(2) with a b = = 1. Throughout the paper, without loss of generality, we employ equation (1) as an example to simplify the discussion.
The phases of collective motions are quantified by two order parameters, i.e., the migratory motion order L is selected so that the density N L 2 equals 2. A sufficiently large density is selected to ensure that all agents are connected with each other directly or indirectly through the proximity network. Otherwise, the group will be split into several clusters with different sizes and moving directions, and hence the ordered state might not be attainable. When α increases, three remarkable transitions occur from gas (figure 1(a)) to crystal (figure 1(b)), and then to liquid (figures 1(c)-(k)), followed by mill-liquid coexistence patterns (phases). In the gaseous pattern, the particles are evenly diffused in the entire movement space with relatively large inter-particle distances, whereas in the crystalline pattern, the particles form a crystalshaped congregation with a strictly identical shorter inter-particle distance, which resembles the grid distribution of molecules in the crystalline phase. By contrast, in the liquid pattern, the particles form dropshaped congregations where the inter-particle distance typically decreases from the surface to the kernel. The specific 'drop' shape of the liquid pattern is defined by the 'surface tension' of the particle cluster. In the millliquid coexistence phase, a milling cluster or a liquid cluster may be formed given different initial conditions. In the milling pattern, particles rotate around a common center with one or several circular trajectories. Figures 2(a)-(b) show the migratory and circular order parameters, respectively. Evidently, the gaseous/ crystalline/liquid phases correspond to the migratory motion V m =1 and V c =0 (a > 1.6), and the milling phase suggests the stationary motion, V m approximates to 0 and V c =1 (a > 1.6). Both V m and V c curves bifurcate into the liquid phase (red line) and the milling phase (blue line) above a = 1.4, which verifies our observations in figures 2(a)-(b). In order to illustrate the phase transitions between the gaseous, crystalline, liquid, and mill-liquid coexistence phases, we show the evolution of the normalized inter-particle distance  figure 2(c). More precisely, as α ascends, D e decreases slightly from 0.9 around a = 1.2, thereby implying that more inter-particle connections are formed and that the liquid phase is entered. Next, the D e curve falls suddenly from 0.9 to 0.3 at around a = 1.4 as the interparticle distances decrease abruptly. The D e curve bifurcates into the liquid phase (red line) and the milling phase  . D r is a ratio between the inner and outer radius over all particles of  i (the neighborhood of particle i). The milling particle clusters can be further classified into three sub-phases, D r =0 indicates a disk (solid mill with multiple loops, see figure 1(g)), ( ) Î D 0, 1 r refers to an annulus (hollow mill with multiple loops, see figure 1(h)), and D r =1 is a ring (single loop, see figure 1(i)). Phase transition between the ring and the annulus phase is at b = 0.7 and that between the annulus and the disk phase is at b = 1.3, while α is greater than 1.4. To quantitatively investigate the influence of β on V m , V c , D e , and D r , we show the results along with both parameters α and β in figure S1 of supporting information (SI). It is observed that, with a fixed α, V m and V c almost keep the same level with increasing β. Meanwhile, D e and D r grow larger along with increasing β, which supports the role of parameter β on the group pattern.

Results
We present the entire phase diagram for the self-propelled particles in figure 3, where the gaseous, crystalline, liquid, and mill-liquid coexistence phases correspond to , respectively. It should be noted that in the mill-liquid coexistence phase, the liquid pattern has similar positional distribution to the corresponding milling pattern shown in figures 1(c)-(f), but they differ because their velocities are synchronized (migratory) rather than being rotational (stationary), as in the mill phase (see figures 1(g)-(k)). To show the phase transition more vividly, we present a heat map diagram in figure 4 that illustrates the phase transition of the MAS, which can be used to predict the emergent phase obtained according to the parameter pairs of α and β. Given a sufficiently large constant particle distribution density r = = N L 2 2 , the heat maps for different particle numbers, = N 30, 50, 100, and 400, are identical. Some substable phenomena still occur in the mill-liquid coexistence phase. For example, in the milling pattern, several smaller mills with the same or different rotational directions are observed, because the connectivity of particles is not guaranteed due to the limited sensing/interacting range. Reference [21]claimed that general random initial conditions typically yield single rotational structure. However, double rotational ring pattern is also found given random initial conditions with velocity = v 0.5 0 , which is a supplement to [21]. Another sub-phase is observed on the boundary between the liquid phase and the ring sub-phase ( [ ] a Î 1.4, 1.5 and [ ] b Î 0.2, 0.8 ), where clumps of particles rotate around a common center. Although [20] showed the existence of both the migratory and stationary solutions to self-driven particles, the coexistence of liquid (migratory) and mill (stationary) does not always happen. For some settings, only mills are formed regardless of the employed initial conditions. The possibility of forming a liquid or a mill phase may be affected by the parameters α, β, and N, which requires further investigation.
Next, we investigate the phase transitions among the gaseous, crystalline, and liquid phases. As shown in figure 3, the particle density in the three phases increases along with the vision range α. In particular, we observe that the gaseous and crystalline phase are yielded with a small vision range α (  a 1.2), whereas the liquid phase is formed when a > 1.2. For a fixed equilibrium pairwise distance d, the attractive force range rises with increasing α. It is expected that, with a smaller effective attractive force range, the potential energy accumulates more slowly with increasing particle density. In particular, let U be the potential energy, which is referred to as H-stable if a constant  g 0 exists such that  g -U N [2,31]. By extensive numerical simulations, it is observed that the particles in a system with H-stable potential energy form the gaseous and crystal phases. With greater potential energy accumulation due to a largerα, the self-propelled particle system undergoes a phase transition to the liquid phase and then the liquid-mill coexistence phase, of which both are not H-stable.
The formation of different milling sub-phases emerge with various short-range interactions among particles. As shown in figure 5, the absolute value of the inter-particle force ( )  F x ij grows stronger along with the increasing exponent β (see equation (1)) in the short inter-particle distance range of ( )    Î x 0.5, 1 ij . Thus, each particle is more likely to capture neighboring particles to form an annulus or a disk. Particles with a small value of β balances the attractive and repulsive forces with smaller average inter-particle distances, and thus rings emerge. However, the inter-particle repulsive force grows much faster with larger values of β when    < x d ij . When β ascends beyond a threshold b (e.g., b = 3), the attractive/repulsive force balance cannot be  The range of repulsive force is fixed to d, α is the vision range for attractive force and β denotes a particle's tendency of obstacle avoidance. The gas, crystal, and liquid phases are distinguished by D e . The co-existence phases are measured by D r . maintained in a single ring because a larger inter-particle distance is required. As a result, annuli with multiple rings are formed. The parameter β represents a particle's tendency of obstacle avoidance.
The radius R of a milling ring can be calculated by balancing the centrifugal and centripetal forces in the selfpropelled particles model(1) with   = v v i 0 , thereby leading to: With the same parameter settings, the radius of liquid (migrating) rings is different from that of milling (circling) rings due to the different motional dynamics. The mathematical proof of equations (3) and (4) is given in SI. The phases and their transitions discussed above are noise-free (s = 0). Next, we conduct investigation to the stochastic model (1) in the presence of noise (s > 0). All the initial conditions and parameters are the same as figure 2 except that a = 2, b = 1, L=4, and  s 0. Since V c is fluctuating in the presence of noise, we measure V c instead of V c , which is the average V c over the second half of each simulation run. converges more slowly with increasing external noise intensity. Once the model converges, the formation patterns with noise are the same as the noise-free case, as shown in figures 1(e) and (i). When σ is in the range of [ ) 0.06, 0.07 , the index G drops; bi-model [27] is observed as the particles transits between the liquid and the milling annuli. When σ is in the range [ ] 0.07, 0.1 , the index G=0 and particles still form a cluster but converge to neither a liquid nor a milling ring. When σ grows lager than 1, noise dominates the collective behavior, i.e., particles span the whole space and move randomly. The universality class of the self-propelled system should satisfy the following conditions: (i) individual stable dynamics; (ii) inter-particle attractive and repulsive forces; (iii) the union of the proximity networks of consequential time instants keeps connected. Generally speaking, condition (ii) can be described as a K-classed function, and condition (iii) is fulfilled by sufficiently large particle density ρ and periodical/bouncing boundary conditions. External noises are allowable in the systems as well.

Conclusion
In summary, this study refines our understanding of the pattern transitions in self-propelled particles as well as natural biological collective behaviors. Multiple formation patterns are observed merely by varying two factors: vision range α and tendency of obstacle avoidance β. The phase transitions that occur between four different MAS phases might improve the flexibility of unmanned vehicles or multi-robot coordination technology. Vehicles/robots can be programmed to cross over from stable gas/crystal/liquid to mill behaviors, thereby shifting from an undesirable mode to a specific beneficial pattern solely by slightly changing a few dynamic parameters. The mechanism that underlies these interesting phenomena can be explained by H-stability analysis of Hamiltonian systems. This phase transition can be expected to scale up to more complex large groups with inter-particle attractions and velocity alignments. From an interdisciplinary biological/physical perspective, this study helps bridge the aggregation phenomena that occur in physical particles and the swarming of organisms in nature by identifying their common underlying mechanism. From the perspective of industrial applications, the control of MAS in engineering may be adjusted more efficiently by the phase transition method because it minimizes the cost for intro-group communication.