A Quantitative Kinetic Theory of Flocking with Three-Particle-Closure
AA Quantitative Kinetic Theory of Flocking with Three-Particle-Closure
Rüdiger Kürsten and Thomas Ihle
Institut für Physik, Universität Greifswald, Felix-Hausdorff-Str. 6, 17489 Greifswald, Germany (Dated: February 25, 2021)We consider aligning self-propelled particles in two dimensions. Their motion is given by general-ized Langevin equations and includes non-additive N -particle interactions. The qualitative behavioris as for the famous Vicsek model. We develop a kinetic theory of flocking beyond mean field. Inparticular, we self-consistently take into account the full pair correlation function. We find excellentquantitative agreement of the pair correlations with direct agent-based simulations within the dis-ordered regime. Furthermore we use a closure relation to incorporate spatial correlations of threeparticles. In that way we achieve good quantitative agreement of the onset of flocking with directsimulations. Compared to mean field theory, the flocking transition is shifted significantly towardslower noise because directional correlations favor disorder. We compare our theory with a recentlydeveloped Landau- kinetic theory. I. INTRODUCTION
Entities equipped with a propulsion mechanism, that isactive matter, transfer free energy into directed motion.Such entities exist from micro or even nano length scalesup to macroscopic size. Examples are man-made microswimmers, bacteria, insects, fish, mammals and robots.Large groups of interacting active particles can exhibitcomplex emergent collective phenomena that differ fun-damentally from their equilibrium counter parts. Someprominent examples are flocking [1, 2], motility inducedphase separation [3, 4], and bacterial turbulence [5]. Anoverview on the rapidly developing field of active matterand its application can be found e.g. in the followingreviews [6–14].The Vicsek model [1] is historically one of the first andcomputationally one of the simplest active models thatexhibits a flocking transition. Therefore it is consideredas one of the prototype models of active matter. In thismodel, point particles move in two dimensions at con-stant speed v in individual directions given by the polarangles φ i , i ∈ { , . . . , N } . At points in time τ, τ, τ, . . . all particle directions φ i are changed due to interactionsin the following way: Each particle takes the directionof the average velocity of all particles within distance R (including the particle itself), disturbed by a noise term.Thus, the interactions favor a local alignment of the par-ticle velocities.For large systems and periodic boundary conditionsthe Vicsek model is known in four phases [1, 12, 15–18].For small noise or large particle densities, on average allparticles move in a similar direction (i) [1].Hence we call the system polarly ordered. Further-more, particles tend to cluster together locally. However,the clusters are distributed equally over space and thusthe particle distribution is homogeneous on larger lengthscales. By increasing the noise strength or decreasing theparticle density the system arranges in a ’cross sea’ phase(ii) [18], where high density regions that look like cross- ing wave fronts are formed. This pattern moves throughthe system and the low density regions have almost nopolar order. For even larger noise strength or smallerdensities the system arranges in parallel non-intersectingtraveling high density bands (iii) [15]. Again, there isalmost no polar order in the low density regions betweenthe high density bands. In phases (ii) and (iii) there isstill an average polar order. The main contributions tothe polar order come from the high density regions. Forvery large noise strength or very small particle densitiesthere is no polar order (iv). That is, for large systemsthere is no motion of the center of mass and the particlesare distributed homogeneously.There have been qualitative descriptions of some ofthe aforementioned phases by field- or kinetic theories[2, 19–22]. Those theories either do not reach quantita-tive agreement with direct simulations or only for veryspecial parameters. In addition, most theories rely onthe mean field assumption. There are also kinetic theo-ries of active systems that consider weak pair correlationsto some extent, see e.g. [23–26]. However, for the Vic-sek model, recently it has been shown quantitatively, thatcorrelations of two and more particles are important evenin large parts of the disordered phase (iv) [27].The aim of this paper is to provide a more quantita-tive theoretical description of the model. For technicalreasons we are not studying the Vicsek model itself buta similar model that is believed to behave qualitativelyequivalent to the Vicsek model. Major problems in ki-netic theories of the standard Vicsek model are difficultiesrelated to the finite time step as well as the presence ofmulti-particle collision integrals that are not analyticallysolvable, see e.g. [20]. The latter was circumvented in[23] by using binary interactions with randomly selectedinteraction partners.In this paper we develop a ring-kinetic theory. Thatmeans, we consider the dynamical equations for the one-particle distribution and for the two-particle correlationfunction explicitly. Higher order correlations are ne-glected in the first step. However, the effects of the pair a r X i v : . [ c ond - m a t . s o f t ] F e b correlations on the one-particle distribution as well as onthe pair correlation function itself, are fully taken intoaccount. This concept has been applied in various fields,see e.g. [28–31]. In this paper, we set up the ring-kineticequations for N -particle interactions, a significant dif-ference to most previous applications. That means theforces do not depend on the state of only two particlesbut on the positions of all particles.Within the disordered parameter regime, we find ex-cellent quantitative agreement between the ring-kinetictheory and agent-based simulations. For a moderate par-ticle density, we find the threshold noise for the onset offlocking within the ring-kinetic theory to be very closeto the value measured in agent-based simulations. Inparticular, the results are much more precise than thepredictions of mean field theory.For larger particle densities, at the onset of flocking,higher order correlations are more important. We extendthe ring-kinetic theory and incorporate also spatial three-particle correlations via a closure ansatz. In that way,the flocking transition can be described also for largerdensities, and the agreement between theory and simu-lations is further improved at moderate densities, wherewe reach quantitative agreement within the consideredresolution. In addition, for higher densities, our theoryis a significant improvement over mean field.The paper is organized as follows. In Sec. II we de-fine and discuss the Vicsek-like model studied here. InSec. III we develop a homogeneous mean field theory andcalculate the critical noise strength of the flocking transi-tion. In Sec. IV we introduce the notion of many particlecorrelations. In Sec. V we develop the full ring-kinetictheory, that is the dynamical equations for the one par-ticle distribution and the pair correlation function ne-glecting higher order correlations. In Sec. VI we presentextensive quantitative comparisons between the solutionsof the ring-kinetic equations and direct agent-based sim-ulations. We identify the parameter region where thering-kinetic theory is applicable. In Sec. VII we employ aclosure relation to take into account spatial three particlecorrelations. In that way we can significantly enlarge theapplicability domain of the kinetic theory. In Sec. VIII wediscuss our results and give an outlook to possible exten-sions of the method. In Appendix A we give all relevantequations in Fourier space that have been used to eval-uate the kinetic equations numerically. In Appendix Bwe compare the results of the kinetic theory presentedin this work with the simplified Landau kinetic theory ofPatelli [26]. II. MODEL
We consider a Vicsek-like model in continuous timethat was investigated as here, or in a similar form in many studies, see e.g. [32–35]. The model is given by ˙ x i = v cos( φ i )˙ y i = v sin( φ i )˙ φ i = w ( | Ω( i ) | ) (cid:88) j ∈ Ω( i ) sin( φ j − φ i ) + σξ i , i = 1 , . . . , N, (1)where r i ( t ) = ( x i ( t ) , y i ( t )) denote the particles positionsand φ i the directions of the particles velocities. The set Ω i := { j ∈ { , . . . , N } : | r j − r i | ≤ R } (2)contains all particles that are within distance R to par-ticle i . The ξ i ( t ) are independent Gaussian white noiseterms satisfying (cid:104) ξ i ( t ) ξ j ( s ) (cid:105) = δ ij δ ( t − s ) . (3)The noise strength is given by σ , and w ( n ) is an inter-action weight function that depends on the number ofneighbors of particle i (particles that are within distance R including particle i itself).We consider the two-dimensional motion of N particlesthat move at constant speed v in individual directions φ i . The directions of neighboring particles tend to align,however, they are disturbed by noise. In the following,we discuss the cases w ( n ) = 1 = const. and w ( n ) = 1 /n .From a technical point of view, the case w ( n ) = 1 = const. is the most desirable since in this case, the modelincludes only pair interactions. This has the advantagethat, as in the regular BBGKY-hierarchy, three-particlecorrelations can be produced only from previously exist-ing pair correlations, four-particle correlations can onlybe produced from three-particle correlations and so on.In contrast, if w is a function of n as e.g. w ( n ) = 1 /n the interactions are in fact N -particle interactions. Thatmeans, all orders of correlations are produced immedi-ately even if the model evolves from uncorrelated initialconditions.Intuitively, one might think that such tiny details of themodel are not that important and might lead to qualita-tively equivalent results. The average number of neigh-bors is a constant anyway and one might hope that thefluctuations of the number of neighbors do not play amajor role. Surprisingly, that is not the case and the twomodels differ qualitatively, see also [36] for a detailedanalysis. For w ( n ) = 1 = const. one still finds a ho-mogeneous phase at large noise and a polarly orderedphase at smaller noise for systems of finite size. How-ever, in the case of polar order, particles do not arrangein high density bands or cross sea patterns but they formhigh density clusters that contain almost all particles, seeFig. 1.It is not even clear if there is a disordered phase at all inthe thermodynamic limit. An alternative hypothesis is as σ . . . . . . | p | (a) (b)(c) (d) FIG. 1. Numerical results for the system (1) with w ( n ) = 1 = const. , simulated with an Euler-Maruyama-scheme [38] with ∆ t = 0 . . Parameters are N = 10 , L ≈ . , v = R = 1 ,thus C := ρ πR = 1 . (a) Polar order parameter as a func-tion of noise strength. (b-d) Snapshots after a thermalizationtime T = 1000 in the disordered phase (b) σ = 2 . , and inthe polarly ordered phase (c) σ = 1 . and (d) σ = 0 . . follows. For every finite but possibly large noise strengthone finds polar ordered clusters of very high densities forlarge enough systems. Due to the additive nature of thealignment interactions, such a cluster can remain stableeven for large noise strength if the density is large enough.It could be that such a clustered state is the steady statefor large systems at any noise strength. However, so-phisticated investigations are necessary to answer thisquestion. Differences between additive and non-additiveinteractions have been noticed already in [37], they havebeen studied in more detail in a very recent work [36].Although the model with additive interactions is very in-teresting as well, we are mainly interested in the studyof models behaving qualitatively like the Vicsek model.Hence we focus on the case w ( n ) = 1 /n (4)in this paper. III. KINETIC THEORY
Instead of investigating the Langevin equations,Eq. (1), we can equivalently study the associated N - particle Fokker-Planck equation [39] ∂ t P N (1 , , . . . , N, t ) = − N (cid:88) i =1 ∂ φ i [ N (cid:88) j =1 w ( | Ω( i ) | ) θ ij × sin( φ j − φ i ) P N (1 , , . . . , N, t )] − N (cid:88) i =1 v cos( φ i ) ∂ x i P N (1 , , . . . , N, t ) − N (cid:88) i =1 v sin( φ i ) ∂ y i P N (1 , , . . . , N, t )+ N (cid:88) i =1 σ ∂ φ i P N (1 , , . . . , N, t ) , (5) θ ij := θ ( R − (cid:113) ( x j − x i ) + ( y j − y i ) ) , (6)where we used the abbreviation of writing just insteadof x , y , φ and similar for , , . . . . The probability den-sity of finding the system in a phase space volume aroundthe coordinates given by , , . . . , N at time t is denotedby P N (1 , , . . . , N, t ) . This equation is exact.We assume that the probability density is initially andhence at all times symmetric with respect to the permu-tation of coordinates i ↔ j . We denote the marginalizedprobability density by P k (1 , , . . . , k, t ) := (cid:90) P N (1 , , . . . , N, t )d( k + 1) . . . d N, (7)where d l denotes d x l d y l d φ l and integration is performedover the interval [0 , L ) for x l and y l and over the interval [0 , π ) for φ l .Using the Fokker-Planck equation (5) and the symme-try of permutations of coordinates we obtain the timeevolution of the marginalized density ∂ t P (1 , t ) = − ∂ φ [( N − (cid:90) w ( | Ω(1) | )d2 . . . d N θ × sin( φ − φ ) P N (1 , , . . . , N, t )] − v cos( φ ) ∂ x P (1 , t ) − v sin( φ ) ∂ y P (1 , t )+ σ ∂ φ P (1 , t ) . (8)We observe that the time evolution of the one particleprobability depends on the complete N particle proba-bility due to the fact that the number of neighbors | Ω(1) | of the first particle depends on the positions of all theother particles.For spatially homogeneous solutions, P depends onlyon φ and t , hence we can define p φ ( φ , t ) := L P (1 , t ) (9)such that Eq. (8) simplifies to ∂ t p φ ( φ , t ) = − ∂ φ [( N − (cid:90) w ( | Ω(1) | )d1 . . . d N θ × sin( φ − φ ) P N (1 , , . . . , N, t )]+ σ ∂ φ p φ ( φ , t ) , (10)Considering the noise strength σ as control parameter,in the thermodynamic limit N → ∞ , NL = const. , thesystem exhibits a transition from disorder to polar order.In the disordered phase, all particles directions are dis-tributed randomly. In the polar ordered phase there is aglobal preferred direction of motion.In the next section we calculate the threshold noisestrength where the transition takes place, within themean field theory. A. Mean Field Theory
Solving Eq. (10) exactly is very difficult or even impos-sible. Therefore we are looking for the stationary solutionunder the mean field assumption P N (1 , , . . . , N ) = P (1) P (2) . . . P ( N ) , (11)which simplifies with Eq. (9) to P N (1 , , . . . , N ) = 1 L N p φ ( φ ) p φ ( φ ) . . . p φ ( φ N ) (12)In order to simplify Eq. (10) it is reasonable to split thespatial integration domain into pieces such that the num-ber of neighbors of particle one is constant on each piece.Inserting this assumption into Eq. (10) we obtain ∂ t p φ ( φ , t ) = − ∂ φ (cid:20) N − L N − N − (cid:88) n =1 (cid:90) w ( n + 1)d2d r . . . d r N θ . . . θ n +1) × (1 − θ n +2) ) . . . (1 − θ N ) sin( φ − φ ) × p φ ( φ , t ) p φ ( φ , t ) + permutations (cid:21) + σ ∂ φ p φ ( φ , t ) , (13)where r i = ( x i , y i ) . In the above integral r , . . . , r n +1 are integrated over the inside of the disk with radius R around particle one and r n +2 , . . . , r N are integrated overthe outside of the disk around particle one. That means n is the number of neighbors of particle one. By permu-tations we mean all other choices of n − particles from { , . . . , N } to be integrated over the inside. All this per-mutations give the same contribution, hence we can takecare of them by a combinatorial factor to arrive at ∂ t p φ ( φ , t ) = − ∂ φ (cid:20) N − L N − N − (cid:88) n =1 (cid:18) N − n − (cid:19) (cid:90) w ( n + 1)d2d r . . . d r N × θ . . . θ n +1) (1 − θ n +2) ) . . . (1 − θ N ) sin( φ − φ ) × p φ ( φ , t ) p φ ( φ , t ) (cid:21) + σ ∂ φ p φ ( φ , t ) . (14)Here, in the two dimensional case, the C coefficient iscalculated as C := NL (cid:90) θ d r , (15)that is the expectation value of the number of particleswithin a circle of radius R . Since we are interested in thethermodynamic limit N → ∞ it suffices to approximate ( N − (cid:18) N − n − (cid:19) ≈ N n ( n − . (16)Inserting Eqs. (15) and (16) into Eq. (14) we obtain ∂ t p φ ( φ , t ) = − N − (cid:88) n =1 C n ( n − (cid:0) − C N (cid:1) N − − n w ( n + 1) × ∂ φ (cid:20) (cid:90) d φ sin( φ − φ ) p φ ( φ , t ) p φ ( φ , t ) (cid:21) + σ ∂ φ p φ ( φ , t ) . (17)Taking the limit N → ∞ and substituting k = n − wearrive at ∂ t p φ ( φ , t ) = − C ∞ (cid:88) k =0 C k k ! exp( − C ) w ( k + 2) × ∂ φ (cid:20) (cid:90) d φ sin( φ − φ ) p φ ( φ , t ) p φ ( φ , t ) (cid:21) + σ ∂ φ p φ ( φ , t ) . (18)One way of treating Eq. (18) is to study its full Fouriertransform. Since here, we are only interested in the crit-ical noise strength it suffices to consider only the zerothand the first order in the Fourier transform. Rotatingthe coordinate system in an appropriate way we can usethe ansatz p φ ( φ ) = 12 π + ε cos( φ ) . (19)The isotropic distribution corresponding to ε = 0 is al-ways a solution of Eq. (18). At the critical noise strengthit changes stability.Inserting the ansatz (19) into Eq. (18), multiplying by cos( φ ) /π integrating over φ and neglecting terms oforder ε we obtain by partial integration ∂ t ε = − π C ∞ (cid:88) k =0 C k k ! exp( − C ) w ( k + 2) × (cid:26) (cid:90) π d φ d φ sin( φ ) sin( φ − φ ) × (cid:20) π + ε cos( φ )2 π + ε cos( φ )2 π (cid:21)(cid:27) − σ ε = (cid:20) C ∞ (cid:88) k =0 C k k ! exp( − C ) w ( k + 2) − σ (cid:21) ε. (20)Hence the mean field critical noise strength is σ c = (cid:118)(cid:117)(cid:117)(cid:116) C ∞ (cid:88) k =0 C k k ! exp( − C ) w ( k + 2) . (21)In case of the constant weight function w ( n ) = 1 thisreduces to σ c = (cid:112) C . (22)For a Poisson distribution we have (cid:28) k + 1 (cid:29) = ∞ (cid:88) k =0 k + 1 C k k ! exp( − C )= 1 C ∞ (cid:88) k =0 C k +11 ( k + 1)! exp( − C )= 1 C ∞ (cid:88) ˜ k =1 C ˜ k ˜ k ! exp( − C )= 1 C (cid:20) ∞ (cid:88) ˜ k =0 C ˜ k ˜ k ! exp( − C ) − exp( − C ) (cid:21) = 1 C [1 − exp( − C )] (23)and similar (cid:28) k + 2 (cid:29) = (cid:28) k + 1 (cid:29) − (cid:28) k + 1)( k + 2) (cid:29) = (cid:28) k + 1 (cid:29) − C ∞ (cid:88) ˜ k =2 C ˜ k ˜ k ! exp( − C )= (cid:28) k + 1 (cid:29) − C [1 − (1 + C ) exp( − C )] . (24)Hence, for the weight function w ( n ) = 1 /n we obtain thecritical noise strength using Eqs. (21), (23) and (24) σ c = (cid:114) { − exp( − C ) − C [1 − (1 + C ) exp( − C )] } = (cid:114) − C [1 − exp( − C )] . (25) IV. CORRELATIONS
In general, the mean field assumption (11) is not validbut must be corrected by terms containing correlations ofvarious order. The correlation functions G k are definedrecursively by [28, 40, 41] G ≡ P , (26) G k (1 , · · · , k ) := P k (1 , · · · , k ) (27) − (cid:88) σ k − (cid:88) l =1 l − k − l )! G l (1 , σ (2) , · · · , σ ( l )) × P k − l ( σ ( l + 1) , · · · , σ ( k )) , where (cid:80) σ denotes the sum over all permutations ofthe elements { , · · · , k } . We can rewrite Eq. (27) as P l ( . . . ) = G l ( . . . ) + . . . and insert it recursively for all P l of order l < k into Eq. (27) and eventually replace P by G . Performing such an expansion, only P k and G -functions remain on the right hand side of Eq. (27). Itfollows inductively from Eq. (27) that the indexes of allcorrelation functions on the right hand side are ordered,that is for each term G l ( i , . . . , i l ) appearing in the ex-pansion of Eq. (27) it holds i < i < · · · < i l . Thus,instead of Eq. (27) we might alternatively write G k (1 , . . . , k ) = P k (1 , . . . , k ) − (cid:88) { over all possibleproducts of G -functions such that each of the arguments , . . . , k appears exactly once and for each G -functionthe arguments are ordered. } (28)For example, the two-, three- and four- particle correla-tion functions are given explicitly as G (1 ,
2) = P (1 , − P (1) P (2) (29) G (1 , ,
3) = P (1 , , − P (1) P (2) P (3) − P (1) G (2 , − P (2) G (1 , − P (3) G (1 , , (30) G (1 , , ,
4) = P (1 , , , − P (1) P (2) P (3) P (4) − G (1 , P (3) P (4) − G (1 , P (2) P (4) − G (1 , P (2) P (3) − G (2 , P (1) P (4) − G (2 , P (1) P (3) − G (3 , P (1) P (2) − G (1 , G (3 , − G (1 , G (2 , − G (1 , G (2 , − G (1 , , P (4) − G (1 , , P (3) − G (1 , , P (2) − G (2 , , P (1) . (31)We want to mention some important properties of thecorrelation functions G k . Since we assume that the N -particle probability distribution is symmetric with re-spect to particle exchanges, also the correlation functionsfollow the same symmetry. Furthermore, for k ≥ itholds (cid:90) G k (1 , . . . , k )d k = 0 , (32)where the integration is performed over the whole spaceand all possible orientations of particle k . This prop-erty follows by induction over k from Eq. (27). It showsthat G k contains only information about k -particle cor-relations and not about lower order correlations. As asimple consequence it follows that (cid:90) X G k (1 , . . . , k )d k = − (cid:90) X C G k (1 , . . . , k )d k, (33)where X is some subset of the set of all possible con-figurations of particle k and X C is the complement of X . V. RING-KINETIC THEORY
The mean field assumption (11) is equivalent to vanish-ing correlation functions G i ≡ for i ≥ . In this paperwe go one step further and take two particle correlationsinto account. However, we still assume that higher ordercorrelations can be neglected. That is, we consider P (1 , t ) ,G (1 , , t ) := P (1 , , t ) − P (1 , t ) P (2 , t ) ,G i (1 , · · · , i, t ) ≡ for i ≥ . (34)Similar assumptions have been made in Ref. [23] butalso in a very recent kinetic Landau theory [26] on themodel with additive interactions, that is with w ( n ) ≡ .Therein, it was assumed additionally that G is smalland furthermore that also noise and coupling are small.It should be mentioned that according to [27] and to theresults presented in this work, it follows that in the vicin-ity of the flocking transition and even in a considerablepart of the disordered phase, the assumption of negligible G and small G is not justifiable. However, far in thedisordered regime it is a reasonable approximation. InAppendix B we present a detailed comparison of the Lan-dau theory of [26] applied to the model of non-additiveinteractions with our full ring-kinetic theory.In this paper we impose no restrictions on G , thatmeans we allow large pair correlations. As we assumespatial homogeneity, the one-particle probability density P is independent on the position ( x , y ) and it sufficesto consider the angular dependence p φ given by Eq. (9).Similarly, the two-particle probability density P dependsonly on the angles φ , φ and on the distance betweenthe spatial coordinates ∆ x := x − x , ∆ y := y − y , ∆ = (∆ x , ∆ y ) . Hence it is reasonable to introduce thereduced probability density p ( φ , φ , ∆) := L (cid:90) P ( φ , φ , x , x , y , y ) × δ (∆ x − ( x − x )) δ (∆ y − ( y − y ))d x d x d y d y = L P ( φ , φ , x , x = x + ∆ x , y , y = y + ∆ y ) (35) and correlation functions g ( φ , φ , ∆) := L (cid:90) G ( φ , φ , x , x , y , y ) × δ (∆ x − ( x − x )) δ (∆ y − ( y − y ))d x d x d y d y = L G ( φ , φ , x , x = x + ∆ x , y , y = y + ∆ y ) . (36)We obtain the time evolution equations of the reducedone- and two-particle probability density functions, thatare the first two equations of the BBGKY-hierarchy, fromthe Fokker-Planck equation (5) using the marginalization(7), plugging in the cluster expansion (27) and using thering-kinetic ansatz (34) as well as properties (32) and(33) of the correlation functions. We use the followingdiagramatic representation of the appearing collision in-tegrals1 2 = sin( φ − φ ) θ ( R − | r − r | ) , (37)1 2 = θ ( R − | r − r | ) , (38)1 2 = 1 − θ ( R − | r − r | ) , (39)1 2 = g (1 , , (40)k = p ( φ k ) , if k is not connected (41)via a coiled line,1 = ∂ φ . (42)All terms representing a symbol in a diagram are meantto be multiplied and integrated over all degrees of free-dom that do not occur on the left hand side of the equa-tion in which the diagram appears. As an example, wegive the meaning of the following integral1 2 = (cid:90) ∂ φ g ( φ , φ , ∆) sin( φ − φ ) × θ ( R − | ∆ | )d φ d∆ . (43)For the one particle angular distribution we find withinour diagramatic notation ∂ t p ( φ ) = − ρw (cid:20) + (cid:21) − ρ ( w − w ) (cid:20) + (cid:21) − ρ ( w − w + w ) + σ ∂ φ p ( φ ) (44)Here, w k denotes the expectation value of the weightfunction for k + l neighbors with respect to the numberof neighbor distribution q n ( l ) , that is w k := (cid:104) w ( k + l ) (cid:105) = ∞ (cid:88) l =0 w ( k + l ) q n ( l ) . (45)The number of neighbor distribution for correlated par-ticles was recently derived in [27]. In case that only two-particle correlations are present this distribution dependson the two parameters C := N (cid:90) G (1 , θ θ d1d2 ,D := N (cid:90) G (1 , θ d1d2 , (46)where θ i := θ ( R − (cid:113) x i + y i ) (47)and θ ij is given by Eq. (6) with the Heaviside function θ .The number of neighbor distribution is given in [27] as q n ( l ) = q ( l ) − [ q ( l − − q ( l )]( D − C , (48)where q ( l ) is the following infinite series q ( l ) =( C − C ) l exp( C / − C ) × ∞ (cid:88) k =0 (cid:20) C C − C ) (cid:21) k k !( l − k )! . (49)As before, C = NL πR is just the average number ofneighbors. Alternatively to the representation by an in-finite sum (49) one can give a simple expression for thecharacteristic function of q ( l ) [27] χ ( u ) = exp (cid:20) (cid:88) l =1 l (cid:88) t =0 ( − l + t C l l ! (cid:18) lt (cid:19) exp( itu ) (cid:21) . (50) Coming back to Eq. (44), we give a guide to explainall interaction integrals. Those are all possible diagramswhere each particle is connected to via a solid line oran arrow (that means all particles are neighbors of par-ticle ) such that every point is connected to also viaa path consisting only of arrows and/or coils and thereis exactly one arrow that starts from . Particles thatare not connected to via a path of arrows and/or coilsdo not directly take part in the interactions and thusdo not appear in the diagrams. Because there is onlyone arrow in the diagrams and each point can be con-nected to only one coil (that represents g ) the diagramscan have no more than four points. Each diagram hasa prefactor of − ρ k − , where k is the number of particlesinvolved in the diagram. This prefactor is a combinationof a combinatorial factor and /L k from Eqs. (9) or (36).Furthermore, there is another prefactor of w k that takesinto account the expectation value of the weight functionwhen an average over all particles not involved in thediagram is performed. In principle, we would have addi-tional diagrams where points not connected via an arroware connected to via dashed lines, that means thoseparticles are not a neighbor of particle . However, wecan replace each dashed line by a solid line and multiplywith − to compensate it due to property (33). However,for these additional diagrams we have different prefactorsof w k − or w k − if one or two dashed lines are involved,respectively, due to the different number of neighbors ofparticle . Considering all the aforementioned interac-tion terms, angular diffusion and streaming we arrive atEq. (44) Note, that for pair interactions, w ( n ) = const. ,we have w = w = w and only the first two interactionintegrals remain.Analogously to Eq. (44), we find starting from Eq. (5),the time evolution equation for the two particle proba-bility distribution ∂ t p ( φ , φ , ∆) = − w (cid:20) + (cid:21) − ρw (cid:20) + + + (cid:21) − ρ ( w − w ) (cid:20) + + + + + (cid:21) − ρ ( w − w + w ) (cid:20) + + (cid:21) − ρ ( w − w + 3 w − w ) (cid:20) (cid:21) − ρw (cid:20) + + + (cid:21) − ρ ( w − w ) (cid:20) + + + + + (cid:21) − ρ ( w − w + w ) (cid:20) + + (cid:21) − ρ ( w − w + 3 w − w ) (cid:20) (cid:21) + σ ∂ φ [ p ( φ ) p ( φ ) + g ( φ , φ , ∆)] + 1 ↔ v (cos φ − cos φ ) ∂ ∆ x g ( φ , φ , ∆) + v (sin φ − sin φ ) ∂ ∆ y g ( φ , φ , ∆) , (51)where ↔ is an abbreviation for all previous termswith particles one and two interchanged.Employing Eq. (29), in our notation (9), (35) and (36), g ( φ , φ , ∆) = p ( φ , φ , ∆) − L p φ ( φ ) p φ ( φ ) , (52) we obtain the time evolution equation of the pair corre-lation function ∂ t g ( φ , φ , ∆) = ∂ t p ( φ , φ , ∆) − L p φ ( φ ) ∂ t p φ ( φ ) − L p φ ( φ ) ∂ t p φ ( φ ) . (53)Inserting Eqs. (44) and (51) yields ∂ t g ( φ , φ , ∆) = − w (cid:20) + (cid:21) − ρw (cid:20) + (cid:21) − ρ ( w − w ) (cid:20) + (cid:21) − ρ ( w − w ) (cid:20) + + + (cid:21) − ρ ( w − w + w ) (cid:20) + (cid:21) − ρ ( w − w + w ) (cid:20) + (cid:21) − ρ ( w − w + 3 w − w ) (cid:20) (cid:21) − ρ ( w − w + 3 w − w ) (cid:20) (cid:21) (54) − ρw (cid:20) + (cid:21) − ρ ( w − w ) (cid:20) + + + (cid:21) − ρ ( w − w + w ) (cid:20) + (cid:21) − ρ ( w − w + 3 w − w ) (cid:20) (cid:21) + σ ∂ φ g ( φ , φ , ∆) + 1 ↔ v (cos φ − cos φ ) ∂ ∆ x g ( φ , φ , ∆) + v (sin φ − sin φ ) ∂ ∆ y g ( φ , φ , ∆) . (55)Note that, for example, the multiplication [ ∂ t p φ ( φ )] × p φ ( φ ) =[ ∂ t p φ ( φ )] + [ ∂ t p φ ( φ )] ∂ t p and − L p φ ( φ ) ∂ t p φ ( φ ) or − L p φ ( φ ) ∂ t p φ ( φ ) cancel such as for example 1 23 . However, theircounter parts with particles one and two being neighbors,such as for example 1 23 do not cancel. This isan effect of the N − particle interactions. As a result theseterms appear with different sign but also with differentprefactors and thus they do not cancel. This effect isresponsible for a number of terms that are not presentfor pair interactions, w ( n ) = const. . VI. COMPARISON OF RING-KINETICTHEORY AND AGENT-BASED SIMULATIONS
Recently, a quantitative numerical method to predictthe minimal required order of correlations has been in-troduced in [27]. It is based on the measurement of thenumber of neighbors distribution (neighbors are particlesthat are closer than R ) of a randomly selected particle.Besides this measurement, the number of neighbors dis-tribution is also calculated under the assumption thatonly correlations up to a given order l max are present,that is G l ≡ for all l > l max . If both distributions agreereasonably well, the correlation order l max is consideredto be sufficient, see [27] for details.In Fig. 2 we show a correlation map that we obtainedby this method for large systems of N = 10 particles. Itshows, depending on parameters, which order of correla-tions is required. We expect excellent quantitative pre-dictions of the ring-kinetic theory in a parameter regimewhere pair correlations are sufficient. Strictly speaking,this is the case in the disordered phase only. However,0 C = ρ πR σ c o rr e l a t i o n o r d e r none ≥ FIG. 2. Correlation map of the Vicsek-like model given byEq. (1) with w ( n ) = 1 /n . Colors encode the necessary or-der of correlations (e.g. cyan ∼ ∼ three particle correlations,etc.) for a quantitative description of the system. The corre-lation map was obtained by a recently introduced method [27]comparing the measured number of neighbor distribution toa theoretically predicted one. One considers the distributionsto be acceptably close if the corresponding Kullback-Leiblerdivergence is less than − . We simulated realizationsfor each noise strength σ ∈ { . , . , . , . . . , . } and eachparticle density C = ρπR ∈ { . , . , . , . . . , . } with pa-rameters v = R = 1 , N = 10 . The correlation parametershave been measured for a time interval of after a thermal-ization time of . Eq. (1) was integrated using an Euler-Maruyama scheme with time step ∆ t = 10 − . The analysishas been done as in [27]. The red curve shows the mean fieldtransition line given by Eq. (25). we find that a reasonable ring-kinetic description of thesystem is still possible if the influence of higher order cor-relations is already measurable but still not dominant.We solve the time evolution equations for p φ and g ,Eqs. (44) and (55), numerically in Fourier space, see Ap-pendix A. As a measurable quantity we consider the stan-dard radial distribution function g ( r ) that is defined by g ( r ) = L N ( N − (cid:88) i (cid:54) = j πr (cid:104) δ ( r − | r i − r j | ) (cid:105) . (57)It is related to the correlation function g given inEq. (36) by g ( r ) = 1+ (cid:90) π d φ (cid:90) π d φ (cid:90) L d∆ x (cid:90) L d∆ y × g ( φ , φ , ∆) 12 πr δ ( r − | ∆ | ) (58)and can be calculated analytically from the Fourier rep-resentation, see Appendix A.In Fig. 3 we compare the radial distribution func-tion measured in agent-based simulations with the re-sults from ring-kinetic theory. Overall they agree very r . . . . . g ( r ) FIG. 3. Radial distribution function g ( r ) obtained from ring-kinetic theory (blue solid line) and sampled directly fromagent-based simulations according to Eq. (57) (black solidline). For the ring-kinetic theory we integrate the time evolu-tion equations (44) and (55) in Fourier space (see AppendixA) with an Euler scheme with step size ∆ t = 10 − for anabsolute time of , starting initially uncorrelated and dis-ordered. For the agent based simulation we integrate Eq. (1)with an Euler-Maruyama scheme with step size ∆ t = 10 − .After a thermalization time of we average over a timeof and over realizations. Parameters are σ = 1 . , v = R = 1 , N = 509 , C = ρπR = 1 , L ≈ . Thering-kinetic theory uses Fourier modes F klmn according toEq. (A2) with angular indexes k, l ∈ {− , . . . , } and spatialindexes m, n ∈ {− , . . . , } . Performing a spatial Fouriertransform of the measured curve (black line) with indexesfrom the same range m, n ∈ {− , . . . , } and transformingback into real space we obtain the dashed green line. Thedashed solid line at y = 1 serves as a guide to the eyes. well, however, there are minimal deviations. Those de-viations can be explained by the finite resolution usedin the Fourier transform. We evaluated the ring-kineticequations in Fourier space using a minimal spatial wavelength of . . Since the radial distribution functions g ( r ) drops rapidly at about r ≈ it is not perfectly re-solved, causing small deviations. To test this hypothesis,we Fourier transform the pair correlation function cor-responding to the measured g ( r ) with the same spatialresolution and Fourier transform it back into real space.The resulting curve deviates by less than one percentfrom the ring-kinetic result, see Fig. 3.The ring-kinetic theory predicts the full two particlecorrelation function g ( φ , φ , ∆) and not only the radialdistribution function g ( r ) that is an integral of it, seeEq. (58). Without spontaneous symmetry breaking, thatis in the disordered phase, the system is isotropic. Inthat case, the pair correlation function depends only onthree independent arguments and not on four, such as φ , φ , ∆ x and ∆ y . We choose those three degrees offreedom as the length of the vector ∆ , r := | ∆ | , the1
11 22 11
11 22
FIG. 4. Sketch of the three arguments r , α and ∆ φ of the paircorrelation function h for homogeneous and isotropic systems,see Eq. (59) for the definition of h ( r, α, ∆ φ ) . The blue pointsdenote the positions of particles one and two. Their distanceis denoted by r and their direction of motion (red arrows)is given by the angles φ and φ . The difference betweenthose two angles is denoted by ∆ φ . The difference betweenthe polar angle γ of the vector pointing from particle one toparticle two and φ is denoted by α . difference between the orientations of the two particles, ∆ φ := φ − φ , and the difference of the polar angleof the vector ∆ and the orientation of the first particle α := γ − φ , where γ is the polar angle of the vector ∆ ,see Fig. 4. Depending one those arguments we define thefunction h ( r, α, ∆ φ ) :=2 π (cid:90) π d φ (cid:90) π d φ (cid:90) π d γ × g ( φ , φ , ∆ x = r cos γ, ∆ y = r sin γ ) × δ ( α − γ + φ ) δ (∆ φ − φ + φ ) . (59)For translational and rotational invariant systems it con-tains all information about two particle correlations. Itcan be directly sampled from numerical data accordingto h ( r, α, ∆ φ ) = L (2 π ) N ( N −
1) 12 πr (cid:88) i (cid:54) = j (cid:104) δ ( r − | r i − r j | ) × δ ( α − γ ij + φ i ) δ (∆ φ − φ j + φ i ) (cid:105) − . (60) A. Dependence of h ( r, α, ∆ φ ) on ∆ φ We consider the dependence of the correlation function h ( r, α, ∆ φ ) on the difference in orientation of the twoparticles ∆ φ for fixed r and α . In Fig. 5 we fix r = 0 . and α = 0 , π/ , π, π/ . That means that the two particlesare in the interior of each others interaction circle. We seethat the correlations are highest for ∆ φ = 0 showing that nearby particles align. For ∆ φ = π the correlations areeven slightly negative. Rotations of one particle aroundthe other seem to be of no particular importance at theconsidered distance as we see almost the same picture fordifferent values of α . . . α ∗ = 0 α ∗ = π/ φ . . α ∗ = π − π π φα ∗ = 3 π/ − π π h ( r = . , α = α ∗ , ∆ φ ) FIG. 5. Dependence of the correlation function h ( r, α, ∆ φ ) ondifference of velocity directions ∆ φ for fixed values of r = 0 . and α (value given in the plot). Results of the ring-kinetictheory (blue line) are compared to direct measurements ofagent-based simulations (black line). Parameters as in Fig. 3. In Fig. 6 we fix r = 1 . . That means both particlesare exactly at the boundary of each others interactionregion. For α = 0 and α = π , that means if particle twois in the front or in the back of particle one (looking inthe direction of motion of particle one), the distributionof h (∆ φ ) is still symmetric with a maximum at ∆ φ = 0 .If particle two is on the left or on the right of particleone it is slightly more likely that particle two points abit outward the interaction region than to be perfectlyaligned with particle one. This is reasonable due to thefollowing argument. When particle two is placed at dis-tance r = 1 at a given time, it is more likely that it wasat r < than that it was at r > shortly before, becausethe pair correlations are much higher for r < .In Fig. 7 we fix r = 1 . . Thus, the two particles arenot directly interacting with each other. We see thatthe strength of the correlations is much smaller in thiscase. Also the fluctuations in the measurements of agent-based simulations are larger, because there are less eventssampled. Quantitatively, the correlations are similar tothe case of r = 1 . . B. Dependence of h ( r, α, ∆ φ ) on α In Fig. 8 we consider the α -dependence of the corre-lation function h for r = 0 . and different values of ∆ φ .2 . . α ∗ = 0 α ∗ = π/ φ . . α ∗ = π − π π φα ∗ =3 π/ − π π h ( r = . , α = α ∗ , ∆ φ ) FIG. 6. Two particle correlation function h as a function ofthe difference of the velocity directions of both particles asin Fig. 5, but here for r = 1 . The black bullets sketch therelative positions of particle one (in the center of the circle)and particle two (on the circumference of the circle). Thered arrows indicate the direction of the particle velocities at ∆ φ = 0 . . . α ∗ = 0 α ∗ = π/ φ . . α ∗ = π − π π φα ∗ = 3 π/ − π π h ( r = . , α = α ∗ , ∆ φ ) FIG. 7. Two particle correlation function h as a function ofthe difference of the velocity directions of both particles as inFig. 5, but here for r = 1 . . As discussed in the previous subsection, there is almostno α -dependence. The value of h is largest for ∆ φ = 0 and smallest and even negative for ∆ φ = π which is inagreement with Fig. 5. . . φ ∗ = 0 φ ∗ = π/ α . . φ ∗ = π − π π αφ ∗ = 3 π/ − π π h ( r = . , α , ∆ φ = φ ∗ ) FIG. 8. Dependence of the correlation function h ( r, α, ∆ φ ) on the angle α giving the rotation of particle two around par-ticle with respect to the direction of motion of particle one.Fixed values of r = 0 . and ∆ φ (value given in the plot) areconsidered. Results of the ring-kinetic theory (blue line) arecompared to direct measurements of agent-based simulations(black line). Parameters as in Fig. 3. In Fig. 9 we fix r = 1 . For aligned particles, ∆ φ = 0 ,the α -dependence of the correlation function is still verysmall with a slight preference of particle two being rightor left of particle one compared to a placement in thefront or back of particle one. For anti-aligned particles, ∆ φ = π , it is unlikely that particle two is in the backof particle one. In that case, the particles would moveaway from each other. That means that they have beencloser to each other in the past which makes it unlikelythat they are anti-aligned. For ∆ = π/ and ∆ φ =3 π/ , there is a complex α -dependence with maximalcorrelations at α ≈ . π and α ≈ − . π , respectively.The value of α with minimal correlations is shifted by π with respect to the values of maximal correlations.In Fig. 10 we fix r = 1 . . The correlation functionbehaves qualitatively similar to the case of r = 1 , howeverhere, the correlations are much smaller. C. Dependence of h ( r, α, ∆ φ ) on r In this subsection we consider the r -dependence of thecorrelation function h for fixed values of α and ∆ φ . Wefind the same qualitative behavior for all considered val-ues of α , see Figs. 11-14. For ∆ φ = 0 , π/ and π/ , forsmall r the function h shows a plateau of high correlationsand decreases relatively fast at about r = 1 towards zero.For ∆ φ = π , the plateau for small r is negative and thedecay towards zero at about r = 1 is still present. Thisshows that nearby anti-aligned particles are unlikely.3 . . φ ∗ = 0 φ ∗ = π/ α . . φ ∗ = π − π π αφ ∗ = 3 π/ − π π h ( r = . , α , ∆ φ = φ ∗ ) FIG. 9. Two particle correlation function h as a function ofthe relative orientation α as in Fig. 8, but here for r = 1 . Theblack bullets sketch the relative positions of particle one (inthe center of the circle) and particle two (on the circumferenceof the circle) at α = 0 . The red arrows indicate the directionof the particle velocities. . . φ ∗ = 0 φ ∗ = π/ α . . φ ∗ = π − π π αφ ∗ = 3 π/ − π π h ( r = . , α , ∆ φ = φ ∗ ) FIG. 10. Two particle correlation function h as a function ofthe relative orientation α as in Fig. 8, but here for r = 1 . . . . φ ∗ = 0 φ ∗ = π/
20 1 2 3 r . . φ ∗ = π rφ ∗ = 3 π/ h ( r , α = , ∆ φ = φ ∗ ) FIG. 11. Dependence of the correlation function h ( r, α, ∆ φ ) on the distance r for fixed values of α = 0 and ∆ φ (value givenin the plot). Results of the ring-kinetic theory (blue line) arecompared to direct measurements of agent-based simulations(black line). Parameters as in Fig. 3. . . φ ∗ = 0 φ ∗ = π/
20 1 2 3 r . . φ ∗ = π rφ ∗ = 3 π/ h ( r , α = π / , ∆ φ = φ ∗ ) FIG. 12. Two particle correlation function h as a function ofdistance r as in Fig. 11, but here for α = π/ . . . φ ∗ = 0 φ ∗ = π/
20 1 2 3 r . . φ ∗ = π rφ ∗ = 3 π/ h ( r , α = π , ∆ φ = φ ∗ ) FIG. 13. Two particle correlation function h as a function ofdistance r as in Fig. 11, but here for α = π . . . φ ∗ = 0 φ ∗ = π/
20 1 2 3 r . . φ ∗ = π rφ ∗ = 3 π/ h ( r , α = π / , ∆ φ = φ ∗ ) FIG. 14. Two particle correlation function h as a function ofdistance r as in Fig. 11, but here for α = 3 π/ . D. Applicability of the Ring-Kinetic Theory
In the previous subsections we have seen that forthe considered parameters, the ring-kinetic theory agreesvery well with direct agent-based simulations. There areonly very small deviations partially caused by the finiteresolution in the Fourier transform in the numerical im-plementation of the ring-kinetic equations. In this sub-section we study the applicability of the ring-kinetic the-ory depending on parameters. We have seen in the pre-vious subsections that the spatial pair correlations decayrather rapidly for r (cid:46) R . Therefore, we compare the in-tegrated spatial pair correlations C and D defined inEq. (46) between ring-kinetic theory and direct simula-tions in Fig. 15. We see very good agreement between ring-kinetic theory and simulations for noise strengths σ (cid:38) . . For smaller noise strengths we find seriousdeviations. . . . . σ . . . . . . C (a) . . . . σ . . . . . . D (b) FIG. 15. Local spatial pair correlation parameters C (a)and D (b) as defined in Eq. (46) compared for direct sim-ulations (black line), ring-kinetic theory (blue line), and akinetic theory including a g -closure (red line). The dashedvertical lines shows the onset of flocking measured in agent-based simulations at σ c,sim = 0 . (black), in the ring-kinetic theory at σ c,rk = 0 . (blue), ring-kinetic theorywith closure at σ c,closure = 0 . (red) and in mean fieldtheory at σ c,mf = 0 . (purple). System parameters are: C = ρπR = 1 , R = v = 1 , N = 509 . For the considered parameter set, homogeneous meanfield theory, Eq. (25), predicts the onset of collectivemotion at σ c ≈ . . In the direct simulations wedetermined the onset of flocking at σ c, sim = 0 . from fluctuations of the polar order parameter p = | p | = (cid:113) p x + p y , with p x = 1 /N (cid:80) Ni =1 cos φ i , p y =1 /N (cid:80) Ni =1 sin φ i , see Fig. 16.Returning to the deviations of the spatial correlationspredicted by ring-kinetic theory in Fig. 15, we find thatthey are significant already before the onset of flock-5 . . . . σ . . . . . . . v a r i a n ce ( p ) FIG. 16. Fluctuations of the polar order parameter p , mea-sured in agent-based simulations. The maximum at σ = 0 . indicates the onset of collective motion. Parameters are as inFig. 15 ing. Nevertheless, the predicted correlations are not com-pletely wrong, but have the correct order of magnitude.Therefore, we can still obtain the onset of collective mo-tion using the ring-kinetic theory.In order to analyze fluctuations or susceptibilities wewould require to analyze stationary states also within theordered phase. However, the solutions of the ring-kineticequations require a very long time to become stationaryin the ordered regime (when started with only a minimalpolar order). Therefore, we determine the onset of flock-ing within the ring-kinetic theory by a numerical stabilityanalysis of the disordered state. In the numerical solutionof the ring-kinetic equations we start with some (verysmall) initial polar order because the disordered statesare always (for small noise unstable) stationary solutions.We started with an initial polar order of p = 2 π × − and considered the state as polarly ordered if the final po-lar order (after a total time of or ) is larger than p and disordered if p < p in the final state, see Fig. 17.Typically we find also a decreasing trend of p at the endof the observation time if p < p and an increasing trendif p > p . In that way, the ring-kinetic theory predictsthe flocking transition at σ c, rk = 0 . . Surprisingly,that result is already pretty close to the value measuredin direct simulations, even though the pair correlationsare not quantitatively correct. Apparently, consideringcorrelations of the correct order of magnitude, improvesmean field theory significantly.Similar to phase transitions in equilibrium spin sys-tems, correlations shift the flocking transition towardssmaller noise, compared to mean field theory. Thatmeans correlations favor disorder. This can be under-stood qualitatively as follows. If we consider a particlethat moves not in the direction of the majority, withoutcorrelations it would be convinced to join the majority .
45 0 .
50 0 .
55 0 .
60 0 .
65 0 .
70 0 .
75 0 . σ . . . . p p = 2 π × − FIG. 17. Polar order parameter p obtained with ring-kinetictheory. The system was initiated without correlations andwith a minimal polar order of p = 2 π × − ≈ . × − .The graph shows the polar order after a waiting time of T = 250 . Far in the disordered regime the system is alreadyin steady state after this time, that means the correlationsare stationary. However, close to the transition and also inthe polarly ordered regime the systems needs much longer tobecome stationary. In particular at small noise, the polar or-der is still much below its steady state value. Therefore, weestimate the transition by a numerical stability analysis of thedisordered state: we consider the system as disordered if thepolar order parameter after time T is below its initial value p .In that case we typically observe that the polar order param-eter is still decreasing at the end of the observation time. Onthe other hand we consider the system as ordered when thefinal value p is larger than p . In that case we also typicallyobserve that p is still increasing at the end of the observa-tion time. The transition noise strength between polar orderand disorder is displayed as the blue vertical dashed line at σ = 0 . . The black and purple vertical dashed lines showthe transition noise strengths obtained in agent-based simu-lations and in mean field theory, respectively. Parameters areas in Fig. 15 due to interactions with other particles (that on averagebehave as the majority). If correlations are present, theinteraction partners can have the same direction as theconsidered particle (different from the majority) due tocorrelations. That means a particle that is oriented dif-ferently from the global average, is likely to be accompa-nied by other particles that differ from the global averageif strong angular correlations are present. This correla-tion effect weakens the alignment mechanism comparedto mean field.It should be mentioned that the considered systemis relatively small ( N = 509 ). The correlation map inFig. 2 is based on simulations of much larger systems( N = 10 ). In those larger systems, the flocking tran-sition becomes discontinuous and higher order correla-tions or inhomogeneous solutions might become moreimportant. In fact, for larger systems, the transition6noise strength is shifted towards larger noise strengths σ c, sim, large = 0 . [27]. This effect can be understoodby density fluctuations. Because the transition noisestrength of the homogeneous solutions is highly densitydependent, polar order is achieved locally due to an ac-cumulation of particles in a band already before the ho-mogeneous solution becomes polarly ordered. We expect,that the system size dependence of the transition noisebecomes less pronounced for higher densities, because inthis case, the density fluctuations are less important.In order to confirm the applicability of the ring-kinetictheory, we apply it also at different densities. We for-mally assumed the limit N → ∞ in the derivation of thekinetic equations, in particular in the number of neighbordistribution. For small densities we require large sys-tems in order to obtain reasonably large particle num-bers N . However, for those large systems the numericalsolution of the ring-kinetic equations becomes computa-tionally to expensive with the Fourier techniques we areusing. Therefore we focus on larger densities C = 3 and C = 5 . For those parameters, the ring-kinetic theorybecomes unstable when approaching the flocking transi-tion. The reason is, that the number of neighbor distri-bution becomes unphysical, producing negative probabil-ities. This artefact is caused by the fact that in reality,higher order correlations have to be taken into accountin order to predict the correct number of neighbor dis-tribution. Ignoring them can lead to negative prefactors w k in the ring-kinetic equations (44) and (55) that makethe equation unstable. In order to fix this problem we re-quire some information about higher order correlations,such that a reasonable (physical) number of neighbor dis-tribution can be used. In particular we need the localintegrals over the correlation functions: C k := N k (cid:90) G k (1 , ..., k ) k (cid:89) l =1 θ ( R − | r l | )d l, (61) D k := N k − (cid:90) G k (1 , , . . . , k )d1 k (cid:89) l =2 θ ( R − | r − r l | )d l, (62)see [27] for details. We are estimating the next ordercoefficients C and D using a closure ansatz presentedin the next section. VII. CLOSURE RELATION
In order to obtain an estimate for the three particlecorrelation function G we employ the closure ansatz P (1 , ,
3) =Ψ(1 , ,
3) + Ψ(1 , , , , , (63)where Ψ(1 , is a symmetric function Ψ(1 ,
2) = Ψ(2 , (64) that satisfies (cid:90) d2Ψ(1 ,
2) = Γ = const. (65)We assume furthermore that
Ψ(1 , is translational in-variant. Similar closures have been used in astronomyor plasma physics, see e.g. [42–44]. Here, the ansatzis motivated by the limit of small noise. There, nearbyparticles are strongly aligned. If particle one is alignedwith particle two and particle two is aligned with particlethree, then particle one will be also aligned with particlethree. For finite noise however, the quality of the ansatzis not evident a priori. It can be justified only when itleads to reasonable results.Integrating Eq. (63) over all degrees of freedom fixes Γ as Γ = 1 L √ π . (66)Integrating Eq. (63) over the degrees of freedom of par-ticle three leads to Ψ(1 ,
2) = L √ π P (1 , − (cid:90) d3Ψ(1 , , . (67)For a given P we can iterate this equation in order tosolve for Ψ . Starting with the initial guess Ψ = const. = √ πL , the iterative procedure converges very fast, usu-ally within one or two iterations. Having solved for Ψ ,we obtain P from the ansatz (63) and G according toEq. (30). In principle, we could calculate an additionalcollision term with G in the kinetic equation of g (55).Here however, we neglect this term, assuming that theangular dependence of G is sufficiently small. We con-sider only the effect of G on the number of neighbordistribution. That means we calculate C and D ac-cording to Eqs. (61) and (62) and use the number ofneighbor distribution of [27] to calculate the prefactorsin the kinetic equations (44) and (55).Considering the parameters of section VI, with den-sity C = 1 we compare the measured values of the threeparticle correlation parameters C and D with the re-sults of the kinetic theory in Fig. 18. For large noise,the relative difference between theory and simulation islarge. However, in that case the influence of three parti-cle correlations is negligible any way. For smaller noise,coming closer to the flocking transition, the closure be-comes better and agrees very well with the direct simu-lations. Even closer to the transition, the closure clearlyunderestimates the correlations. However, the agreementis satisfactory almost up to the transition. The introduc-tion of the closure relations also improves the agreementof pair correlations with the simulation, cf. Fig. 15. Fur-thermore, the prediction of the transition noise strengthis further improved to σ c, closure = 0 . . It agrees withthe value from direct simulations within the reached un-certainty.7 . . . . σ C (a) . . . . σ . . . . . . D (b) FIG. 18. Local spatial three particle correlation parameters C (a) and D (b) as defined in Eq. (46) compared for directsimulations (black line) and a kinetic theory including a g -closure (red line). The dashed vertical lines shows the onset offlocking measured in direct simulations (black), in the ring-kinetic theory (blue), ring-kinetic theory with closure (red)and in mean field theory (purple). System parameters are asin Fig.15. A. Density Dependence
For larger densities C = 3 and C = 5 we find similarresults. The pair correlations predicted by the kinetictheory agree within about 10 % with the direct simula-tions up to the flocking transition, cf. Fig. 19a and band Fig. 20a and b. Also the three particle correlationspredicted by the closure ansatz have the correct order ofmagnitude, see Fig. 19c and d and Fig. 20c and d.For C = 3 we find the transition noise strength at σ c, closure = 0 . compared to the value measured indirect simulations σ c, sim = 0 . and the mean fieldvalue σ c, mf = 0 . . For C = 5 the obtained transitionnoise strengths are σ c, closure = 0 . , σ c, sim = 0 . and σ c, mf = 0 . . In summary, we find that the role of . . . σ C (a) . . . σ D (b) . . . σ C (c) . . . σ D (d) FIG. 19. Spatial pair correlation parameters C (a) and D (b) and three particle correlation parameters C (c) and D (d) measured in agent-based simulations (black line), withinpure ring-kinetic theory (blue line) and within kinetic the-ory including a three particle closure (red line). The verticaldashed lines display the onset of flocking in kinetic theory withclosure at σ c,closure = 0 . (red), agent-based simulations at σ c,sim = 0 . (black) and mean field theory at σ c,mf = 0 . (purple). Sytem parameters are C = ρπR = 3 , v = R = 1 , N = 300 . three particle correlations increases for larger densities atthe flocking transition. Therefore one expects that alsohigher order correlations become important. This ex-plains that deviations between theory and agent-basedsimulations increase for large densities. Nevertheless,we still find satisfactory quantitative agreement with thesimulations up to C = 5 . The introduction of a closurerelation for spatial three particle correlations enlargedthe parameter region of applicability of the kinetic the-ory significantly compared to the pure ring-kinetic theorywithout closure. B. Velocity Dependence
We study the parameter set of Fig. 15 at C = ρπR =1 for different velocities. For a smaller velocity of v = 0 . we find larger correlations, see Fig.21. This is expectedbecause particles have more time to interact and alignat smaller velocities. The predictions of the pure ring-kinetic theory do not agree as well as for v = 1 similaras for larger densities. These deviations are very likelycaused by higher order correlations that need to be takeninto account. As for larger densities, the ring-kinetic the-ory is unstable for small noise strengths and the descrip-tion breaks down already above the flocking transition.Including the closure improves the kinetic theory signifi-cantly. However, there are significant deviations of pair-8 σ C (a) 0.7 0.8 0.9 σ D (b)0.7 0.8 0.9 σ C (c) 0.7 0.8 0.9 σ D (d) FIG. 20. Spatial pair correlation parameters C (a) and D (b) and three particle correlation parameters C (c) and D (d) measured in agent-based simulations (black line), withinpure ring-kinetic theory (blue line) and within kinetic the-ory including a three particle closure (red line). The verticaldashed lines display the onset of flocking in kinetic theory withclosure at σ c,closure = 0 . (red), agent-based simulations at σ c,sim = 0 . (black) and mean field theory at σ c,mf = 0 . (purple). Sytem parameters are C = ρπR = 5 , v = R = 1 , N = 500 . and three particle correlations when the flocking transi-tion is approached. Nevertheless, the flocking transitionpredicted by our kinetic theory is a clear improvementcompared to mean field theory.For a larger velocity, v = 1 . we find smaller spatialcorrelations as expected, see Fig. 22. The spatial paircorrelations agree also very well between ring-kinetic the-ory and agent-based simulations almost until the onsetof flocking. Here, the ring-kinetic theory is stable andpredicts the flocking the transition at slightly too smallnoise. Thus, the closure improves the results slightly.In principle, we would expect better agreement due tothe smaller correlations. However, only the pure spatialcorrelations decrease for larger velocities, but the angularcorrelations increase because particles can only stay closeto each other for a long time if they are aligned at highvelocities. For even larger velocity (e.g. v = 2 ) we ob-serve very strong angular correlations that suppress theonset of flocking. That means we find no onset of flockingwithin the ring-kinetic theory (with or without closure)although the spatial correlations are predicted reasonablywell. We believe that higher order angular correlationsneed to be taken into account for high velocities. Thenext step in this direction will be the incorporation of g -collision integrals into Eq. (55), where g is still deter-mined by the closure (63). . . . . . σ C (a) . . . . . σ D (c) . . . . . σ C (b) . . . . . σ D (d) FIG. 21. Spatial pair correlation parameters C (a) and D (b) and three particle correlation parameters C (c) and D (d) measured in agent-based simulations (black line), withinpure ring-kinetic theory (blue line) and within kinetic the-ory including a three particle closure (red line). The verticaldashed lines display the onset of flocking in kinetic theory withclosure at σ c,closure = 0 . (red), agent-based simulations at σ c,sim = 0 . (black) and mean field theory at σ c,mf = 0 . (purple). Sytem parameters are C = ρπR = 1 , v = 0 . , R = 1 , N = 509 . VIII. DISCUSSION
We consider polarly aligning self-propelled point par-ticles in two dimensions. The investigated models are inthe spirit of the famous Vicsek models. However, theyfollow a continuous time dynamics given by a system ofgeneralized Langevin equations. For technical reasonsit would be desirable to introduce only pair interactions.But in that case the model behaves qualitatively differentfrom the Vicsek model: in the ordered phase, no bands orcross sea patterns are observed. Instead, strongly alignedhigh density clusters are formed.In order to observe a behavior qualitatively equivalentto the Vicsek model, we need to introduce N -particle-,that is non-additive interactions as pointed out recentlyin [36, 37].The presence of N -particle interactions seriously com-plicates matters. As a consequence, the kinetic equationscontain weight factors that are expectation values of thenumber of neighbor distribution. It was discovered re-cently, how this distribution can be calculated exactlyeven in the presence of many particle correlations [27].Employing this predecessor work, we were able to handlethe complications arising from the non-additive interac-tions. From a technical point of view, the incorporationof N -particle interactions into ring-kinetic theories is asignificant development that can be of interest also for9 . . . . . σ . . . . C (a) . . . . . σ . . . . . D (c) . . . . . σ C (b) . . . . . σ D (d) FIG. 22. Spatial pair correlation parameters C (a) and D (b) and three particle correlation parameters C (c) and D (d) measured in agent-based simulations (black line), withinpure ring-kinetic theory (blue line) and within kinetic the-ory including a three particle closure (red line). The verti-cal dashed lines display the onset of flocking in pure ring-kinetic theory at σ c,rk = 0 . (blue), kinetic theory withclosure at σ c,closure = 0 . (red), agent-based simulations at σ c,sim = 0 . (black) and mean field theory at σ c,mf = 0 . (purple). System parameters are C = ρπR = 1 , v = 1 . , R = 1 , N = 509 . other (possibly passive) systems with interactions of thistype.Truncating the BBGKY-hierarchy after the secondequation we obtain the time evolution equations for theone particle distribution and the pair correlation func-tion, neglecting higher order correlations. We solve thoseequations numerically, transforming both, spatial and an-gular dependence, to Fourier space. We compare the re-sulting steady states with direct, agent-based simulationsof the Langevin equations. In the disordered phase andnot too close to the onset of flocking, we find excellentagreement between ring-kinetic theory and direct simu-lations.Reducing the noise strength, we find the onset of flock-ing in the ring-kinetic theory. Compared to the homo-geneous mean field theory, the transition is shifted to-wards smaller noise. This shift is caused by positiveangular pair correlations that stabilize the disorderedphase. The same effect is seen in direct simulations. For C = ρπR = 1 we find good quantitative agreementof the onset of flocking between ring-kinetic theory andagent-based simulations. For larger densities, close tothe onset of flocking, the ring-kinetic equations becomeunstable because the weights depending on the numberof neighbor distribution are predicted wrong by the ring-kinetic theory. This is not too surprising as it is knownthat higher order correlations are required in order to predict the correct number of neighbor distribution andthus the correct weights within the ring-kinetic equations[27].In order to enlarge the applicability of the kinetic the-ory we introduce a closure ansatz. With that relation wecalculate the three particle correlation function g fromthe single particle distribution and the pair correlations.In this work, we do not consider the full effect of the threeparticle correlations, but only its influence on the num-ber of neighbor distribution and hence on the weights inthe kinetic equations. Phrased differently, we consideronly the effect of three particle correlations on the spa-tial distribution and we neglect angular three particlecorrelations that are entering the kinetic equations byadditional collision integrals.Far in the disordered phase, for large noise, the closureansatz does not agree very well with direct simulations.However, in that region three particle correlations canbe ignored as discussed above. At about the onset offlocking, the closure ansatz for the spatial three particlecorrelations agrees within about with direct simula-tions. That means it is a considerable improvement com-pared to the neglect of those correlations. Furthermore,the introduction of the closure significantly enlarges theparameter regime where spatial pair correlations agreequantitatively with direct simulations.Employing the closure ansatz we are able to describethe onset of flocking within kinetic theory also for largerdensities such as C = ρπR = 3 , with deviations of theflocking noise strength of about and . In compar-ison, the deviations of the transition noise in mean fieldtheory are and , respectively.We also consider the dependence on the particle veloc-ity and find good agreement between ring-kinetic theoryand simulations at large noise for all considered velocities.Close to the onset of flocking, small velocities increase inparticular spatial correlation whereas high velocities sup-press spatial correlations but favor angular correlations.This is intuitive because at high speed, particles remainclose to each other for long times only if they have roughlythe same velocity, and for small velocities, particles thatinteract have a long time to align. Both effects lead todeviations of the predicted onset of flocking in the ki-netic theory with three particle closure. However, theresults of the kinetic theory are still a major improve-ment over mean field theory. There is further potentialto improve the theory for high velocities by incorporating g collision integrals in the g equation by means of thepresented closure ansatz.In general, it can be seen that the ring-kinetic theory(possibly extended by a closure relation estimating g )gives quantitatively good results as long as higher ordercorrelations are not too large. Far in the disordered phase(for large noise), this is always the case. Close to theflocking transition we observe larger spatial correlations0for higher particle densities or small velocities and largeangular correlations for high velocities.Note, in this work we considered relatively small sys-tems consisting of only a few hundred particles. Thereason is that the numerical solution of the kinetic equa-tions is too time consuming for much larger systems be-cause the number of necessary spatial Fourier modes isproportional to the area of the system. Our analysis pre-dicts that the correlations are mainly localized withinthe interaction region. That suggests the use of differentspatial decompositions like e.g. into Hermite functionsinstead of Fourier modes. ACKNOWLEDGMENTS
The authors gratefully acknowledge the GWK supportfor funding this project by providing computing timethrough the Center for Information Services and HPC(ZIH) at TU Dresden on the HRSK-II. The authors grate-fully acknowledge the Universitätsrechenzentrum Greif-swald for providing computing. We thank Aurelio Patelli,Fernando Peruani and Sven Stroteich for valuable discus-sions.
Appendix A: Fourier Transform1. Time Evolution Equations
We solve the time evolution equations (44) and (51) nu-merically in Fourier space, transforming both, spatial andangular coordinates. This has the huge advantage thatall appearing integrals (also high dimensional ones) canbe solved analytically. We employ the following Fourieransatz for the one particle distribution and the pair cor-relation function. p ( φ ) = (cid:88) k A k exp( ikφ ) , (A1) g ( φ , φ , ∆) = (cid:88) k,l,m,n F k,l,m,n exp( ikφ ) exp( ilφ ) × exp( im ∆ x π/L ) exp( in ∆ y π/L ) . (A2)The number of neighbor distribution that determines theweights in Eqs. (44) and (51) is known analytically up toa discrete Fourier transform [27]. It depends only on theintegrals (46) when three particle and higher order corre-lations are neglected. Those integrals can be evaluated in Fourier space. It is useful to introduce the abbreviations K m,n := 1 L (cid:90) L/ − L/ (cid:90) L/ − L/ e im πL x e in πL y × θ ( R − | x | )d x d y = RL √ m + n × J (cid:16) πL R (cid:112) m + n (cid:17) , (A3) Q ± klmn := (cid:88) s,t F klst K m ± s,n ± t (A4)and R ± klmn := (cid:88) s,t F klst K m ± s,n ± t C s,t , (A5)where J is the Bessel function of the first kind. Withthose abbreviations one obtains C =4 π N R +0000 , (A6) D =4 π N Q +0000 . (A7)Thus, with the results of [27], the weights in Eqs. (44) and(51) can be expressed explicitly in terms of the Fouriermodes F klmn . Fourier transforming those time evolutionequations results in ∂ t A k = N w kπ ( Q + k − , , , − Q + k +1 , − , , )+ C w kπ ( A k − A − A k +1 A − )+ N π k ( w − w )( A k − R +1 , , , − A k +1 R + − , , , )+ C N π k ( w − w )( A Q + k − , , , − A − Q + k +1 , , , )+ N π k ( w − w + w ) × ( Q + k − , , , R +1 , , , − Q + k +1 , , , R + − , , , ) − σ k A k . (A8)and ∂ t F klmn = (cid:88) i =1 (cid:13) i + { k ↔ l, m ↔ − m, n ↔ − n } , (A9)where the terms (cid:13) i are given by (cid:13) w k Q − k − ,l +1 ,m,n − Q − k +1 ,l − ,m,n ) , (A10) (cid:13) w k K mn ( A k − A l +1 − A k +1 A l − ) , (A11) (cid:13) N kw π ( A k − R + l, ,m,n − A k +1 R + l, − ,m,n ) , (A12)1 (cid:13) C πkw ( Q − k − ,l,m,n A − Q − k +1 ,l,m,n A − ) , (A13) (cid:13) N πk ( w − w ) K mn A l ( Q + k − , , , − Q + k +1 , − , , ) , (A14) (cid:13) N ( w − w ) k π × ( Q − k − ,l,m,n R +1 , , , − Q − k +1 ,l,m,n R + − , , , ) , (A15) (cid:13) N ( w − w )2 π k × ( R + l, ,m,n Q + k − , , , − R + l, − ,m,n Q + k +1 , , , ) , (A16) (cid:13) N ( w − w ) k π × R + l, ,m,n ( Q + k − , , , − Q + k +1 , − , , ) , (A17) (cid:13) C N ( w − w ) k π × R + l, ,m,n ( A k − A − A k +1 A − ) , (A18) (cid:13)
10 = N ( w − w + w ) k π A l K mn × ( A k − R +1 , , , − A k +1 R + − , , , ) , (A19) (cid:13)
11 = C N k π ( w − w + w ) A l K mn × ( Q + k − , , , A − Q + k +1 , , , A − ) , (A20) (cid:13)
12 = N C π k ( w − w + w ) R + l, ,m,n × ( Q + k − , , , A − Q + k +1 , , , A − ) , (A21) (cid:13)
13 = N k π ( w − w + w ) R + l, ,m,n × ( A k − R +1 , , , − A k +1 R + − , , , ) , (A22) (cid:13)
14 = N k π ( w − w + 3 w − w K mn A l × ( Q + k − , , , R +1 , , , − Q + k +1 , , , R + − , , , , (A23) (cid:13)
15 = N π k ( w − w + 3 w − w ) R + l, ,m,n × ( Q + k − , , , R +1 , , , − Q + k +1 , , , R + − , , , ) , (A24) (cid:13)
16 =
N kπw K mn × ( A k − F l, , − m, − n − A k +1 F l, − , − m, − n ) , (A25) (cid:13)
17 = − N kπw ( A k − R + l, ,m,n − A k +1 R + l, − ,m,n ) (A26) (cid:13)
18 = C w kπ ( F k − ,l,m,n A − F k +1 ,l,m,n A − ) , (A27) (cid:13)
19 = − C w kπ ( Q − k − ,l,m,n A − Q − k +1 ,l,m,n A − ) , (A28) (cid:13)
20 = N π k ( w − w ) × ( R +1 , , , F k − ,l,m,n − R + − , , , F k +1 ,l,m,n ) , (A29) (cid:13)
21 = − N π k ( w − w ) × ( Q − k − ,l,m,n R +1 , , , − Q − k +1 ,l,m,n R + − , , , ) , (A30) (cid:13)
22 = N k π ( w − w ) K mn × ( Q + k − , , , F l, , − m, − n − Q + k +1 , , , F l, − , − m, − n ) , (A31) (cid:13)
23 = − N k π ( w − w ) × ( Q + k − , , , R + l, ,m,n − Q + k +1 , , , R + l, − ,m,n ) , (A32) (cid:13)
24 = N k π ( w − w ) K mn × F l, , − m, − n ( Q + k − , , , − Q + k +1 , − , , ) , (A33) (cid:13)
25 = − N k π ( w − w ) R + l, ,m,n × ( Q + k − , , , − Q + k +1 , − , , ) , (A34) (cid:13)
26 = C N k π ( w − w ) K mn × F l, , − m, − n ( A k − A − A k +1 A − ) , (A35) (cid:13)
27 = − C N k π ( w − w ) × R + l, ,m,n ( A k − A − A k +1 A − ) , (A36) (cid:13)
28 = C N k π ( w − w + w ) K mn × F l, , − m, − n ( Q + k − , , , A − Q + k +1 , , , A − ) , (A37) (cid:13) − C N k π ( w − w + w ) × R + l, ,m,n ( Q + k − , , , A − Q + k +1 , , , A − ) , (A38)2 (cid:13)
30 = N k π ( w − w + w ) K mn × F l, , − m, − n ( A k − R +1 , , , − A k +1 R + − , , , ) , (A39) (cid:13)
31 = − N k π ( w − w + w ) × R + l, ,m,n ( A k − R +1 , , , − A k +1 R + − , , , ) , (A40) (cid:13)
32 = N k π ( w − w + 3 w − w ) F l, , − m, − n × K mn ( Q + k − , , , R +1 , , , − Q + k +1 , , , R + − , , , ) , (A41) (cid:13)
33 = − N k π ( w − w + 3 w − w ) R + l, ,m,n × ( Q + k − , , , R +1 , , , − Q + k +1 , , , R + − , , , ) , (A42) (cid:13)
34 = vm iπL ( F k − ,l,m,n + F k +1 ,l,m,n ) , (A43) (cid:13)
35 = vn πL ( F k − ,l,m,n − F k +1 ,l,m,n ) , (A44) (cid:13)
36 = − σ k F klmn (A45)and { k ↔ l, m ↔ − m, n ↔ − n } refers to the same termsbut with k and l , m and − m , n and − n interchanged,respectively.
2. Closure Relation
In order to evaluate the closure ansatz (63) we considerthe function Ψ in Fourier space Ψ( φ , φ , ∆) = (cid:88) klmn U klmn exp( ikφ ) × exp( ilφ ) exp( im ∆ x π/L ) exp( in ∆ y π/L ) . (A46)It is reasonable to consider scaled modes u stuv = L U stuv . (A47)From Eqs. (65) and (66) it follows that u s = δ s, u = δ s, L Γ2 π = δ s, π √ π . (A48)Fourier transforming the iteration equation (67) we ob-tain u klmn = 12 u (cid:20) π δ k, δ l, δ m, δ n, − (cid:88) s u ksmn u − slmn + 12 π F klmn (cid:21) . (A49) Note that this equation is compatible with the nor-malization condition (A48). In practice, we start with u klmn = 0 except for u that is given by Eq. (A48).Then we iterate Eq. (A49) to calculate u klmn , that is theclosure ansatz function Ψ . In the parameter regime weinvestigated, the recursion converges very fast. It almostreaches its fixed point already after two iterations. Weused five iterations after which we reach perfect conver-gence within our numerical accuracy.Once Ψ is known, the spatial three particle correlationparameters C and D can be calculated via C + 3 C C + C = N (cid:90) P (1 , , N (cid:90) Ψ(1 , , θ θ θ d1d2d3= 3 ρ (cid:90) L Ψ(1 , , θ θ θ d1d2d3= 3 ρ (cid:90) (cid:88) stuv exp( isφ ) exp( itφ ) exp( iu ( x − x )2 π/L ) × exp( i ( y − y )2 π/L ) (cid:88) klmn exp( ikφ ) exp( ilφ ) × exp( im ( x − x )2 π/L ) × exp( in ( y − y )2 π/L ) θ θ θ d1d2d3= 3 N (2 π ) (cid:88) tuvmn u ,t,u,v u − t, ,m,n K u,v K u − m,v − n K m,n (A50)and similar D + 2 D C + C + C = N (cid:90) θ θ G (1 , , π ) N (cid:88) lmnuv u lmn u − l uv K m − u,n − v K u,v + (2 π ) N (cid:88) lmnuv u lmn u − luv K m + u,n + v K u,v + (2 π ) N (cid:88) lmnuv u − l mn u l uv K m,n K u,v . (A51)Note that C := NL πR and D := 1 per definition.In sumary, the numerical time evolution works as fol-lows. We calculate the pair correlation coefficients C and D according to Eq. (A7). Next, we calculate theclosure ansatz function with the iteration (A49) and thenthe three particle correlation coefficients C and D ac-cording to Eqs. (A50) and (A51). Having calculatedthe correlation coefficients we compute the number ofneighbor distribution as described in detail in [27]. Withthis distribution we obtain the weights w k according toEq. (45). Given the weights we can eventually time evolvethe time evolution equations (A8) and (A9) with a simpleEuler scheme.3
3. Comparison to Agent-based Simulations
The correlation functions g ( r ) and h ( r, α, ∆ φ ) consid-ered in Sec. VI can be sampled in agent-based simulationsaccording to Eqs. (57) and (60). On the other hand, thesefunctions can be calculated from the Fourier modes of g according to g ( r ) = 1 + 4 π (cid:88) m,n F mn J ( 2 πrL (cid:112) m + n ) , (A52) h ( r, α, ∆ φ ) =(2 π ) (cid:88) klmn exp( il ∆ φ ) exp[ − i ( k + l ) α ] × (cid:18) n − im √ n + m (cid:19) k + l J − ( k + l ) (cid:18) πL r (cid:112) m + n (cid:19) . (A53)Assuming rotational symmetry, it is also possible to re-vert Eq. (A52) that is to calculate the spatial Fouriermodes F mn from the function g ( r ) according to F mn = 1 L π (cid:90) [ g ( r ) − rJ (cid:18) πrL (cid:112) m + n (cid:19) d r. (A54)We use this relation to see the impact of the resolutionin the Fourier transform on g ( r ) by Fourier transformingthe measured function according to Eq. (A54) and trans-forming it back according to Eq. (A52). In that way weachieved the green dashed line in Fig. 3. The solid blackline in this figure is a histogram of agent-based simula-tions according to Eq. (57) and the solid blue line is cal-culated according to Eq. (A52), where the Fourier modes F klmn are the result of the ring-kinetic theory. Appendix B: Comparison between full Ring-kineticTheory and Landau Kinetic Theory
The Landau theory of [26] considers the model withadditive interactions, w ( n ) ≡ . Here, we adopt the the-ory to the case of non-additive interactions, w ( n ) = 1 /n in order to compare it with the present approach. Theequivalent of Eq. (17) of Ref. [26] is in our notation ∂ t g ( φ , φ , ∆) = − w (cid:20) + (cid:21) + v (cos φ − cos φ ) ∂ ∆ x g ( φ , φ , ∆)+ v (sin φ − sin φ ) ∂ ∆ y g ( φ , φ , ∆) . (B1)That is, compared to Eq. (55), all collision integrals onthe right hand side that contain g are neglected be-cause g and the coupling constant are both assumedto be small. Furthermore, the angular diffusion term σ ( ∂ φ + ∂ φ ) g ( φ , φ , ∆) is neglected because both, σ and g , are assumed to be small. However, we could notmake sense out of Eq. (B1) without the angular diffusionterm. Intuitively, it is clear that the correlations are go-ing to diverge according to Eq. (B1) because there is asource term but no damping term. In fact, we find in-deed a diverging g if we integrate Eq. (B1) numericallyin Fourier space. Furthermore, we do not see how theprincipal problem of diverging g can be prevented in anyparameter limit if the noise term is neglected. Thereforewe did not neglect the angular diffusion and consideredinstead the time evolution equation ∂ t g ( φ , φ , ∆) = − w (cid:20) + (cid:21) + v (cos φ − cos φ ) ∂ ∆ x g ( φ , φ , ∆)+ v (sin φ − sin φ ) ∂ ∆ y g ( φ , φ , ∆)+ σ ∂ φ + ∂ φ ) g ( φ , φ , ∆) . (B2)In Fig.23a we show the pair correlation function g ( r ) ob-tained from the Landau kinetic theory (purple line) com-pared to the results of the full ring-kinetic theory (blueline) and direct agent-based simulations (black line) forthe parameters of Fig.3 far in the disordered phase at σ = 1 . (the flocking transition occurs at σ = 0 . ).We find that the qualitative behavior of the pair corre-lation function is predicted correctly, but it is too smallby about . We expect that the results of the Landaukinetic theory improve if much larger noise strengths areconsidered. In that case however, the correlations de-crease and become less important.In Fig.23b we display the pair correlation function forthe same parameters but smaller noise, σ = 0 . . Thesystem is still disordered but closer to the flocking tran-sition. We added also the curve of our kinetic theoryincluding a three particle closure (red line). Here, we seethat the Landau theory is off by more than , the ring-kinetic theory is much better but also shows significantdeviations and the theory including spatial three parti-cle correlations by a closure ansatz agrees quantitativelyvery well with the agent-based simulations. We concludethat the simplified Landau kinetic theory of [26] is a goodapproximation far in the disordered regime. However, itis not suitable to describe the system in the vicinity ofthe flocking transition.4 r . . . . . g ( r ) (a) r . . . . . . g ( r ) (b) FIG. 23. Radial distribution function g ( r ) obtained fromagent-based simulations (black line), ring-kinetic theory (blueline), kinetic theory with closure (red line) and Landau-kinetictheory according to Eq. (B2) (purple line). In (a) we used theparameters of Fig.(3), in (b) the noise is smaller σ = 0 . butstill in the disordered phase, other parameters are the same.[1] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, andO. Shochet, Novel type of phase transition in a system ofself-driven particles, Phys. Rev. Lett. , 1226 (1995).[2] J. Toner and Y. Tu, Long-range order in a two-dimensional dynamical xy model: how birds fly together,Phys. Rev. Lett. , 4326 (1995).[3] Y. Fily and M. C. Marchetti, Athermal phase separationof self-propelled particles with no alignment, Phys. Rev.Lett. , 235702 (2012).[4] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen,C. Bechinger, and T. Speck, Dynamical clustering andphase separation in suspensions of self-propelled colloidalparticles, Phys. Rev. Lett. , 238301 (2013).[5] C. W. Wolgemuth, Collective swimming and the dynam-ics of bacterial turbulence, Biophys. J. , 1564 (2008).[6] J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics and phases of flocks, Ann. Phys. , 170 (2005).[7] S. Ramaswamy, The mechanics and statistics of activematter, Annu. Rev. Condens. Matter Phys. , 323 (2010).[8] T. Vicsek and A. Zafeiris, Collective motion, Phys. Rep. , 71 (2012).[9] P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, andL. Schimansky-Geier, Active brownian particles, Eur.Phys. J. - Spec. Top. , 1 (2012).[10] D. Klotsa, As above, so below, and also in between:mesoscale active matter in fluids, Soft Matter , 8946(2019).[11] O. Dauchot and H. Löwen, Chemical physics of activematter, J. Chem. Phys. , 114901 (2019).[12] H. Chaté, Dry aligning dilute active matter, Annu. Rev.Condens. Matter Phys. , 189 (2020).[13] M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gomp- per, and H. Rieger, Computational models for active mat-ter, Nature Rev. Phys. , 181 (2020).[14] M. Bär, R. Großmann, S. Heidenreich, and F. Peru-ani, Self-propelled rods: Insights and perspectives for ac-tive matter, Annu. Rev. Condens. Matter Phys. , 441(2020).[15] G. Grégoire and H. Chaté, Onset of collective and cohe-sive motion, Phys. Rev. Lett. , 025702 (2004).[16] H. Chaté, F. Ginelli, G. Grégoire, and F. Raynaud, Col-lective motion of self-propelled particles interacting with-out cohesion, Phys. Rev. E , 046113 (2008).[17] A. P. Solon, H. Chaté, and J. Tailleur, From phase tomicrophase separation in flocking models: The essentialrole of nonequilibrium fluctuations, Phys. Rev. Lett. ,068101 (2015).[18] R. Kürsten and T. Ihle, Dry active matter exhibits a self-organized cross sea phase, Phys. Rev. Lett. , 188003(2020).[19] J. Toner and Y. Tu, Flocks, herds, and schools: A quan-titative theory of flocking, Phys. Rev. E , 4828 (1998).[20] T. Ihle, Kinetic theory of flocking: Derivation of hydro-dynamic equations, Phys. Rev. E , 030901 (2011).[21] T. Ihle, Invasion-wave-induced first-order phase transi-tion in systems of active particles, Phys. Rev. E ,040303 (2013).[22] T. Ihle, Chapman–enskog expansion for the vicsek modelof self-propelled particles, J. Stat. Mech. - Theory E. , 083205 (2016).[23] Y.-L. Chou and T. Ihle, Active matter beyond mean-field: Ring-kinetic theory for self-propelled particles,Phys. Rev. E , 022103 (2015).[24] J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo,and A. Morozov, Role of correlations in the collectivebehavior of microswimmer suspensions, Phys. Rev. Lett. , 028005 (2017).[25] V. Škultéty, C. Nardini, J. Stenhammar, D. Marenduzzo,and A. Morozov, Swimming suppresses correlations indilute suspensions of pusher microorganisms, Phys. Rev.X , 031059 (2020).[26] A. Patelli, Landau kinetic equation for dry aligning activemodels, arXiv preprint arXiv:2010.12213 (2020).[27] R. Kürsten, S. Stroteich, M. Zumaya-Hernández, andT. Ihle, Multiple particle correlation analysis of many-particle systems: Formalism and application to activematter, Phys. Rev. Lett. , 088002 (2020).[28] M. Ernst and E. Cohen, Nonequilibrium fluctuations in µ space, J. Stat. Phys. , 153 (1981). [29] E. Leutheusser, Self-consistent kinetic theory for thelorentz gas, Phys. Rev. A , 1762 (1983).[30] H. Bussemaker, M. Ernst, and J. Dufty, Generalizedboltzmann equation for lattice gas automata, J. Stat.Phys. , 1521 (1995).[31] T. Van Noije, M. Ernst, and R. Brito, Ring kinetic theoryfor an idealized granular gas, Physica A , 266 (1998).[32] F. Peruani, A. Deutsch, and M. Bär, A mean-field theoryfor self-propelled particles interacting by velocity align-ment mechanisms, Eur. Phys. J. - Spec. Top. , 111(2008).[33] F. Peruani, L. Schimansky-Geier, and M. Bär, Clusterdynamics and cluster size distributions in systems of self-propelled particles, Eur. Phys. J. - Spec. Top. , 173(2010).[34] F. D. C. Farrell, M. C. Marchetti, D. Marenduzzo, andJ. Tailleur, Pattern formation in self-propelled particleswith density-dependent motility, Phys. Rev. Lett. ,248101 (2012).[35] B. Liebchen and D. Levis, Collective behavior of chiralactive matter: Pattern formation and enhanced flocking,Phys. Rev. Lett. , 058002 (2017).[36] O. Chepizhko, D. Saintillan, and F. Peruani, Revisitingthe emergence of order in active matter, Soft Matter ,(2021).[37] S. Stroteich, Linear response and closure methods: acomputer-assisted search for a quantitative theory of cor-related active particle systems, Masterthesis, Universityof Greifswald (2019).[38] P. E. Kloeden and E. Platen, Numerical Solutionof Stochastic Differential Equations (Springer-Verlag,Berlin, 1992).[39] H. Risken,
The Fokker-Planck Equation , 2nd ed.(Springer-Verlag, Berlin, 1989).[40] H. D. Ursell, The evaluation of gibbs’ phase-integral forimperfect gases, Math. Proc. Cambridge , 685–697(1927).[41] J. E. Mayer and E. Montroll, Molecular distribution, J.Chem. Phys. , 2 (1941).[42] M. Davis and P. Peebles, On the integration of the bbgkyequations for the development of strongly nonlinear clus-tering in an expanding universe, Astrophys. J. Suppl. S. , 425 (1977).[43] S. D. M. White, The hierarchy of correlation functionsand its relation to other measures of galaxy clustering,Mon. Not. R. Astr. Soc. , 145 (1979).[44] T. O’Neil and N. Rostoker, Triplet correlation for aplasma, Phys. Fluids8