Large density expansion of a hydrodynamic theory for self-propelled particles

Large density expansion of a hydrodynamic theory for self-propelled particles

Thomas Ihle Department of Physics, North Dakota State University, Fargo, ND 58108-6050, USA

Recently, an Enskog-type kinetic theory for Vicsek-type models for self-propelled particles has been proposed [T. Ihle, Phys. Rev. E 83, 030901 (2011)]. This theory is based on an exact equation for a Markov chain in phase space and is not limited to small density. Previously, the hydrodynamic equations were derived from this theory and its transport coefficients were given in terms of infinite series. Here, I show that the transport coefficients take a simple form in the large density limit. This allows me to analytically evaluate the well-known density instability of the polarly ordered phase near the flocking threshold at moderate and large densities. The growth rate of a longitudinal perturbation is calculated and several scaling regimes, including three different power laws, are identified. It is shown that at large densities, the restabilization of the ordered phase at smaller noise is analytically accessible within the range of validity of the hydrodynamic theory. Analytical predictions for the width of the unstable band, the maximum growth rate and for the wave number below which the instability occurs are given. In particular, the system size below which spatial perturbations of the homogeneous ordered state are stable is predicted to scale with where is the average number of collision partners. The typical time scale until the instability becomes visible is calculated and is proportional to .

1 Introduction

The emergence of collective motion of living things, such as insects, birds, slime molds, bacteria and fish is a fascinating far-from-equilibrium phenomenon which has attracted a great deal of cross-disciplinary attention vicsek_zafeiris_12 (). The study of these systems falls under the broader umbrella of “active matter”. This term stands for collections of agents that are able to extract and dissipate energy from their surrounding to produce systematic motion ramaswamy_10 (); marchetti_13 (). Non-living examples of active matter are chemically powered nanorods rueckner_07 (), networks of actin fibers driven by molecular motors humphrey_02 () and swarms of interacting robots rubenstein_14 (); wilson_14 (). In 1995, a minimal computational model that captures the essentials of collective motion without entering too much detail was introduced by Vicsek vicsek_95_97 (). In this so-called Vicsek model (VM), agents are represented as point particles traveling at constant speed. The agents follow noisy interaction rules which aim to align the velocity vector of any given particle with its neighbors. Later, more sophisticated models that include additional interactions such as attraction and short-range repulsion were introduced levine_00 (); couzin_02 (); romenskyy_13 (); grossmann_13 (). On the theoretical side, coarse-grained descriptions of the dynamics in terms of phenomenological hydrodynamic equations have been proposed on the basis of symmetry and conservation law arguments.Well-knowns examples are the Toner-Tu theory for polar active matter toner_98 (); toner_12 () and the theory by Kruse et al. for active polar gels kruse_05 (). While being very successful, a disadvantage of these approaches is that no link between microscopic collision rules and the coefficients of the macroscopic equations is provided. In particular, since all coefficients must be related to only a few microscopic parameters, they cannot be varied independently and the actual parameter space should be more restricted than the hydrodynamic equations might suggest. Furthermore, by just postulating equations there is the immanent danger that some terms or even entire equations might have been ommitted, which could become relevant in new, previously untested, situations. To address the missing-link issue, several groups have derived hydrodynamic equations from the underlying microscopic rules and provided expressions for the transport coefficients in terms of microscopic parameters aranson_05 (); bertin_06 (); mishra_10 (); ihle_11 (); menzel_12 (); romanczuk_12 (); peshkov_12 (). One of the first attempts to put the Toner-Tu theory on a microscopic basis was published by Bertin et al., who studied a Vicsek-type model with continuous time-dynamics and binary interactions by means of a Boltzmann approach bertin_06 (); bertin_09 (). While the authors recently clarified bertin_14A () that even at low density their underlying microscopic model is not identical to the Vicsek-model, many predictions agree qualitatively with the ones for the VM. This Boltzmann approach was later extended to systems of self-propelled particles with nematic and metric-free interactions, and hydrodynamic equations were drived peshkov_12 (); ngo_12 (). Very recently, the very basis of these derivations – the Boltzmann equation and its validity in active matter – has been critically assessed thueroff_13 (); ihle_14 (); ihle_14A (). In particular, it has been shown that, at least near the threshold to collective motion, the binary collision assumption is not valid in the VM at realistic particle speeds and even at very low densities ihle_14 (). Furthermore, it has been demonstrated that the mean-field assumption of molecular chaos is not justified near the threshold to collective motion in soft active colloids and that the Boltzmann theory must be amended by pre-collisional correlations hanke_13 (). A first attempt to rigorously include correlations by means of a ring-kinetic theory for Vicsek-type models was put forward by Chou et al. chou_14 (). This theory is very complicated; work is still in progress to simplify it ihle_15 ().

Nevertheless, arguments from ring-kinetic theory can be used to confirm the plausible presumption that mean-field theories become reliable in the VM at large particle speeds (or time steps) and/or at large particle densities FOOT2 (). Here, large density means that the average number of collision partners of a given agent is much larger than one. It is known from equilibrium statistical mechanics that the behavior of spin-models become mean-field or more mean-field-like if the number of interaction partners is increased which can either be achieved by extending the range of interaction or by increasing the spatial dimension. It appears reasonable to assume a similar tendency in the VM which consists of moving “spins”.

The large density limit is beyond the capability of Boltzmann approaches because Boltzmann equations are restricted to binary interactions. However, a recently proposed Enskog-like theory ihle_11 (); ihle_13 (); ihle_14 () has no limitation on density and has already be applied to a model with and metric-free interactions chou_12 (). This Enskog-like theory is based on an exact equation for a Markov chain in phase space. In a previous paper ihle_11 (), hydrodynamic equations were derived from this theory and its transport coefficients were given in terms of infinite series. In this paper, I analyze the transport coefficients in the large density limit, , and show that they take a simple form. This allows me to analytically evaluate the well-known density instability of the polarly ordered phase near the flocking threshold chate_04_08 (); bertin_06 (), and to obtain simple formulas and scaling laws for the dispersion relation. Note that a similar analysis for the opposite limit has already been performed by Bertin et al. bertin_09 ().

The large density expansion performed in this paper is supposed to be benefitial for several reasons. First, one does not have to worry about the validity of the underlying mean-field assumption of Molecular Chaos: As discussed above, this assumption is expected to be asymptotically valid for . Thus, the obtained relations should be quantitatively correct in this limit. As will be shown below, at large density, the restabilization of the homogeneous ordered phase happens closer to the threshold. Hence, this effect occurs within the validity domain of the hydrodynamic theory, and, aiming at quantitative agreement, one is not forced to evaluate the kinetic eqution directly as proposed in Ref. bertin_09 (). Second, biologically realistic models of fish schools tegeder_95 () and bird flocks ballerini_08 () assume that is between 2 and 7, thus an approach for large seems promising. Third, similar to the opposite limit, , the large density approach leads to a simplification of the expressions for the transport coefficients and allows a physical interpretation of the instability close to the threshold of collective motion.

An obvious concern about the large density limit is how realistic it is for active matter in general. Of course, for example, for systems of active colloids with short-ranged excluded volume interactions it will be problematic. However, it does make sense for systems with long-ranged alignment interactions. Speaking in terms of the terminology put forward by Couzin et al. couzin_02 (), the condition can occur in a model with a large “zone of orientations” and a very small “zone of repulsion” Furthermore, in the process of doing the actual calculations and performing expansions in powers of , one realizes that error terms are often proportional to . For example, for , this is already a reasonably small number. Therefore, it is possible that the large density formulas derived in this paper could still be valid at slightly larger than one without loosing too much accuracy.

To the best of my knowledge, almost no analytical results exist for the Vicsek-model at high particle density. This paper is intended to fill this void. One of the main results is the prediction that the homogeneous ordered phase remains stable to small perturbations if the linear system size is chosen to be smaller than . Another main result is that in systems larger than the number of time steps to develop an instability is equal to or larger than .

2 Vicsek model

In the two-dimensional Vicsek-model (VM) vicsek_95_97 (); nagy_07 (), particles with average number density move at constant speed . The system of area evolves in discrete time intervals of size . The particles are assumed to have zero volume. The evolution takes place in of two steps: streaming and collision. During streaming, the particle positions are updated in parallel according to


Because the particle speeds remain the same at all times, the velocities, , are parametrized by the “flying” angles, , . In the collision step, the directions are modified in the following way: a circle of radius is drawn around the focal particle , and the average direction of motion of the particles within the circle is determined according to


where the sum goes over all particles inside the corresponding circle (including particle ). The new flying angles are then given as,


where is a random number uniformly distributed in the interval . In this paper, a forward-upating rule baglietto_08_09 (), corresponding to the so-called standard Vicsek-model, is used: the updated positions instead of the “old” positions enter the calculations of the average directions . Important dimensionless parameters of the model are (i) the noise strength , (ii) the average number of particles in a collision circle, , which is basically a normalized density, and (iii) the ratio of the range of the interaction to the mean displacement during streaming which can be interpreted as a mean free path (mfp).

3 Hydrodynamic theory

3.1 Macroscopic equations for the Vicsek model

Using Boltzmann’s assumption of molecular chaos, an Enskog-type kinetic equation for the one-particle density was derived from the exact evolution equation of the -particle probability. This was done for the standard Vicsek-model (VM) ihle_11 () as well as for other microscopic models ihle_09 (); chou_12 (); romensky_14 (). By means of a Chapman-Enskog expansion and an angular Fourier transformation of the one-particle density,


hydrodynamic equations were obtained for the dimensionless particle density and momentum density of the VM ihle_11 (). In this formulation, all times are scaled by the time step , all lengths are scaled by the mean free path (mfp) , , and all other quantities are scaled accordingly. For example, the actual two-dimensional number density and momentum density can be obtained from their scaled versions as and . To emphasize rotational invariance, the hydrodynamic equations are given in terms of the tensors , , and which depend on spatial derivatives of and :


with . Here, is the amplification factor for the first Fourier mode, which also determines the linear growth rate of the momentum density. This is because density and momentum density are the zeroth and first order moments of the one-particle distribution function, respectively, and therefore, the components of the momentum density are proportional to the first Fourier coefficients, , . In Ref. ihle_11 (), in the thermodynamic limit , an infinite series was found for ,


where is equal to the average angle defined in Eq. (2). Here, is the local average number of particles within a collision circle centered around the position ,


For a homogeneous system, one has . If is larger than one, the disordered state of zero average motion becomes unstable and a polarly ordered state with non-zero global momentum develops. Thus, setting defines the threshold noise of the flocking transition, at the level of homogeneous mean field theory. The hydrodynamic equations, Eqs. (5) and (6), are only valid in the vicinity of the transition to polar order, that is, for . This condition was essential to obtain a closed set of just two equations, see Refs. ihle_11 (); ihle_14 (). Physically, this assumption means that the first order Fourier coefficients and are evolving so slowly that the higher order coefficients, , become enslaved to them FOOT1 (). Therefore, no additional equations for the temporal evolution of the higher order coefficients are needed.

The momentum flux tensor and the tensors , ,


are given in terms of five symmetric traceless tensors ,


The tensor is the viscous stress tensor of a two-dimensional fluid. The transport coefficients in Eq. (9) were obtained in the additional limit of large mean free path, , and are given in Table 1. This means, contributions from the so-called collisional momentum transfer, which are supposed to be relevant at small mean free paths, see Ref. ihle_09 (), are neglected here. The rationale for this choice was that at small mean free paths and not too large density, the mean-field assumption (on which our theory is based) is supposed to break down anyway.

The transport coefficients depend on the following variables,

Table 1: The transport coefficients , and , defined in Eq. (9), are expressed as functions of , , , , see Eq. (11).

where is the ratio of the interaction radius to the mfp, . Note that these coefficients have a non-local dependence on position through the normalized density . The transport coefficients contain the following four types of dimensional integrals,


where is given by , ,
, and . The average angle is a function of the angles , see Eq. (2).

The Navier-Stokes-type equation, Eq. (6) is consistent with the one from Toner-Tu theory toner_98 (); toner_12 () but contains additional gradient terms. The additional terms are a result of a Chapman-Enskog expansion that includes all terms up to third order in a formal expansion parameter. Discussions of higher order expansions can be found in Refs. ihle_14 (); peshkov_14 (); ihle_14A (); bertin_14A (). Eq. (6) has a simple homogeneous flocking solution: and . The amplitude of the flow is given by


3.2 Large expansion

Using an analogy to the freely-jointed chain model for polymers, the angular integrals, Eqs. (12), can be asymptotially evaluated in the limit ihle_15 (). To leading order, the results are ihle_10_V1 (),


The angular integral in the definition of , Eq. (7), is evaluated similarly,


Even with these approximations, the infinite series in Eqs. (7) and (11), still cannot be calculated exactly. However, we note that they can be written as an average over a Poisson distribution. This is a consequence of the mean-field assumption of molecular chaos which predicts that the density fluctuations are equal to the ones of an ideal gas, e.g. are Poisson distributed. For example, using Eq. (15) we can rewrite the amplification coefficient from Eq. (7) as


where the average of an arbitrary function is defined as


Expanding around the mean value of the distribution one obtains,


where is the derivative of evaluated at , and are the so-called central moments,


which are known for the Poisson distribution. For example, one has


Thus, the following approximation is obtained,


To evaluate transport coefficients, averages over polynomials in are needed. Therefore, we investigate a generic example, with some power and a positive constant . The approximation, Eq. (22), then gives for small


This approximation becomes accurate for . For very large it suffices to replace by . Thus, replacing by in Eq. (16), leads to the leading order approximation for the amplification factor ,


Thus, fulfilling the critical noise condition in the limit of infinite density requires that the critical noise becomes equal to . Expanding around its infinite density limit gives,


Approximation (25) leads to


Introducing the relative noise distance to the threshold,


we consider the difference of the amplification factor to its value at the threshold, and find,


As a preliminary to calculate transport coefficients, we express the variables , , and from Eqs. (11), in terms of averages over the Poisson distribution,


and evaluate the averages to leading order by means of Eqs. (14). Finally, using Eq. (26), at the threshold one obtains


These expressions, Eqs. (30) are used to evaluate the transport coefficients from Table 1 in the limit of large and at the flocking threshold. Only the leading order contributions in the limit of large are given. For example, the product in coefficients such as , and so on is neglected because decays less strongly in the infinite density limit. The results are given in Table 2.

Table 2: The transport coefficients , and , defined in Eq. (9), evaluated at the threshold , in the limit .

It is instructive to evaluate the momentum density of a homogeneous, ordered system, close to the threshold, at , According to Eqs. (13), (28), and Table 2, we find


Even if the distance to the flocking threshold , is tiny, the amplitude eventually diverges in the infinite density limit. It also has a strong dependence on the mean free path ratio and diverges for infinite mean free path . This seemingly unphysical behavior is just a consequence of the particular choice of dimensionless variables where all lengths are measured in units of the mean free path , and all times in units of the time step . The amplitude is a measure of global polar order but it is not a convenient variable in the large density limit. Instead, I “renormalize” be relating it to the polar order parameter , which is defined as the average speed of a particle divided by the individual particle speed . This way, will always be between zero and one, and a small is expected near the order-disorder transition. One has


where is the dimensionless particle density. It is related to the regular number density of a two-dimensional system, by Using the definition, we find and an expression for is obtained,


Inserting from Eq. (31), the order parameter near threshold follows as


In contrast to , the order parameter does not depend on particle speed. This is the expected result because, (i) we work in the mean-field approximation, and (ii) the collisional contribution to the transport coefficients was neglected earlier.

3.3 Stability analysis

A long wavelength instability has been reported in Vicsek-type models bertin_06 (); bertin_09 (); ihle_11 (); ihle_13 (). This instability is strongest for longitudinal sound-wave-like perturbations and exists in the ordered phase at a range of noise values , next to the flocking threshold. In this chapter, I reanalyse this instability in the large -limit for a longitudinal perturbation with wave vector . Without loss of generality, the collective motion is assumed to go into the x-direction, . As small perturbation of amplitude is imposed on the homogeneous particle density . A similar perturbation is added to the x-component of the homogeneous momentum density , where is given in Eqs. (13) and (31),


and is the complex growth rate of the perturbation. From the continuity equation, Eq. (5), a simple relation between and follows,


The hydrodynamic equation for the momentum density, Eq. (6), in linear order of the perturbations leads to


Here, all transport coefficients are evaluated for a homogeneous system with constant density . The density derivatives of the coefficients , and are also needed and, for example, are given as


Using Eq. (24), one finds


This expression is then evaluated at the threshold with the help of Eq. (26), leading to


Similiarly, one obtains


Substituting Eq. (36) into (37), a quadratic equation for is obtained,


with complex coefficients , . The real and imaginary parts of and are found as,


where Eq. (13) has been used to eliminate . Eq. (42) has two solutions, , where is the one with the larger real part. Solving for the real part, one finds,


with the abbreviations


Instead of directly analysing this complicated expression, a different approach is pursued: Splitting the growth rate into its real and imaginary part, , Eq. (42) can be rewritten as two equations for and ,


The sign of the real part determines whether the system is linearly stable. We therefore eliminate from Eq.(46) and obtain a quartic equation with real coefficients for the real part of ,


with coefficients


This equation can be solved exactly by radicals as given in Eq. (44). However, in the high density limit and near threshold , the terms in the quartic equation can differ by several orders of magnitude and the analytical expressions become difficult to evaluate. Furthermore, our aim is to gain a better physical understanding and to obtain simple expressions for relevant quantities like the maximum growth rate of a perturbation of the ordered state. Therefore, in the next section, instead of relying on the exact solution of the quartic equation, various scaling regimes of are identified by different Ansatzes to solve Eq. (47) approximately.

3.4 Analysis of the dispersion relation

3.4.1 The limit of vanishing wave number

To see whether there is an instability at all, the limit of zero wave number is investigated by inserting the “hydrodynamic” Ansatz


into Eq. (47) and expanding the coefficients from Eq. (48) in powers of ,


with new coefficients that do not depend on , for example,


Collecting terms of order and , we find


To simplify the analysis we evaluate the dispersion in the large density limit and replace , using Eq. (33) by . In order to trace the origin of stabilizing and destabilizing effects, mathematical “markers” are put on the transport coefficients. This is done by taking the asymptotic expressions from Table 2 and multiplying them with dummy variables that become equal to one in the limit . For example, we have


where and the “marker”-variables are given by , , and so on. Inserting these “marked” coefficients into Eqs. (52) and (53), one finds that for vanishing order, , the coefficient is always positive but becomes negative for a strongly ordered state, . To obtain quantitative results and to investigate the sign change of , the dominant positive and negative terms must be balanced. Balance is achieved by assuming that the order parameter is of order or smaller. Mathematically, this is imposed by introducing the scaled order parameter and assuming that is at most of order one. Replacing in the equations for the coefficients and , Eqs. (48) and (52), and neglecting terms of higher order in , the coefficients dramatically simplify to,