Role of particle conservation in self-propelled particle systems

Actively propelled particles undergoing dissipative collisions are known to develop a state of spatially distributed coherently moving clusters. For densities larger than a characteristic value clusters grow in time and form a stationary well-ordered state of coherent macroscopic motion. In this work we address two questions: (i) What is the role of the particles' aspect ratio in the context of cluster formation, and does the particle shape affect the system's behavior on hydrodynamic scales? (ii) To what extent does particle conservation influence pattern formation? To answer these questions we suggest a simple kinetic model permitting to depict some of the interaction properties between freely moving particles and particles integrated in clusters. To this end, we introduce two particle species: single and cluster particles. Specifically, we account for coalescence of clusters from single particles, assembly of single particles on existing clusters, collisions between clusters, and cluster disassembly. Coarse-graining our kinetic model, (i) we demonstrate that particle shape (i.e. aspect ratio) shifts the scale of the transition density, but does not impact the instabilities at the ordering threshold. (ii) We show that the validity of particle conservation determines the existence of a longitudinal instability, which tends to amplify density heterogeneities locally, and in turn triggers a wave pattern with wave vectors parallel to the axis of macroscopic order. If the system is in contact with a particle reservoir this instability vanishes due to a compensation of density heterogeneities.

Theoretically, the emergence of collective motion has mostly been studied in the context of particle conserving systems. There are, however, a number of experimental systems, in which the assumption of particle conservation is questionable. In typical gliding assays [4-6, 9, 10], for instance, collective motion of filaments is observed on a two-dimensional "motor carpet" which itself is in contact with a three-dimensional bulk reservoir of filaments. However, the impact of particle conservation on the formation of patterns of collective motion remains largely elusive.
Here, we address the significance of constraints for particle number by highlighting the differences in the collective properties between particle conserving systems and those in contact with a particle reservoir. Our focus will be on the comparison of two archetypical scenarios, which we will refer to as the canonical (particle conserving), and the grand canonical (violating particle conservation) scenario, respectively.
To this end, we will resort to a kinetic approach, which has been set up previously by Aranson et al. [16] to describe pattern formation in a system of interacting microtubules, and which has been extended to the case of self-propelled spheres by Bertin et al. [17,22]. In the following we will extend this description in accordance with a physical picture of collective motion that has been developed over the last decade based on observations in agent-based simulations of locally interacting, particle conserving systems [31,32,37,39]. Among the most pertinent phenomena that have been reported in the context of these studies is the formation of intricate local structures pervading these systems in the vicinity of the ordering transition: Densely packed cohorts of coherently moving particles-subsequently referred to as clusters-incessantly "nucleate" and "evaporate" on local scales, even below threshold, stable longitudinally unstable olutions to the full canonical hydrodynamic equations (28a) to (28c). In the ing we are going to investigate the linear stability of these solutions against waveerturbations, employing the following ansatz: ⌘(x, t) = ⌘ ⇤ + ⌘(x, t), g(x, t) = g 0 + g(x, t), g(x, t) = g 0 e st+iq·x . (44c) e equations above, q denotes the wave vector and s is the growth rate. Inserting nsatz into equations (28a) to (28c), we find for the canonical model: iq cos( ) g x,0 iq sin( ) g y,0 , where ⇢(x, t) = ⇢ 0 e st+iq·x , ⌘(x, t) = ⌘ 0 e st+iq·x , g(x, t) = g 0 e st+iq·x .
In the equations above, q denotes the wave this ansatz into equations (28a) to (28c), we s ⇢ 0 = iq cos( ) g x,0 iq s state (left), the cluster particles density (blue arrows) constitutes the system's macroscopic net momentum g 0 , while some fraction of the system's particles, the single particles (orange dots) exhibit zero net momentum. Spatial perturbations of both density fields lead to two fundamentally different outcomes: (i) In case of a closed system obeying total particle conservation (single particles+cluster particles), termed as the canonical model, the homogeneously polarized state is longitudinally unstable, with a wave vector q parallel to the polarized state g 0 , potentially enforcing a wave-like pattern. (ii) In contrast, open systems turn out to be stable against this kind of density fluctuations.
rendering the system isotropic and homogeneous only in the limit of macroscopic length scales. Individual particles exhibit superdiffusive behavior in this regime, performing quasiballistic "flights" as long as they are part of a cluster, and conventional particle diffusion if they are not. Above threshold, collective motion manifests itself on macroscopic scales in the form of coherently moving and dense bands, which are submersed in an isotropic lowdensity "particle sea". Spatially homogeneous flowing states, in contrast, are observed only well beyond the ordering threshold [31]. Moreover, particle geometry was demonstrated to play an essential role in the context of clustering dynamics, with higher aspect ratios facilitating the formation of clusters of coherently moving particles [37].
In the light of the above, we suggest a simplified modeling framework to incorporate the intricate role of clusters on the ordering behavior, which will be presented in greater technical detail in the following section: Particles interact via binary collisions with a scattering cross section which is explicitly derived as a function of particle shape. Depending on whether a given particle is part of a cluster or not, it will be associated with one of two distinct particle classes, which we will refer to as the class of cluster particles and the class of single particles, respectively. Single particles are "converted" to cluster particles by "condensation" every time a single particle collides with a cluster. Conversely, cluster particles are "converted" back to single particles by an "evaporation" process which we assume to occur at some constant (possibly particle shape dependent [37]) rate. Moreover, in the absence of interactions, cluster particles will be assumed to move ballistically, whereas single particles will be assumed to perform random walks. Taken together, the conversion dynamics and the class-specificity of particle motion provide a simple way to implement the typical superdiffusive behavior of individual particles, that was alluded to above. To assess the importance of particle conservation in the context of pattern formation, we will analyze two variants of this model: Firstly, we study closed systems in which the total number of particles is conserved (canonical scenario) and where, consequently, the denser cluster phase grows at the expense of the single phase. Secondly, we examine open systems in contact with a particle reservoir (grand canonical scenario), where the particle current out of the single phase is compensated so as to retain the density of the isotropic sea of single particles at a constant level; cf. figure 1.
Our work is structured as follows: In II the modeling framework for the canonical and grand canonical model is introduced and the model equations are discussed in detail. The corresponding hydrodynamic equations are derived in III by means of an appropriate truncation scheme in Fourier space. Therein, we also give explicit expressions of the kinetic coefficients as a function of the particles' aspect ratio and velocity, noise level and density for single particles and cluster particles. IV is devoted to the analysis of the homogeneous equations. The dynamic's stationary fixed points are determined and the phase boundary between the isotropic and homogeneous state is calculated. V deals with the implications of the inhomogeneous equations in the framework of a linear stability analysis, which are concluded in VI.

II. COARSE-GRAINED KINETIC MODEL
We consider rod-like particles of length L and diameter d moving in two dimensions with a constant velocity v. A particle's state is determined by its position x and the orientation θ of its velocity vector. To describe the time evolution of the system, we adopt a kinetic approach [16][17][18]22].
On mesoscopic scales, the system's spatio-temporal evolution is then governed by Boltzmann-like equations for the one-particle distribution functions within the classes of single particles and cluster particles, respectively. Interactions enter this description by means of collision integrals. The kernel of these integrals involves both, a measure for the rate of collisions, as well as a "collision rule" implementing a mapping between pre-and post-collisional directions θ and θ of each of the two partaking particles. Here, we are led to consider a simplified model of binary particle interactions, which builds on the distinction between single particles and cluster particles. The details of this model will be described in the following section.

A. Reaction equations
Let S(θ) and C(θ) refer to a particle moving in the direction of θ and being associated with the class of single particles or cluster particles, respectively. In the absence of interactions, single particles are assumed to perform a persistent random walk, which we model as a succession of ballistic straight flights, interspersed by self-diffusion ("tumble") events. These tumble events are assumed to occur at a constant rate λ and reorient the particle's orientation θ by a random amount ϑ 0 : For simplicity we assume ϑ 0 to be Gaussian-distributed, with σ 0 denoting the standard deviation. On times scales much larger than λ −1 , this tumbling behavior can be described as conventional particle diffusion, with the particles' diffusion constant being a function of λ and σ 0 [40].  2: (a) Illustration of two single particle species (light orange) with a pre-collisional relative angle of θ 12 , colliding such that they align collinear to the average angleθ. Both particles become a cluster species after the collision. (b) Right: Illustration of a possible scenario, where a single particle joins a cluster by perfectly aligning to the cluster particles (blue). Left: A particle leaves the cluster by a random change of its direction at a characteristic rate .
When two single particles S(θ 1 ) and S(θ 2 ) collide they are assumed to assemble a cluster, i.e. each of the two particles becomes a cluster particle (see figure 2a): where [42]θ denotes the average of both pre-collisional angles θ 1 and θ 2 , and where ϑ is a random variable which we, again, assume to be Gaussian-distributed: The rate of binary collisions, such as equation (3), are determined by a particle-shape dependent differential scattering cross section, which will be discussed below; see II B and A.
Collisions involving cluster particles are distinct from single particle events. Due to the close spatial proximity of particles within each cluster these collisions correspond to manyparticle interactions. Needless to say, a detailed description of cluster formation and the ensuing particle dynamics represents a highly complex matter, requiring explicit consideration of such many-particle interactions. For simplicity, we will resort to the following simplified interaction picture: We assume that (binary) collisions between single particles and cluster particles lead to a condensation process during which the single particle aligns to the cluster particle without changing the direction of the cluster as a whole: Eq. (6) thus captures the net effect of collisions between single particles and cluster particles, during which multiple collisions, involving neighboring particles belonging to the same cluster, stabilize the cluster's direction; cf. figure 2 b for an illustration.
Collisions among cluster particles is an even more intricate process, since they actually depend on size and shape of both colliding clusters, and in general involve multi-particle interactions. In the framework of a Boltzmann-like description, correlations in the particle distribution are neglected and only binary interactions are considered. The frequency of interactions are determined by a geometrical construction called the "Boltzmann cylinder", assuming that particle positions are homogeneously distributed on local scales. With regard to many-particle interactions during collisions among cluster particles, we thus have to resort to some kind of simplified, binary collision picture. Since our kinetic model lacks any direct notion of cluster size or shape, we will stick to the assumption that, on average, collisions between cluster particles are devoid of any directional bias, leading to the same type of collision rule as for single particles: Again, ϑ constitutes a Gaussian-distributed random variable given in (5). Moreover, due to external (e.g. thermal background) and internal (e.g. noisy propelling mechanism) noise, cluster particles evaporate to become single particles. In analogy to the self-diffusion of single particles, we thus introduce a rate [43] characterizing the following evaporation process: Also in this case, the strength of the angular changes are Gaussian-distributed according to (2), and for simplicity we use the same standard deviation σ 0 as for the single particles' persistent random walk. As discussed above, cluster particles are strongly caged due to their close proximity to neighboring, collinearly moving particles. Reorientations of cluster particles due to noise are therefore strongly counteracted by realigning particle collisions, rendering cluster particles considerably less susceptible to random fluctuations than single particles. Hence, we assume λ, which is consistent with the observations in agent-based simulations slightly below the ordering transition [31], finding coherently moving clusters in an unpolarized background of randomly moving particles. In this regime individual particles exhibit superdiffusive behavior, performing quasi-ballistic "flights" as long as they are part of a cluster, and conventional particle diffusion if they are not.

B. Constitutive equations
Building on the modeling framework defined above, we now set up a kinetic description for the canonical model. We denote by s(θ, x, t) and c(θ, x, t) the one-particle distribution functions within the class of single particles and cluster particles, respectively, i.e. s(θ, x, t) dθ d 2 x gives the number of single particles located in an infinitesimal region [x, x + dx] with orientations in the interval [θ, θ + dθ] (and likewise for c(θ, x, t) dθ d 2 x). Both one-particle distribution functions are subject to convection due to the propelling velocity v of each particle. Moreover, local fluctuations in the one-particle distribution functions due to selfdiffusion and collision events are to be accounted for. We thus arrive at the following set of Boltzmann-like equations for the canonical model: where the source termsṡ(θ, x, t) andċ(θ, x, t) reaḋ They give the net number of single particles and cluster particles entering the phase space (i) Self-diffusion and evaporation. In these cases the source terms are products of the corresponding rates and probability densities, with (12) denoting the probability density for a particular species f to have a certain angle θ, and denoting the transition probability from θ = θ − ϑ 0 to θ averaged over all ϑ 0 with respect to the Gaussian weight (2). Note that, here and in the following, the argument of f is understood modulo 2π.
(ii) Collisions within the same class of particles. The collision integrals, representing the processes defined in equations (3) and (7), are given by standard expressions [16][17][18]22] C (+) Here ... denotes an average over ϑ ∈ (−∞, ∞) with respect to the Gaussian weight (5) and the average angleθ is given in Eq. (4). The integral measure contains the differential scattering cross section characterizing the frequency of collisions (i.e. hard-core interactions) between rod-like particles. The scattering function Γ itself carries all information concerning the shape of the particles and is a function of the relative orientation of the colliding particles. Reminiscent of the Boltzmann scattering cylinder, Γ can be derived on the basis of purely geometric considerations assuming that all spatial coordinates within the cylinder are equally probable; for details see A.
(iii) Assembly events of a single particle joining a cluster. These events, represented by (6), occur through binary collisions between single particles and cluster particles and are thus represented by an analogous integral expression:

III. DERIVATION OF HYDRODYNAMIC EQUATIONS
In order to reduce our kinetic description to a set of hydrodynamic equations valid on large length and time scales, we follow the well-established procedure of Aranson et al. [16] and Bertin et al. [17,22], and analyze the angular dependence of equations (10a) and (10b) in Fourier space. Due to the 2π-periodicity in θ, the one particle distribution functions can be expanded in Fourier series where , the zeroth and first Fourier modes are directly connected to the hydrodynamic densities ρ s (single particle density) and ρ c (cluster particle density), and the corresponding current density g s and g c , i.e.
In equations (20c) and (20d), the "=" signs indicate identification of vectors and complex numbers. The quantities u s/c denote the velocities of the macroscopic flow fields established by single particles and cluster particles, respectively. Also note that the second Fourier components are proportional to the nematic order parameter within the respective class of particles (as reflected by the symmetry of e i2θ under θ → θ + π). Using equations (18a) and (18b), the Boltzmann-like equations (10a) and (10b) transform to where the collision integrals I n,k are defined as follows: Note, in particular, that I 0,0 gives the total scattering cross section.

A. Truncation scheme
Equations (21a) and (21b) constitute an infinite set of coupled equations in Fourier space, which are fully equivalent to the Boltzmann-like equations (10a) and (10b). To derive a closed set of hydrodynamic equations, we need to consider some additional assumptions, allowing us to truncate this infinite Fourier space representation.
Here, our focus will be on virtually isotropic systems in the vicinity of an ordering transition breaking rotational symmetry. In this case, deviations of the one-particle distribution functions from the constant distribution ∼ 1/2π are small and contributions from large wave numbers in the Fourier series (21a) and (21b) are negligible. We further consider sufficiently dilute systems, in which the number of (binary) particle collisions per unit time and area [∼ (ρ c + ρ s ) 2 I 0,0 ] is much smaller than the corresponding number of single particle diffusion events [∼ λ ρ s ]. Together with λ [equation (9)], stating that disassembly from a cluster is strongly hindered by particle caging, allows us to treat single particle diffusion as a fast process. The single particle phase thus acts as an isotropic sea of particles where particle orientations (but not necessarily particle densities) are equilibrated, and hence the In summary, we resort to the following truncation scheme, leading to a set of hydrodynamic equations, valid near the onset of the ordering transition:

B. Derivation of the hydrodynamic equations
With the above truncation scheme, (21a) and (21b) reduce to and where Re(a) [Im(a)] denote the real [imaginary] part of a. Moreover, as can be seen from the definition in (22), the collision integrals I n,k only depend on the value |n − k/2|, whence only five of the collision integrals appearing in the above equations are independent. These integrals as a function of the particle's aspect ratio are evaluated and summarized in table I. Also note that the entire set of equations (24a) -(24d) is independent of the fast single particle diffusion time scale λ −1 (and, hence, also of the diffusion noise parameter σ 0 ). In our present approach, λ has only a conceptual meaning in maintaining a well-mixed particle bath within the class of single particles.
For given particle densities, the time scales governing the dynamics of the polar and nematic order parameter fields, represented by c 1 and c 2 , are given by the linear coefficients in the second line of (24c) and (24d), respectively. As will be detailed in IV B, the onset of collective motion is hallmarked by a change in sign of the linear coefficient in (24c), implying a diverging time scale for the dynamics of the polarity field. On the other hand the time scale for c 2 is finite for all densities which implies that the relaxation of the nematic order parameter field is fast compared to the polarity field. This allows us to set ∂ t c 2 ≈ 0 in (24d).
In the following it will be convenient to write down equations in dimensionless form. To this end, we construct the following characteristic scales: Time and space will be measured in units of the cluster evaporation time and length scalê From the cluster evaporation time scaleτ e and the total scattering cross section I 0,0 , we can construct the characteristic density scalê The single particle and cluster particle phase constantly exchange particles at rates that are determined by cluster evaporation ( ) on the one hand (cluster particles → single particles), and cluster nucleation due to particle collisions on the other hand (single particles → cluster particles), which occur with a rate ∼ ρ I 0,0 . Therefore the characteristic density scaleρ b marks the particle density, where both rates balance. In particular, gives the rate of inter-particle collisions relative to cluster evaporation events. Thus, the numerical quantity ρ/ρ b provides a direct measure expressing the competition between the randomizing effects of noise and the order creating effects of particle collisions, hallmarking the onset (and maintenance) of collective motion [29].
We thus arrive at the following rescaling scheme x where the characteristic scales for momentum (g) and scattering cross section (I n,k ) have been constructed from those of time, space, and density. In this rescaling the momentum current density is equal to one if the corresponding fluid element with a characteristic densitŷ ρ b , for which cluster evaporation and nucleation balance, is convected with the particle velocity v.
Then, upon eliminating c 2 from (24c) as discussed above, and using the relations between Fourier modes and hydrodynamic fields (for details see B), (20a) -(20d), equations (24a) -(24d) give rise to the hydrodynamic equations corresponding to the canonical model. In rescaled variables they read: where and where we have introduced the following abbreviations: Equations (28a) -(28c) capture the evolution of our canonical model system on a hydrodynamic level. More specifically, (28a) and (28b) describe the spatio-temporal evolution of the particle densities ρ s and ρ c . Since, by the assumptions underlying our model, no macroscopic flow of single particles can build up, only the density of cluster particles (ρ c ) is subject to convection. This implies that the genuine hydrodynamic momentum field g = g c + g s ≡ g c , is carried solely by the subset of cluster particles. Therefore we omit where L and d denote particle length and diameter, and where v is the particle velocity. The quantities I n,k /I 0,0 depend only weakly on the aspect ratio ξ. In particular, the signs of I n,k /I 0,0 do not change with ξ, leaving all our present conclusions made on the basis of the kinetic coefficients qualitatively unchanged.
the subscript c in (28b) and(28c) and denote g ≡ g c . The dynamics of both densities is, moreover, driven by source terms, as determined by the reactions discussed in section II A.
The gain and loss parts in these source terms of ρ c and ρ s are exactly balanced, such that the total density ρ = ρ c + ρ s is conserved. As an aside we note that any distinction between single particles and cluster particles is a purely conceptual matter. Experimentally, only the total density ρ and the momentum field g are accessible.
(28c), governing the evolution of the current density g, can be interpreted as a generalization of the Navier-Stokes equation to active systems. The terms on the right hand side of (28c) can be given the following interpretation: In the first line, the first two terms account for the local dynamics of g. They play a crucial role in establishing and maintaining a state of macroscopic flow, as will be detailed below. The Navier-Stokes equation itself, which conserves momentum, is devoid of these terms. In formal analogy to the Navier-Stokes equation, the density gradient in the first line together with the last term in the second line can be interpreted as a pressure gradient. This effective pressure is given by 1 2 ρ c + ζ − ν 2 g 2 , when neglecting the density-dependence of ν 2 . The last term in the first line is analogous to the shear stress term in the Navier-Stokes equation, with a kinematic viscosity ∼ ν −1 2 . The second line in (28c) is a generalization of the convection term to systems not obeying Galilean invariance, where all combinations of ∇ and factors second order in g transforming as vectors are allowed [41]. Finally, the last two lines describe couplings of the current density g and gradients thereof to density gradients. Note that the density gradients in these coupling terms are all of the same generic structure (29).
As already noted, the canonical model equations (28a) -(28c) conserve the total number of particles. To make this explicit, we define where ρ denotes the overall particle density, and η measures density difference between the two particle classes. The canonical model equations then attain the following form: The equation governing ρ expresses the overal conservation of particle number, whereas the source terms of equations (28a) and (28b) combine to determine the local dynamics of the relative density η in (32b).
Now we turn to the grand canonical model, where the single particle phase is coupled to a particle reservoir, resulting in a situation where single particles constitute an isotropic sea of particles which is maintained at a constant density ρ 0 s . Particle number conservation is now violated, and the only non-trivial density dynamics takes place within the phase of cluster particles. The hydrodynamic equations corresponding to the grand canonical model can be obtained immediately by setting in (28a) -(28c) the density of single particles to a constant value, yielding: One final remark is in order: The rescaling scheme introduced in equations (27a) -(27e) renders both, the canonical and grand canonical model equations virtually independent of particle shape. While these equations exhibit a weak dependence on the particles' aspect ratio L/d (via the rescaled collision integrals I n,k ), this dependence introduces only minor quantitative effects, which are negligible for all present purposes. To a good approximation we can thus set L/d = 1 while working with dimensionless variables, and assess the effects entailed by particle shape by restoring original units. Within our present approach, the effects of particle shape are purely quantitative, causing a numerical shift in the characteristic scales, but leaving the qualitative features of the problem unaffected. Deep within the ordered phase, i.e. for large densities, we indeed find a qualitative change of the ensuing hydrodynamic instability, as detailed in V. Nevertheless, this statement has to be taken with a grain of salt because corresponding threshold densities are far beyond the validity of the hydrodynamic equations.

IV. SPATIALLY HOMOGENEOUS SYSTEMS
To investigate the implications of the hydrodynamic equations, we start with the simplest case by analyzing spatially homogeneous solutions. These considerations will provide the basis for the study of spatially inhomogeneous systems, which will be the subject of section V. Dropping all gradients, the hydrodynamic equations for spatially homogeneous systems for the canonical model read For the grand canonical model we get In both cases, the density dynamics decouples from the momentum current dynamics and can be addressed separately.
In this section, our focus is on the stationary properties of the canonical and grand canonical model, respectively. While the dynamical approach to the stationary state is model dependent, the system's composition in terms of single particles and cluster particles, for given total density ρ, in the limit t → ∞ is identical in both cases [refer to (28a) -(28b) and (33a) -(33b)]. Since, moreover, the momentum current densities g obey identical dynamical equations, the ensuing analysis of the stationary state is equal for both models.
A. Crossover to clustering To assess the density difference between the cluster particle and the single particle phase η, we calculate the dynamical fixed point η * of (34b), attracting the dynamics of η(t) in the long time limit t → ∞: The defining equations (31a) -(31b) can be used to determine the corresponding (stationary) fixed point densities of single particles (ρ * s ) and cluster particles (ρ * c ) as a function of the total density ρ. figure 3 summarizes these findings: Upon increasing the total density ρ, the ratio η * /ρ continuously grows from η * /ρ = −1 at ρ = 0, asymptotically approaching η * /ρ = 1 as ρ → ∞. Based on the sign of η * , two density regimes can be distinguished: In the low density regime (ρ 1, η * < 0) particle collisions, underlying the formation of clusters, occur at much smaller rates than cluster evaporation events. Only a small fraction The stationary relative density η * /ρ, as well as the stationary cluster and single particle density, ρ * c /ρ and ρ * s /ρ, respectively. The larger ρ, the more cluster particles exist in the system. The of all particles organize themselves in clusters leading to a relatively dense population of single particles and correspondingly small density of cluster particles. In the high density regime (ρ 1, η * > 0), the situation is reversed: Large overall densities imply frequent particle collisions and, consequently, cluster formation and cluster growth dominate over cluster evaporation. In this regime, the number of cluster particles exceeds the number of single particles.
The crossover between the single particle dominated low density regime and the cluster particle dominated high density regime occurs at the crossover densityρ =ρ b = 1, where both, the single particle and the cluster particle populations are of equal size [i.e. η * (ρ) = 0]. The relation between the crossover density to clustering, and the geometrical shape of the constituent particles has been addressed previously in Ref. [37], based on agent-based simulations and a mean-field type analytical analysis. Using our definition of the crossover densityρ we can establish the corresponding relation simply by restoring original units [equation (26)]. Using packing fractionp ρ L d instead of particle density, and assuming for the sake of simplicity L/d 1, which allows us to estimate the particle surface A 0 Ld, we find:p which correctly reproduces the findings of Ref. [37] (taking into account that the cluster evaporation rate is assumed to be proportional to the inverse particle length, ∝ L −1 ). For the sake of completeness, we note that the definition of the clustering crossover density in reference [37] is based on the cluster size distribution, and thus does not necessarily coincide with our definition. We stress, however, that in our description the scaling structure in equation (37) is completely generic. It is an immediate consequence of the characteristic scales of our model and of the fact that the rescaled hydrodynamic model equations are (virtually) independent of particle shape. The structure of equation (37) is thus robust under an arbitrary redefinition of the (rescaled) crossover densityρ.

B. Homogeneous equations for momentum current density
Having examined the composition of the system in terms of single particle and cluster particle densities, we now turn to a discussion of the spatially homogeneous solutions for the momentum current density g. Due to rotational invariance of (34c), only the magnitude g = |g| of the momentum current density, but not its direction, evolves in time. We can thus concentrate on the scalar equation which leads to the following fixed points g * as the attractor of the dynamics of g in the limit of long times: It can be shown, that the coefficient in front of the cubic term in (38) is indeed strictly positive for all control parameters of density ρ and noise σ consistent with ν 1 < 0, ensuring the existence of the non-trivial fixed point in the second line of (39).
Depending on the sign of the linear coefficient ν 1 , two parameter regimes can thus be distinguished: Parameters leading to ν 1 > 0 render stable an overall homogeneous and isotropic state with vanishing macroscopic flow g = 0. Upon crossing the phase boundary in parameter space, the isotropic solution gets unstable and a macroscopic current density of non-zero amplitude builds up. In equation (40) we used that the density difference η, in the stationary limit, is a function of the total density ρ; cf. equation (36). Hence, in the limit of long times, ν 1 is a function of the total density ρ and the noise parameter σ, only.
Using the definition of the coefficient ν 1 , equation (30a), we can readily calculate the shape of the phase boundary in the σρ -plane: where we used I 1,0 = − 1 3 and I 1,1 = 1 2 . The corresponding phase diagram is shown in figure  4.
To conclude this section, we note that the analysis of spatially homogeneous systems corroborates the general physical picture of active systems, that was alluded to in the introduction (e.g. cf. reference [31]): Even in the absence of noise, σ = 0, for which the threshold density ρ (c) is lowest, the fully isotropic state g = 0 remains stable up to a critical density ρ (c) (σ = 0) = √ 3ρ, which lies well beyond the densityρ indicating the crossover to clustering. We thus extract the following physical picture; cf. figure 4: For low densities, ρ <ρ, cluster evaporation dominates over cluster assembly via particle collisions and clusters form only transiently. The system most closely resembles a structureless, isotropic "sea of particles". At intermediate densities,ρ < ρ < ρ (c) , particle collisions are more frequent. The emergence of clusters is now a virtually persistent phenomenon, with cluster evaporation occurring at a lower rate than cluster formation and growth. Yet, the collision rates between clusters (i.e. collisions among cluster particles) are still too low to orchestrate macroscopic order, leading to an overall isotropic "sea of clusters". Finally, for large densities, ρ > ρ (c) , the frequency of collisions among clusters is high enough to establish collective motion even on macroscopic scales.

V. STABILITY OF INHOMOGENEOUS HYDRODYNAMIC EQUATIONS
From our hitherto discussions, we have ascertained that the isotropic, homogeneous state (ρ = const. and g = 0) becomes unstable for sufficiently large densities. Yet, from a purely homogeneous analysis we cannot tell anything about the spatial structure of such a macroscopic broken-symmetry state. Nor can we be sure that the isotropic and homogeneous solution for ρ < ρ c (σ) is indeed stable with respect to spatially inhomogeneous perturbations.
In this section, we therefore test the linear stability of the homogeneous isotropic and nonisotropic base states with respect to wavelike perturbations of arbitrary wave number. Unlike the homogeneous model equations, the full hydrodynamic model equations are different for both, the canonical and grand canonical model, implying different dispersion relations describing the growth of such wave-like perturbations. We will thus analyze both models separately, and show that particle conservation does indeed influence pattern formation in essential respects.
A. Linearization about stationary, spatially homogeneous base states We start by linearizing the hydrodynamic equations for the canonical model. In the canonical model, the total number of particles is conserved, and the appropriate base state reads (cf. section IV) whereê g denotes the unit vector in the direction of the homogeneous polarization, and where all fields of the base states are assumed to be constant both in space and time. We are going to investigate the linear stability of the solutions (42a) -(42c) against wave-like perturbations, employing the following ansatz: where the perturbations are plane waves δη(x, t) = δη 0 e st+iq·x , (44b) In the equations above, q denotes the wave vector and s is the growth rate. Inserting this ansatz into the hydrodynamic equations (32a) -(32c), we obtain the following eigenvalue problem: Unlike the canonical model, the grand canonical model conserves the number of single particles, but not the total number of particles. The appropriate base state in this case reads: We investigate the linear stability of these solutions, using a perturbation ansatz analogous to equations. (43a) -(44c): δg(x, t) = δg 0 e st+iq·x .
Inserting this ansatz into equations (33a) -(33c), we obtain: B. Stability of the disordered state g 0 = 0 We start by considering the homogeneous and isotropic base state, which was shown to be stable against spatially homogeneous perturbations for ρ < ρ (c) (σ); cf. section IV B.
To assess the stability of this state with respect to perturbations of arbitrary (non-zero) The corresponding eigenvalue problem for the grand canonical model is found from equations (49a) -(49b), and attains the following form: For g h = 0, (45c) or (49b), respectively, implies q || δg 0 , allowing to replace the vectors q and δg 0 by their respective magnitudes q and δg 0 . We solved both eigenvalue problems numerically, for arbitrary wavenumbers q > 0 with the results shown in figure 5. Note that in both models the real parts of all eigenvalues are negative for all wavenumbers q > 0, provided the particle density ρ < ρ (c) . The spatially homogeneous, isotropic state is thus stable against small perturbations with arbitrary wavevectors.
For densities ρ > ρ (c) , in contrast, a narrow band of positive eigenvalues emerges in both models, located at wavenumbers q 1. Equations (50) and (51) evaluated at q = 0 return nothing but the linearized versions of the homogeneous hydrodynamic equations, (34a) -(34c) and (35a) -(35b). To gain new insights, we will therefore examine the limit q → 0 and consider the eigenvalues of the above coefficient matrices to leading order in the wavenumber q.
In this limit of small wavenumbers, the grand canonical coefficient matrix, given in (51) To sum up, the study of the linear stability of the homogeneous, isotropic state against spatially inhomogeneous perturbations of arbitrary wave vectors strongly suggests that particle conservation plays a vital role in the context of pattern formation. Both models exhibit spontaneous symmetry breaking by establishing a state of macroscopic collective motion. In the canonical model, in addition, conservation of total particle number entails a marginally stable density mode at q = 0 which is absent in the grand canonical model. This mode, in turn, gives rise to a density instability at small, non-zero wavenumbers, accompanying the spontaneous symmetry breaking event for ρ > ρ (c) . We note, however, that, at this point of the discussions, the existence of a narrow band of unstable modes at small wavenumbers does not allow for any conclusions concerning the structure of the macroscopic density and momentum current density for ρ > ρ (c) . We will address this issue in greater detail in the following section.
C. Stability of the broken symmetry state g 0 > 0 Both, the canonical and grand canonical model exhibit spontaneous symmetry breaking for overall densities ρ > ρ (c) (σ). To illuminate the spatial structure of this broken symmetry state, we start from the most simple case of a spatially homogeneous state of collective motion, and examine its stability with respect to wavelike perturbations in the hydrodynamic particle and momentum current densities. Without loss of generality, we assume the direction of the macroscopic momentum current density to coincide with the x-direction and choose g h = g 0êx . The wave vector q of the perturbation fields is assumed to make an angle ψ with the macroscopic momentum current density g h , yielding ψ = ∠(q, e x ) and q = q (cos (ψ), sin (ψ)) with q = |q|. The linearized canonical model equations (45a) -(45c) then attain the following form: The corresponding equations for the grand canonical model read: In equations (52a) -(53c) we used −ν 1 − µκ ν 2 g 2 0 = 0, which directly follows from the definition of g 0 , given in equation (39). We numerically solved both eigenvalue problems in the immediate vicinity of the ordering transition line ρ = ρ (c) (σ).
In the case of the canonical model, we find that the most unstable mode occurs for longitudinal perturbations, i.e. perturbations with wave vectors parallel to the direction of macroscopic motion, q || g 0 (ψ = 0). 6a shows the corresponding eigenvalues as functions of the wavenumber q for a set of density values slightly beyond ρ = ρ (c) . Further inspection of the coupling coefficients in equations (52a) -(52d) reveals that this longitudinal instability only affects the amplitude of g leaving the direction unchanged: For ψ = 0, the dynamics of δg y,0 decouples and momentum current density fluctuations perpendicular to the direction of macroscopic motion decay exponentially, with a rate which approaches zero for q → 0, as expected for a broken symmetry variable. To assess the nature of the instability in greater detail, we calculated the eigenvector corresponding to the most unstable longitudinal mode (evaluated at the most unstable wavenumber). It turns out that this eigenvector has approximately equally large components along the remaining fluctuation amplitudes δg x,0 , δρ 0 and δη 0 . This is consistent with our previous findings, indicating that the density mode, which was alluded to in section V B and which turns unstable at ρ = ρ c , renders the state of homogeneous collective motion unstable to fluctuations of the magnitude of the momentum current density. We further note that this picture is in agreement with previous numerical [31] and analytical [22] results (cf. 1).
The stability regions of the grand canonical model strongly deviate from the above picture. Setting ψ = 0 (longitudinal perturbations), we calculated the largest eigenvalue s (GC) max of the linear system of equations (53a) -(53c): which is always negative since ρ s < 1. In contrast to the canonical model, longitudinal perturbations thus always decay exponentially fast in the grand canonical model. For perturbations in transverse directions, in contrast, a positive eigenvalue can be found for sufficiently low noise levels σ < σ r , with the fastest growing modes posessing wavevectors q ⊥ g 0 . 6b shows the eigenvalue of the most unstable modes, which occur for σ = 0. To assess the implications of this instability for the dynamics of the various fluctuation amplitudes, we numerically examined the eigenvector corresponding to the positive eigenvalue, evaluated at the most unstable wavenumber. For densities in the vicinity of the ordering transition, we find that this eigenvector has approximately equal components in both momentum current density fluctuation amplitudes, δg x,0 and δg y,0 , but an essentially vanishing component along the direction of density fluctuations δρ 0 c . The corresponding instability can thus be classified as a hybrid shear/splay instability, leaving the spatially homogeneously distributed particle density virtually unaffected.
Three remarks are in order: First of all, for noise values σ c (ρ → ∞) > σ > σ r , the state of homogeneous collective motion becomes linearly stable with respect to arbitrary perturbations, including transverse perturbations. Secondly, even the most unstable eigenvalues "restabilize" for densities which are in the vicinity of the ordering transition threshold, which is depicted in 7. Finally, a restabilization can also be "observed" for the longitudinal instability in the canonical model. In this case, however, the restabilization occurs for relatively large densities and thus lies outside the range of validity of the linearized equations (52a) -(52d).

VI. DISCUSSION AND CONCLUSION
To conclude, we discuss and summarize our main findings. To study the onset of collective motion in active media, we started out with a simplified model for a system of self-propelled rod-like particles of variable aspect ratio. Collective motion was assumed to be established in a completely self-organized fashion solely by means of interactions among the constituent particles and in the absence of any external alignment fields. These interactions were assumed to occur via binary, inelastic particle collisions during which the rods align their direction of motion. Moreover, interactions were assumed to be subject to noise, which we controlled by a single model parameter σ. To assess some of the structural properties of such systems, we associated each of the particles with one of two classes: single particles and cluster particles, each with the corresponding density fields denoted as ρ s and ρ c . The class of cluster particles hosts all particles belonging to some coherently moving group of particles, which we referred to as cluster. The rest of the particles can be imagined to make up an isotropic sea of particles and are associated with the class of single particles. Using this classification scheme, we implemented simple interaction rules, representing cluster nucleation, cluster growth and cluster evaporation; the latter is assumed to occur at some fixed rate .
To illuminate the self-organization of collective motion, we set up an analytical, kinetic description of such systems, focusing on two archetypical modeling frameworks. Firstly, we considered isolated systems in which the total number of constituent particles is a conserved quantity. This case was referred to as the canonical model. Secondly, we examined open systems, which we referred to as the grand canonical model. Open systems are in contact with a particle reservoir which keeps the density of single particles at a constant level.
Inspecting the corresponding hydrodynamic equations, we were able to establish the following physical picture, portraying the formation of collective motion via dissipative particle interactions: For both, the canonical and the grand canonical model, we identified two characteristic density scalesρ and ρ (c) (σ), with ρ (c) (σ) >ρ, which allowed us to distinguish three density regimes.
For low densities, ρ <ρ, the rate at which particles collide is much smaller than the rate at which clusters disassemble. In terms of a particle based picture, this regime corresponds to a situation, where particle clusters are unstable, evaporating shortly after their nucleation.
In the stationary state, the vast majority of particles populates the single particle phase, rendering the system homogeneous and isotropic even on mesosopic scales. This low density regime terminates at the characteristic densityρ, where both classes exchange particles at equal rates.
In the contiguous regime of intermediate densities,ρ < ρ < ρ (c) , the overall rate of cluster formation and growth outstrips the rate at which clusters evaporate, and the majority of particles becomes organized in clusters. Translated to a particle based notion, clusters grow to finite sizes and persist over macroscopic time scales. Clusters of coherently moving particles now dominate the physical picture on mesoscopic scales. Yet, interactions among clusters are too rare to establish a macroscopic state of collective motion. On hydrodynamic length scales, the system can be viewed as a homogeneous and isotropic sea of clusters.
For densities exceeding the critical density, ρ > ρ (c) (σ), collisions within the cluster phase occur at sufficiently high rates, and macroscopic collective motion emerges. The homogeneous and isotropic state, which has been shown to be stable within the two preceding regimes, thus gets unstable and rotational symmetry is spontaneously broken. Resorting to a particle based image, we can imagine the mean cluster size to reach a "percolation threshold", leading to coagulation and net alignment between clusters. This is in agreement with previous analytical [22] and numerical [31] results for particle conserving systems, where the emergence of solitary wave structures has been reported in the vicinity of the ordering transition ρ ρ (c) (σ). The longitudinal instability thus seems to be a quite generic feature of particle conserving systems with hard core interactions. For an interesting counter example we refer the reader to Ref. [26], where a particle conserving system with topological interactions has been studied.
We can now combine our findings for both, the canonical and the grand canonical model, to offer the following mechanistic explanation concerning the emergence of the longitudinal instability. The prerequisite, underlying the establishment of coherent motion, is embodied by two basic processes: Cluster nucleation by collisions among single particles, and cluster growth by alignment of single particles to clusters. Only if, by virtue of these processes, the concentration of cluster particles grows sufficiently large, clusters are able to synchronize their movements by coagulation and macroscopic collective motion emerges.
Now consider the effect of a density fluctuation in an otherwise homogeneous state of macroscopic collective motion. In the grand canonical model, where the density of single particles is kept fixed by virtue of a particle reservoir, this fluctuation occurs within the class of cluster particles. We can use the right hand side of equation (35a), to assess the implications of such a fluctuation on the local composition of the system in terms of cluster particles and single particles: Note that this equation captures the balance of the two particle currents between the single particle and the cluster particle phase in the stationary limit. As can be seen from this equa-tion, locally enhancing the density of cluster particles implies a net current from the cluster particle phase into the single particle phase, thus counteracting the effect of the original density fluctuation. Conversely, locally diminishing the density of cluster particles leads to the opposite effect. Density fluctuations are thus damped in the grand canonical model and do not impact the macroscopic velocity field, which is set up by the cluster particles.
Exactly the opposite happens in the particle conserving canonical model.
where the left hand side describes cluster nucleation and condensation, and the right hand side corresponds to cluster evaporation. This can be seen by using the definitions of the relative density η = ρ c − ρ s , and the total particle density ρ = ρ s + ρ c . Now, consider a fluctuation in the total density ρ, where, for the sake of simplicity, we assume the relative density η to remain constant. In regions, where the fluctuation leads to an increase in the total density by a factor k > 1 we have Hence, the particle current into the cluster particle phase grows. As a consequence, the local value of the momentum current density increases, since the cluster particles are the "carriers" of the macroscopic momentum. In contrast, in regions, where the fluctuation decreases the total density by a factor k < 1 we have There the cluster particle phase gets depleted and the local magnitude of the momentum current density declines. As a result, high density regions move at faster speeds than low density regions, gathering more and more particles on their way through the system. Conversely, lower density regions continually lose particles to the faster high density structures.
In particle conserving systems, every density fluctuation thus automatically triggers a corresponding fluctuation in the momentum current density, which in turn amplifies the density fluctuation. As a result of this process, high density bands of collectively moving cluster particles might emerge [31]. These bands being interspersed by regions where the particle density has fallen below the critical density ρ (c) (and possibly belowρ), leading to local destruction of clusters and collective motion.
We close by adding some remarks on the importance of the particles' shape on the establishment of collective motion on hydrodynamic scales. We found that the impact of particle shape on the macroscopic properties of such systems is purely quantitative in the framework of our present study: Varying the particles' aspect ratio results in a shift of the characteristic density scalesρ and ρ (c) (σ), which we quantified in equation (37). Qualitatively, our conclusions concerning the macroscopic properties of these systems remain unaffected by a change in the particles' aspect ratio. Note that, in our approach, the aspect ratio basically determines the total scattering cross section and thus "merely" impacts the rate at which particles collide. We stress, however, that in real systems particle shape is likely to have a profound impact on the entire physical picture of particle interactions, and not just on their rate. The study of those effects lies outside the scope of our present work and would be an interesting topic for future research.

Appendix A: Derivation of Boltzmann collision cylinder for driven rods
In the framework of our Boltzmann-like description, binary collisions, such as equations (3), (6) and (7), occur with a certain rate Γ depending on particle shape (L and d), relative angle of both collision partners θ 12 = |θ 1 − θ 2 | and the constant velocity v. The quantity Γ(L, d, θ 12 ) characterizes the collision area per unit time -more commonly referred to as Boltzmann collision cylinder. On the scale of the Boltzmann equation, binary collisions occur locally, say in an infinitesimal volume element centered at r. Assume that particle 1 has an orientation θ 1 . Then, Γ(L, d, θ 12 ) dt gives the area around particle 1 in which every particle with orientation θ 2 will collide during a time interval [t, t + dt] with particle 1 . As a consequence Γ(L, d, θ 12 ) f (r, θ 1 , t) f (r, θ 2 , t)dθ 1 dθ 2 equals the number of collisions per unit time and unit area at time t, with f (r, θ, t) denoting the one-particle distribution function.
To determine Γ(L, d, θ 12 ), we take a microscopic point of view. Since the model employed in this work assigns to each particle a velocity vector pointing along its rod axis, we can distinguish "head" and "tail". Referring to figure 8, without loss of generality we assume π − θ 12 ≡ θ ∈ [0, π] (negative relative angles lead to the same result), and consider the blue rod, with the position of its head indicated by the blue dot. All rods of relative orientation θ 12 = θ 1 − θ 2 , and with their heads lying in the area S = A ∪ S 1 ∪ S 2 at time t, will collide with the blue rod during the time interval [t, t + dt]. Since A, S 1 and S 2 are disjoint, where |X| denotes the area of the region X. The respective areas are given by: and Returning to the laboratory frame we have v rel = v|ê(θ 1 ) −ê(θ 2 )| = 2 v| sin(θ 12 /2)|. Noting that Γ = |S|/dt (cf. eq. (A1)), we find: Γ(L, d, θ 12 ) = 4v d sin In figure 9, Γ(L, d, θ 12 ) is shown as a function of relative angle θ 12 for different particle lengths, whereby the particle width d is kept fixed. Increasing L/d shifts the most probable collision from θ 12 = π for L/d = 1 (the case of a sphere; θ 12 = π leads to the largest value of the relative velocity), towards θ 12 = π/2 for L/d → ∞ (limiting case of a needle; largest target area for θ 12 = π/2). Here we collect all such (complex) gradient terms and give a brief derivation of their vector-analytic counterparts. As in the main text, we identify C and R 2 , i.e.
To distinguish (genuinely) complex from purely real quantities, we assume f ∈ C and ρ ∈ R in the following. For the figure, we chose for particle width d = 1 and for particle velocity v = 1. Increasing the aspect ratio L/d, the most probable collision approaches θ 12 = π/2, whereas for L/d = 1 the most probable collision is the head-head collision with θ 12 = π.
(∂ x + i∂ y )ρ Using (B1) we immediately obtain By straightforward expansion we find Decomposing f into real and imaginary part and expanding, we find +i ∂ y f 2 y − ∂ y f 2 x + 2∂ x (f x f y ) where e j denotes the j-th Cartesian unit vector.
Expanding and collecting real and imaginary parts, we find Thus, We find Hence,