Swirlonic state of active matter

We report a novel state of active matter—a swirlonic state. It is comprised of swirlons, formed by groups of active particles orbiting their common center of mass. These quasi-particles demonstrate a surprising behavior: In response to an external load they move with a constant velocity proportional to the applied force, just as objects in viscous media. The swirlons attract each other and coalesce forming a larger, joint swirlon. The coalescence is extremely slow, decelerating process, resulting in a rarified state of immobile quasi-particles. In addition to the swirlonic state, we observe gaseous, liquid and solid states, depending on the inter-particle and self-driving forces. Interestingly, in contrast to molecular systems, liquid and gaseous states of active matter do not coexist. We explain this unusual phenomenon by the lack of fast particles in active matter. We perform extensive numerical simulations and theoretical analysis. The predictions of the theory agree qualitatively and quantitatively with the simulation results.

Active matter is a substance comprised of active particle that demonstrate the motility, that is "the ability to exhibit motion and to perform mechanical work at the expense of metabolic energy" 1 . The active particles consume energy from the environment and drive themselves far from equilibrium 2,3 . Living matter (biological systems) provides an uncountable amount of examples of active particles systems 4 . Other examples refer to robotics 5 , biomedicine 6 , and social science 7,8 ; an interesting laboratory realization of artificial active particles are the so-called vibrots [9][10][11][12][13] . Similarly to common molecular systems active matter undergo different phase transitions -separation into dense and dilute phase 9,10 and crystallization [14][15][16] ; these phenomena my be described within the framework of conventional thermodynamics [14][15][16][17][18] .
The most prominent feature of systems of active (self-propelled) particle is the formation of self-organized coherent structures, see e.g. 2,3,[19][20][21][22][23][24][25][26][27][28][29] . Among these are intriguing milling patterns emerging in circular motion, when a group of individuals follow one another around an empty core. Such swirling patterns have been observed for animals at different evolution stages, ranging from plant-animal worms and insects to fish [30][31][32][33][34][35][36][37][38][39][40] . To describe the swirling motion several models have been proposed, including the celebrated Vicsek [41][42][43] and Vicsek-like models [44][45][46][47] , as well as models based on the inter-particle interaction potential 20,29,[48][49][50][51][52] . Although the inter-particle potential between active particles (agents) does not exists in reality, it mimics an intention of an agent to change its velocity according to some rules. It is convenient to put these rules in the form of Newtonian forces acting between the agents. Alternatively, the rules may be explicitly written in terms of velocity of an agent at the next time step, that is, determined by the velocities and coordinates of the (usually the nearest) neighbors. Below we will use the former approach, based on the inter-particle potential. Note that here we address the swirling motion of active particles; swirling also exists in a system of agitated passive particles, which have been studied in 53 .
The milling motion reported in Refs. 20,29,[48][49][50][51][52][53] referred to a single pattern where particles orbited around a common center. The number of particles was not large-it did not exceed a few hundreds. Formation of multiple milling patterns has not been reported. A possible reason for this could be a small system size due to the computational constrains. Indeed, the system of dynamic equations for active particles forms a set of very stiff ordinary differential equations (ODEs). An accurate convergent solution requires in this case an extremely small time step 54 . Hence it seems interesting to simulate much larger systems of active particles to explore the formation of multiple structures and their evolution.
In the present study we investigate large systems of active particles and report the formation of the conventional phases-gaseous and liquid phase, as well as solid phase, but also a new state of active matter-"swirlonic" state. In a gaseous state the active matter tends to occupy all the available space, while in a liquid state its volume is limited. However, in contrast to common molecular systems, we did not observe the gas-liquid coexistence, but only the presence of a single phase. In the swirlonic state multiple milling patterns-swirlons are formed. The swirlons behave like individual super-particles with surprising properties. They attract each other and coalesce www.nature.com/scientificreports/ upon collision, forming swirlons of a larger mass. When an external force is applied to a swirlon, it moves in the direction of the force, with a constant velocity, proportional to the applied force, as a particle in viscous medium. The steady velocity of the swirlon is inversely proportional to the mean-square velocity of the orbital motion of particles comprising the swirlon. The larger the swirlon, the larger the mean-square velocity and hence smaller the steady velocity. In other words, the mobility of swirlons decays with its mass. The swirlonic state is not a stable but transient state, since eventually all swirlons should collapse into a single milling structure. However the mobility of heavy swirlons drastically decreases, which entails, for a large system, a diverging transient time from the swirlonic state to the single milling state. This makes the swirlonic state of a special interest.
To study large systems of active particles we slightly modify the repulsion part of the commonly used Morse interaction potential. Namely, we apply the Gaussian dependence of the potential on the inter-particle distance, instead of previously used exponential dependence. This simple modification leads to a drastic reduction of the stiffness of the ODEs describing the system and keeps, at the same time, qualitatively the same behavior. The reduced ODE stiffness allowed to release very severe restrictions for the size of the computation time step and simulate relatively large systems, up to few tens of thousand particles. Simulation of such large system gives a new insight into the properties of active matter.

Model
We consider a two-dimensional system comprising N identical self-driven particles of mass m = 1 . Then the equation of motion of ith particle reads, where r i and v i are the coordinate and velocity of the i-th particle, f i is a force acting from all other particles, g i is an external force and ξ i (t) is a random force, modeled by a Gaussian white noise, where α and β denote the Cartesian components of the vector ξ i The self-driving force (b − v 2 )v describes the motion of the active particle with a constant velocity v = √ b in the absence of all other interactions. We use the following model for the interaction with other particles: U r (r) and U a (r) refer respectively to the repulsive and attractive part of the interaction potential with the according constant C r and C a characterizing the interaction strength. l r and l a describe the characteristic lengths of U r (r) and U a (r) . As we have already mentioned, we use the Gaussian dependence on distance for the repulsive part of the potential, instead of exponential dependence of the Morse potential. Eq. (2) has the form of stochastic Langevin equation. The corresponding Fokker-Planck equation for the distribution function �(r, v) is given in the Appendix A; it can be derived using the standard techniques, see e.g. 55,56 .
The dynamic of the system is described by N Eqs. (1)-(4) for i = 1, . . . N along with Eqs. (4)-(5) for the interaction potential. We simulate two-dimensional systems. In simulations we use up to 38,000 ODE, different sets of parameters (C a , l a , C r , l r ) and different initial densities of particles and system size. We start with random configurations and random velocities, randomized both in magnitude and direction. Depending on the above parameters one obtains rather different behavior of the system, which is qualitatively the same as reported in Refs. 29,52 , that is, flocking, swarming and milling. However, additionally, a novel state with many milling structures has been detected, which we call "swirlonic state". We have checked that this state emerges also for other potentials, with attractive and repulsive parts, including Morse potential. Still an accurate modeling of such states with multiple swirlons is computationally very challenging. Here we report our findings for the soft (Gaussian) form of the repulsive potential, see Eq. (5). The results have been obtained for a system in liquid, gaseous and swirlonic state confined in a circular region.

Phases of active matter
Starting from a uniformly distributed particles, confined in some area, Fig. 1a different phases may be formed depending on the density and parameters of interaction potential. If the density is low, so that the average distance between the particles exceeds the characteristic length of the interaction potential, the interactions are weak and compact patterns do not emerge. The particles run away from their initial position and occupy all the available space. They move chaotically, and thus form a "gaseous" phase, Fig. 1b. In contrast, for larger density, when many particles are located within the interaction distance, the attraction is stronger, and a state with a fixed density emerges, Fig. 1c. In this case the active matter does not occupy all the available volume, but only part of it. This is a "liquid" state. Interestingly, we did not observe coexistence of two phases-gaseous and liquid as in molecular systems. This may be explained by the fact that in contrast to molecular systems, where the distribution function Hence, while the fast particles are abundant in molecular systems, they are practically lacking in active matter. These fast particle can overcome the attraction forces of the molecular surrounding and give rise to the gaseous phase; the lack of fast particles in active matter results in the lack of the gas, coexisting with the liquid phase. If the particle density is still larger and the repulsive forces are strong, the active matter falls into a solid phase, Fig. 1d, where the particles mobility is suppressed. In some range of parameters a novel swirlonic phase is formed, Fig. 1e. It is comprised of swirlons-"superparticles" with many astonishing properties which we address below. Individual active particles in swirlons perform a swirling motion around their common center. As we demonstrate in what follows, the swirlons are formed when local fluctuating force exceeds the critical force, which characterizes the ability of an active particle to move against an applied load.
Finally, for large density and very strong attractive forces a collapsed state is observed, Fig. 1f. In the collapsed state active particle also move around the common center, however the character of the motion is rather irregular.
In the present study we will focus on the gaseous, liquid and especially on the swirlonic phase.
Qualitative theory of phase behavior critical force. To understand the conditions for the formation of different pases we consider first a single active particle under an action of a constant force. It obeys the equation (recall that m = 1):  www.nature.com/scientificreports/ For simplicity we assume that the direction of the force and velocity coincide, then Eq. (7) is a scalar equation. In Fig. 2a the dependence dv/dt = 0 on v is demonstrated for different f. It shows that the system evolves to states with dv/dt = 0 , which are depicted as stable points in the figure. These are the roots of the equation, bv − v 3 + f = 0 . For |f | ≤ |f c | there are roots, which correspond to the steady motion in the direction of the force and in the opposite direction. For |f | ≥ |f c | there is only one root, corresponding to the steady motion in the direction of the force. The critical force reads, Fig. 2a. One can interpret this results as follows: If the applied force f does not exceed the critical force f c an active particle can move against the force, that is, such force would not prevent the escape of the particle from some area, where such force is acting. In contrast, if f is larger than the critical force, an active particle can not escape from the region where this force is acting. The above estimate of the critical force will be used to estimate the conditions for different phase formation.
Conditions for "gaseous" phase. Now we estimate the average force acting on active particles at different locations. Let N active particles are uniformly distributed within a circle of radius R, see Fig. 1c (note that R may differ from the initial value R 0 ). If l a , l r ≪ R the average force, acting on a particle from its neighbors, in the internal part of the circle is zero. At the same time particles near the border experience uncompensated average force f av directed towards the circle center, see Fig. 2b. Consider a particle located exactly at the border, then the average force may be estimated as Here n = N/(πR 2 ) is the number density of particles, nrdrdφ gives the amount of particles within a small surface element rdrdφ which act on the chosen particle at the border with the force f (r) = (C a /l a )e −r/l a − 2(rC r /l 2 r )e −r 2 /l 2 r . The factor sin φ accounts for the fact that only the component of the above force, directed towards the circle center, contributes to the average force; the component normal to this direction does not contribute due to the symmetry, see Fig. 2b. As it follows from Eq. (9) both the repulsive and attractive forces contribute to f av , but these act in opposite directions. Since we focus in our study on the swirlonic states, which require that l a > l r (see also the classification of Ref. 52 ) we can neglect the repulsive force to find an upper estimate for f av . Then the integration over r in Eq. (9) yields, Writing the second line in the Eq. (10) we assume that R/l a ≫ 1 , which implies that only φ ≃ sin φ ≪ 1 contribute from the part of the integrand with the exponential factor; this also allows to change the upper limit π/2 to ∞ . Although the above estimate is justified for R/l a ≫ 1 , it works, in practice, rather well (with the accuracy better than 17%) already for R ≥ l a . We assume that a particle remains in a condensed phase if the average force is larger than the critical force. Then from Eqs. (10) and (8) (11) is not fulfilled, the "evaporation" takes place and the system evolves to the "gaseous" state, occupying all available space, see Fig. 1b. The evaporation is facilitated for small systems, at low density, with a weak, short-range potential. The evaporation through the flat boundary, R → ∞ , occurs at b 2/3 ∼ αC a nl a . With Eq. (8) for the critical force, the condition for condensed state may be written as where f * c = f c l a /C a is the reduced critical force and φ a = πnl 2 a is the effective "packing fraction", associated with the length of attractive interactions l a and particle number density n = N/� , where � = πR 2 is the area occupied by the active matter. The latter can differ from the initial area as well as from the total available area, see Fig. 1a and c.
Conditions for liquid and swirlonic state. Now we need to discriminate between two condensed statesliquid, with a uniform density, occupying a fixed volume, Fig. 1c, and swirlonic, with multiple milling structures, (11) C a nl a 1 − 3 4 www.nature.com/scientificreports/ Fig. 1e. That is, we need to find the condition, when a uniform state becomes unstable. We notice that swirlons emerge due to a "local capture" of active particles by their surrounding. Indeed, these structures emerge, when particles remain localized in the vicinity of their initial positions. In other words, it is expected that there exist local forces with the strength exceeding the critical force f c . Since the average force is zero for uniformly distributed particles, the arising forces have fluctuating nature. Let us estimate these fluctuating forces. Consider a particle in the internal part of the system. The number of neighbors that act on the particle is about ∼ nl 2 a and they act on the particle with the force of the order ∼ (C a /l a ) (again we neglect the repulsive forces). The average force is zero, but due to the density fluctuations, which may be estimated as ∼ nl 2 a a non-compensated fluctuating force of the order of ∼ (C a /l a ) nl 2 a ∼ C a n 1/2 will arise. A more rigorous analysis presented in Appendix B confirms the above estimates. Namely, it may be shown that the distribution of the local fluctuating force F in a space uniform system with the interaction potential U = C a e −r/l a reads as shown in Fig. 3 (see the "Appendix B"), The average force may be obtained from the above distribution (13) as This average force is to be compared with the critical force f c , (8), to escape from the local surroundings. Hence we can write the estimate for the condition that an active particle is locally caught by the fluctuating force and participate at the formation of a swirlon: where α 1 = O(1) is the constant of the order of one. The above condition (15) does not, however, discriminate between swirlonic and collapsed state, but rather demarcates liquid and swirlonic states. With the above notations the condition of liquid phase may be put into the form, The three phases of active matter-gaseous, liquid and swirlonic, obtained in simulations are demonstrated in Fig. 4a. If we assume that α 1 is considerably larger than α , we conclude that the observed phase diagram is in a qualitative agreement with the theoretical prediction (16).
conditions for collapsed state. Here we call a state with a single structure formed as a "collapsed" state.
This may happen if a considerable amount of particles does not experience an action of a local fluctuating forces, that cause the emergence of local swirlons, but rather experience a "global" force, directed towards the center of the occupied area. In other words, it is expected that collapsed structures would be formed, provided that number of "boundary" particles, located in a layer with the distance of l a from the boundary is comparable to the number of "internal" particles. This condition may be written as follows, If we neglect l 2 a as compared to Rl a we can write approximate condition for the formation of a single structure: As one can see from the above equation, the collapsed state emerges in small systems, with a long range attraction. This tendency has been confirmed in our simulations.
(15) C a n 1/2 > α 1 b 2/3 , (17) πR 2 n − π(R − l a ) 2 n ≈ π(R − l a ) 2 n.   Fig. 4b; still, they move coherently. Namely, in our simulations of active matter we observe that swirlons (in the swirlonic phase) move under an action of external force, as material particles in viscous medium. That is, their velocity linearly increases with the applied force. The mobility of a swirlon depends on its mass, that is, on the number of active particles comprising the swirlon. The larger the mass of a swirlon the smaller the mobility. Below we explain this astonishing property of swirlons. Let a constant external force g is directed along the x-axis. Let a swirlon be in a steady motion with a constant velocity v 0 , also along the x-axis. Consider a motions of an active particle inside the swirlon. Let F be the resulting total force that acts on the particle from the other particles of the swirlon. It supports the circular motion; the force and velocity components may be written as: where ω is the angular frequency and v 1 is the linear velocity of the circular motion . The x-component of the equation of motion then reads, Now we take into account that and average over time T which is large, as compared with the rotational period T swir = 2π/ω . Introducing the notation, we obtain,  22) and (23), which is a cubic equation for v 0 . Almost in all cases the constant velocity v 0 is much smaller than the orbiting velocity v 1 (that is, v 0 ≪ v 1 ), which implies the approximation, v 2 where v 2 is the mean square velocity of particles in the swirlon. Then the solution of Eq. (24) may be written as Note that the orbiting velocity v 1 may strongly differ from the "self-velocity" b 1/2 . Indeed, when a swirlon is formed, the initial potential energy of particles may transform into the kinetic energy associated with the orbiting motion around the center of mass of the swirlon. It is not easy, however, to obtain the value of v 2 , due to the lack of conservation laws for energy, momentum and angular momentum for active particles. In Fig. 5 the simulation results are compared to the theoretical prediction, Eq. (25). The theoretical results are in a very good agreement with the numerical data. coalescence of swirlons. In Fig. 6 we demonstrate the main stages of the swirlon emergence and evolution. One can observe the coalescent dynamics in the swirlonic phase. Initially the aggregation is fast, since the swirlons attract each other and have large mobility. In a course of time the average mass of swirlons increases and the mobility drops down (see the discussion below). This results in the slowing down of the coalescence. Eventually a small amount of very massive swirlons remain that have vanishing mobility and the aggregation practically ceases.
Consider some particular stage of the swirlon evolution. Let N s swirlons emerge in a system with the area S and the initial number density of active particles be n. Then the average "mass" of a swirlon (that is, the number of active particles forming a swirlon) reads, m s = nS/N s and the number density of the swirlons is n s = N s /S = n/m s . The average distance between the swirlons may be estimated as l ∼ n −1/2 s . The coalescence of these objects occurs as a process that symbolically may be written as A + A → A , where A denotes a single swirlon. We estimate the rate of this process. Let two swirlons were initially at the average distance r(0) =l and after time T coalesce, r(T) ≈ 0 . Their relative motion obeys the equation, where µ accounts for the mobility, which we approximate as size-independent and the factor m 2 s accounts that in each swirlon there are m s monomers, each attracting with the force (C a /l a )e −r/l a . Solving Eq. (26) with the above boundary conditions we find:  www.nature.com/scientificreports/ Since e¯l /l a ≫ 1 we conclude that The quantity T −1 may be treated as the reaction rate between two swirlons. Therefore the equation for the reaction kinetics A + A → A reads where we use the expressions m s = n/n s and l = 1/n 1/2 s . Solving the above equation we obtain for large t and small n s l 2 a → 0, where x = n s l 2 a or, asymptotically, for t → ∞, Eq. (28) predicts an extremely slow decay of the swirlons concentration, which agrees qualitatively with the simulation results, as it is illustrated in Fig. 7.
It is interesting to compare the coalescence of swirlons with the aggregation of vertexes in two-dimensional turbulence, although the background physics of these two phenomena is rather different. The vertex aggregation results in a power-law decay of a number of vertexes, n v ∼ t −ξ with ξ ≈ 0.70 − 0.75 57 , which is much faster than the inverse square logarithm of time in Eq. (28). This is explained by the long-range interaction of the vertexes, while the swirlon interactions decay exponentially with distance.

Conclusion
We reveal a new state of active matter, which we call a swirlonic state, where all active particles belong to swirlons, comprised of particles orbiting around common center of mass. Although such milling structures have been reported in many studies, the state of matter constituted by quasi-particles-swirlons was not investigated. The dimension of swirlons is almost independent on the number of particles inside it and determined by the intensity of the self-driven and inter-particle forces. These quasi-particles demonstrate an astonishing property-under an applied external force they move with a constant velocity, proportional to the force, just as the objects in viscous medium. We observe that the mobility of a swirlon is inversely proportional to the square average velocity of particles orbiting the swirlon's center of mass. Swirlons attract to each other and coalesce upon a collision, forming a joint swirlon. Therefore evolution of a swirlonic state of matter corresponds to continuous aggregation of particles when fewer and fewer swirlons of larger and larger mass are formed. This process being fast in the beginning, becomes extremely slow afterwards, leading to the formation of a rarified state of massive immobile swirlons; the evolution at this stage is practically frozen.
We also observe, that depending on the parameters of the inter-particle potential and characteristics of the self-driving force of active particles, a few other states may be formed. These correspond to the gaseous, liquid and solid states of matter, as it is known for molecular systems. While in the gaseous state the active matter occupies all the available space, in the liquid state only a part of the space is occupied. Moreover, in the liquid state the active matter possesses a surface tension. Surprisingly, we did not observe the coexistence between liquid and gaseous phase, as it is common for molecular systems. We explain this feature of the active matter by the lack of fast particles, which can escape from the liquid phase and form the gaseous phase. While the abundance of fast particles is secured in molecular systems by the Maxwellian velocity distribution, the velocity distribution for fast particles in active matter has much steeper decay with increasing velocity.
For high density of active particles with a small self-driving force and relatively strong repulsive interactions, a solid state is formed. It is characterized by a complete suppression of the long-range motion (diffusion) of the particles. In the solid state the active particles perform a very limited motion around a center of a potential well formed by their nearest neighbors. Finally, for the case of a low self-driving force and very long-range attractive interactions the active matter forms a collapsed state. This is a dense state where we have observed quite irregular motion of the particles.
We performed both numerical and theoretical analysis of the system. To make the set of ODEs, that describe the particles motion, less stiff, we use a modified Morse potential, which has the gaussian form of its repulsive part. This allows to simulate a few tens of thousands of equations, which suffices to investigate different macroscopic states of the system. The qualitative and quantitative predictions of our theory are in a good agreement with the simulation results. µC a n 2 l 2 a e −1/(n s l 2 a ) 1/2 , where as previously g is the external force. We consider the stationary distribution for a space uniform system in the absence of external forces. In this case the l.h.s of Eq. (29) vanishes, yielding, For uniform and isotropic system the distribution function depends only on velocity module, that is, � = �(v) . Changing to the polar coordinates for the velocity, v = (v, φ) , we obtain, using the independence of on the angular variable φ: Integration this equation and using the boundary condition at v = 0 , we finally obtain: Derivation of the force distribution. Following the derivation of the Holtsmark distribution given for the gravitational force f (r) ∼ 1/r 2 (see e.g. 59 ), we write for the distribution of the local force F: where is the system volume and N is the number of particles. Now we take a Fourier transform of P(F) to obtain: In the limit → ∞ , N → ∞ , n = N/� = const. the last result takes the form 59 , which defines �(k) . Hence we obtain for the distribution function, (32) � ∼ e bv 2 /D−v 4 /2D , www.nature.com/scientificreports/ To compute �(k) we write the force as f (r) = (r/r)f 1 (r) with f 1 (r) = (C a /l a )re −r/l a . Then one obtains for the two-dimensional system: Integration over φ yields, where J 0 (x) is the zero-order Bessel function. Integration over r in (36) may not be performed analytically, hence we apply approximations. To analyze capturing of particles by the fluctuating force we are interested in the force distribution for large forces, that exceed the characteristic force of the order of ∼ C a /l a . In the Fourier space this would correspond to values of k ∼ 1/F smaller than ∼ (C a /l a ) −1 . In other words, we are interested in the values of k satisfying the relation, kC a /l a < 1 . This motivate us to apply the approximation, J 0 (x) ≃ 1 − 1 4 x 2 with x = (kC a /l a )e −r/l a , which yields, Using Eq. (37) we find the distribution of the fluctuating forces: where a = πC 2 a /8 . It may be easily checked that P(F) is normalized and that �F� = dFP(F)F = 0.