Self-propelled particles with selective attraction–repulsion interaction: from microscopic dynamics to coarse-grained theories

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-) micrometre 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 self-propelled 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–12], where the spatial dimensions can be in the order of kilometres.

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 their 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 [13]. 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 [14, 15]. 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 [16, 17].

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 self-propelled particles with velocity-alignment [18–21] and self-propelled aligning rods [22]. Furthermore, coarse-grained descriptions for active particles with variable speeds and velocity-alignment were derived in [23–27].

In this paper, we will derive a kinetic and a hydrodynamic description for self-propelled 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 [28, 29].

A fundamental difference between our model and the Vicsek model is its formulation in continuous time using stochastic differential equations. The interaction of individuals is modelled by superposition of (binary) interaction forces. This allows, on the one hand, a straightforward 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 (FPE) for the corresponding probability density functions (PDFs), 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 FPE for the one-particle density function, which allows efficient numerical analysis of the stability of solutions of the FPE in Fourier space. Further on, we will derive a hydrodynamic theory for self-propelled particles with selective attraction–repulsion 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 PDF, 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 expansions. Finally, we compare the results of the kinetic and hydrodynamic theory with direct numerical simulations of the individual-based model.

We consider N self-propelled particles of mass m = 1 in two spatial dimensions, so-called agents. Each individual moves at a constant speed s₀, 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

$$\begin{equation} \frac{{\rm d} \boldsymbol{r}_i}{{\rm d} t} = s_0 \,\boldsymbol{e}_{i,h}(t) = s_0 \! \left( \begin{array}{@{}c@{}} \cos\varphi_i(t) \\ \sin\varphi_i(t) \\ \end{array}\right), \end{equation} \tag{ 2.1a }$$

$$\begin{equation} \frac{{\rm d} \varphi_i}{{\rm d} t} = \frac{1}{s_0}\left( F_{i,\varphi} + \sqrt{2 D_\varphi} \, \xi_i (t)\right) . \end{equation} \tag{ 2.1b }$$

Here, F_i,φ = F_i·e_i,φ is a projection of an effective social force vector F_i on the angular degree of freedom φ_i. It corresponds to a torque, which induces a turning behaviour of the focal individual due to social interactions. e_i,φ = (− sin φ_i,cos φ_i)^T is the angular unit vector perpendicular to e_i,h. In the following, despite considering here only agents moving with constant speed, we will use for simplicity the term 'forces' instead of 'torques' when referring to the social interactions. This takes into account that the general social forces introduced here naturally extend to the variable velocity case, where the social forces induce not only turning behaviour, but also accelerations and deceleration of individuals.

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. $\left < \xi _i(t) \xi _j(t+\tau ) \right \rangle = \delta _{ij} \, \delta (\tau )$ (Gaussian white noise).

The total social force is given by a sum of three components:

$$\begin{eqnarray} &&\boldsymbol{ F}_i = \boldsymbol{ f}_{i,r}+\boldsymbol{ f}_{i,m} +\boldsymbol{ f}_{i,a}. \end{eqnarray} \tag{ 2.2 }$$

The first term represents a short-ranged repulsion allowing for finite sized agents. It reads

$$\begin{eqnarray} \boldsymbol{ f}_{i,r} & = -\sum_{j=1j\neq i}^N \mu_r({r}_{ji}) \hat{\boldsymbol{r}}_{ji} \theta(l_{\rm c} - r_{ji}) \end{eqnarray} \tag{ 2.3 }$$

with μ_r(r_ji) ⩾ 0 being a repulsive turning rate, which, in general, depends on distance $r_{ji} = \left | \boldsymbol {r}_{ji} \right |=\left | \boldsymbol {r}_j - \boldsymbol {r}_i \right |$ 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

$$\begin{equation} \boldsymbol{ f}_{i,m} = \sum_{j=1j\neq i}^N \mu_{\rm m}({r}_{ji}) \hat{\boldsymbol{ r}}_{ji} \theta(l_{\rm s} - r_{ji}) \theta(r_{ji} - l_{\rm c}) \left | \tilde v_{ji} \right | \theta(+\tilde v_{ji}), \end{equation} \tag{ 2.4a }$$

$$\begin{equation} \boldsymbol{ f}_{i,a} = \sum_{j=1j\neq i}^N \mu_{\rm a}({r}_{ji}) \hat{\boldsymbol{ r}}_{ji} \theta(l_{\rm s} - r_{ji}) \theta(r_{ji} - l_{\rm c}) \left | \tilde v_{ji} \right | \theta(-\tilde v_{ji}). \end{equation} \tag{ 2.4b }$$

These forces can be considered as a sum over binary interactions, which act always along the unit vector $\hat {\boldsymbol { r}}_{ji}=(\boldsymbol { r}_j-\boldsymbol { r}_i)/r_{ji}$ pointing towards the centre 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 $\tilde v_{ji}$ defined by the temporal derivative of the distance ${\dot r}_{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 vector $\hat {\boldsymbol {r}}_{ji}$: $\tilde v_{ji}={\bf{ v}}_{ji}\cdot \hat {\boldsymbol {r}}_{ji}=\tilde v_{ij}$. Hence f_i,a is the response to approaching individuals characterized by a negative relative velocity $\tilde v_{ji} {<} 0$, whereas f_i,m is the corresponding response to moving away (or receding) individuals characterized by positive relative velocity $\tilde v_{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 $\tilde {v}_{ji}=\dot {r}_{ji}$, one obtains the relative velocity

$$\begin{eqnarray} &&\tilde v_{ji} = \tilde v_{ji}(\varphi_i,\varphi_j,\alpha_{ji})=2 {s}_0 \sin\left(\frac{\varphi_j+\varphi_i}{2}-\alpha_{ji}\right)\,\sin\left(\frac{\varphi_i-\varphi_j}{2}\right) \end{eqnarray} \tag{ 2.5 }$$

as a function of the velocity angles φ_i, φ_j and the angle α_ji of the distance vector r_ji, with $\boldsymbol {r}_{ji}= r_{ji} \left ( \cos \alpha _{ji}, \sin \alpha _{ji} \right )$. 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π:

$$\begin{eqnarray} \begin{array}{l l l } \tilde{v}_{ji} > 0 & \mbox{for} & \displaystyle\frac{\varphi_j+\varphi_i}{2} < \alpha_{ji} < \displaystyle\frac{\varphi_j+\varphi_i}{2}+\pi \ \ \ \mbox{moving away,}\\ \tilde{v}_{ji} < 0 & \mbox{for} & \displaystyle\frac{\varphi_j+\varphi_i}{2}+\pi < \alpha_{ji} < \displaystyle\frac{\varphi_j+\varphi_i}{2}+2\pi \ \ \ \mbox{approaching.} \end{array} \end{eqnarray} \tag{ 2.6 }$$

Hence, a particle located in the half-sphere in a clockwise direction from the mean angle (φ_j + φ_i)/2 approaches the focal particle i, whereas particles located in the other half-sphere anticlockwise from the mean angle are moving away (see figure 1). Please note that for φ_j = φ_i the social force vanishes as $\tilde v_{ji}=0$ according to (2.5).

Zoom In Zoom Out Reset image size

Four different regions in the (μ_m,μ_a)-parameter space are distinguished [29]:

• pure attraction: μ_m > 0, μ_a > 0;
• effective alignment: μ_m > 0, μ_a < 0, i.e. attraction to particles moving away, repulsion from particles coming closer;
• pure repulsion: μ_m < 0, μ_a < 0;
• effective anti-alignment: μ_m < 0, μ_a > 0, i.e. attraction to particles coming closer, repulsion from particles moving away.

A schematic visualization of the turning behaviour of interacting agents in the different regimes is shown in figure 2. Some typical spatial snapshots obtained from individual-based simulations of (2.1) are shown in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Please note that in previous publications [28, 29], 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 [28] and the respective response of a focal agent to neighbours irrespective of their behaviour. Furthermore, an asymmetric alignment response for μ_m > 0, μ_a < 0, 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 [28, 29] 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.

3.1. 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 PDF

$$\begin{eqnarray} &&P_N(\boldsymbol{r}_1,\varphi_1; \boldsymbol{r}_2,\varphi_2; \ldots ; \boldsymbol{r}_N,\varphi_N ; t) \end{eqnarray} \tag{ 3.1 }$$

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:

$$\begin{eqnarray} &&\left [\; \prod_{j = 1 }^N \left ( \int {\rm d}^2 r_j \int_0^{2 \pi} {\rm d} \varphi_j \right) \right ] P_N(\boldsymbol{r}_1,\varphi_1; \boldsymbol{r}_2,\varphi_2; \ldots ; \boldsymbol{r}_N,\varphi_N ; t) = 1. \end{eqnarray} \tag{ 3.2 }$$

In agreement with (2.1), we can write down the FPE for the dynamics of the PDF as follows:

$$\begin{eqnarray} &&\frac{\partial P_N}{\partial t} = - s_0 \sum_{i=1}^N \nabla_{\boldsymbol{r}_i} \cdot \left ( \boldsymbol{e}_{i,h} P_N \right) - \sum_{i=1}^N \frac{\partial}{\partial \varphi_i} \left ( \frac{ F_{i,\varphi} \, P_N }{s_0} \right) + \frac{D_{\varphi}}{s_0^2} \sum_{i=1}^N \frac{\partial^2 P_N}{\partial \varphi_i^2}. \end{eqnarray} \tag{ 3.3 }$$

From the linear FPE above one can derive an evolution equation for the marginal PDF

$$\begin{eqnarray} &&\fl P(\boldsymbol{r}_i,\varphi_i,t) = \left [ \; \prod_{j \neq i}^N \left ( \int {\rm d}^2 r_j \int_0^{2 \pi} {\rm d} \varphi_j \right) \right ] P_N(\boldsymbol{r}_1,\varphi_1; \boldsymbol{r}_2,\varphi_2; \ldots ; \boldsymbol{r}_N,\varphi_N ; t) \end{eqnarray} \tag{ 3.4 }$$

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.

$$\begin{eqnarray} &&P_N(\boldsymbol{r}_1,\varphi_1; \boldsymbol{r}_2,\varphi_2; \ldots ; \boldsymbol{r}_N,\varphi_N ; t) = \prod_{i=1}^N P(\boldsymbol{r}_i,\varphi_i,t). \end{eqnarray} \tag{ 3.5 }$$

By this means, one obtains an effective one-particle description

$$\begin{eqnarray} &&\frac{\partial P (\boldsymbol{r}_i,\varphi_i,t)}{\partial t} = - s_0 \nabla_{\boldsymbol{r}_i} \cdot \left (\boldsymbol{e}_{i,h} P \right) - \frac{N-1}{s_0} \frac{\partial }{\partial \varphi_i} \left ( F_{\varphi}\, P \right) + \frac{D_\varphi}{s_0^2} \frac{\partial^2 P}{\partial \varphi_i^2}, \end{eqnarray} \tag{ 3.6 }$$

where the force F_φ is given by the following integral over P(r_i,φ_i,t):

$$\begin{eqnarray} &&\fl F_{\varphi}(\boldsymbol{r}_i,\varphi_i,t) = 2 s_0 \int_{\varphi_i}^{\varphi_i + 2 \pi} \mbox{d}\varphi_j \, \sin \left ( \frac{\varphi_i - \varphi_j}{2} \right)\int_{l_{\rm c}}^{l_{\rm s}} \mbox{d}r_{ji} \, r_{ji} \nonumber \\ &&\times \left [ \mu_{\rm m}(r_{ji}) \int_{\frac{\varphi_i + \varphi_j}{2}}^{\frac{\varphi_i + \varphi_j}{2} + \pi} \mbox{d}\alpha_{ji} \, \sin(\alpha_{ji} - \varphi_i) \sin \left ( \frac{\varphi_i + \varphi_j}{2} - \alpha_{ji} \right) P(\boldsymbol{r}_i+\boldsymbol{r}_{ji},\varphi_j,t) \right. \nonumber \\ &&\left. - \mu_{\rm a}(r_{ji}) \int_{\frac{\varphi_i + \varphi_j}{2}+\pi}^{\frac{\varphi_i + \varphi_j}{2} + 2 \pi} \mbox{d}\alpha_{ji} \, \sin(\alpha_{ji} - \varphi_i) \sin \left ( \frac{\varphi_i + \varphi_j}{2} - \alpha_{ji} \right) P(\boldsymbol{r}_i+\boldsymbol{r}_{ji},\varphi_j,t) \right] \nonumber \\ &&- \int_{\varphi_i}^{\varphi_i+2\pi}\mbox{d}\varphi_j \int_{0}^{l_{\rm c}} \mbox{d}r_{ji} \, r_{ji} \int_{0}^{2\pi} \mbox{d}\alpha_{ji} \, \mu_r(r_{ji}) \sin(\alpha_{ji} - \varphi_i) P(\boldsymbol{r}_i+\boldsymbol{r}_{ji},\varphi_j,t). \end{eqnarray} \tag{ 3.7 }$$

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 FPE

$$\begin{eqnarray} &&\frac{\partial p (\boldsymbol{r},\varphi,t)}{\partial t} = - s_0 \nabla_{\boldsymbol{r}} \cdot \left (\boldsymbol{e}_h \, p \right) - \frac{\partial }{\partial \varphi} \frac{ F_{\varphi} (\boldsymbol{r},\varphi,t) \, p }{s_0} + \frac{D_\varphi}{s_0^2} \frac{\partial^2 p}{\partial \varphi^2}. \end{eqnarray} \tag{ 3.8 }$$

First, we assume a spatially homogeneous situation where p(r,φ,t) = p(φ,t) holds. In this case, (3.8) is reduced to

$$\begin{eqnarray} &&\fl \frac{\partial p (\varphi,t)}{\partial t} = - \frac{D_{\varphi}}{s_0^2} \frac{\partial }{\partial \varphi} \left [ \frac{\kappa}{\rho_0} \, \, \int_{\varphi}^{\varphi + 2 \pi} {\rm d} \varphi_j \, \sin (\varphi_j - \varphi) p(\varphi_j,t) \, p(\varphi,t) - \frac{\partial p}{\partial \varphi} \right], \end{eqnarray} \tag{ 3.9 }$$

where the dimensionless coupling parameter

$$\begin{eqnarray} &&\kappa = \frac{\pi \rho_0 s_0^2}{2 D_{\varphi}} \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji} \left ( \mu_{\rm m}(r_{ji}) - \mu_{\rm a}(r_{ji}) \right) \end{eqnarray} \tag{ 3.10 }$$

is introduced. Here, ρ₀ = N/L² is the homogeneous particle density, with L being the linear spatial dimension of the two-dimensional square-shaped system.

One obtains the same FPE 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 velocity-alignment as considered in [38] (continuum time form of the original Vicsek model [13]), 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 2(D)).

The spatially homogeneous, time-independent particle density

$$\begin{eqnarray} &&p^{(\Phi)}(\varphi) = \frac{\rho_0}{2 \pi} \, \frac{\exp \left ( \kappa \Phi \cos (\varphi - \varphi_0) \right)}{I_0(\kappa \Phi)} \end{eqnarray} \tag{ 3.11 }$$

is a solution to both nonlinear FPEs (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 φ₀ defines the x-coordinate, i.e. φ₀ = 0. From the definition of the polar order parameter

$$\begin{eqnarray} &&\Phi = \left | \left \langle {\rm e}^{{\rm i} \varphi} \right \rangle \right | = \left | \int_0^{2\pi}\,{\rm d}\varphi \, p^{(\Phi)}(\varphi) {\rm e}^{{\rm i} \varphi} \right |, \end{eqnarray} \tag{ 3.12 }$$

we can determine the order parameter by solving the transcendental equation

$$\begin{eqnarray} &&\Phi = \frac{I_1(\kappa \Phi)}{I_0(\kappa \Phi)}. \end{eqnarray} \tag{ 3.13 }$$

Equation (3.13) is always fulfilled for Φ = 0. In that case, (3.11) yields the uniform distribution (spatially homogeneous, disordered state), i.e.

$$\begin{eqnarray} &&p^{(0)} = \frac{\rho_0}{2 \pi}. \end{eqnarray} \tag{ 3.14 }$$

If κ > 2 (high effective-alignment strength μ_m − μ_a and low noise D_φ, respectively), there is a second non-trivial 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⁽⁰⁾ is stable for κ < 2 [36, 37].

The critical line κ = 2 for the order–disorder transition is translated to microscopic model parameters as follows:

$$\begin{eqnarray} &&\int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji} \left( \mu_{\rm m}(r_{ji}) - \mu_{\rm a}(r_{ji}) \right) = \frac{4 D_{\varphi}}{\pi \rho_0 s_0^2}. \end{eqnarray} \tag{ 3.15 }$$

The above is a generalization of the previous result obtained in [29]. Since the order–disorder 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:

$$\begin{eqnarray} &&\Phi(\kappa) \simeq \sqrt{\kappa - 2} + \mathcal{O}\left ((k-2)^{3/2} \right ). \end{eqnarray} \tag{ 3.16 }$$

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.

3.2. Stability analysis in Fourier space

In this section it is shown how the stability of solutions of the nonlinear FPE (3.8) can be analysed in Fourier space. A similar approach was used by Chou et al [30] in order to achieve a kinetic theory for self-propelled particles with metric-free interactions.

In order to analyse the FPE (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

$$\begin{eqnarray} &&{\hat{g}}_{n}(\boldsymbol{k},t) = \int {\rm d}^2 r \int_0^{2 \pi} {\rm d} \varphi \, p(\boldsymbol{r},\varphi,t) {\rm e}^{{\rm i} \boldsymbol{k}\cdot\boldsymbol{r} + \mathrm{i} n \varphi } \end{eqnarray} \tag{ 3.17 }$$

reads

$$\begin{eqnarray} &&\fl \frac{\partial {\hat{g}}_{n}(\boldsymbol{k},t)}{\partial t} - \frac{{\rm i} k_x s_0}{2} \left ( {\hat{g}}_{n+1} + {\hat{g}}_{n-1} \right) - \frac{k_y s_0}{2} \left ({\hat{g}}_{n+1} - {\hat{g}}_{n-1} \right) = -n^2 \frac{D_{\varphi}}{s_0^2} {\hat{g}}_{n} \nonumber \\ &&+\, \frac{{\rm i} n }{(2 \pi)^3 s_0} \sum_{j=-\infty}^{\infty} \int {\rm d}^2 q \, {\hat{g}}_{j}(\boldsymbol{q},t) {\hat{K}}_{n-j}(\boldsymbol{k}-\boldsymbol{q},t), \end{eqnarray} \tag{ 3.18 }$$

where the Fourier coefficients of the force (3.7) are defined as

$$\begin{eqnarray} &&{\hat{K}}_{n}(\boldsymbol{k},t) = \int {\rm d}^2 r \int_0^{2 \pi} {\rm d} \varphi \, F_{\varphi}(\boldsymbol{r},\varphi,t) {\rm e}^{{\rm i} \boldsymbol{k}\cdot\boldsymbol{r} + i n \varphi}. \end{eqnarray} \tag{ 3.19 }$$

Obviously, the FPE turns into an infinite hierarchy of equations (3.18) in Fourier space.

For the stability analysis, the Fourier coefficients ${\hat {K}}_{n}(\boldsymbol {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 coefficients ${\hat {K}}_{n}(\boldsymbol {k},t)$ of the force can be written as a linear combination of the Fourier coefficients ${\hat {g}}_{r}(\boldsymbol {k},t)$

$$\begin{eqnarray} &&{\hat{K}}_{n}(\boldsymbol{k},t) = \sum_{r=-\infty}^{\infty} \tilde{K}_{n,r}(\boldsymbol{k}) \, {\hat{g}}_{r}(\boldsymbol{k},t) \end{eqnarray} \tag{ 3.20 }$$

using an infinite-dimensional matrix $\tilde {K}_{n,r}(\boldsymbol {k})$. Details on the calculation of the matrix elements are given in appendix A. Consequently, the dynamics of the Fourier coefficients ${\hat {g}}_{n}(\boldsymbol {k},t)$ read

$$\begin{eqnarray} &&\fl \frac{\partial {\hat{g}}_{n}(\boldsymbol{k},t)}{\partial t} - \frac{{\rm i} k_x s_0}{2} \left ( {\hat{g}}_{n+1} + {\hat{g}}_{n-1} \right) - \frac{k_y s_0}{2} \left ({\hat{g}}_{n+1} - {\hat{g}}_{n-1} \right) = -n^2 \frac{D_{\varphi}}{s_0^2} {\hat{g}}_{n} \nonumber \\ &&+\, \frac{{{\rm i}} n }{(2 \pi)^3 s_0} \sum_{j=-\infty}^{\infty} \sum_{r=-\infty}^{\infty} \int {\rm d}^2 q \, {\hat{g}}_{j}(\boldsymbol{q},t) {\hat{g}}_{r}(\boldsymbol{k} - \boldsymbol{q},t) \tilde{K}_{n-j,r}(\boldsymbol{k} - \boldsymbol{q}). \end{eqnarray} \tag{ 3.21 }$$

Let ${\hat {g}}_{n}^{(0)}(\boldsymbol {k},t)$ be a solution of the equation above. The dynamics of a small perturbation $\delta {\hat {g}}_{n}(\boldsymbol {k},t) = {\hat {g}}_{n}(\boldsymbol {k},t) - {\hat {g}}_{n}^{(0)}(\boldsymbol {k},t)$ is determined by the linearized equation

$$\begin{eqnarray} &&\fl \frac{\partial \delta {\hat{g}}_{n}(\boldsymbol{k},t)}{\partial t} - \frac{{\rm i} k_x s_0}{2} \left ( \delta {\hat{g}}_{n+1} + \delta {\hat{g}}_{n-1} \right) - \frac{k_y s_0}{2} \left (\delta {\hat{g}}_{n+1} - \delta {\hat{g}}_{n-1} \right) = -n^2 \frac{D_{\varphi}}{s_0^2} \delta {\hat{g}}_{n} + \frac{{\rm i} n }{(2 \pi)^3 s_0} \nonumber \\ &&\times \sum_{j=-\infty}^{\infty} \sum_{r=-\infty}^{\infty} \int {\rm d}^2 q \, \left (\delta {\hat{g}}_{j}(\boldsymbol{q},t) {\hat{g}}_{r}^{(0)}(\boldsymbol{k}-\boldsymbol{q},t) + {\hat{g}}_{j}^{(0)}(\boldsymbol{q},t) \delta {\hat{g}}_{r}(\boldsymbol{k}-\boldsymbol{q},t) \right )\nonumber\\ &&\times \tilde{K}_{n-j,r}(\boldsymbol{k} - \boldsymbol{q}). \end{eqnarray} \tag{ 3.22 }$$

The stability analysis of the solutions

$$\begin{eqnarray} &&p^{(\Phi)}(\varphi) = \frac{\rho_0}{2 \pi} \, \frac{\exp \left ( \kappa \Phi \cos \varphi \right)}{I_0(\kappa \Phi)} \end{eqnarray} \tag{ 3.23 }$$

requires the corresponding Fourier coefficients

$$\begin{eqnarray} &&{\hat{g}}_{n}(\boldsymbol{k}; \Phi) =(2 \pi)^2 \rho_0 \, \frac{I_n(\kappa \Phi)}{I_0(\kappa \Phi)} \, \delta (\boldsymbol{k}). \end{eqnarray} \tag{ 3.24 }$$

Inserting those in (3.22) yields the following linearized ordinary differential equations:

$$\begin{eqnarray} &&\frac{\partial \delta {\hat{g}}_{n}(\boldsymbol{k},t)}{\partial t} = \sum_{r=-\infty}^{\infty} \tilde{M}_{n,r}(\boldsymbol{k}) \, \delta {\hat{g}}_{r}(\boldsymbol{k},t) \end{eqnarray} \tag{ 3.25 }$$

with the matrix

$$\begin{eqnarray} &&\fl \tilde{M}_{n,r}(\boldsymbol{k}) = \frac{{\rm i} k_x s_0}{2} \left ( \delta_{n+1,r} + \delta_{n-1,r} \right) + \frac{k_y s_0}{2} \left ( \delta_{n+1,r} - \delta_{n-1,r} \right) - n^2 \frac{D_{\varphi}}{s_0^2} \, \delta_{n,r} \nonumber \\ &&+\, \frac{{\rm i} n \rho_0}{2 \pi s_0} \sum_{j=-\infty}^{\infty} \frac{I_j(\kappa \Phi)}{I_0(\kappa \Phi)} \left ( \tilde{K}_{n-r,j}(0) + \tilde{K}_{n-j,r}(\boldsymbol{k}) \right). \end{eqnarray} \tag{ 3.26 }$$

The eigenvalues σ of the matrix above, given by the solutions of $\det (\tilde M -\sigma \mathbb {I})=0$, 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 $\tilde {n}{>}0$ exists, such that $\delta {\hat {g}}_{n} = 0$ for $\left | n \right |{>}\tilde {n}$ holds approximatively (cf [30]). That assumption is reasonable, because the nth Fourier coefficient is strongly damped, since

$$\begin{eqnarray} &&\frac{\partial \delta {\hat{g}}_{n}(\boldsymbol{k},t)}{\partial t} \propto - n^2 \frac{D_{\varphi}}{s_0^2} \, \delta {\hat{g}}_{n}(\boldsymbol{k},t). \end{eqnarray} \tag{ 3.27 }$$

For the numerical analysis, $\tilde {n}=50$ is assumed. Using this approximation, the stability of the spatially homogeneous states can be studied. Please note that only two approximations were used: firstly, correlations between particles are neglected (3.5) (mean-field); secondly, the truncation of the hierarchy of equations (3.25). In principle, it is possible to analyse 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. Thus, possible instabilities in the ordered state correspond to the emergence of band-like structures perpendicular to the direction of motion, commonly observed in related systems (see figure 2(B)).

The illustrated method allows us 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 p⁽⁰⁾ = ρ₀/(2π) is investigated (Φ = 0), whereas the solution p^(Φ)(φ) (Φ > 0) describing polar order is inserted in (3.26) for κ > 2. Figures 4 and 5 show the results of the stability analysis in the (μ_m,μ_a)-parameter space and the dispersion relations for specific points (named (A)–(F) in figures 4 and 5), respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The kinetic approach reveals instabilities nearby the points (A) and (B). In these parameter regions, an interval of wavenumbers becomes unstable (cf figure 5). 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⁽⁰⁾ is found in (C). According to the individual-based simulations, a clustering phase emerges nearby the point (C) (cf figure 3(E) 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 5(D). According to figures 5(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 4, 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 elements $\tilde {K}_{n,r}(\boldsymbol {k})$ involve Bessel functions of the first kind of k (see appendix A and [30]), which are oscillating functions with alternating Taylor coefficients [31]. The approximation

$$\begin{eqnarray} &&J_{\nu,N}(k x) \approx \sum_{m=0}^N \frac{(-1)^m}{m! \, \Gamma(m+\nu+1)} {\left( \frac{k x}{2} \right)}^{2m+\nu} \end{eqnarray} \tag{ 3.28 }$$

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.

4.1. 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:

$$\begin{equation} p(\boldsymbol{r},\varphi,t) = \frac{1}{2 \pi} \sum_{n=-\infty}^{\infty} {\skew6\hat{f}}_{n} (\boldsymbol{r},t) \, {\rm e}^{-{\rm i} n \varphi}, \end{equation} \tag{ 4.1a }$$

$$\begin{equation} F_{\varphi}(\boldsymbol{r},\varphi,t) = \frac{1}{2 \pi} \sum_{n=-\infty}^{\infty} {\hat{Q}}_{n} (\boldsymbol{r},t) \, {\rm e}^{-{\rm i} n \varphi}. \end{equation} \tag{ 4.1b }$$

The FPE in Fourier space is obtained by multiplying (3.8) with e^inφ and integration over $\varphi \in \left [0,2 \pi \right ]$:

$$\begin{eqnarray} &&\frac{\partial {\skew6\hat{f}}_{n}}{\partial t} + s_0 \left ( {\underset{\sim}{\nabla}} {\skew6\hat{f}}_{n-1} + {\underset{\sim}{\nabla}}^{*} {\skew6\hat{f}}_{n+1} \right) = \frac{1}{2 \pi} \frac{{\rm i} n}{s_0} \sum_{j=-\infty}^{\infty} {\skew6\hat{f}}_{j} \, {\hat{Q}}_{n-j} - n^2 \frac{D_{\varphi}}{s_0^2} {\skew6\hat{f}}_{n} \, \mbox{. } \end{eqnarray} \tag{ 4.2 }$$

In (4.2), the complex derivative ${\underset{\sim}{\nabla}} = 0.5\left (\partial _x + \mathrm {i} \partial _y \right )$ was introduced. The Fourier coefficients

$$\begin{eqnarray} &&{\skew6\hat{f}}_{n}(\boldsymbol{r},t) = \int_{0}^{2 \pi} {\rm d} \varphi \, p(\boldsymbol{r},\varphi,t) {\rm e}^{{\rm i} n \varphi } = \rho(\boldsymbol{r},t) \left \langle {\rm e}^{{\rm i} n \varphi } \right \rangle \end{eqnarray} \tag{ 4.3 }$$

obey the symmetry ${\skew6\hat{f}}_{n}^{*}(\boldsymbol {r},t) = {\skew6\hat{f}}_{-n}(\boldsymbol {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

$$\begin{equation} \rho(\boldsymbol{r},t) = {\skew6\hat{f}}_{0}(\boldsymbol{r},t), \end{equation} \tag{ 4.4a }$$

$$\begin{equation} w_x(\boldsymbol{r},t) = \frac{s_0}{2} \left ( {\skew6\hat{f}}_{1}(\boldsymbol{r},t) + {\skew6\hat{f}}_{-1}(\boldsymbol{r},t) \right), \end{equation} \tag{ 4.4b }$$

$$\begin{equation} w_y(\boldsymbol{r},t) = \frac{s_0}{2i} \left ( {\skew6\hat{f}}_{1}(\boldsymbol{r},t) - {\skew6\hat{f}}_{-1}(\boldsymbol{r},t) \right). \end{equation} \tag{ 4.4c }$$

The second Fourier coefficients are related to the nematic order parameter $\left | \left \langle {\mathrm { e}}^{2\mathrm {i}\varphi } \right \rangle \right |$ and the symmetric temperature tensor, as defined in [24, 26, 29], which is a measure for the width of the velocity distribution (see appendix C 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 force ${\hat {Q}}_{n}(\boldsymbol {r},t)$ are approximatively calculated by expanding p(r + r_ji,φ,t) into a multidimensional Taylor series for small $\left | \boldsymbol {r}_{ji} \right |$ [32], 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 FPE 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. $\left | \boldsymbol {w} \right |/(s_0 \rho _0) \propto \epsilon$, 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 $\left | {\skew6\hat{f}}_{1} \right |$ is larger than $\left | {\skew6\hat{f}}_{n} \right |$ with $\left | n \right |{>}1$, because of the damping term in (4.2) proportional to n², leading to the scaling relations

$$\begin{eqnarray} &&\rho(\boldsymbol{r},t) - \rho_0 \propto \epsilon, \quad {\underset{\sim}{\nabla}} \propto \epsilon, \quad \partial_t \propto \epsilon, \quad {\skew6\hat{f}}_{n} \propto \epsilon^{\left | n \right |}. \end{eqnarray} \tag{ 4.5 }$$

In [21] 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 [18, 21, 23, 32] and in [22] 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 ³.

$$\begin{equation} \fl \frac{\partial \rho}{\partial t} = - s_0 \left ( {\underset{\sim}{\nabla}} {\skew6\hat{f}}_{1}^* + {\underset{\sim}{\nabla}}^{*} {\skew6\hat{f}}_{1} \right), \end{equation} \tag{ 4.6a }$$

$$\begin{eqnarray} &&\fl \frac{\partial {\skew6\hat{f}}_{1}}{\partial t} = \left (\xi_1 \rho - \xi_5 \right) {\skew6\hat{f}}_{1} - \xi_1 {\skew6\hat{f}}_{2} {\skew6\hat{f}}_{1}^{*} + \left ( \frac{32}{45} \xi_3 \rho_0 - s_0 \right) {\underset{\sim}{\nabla}}^* {\skew6\hat{f}}_{2} + \left ( \frac{32}{9} \xi_3 \rho - \xi_2 \rho - s_0 \right) {\underset{\sim}{\nabla}} \rho \nonumber \\ &&- \frac{64}{45} \xi_3 {\skew6\hat{f}}_{1}^* {\underset{\sim}{\nabla}} {\skew6\hat{f}}_{1} + \frac{\xi_4 \rho_0}{2} \left ( \Delta {\skew6\hat{f}}_{1} + 2 {\underset{\sim}{\nabla}}^2 {\skew6\hat{f}}_{1}^* \right), \end{eqnarray} \tag{ 4.6b }$$

$$\begin{equation} \fl \frac{\partial {\skew6\hat{f}}_{2}}{\partial t} = 2 \left [ -2 \xi_5 {\skew6\hat{f}}_{2} + \xi_1 {\skew6\hat{f}}_{1}^2 - \left ( \frac{s_0}{2} + \frac{64}{45} \xi_3 \rho \right) {\underset{\sim}{\nabla}} {\skew6\hat{f}}_{1} + \left ( \frac{32}{9} \xi_3 - \xi_2 \right) {\skew6\hat{f}}_{1} {\underset{\sim}{\nabla}} \rho - \xi_4 \rho_0 {\underset{\sim}{\nabla}}^2 \rho \right ]. \end{equation} \tag{ 4.6c }$$

The following coefficients are introduced:

$$\begin{equation} \xi_1 = \frac{\pi}{4} \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji} \left [ \mu_{\rm m}(r_{ji}) - \mu_{\rm a}(r_{ji}) \right ], \end{equation} \tag{ 4.7a }$$

$$\begin{equation} \xi_2 = \frac{\pi}{s_0} \int_{0}^{l_{\rm c}} {\rm d} r_{ji} \, r_{ji}^2 \, \mu_r(r_{ji}), \end{equation} \tag{ 4.7b }$$

$$\begin{equation} \xi_3 = \frac{1}{\pi} \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji}^2 \left [ \mu_{\rm m}(r_{ji}) + \mu_{\rm a}(r_{ji}) \right ], \end{equation} \tag{ 4.7c }$$

$$\begin{equation} \xi_4 = \frac{\pi}{8} \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji}^3 \left [ \mu_{\rm m}(r_{ji}) - \mu_{\rm a}(r_{ji}) \right ], \end{equation} \tag{ 4.7d }$$

$$\begin{equation} \xi_5 = \frac{D_{\varphi}}{s_0^2}. \end{equation} \tag{ 4.7e }$$

Equation (4.6a) is simply the continuity equation representing the conservation of the particle number N. Please note the difference in dependence of the transport coefficients in (4.6), due to a consistent application of the scaling ansatz ρ(r,t) − ρ₀∝. Some depend on the global density ρ₀ whereas others depend on the local density ρ(r,t).

4.2. Hydrodynamic limit

Another simplification of the system (4.6) is reached by adiabatically eliminating ${\skew6\hat{f}}_{2}$, i.e. $\partial _t {\skew6\hat{f}}_{2} \approx 0$, where

$$\begin{eqnarray} &&{\skew6\hat{f}}_{2} \approx \frac{1}{2 \xi_5} \left [ \xi_1 {\skew6\hat{f}}_{1}^2 - \left ( \frac{s_0}{2} + \frac{64}{45} \xi_3 \rho_0 \right) {\underset{\sim}{\nabla}} {\skew6\hat{f}}_{1} \right ] \end{eqnarray} \tag{ 4.8 }$$

follows. This approximation is only valid if the first Fourier mode relaxes considerably more slowly than ${\skew6\hat{f}}_{2}$. The assumption of a time scale separation between the first and second Fourier modes can only be justified if $\left | \xi _1 \rho -\xi _5 \right |\ll 4\xi _5$, which is equivalent to −6 ≪ κ ≪ 10. Furthermore, ξ₅∝D_φ > 0 has to hold.

Please note that all terms of the order ³ 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:

$$\begin{equation} \fl\frac{\partial \rho}{\partial t} = - \nabla \cdot \boldsymbol{w} \end{equation} \tag{ 4.9a }$$

$$\begin{eqnarray} &&\fl \frac{\partial \boldsymbol{w}}{\partial t} = \left ( \lambda_1 - \eta_1 \left | \boldsymbol{w} \right |^2 \right)\boldsymbol{w} + \lambda_2 \nabla \rho + \lambda_4 \, \Delta \boldsymbol{w} + \lambda_5 \nabla (\nabla \cdot \boldsymbol{w}) + \eta_2 (\boldsymbol{w} \cdot \nabla) \, \boldsymbol{w} + \eta_3 ( \textstyle\frac{1}{2} \nabla \left | \boldsymbol{w} \right |^2 - \boldsymbol{w} \, (\nabla \cdot \boldsymbol{w}) ).\nonumber\\ \end{eqnarray} \tag{ 4.9b }$$

The transport coefficients {λ_i,η_i} as functions of ξ_i (4.7) read

$$\begin{equation} \lambda_1 = \xi_1 \rho - \xi_5, \end{equation} \tag{ 4.10a }$$

$$\begin{equation} \lambda_2 = \frac{s_0}{2} \left ( \frac{32}{9} \xi_3 \rho - \xi_2 \rho - s_0 \right), \end{equation} \tag{ 4.10b }$$

$$\begin{equation} \lambda_4 = \frac{1}{4} \left [ \frac{s_0 - \frac{32}{45} \xi_3 \rho_0}{2 \xi_5} \left ( \frac{s_0}{2} + \frac{64}{45} \xi_3 \rho_0 \right) + \xi_4 \rho_0 \right ], \end{equation} \tag{ 4.10c }$$

$$\begin{equation} \lambda_5 = \frac{\xi_4 \rho_0}{2}, \end{equation} \tag{ 4.10d }$$

$$\begin{equation} \eta_1 = \frac{\xi_1^2}{2 \xi_5 s_0^2}, \end{equation} \tag{ 4.10e }$$

$$\begin{equation} \eta_2 = \frac{3 \xi_1}{8 s_0 \xi_5} \left ( \frac{256 \xi_3 \rho_0}{135} - s_0 \right) - \frac{32 \xi_3}{45 s_0}, \end{equation} \tag{ 4.10f }$$

$$\begin{equation} \eta_3 = \frac{5 \xi_1}{8 \xi_5} - \frac{32 \xi_3}{45 s_0}. \end{equation} \tag{ 4.10g }$$

The hydrodynamic equations (4.9) are comparable in structure to the theory of Toner and Tu [14, 15].

In the theory of Toner and Tu, the coefficients η₂ and η₃ are assumed to be non-zero phenomenological parameters due to the absence of Galilean invariance [15]. Indeed, in our model, η₂ and η₃ can take both negative and positive values. Likewise, λ₂ can be positive and negative. A positive λ₂-term leads to a flow into the direction of higher densities and consequently to clustering, as shown below. The damping coefficient η₁ is positive for all microscopic parameters, consistently. The sign of the different transport coefficients is discussed in appendix B in detail.

The stability of the hydrodynamic equations (4.9) requires that λ₄ and λ₅ 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 λ₄ and λ₅ may change signs. As shown in [3], a negative λ₄ 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 separation of the first and second Fourier modes is only justified for λ₁ ≈ 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 analysing the fixed points of

$$\begin{eqnarray} &&\frac{{\rm d} \boldsymbol{w}}{{\rm d} t} = \left ( \lambda_1 - \eta_1 \left | \boldsymbol{w} \right |^2 \right)\boldsymbol{w}. \end{eqnarray} \tag{ 4.11 }$$

That is the standard symmetry-breaking term present in all models of collective motion [3, 15, 19–21, 23, 32]. Apparently, two spatially homogeneous solutions exist: Firstly, w₀ = 0 corresponding to a disordered phase and vanishing centre of mass velocity. Secondly,

$$\begin{eqnarray} &&\boldsymbol{w}_1 = \sqrt{\frac{\lambda_1}{\eta_1}} \, \boldsymbol{e}, \end{eqnarray} \tag{ 4.12 }$$

where e denotes an arbitrary unit vector reflecting the isotropy of the system. Without loss of generality we set e = e_x. For λ₁ < 0, only the fixed point w₀ exists and is stable against spatially homogeneous perturbations. For λ₁ ⩾ 0, w₁ is a second fixed point corresponding to an ordered phase with non-zero mean-speed (swarming phase).

The supercritical pitchfork bifurcation in λ₁ = 0, i.e.

$$\begin{eqnarray} &&\int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji} \left( \mu_{\rm m}(r_{ji}) - \mu_{\rm a}(r_{ji}) \right) = \frac{4 D_{\varphi}}{\pi \rho_0 s_0^2} \, \mbox{ } \end{eqnarray} \tag{ 4.13 }$$

corresponds to the order–disorder transition (3.16) described in section 3.1. By introducing the dimensionless coupling strength κ (3.10) again,

$$\begin{eqnarray} &&\kappa = \frac{\pi \rho_0 s_0^2}{2 D_{\varphi}} \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji} \left ( \mu_{\rm m}(r_{ji}) - \mu_{\rm a}(r_{ji}) \right), \end{eqnarray} \tag{ 4.14 }$$

the polar order parameter $\Phi _{\mathrm { H}} = \left | \boldsymbol {w}_1 \right |/(s_0 \rho _0)$ as described by the hydrodynamic theory can be written as follows⁵:

$$\begin{eqnarray} &&\Phi_{\rm H} = 2 \, \frac{\sqrt{\kappa -2}}{\kappa}. \end{eqnarray} \tag{ 4.15 }$$

Close to the critical point κ_c = 2, one recovers the correct scaling behaviour $\Phi _{\mathrm { H}} \simeq \sqrt {\kappa -2}$ (3.16). However, for large κ the polar order parameter Φ_H tends to zero according to (4.15) as shown in figure 6. 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.

Zoom In Zoom Out Reset image size

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.

4.3. Stability of the spatially homogeneous solutions

We proceed with the analysis of homogeneous solutions against spatial perturbations. At first, let us consider the stability of the disordered, homogeneous solution with w₀ = 0. In order to investigate the stability of w₀, the ansatz

$$\begin{equation} \rho(\boldsymbol{r},t) = \rho_0 + \delta \rho(\boldsymbol{r},t), \end{equation} \tag{ 4.16a }$$

$$\begin{equation} \boldsymbol{w}(\boldsymbol{r},t) = \delta \boldsymbol{w}(\boldsymbol{r},t) \, \end{equation} \tag{ 4.16b }$$

is inserted in (4.9) and the resulting equations are linearized in the perturbations which are assumed to be small. The linearized equations read

$$\begin{equation} \frac{\partial \delta \rho}{\partial t} = - \nabla \cdot \delta \boldsymbol{w}, \end{equation} \tag{ 4.17a }$$

$$\begin{equation} \frac{\partial \delta \boldsymbol{w}}{\partial t} = \lambda_1 \delta \boldsymbol{w} + \lambda_2 \nabla \delta \rho + \lambda_4 \, \Delta \delta \boldsymbol{w} + \lambda_5 \nabla (\nabla \cdot \delta \boldsymbol{w}). \end{equation} \tag{ 4.17b }$$

Inserting the exponential ansatz

$$\begin{equation} \delta \rho(\boldsymbol{r},t) = \zeta_{\rho} \exp \left (\sigma t + {\rm i} \boldsymbol{k}\cdot \boldsymbol{r} \right), \end{equation} \tag{ 4.18a }$$

$$\begin{equation} \delta \boldsymbol{w}(\boldsymbol{r},t) = \boldsymbol{\zeta}_{\boldsymbol{w}} \exp \left (\sigma t + {\rm i} \boldsymbol{k}\cdot \boldsymbol{r} \right) \end{equation} \tag{ 4.18b }$$

yields the dispersion relation σ(k). Since w₀ is isotropic, the eigenvalues σ are independent of the direction of the wavevector:

$$\begin{equation} \sigma_1 = \lambda_1 - \lambda_4 k^2, \end{equation} \tag{ 4.19a }$$

$$\begin{equation} \sigma_{2,3} = \frac{1}{2} \left [ \lambda_1 - \left ( \lambda_4 + \lambda_5 \right) k^2 \pm \sqrt{\left [\lambda_1 - (\lambda_4+\lambda_5) \, k^2 \right]^2 + 4 k^2 \lambda_2 } \right ]. \end{equation} \tag{ 4.19b }$$

Provided all eigenvalues are less than zero, the disordered spatially homogeneous state is stable. The first eigenvalue σ₁ is equal to λ₁ for k = 0, i.e. the stability of w₀ is determined by the sign of λ₁. That is the destabilization of the homogeneous state as discussed before. For λ₄ > 0, the first eigenvalue σ₁ is monotonically decreasing. However, the eigenvalue σ₁ is monotonically increasing for λ₄ < 0. In this case, the hydrodynamic equations (4.9) lose their validity. This unphysical behaviour is due to the restriction to second derivatives when the integral (3.7) was calculated. Furthermore, it is problematic 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₀ by analysing those eigenvalues, which tend to zero for small wavenumbers (hydrodynamic modes [15]). For small wavenumbers, the eigenvalue σ₂ (σ₂ > σ₃) is expanded in a Taylor series:

$$\begin{equation} \sigma_2 \simeq - \frac{\lambda_2}{\lambda_1} \, k^2 + \frac{\lambda_2^2 - \lambda_1 \lambda_2 \left ( \lambda_4 + \lambda_5 \right) }{\lambda_1^3} \, k^4 + \mathcal{O}(k^6). \end{equation} \tag{ 4.20 }$$

Close to the order–disorder transition line, λ₁ is approximately zero. Therefore it is sufficient to keep the leading order terms in 1/λ₁ [19]:

$$\begin{equation} \sigma_2 \simeq - \frac{\lambda_2}{\lambda_1} \, k^2 + \frac{\lambda_2^2}{\lambda_1^3} \, k^4 + \mathcal{O}(k^6). \end{equation} \tag{ 4.21 }$$

Suppose w₀ is stable against spatially homogeneous perturbations, i.e. λ₁ < 0. In that case, a long-wavelength instability emerges for positive λ₂. Considering (4.9), this instability is reasonable: for λ₂ > 0, inhomogeneities in the particle density will lead to a flow in the direction of the density gradients while for λ₂ < 0 density inhomogeneities decay due to inverse flows. The amplification of density gradients will lead to an agglomeration of particles (cf figure 3(E)). This hypothesis is supported by numerical simulations of the microscopic model, as shown in the next section and in [29]. A similar instability due to the gradient term in (4.9) is also known from other active matter systems [23]. Whereas in [23] 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 λ₂ = 0, i.e.

$$\begin{eqnarray} &&\fl \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji}\, r_{ji}^2 \left ( \mu_{\rm m}(r_{ji}) + \mu_{\rm a}(r_{ji}) \right) = \frac{9 \pi}{16 s_0 \rho_0} \left ( \frac{s_0^2}{2} + \frac{\rho_0 \pi}{2} \int_0^{l_{\rm c}} {\rm d} r_{ji} \, r_{ji}^2 \mu_r(r_{ji}) \right). \end{eqnarray} \tag{ 4.22 }$$

The corresponding analysis of the stability of the spatially homogeneous ordered state for |w₁| > 0 using the full hydrodynamic equations is much more complicated. It is done in great detail for the original Vicsek model in [19] 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 λ₂, η₂ and η₃ may change their algebraic signs. Inserting the ansatz

$$\begin{equation} \rho(\boldsymbol{r},t) = \rho_0 + \delta \rho(\boldsymbol{r},t), \end{equation} \tag{ 4.23a }$$

$$\begin{equation} \boldsymbol{w}(\boldsymbol{r},t) = \boldsymbol{w}_1 + \delta \boldsymbol{w}(\boldsymbol{r},t) \, \end{equation} \tag{ 4.23b }$$

into (4.9) and linearization yields

$$\begin{equation} \fl \frac{\partial \delta \rho}{\partial t} = - \nabla \cdot \delta \boldsymbol{w}, \end{equation} \tag{ 4.24a }$$

$$\begin{eqnarray} &&\fl \frac{\partial \delta \boldsymbol{w}}{\partial t} = \left. \frac{\partial \lambda_1}{\partial \rho}\right |_{\rho_0} \delta \rho \, \boldsymbol{w}_1 - 2 \eta_1 \left (\boldsymbol{w}_1 \cdot \delta \boldsymbol{w} \right) \boldsymbol{w}_1 + \lambda_2 \nabla \delta \rho + \lambda_4 \Delta \delta \boldsymbol{w} + \lambda_5 \nabla (\nabla \cdot \delta \boldsymbol{w}) + \eta_2 \left ( \boldsymbol{w}_1 \cdot \nabla \right) \delta \boldsymbol{w} \nonumber \\ &&+ \eta_3 \left ( \nabla \left ( \boldsymbol{w}_1 \cdot \delta \boldsymbol{w} \right) - \boldsymbol{w}_1 \left ( \nabla \cdot \delta \boldsymbol{w} \right) \right). \end{eqnarray} \tag{ 4.24b }$$

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₁ are pointing in the same direction, because we are interested in the emergence of collectively moving bands (cf figure 3). The growth rate reads

$$\begin{eqnarray} &&\sigma_{\pm} = \textstyle\frac{1}{2} \left [ -2 \lambda_1 - \left (\lambda_4 + \lambda_5 \right) k^2 + {\rm i} \eta_2 \left | \boldsymbol{w}_1 \right | k \pm \sqrt{\alpha_1+{\rm i} \alpha_2} \, \right ], \end{eqnarray} \tag{ 4.25 }$$

where the coefficients

$$\begin{equation} \alpha_1 = \left [ \left ( \lambda_4 + \lambda_5 \right) k^2 + 2 \lambda_1 \right]^2 - k^2 \left [ \left (\eta_2 \left | \boldsymbol{w}_1 \right | \right)^2 - 4 \lambda_2 \right], \end{equation} \tag{ 4.26a }$$

$$\begin{equation} \alpha_2 = - 2 \left | \boldsymbol{w}_1 \right | k \left \{ 2 \left. \frac{\partial \lambda_1}{\partial \rho}\right |_{\rho_0} + \eta_2 \left [ \left ( \lambda_4 + \lambda_5 \right) k^2 + 2 \lambda_1 \right ] \right \} \end{equation} \tag{ 4.26b }$$

were introduced. The real part of the largest eigenvalue reads as follows:

$$\begin{eqnarray} \mbox{Re}\,(\sigma_{+})(k^2) &= -\lambda_1 - \frac{\lambda_4 + \lambda_5}{2} k^2 + \sqrt{\frac{\alpha_1 + \sqrt{\alpha_1^2 + \alpha_2^2}}{8}}. \end{eqnarray} \tag{ 4.27 }$$

Close to the order–disorder transition, λ₁ is small, with the result that it is sufficient to keep the lowest order in 1/λ₁. The following expression is obtained for small wavenumbers [19]:

$$\begin{eqnarray} &&\mbox{Re}\,(\sigma_{+})(k^2) = \frac{\xi_1^2}{8\eta_1 \lambda_1^2} k^2 - \frac{5 \xi_1^4}{128 \eta_1^2 \lambda_1^5} k^4 + \mathcal{O}(k^6). \end{eqnarray} \tag{ 4.28 }$$

That implies that for λ₁ ≳ 0, a long-wavelength instability emerges. By expanding (4.27) in a Taylor series for small k (not assuming that λ₁ ≃ 0), one obtains

$$\begin{eqnarray} &&\mbox{Re}\,(\sigma_{+})(k^2) = \left ( \frac{1}{\eta_1} \left ( \frac{\xi_1}{\lambda_1} + \eta_2 \right)^2 - \frac{\eta_2^2}{\eta_1} + \frac{4 \lambda_2}{\lambda_1} \right) \frac{k^2}{8} + \mathcal{O}(k^4). \end{eqnarray} \tag{ 4.29 }$$

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 7)

$$\begin{equation} \mu_{\rm m} + \mu_{\rm a} = \frac{135 \pi}{64 \rho_0 s_0 \left (l_{\rm s}^3 - l_{\rm c}^3 \right)} \left ( \frac{s_0^2}{2} + \frac{2 \pi \rho_0 l_{\rm c}^3 \mu_r}{21} \right) , \end{equation} \tag{ 4.30a }$$

$$\begin{equation} \mu_{\rm m} - \mu_{\rm a} = \frac{8 D_{\varphi}}{\pi \rho_0 s_0^2 \left (l_{\rm s}^2 - l_{\rm c}^2 \right)} . \end{equation} \tag{ 4.30b }$$

The instability of w₁ for

$$\begin{eqnarray} &&\mu_{\rm m} + \mu_{\rm a} > \frac{135 \pi}{64 \rho_0 s_0 \left (l_{\rm s}^3 - l_{\rm c}^3 \right)} \left ( \frac{s_0^2}{2} + \frac{2 \pi \rho_0 l_{\rm c}^3 \mu_r}{21} \right) \end{eqnarray} \tag{ 4.31 }$$

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 λ₁ = 0 (4.13).

Zoom In Zoom Out Reset image size

As already mentioned, negative values of λ₄ and λ₅ 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 λ₄ and λ₅. 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 with the predictions of the stability analysis using the kinetic description in section 3.2 in figure 7(b). The long-wavelength instability of w₀, leading to clustering for λ₂ > 0 as well as the long-wavelength instability close to the critical point (κ = 2, λ₁ = 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 7(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₁ is not sufficiently well described by the hydrodynamic equations. Please note that the hydrodynamic theory is only valid for λ₄ > 0 (guarantees the stability for large wavenumbers) and 0 ⩽ κ ⩽ 2.617, where the lower bound is equivalent to λ₅ > 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 analysed analytically and it allows us 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 λ₂ and η₂. Furthermore, the structure of the theory suggests that all predictions of the theory by Toner and Tu [15] 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.

4.4. 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 particles within the interaction area is kept fixed: ρ l²_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 [19, 21] 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

$$\begin{equation} \fl \frac{\partial \rho}{\partial t} = - \nabla \cdot \boldsymbol{w}, \end{equation} \tag{ 4.32a }$$

$$\begin{eqnarray} &&\fl \frac{\partial \boldsymbol{w}}{\partial t} + \frac{3 s_0^2 \, \xi_1 }{8 D_{\varphi}} \, (\boldsymbol{w} \cdot \nabla) \, \boldsymbol{w} = \left [ \left ( \xi_1 \rho - \frac{D_{\varphi}}{s_0^2} \right) - \frac{\xi_1^2}{2 D_{\varphi}} \left | \boldsymbol{w} \right |^2 \right]\boldsymbol{w} - \frac{s_0^2}{2} \nabla \rho + \frac{s_0^4}{16 D_{\varphi}} \, \Delta \boldsymbol{w} \nonumber \\ &&+ \frac{5 s_0^2 \, \xi_1}{8 D_{\varphi}} \left ( \frac{1}{2} \nabla \left | \boldsymbol{w} \right |^2 - \boldsymbol{w} \, (\nabla \cdot \boldsymbol{w}) \right). \end{eqnarray} \tag{ 4.32b }$$

Since D_φ is always positive, the dynamics is stable for all microscopic parameters. The coefficient ξ₁ 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' $\nabla \left | \boldsymbol {w} \right |^2$ and w (∇·w). Microscopically, the transition from positive ξ₁ to negative ξ₁ 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 2). Similar conclusions could possibly be drawn from [23, 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.

In order to analyse 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

$$\begin{equation} \left \langle \Phi \right \rangle_t = \left \langle \left | \frac{1}{N s_0} \sum_{i=1}^N \boldsymbol{v}_i(t) \right | \right \rangle_{\! t}. \end{equation} \tag{ 5.1 }$$

Here 〈·〉_t represents a temporal average. $\left \langle \Phi \right \rangle _t=1$ corresponds to perfect orientational order with all agents moving in the same direction, whereas a vanishing $\left \langle \Phi \right \rangle _t$ corresponds to a completely disordered system. Please note that $\left \langle \Phi \right \rangle _t=0$ can only be observed in the thermodynamic limit (N → ∞). In a finite, disordered system we will measure a small, but finite $\left \langle \Phi \right \rangle _t\gtrapprox 0$ due to finite-size fluctuations of the order $1/\sqrt {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

$$\begin{eqnarray} &&S=-\sum_{j, n_j \neq 0} p_j \log p_j = -\sum_{j, n_j \neq 0} \frac{n_j}{N}\log\frac{n_j}{N}, \end{eqnarray} \tag{ 5.2 }$$

where the summation occurs over all occupied cells of the grid with n_j ⩾ 1 (n_j is the number of particles in the jth cell). This allows us to define the following spatial order parameter:

$$\begin{equation} \left \langle \Psi \right \rangle_t = \left \langle 1-\frac{S}{S_{\rm max}} \right \rangle_t \end{equation} \tag{ 5.3 }$$

with S_max being the maximal value of the spatial entropy corresponding to a homogeneous distribution of particles. $\left \langle \Psi \right \rangle _t=0$ corresponds to a perfectly disordered state, whereas $\left \langle \Psi \right \rangle _t{>}0$ indicates a spatially inhomogeneous distribution of particles (clusters, bands), where $\left \langle \Psi \right \rangle _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 $\left \langle \Phi \right \rangle _t$ and $\left \langle \Psi \right \rangle _t$. Thus, we set the upper limit of the colour bar for both order parameters in figures 8 and 9 to 0.2. In consequence, all regions of parameter space, where $\left \langle \Phi \right \rangle _t \geqslant 0.2$ or $\left \langle \Psi \right \rangle _t\geqslant 0.2$, will appear white in figures 8 and 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 8). The first one (μ_a ∼ μ_m) corresponds to the orientational instability and onset of collective motion, cf (4.13) and line (1) in figure 8, 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 8. 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 $\left (\mu _{\mathrm { m}},\mu _{\mathrm { a}} \right ) \in \left \{ \left (\mu _{\mathrm { m}},\mu _{\mathrm { a}} \right ) : |\mu _{\mathrm { m}}|>|\mu _{\mathrm { a}}|, \mu _{\mathrm { m}}>0 \right \}$). This also confirms our previous results of comprehensive numerical investigations of related individual-based models [29, 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 one-particle 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 9. This instability was missed in the previous study of the model [29]. 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 move in either direction along the filament (see figure 10).

Zoom In Zoom Out Reset image size

The hydrodynamic theory derived in section 4 is not able to account for this kind of structure 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 [22, 34, 35, 39, 40].

The self-propelled particle model with selective attraction–repulsion interaction shows a rich dynamical behaviour. In the 'effective alignment' regime, characterized by repulsion from approaching agents and attraction to moving away agents, the model behaviour is closely related to the well-known Vicsek model, in particular for |μ_a| = |μ_m|, μ_a < 0, μ_m > 0. In general, for repulsion from approaching particles equal to, or stronger than, the attraction to particles moving away (|μ_a| ⩾ |μ_m|, μ_a < 0, μ_m > 0), we typically observe the homogeneous flocking phase with giant number fluctuations predicted by Toner and Tu (see e.g. figure 3(B)). 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 PDF. 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 FPE. We have analysed 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 'small wavenumber 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 individual-based model and the macroscopic parameters governing the behaviour of the coarse-grained hydrodynamic variables (density and momentum fields). Interestingly, the hydrodynamic parameters, as for example η₃, which relates to a splay elasticity in corresponding equilibrium systems, or λ₄ and λ₅ which govern the relaxation of splay and bend fluctuations [43], 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 ξ₁, which indicates the transition from an overall aligning effect of the social interaction to the opposite case where the interaction tends 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, such 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 one 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 [18, 20]. 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 one-particle FPE. Similar arguments have been put forward also in [23]. The impact of individual (angular) noise on the mean-field assumption may also be not straightforward. 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 self-propelled 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 also offers an interesting playground for the study of self-organization, pattern formation and phase transitions at far-from-equilibrium conditions.

In this section, it is sketched how the matrix elements $\tilde {K}_{p,n}(\boldsymbol {k})$, defined by

$$\begin{equation} {\hat{K}}_{p}(\boldsymbol{k},t) = \int {\rm d}^2 r \int_0^{2 \pi} {\rm d} \varphi \, F(\boldsymbol{r},\varphi,t) {\rm e}^{{\rm i} \boldsymbol{k}\cdot \boldsymbol{r} + {\rm i} p \varphi}, \end{equation} \tag{ A.1a }$$

$$\begin{equation} {\hat{K}}_{p}(\boldsymbol{k},t) = \sum_{r=-\infty}^{\infty} \tilde{K}_{p,r}(\boldsymbol{k}) {\hat{g}}_{r}(\boldsymbol{k},t) \end{equation} \tag{ A.1b }$$

and used for the stability analysis of the spatial homogeneous solutions (3.11) of the FPE (3.8) in section 3.2, are calculated. $\tilde {K}_{p,n}$ is obtained by inserting the inverse Fourier transform of the one-particle PDF

$$\begin{eqnarray} &&p(\boldsymbol{r},\varphi,t) = \frac{1}{(2 \pi)^3} \sum_{n=-\infty}^{\infty} \int {\rm d}^2 k \, {\hat{g}}_{n}(\boldsymbol{k},t) {\rm e}^{-{\rm i} \boldsymbol{k}\cdot \boldsymbol{r} - {\rm i} n \varphi} \end{eqnarray} \tag{ A.2 }$$

into (3.7) and solving the remaining integral:

$$\begin{eqnarray} &&\fl F_{\varphi}(\boldsymbol{r},\varphi,t) = \frac{1}{(2 \pi)^3} \sum_{p=-\infty}^{\infty} \int \mbox{d}^2 k \, {\hat{g}}_{p}(\boldsymbol{k},t) e^{-i \boldsymbol{k}\cdot \boldsymbol{r}} \int_{\varphi}^{\varphi + 2 \pi} \mbox{d}\varphi_j\, {\rm e}^{-{\rm i} p \varphi_j} \left \{ 2 s_0 \sin \left ( \frac{\varphi - \varphi_j}{2} \right)\right. \, \nonumber \\ &&\times \int_{l_{\rm c}}^{l_{\rm s}} \mbox{d}r_{ji} \, r_{ji} \left [ \mu_{\rm m}(r_{ji}) \int_{\frac{\varphi + \varphi_j}{2}}^{\frac{\varphi + \varphi_j}{2} + \pi} \mbox{d}\alpha \, {\rm e}^{-{\rm i} \boldsymbol{k}\cdot \boldsymbol{r}_{ji}} \sin(\alpha - \varphi) \sin \left ( \frac{\varphi+ \varphi_j}{2} - \alpha \right) \right. \nonumber \\ &&\left. - \mu_{\rm a}(r_{ji}) \int_{\frac{\varphi + \varphi_j}{2}+\pi}^{\frac{\varphi + \varphi_j}{2} + 2 \pi} \mbox{d}\alpha \, {\rm e}^{-{\rm i} \boldsymbol{k}\cdot \boldsymbol{r}_{ji}} \sin(\alpha - \varphi) \sin \left ( \frac{\varphi+ \varphi_j}{2} - \alpha \right) \right] \nonumber \\ &&\left. - \int_0^{l_{\rm c}} {\rm d} r_{ji} \, r_{ji} \, \mu_r(r_{ji}) \int_0^{2 \pi} {\rm d} \alpha \, \sin (\alpha - \varphi) {\rm e}^{-{\rm i} \boldsymbol{k}\cdot \boldsymbol{r}_{ji}} \right \}. \end{eqnarray} \tag{ A.3 }$$

In order to perform the integration over α and φ_j, the following identity is used to express the exponential function:

$$\begin{eqnarray} &&{\rm e}^{-{\rm i} \boldsymbol{k}\cdot \boldsymbol{r}_{ji}} = J_0(k \, r_{ji}) + 2 \sum_{s=1}^{\infty} (-{\rm i})^s J_{\rm s}(k \, r_{ji}) \cos \left (s \left ( \alpha-\chi \right) \right), \end{eqnarray} \tag{ A.4 }$$

where $\boldsymbol {k}=k\left (\cos \chi , \sin \chi \right )$ and $\boldsymbol {r}_{ji} = r_{ji} \left ( \cos \alpha , \sin \alpha \right )$. 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 (A.1b) of the Fourier coefficients ${\hat {g}}_{r}(\boldsymbol {k},t)$, where one can read off the elements of the infinite-dimensional matrix $\tilde {K}_{p,n}(\boldsymbol {k})$:

$$\begin{eqnarray} &&\fl \tilde{K}_{p,n}(\boldsymbol{k}) = \frac{s_0 \pi^2}{2 i} \left \{ \left (\tilde{\mu}_{m,0}(k)-\tilde{\mu}_{a,0}(k) \right) \left ( \delta_{p,1} \delta_{n,1} - \delta_{p,-1} \delta_{n,-1} \right) + \left (\tilde{\mu}_{m,2}(k)-\tilde{\mu}_{a,2}(k) \right) \right. \nonumber \\ &&\times \left. \left [ {\rm e}^{-2i\chi} \left ( \delta_{p,-1} \, \delta_{n,1} - \delta_{p,-2} \, \delta_{n,0} \right) - {\rm e}^{2 i \chi} \left ( \delta_{p,1} \, \delta_{n,-1} - \delta_{p,2} \, \delta_{n,0} \right) \right ] \right \} \nonumber \\ &&+\, 2 \pi^2 \tilde{\mu}_r(k) \, \delta_{n,0} \left ( {\rm e}^{{\rm i} \chi} \, \delta_{p,1} - e^{-i \chi} \, \delta_{p,-1} \right) + 2 s_0 \sum_{s=0}^{\infty} (-1)^s \left ( \tilde{\mu}_{m,2s+1}(k) + \tilde{\mu}_{a,2s+1}(k) \right) \nonumber \\ &&\times \left [ \Lambda_{p,s} \, {\rm e}^{-{\rm i}(2s+1)\chi} \, \delta_{p+2s+1,n} - \Gamma_{p,s} \, {\rm e}^{{\rm i}(2s+1)\chi} \, \delta_{p-(2s+1),n} \right ]. \end{eqnarray} \tag{ A.5 }$$

The following coefficients were introduced:

$$\begin{equation} \fl \tilde{\mu}_{m,\nu}(k) = \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji} \, \mu_{\rm m}(r_{ji}) J_{\nu}(k \, r_{ji}), \end{equation} \tag{ A.6a }$$

$$\begin{equation} \fl \tilde{\mu}_{a,\nu}(k) = \int_{l_{\rm c}}^{l_{\rm s}} {\rm d} r_{ji} \, r_{ji} \, \mu_{\rm a}(r_{ji}) J_{\nu}(k \, r_{ji}), \end{equation} \tag{ A.6b }$$

$$\begin{equation} \fl \tilde{\mu}_r(k) = \int_0^{l_{\rm c}} {\rm d} r_{ji} \, r_{ji} \, \mu_r(r_{ji}) J_1(k \, r_{ji}), \end{equation} \tag{ A.6c }$$

$$\begin{equation} \fl \Lambda_{p,s} = \frac{32 p}{(2s-1)(1+2s)(3+2s)(-1+2p+2s)(1+2p+2s)(3+2p+2s)}, \end{equation} \tag{ A.6d }$$

$$\begin{equation} \fl \Gamma_{p,s} = \frac{32 p}{(2s-1)(1+2s)(3+2s)(-1+2p-2s)(1+2p-2s)(-3+2p-2s)}. \end{equation} \tag{ A.6e }$$

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.

Besides the coefficient η₁, which is always positive, all transport coefficients may be of either sign, positive or negative. In the following, the sign of every coefficient is discussed in the (μ_m,μ_a)-parameter plane. For simplicity, it is assumed that the interaction strengths μ_a and μ_m are constant for l_c < r_ji < l_s and zero otherwise.

λ₁ is positive above a critical line (positive slope and ordinate-intercept γ₁) which is given by

$$\begin{eqnarray} &&\mu_{\rm m} - \mu_{\rm a} = \frac{8 D_{\varphi}}{\pi \rho_0 s_0^2 (l_{\rm s}^2 - l_{\rm c}^2)} \; {=:} \; \gamma_1 > 0, \end{eqnarray} \tag{ B.1 }$$

whereas λ₂ is positive above a critical line (negative slope and positive ordinate-intercept γ₂), given by

$$\begin{eqnarray} &&\mu_{\rm m} + \mu_{\rm a} = \frac{27 \pi}{16 s_0 \rho_0 (l_{\rm s}^3 - l_{\rm c}^3)} \left ( \frac{s_0^2}{2} + \frac{\pi \rho_0 \mu_r l_{\rm c}^3}{6} \right) \; {=:} \; \gamma_2 \geqslant \frac{9 \pi}{32 s_0} \frac{l_{\rm c}^3}{l_{\rm s}^3 - l_{\rm c}^3} \, \mu_r > 0. \end{eqnarray} \tag{ B.2 }$$

The analysis of the coefficient λ₄ is somewhat more involved. λ₄ is greater than zero inside a parabola, which is oriented along the line

$$\begin{eqnarray} &&\mu_{\rm m} + \mu_{\rm a} = \frac{405 \pi s_0}{256 \rho_0 \left ( l_{\rm s}^3 - l_{\rm c}^3 \right)} \; {=:} \; \gamma_3 > 0. \end{eqnarray} \tag{ B.3 }$$

The width of the parabola is proportional to the noise strength D_φ and inverse proportional to the density ρ₀. The vertex of the parabola has the coordinates

$$\begin{equation} \mu_{\rm a}^{(v)} = \frac{1}{2} \left ( \frac{405 \pi s_0}{256 \rho_0 (l_{\rm s}^3 - l_{\rm c}^3)} + \frac{25 s_0^4}{2 \pi \rho_0 (l_{\rm s}^4-l_{\rm c}^4) D_{\varphi}} \right ), \end{equation} \tag{ B.4a }$$

$$\begin{equation} \mu_{\rm m}^{(v)} = \frac{1}{2} \left ( \frac{405 \pi s_0}{256 \rho_0 (l_{\rm s}^3 - l_{\rm c}^3)} - \frac{25 s_0^4}{2 \pi \rho_0 (l_{\rm s}^4-l_{\rm c}^4) D_{\varphi}} \right ). \end{equation} \tag{ B.4b }$$

The parabola itself is given by⁶

$$\begin{eqnarray} &&\fl \mu_{\rm m} - \mu_{\rm a} = \frac{8 s_0^4}{\pi \rho_0 (l_{\rm s}^4 - l_{\rm c}^4) D_{\varphi}} \left [ \left ( \frac{64 \rho_0 (l_{\rm s}^3 - l_{\rm c}^3) }{135 \pi s_0} \right)^2 (\mu_{\rm m} + \mu_{\rm a})^2 - \frac{32 \rho_0 (l_{\rm s}^3 - l_{\rm c}^3) }{45 \pi s_0} \, (\mu_{\rm m} + \mu_{\rm a}) - 1 \right ].\nonumber\\ \end{eqnarray} \tag{ B.5 }$$

The coefficient λ₅ is positive above the main diagonal, i.e. μ_m > μ_a. Transport coefficient η₂ is less than zero in between two hyperbolas, which are bounded by the asymptotes (B.1) and (B.3). The hyperbolas are determined by

$$\begin{eqnarray} &&\mu_{\rm m}-\mu_{\rm a} = \gamma_3 \, \frac{\left (\mu_{\rm m} + \mu_{\rm a}\right)}{\left (\mu_{\rm m} + \mu_{\rm a}\right) - \gamma_1} \end{eqnarray} \tag{ B.6 }$$

or rather

$$\begin{eqnarray} &&\mu_{\rm m} = \frac{1}{2} \left (\gamma_1 + \gamma_3 \pm \sqrt{4 \mu_{\rm a}^2 + 4 \mu_{\rm a} (\gamma_1 - \gamma_3) + (\gamma_1 + \gamma_3)^2} \right). \end{eqnarray} \tag{ B.7 }$$

The extrema of the hyperbolas are located at

$$\begin{equation} \left (\mu_{\rm a}^{({\rm ext})},\mu_{\rm m}^{(\rm min)} \right ) = \left ( \frac{\gamma_3 - \gamma_1}{2}, \frac{\left ( \sqrt{\gamma_1} - \sqrt{\gamma_3} \right)^2}{2} \right), \end{equation} \tag{ B.8a }$$

$$\begin{equation} \left (\mu_{\rm a}^{({\rm ext})},\mu_{\rm m}^{({\rm max})} \right ) = \left ( \frac{\gamma_3 - \gamma_1}{2}, \frac{\left ( \sqrt{\gamma_1} + \sqrt{\gamma_3} \right)^2}{2} \right). \end{equation} \tag{ B.8b }$$

The critical line for η₃ = 0 depends on the noise strength, whereby the critical noise strength reads

$$\begin{eqnarray} &&D_{\varphi}^{(\mathrm{c})} = \frac{675 \pi^2 s_0^3 (l_{\rm s}^2-l_{\rm c}^2)}{2048 (l_{\rm s}^3-l_{\rm c}^3)}. \end{eqnarray} \tag{ B.9 }$$

η₃ is positive, if

$$\begin{eqnarray} \mu_{{\rm m}} \left\{\begin{array}{@{}lc} > \mu_{{\rm a}} \displaystyle\frac{D_{\varphi}^{(\mathrm{c})} + D_{\varphi}}{D_{\varphi}^{(\mathrm{c})} - D_{\varphi}}, & D_{\varphi} < D_{\varphi}^{(\mathrm{c})}, \\ < \mu_{{\rm a}} \displaystyle\frac{D_{\varphi}^{(\mathrm{c})} + D_{\varphi}}{D_{\varphi}^{(\mathrm{c})} - D_{\varphi}}, & D_{\varphi} > D_{\varphi}^{(\mathrm{c})}. \end{array} \right. \end{eqnarray} \tag{ B.10 }$$

Please note that the slope of the lines is bounded by one (D_φ = 0) and minus one (D_φ → ∞), respectively. Hence, η₃ is always positive for $\left \{ (\mu _{\mathrm { a}},\mu _{\mathrm { m}}): \mu _{\mathrm { a}}<0, -\left | \mu _{\mathrm { a}} \right | < \mu _{\mathrm { m}} < \left | \mu _{\mathrm { a}} \right | \right \}$ and η₃ always negative for $\left \{ (\mu _{\mathrm { a}},\mu _{\mathrm { m}}): \mu _{\mathrm { a}}>0, -\mu _{\mathrm { a}} < \mu _{\mathrm { m}} < \mu _{\mathrm { a}} \right \}$. All criteria discussed above are summarized in figure B.1.

Zoom In Zoom Out Reset image size

The stability of the hydrodynamic equations (4.9) requires that both λ₄ and λ₅ 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.

In [24, 26, 29], the symmetric temperature tensor

$$\begin{eqnarray} \hat{T} = \left ( \begin{array}{@{}cc@{}} \left \langle \cos \varphi^2 \right \rangle - \left \langle \cos \varphi \right \rangle^2 & \left \langle \sin \varphi \cos \varphi \right \rangle - \left \langle \sin \varphi \right \rangle\left \langle \cos \varphi \right \rangle\\ \left \langle \sin \varphi \cos \varphi \right \rangle - \left \langle \sin \varphi \right \rangle\left \langle \cos \varphi \right \rangle & \left \langle \sin \varphi^2 \right \rangle - \left \langle \sin \varphi \right \rangle^2 \end{array} \right ) \end{eqnarray} \tag{ C.1 }$$

is defined. It is a measure for the width of the velocity distribution (fluctuations around the mean velocity). The temperature is basically related to the second Fourier coefficients:

$$\begin{equation} T_{xx} = \frac{1}{2}\left (1 + \frac{{\skew6\hat{f}}_{2}(\boldsymbol{r},t)+{\skew6\hat{f}}_{-2}(\boldsymbol{r},t)}{2 \rho(\boldsymbol{r},t)} \right) - \left ( \frac{{\skew6\hat{f}}_{1}(\boldsymbol{r},t) + {\skew6\hat{f}}_{-1}(\boldsymbol{r},t)}{2 \rho(\boldsymbol{r},t)} \right)^2, \end{equation} \tag{ C.2a }$$

$$\begin{equation} T_{xy} = \frac{{\skew6\hat{f}}_{2}(\boldsymbol{r},t)-{\skew6\hat{f}}_{-2}(\boldsymbol{r},t)}{4 i \rho(\boldsymbol{r},t)} - \frac{\left ({\skew6\hat{f}}_{1}(\boldsymbol{r},t)\right)^2 - \left ({\skew6\hat{f}}_{-1}(\boldsymbol{r},t)\right)^2}{4 i \left (\rho(\boldsymbol{r},t)\right)^2}, \end{equation} \tag{ C.2b }$$

$$\begin{equation} T_{yy} = \frac{1}{2}\left (1 - \frac{{\skew6\hat{f}}_{2}(\boldsymbol{r},t)+{\skew6\hat{f}}_{-2}(\boldsymbol{r},t)}{2 \rho(\boldsymbol{r},t)} \right) + \left ( \frac{{\skew6\hat{f}}_{1}(\boldsymbol{r},t) - {\skew6\hat{f}}_{-1}(\boldsymbol{r},t)}{2 \rho(\boldsymbol{r},t)} \right)^2. \end{equation} \tag{ C.2c }$$

The second Fourier coefficients are assumed to be fast variables in section 4.2, so that the temperature is implicitly contained in the hydrodynamic description of the system, even though the dynamics of the temperature tensor is not derived explicitly in this context.