Self-propelled particles with selective attraction-repulsion interaction - From microscopic dynamics to coarse-grained theories

In this work we derive and analyze coarse-grained descriptions of self-propelled particles with selective attraction-repulsion interaction, where individuals may respond differently to their neighbours depending on their relative state of motion (approach versus movement away). Based on the formulation of a nonlinear Fokker-Planck equation, we derive a kinetic description of the system dynamics in terms of equations for the Fourier modes of a one-particle density function. This approach allows effective numerical investigation of the stability of possible solutions of the system. The detailed analysis of the interaction integrals entering the equations demonstrates that divergences at small wavelengths can appear at arbitrary expansion orders. Further on, we also derive a hydrodynamic theory by performing a closure at the level of the second Fourier mode of the one-particle density function. We show that the general form of equations is in agreement with the theory formulated by Toner and Tu. Finally, we compare our analytical predictions on the stability of the disordered homogeneous solution with results of individual-based simulations. They show good agreement for sufficiently large densities and non-negligible short-ranged repulsion. Disagreements of numerical results and the hydrodynamic theory for weak short-ranged repulsion reveal the existence of a previously unknown phase of the model consisting of dense, nematically aligned filaments, which cannot be accounted for by the present Toner and Tu type theory of polar active matter.


Introduction
In recent decades, there has been an increased research focus on far-from-equilibrium systems in biology and physics which is referred to as "active matter". The relevant length scales of such systems span several orders of magnitude. They range from the (sub-)micrometer scale governing the dynamics of individual active units in motility assays in-vitro [1] and the actin cortex in-vivo [2], via the mesoscopic length scales of interest in collective dynamics of large bacterial ensembles [3,4] and artificial selfpropelled particles [5], up to the macroscopic scales of driven granular matter [6,7] or flocks of birds [8], schools of fish [9], and swarms of insects [10,11], where the spatial dimensions can be in the order of kilometers.
Despite the apparent variety, all these systems share the fundamental property of local uptake and/or conversion of internal energy into kinetic energy of motion by its individual units. This -together with additional interactions between those units -distinguishes these systems from related equilibrium systems and yields fascinating examples of self-organization and collective dynamics.
The question of universal properties of such active matter systems from the statistical physics point of view is a vibrant research field. To a large extent, it was initiated by the numerical study of a minimal, individual-based model of active matter published by Vicsek et al in 1995 [12]. Shortly after this publication, Toner and Tu made a seminal contribution by formulating the hydrodynamic equations of polar active matter at largest relevant length and time scales purely based on symmetry arguments [13,14]. The analysis of these generic equations, as well as their counterparts for nematic order, improved our understanding of the fundamental properties of active matter, such as the existence of long range order or giant number fluctuations [15,16].
However, the direct derivation of a hydrodynamic theory of the Toner and Tu type from microscopic models of active matter was a long standing problem. Only recently, such a link between microscopic parameters determining the dynamics of individual active units and parameters governing the macroscopic flow of active matter was established by formulating kinetic equations for minimal models of selfpropelled particles with velocity-alignment [17][18][19][20] and self-propelled aligning rods [21]. Furthermore, coarse-grained descriptions for active particles with variable speeds and velocity-alignment were derived in [22][23][24][25][26].
In this paper, we will derive a kinetic and a hydrodynamic description for selfpropelled agents interacting via a selective attraction and repulsion interaction. A corresponding model was recently introduced to describe the onset of collective motion in insect swarms and is directly motivated by response of individual agents to looming visual stimuli and in particular the distinction of approaching and moving away objects [27,28].
A fundamental difference of our model and the Vicsek model is its formulation in continuous time using stochastic differential equations. The interaction of individuals are modelled by superposition of (binary) interaction forces. This allows on the one hand a straight forward generalization of the model, e.g. towards variable speed of individuals, and on the other hand can be directly coarse-grained by deriving the Fokker-Planck equation for the corresponding probability density functions. The latter being the main focus of this work.
The model shows different phases including large scale collective motion, a disordered clustering phase and a nematic phase despite the absence of an explicit velocity-alignment interaction.
First, we will introduce the microscopic, individual-based model in terms of stochastic differential equations. Then we will proceed with the discussion of a kinetic description of the collective dynamics based on the nonlinear Fokker-Planck equation for the one-particle density function, which allows efficient numerical analysis of the stability of solutions of the Fokker-Planck equation in Fourier space. Further on, we will derive a hydrodynamic theory for self-propelled particles with selective attractionrepulsion interaction, which yields hydrodynamic equations in agreement with the theory by Toner and Tu. A direct comparison of the kinetic approach, which in principle can be considered up to arbitrary accuracy, to the hydrodynamic theory, which corresponds to a closure at the level of the second Fourier mode of the probability density function, reveals the range of validity of the hydrodynamic equations at large wavenumbers. Moreover, the analysis provides insights into the origin of unphysical divergences at large wavenumbers related to the approximations used namely the usage of Taylor polynomials. Finally, we compare the results of the kinetic and hydrodynamic theory with direct numerical simulations of the individual-based model.

Microscopic Model
We consider N self-propelled particles of mass m = 1 in two spatial dimensions, socalled agents. Each individual moves at a constant speed s 0 , thus the velocity vector of each agent is determined by its polar orientation angle ϕ i . The equations of motion for the positions r i and the polar orientation angles ϕ i read: Here, F i,ϕ = F i · e i,ϕ is the projection of an effective social force vector F i on the angular degree of freedom ϕ i , which induces a turning behaviour of the focal individual due to interaction with others. e i,ϕ = (− sin ϕ i , cos ϕ i ) T is the angular unit vector perpendicular to e i,h . The second term within the brackets in (2.1b) stands for random angular noise with intensity D ϕ . ξ i (t) are independent, Gaussian random processes with vanishing mean and temporal δ-correlations, i.e. ξ i (t)ξ j (t + τ ) = δ ij δ(τ ) (Gaussian white noise).
The total social force is given by a sum of three components: (2. 2) The first term represents a short-ranged repulsion allowing for finite sized agents. It reads with µ r (r ji ) ≥ 0 being a repulsive turning rate, which in general depends on distance r ji = |r ji | = |r j − r i | between two particles. We assume that this repulsive interaction is strictly short-ranged, i.e. it vanishes above a finite repulsive radius l c , as indicated by the Heaviside (unit-step) function θ(x). The other two forces read: These forces can be considered as a sum over binary interactions, which act always along the unit vectorr ji = (r j − r i )/r ji pointing towards the center of mass of the neighbouring particle j. The corresponding response strengths µ a,m (r ji ) are distance dependent and may in general be both positive (attraction) and negative (repulsion). Furthermore, the overall response to other individuals is assumed to vanish above a finite sensory range l s : µ a,m (r ji > l s ) = 0, whereby l s ≥ l c . The decisive factor in the distinction of the two forces is the sign of the relative velocityṽ ji defined by the temporal derivative of the distanceṙ ji between particles i and j, and equals the projection of the velocity difference v ji = v j − v i of neighbour j and the focal individual i on the relative position unit vectorr ji :ṽ ji = v ji ·r ji =ṽ ij . Hence f i,a is the response to approaching individuals characterized by a negative relative velocityṽ ji < 0, whereas f i,m is the corresponding response to moving away (or receding) individuals characterized by positive relative velocityṽ ji > 0. Both force terms are proportional to the absolute value of the relative velocity, leading to stronger responses to faster approaching or receding individuals. Using the definitionṽ ji =ṙ ji , one obtains the relative velocitỹ as a function of the velocity angles ϕ i , ϕ j and the angle α ji of the distance vector r ji , with r ji = r ji (cos α ji , sin α ji ). From (2.5), one finds the two different spatial regions of approaching and receding particles, respectively, where the relative velocity has a different sign for fixed ϕ i and ϕ j in the interval ϕ i ≤ ϕ j < ϕ i + 2π: Figure 1. Visualization of the spatial regions for approaching and receding selfpropelled particles for a binary interaction of the ith particle with velocity vector v i (heading angle ϕ i ), and a neighbouring particle j with velocity v j (heading angle ϕ j ) within its sensory range l s . The center of the circle corresponds to the position of the focal particle i. The dotted line, determined by the mean angle (ϕ i + ϕ j )/2, represents the border between the two distinct spatial regions (half-discs) corresponding to approaching and moving away of individual j. If the jth particle is located above the dotted line in the red region, the two particles are coming closer (approaching). Otherwise, if the neighbouring particle is located in the blue region, the two particles move away from each other (as shown in this example).
Some typical spatial snapshots obtained from individual-based simulations of (2.1) are shown in figure 2. A schematic visualisation of the turning behavior of interacting agents in the different regimes is shown in Appendix A. Please note that in previous publications [27,28], the "effective alignment" regime was referred to as "escape and pursuit", whereas "effective anti-alignment" was labelled "head on head". These previous labels had their origin in the context of heterogeneous agents and the respective response of individual agents to neighbors irrespective on their behavior. Furthermore, an asymmetric alignment response for µ m > 0, µ a < 0, Spatial snapshots of the microscopic model for different interaction strengths after the system relaxes towards a steady-state: (A) disordered, homogeneous state for µ a = µ m = −0.6 (pure repulsion), (B) diffuse collectively moving bands for µ a = −0.6, µ m = 0 (repulsion from approaching individuals), (C) collectively moving bands for µ a = −0.6, µ m = 0.6 (effective alignment), (D) dense, collectively moving cluster for µ m = 0.6, µ a = 0 (attraction to individuals moving away), (E) cluster state for µ m = µ a = 0.6 (pure attraction), (F) disordered, homogeneous state for µ m = −0.6, µ a = 0.6 (effective anti-alignment). The velocity vectors of individual particles are coloured corresponding to their polar direction of motion according to the inset of (A). The inset in (D) indicates the positions of the snapshots in the (µ m , µ a )-parameter plane. Other parameters: N = 4000, L = 40, D ϕ = 0.02, l s = 1, l c = 0.2, µ r = 5, s 0 = 0. 25. yields in the present model the same qualitative dependence of the observed spatial patterns on the interaction strengths as in the original "escape and pursuit" model [33]. Nevertheless, for self-propelled agents with constant speed, the above labels appear more appropriate.
One could argue, that in the effective alignment regime the macroscopic dynamics is essentially identical to other models, as for example the minimal Vicsek-model. However, for strong attraction to particles moving away in comparison to the repulsion to approaching particles, we observe the emergence of dense, collectively moving clusters, which do not appear in simple alignment models.
Our model corresponds to the ones studied in [27,28], if the distance dependent interaction strengths µ r (r ji ) and µ a,m (r ji ) are assumed to be constant and the forces are rescaled by the number of particles within the interaction area of the focal particle.

Mean-Field Theory -Derivation and Analysis of a Nonlinear Fokker-Planck Equation
In this section, we derive a kinetic description for the above individual-based model (2.1). For this purpose, we introduce the N-particle probability density function (PDF) which determines the probability to find particle i at position r i moving into the direction ϕ i (i = 1, 2, . . . , N) at time t. It is normalized with respect to integration over all positions and angles: In agreement with (2.1), we can write down the Fokker-Planck equation (FPE) for the dynamics of the PDF as follows.
From the linear FPE above one can derive an evolution equation for the marginal PDF by integrating (3.3) over the degrees of freedom of particles j = i. In what follows, we assume that correlations between particles can be neglected. Therefore, the N-particle PDF factorizes, i.e.
By this means, one obtains an effective one-particle description where the force F ϕ is given by the following integral over P (r i , ϕ i , t): The factor N −1 in (3.6) arises, because (3.3) was integrated over the positions and angles of N − 1 identical particles. In other words, the focal particle can interact with N − 1 identical neighbours. For the same reason, the particle index i is omitted henceforward. The FPE (3.6), which was derived from the linear FPE (3.3) governing the dynamics of P N , is nonlinear since the force terms F ϕ depend on P (r, ϕ, t). In this sense, assumption (3.5) is a mean-field approximation: a single particle is affected by a force due to its own PDF (3.7).
By introducing the one-particle density function p(r, ϕ, t) = NP (r, ϕ, t) we can eliminate the factor (N − 1) ≈ N from (3.6). Accordingly, p(r, ϕ, t) is interpreted as particle density obeying the nonlinear Fokker-Planck equation First, we assume a spatially homogeneous situation where p(r, ϕ, t) = p(ϕ, t) holds. In this case, (3.8) is reduced to where the dimensionless coupling parameter κ = πρ 0 s 2 0 2D ϕ ls lc dr ji r ji (µ m (r ji ) − µ a (r ji )) (3.10) is introduced, where ρ 0 = N/L 2 is the homogeneous particle density, with L being the linear spatial dimension of the two-dimensional square-shaped system. One obtains the same Fokker-Planck equation as (3.9) for self-propelled particles interacting via a velocity-alignment mechanism or globally coupled Kuramoto oscillators [36,37]. Hence, the selective attraction-repulsion interaction reduces to effective velocityalignment as considered in [38] (continuum time form of the original Vicsek model [12]), if spatial inhomogeneities are neglected. Furthermore, the effective alignment strength is proportional to µ m −µ a and may be negative, leading to anti-alignment as a consequence of repulsion from particles moving away and attraction to particles coming closer, respectively (see also figure A1(D)).
The spatially homogeneous, time-independent particle density is a solution to both nonlinear Fokker-Planck equations (3.8) and (3.9), as can be easily shown by inserting (3.11) into (3.9). I ν (x) denotes the modified Bessel-function of the first kind. Without loss of generality, we choose the frame of reference such that the direction of collective motion ϕ 0 defines the x-coordinate, i.e. ϕ 0 = 0. From the definition of the polar order parameter we can determine the order parameter by solving the transcendental equation Equation (3.13) is always fulfilled for Φ = 0. In that case, (3.11) yields the uniform distribution (spatially homogeneous, disordered state), i.e.
If κ > 2 (high effective alignment strength µ m − µ a and low noise D ϕ , respectively), there is a second nontrivial solution Φ > 0 to (3.13) describing a spatially homogeneous state with polar order (collective motion; swarming phase). In this parameter region, the distribution p (Φ) describing polar order (Φ > 0) is stable with respect to spatially homogeneous perturbations δp = δp(ϕ, t), whereas the homogeneous distribution p (0) is stable for κ < 2 [36,37]. The critical line κ = 2 for the order-disorder transition is translated to microscopic model parameters as follows: The above is a generalization of the previous result obtained in [28]. Since the orderdisorder transition described by (3.13) is continuous, the right hand side of (3.13) can be expanded in a Taylor series so that the dependence Φ = Φ(κ) is known close to the critical point κ c = 2 analytically: The stability analysis of the spatially homogeneous solutions p (Φ) with respect to arbitrary perturbations δp = δp(r, ϕ, t) is carried out in the next section.

Stability Analysis in Fourier Space
In this section it is shown, how the stability of solutions of the nonlinear Fokker-Planck equation (3.8) can be analyzed in Fourier space. A similar approach was used by Chou et al in order to achieve a kinetic theory for self-propelled particles with metric-free interactions [29]. In order to analyze the Fokker-Planck equation (3.8) analytically and especially to solve the integral (3.7), it is convenient to work in Fourier space with respect to the spatial coordinates r and the angular variable ϕ. The dynamics of the Fourier coefficientsĝ where the Fourier coefficients of the force (3.7) are defined aŝ Obviously, the Fokker-Planck equation turns into an infinite hierarchy of equations (3.18) in Fourier space.
For the stability analysis, the Fourier coefficientsK n (k, t) of the force (3.7) are required. Due to the fact, that the force F ϕ (r, ϕ, t) depends on the integral over p(r, ϕ, t) (3.7), the Fourier coefficientsK n (k, t) of the force can be written as a linear combination of the Fourier coefficientsĝ r (k, t) using an infinite dimensional matrixK n,r (k). Details on the calculation of the matrix elements are given in Appendix B. Consequently, the dynamics of the Fourier coefficientŝ n (k, t) be a solution of the equation above. The dynamics of a small perturbation δĝ n (k, t) =ĝ n (k, t) −ĝ The stability analysis of the solutions requires the corresponding Fourier coefficientŝ Inserting those in (3.22) yields the following linearized ordinary differential equations with the matrix The eigenvalues σ of the matrix above determine the stability of p (Φ) . In order to calculate the eigenvalues numerically as a function of the wavevector k (dispersion relation σ(k)), it is necessary to find an appropriate closure of the infinite dimensional system of linear equations. Here we assume, that a criticalñ > 0 exists, such that δĝ n = 0 for |n| >ñ holds approximatively (cf. [29]). That assumption is reasonable, because the nth Fourier coefficient is strongly damped, since For the numerical analysis,ñ = 50 is assumed. Using this approximation, the stability of the spatially homogeneous states can be studied. Please note, that only two approximations were used: First, correlations between particles are neglected (3.5) (mean-field); Second, the truncation of the hierarchy of equations (3.25). In principle, it is possible to analyze the stability for arbitrary wavevectors k, notably, there is no restriction to small wavenumbers k. The procedure described above implies arbitrary directions of perturbations (with respect to the direction of collective motion) of the homogeneous distribution p (Φ) . However, the wavevector k = (k, 0) T is chosen parallel with respect to the direction of collective motion. The illustrated method allows to explore the structure of the phase space of the dynamical system. In particular, we are interested in the stability of the solutions p (Φ) in the (µ m , µ a )-parameter space. For κ < 2, the stability of The kinetic approach reveals instabilities nearby the points (A) and (B). In these parameter regions, an interval of wavenumbers becomes unstable (cf. figure 4). In (A), nematic filaments are found in the individual based simulations (see section 5). However, a long-wavelength instability of the spatially homogeneous, disordered state p (0) is found in (C). According to the individual based simulations, a clustering phase emerges nearby the point (C) (cf. figure 2E for a snapshot of a related simulation).
Close to the order-disorder transition line in the regime of collective motion (κ > 2), a long-wavelength instability emerges as can be seen in figure 4(D). According to figures 4(D)-(F), two hydrodynamic modes exist, which satisfy σ(k → 0) → 0, corresponding to transversal and longitudinal perturbations with respect to the direction of collective motion. The destabilization of the spatially homogeneous, ordered state p (Φ) is always determined by these hydrodynamic modes, leading to long-wavelength instabilities.
Furthermore, there exists a region in the parameter space, where p (Φ) is stable against spatially dependent perturbations. According to figure 3, this region is approximately bounded by the critical line κ = 2 (3.15) and the secondary diagonal µ m = −µ a , which is equally confirmed by the individual based simulation (see section 5).
In summary, it can be stated that the kinetic approach enables the prediction of the structure of the phase space as well as the dispersion relation for any microscopic model parameters. However, the analysis involves numerical methods. Unfortunately, it is not possible to deduce further analytical criteria for the stability of the solutions considered. For that purpose, additional assumptions are necessary, namely the restriction to the lowest Fourier coefficients and small wavenumbers, as shown in section 4. The hydrodynamic theory supplies analytical criteria for the stability of the solutions considered, such as critical microscopic model parameters.
In this context, we would like to emphasize one difficulty which occurs if the focus is put on the dynamics of the lowest Fourier coefficients on large length scales, i.e. small wavenumbers k: The matrix elementsK n,r (k) involve Bessel functions of the first kind of k (see Appendix B and [29]), which are oscillating functions with alternating Taylor coefficients [30]. The approximation will only be valid in the vicinity of k = 0. For large k, J ν,N (kx) tends to either plus or minus infinity, but never to zero. This issue will be important for the derivation of hydrodynamic equations in the next section and is restricting their validity to small wavenumbers k and large length scales, respectively.

Expansion of the Fokker-Planck Equation in the Angular Fourier Domain
In this section, a hydrodynamic description of the many-particle system is derived directly from the one-particle FPE (3.8). For this purpose, it is convenient to work in Fourier space with respect to the angular variable ϕ. Both, the particle density p(r, ϕ, t) and the force F ϕ (r, ϕ, t) (3.7) are 2π-periodic functions in ϕ. Hence, they can be expanded in a Fourier series as follows: The Fokker-Planck equation in Fourier space is obtained by multiplying (3.8) with e inϕ and integration over ϕ ∈ [0, 2π]: In (4.2), the complex derivative ∇ = 0.5 (∂ x + i∂ y ) was introduced. The Fourier coefficientsf n (r, t) = obey the symmetryf * n (r, t) =f −n (r, t), since p(r, ϕ, t) is real-valued. The lowest Fourier coefficients are related to macroscopic physical quantities, namely the marginal particle density ρ(r, t) and also the momentum field w(r, t) = (w x , w y ) via ρ(r, t) =f 0 (r, t) , The second Fourier coefficients are related to the nematic order parameter | e 2iϕ | and the symmetric temperature tensor, as defined in [23,25,28], which is a measure for the width of the velocity distribution (see Appendix D for details). The dynamics of the temperature tensor is not derived explicitly in this context. In the following, the dynamics of the density and the momentum field is deduced from (4.2). Therefore, the Fourier coefficients of the forceQ n (r, t) are approximatively calculated by expanding p(r + r ji , ϕ, t) into a multidimensional Taylor series for small |r ji | [31], substituting (4.1a) into (3.7) and evaluating the remaining integral. Here, we consider all terms up to the second order of the Taylor expansion. This approximation holds for spatially slowly varying densities, i.e. large system size compared to the interaction radius l s . Again, the Fokker-Planck equation turns into an infinite hierarchy of equations (4.2) in Fourier space. An appropriate closure scheme is used in order to consider only the first Fourier coefficients. Close to the order-disorder transition, the degree of polar order is assumed to be small, i.e. |w| /(s 0 ρ) ∝ ǫ, where ǫ is a small number and the particle density is locally close to the homogeneous distribution. Furthermore, the dynamics of the Fourier coefficients (4.2) suggests, that f 1 is larger than f n with |n| > 1, because of the damping term in (4.2) proportional to n 2 , leading to the scaling relations In [20] it is argued, that the scaling ansatz of the temporal and spatial derivatives reflects the propagating nature of the system. This closure scheme described above was already used in [17,20,22,31] and in [21] in a modified way for nematic particles. A system of three nonlinear partial differential equations is obtained by keeping all terms up to the order ǫ 3 .
The following coefficients are introduced: 7c) (4.7e) Equation (4.6a) is simply the continuity equation representing the conservation of the particle number N.

Hydrodynamic Limit
Another simplification of the system (4.6) is reached by adiabatically eliminatingf 2 , i.e. ∂ tf2 ≈ 0, wherê follows. This approximation is only valid if the first Fourier mode relaxes considerably slower thanf 2 . The assumption of a time scale separation between the first and second Fourier mode can be only be justified if |ξ 1 ρ − ξ 5 | ≪ 4ξ 5 , which is equivalent to −6 ≪ κ ≪ 10. Furthermore ξ 5 ∝ D ϕ > 0 has to hold. Please note, that all terms of the order ǫ 3 were dropped in (4.8) according to the scaling relations (4.5).
By identifying Fourier coefficients with physical quantities (4.4), familiar hydrodynamic equations are obtained.
The transport coefficients {λ i , η i } as functions of ξ i (4.7) read The hydrodynamic equations (4.9) are comparable in structure to the theory of Toner and Tu [13,14].
In the theory of Toner and Tu, the coefficients η 2 and η 3 are assumed to be nonzero phenomenological parameters due to the absence of Galilean invariance [14]. Indeed, in our model, η 2 and η 3 can take both, negative and positive values. Likewise, λ 2 can be positive and negative: A positive λ 2 -term leads to a flow into the direction of higher densities and consequently to clustering, as shown below. The damping coefficient η 1 is positive for all microscopic parameters, consistently. The sign of the different transport coefficients is discussed in Appendix C in detail.
The stability of the hydrodynamic equations (4.9) requires that λ 4 and λ 5 are greater than zero. In this region, the hydrodynamic equations describe collectively moving bands and the theory of Toner and Tu applies. In contrast to the original Vicsek model, in our model the coefficients λ 4 and λ 5 may change signs. As shown in [3], a negative λ 4 may lead to an instability of spatially homogeneous states. In this case, higher order derivatives are needed in (4.9b) in order to guarantee the stability of the dynamics + . Moreover, the time scale seperation of the first and second Fourier mode is only justified for ξ 1 ρ ≈ ξ 5 , i.e. λ 1 ≈ 0. For the remaining analysis of the hydrodynamic equations, we will assume the stability and validity of our hydrodynamic theory (4.9).
Let us begin the analysis of (4.9) by neglecting the spatial derivatives and analyzing the fixed points of dw dt = λ 1 − η 1 |w| 2 w . (4.11) That is the standard symmetry-breaking term present in all models of collective motion [3, 14, 18-20, 22, 31]. Apparently, two spatially homogeneous solutions exist: First, w 0 = 0 corresponding to a disordered phase and vanishing center of mass velocity. Second, where e denotes an arbitrary unit vector reflecting the isotropy of the system. Without loss of generality we set e = e x . For λ 1 < 0, only the fixed point w 0 exists and is stable against spatially homogeneous perturbations. For λ 1 ≥ 0, w 1 is a second fixed point corresponding to an ordered phase with nonzero mean-speed (swarming phase). The supercritical pitchfork bifurcation in λ 1 = 0, i.e. introducing the dimensionless coupling strength κ (3.10) again, κ = πρ 0 s 2 0 2D ϕ ls lc dr ji r ji (µ m (r ji ) − µ a (r ji )) , (4.14) the polar order parameter Φ H = |w 1 | /(s 0 ρ 0 ) as described by the hydrodynamic theory can be written as follows * : Close to the critical point κ c = 2, one recovers the correct scaling behaviour Φ H ≃ √ κ − 2 (3.16). However, for large κ the polar order parameter Φ H tends to zero according to (4.15) as shown in figure 5. This is a consequence of the approximations made, namely the elimination of higher Fourier modes. If one allows a maximal relative error (Φ − Φ H )/Φ of 5%, an upper bound for κ is found numerically: κ ≤ 2.617. Equation (4.13) defines the transition line in the (µ m ,µ a )-parameter space. Interestingly, only the integrated interaction strengths are important, the exact functional dependence on the distance between two particles µ m,a = µ m,a (r ji ) does not matter. Furthermore, collective motion is possible in both situations, pure attraction to particles moving away (µ m > 0, µ a = 0) and pure repulsion from approaching particles (µ m = 0, µ a < 0), respectively. Since the critical coupling parameter κ c = 2 is positive, polar order cannot be found in the anti-alignment regime.
Provided all eigenvalues are less than zero, the disordered spatially homogeneous state is stable. The first eigenvalue σ 1 is equal to λ 1 for k = 0, i.e. the stability of w 0 is determined by the sign of λ 1 . That is the destabilization of the homogeneous state as discussed before. For λ 4 > 0, the first eigenvalue σ 1 is monotonically decreasing. However, the eigenvalue σ 1 is monotonically increasing for λ 4 < 0. In this case, the hydrodynamic equations (4.9) loose their validity. This unphysical behaviour is due to the restriction to second derivatives when the integral (3.7) was calculated. Furthermore, it is problematically that the Taylor coefficients have alternating signs so that Taylor polynomials converge slowly as argued at the end of section 3.2. To predict the behaviour of the system in that case, one needs to consider higher order terms. Unfortunately, the number of terms and Fourier coefficients that has to be considered to be consistent with the scaling ansatz (4.5), grows rapidly. Nevertheless, it is possible to predict the stability of w 0 by analyzing those eigenvalues, which tend to zero for small wavenumbers (hydrodynamic modes [14]). For small wavenumbers, the eigenvalue σ 2 (σ 2 > σ 3 ) is expanded in a Taylor series: Close to the order-disorder transition line, λ 1 is approximately zero. Therefore it is sufficient to keep the leading order terms in 1/λ 1 [18]: Suppose, w 0 is stable against spatially homogeneous perturbations, i.e. λ 1 < 0. In that case, a long-wavelength instability emerges for positive λ 2 . Considering (4.9), this instability is reasonable: for λ 2 > 0, inhomogeneities in the particle density will lead to a flow in the direction of the density gradients while for λ 2 < 0 density inhomogeneities decay due to inverse flows. The amplification of density gradients will lead to an agglomeration of particles (cf. figure 2E). This hypothesis is supported by numerical simulations of the microscopic model, as shown in the next section and in [28]. A similar instability due to the gradient term in (4.9) is also known from other active matter systems [22]. Whereas in [22] the instability is due to the density-dependent motility of the particles, in our case the instability is referable to the selective interaction. The critical line is given by λ 2 = 0, i.e.  The corresponding analysis of the stability of the spatially homogeneous ordered state for |w 1 | > 0 using the full hydrodynamic equations is much more complicated. It is done in great detail for the original Vicsek model in [18] for structurally identical hydrodynamic equations. In the following, the key points of the analysis are summarized. Moreover, an additional instability is found that is not present in the Vicsek model because in our model the transport coefficients λ 2 , η 2 and η 3 may change their algebraic signs. Inserting the ansatz ρ(r, t) = ρ 0 + δρ(r, t) , into (4.9) and linearization yields Again, inserting (4.18) yields the growth rate σ dependent on the wavevector k. We restrict ourselves to longitudinal perturbations, meaning that the wavevector k, the perturbation δw and the direction of collective motion w 0 are pointing in the same direction, because we are interested in the emergence of collectively moving bands (cf. figure 2). The growth rate reads where the coefficients were introduced. The real part of the largest eigenvalue reads as follows: Close to the order-disorder transition, λ 1 is small, with the result that it is sufficient to keep the lowest order in 1/λ 1 . The following expression is obtained for small wavenumbers [18]: That implies, that for λ 1 0, a long-wavelength instability emerges. By expanding (4.27) in a Taylor series for small k (not assuming that λ 1 ≃ 0), one obtains In general, an instability appears, if the first Taylor coefficient becomes positive, because σ + (k → 0) → 0. Unfortunately, it is not illuminating to derive the critical line µ m = f (µ a ), where the relevant Taylor coefficient is zero, analytically, because the resulting expressions are quite complicated. The numerical analysis reveals that the spatially homogeneous, ordered state is unstable beyond a hyperbola which is bounded by the two asymptotes (see also figure 6) . (4.30b) The instability of w 1 for is due to the selective interaction, therefore it is not present in the original Vicsek model. The two asymptotes (4.30) were derived using the approximation that the interaction strengths do not depend on the distance. The second line (4.30b) is equal to the critical line λ 1 = 0 (4.13).
As already mentioned, negative values of λ 4 and λ 5 may lead to additional instabilities, which go beyond the scope of this work, because the hydrodynamic equations are restricted to second derivatives. However, neither (4.21), (4.28) nor (4.29) depend on λ 4 and λ 5 . Hence, the spatially homogeneous states and instabilities described above are always present in our model.
The predictions of the hydrodynamic theory, derived in this section, are compared to the predictions of the stability analysis using the kinetic description in section 3.2 in figure 6(b). The long-wavelength instability of w 0 , leading to clustering for λ 2 > 0 as well as the long-wavelength instability close to the critical point (κ = 2, λ 1 = 0) are confirmed by the kinetic approach.
Within the range of validity of the hydrodynamic theory, no contradictions are found with the kinetic theory. However, figure 6(b) indicates that the range of validity of the hydrodynamic theory does not cover all relevant parts of the parameter space. In particular, the destabilization of the ordered state w 1 is not sufficiently well described by the hydrodynamic equations. Please note, that the hydrodynamic theory is only valid for λ 4 > 0 (guarantees the stability for large wavenumbers) and 0 ≤ κ ≤ 2.617, where the lower bound is equivalent to λ 5 > 0 and the upper bound ensures small deviations Φ − Φ H . Due to the last limitation, it is impossible to study the dynamics of the system for high order parameters, i.e. far away from the critical point.
However, one can show numerically that the dispersion relations of the hydrodynamic theory (4.9) and the full system (3.18) coincide for small wavenumbers in the parameter range, where the hydrodynamic equations are valid. The hydrodynamic theory yields only three eigenvalues, which coincide with three eigenvalues given by the kinetic theory in the limit k → 0. For large k, deviations occur due to the restriction to second derivatives (cf. scaling ansatz (4.5)).
Nevertheless, the hydrodynamic theory yields important insights on the macroscopic behaviour of the system. The symmetry breaking, i.e. the existence of a collective motion mode is predicted. Moreover, the stability of the homogeneous solutions to the hydrodynamic equations can be analyzed analytically and it allows to obtain analytical results on the order-disorder transition line as well as the occurrence of clustering. Thus, the structure of the phase space is roughly estimated by the hydrodynamic theory. New instabilities are found due to changes in the algebraic signs of the coefficients λ 2 and η 2 . Furthermore, the structure of the theory suggests, that all predictions of the theory by Toner and Tu [14] are expected to arise in our system, such as giant number fluctuations. Within the range of validity of the hydrodynamic theory, the predictions of the kinetic theory, which is based on the mean-field assumption only (see section 3 for details), are in agreement with the hydrodynamic theory, whose range of validity is known quantitatively.

Large-System Limit
The hydrodynamic theory is simplified further, if the limit of large systems, i.e. large particle densities (N → ∞) and small interaction regions (l s → 0) is considered, such that the number of particle within the interaction area is kept fixed: ρ l 2 s ≈ const. This approximation corresponds to a zeroth order Taylor expansion of the force (3.7): p(r + r ji , ϕ j , t) ≈ p(r, ϕ j , t). This approximation of "ultra-local" coupling is used in [18,20] within a Boltzmann approach, where the collision integral is an integral over the particles orientation but not over relative distances. In that case, the hydrodynamic equations read ∂ρ ∂t = −∇ · w , (4.32a) Since D ϕ is always positive, the dynamics is stable for all microscopic parameters. The coefficient ξ 1 in front of the hydrodynamic term (w · ∇) w, which corresponds to the convective derivative in the equilibrium case, is positive for µ m > µ a and is negative for µ m < µ a . The same holds true for the other "hydrodynamic terms" ∇ |w| 2 and w (∇ · w). Microscopically, the transition from positive ξ 1 to negative ξ 1 means, that the selective interaction turns from an effective alignment to an effective anti-aligning interaction, meaning that a focal particle aligns its velocity in an anti-parallel fashion with respect to the local mean-velocity (see figure A1). Similar conclusion could possibly be drawn from [22,38], if one allows the alignment strength to be negative. In that part of the parameter space, the nematic structures emerge in our model according to the individual based simulations discussed in the next section.

Comparison with Numerical Simulations
In order to analyze the stability of the disordered, homogeneous state, we have performed systematic numerical simulations of the individual-based model (2.1). The degree of collective motion was measured using the time averaged polar order parameter Here · t represents a temporal average. Φ t = 1 corresponds to perfect orientational order with all agents moving in the same direction, whereas a vanishing Φ t corresponds to a completely disordered system. Please note that Φ t = 0 can only be observed in the thermodynamic limit (N → ∞). In a finite, disordered system we will measure a small, but finite Φ t 0 due to finite size fluctuations of the order 1/ √ N. In order to measure the deviations from a spatially homogeneous state, we have subdivided the simulation domain into square cells of size l s × l s set by the sensory range. We used this spatial subdivision to calculate the spatial entropy function where the summation occurs over all occupied cells of the grid with n j ≥ 1 (n j is the number of particles in the jths cell). This allows us to define the following spatial order with S max being the maximal value of the spatial entropy corresponding to a homogeneous distribution of particles. Ψ t = 0 corresponds to a perfectly disordered state, whereas Ψ t > 0 indicates a spatially inhomogeneous distribution of particles (clusters, bands), where Ψ t = 1 corresponds to the extreme situation where all particles are located in a single cell. In order to test the stability of the disordered state, we have averaged the two order parameters over a time interval ∆t = 1000 after an initial time t ini = 1000. Please note, that for certain parameter values it is possible that the system has not reached a stationary state after t = 1000. However, this initial time is sufficient to account for deviations from the homogeneous disordered state, which we are interested in. Accordingly, we are not interested in the actual stationary values of Φ t and Ψ t . Thus, we set the upper limit of the colour bar for both order parameters in figures 7 and 8 to 0.2. In consequence, all regions of parameter space, where Φ t ≥ 0.2 or Ψ t ≥ 0.2 will appear white in figures 7 & 8.
The stochastic differential equations of the microscopic model were integrated with periodic boundary conditions using the stochastic version of the Euler algorithm with a numerical time step dt = 0.01. This time steps is two orders of magnitude smaller than the average time scale of turning of individuals due to binary interactions for the range of interaction strengths in the simulations. Additional sample runs with smaller time steps did not show any significant differences in the simulation results. The initial condition for all simulations was the disordered, spatially homogeneous state.
For simplicity, we consider only constant interaction strengths (independent of the distance): µ m,a (r ji ) = µ m,a = const. Here, we focus on the analysis of the stability of the disordered, homogeneous solution with respect to variations of the interaction parameters µ a and µ m . The analytical predictions of the hydrodynamic theory for the onset of instability of the disordered solution are represented by two intersecting critical lines in the (µ a , µ m )-plane perpendicular to each other (figure 7). The first one (µ a ∼ µ m ) corresponds to the orientational instability and onset of collective motion, cf. (4.13) and line (1) in figure 7, whereas the second one (µ a ∼ −µ m ) corresponds to the density instability associated with structure formation due to effective attraction between particles, cf. (4.22) and line (2) in figure 7. The homogeneous, disordered solution is predicted to be linearly stable only "below" both critical lines.
The emergence of collective motion as well as the destabilization of the homogeneous density distribution is mostly confirmed by the numerical individual-based simulations. In specific parameter regions, deviations from the hydrodynamic theory are observed.
In particular, we observe a high degree of collective motion in the effective-alignment region and strong deviations from the homogeneous spatial distribution of particles in the pure attraction regime as well as in the effective-alignment regime (in particular in the regime (µ m , µ a ) ∈ {(µ m , µ a ) : |µ m | > |µ a |, µ m > 0}). This also confirms our previous results of comprehensive numerical investigations of related individual-based models [28,33].
Disagreements between simulations and the theoretical prediction (non-vanishing orientational order and/or clustering in simulations below the two critical lines) appear predominantly close to the intersection of the two critical lines and are stronger at low densities. They might be associated with the mean-field assumptions used in order to derive the hydrodynamic theory. On the one hand, at low densities the assumption of a continuous density of neighbours is strongly violated. On the other hand, correlations between particles may play an important role in the respective parameter range (likewise for high densities). Thus, the factorization of the N-particle PDF into a product of oneparticle PDFs (3.5) leads to a questionable approximation in that case.
Another possible explanation for disagreement between theory and simulation is a breakdown of the homogeneous, disordered solution due to finite amplitude instabilities at parameters where this solution is still linearly stable.
For weak (or vanishing) short-ranged repulsion (l s ≫ l c or µ r ≪ |µ m,a |), we observe inhomogeneous states without polar order far in the effective anti-alignment regime (µ m < 0, µ a > 0), clearly below the critical line predicted by the hydrodynamic theory as shown in figure 8. This instability was missed in the previous study of the model [28]. We were able to confirm it, using our kinetic approach (see section 3) by positive eigenvalues of the matrix determining the stability of the disordered solution at the respective parameter values. A close inspection of the dynamics of the individual-based model in this parameter region reveals the emergence of dense nematic filaments with particles moving in an anti-parallel fashion within the filament, whereby approximately half of the particles moves in either direction along the filament (see figure 9).
The hydrodynamic theory derived in section 4 is not able to account for this kind of structures as it considers only the (polar) momentum field and not the nematic director field as a coarse-grained variable. An extension of the hydrodynamic theory by including an equation for the nematic director field may allow for identifying instabilities due to the onset of nematic order, however it goes beyond the scope of this work. We refer the reader to some recent theoretical works on active nematics [21,34,35,39,40].

Discussion
The self-propelled particle model with selective attraction-repulsion interaction shows a rich dynamical behavior. In the "effective alignment" regime, characterized by repulsion from approaching agents and attraction to moving away agents, the model behavior is closely related to the well-known Vicsek model: For repulsion from approaching particles equal to, or stronger than, the attraction to particles moving away, we typically observe the homogeneous flocking phase with giant number fluctuations predicted by Toner and Tu. However, in the opposite case, if the attraction to particles moving away starts to dominate, we observe strong clustering and effective phase separation, which is not observed in minimal models with constant speed and velocity alignment. Furthermore, by continuously varying model parameters, we can observe other phases such as unpolarized clusters for pure attraction or dense filamentous nematic structures in the effective anti-alignment regime.
In order to obtain a fundamental understanding of the phase diagram of the model, we derived first a kinetic description of the system based on the Fourier transform of the probability density function. The corresponding system of equations for the successive Fourier modes can be used for efficient numerical analysis of the linear stability of solutions of the nonlinear Fokker-Planck equation. We have analyzed the stability of the spatially homogeneous solutions. In addition, we have shown that the integration over the social forces yields Bessel functions of the first kind, which enter the matrix elements of the corresponding linearized system of differential equations in Fourier space. Due to the alternating Taylor coefficients of the corresponding Bessel functions, a closure approximation, corresponding to a finite order expansion, immediately leads to unphysical divergences at large wavenumbers k.
Furthermore, we have derived a hydrodynamic theory by truncating the "smallwavenumber-expansion" at the second order. The resulting hydrodynamic equations are in agreement with the generic Toner and Tu theory of active matter. Our work establishes a direct link between the microscopic parameters of the individualbased model and the macroscopic parameters governing the behaviour of the coarsegrained hydrodynamic variables (density and momentum fields). Interestingly, the hydrodynamic parameters, as for example η 3 , which relates to a splay elasticity in corresponding equilibrium systems, or λ 4 and λ 5 which govern the relaxation of splay and bend fluctuations [?], may change their sign in dependence on microscopic model parameters. Using the hydrodynamic theory, we can track down such sign changes to corresponding microscopic dynamics. Considering a limit of large system sizes, we reveal the importance of the change of sign of σ 1 , which indicates the transition from an overall aligning effect of the social interaction to the opposite case where the interaction tend to anti-align interacting particles. Finally, a comparison between the eigenvalues obtained from hydrodynamic theory and the kinetic description, which takes higher orders into account, allows us to assess the validity of the hydrodynamic theory.
We performed extensive simulations of the microscopic model based on stochastic differential equations focusing on the (µ m , µ a )-plane in the parameter space. The critical lines obtained from the hydrodynamic theory, where the disordered, spatially homogeneous solution becomes unstable, show good agreement with the numerical results at sufficiently high densities. However, at certain parameters clear deviations appear, as for example in the vicinity of the crossing of the two critical lines at low densities.
This leads us to the important question of the validity of the approximations made during the coarse-graining. One common assumption is that of molecular chaos, which allows to factorize the N-particle PDF into a product of N one-particle PDFs. In models with collision-like interactions, the approximation is supposed to work best at low densities, where the mean-free path of particles between interactions is large [17,19]. However, in our case, we observe large deviations at low densities. This contradicting effect may be due to an approximation required to evaluate the integrals over the social forces. It relies on a Taylor expansion of the one-particle density function around the position of a focal individual. However at low densities, the interaction range l s is not the suitable coarse-graining scale as assumed implicitly in the formulation of the oneparticle Fokker-Planck equation. Similar arguments have been put forward also in [22]. The impact of individual (angular) noise on the mean-field assumption may also be not straight forward. Intuitively, one would argue that uncorrelated individual noise terms always decrease correlations between interacting particles. However, for self-propelled particles, angular noise leads to a stronger localization of particles [41,42], thus in principle also to a prolonged interaction between neighbours, which may in principle enhance multi-particle correlations.
Furthermore, the systematic comparison of the prediction of the kinetic theory with the predictions of the hydrodynamic theory and numerical simulations revealed an unexpected additional instability. It corresponds to the emergence of dense filamentous structures with nematic order. Onset of nematic order is linked to the dynamics of the second Fourier amplitude, which was adiabatically eliminated in order to derive the hydrodynamic theory. Thus, the presented hydrodynamic theory cannot account for this instability, but it can be well traced by numerical evaluation of the linearized kinetic equations derived in section 3.
So far, hydrodynamic equations of active matter were derived directly from minimal microscopic models of self-propelled particles with velocity-alignment. Here, we show that it is also possible to derive such equations for a more complex model of selfpropelled particles with selective attraction-repulsion interaction and establish a direct link between the microscopic and macroscopic level of description. The model exhibits a large variety of different phases and we believe, it might be not only of interest from the biological point of view, as an alternative to models including explicit alignment of individual agents, but that it offers also interesting playground for the study of self-organization, pattern formation and phase transitions at far-from-equilibrium conditions. into (3.7) and solving the remaining integral: In order to perform the integration over α and ϕ j , the following identity is used to express the exponential function e −ik·r ji = J 0 (k r ji ) + 2 where k = k (cos χ, sin χ) and r ji = r ji (cos α, sin α). The Bessel functions of the first kind are denoted by J ν (x). Performing the integration yields the force expanded into a Fourier series and as a linear combination (B.1b) of the Fourier coefficientsĝ r (k, t), where one can read off the elements of the infinite dimensional matrixK p,n (k): K p,n (k) = s 0 π 2 2i {(μ m,0 (k) −μ a,0 (k)) (δ p,1 δ n,1 − δ p,−1 δ n,−1 ) + (μ m,2 (k) −μ a,2 (k)) e −2iχ (δ p,−1 δ n,1 − δ p,−2 δ n,0 ) − e 2iχ (δ p,1 δ n,−1 − δ p,2 δ n,0 ) ... + 2π 2μ r (k) δ n,0 e iχ δ p,1 − e −iχ δ p,−1 + 2s 0 ∞ s=0 (−1) s (μ m,2s+1 (k) +μ a,2s+1 (k)) Λ p,s e −i(2s+1)χ δ p+2s+1,n − Γ p,s e i(2s+1)χ δ p−(2s+1),n . (B.5) The following coefficients were introduced: The Bessel functions J ν (x) are oscillating functions with alternating Taylor coefficients. A truncation of the Taylor series (as done in section 4) at a finite order may therefore lead to unphysical behaviour for large wavenumbers k.
The extrema of the hyperbolas are located at The critical line for η 3 = 0 depends on the noise strength, whereby the critical noise strength reads . (C.9) ϕ . (C.10) Please note, that the slope of the lines is bounded by one (D ϕ = 0) and minus one (D ϕ → ∞), respectively. Hence, η 3 is always positive for {(µ a , µ m ) : µ a < 0 , − |µ a | < µ m < |µ a |} and η 3 always negative for {(µ a , µ m ) : µ a > 0 , −µ a < µ m < µ a }. All criteria discussed above are summarized in figure C1. The stability of the hydrodynamic equations (4.9) requires that both λ 4 and λ 5 are greater than zero. This is in principle true in the effective-alignment region (at least for sufficient high noise and low densities), as well as parts of the pure attraction and pure repulsion region. In this parameter region, the hydrodynamic equations describe collectively moving bands.